subroutine spllay
C ***************************************************************** wjr
C Converts NASIS layered IFC files into 10,40,50,... IFC files
C
C     Edit History
C     07-Feb-01   wjr   created
C
      include 'p1werm.inc'
      include 'wpath.inc'
      include 'm1subr.inc'
      include 'm1sim.inc'
      include 'm1geo.inc'
      include 'm1flag.inc'
      include 'm1dbug.inc'
      include 's1layr.inc'
      include 's1surf.inc'
      include 's1phys.inc'
      include 's1agg.inc'
      include 's1dbh.inc'
      include 's1dbc.inc'
      include 's1sgeo.inc'
      include 'h1hydro.inc'
      include 'h1scs.inc'
      include 'h1db1.inc'
      include 'file.fi'

c     + + + LOCAL COMMON BLOCKS + + +
      include 'main/main.inc'
c     + + + LOCAL VARIABLES + + +
c      integer       lay
C      character     line*256
      integer       isr
	  real			totthk, curdep, tgtthk
	  integer		otnlay(mnsz)
	  real			newthk(mnsz)
	  integer		newnumlay
	  integer		oldcur
	  integer		newcur
	  integer		dolcur
	  integer		ldx

c hot fix
c     set up target layers (use as max for now)
      real          targetthk(mnsz)
      targetthk(1) = 10
      targetthk(2) = 40
      targetthk(3) = 50
      targetthk(4) = 50
      targetthk(5) = 50
      targetthk(6) = 75
      targetthk(7) = 75
      do ldx=8,mnsz
          targetthk(ldx) = 100
       end do

c alternative layering
      targetthk(1) = 10
      do ldx=2,mnsz
          targetthk(ldx) = targetthk(ldx-1) * 1.5
       end do

      do isr = 1,nsubr

      oldcur = 0
      newcur = 0
      totthk = 0
      tgtthk = 0
   10 continue
        oldcur = oldcur + 1
        totthk = totthk + aszlyt(oldcur,isr)
        newcur = newcur + 1
        tgtthk = tgtthk + targetthk(newcur)
        otnlay(newcur) = oldcur
   20   if( tgtthk .lt. totthk ) then
          newcur = newcur + 1
          tgtthk = tgtthk + targetthk(newcur)
          otnlay(newcur) = oldcur
          go to 20
        end if
		if (oldcur .lt. nslay(isr)) goto 10
      nslay(isr) = newcur
      do ldx = 1, nslay(isr)
        aszlyt(ldx, isr) = targetthk(ldx)
      enddo
      go to 1000
c end of hot fix

c     + + + FUNCTION DECLARATIONS + + +
C
c     convert layers by subregion
C
c      do isr = 1,nsubr

        write(*,*) 'ready to spllay:'

		totthk = 0;
		do ldx = 1, nslay(isr)
			totthk = totthk + aszlyt(ldx, isr)
		enddo
C
C	  stop if soil sample is less than 40mm thick
C
		if (totthk .lt. 40.0) then
			write(*,9002) 
9002		format('Total soil thickness too thin for WEPS')
			stop 14
		endif
		  
		oldcur = 1
		newcur = 1
		dolcur = 1
		if (aszlyt(dolcur,isr) .ge. 150) then
		  newthk(1) = 10
		  newthk(2) = 40
		  newthk(3) = 50
		  otnlay(1) = oldcur
		  otnlay(2) = oldcur
		  otnlay(3) = oldcur
		  aszlyt(1, isr) = aszlyt(1, isr) - 100
          newcur = 4
		  oldcur = 1
        elseif (aszlyt(dolcur,isr) .ge. 120) then
		  newthk(1) = aszlyt(1, isr) * 1 / 15
		  newthk(2) = aszlyt(1, isr) * 4 / 15
		  newthk(3) = aszlyt(1, isr) * 5 / 15
		  newthk(4) = aszlyt(1, isr) * 5 / 15
		  otnlay(1) = oldcur
		  otnlay(2) = oldcur
		  otnlay(3) = oldcur
		  otnlay(4) = oldcur
		  newcur = 5
		  oldcur = 2
		elseif (aszlyt(dolcur,isr) .ge. 70) then
		  newthk(1) = aszlyt(1, isr) * 1 / 10
		  newthk(2) = aszlyt(1, isr) * 4 / 10
		  newthk(3) = aszlyt(1, isr) * 5 / 10
		  otnlay(1) = oldcur
		  otnlay(2) = oldcur
		  otnlay(3) = oldcur
		  newcur = 4
		  oldcur = 2
		elseif(aszlyt(dolcur,isr) .ge. 20) then
		  newthk(1) = aszlyt(1, isr) * 1 / 5
		  newthk(2) = aszlyt(1, isr) * 4 / 5
		  otnlay(1) = oldcur
		  otnlay(2) = oldcur
		  newcur = 3
		  oldcur = 2
		else
		  newthk(1) = aszlyt(dolcur, isr) * 1 / 5
		  otnlay(1) = oldcur
		  newcur = 2
		  oldcur = 2
		  if (newthk(1).lt.4) newthk(1) = 4
		endif
		
		dolcur = oldcur
		if (newcur .eq. 2) then
	 	  if (aszlyt(dolcur,isr) .ge. newthk(1) * 14) then
		    newthk(2) = newthk(1) * 4
			newthk(3) = newthk(1) * 5
			otnlay(2) = oldcur
			otnlay(3) = oldcur
			aszlyt(dolcur,isr) = aszlyt(dolcur,isr) - newthk(2) - 
     *        newthk(3)
			newcur = 4
			oldcur = 2
		  elseif(aszlyt(dolcur,isr) .ge. newthk(1) * 11) then
		    newthk(2) = newthk(1) * 4 / 14
			newthk(3) = newthk(1) * 5 / 14
			newthk(4) = newthk(1) * 5 / 14
			otnlay(2) = oldcur
			otnlay(3) = oldcur
			otnlay(4) = oldcur
			newcur = 5
			oldcur = 3
		  elseif(aszlyt(dolcur,isr) .ge. newthk(1) * 6) then
		    newthk(2) = newthk(1) * 4 / 9
			newthk(3) = newthk(1) * 5 / 9
			otnlay(2) = oldcur
			otnlay(3) = oldcur
			newcur = 4
			oldcur = 3
		  else
		    newthk(2) = aszlyt(dolcur,isr)
			otnlay(2) = oldcur
			newcur = 3
			oldcur = 3
          endif
        endif

		dolcur = oldcur
		if(newcur.eq.3) then
		  if (aszlyt(dolcur,isr) .ge. newthk(2)*1.25 + 50) then
		     newthk(3) = newthk(2) * 1.25
			 aszlyt(dolcur,isr) = aszlyt(dolcur,isr) - 50
			 otnlay(3) = oldcur
			 newcur = 4
			 oldcur = 3
		  elseif (aszlyt(dolcur,isr) .ge. newthk(2)*1.25) then
		     newthk(3) = newthk(2) * .5
		     newthk(4) = newthk(2) * .5
			 otnlay(3) = oldcur
			 otnlay(4) = oldcur
			 newcur = 5
			 oldcur = 3
		  else
		     newthk(3) = aszlyt(dolcur,isr)
			 if (newthk(3) .lt. 10) newthk(3) = 10
			 otnlay(3) = oldcur
			 newcur = 5
			 oldcur = 3
		  endif
		endif

		curdep = 0
		do ldx = 1, newcur-1
		  curdep = curdep + newthk(ldx)
		enddo

C
C layers in top 200 mm are targeted to 50 mm thickness
C
		do while (curdep .lt. 200)
		  if (aszlyt(oldcur, isr) .gt. 100) then
		    newthk(newcur) = 50
			aszlyt(oldcur, isr) = aszlyt(oldcur, isr) - 50
			otnlay(newcur) = oldcur
			newcur = newcur + 1
			curdep = curdep + 50
		  else if (aszlyt(oldcur, isr) .gt. 50) then
			newthk(newcur) = aszlyt(oldcur, isr) / 2.0
			newthk(newcur+1) = newthk(newcur)
			otnlay(newcur) = oldcur
			otnlay(newcur+1) = oldcur
			newcur = newcur + 2
			oldcur = oldcur + 1
			if (oldcur .gt. nslay(isr)) exit
		  else
			newthk(newcur) = aszlyt(oldcur, isr)
			curdep = curdep + newthk(newcur)
			otnlay(newcur) = oldcur
			newcur = newcur + 1
			oldcur = oldcur + 1
			if (oldcur .gt. nslay(isr)) exit
		  endif
		enddo
C
C no more input layers left
C
		if (oldcur .gt. nslay(isr)) goto 120	
C
C layers in top 200-400 mm are targeted to 50-100 mm thickness
C   depending on depth
C
		do while (curdep .lt. 400)
		  tgtthk = 50 + ((curdep - 200) / 200) * 50
		  if (aszlyt(oldcur, isr) .gt. tgtthk * 2) then
		    newthk(newcur) = tgtthk
			aszlyt(oldcur, isr) = aszlyt(oldcur, isr) - 50
			otnlay(newcur) = oldcur
			newcur = newcur + 1
			curdep = curdep + 50
		  else if (aszlyt(oldcur, isr) .gt. tgtthk) then
			newthk(newcur) = aszlyt(oldcur, isr) / 2.0
			newthk(newcur+1) = newthk(newcur)
			otnlay(newcur) = oldcur
			otnlay(newcur+1) = oldcur
			newcur = newcur + 2
			oldcur = oldcur + 1
			if (oldcur .gt. nslay(isr)) exit
		  else
			newthk(newcur) = aszlyt(oldcur, isr)
			curdep = curdep + newthk(newcur)
			otnlay(newcur) = oldcur
			newcur = newcur + 1
			oldcur = oldcur + 1
			if (oldcur .gt. nslay(isr)) exit
		  endif
		enddo
C
C  no more input layers left
C
		if (oldcur .gt. nslay(isr)) goto 120	
C
C layers lower than 400 mm are targeted at 100 mm
C
		do while (oldcur .le. nslay(isr))
		  if (aszlyt(oldcur, isr) .gt. 200) then
			newthk(newcur) = 100
			aszlyt(oldcur, isr) = aszlyt(oldcur, isr) - 100
			otnlay(newcur) = oldcur
			newcur = newcur + 1
		  else
			newthk(newcur) = aszlyt(oldcur, isr) / 2
			newthk(newcur+1) = aszlyt(oldcur, isr) / 2
			otnlay(newcur) = oldcur
			otnlay(newcur+1) = oldcur
			newcur = newcur + 2
			oldcur = oldcur + 1
		  endif
		enddo
C
C update number of layers
C
		nslay(isr) = newcur -1
C
C transfer new layer thicknesses to global array
C
		do ldx = 1, nslay(isr)
		  aszlyt(ldx, isr) = newthk(ldx)
		enddo
C
C      update soil properties for new layers
C
c     note!, this is done backward assuming that more layers are created
c     than exist in the original soil layering and the old value
c     will be used before it is overwritten
 1000 continue
	  do ldx = nslay(isr), 2, -1
		asfsan(ldx, isr) = asfsan(otnlay(ldx), isr)
		asfsil(ldx, isr) = asfsil(otnlay(ldx), isr)
		asfcla(ldx, isr) = asfcla(otnlay(ldx), isr)
		asvroc(ldx, isr) = asvroc(otnlay(ldx), isr)
		asfcs(ldx, isr)  = asfcs(otnlay(ldx), isr)
		asfms(ldx, isr)  = asfms(otnlay(ldx), isr)
		asffs(ldx, isr)  = asffs(otnlay(ldx), isr)
		asfvfs(ldx, isr) = asfvfs(otnlay(ldx), isr)
		asfwdc(ldx, isr) = asfwdc(otnlay(ldx), isr)
		asdblk(ldx, isr) = asdblk(otnlay(ldx), isr)
		asdwblk(ldx, isr)= asdwblk(otnlay(ldx), isr)
		aslagm(ldx, isr) = aslagm(otnlay(ldx), isr)
		as0ags(ldx, isr) = as0ags(otnlay(ldx), isr)
		aslagx(ldx, isr) = aslagx(otnlay(ldx), isr)
		aslagn(ldx, isr) = aslagn(otnlay(ldx), isr)
		asdagd(ldx, isr) = asdagd(otnlay(ldx), isr)
		aseags(ldx, isr) = aseags(otnlay(ldx), isr)
		ahrwc(ldx, isr)  = ahrwc(otnlay(ldx), isr)
		ahrwcs(ldx, isr) = ahrwcs(otnlay(ldx), isr)
		ahrwcf(ldx, isr) = ahrwcf(otnlay(ldx), isr)
		ahrwcw(ldx, isr) = ahrwcw(otnlay(ldx), isr)
		ahrwc1(ldx, isr) = ahrwc1(otnlay(ldx), isr)
		ah0cb(ldx, isr)  = ah0cb(otnlay(ldx), isr)
		aheaep(ldx, isr) = aheaep(otnlay(ldx), isr)
		ahrsk(ldx, isr)  = ahrsk(otnlay(ldx), isr)
		asfom(ldx, isr)  = asfom(otnlay(ldx), isr)
		asdsblk(ldx, isr)= asdsblk(otnlay(ldx), isr)
		asdpart(ldx, isr)= asdpart(otnlay(ldx), isr)
		as0ph(ldx, isr)  = as0ph(otnlay(ldx), isr)
		asfcce(ldx, isr) = asfcce(otnlay(ldx), isr)
		asfcec(ldx, isr) = asfcec(otnlay(ldx), isr)
		asfsmb(ldx, isr) = asfsmb(otnlay(ldx), isr)
		as0ec(ldx, isr)  = as0ec(otnlay(ldx), isr)
		asrsar(ldx, isr) = asrsar(otnlay(ldx), isr)
		asftan(ldx, isr) = asftan(otnlay(ldx), isr)
		asftap(ldx, isr) = asftap(otnlay(ldx), isr)
      enddo
 120  continue
      enddo
C
      end