	Subroutine MGAR

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	Written for  	PhD Dissertation                                                c
c	Last Updated	November 2007													c
c	Written by	Leslie Gowdish														c
c			Ag. & Bio. Eng. Department											    c
c			University of Florida													c
c			Email Leslie@Gowdish.com    											c
c																					c
c	Program 1-dimensional Green-Ampt infiltration and redistribution method (MGAR)  c
c	Code is written in Fortran language											    c	
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	
	USE ModuleDefs								
	Implicit None

	Integer j, precipCount,T
	Integer k, i, runFlag
	Double Precision dM,t1
	Double Precision tstart,tend
	Double Precision top(maxNumLayers), bottom(maxNumLayers)

	
	Parameter (RunFlag = 1)

	Open(Unit = 500, file = 
     +	'Soil_In.txt',status = 'old')
     	Open(Unit = 1000, file = 
     +	'Rain_In.txt',status = 'old')
      
	Open(Unit = 200, file = '1DInfilOutput.txt')
	Open(unit = 300, file = '1DInfilOutput2.txt')
	
	Read(500,*) Ksat,Sav,wcs,wci,wcr,psiB,lambda,smax,wcmin
	Read(500,*) NumObsLayers,PrecipCount,deltaT
	Read(500,*) (ProfileDepth(j), j = 1,NumObsLayers)
	Read(500,*) (TopObsLayers(j), j = 1,NumObsLayers)
	Read(1000,*) (precips(j), j = 1,PrecipCount)
	
	Do j = 1, NumObsLayers
		top(j) = (TopObsLayers(j))
		bottom(j) = (TopObsLayers(j) + ProfileDepth(j))
	EndDo
		
	
	IF (RunFlag.EQ.1) THEN
	
		Call Initialize
		tstart = 0.0
		tend = tstart + deltaT
		
		Write(200,8500) 'Time', 'F1', 'Z1', 'WC1', 'F2', 'Z2', 'WC2',
     +		 'F3', 'Z3', 'WC3', 'F4', 'Z4', 'WC4', 'WCO', 
     +        'RWC', 'CummF',(Top(j),'-',Bottom(j),j = 1,numObsLayers)
		Write(200,8600)
!		Write(200,7500) 'Time', 'F1', 'Z1', 'WC1', 'F2', 'Z2', 'WC2',
!     +		 'F3', 'Z3', 'WC3', 'F4', 'Z4', 'WC4', 'WCO', 
!     +        'RWC', 'CummF','0-500mm','500-1000mm','0-1000mm'

	EndIf
	
	Do T = 1, precipCount - 1
		precipRate = precips(T)
		Call CalcInfilRate(tstart,tend)
		tstart = tstart + deltaT
		tend = tend + deltaT
!		Write(200,7000) tstart,WF(1)%Famt,WF(1)%z,WF(1)%WC
!		Write(200,7000) tstart,(LayerAvgWC(i),i=1,3)
	Write(200,8000) tstart, (WF(I)%FAmt, WF(I)%z, WF(I)%WC, i = 1,4),
     +      WCO,RWC,CummInfil,(LayerAvgWC(j),j=1,NumObsLayers)			

	EndDo
7500  Format (A6, 17A11,/,200('-'))
8000	Format(F6.2, 18F11.3)
8500  Format (A6, 15A11,50(F6.0,A1,F6.0))
8600  Format (210('-'))
!7000	Format (F6.2,F10.6,F12.6,F9.6)
!7000	Format (F6.2,3F9.6)
	END


	
	
	
	Subroutine Initialize
	USE ModuleDefs
	Implicit None

	Integer T, L

ccccccc initialize calculated variables ccccccccccccccccccccccccccccccc
   
      pondFlag		= 0
      redistStatus    = 10
      WCO				= 0.0
      RWC				= 0.0
      fp				= 0.0
      Bfp				= 0.0
      timeToPond		= 0.0
      tpp				= 0.0
      precipRate      = 0.0
      cummPrecip      = 0.0
      pondingAmt      = 0.0
      infilRate       = 0.0
      cummInfil       = 0.0
      runoffAmt       = 0.0
    
      numRedist		= 0
      redistTime		= 0.0
      time			= 0.0
      tpAdj			= 0.0
   
cccccc  loop for initializing wetting front data ccccccccccccccccccccccccccccccccccccc
	DO T = 1,4
		WF(T)%z			 = 0.0
		WF(T)%FAmt		 = 0.0
		WF(T)%WC		 = 0.0
		WF(T)%WCHold	 = 0.0
		WF(T)%RedistTime = 0.0
		WF(T)%adjFactor  = 0.0
		WF(T)%NumRedist	 = 0
	ENDDO
  
      WF(1)%WC = WCI
      WF(1)%WCHold = WCI


cccccc  Loop to initialize observation layers and average water content ccccccccccccc
      DO L = 1, NumObsLayers
		LayerBottomZ(L) = TopObsLayers(L) + ProfileDepth(L)
		LayerAvgWC(L)   = 0.0
      ENDDO
    	
	
	END

	Subroutine CalcInfilRate(tstart,tend)
	USE ModuleDefs
	Implicit None
	Double Precision tstart, tend

	Double Precision precip, totalInfil
	Double Precision t1, intertime, deltaTerr
	Double Precision a, b, c
	Integer T, nsteps, i 
	

	precip = precipRate

cccccc update number of redistribution phases and redistribution times ccccccccccccccc
	IF ((redistStatus.EQ.10.OR.redistStatus.EQ.20.OR.redistStatus.EQ.
     +		30.OR.redistStatus.EQ.40).AND.(precip.LT.Ksat).AND.
     +		(pondingAmt.LE.1E-10)) THEN 
		
		numRedist = numRedist + 1
		IF (redistStatus.EQ.10) THEN
			redistTime = 0.0
              WF(1)%redistTime = 0.0
			WF(1)%numRedist = numRedist
		ELSE IF (redistStatus.EQ.20) THEN
			redistTime = 0.0
			WF(2)%redistTime = 0.0
			WF(2)%numRedist = numRedist
		ELSE IF (redistStatus.EQ.30) THEN
			redistTime = 0.0
			WF(3)%redistTime = 0.0
			WF(3)%numRedist = numRedist
		ELSE IF (redistStatus.EQ.40) THEN
		!wf1 and 2 merge; old 1 and 2 are wf1; old3 is wf2; old4 is wf3;
			WF(1)%redistTime = WF(2)%redistTime
			WF(1)%numRedist = WF(2)%numRedist
			WF(2)%redistTime = WF(3)%redistTime
			WF(2)%numRedist = WF(3)%numRedist
			WF(3)%redistTime = 0.0
			WF(3)%numRedist = numRedist
		ENDIF
	ENDIF

	redistTime = redistTime + deltaT
	
	DO T = 1,3
		WF(T)%redistTime = WF(T)%redistTime + 1
	ENDDO
	
	IF ((precip.GT.0.0).OR.(pondingAmt.GT.1E-10)) THEN 
	!BEGIN
cccccccccccccccccccc Start of infiltration through soil cccccccccccccccccccccccccccccc
		IF (pondingAmt.LT.1E-10) THEN
		!BEGIN
		!ccccccc update redistribution status cccccccccccccccccccccccccccccccccccccc
			IF ((redistStatus.EQ.10).AND.(precip.LT.ksat)) THEN
			!BEGIN //wetting front can redistribute
				redistStatus = 11
			!END
			ELSE IF ((redistStatus.EQ.11).AND.(precip.GT.ksat)) THEN
			!BEGIN//two wetting fronts with one redistributing
				redistStatus = 20
				WF(2)%FAmt = cummInfil - WF(1)%FAmt
				WF(2)%WC = WCS
				WF(2)%WCHold = WCS
				fp = 0.0
				tpAdj = 0
			!END
			ELSE IF ((redistStatus.EQ.20).AND.(precip.LT.ksat)) THEN
			!BEGIN //two wetting fronts with both redistributing
				redistStatus = 21
			!END
			ELSE IF ((redistStatus.EQ.21).AND.(precip.GT.ksat)) THEN
			!BEGIN//three wetting fronts with two redistributing
				redistStatus = 30
				WF(3)%FAmt = cummInfil - WF(2)%FAmt - WF(1)%FAmt
                	WF(3)%WC = WCS
				WF(3)%WCHold = WCS
				fp = 0.0
				tpAdj = 0
			!END
			ELSE IF ((redistStatus.EQ.30).AND.(precip.LT.ksat)) THEN
			!BEGIN //three wetting fronts with three redistributing
				redistStatus = 31
			!END
			ELSE IF ((redistStatus.EQ.31).AND.(precip.GT.ksat)) THEN
			!BEGIN//four wetting fronts with 3 redistributing
				redistStatus = 40
				WF(4)%FAmt = cummInfil - WF(3)%FAmt - WF(2)%FAmt - 
     +							WF(1)%FAmt
                  WF(4)%WC = WCS
				WF(4)%WCHold = WCS
				fp = 0.0
				tpAdj = 0
			!END
			ELSE IF ((redistStatus.EQ.40).AND.(precip.LT.ksat)) THEN
			!BEGIN//merge wetting front 1 and 2 - three wetting fronts redistributing
				redistStatus = 31
				WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
				WF(1)%WC = WF(2)%WC
				WF(1)%WCHold = WF(2)%WCHold
				WF(1)%Z = WF(1)%FAmt/(WF(1)%WC - WCI)
				WF(2)%Z = WF(3)%Z
				WF(2)%FAmt = WF(3)%FAmt
				WF(2)%WC = WF(3)%WC
				WF(2)%WCHold = WF(3)%WCHold
				WF(3)%Z = WF(4)%Z
				WF(3)%FAmt = WF(4)%FAmt
				WF(3)%WC = WF(4)%WC
				WF(3)%WCHold = WF(4)%WCHold
				WF(4)%Z = 0.0
				WF(4)%FAmt = 0.0
				WF(4)%WC = 0.0
				WF(4)%WCHold = 0.0
				fp = 0.0
			ENDIF

			IF ((fp.GT.1E-10).AND.(precip.GE.fp)) THEN
cccccccccccccccc ponding starts immediately ccccccccccccccccccccccccccccccccccccccc
cccccccccccccccc calculate cumulative infiltration at ponding and time shift cccccc
			!BEGIN
				timeToPond = tstart
				IF ((redistStatus.EQ.10).OR.(redistStatus.EQ.11)) THEN
				!BEGIN
					Bfp = ksat*Sav*(WF(1)%WC - WCI)/(fp - ksat)
                      Call calcTpp(WCI,Sav)
				!END
				ELSE IF ((redistStatus.EQ.20).OR.
     +						(redistStatus.EQ.21)) THEN
				!BEGIN
					Bfp = ksat*Sav*(WF(2)%WC - WF(1)%WC)/(fp - ksat)
                      Call calcTpp(WF(1)%WC,Sav)
				!END
				ELSE IF ((redistStatus.EQ.30).OR.
     +						(redistStatus.EQ.31)) THEN
				!BEGIN
					Bfp = ksat*Sav*(WF(3)%WC - WF(2)%WC)/(fp - ksat)
                      Call calcTpp(WF(2)%WC,Sav)
				!END
				ELSE
				!BEGIN
					Bfp = ksat*Sav*(WF(4)%WC - WF(3)%WC)/(fp - ksat)
                      Call calcTpp(WF(3)%WC,Sav)
				ENDIF
			!END //if fp > 0
cccccccccccccccccccc ponding does not necessarily start immediately ccccccccccccccccccccccccc
cccccccccccccccccccc calculate potential cumulative infil, time to pond and time shift cccccc
			ELSE IF (precip.GT.ksat) THEN
			!BEGIN
				IF ((redistStatus.EQ.10).OR.(redistStatus.EQ.11)) THEN
				!BEGIN
					Call cummInfilAtPond(WCI,Sav)
					Call calcTimeToPond
					timeToPond = timeToPond + tstart - 
     +								tpAdj*deltaT
					Call calcTpp(WCI,Sav)
				!END
				ELSE IF ((redistStatus.EQ.20).OR.
     +						(redistStatus.EQ.21)) THEN
				!BEGIN
					Call cummInfilAtPond(WF(1)%WC,Sav)
					Call calcTimeToPond
					timeToPond = timeToPond + tstart - 
     +								tpAdj*deltaT
					Call calcTpp(WF(1)%WC,Sav)
				!END
				ELSE IF ((redistStatus.EQ.30).OR.
     +						(redistStatus.EQ.31)) THEN
				!BEGIN
					Call cummInfilAtPond(WF(2)%WC,Sav)
					Call calcTimeToPond
					timeToPond = timeToPond + tstart - 
     +								tpAdj*deltaT
					Call calcTpp(WF(2)%WC,Sav)
				!END
				ELSE
				!BEGIN
					Call cummInfilAtPond(WF(3)%WC,Sav)
					Call calcTimeToPond
					timeToPond = timeToPond + tstart - 
     +								tpAdj*deltaT
					Call calcTpp(WF(3)%WC,Sav)
				ENDIF
			!END //if R > ksat
cccccccccccccccccccccc no ponding will occur ccccccccccccccccccccccccccccccccccccccccc
			ELSE
			!BEGIN
				timeToPond = 99999999
			ENDIF  !//else
			
			IF (timeToPond.GT.tend) THEN
cccccccccccc no ponding occurs during time step; calc infil using calcNoPond only ccccc
			!BEGIN
				timeToPond = 99999999
				tpp = 99999999
				CALL calcNoPond(tstart,tend,precip)
			!END
			ELSE IF (timeToPond.LT.tend) THEN
			!BEGIN
ccccccc ponding will occur during time step; calc infil using both noPond and Pond cccc
ccccccc infil rate from start until ponding at tp ccccccccccccccccccccccccccccccccccccc
				CALL calcNoPond(tstart, timeToPond, precip)

ccccccc infil rate from tp until END cccccccccccccccccccccccccccccccccccccccccccccccccc
				CALL calcPond(TimeToPond,tend,TimeToPond,precip)
					
              ENDIF
		!END  //if ponding < 0
		ELSE IF (pondingAmt.GE.1E-10) THEN
ccccccc ponding from start cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
		!BEGIN
			timeToPond = tstart
			IF (Bfp.EQ.0.0) THEN !//IF ponding on surface before start of simulation
			!BEGIN
				tpp = tstart
			ENDIF
			IF ((redistStatus.EQ.40).AND.(precip.LT.ksat)) THEN
			!BEGIN//merge wetting front 1 and 2 - three wetting fronts
				redistStatus = 30
				WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
				WF(1)%WC = WF(2)%WC
				WF(1)%WCHold = WF(2)%WCHold
				WF(1)%Z = WF(1)%FAmt/(WF(1)%WC-WCI)
				WF(2)%Z = WF(3)%Z
				WF(2)%FAmt = WF(3)%FAmt
				WF(2)%WC = WF(3)%WC
				WF(2)%WCHold = WF(3)%WCHold
				WF(3)%Z = WF(4)%Z
				WF(3)%FAmt = WF(4)%FAmt
				WF(3)%WC = WF(4)%WC
				WF(3)%WCHold = WF(4)%WCHold
				WF(4)%Z = 0.0
				WF(4)%FAmt = 0.0
				WF(4)%WC = 0.0
				WF(4)%WCHold = 0.0
			ENDIF
			IF (precip.LT.fp) THEN !//chance it will not remain ponded the entire time
	!need to find the time all water will infil, ie when no longer ponded
			!BEGIN
				IF (redistStatus.EQ.10) THEN
				!BEGIN
					Bfp = WF(1)%FAmt
					Call calcTpp(WCI,Sav)
					totalInfil = WF(1)%FAmt + pondingAmt
					Call calcTNewton(tstart, tend, timeToPond, 
     +                          precip,totalInfil,WCI,Sav,t1)
				!END
				ELSE IF (redistStatus.EQ.20) THEN
				!BEGIN
					Bfp = WF(2)%FAmt
					Call calcTpp(WF(1)%WC,Sav)
					totalInfil = WF(2)%FAmt + pondingAmt
					Call calcTNewton(tstart, tend, timeToPond,
     +                              precip, totalInfil,WF(1)%WC,Sav,t1)
				!END
				ELSE IF (redistStatus.EQ.30) THEN
				!BEGIN
					Bfp = WF(3)%FAmt
					Call calcTpp(WF(2)%WC,Sav)
					totalInfil = WF(3)%FAmt + pondingAmt
					Call calcTNewton(tstart, tend, timeToPond,
     +                              precip, totalInfil,WF(2)%WC,Sav,t1)
				!END
				ELSE !//redistStatus = 40
				!BEGIN
					Bfp = WF(4)%FAmt
					Call calcTpp(WF(3)%WC,Sav)
					totalInfil = WF(4)%FAmt + pondingAmt
					Call calcTNewton(tstart, tend, timeToPond, 
     +                              precip, totalInfil,WF(3)%WC,Sav,t1)
				ENDIF
			!END
			ELSE !// there is no chance it will loose all ponding
				t1 = tend + 99999999.0 !//set time so that it will remain ponded
			ENDIF
			
			IF (t1.GT.tend) THEN
cccccccccc will remain ponded entire period; infil using pond conditions only ccccccc
			!BEGIN
				CALL calcPond(tstart,tend,timeToPond,precip)
			!END
			ELSE IF (t1.LT.tend) THEN 
cccccccccc ponding will END at t1; calc infil using pond and noPond conditions cccccc
			!BEGIN
ccccccccc ponded portion from tstart to t1 cccccccccccccccccccccccccccccccccccccccccc
				CALL calcPond(tstart,t1,timeToPond,precip)
                                
				IF ((redistStatus.EQ.10).AND.(precip.LT.ksat)) THEN
					redistStatus = 11
				ELSE IF ((redistStatus.EQ.20).AND.
     +							(precip.LT.ksat)) THEN
					redistStatus = 21
				ELSE IF ((redistStatus.EQ.30).AND.
     +							(precip.LT.ksat)) THEN
					redistStatus = 31 
				ENDIF
!/ccccccc no pond portion from t1 to tEND ccccccccccccccccccccccccccccccccccccccccccc
!				calcNoPond(t1,tend,deltaT,precip)
			ENDIF
		ENDIF ! {IF PondingAmt > 0%0}
	!END  {IF precip > 0 or ponding > 0 }
!/cccccccc redistribution only, no infiltration ccccccccccccccccccccccccccccccccccccc
	ELSE
	!BEGIN  //still need to redistribute soil water
		IF (redistStatus.EQ.40) THEN
		!BEGIN//merge wetting front 1 and 2 - three wetting fronts
			redistStatus = 31
			WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
			WF(1)%WC = WF(2)%WC
			WF(1)%WCHold = WF(2)%WCHold
			WF(1)%Z = WF(1)%FAmt/(WF(1)%WC-WCI)
			WF(2)%Z = WF(3)%Z
			WF(2)%FAmt = WF(3)%FAmt
			WF(2)%WC = WF(3)%WC
			WF(2)%WCHold = WF(3)%WCHold
			WF(3)%Z = WF(4)%Z
			WF(3)%FAmt = WF(4)%FAmt
			WF(3)%WC = WF(4)%WC
			WF(3)%WCHold = WF(4)%WCHold
			WF(4)%Z = 0.0
			WF(4)%FAmt = 0.0
			WF(4)%WC = 0.0
			WF(4)%WCHold = 0.0
		!END
		ELSE IF (redistStatus.EQ.30) THEN
			redistStatus = 31
		ELSE IF (redistStatus.EQ.20) THEN
			redistStatus = 21
		ELSE IF (redistStatus.EQ.10) THEN
			redistStatus = 11
		ENDIF

		nsteps = INT((tend - tstart)/deltaT)
		interTime = tstart
		DO T = 0,(nsteps -1)
		!BEGIN
			interTime = interTime + deltaT
			time = interTime
			cummPrecip = cummPrecip
			cummInfil = cummInfil
			pondingAmt = 0.0
			runoffAmt = runoffAmt
			CALL Redistribution(intertime, deltaT)
			CALL CalcWCO
			CALL CalcRWC
			CALL CalcLayerAvgWC
!			Append(SomeTxtFile);
!			Write(SomeTxtFile,intertime10,''1)
			DO I = 1,4 
			!BEGIN
!				Write(SomeTxtFile,WF(I)%FAmt36,''1)
!				Write(SomeTxtFile,WF(I)%z36,''1)
!				Write(SomeTxtFile,WF(I)%WC36,''1)
			ENDDO
!			Write(SomeTxtFile,WCO36,''1)
!			Write(SomeTxtFile,RWC36,''1)
!			Write(SomeTxtFile,cummInfil36,''1)
			DO I = 1,(NumObsLayers-1)
!				WRite(SomeTxtFile,LayerAvgWC(I)36,''1)
			ENDDO
!			WriteLn(SomeTxtFile,LayerAvgWC(NumObsLayers)36,''1)
		ENDDO
		deltaTerr = tend - interTime
		IF (deltaTerr.GT.1E-10) THEN
		!BEGIN
			interTime = tend
			time = interTime
			cummPrecip = cummPrecip
			cummInfil = cummInfil
			pondingAmt = 0.0
			runoffAmt = runoffAmt
			CALL Redistribution(intertime, deltaTerr)
			CALL CalcWCO
			CALL CalcRWC
			CALL CalcLayerAvgWC
!			Append(SomeTxtFile)
!			Write(SomeTxtFile,intertime10,''1)
			DO I = 1,4 
			!BEGIN
!				Write(SomeTxtFile,WF(I)%FAmt36,''1)
!				Write(SomeTxtFile,WF(I)%z36,''1)
!				Write(SomeTxtFile,WF(I)%WC36,''1)
			ENDDO
!			Write(SomeTxtFile,WCO36,''1)
!			Write(SomeTxtFile,RWC36,''1)
!			Write(SomeTxtFile,cummInfil36,''1)
			DO I = 1,(NumObsLayers-1) 
!				WRite(SomeTxtFile,LayerAvgWC(I)36,''1)
			ENDDO
!			WriteLn(SomeTxtFile,LayerAvgWC(NumObsLayers)36,''1)
		ENDIF
	ENDIF
	tpAdj = tpAdj + 1
	End

cccccccccc Subroutines for Green-Ampt calculations ccccccccccccccccccccccccccccccccccccc

ccccccccc	calculates potential cumulative infiltration ccccccccccccccccccccccccccccccccc
	Subroutine CummInfilAtPond(WC,S)
	USE ModuleDefs
	Implicit None
	Double Precision rk, S, WC 

	rk = (precipRate)/(ksat)
	BFp = S *(WCS - WC)/(rk - 1)
	END

cccccccccc calculates time to ponding ccccccccccccccccccccccccccccccccccccccccccccccccccc
	Subroutine CalcTimeToPond
	USE ModuleDefs
  
	TimeToPond = BFp/precipRate
	END

cccccccccc calculates time shift cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	Subroutine CalcTpp(WC,S)
	USE ModuleDefs
	Double Precision WC, S

	Tpp = (BFp - S*(WCS - WC) * log(1 + (BFp/
     +                (S*(WCS-WC)))))/ksat
	END

cccccccc calculates time it will take for remaining ponded water to infiltrate cccccccccc
	SUBROUTINE CalcTNewton(tstart, tend, tp, precip,
     +						totalInfil,WC,S,t1)
      USE ModuleDefs
	Implicit None                  
	Double Precision tstart, tend,WC,S, totalInfil
	Double Precision t1, tp, tshift, precip
	Double Precision tolerance 
	Double Precision t1new 
	Double Precision t1old
	Double Precision F
	Double Precision lnarg
	Double Precision g 
	Double Precision dgdeltaT1
	Double Precision error 
	
	t1 = 0.0
	tolerance = 1.0E-6
	error = 1.0
	t1new = tend
	t1old = t1new
	F = (t1old - tstart)* precip + totalInfil
cccccccccccc the time that occurs is t1 cccccccccccccccccccccccccccccccccccccccccc
	lnarg = 1.0 + (F/(S*(WCS - WC)));
	g = F - (S*(WCS - WC))*(log(lnarg)) - ksat*(t1old - tp + tpp)
	DO WHILE (error.GT.tolerance)
		dgdeltaT1 = precip - precip*(1.0/lnarg) - ksat
		t1new = t1old - (g/dgdeltaT1)
		t1old = t1new
		F = (t1old - tstart)* precip + totalInfil
		lnarg = 1.0 + (F/(S*(WCS - WC)))
		g = F - (S*(WCS - WC))*(log(lnarg)) - ksat*(t1old - tp + tpp)
		error = abs(g)
		t1 = t1new
	ENDDO
	
	END

cccccccccccccc calculates cumulative infiltration using Newton Raphson method cccccccccccc
cccccccccc	Adapted from:  Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1989.
cccccccccc	Numerical Recipes in Pascal, Cambridge University Press: Port Chester, New York.}
	Subroutine CalcFNewton(Fnew, t, tp, WC,S,F)
	USE ModuleDefs
	Implicit None
	Double Precision Fnew, t, tp, WC, S,F
	Double Precision  g 
	Double Precision  dgdf 
	Double Precision  Fold 
	Double Precision  error 
	Double Precision  tolerance 
	
	  F = 0.0
	  tolerance = 1.0E-6
	  error = 1.0
	  Fold = Fnew + 1.0E-6
	  g = Fold - (S*(WCS - WC)) * log(1.0 + (Fold/(S*(WCS - WC))))-
     +					ksat*(t - tp + tpp)
	  DO WHILE (error.GE.tolerance)
	  !BEGIN
	    dgdf = 1.0 - ((S*(WCS - WC))/((S*(WCS - WC)) + Fold))
	    Fnew = Fold - (g/dgdf)
	    Fold = Fnew
	    g = Fold - S*(WCS - WC) * log(1.0 + (Fold/(S*(WCS - WC))))-
     +				ksat*(t - tp + tpp)
	    error = abs(g)
	    F = Fnew
	    ENDDO
	END

ccccccccc calculates approximate cumulative infiltration solution proposed by
ccccccccc	Srivastava, Costello, Edwards. 1996. A Direct, Approximate Solution to the
ccccccccc	Modified Green-Ampt Infiltration Equation.  Transactions of ASAE. 39(4): 1411-1413
	Subroutine CalcFApprox(Fnew, t, tp, WC,S,F)
	USE ModuleDefs
	Implicit None
	Double Precision Fnew, t, tp, WC,S,F
	Double Precision  phi  
	Double Precision  alpha  
	Double Precision  beta  
	Double Precision  delta 
	Double Precision  lambda2  

	
	  F = 0.0
	  alpha = 1.851
	  beta = 0.565
	  delta = 0.004
	  phi = ksat*(t - tp + tpp)/(S*(WCS - WC))
	    IF (phi.GE.0.0001) THEN
		!BEGIN
			IF (phi.LE.0.095) THEN
			!BEGIN
			  alpha = 1.851
			  beta = 0.565
			  delta = 0.004
			!END
			ELSE IF (phi.LE.0.911) THEN
			!BEGIN
			  alpha = 2.137
			  beta = 0.667
			  delta = 0.021
			!END
			ELSE IF (phi.LE.17.00) THEN
			!BEGIN
			  alpha = 2.141
			  beta = 0.689
			  delta = 0.035
			!END
			ELSE
			   PRINT*, 'ERROR: OUT OF RANGE FOR FAPPROX CALC'  
			ENDIF
	    ENDIF
	  lambda2 = (phi**(beta + delta*LOG(phi)))
	  F = alpha*(S*(WCS - WC))*lambda
	END



cccccc calculates infiltration rate, cumulative infiltration, runoff and surface storage
cccccc under nonponded conditions; calls redistribution process for lower layers 
	Subroutine CalcNoPond(tstart, tend, R)
	USE ModuleDefs
	Implicit None

	Double Precision tstart, tend, R
	Integer nsteps 
	Double Precision deltaTerr 
	Double Precision interTime 
	Integer kk, i , j
	Double Precision expon 
	Double Precision deltaInfil  
	Double Precision y 
	Integer nbad, nok 

	infilRate = 0.0

	nsteps = INT((tend - tstart)/deltaT)
	interTime = tstart
	DO kk = 0,(nsteps - 1) 
	!BEGIN
		interTime = interTime + deltaT
		time = interTime
		cummPrecip = cummPrecip + deltaT*R
		cummInfil = cummInfil + R*deltaT
		pondingAmt = 0.0
		runoffAmt = runoffAmt
		infilRate = R
		IF (redistStatus.EQ.40) THEN
		!BEGIN
			WF(4)%FAmt = cummInfil - WF(3)%FAmt -
     +                             WF(2)%FAmt - WF(1)%FAmt
			fp =  Ksat + (Ksat*Sav*(WCS - WF(3)%WC)/WF(4)%FAmt)
                              
		!END
		ELSE IF ((redistStatus.EQ.30).OR.(redistStatus.EQ.31)) THEN
		!BEGIN
			WF(3)%FAmt = cummInfil - WF(2)%FAmt - WF(1)%FAmt
           	fp =  Ksat + (Ksat*Sav*(WCS - WF(2)%WC)/WF(3)%FAmt)
		!END
		ELSE IF ((redistStatus.EQ.20).OR.(redistStatus.EQ.21)) THEN
		!BEGIN
			WF(2)%FAmt = cummInfil - WF(1)%FAmt
			fp =  Ksat + (Ksat*Sav*(WCS - WF(1)%WC)/WF(2)%FAmt)
          !END
		ELSE
		!BEGIN
			WF(1)%FAmt= cummInfil
			fp =  Ksat + (Ksat*Sav*(WCS -WCI)/WF(1)%FAmt)
          ENDIF

cccccccccccccccccccccc redistribute water ccccccccccccccccccccccccccccccccccccccc
ccccccccccccccc calculate generic gamma coefficient ccccccccccccccccccccccccccccc
		DO J = 1,3 
			Call adjFactorCalc(wf(j)%numRedist,
     +						wf(j)%redistTime,wf(j)%adjFactor)                       
		ENDDO
ccccccccccccccc preponding, estimate rising water content cccccccccccccccccccccc
		IF ((pondFlag.EQ.0).AND.(R.GT.0.0).AND.(WF(2)%WC.EQ.0.0)) THEN
		!BEGIN
			expon = 1/(3 + (2/lambda))
			WF(1)%WC = WCI + (WCS-WCI)*((R/fp)**expon)
			WF(1)%WCHold = WF(1)%WC
			Call calcRedistZf(cummInfil,WF(1)%WC,WCI,WF(1)%z)
                                  
		!END !//end pondflag
		ELSE IF (redistStatus.EQ.11) THEN
		!BEGIN
			IF (pondingAmt.EQ.0.0) THEN
			!BEGIN
				redistStatus = 10 !// do not redistribute wf 1
				Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
              !END
			ELSE
				Call Redist11(R,interTime-deltaT,intertime,deltaT)
			ENDIF
		!END
		ELSE IF (redistStatus.EQ.20) THEN
		!BEGIN
			Call Redist20(R,interTime-deltaT,intertime,deltaT)
		!END
		ELSE IF (redistStatus.EQ.21) THEN
		!BEGIN
			IF (pondingAmt.EQ.0.0) THEN
			!BEGIN
				redistStatus = 20
				Call Redist20(R,interTime-deltaT,intertime,deltaT)
			!END
			ELSE
				Call Redist21(R,interTime-deltaT,intertime,deltaT)
			ENDIF
		!END
		ELSE IF (redistStatus.EQ.30) THEN
		!BEGIN
			Call Redist30(R,interTime-deltaT,intertime,deltaT)
		!END
		ELSE IF (redistStatus.EQ.31) THEN
		!BEGIN
			IF (pondingAmt.EQ.0.0) THEN
			!BEGIN
				redistStatus = 30
				Call Redist30(R,interTime-deltaT,intertime,deltaT)
			!END
			ELSE
				Call Redist31(R,interTime-deltaT,intertime,deltaT)
			ENDIF
		!END
		ELSE IF (redistStatus .EQ.40) THEN
		!BEGIN
			Call Redist40(R,interTime-deltaT,intertime,deltaT)
		!END
		ELSE !//redistStatus = 10, do not redistribute , only move wf1 down
			 Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
          ENDIF                         

		IF (WF(1)%WC.GT.WCS) THEN
		!BEGIN
			WF(1)%WC = WCS
			WF(1)%WCHold = WCS
			Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
          ENDIF  !//end wc1 > wsc

		CALL CalcWCO
		Call CalcRWC
		CALL CalcLayerAvgWC
!		Append(SomeTxtFile);
!		Write(SomeTxtFile,intertime10,''1);
		DO I = 1,4 
		!BEGIN
!			Write(SomeTxtFile,WF(I)%FAmt36,''1);
!			Write(SomeTxtFile,WF(I)%z36,''1);
!			Write(SomeTxtFile,WF(I)%WC36,''1);
		ENDDO
!		Write(SomeTxtFile,WCO36,''1);
!		Write(SomeTxtFile,RWC36,''1);
!		Write(SomeTxtFile,cummInfil36,''1);

		DO I = 1,(NumObsLayers-1)
!			WRite(SomeTxtFile,LayerAvgWC(I)36,''1);
!		WriteLn(SomeTxtFile,LayerAvgWC(NumObsLayers)36,''1);
		ENDDO
	ENDDO !//for kk

ccccccccccccccc  Check that really at END of period cccccccccccccccccccccccccc
	deltaTerr = tend - interTime
	IF (deltaTerr.GT.1E-10) THEN
	!BEGIN
		interTime = tend
		time = interTime
		cummPrecip = cummPrecip + deltaTerr*R
		cummInfil = cummInfil + R*deltaTerr
		pondingAmt = 0.0
		runoffAmt = runoffAmt
		infilRate = R
		IF (redistStatus.EQ.40) THEN
		!BEGIN
			WF(4)%FAmt = cummInfil - WF(3)%FAmt -
     +                            WF(2)%FAmt - WF(1)%FAmt
			fp =  Ksat + (Ksat*Sav*(WCS -WF(3)%WC)/WF(4)%FAmt)
                             
		!END
		ELSE IF ((redistStatus.EQ.30).OR.(redistStatus.EQ.31)) THEN
		!BEGIN
			WF(3)%FAmt = cummInfil - WF(2)%FAmt - WF(1)%FAmt
              fp =  Ksat + (Ksat*Sav*(WCS - WF(2)%WC)/WF(3)%FAmt)
		!END
		ELSE IF ((redistStatus.EQ.20).OR.(redistStatus.EQ.21)) THEN
		!BEGIN
			WF(2)%FAmt = cummInfil - WF(1)%FAmt
			fp =  Ksat + (Ksat*Sav*(WCS -WF(1)%WC)/WF(2)%FAmt)
		!END
		ELSE
		!BEGIN
			WF(1)%FAmt = cummInfil
			fp =  Ksat + (Ksat*Sav*(WCS - WCI)/WF(1)%FAmt)                                  
		ENDIF
cccccccccccccccccc redistribute water ccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccc calc generic gamma coefficient cccccccccccccccccccccccccccccc
			DO J = 1, 3 
				 Call adjFactorCalc(wf(j)%numRedist,(wf(j)%redistTime
     +                              -deltaT+deltaTerr),wf(j)%adjFactor)
			ENDDO
cccccccccccccc   preponding, estimate rising water content ccccccccccccccccccccc
		IF ((pondFlag.EQ.0).AND.(R.GT.0.0).AND.(WF(2)%WC.EQ.0.0)) THEN
		!BEGIN
			expon = 1/(3 + (2/lambda))
			WF(1)%WC = WCI + (WCS-WCI)*((R/fp)**expon)
          	WF(1)%WCHold = WF(1)%WC
			Call calcRedistZf(cummInfil,WF(1)%WC,WCI,WF(1)%z)
         	!END //end pondflag
		ELSE IF (redistStatus.EQ.11) THEN
		!BEGIN
			IF (pondingAmt.GT.0.0) THEN
			!BEGIN
				redistStatus = 10 ! // do not redistribute wf 1
				Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
              !END
			ELSE
				Call Redist11(R,interTime-deltaTerr,intertime,
     +									deltaTerr)
			ENDIF
		!END
		ELSE IF (redistStatus.EQ.20) THEN
		!BEGIN
			Call Redist20(R,interTime-deltaTerr,intertime,deltaTerr)
		!END
		ELSE IF (redistStatus.EQ.21) THEN
		!BEGIN
			IF (pondingAmt.GT.0.0) THEN
			!BEGIN
				redistStatus = 20
				Call Redist20(R,interTime-deltaTerr,intertime,
     +								deltaTerr)
			!END
			ELSE
				Call Redist21(R,interTime-deltaTerr,intertime,
     +								deltaTerr)
			ENDIF
		!END
		ELSE IF (redistStatus.EQ.30) THEN
		!BEGIN
			Call Redist30(R,interTime-deltaTerr,intertime,
     +								deltaTerr)
		!END
		ELSE IF (redistStatus.EQ.31) THEN
		!BEGIN
			IF (pondingAmt.GT.0.0) THEN
			!BEGIN
				redistStatus = 30
				Call Redist30(R,interTime-deltaTerr,intertime,
     +								deltaTerr)
			!END
			ELSE
				Call Redist31(R,interTime-deltaTerr,intertime,
     +								deltaTerr)
			ENDIF
		!END
		ELSE IF (redistStatus.EQ.40) THEN
		!BEGIN
			Call Redist40(R,interTime-deltaTerr,intertime,deltaTerr)
		!END
		ELSE !//redistStatus = 10, do not redistribute , only move wf1 down
			Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
		ENDIF
       
		IF (WF(1)%WC.GT.WCS) THEN
		!BEGIN
			WF(1)%WC = WCS
			WF(1)%WCHold = WCS
			Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
          ENDIF  !//end wc1 > wsc
	ENDIF !//end if dterr

	END !//calcNoPond

ccccccc   calculates infiltration rate, cumulative infiltration, runoff and surface storage ccccccccc
ccccccc   under ponded conditions; calls redistribution process for lower layers cccccccccccccccccccc
	Subroutine CalcPond(tstart, tend, tp, R)
	USE ModuleDefs
	Implicit None
	Double Precision tstart, tend,tp,R
	Double Precision  amtWater 
	Double Precision  Fnew 
	Double Precision  Fhold
	Double Precision  t
	Integer  nsteps 
	Integer  kk, i ,j
	Double Precision  interTime 
	Double Precision  deltaTerr 
	Double Precision  deltaCummInfil 
	Double Precision  y 
	
	
	  infilRate = 0.0
	  IF (pondFlag.EQ.0.AND.redistStatus.EQ.10) THEN !//used when preponding rise
	  !BEGIN
	    WF(1)%WCHold = WCS
	    WF(1)%WC = WCS
	  ENDIF
	  pondFlag = 1
	  nsteps = INT((tend - tstart)/deltaT)
	  interTime = tstart
	  DO kk = 0,(nsteps -1)
	  !BEGIN
	    interTime = interTime + deltaT
	    time = interTime
	    cummPrecip = cummPrecip + deltaT*R
	    amtWater = pondingAmt + R*deltaT
	    Fhold = cummInfil
	    t = interTime
	    IF (redistStatus.EQ.10.OR.redistStatus.EQ.11) THEN
	    !BEGIN
			Fnew = cummInfil + amtWater
ccccccccccccccccccccccc  Use calcFNewton function for regular GA cccccccccccccccccccccccccccccccc
ccccccccccccccccccccccc  Use calcFApprox function for approximate GA cccccccccccccccccccccccccccc
			CALL calcFNewton(Fnew,t,tp,WCI,Sav,cummInfil)
!			CALL calcFApprox(Fnew,t,tp,WCI,Sav,cummInfil)
			infilRate =  ksat + ((ksat * Sav*(WF(1)%WC -
     +					  WCI))/cummInfil)
			WF(1)%FAmt = cummInfil
			deltaCummInfil = cummInfil - Fhold
	    !END// end if redist 10 or 11
	    ELSE IF (redistStatus.EQ.20.OR.redistStatus.EQ.21) THEN
	    !BEGIN
cccccccccccccccccccccc  update values since WC1 is changing everytime step cccccccccccccccccccccc
			Bfp = WF(2)%FAmt
			tp = t - deltaT
			CALL calcTpp(WF(1)%WC,Sav)
			tpp = tpp
			Fnew = WF(2)%FAmt + amtWater
cccccccccccccccccccccc Use calcFNewton function for regular GA cccccccccccccccccccccccccccccccccc
cccccccccccccccccccccc Use calcFApprox function for approximate GA cccccccccccccccccccccccccccccc
			CALL calcFNewton(Fnew,t,tp,WF(1)%WC,Sav,WF(2)%FAmt)
!			Call calcFApprox(Fnew,t,tp,WF(1%WC,Sav,WF(2)%FAmt)
			infilRate =  ksat + ((ksat * Sav*(WCS-
     +						WF(1)%WC))/WF(2)%FAmt)
			cummInfil = WF(1)%FAmt + WF(2)%FAmt
			deltaCummInfil = cummInfil - Fhold
	    !END // end redist = 20 or 21
	    ELSE IF (redistStatus.EQ.30.OR.redistStatus.EQ.31) THEN
	    !BEGIN
cccccccccccccccccccccccc  update values since WC1 is changing everytime step cccccccccccccccccc
			Bfp = WF(3)%FAmt
			tp = t - deltaT
			CALL calcTpp(WF(2)%WC,Sav)
			tpp = tpp
			Fnew = WF(3)%FAmt + amtWater
cccccccccccccccccccccccc Use calcFNewton function for regular GA cccccccccccccccccccccccccccccc
cccccccccccccccccccccccc Use calcFApprox function for approximate GA cccccccccccccccccccccccccccc
			Call calcFNewton(Fnew,t,tp,WF(2)%WC,Sav,WF(3)%FAmt)
!			Call calcFApprox(Fnew,t,tp,WF(2)%WC,Sav,WF(3)%FAmt)
			infilRate =  ksat + ((ksat * Sav*(WCS-
     +						  WF(2)%WC))/WF(3)%FAmt)
			cummInfil = WF(1)%FAmt + WF(2)%FAmt +
     +								  WF(3)%FAmt
			deltaCummInfil = cummInfil - Fhold
	    !END // end redist = 30 or 31
	    ELSE !// redist = 40
	    !BEGIN
cccccccccccccccccccccccc  update values since WC1 is changing everytime step cccccccccccccccccccc
			Bfp = WF(4)%FAmt
			tp = t - deltaT
			CALL calcTpp(WF(3)%WC,Sav)
			tpp = tpp
			Fnew = WF(4)%FAmt + amtWater
cccccccccccccccccccccccc Use calcFNewton function for regular GA cccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccc Use calcFApprox function for approximate GA cccccccccccccccccccccccccccccc
			Call calcFNewton(Fnew,t,tp,WF(3)%WC,Sav,WF(4)%FAmt)
!			Call calcFApprox(Fnew,t,tp,WF(3)%WC,Sav,WF(4)%FAmt)
			infilRate =  ksat + ((ksat * Sav*(WCS-
     +							  WF(3)%WC))/WF(4)%FAmt)
			cummInfil = WF(1)%FAmt + WF(2)%FAmt +
     +								  WF(3)%FAmt + WF(4)%FAmt
			deltaCummInfil = cummInfil - Fhold
	    ENDIF  !// end else
	    fp = infilRate
	    IF (amtWater.GT.deltaCummInfil) THEN
	    !BEGIN
			pondingAmt = amtWater - deltaCummInfil
			IF (pondingAmt.GE.smax) THEN
			!BEGIN
				runoffAmt = runoffAmt + (pondingAmt - smax)
				pondingAmt = smax
			!END // end if pond > smax
			ELSE
				runoffAmt = runoffAmt
			ENDIF
	    !END //end if water > deltacumm
	    ELSE
	    !BEGIN
			runoffAmt = runoffAmt
			pondingAmt = 0.0
	    ENDIF  !// end else

cccccccccccccccccccccccccccccccc redistribute water ccccccccccccccccccccccccccccccccccccc
	    !calculate calculate generic gamma coefficient
	    DO J = 1,3 
			 CALL adjFactorCalc(wf(j)%numRedist,wf(j)%redistTime,
     +									wf(j)%adjFactor)
		ENDDO

	    IF (redistStatus.EQ.11) THEN
	    !BEGIN
			IF (pondingAmt.GT.0.0) THEN
			!BEGIN
				redistStatus = 10 !// do not redistribute wf 1
				CALL calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
									  
			!END
			ELSE
				CALL Redist11(R,interTime-deltaT,intertime,deltaT)
			ENDIF
	    !END
	    ELSE IF (redistStatus.EQ.20) THEN
	    !BEGIN
			CALL Redist20(R,interTime-deltaT,intertime,deltaT)
	    !END
	    ELSE IF (redistStatus.EQ.21) THEN
	    !BEGIN
			IF (pondingAmt.GT.0.0) THEN
			!BEGIN
				redistStatus = 20
				CALL Redist20(R,interTime-deltaT,intertime,deltaT)
			!END
			ELSE
				CALL Redist21(R,interTime-deltaT,intertime,deltaT)
			ENDIF
	    !END
	    ELSE IF (redistStatus.EQ.30) THEN
	    !BEGIN
			CALL Redist30(R,interTime-deltaT,intertime,deltaT)
	    !END
	    ELSE IF (redistStatus.EQ.31) THEN
	    !BEGIN
			IF (pondingAmt.GE.0.0) THEN
			!BEGIN
				redistStatus = 30
				CALL Redist30(R,interTime-deltaT,intertime,deltaT)
			!END
			ELSE
				CALL Redist31(R,interTime-deltaT,intertime,deltaT)
			ENDIF
	    !END
	    ELSE IF (redistStatus.EQ.40) THEN
	    !BEGIN
			CALL Redist40(R,interTime-deltaT,intertime,deltaT)
	    !END
	    ELSE !//redistStatus = 10, do not redistribute , only move wf1 down
			CALL calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
     							  
		ENDIF

	    IF (WF(1)%WC.GE.WCS) THEN
	    !BEGIN
			WF(1)%WC = WCS
			WF(1)%WCHold = WCS
			CALL calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
								
	    ENDIF   !//end wc1 > wsc

	    CALL CalcWCO
	    CALL CalcRWC
	    CALL CalcLayerAvgWC
!	    Append(SomeTxtFile)
!	    Write(SomeTxtFile,intertime10,''1)
	    DO I = 1,4 
	    !BEGIN
!			Write(SomeTxtFile,WF(I)%FAmt36,''1)
!			Write(SomeTxtFile,WF(I)%z36,''1)
!			Write(SomeTxtFile,WF(I)%WC36,''1)
	    ENDDO
!	    Write(SomeTxtFile,WCO36,''1)
!	    Write(SomeTxtFile,RWC36,''1)
!	    Write(SomeTxtFile,cummInfil36,''1)

	    DO I = 1,(NumObsLayers-1) 
!			Write(SomeTxtFile,LayerAvgWC(I)36,''1)
		ENDDO
!	    WriteLn(SomeTxtFile,LayerAvgWC(numObsLayers)36,''1)
	  ENDDO !// for kk

ccccccccccccccccccccc Check that at END of period cccccccccccccccccccccccccccccccccccccccc
	  deltaTerr = tend - interTime
	  IF (deltaTerr.GE.1E-10) THEN
	  !BEGIN
	    interTime = tend
	    time = interTime
	    cummPrecip = cummPrecip + deltaTerr*R
	    amtWater = pondingAmt + R*deltaTerr
	    Fhold = cummInfil
	    t = interTime
	    IF (redistStatus.EQ.10.OR.redistStatus.EQ.11) THEN
	    !BEGIN
			Fnew = cummInfil + amtWater
cccccccccccccccccccccc Use calcFNewton function for regular GA cccccccccccccccccccccccccc
cccccccccccccccccccccc Use calcFApprox function for approximate GA ccccccccccccccccccccccc
			Call calcFNewton(Fnew,t,tp,WCI,Sav,cummInfil)
!			Call calcFApprox(Fnew,t,tp,WCI,Sav,cummInfil)
			infilRate =  ksat + ((ksat * Sav*(WF(1)%WC -
     +						WCI))/cummInfil)
			WF(1)%FAmt = cummInfil
			deltaCummInfil = cummInfil - Fhold
	    !END// end if redist 10 or 11
	    ELSE IF (redistStatus.EQ.20.OR.redistStatus.EQ.21) THEN
	    !BEGIN
ccccccccccccccccccccccccc   update values since WC1 is changing everytime step ccccccccccccccc
			Bfp = WF(2)%FAmt
			tp = t - deltaTerr
			CALL calcTpp(WF(1)%WC,Sav)
			tpp = tpp
			Fnew = WF(2)%FAmt + amtWater
cccccccccccccccccccccc Use calcFNewton function for regular GA cccccccccccccccccccccccccc
cccccccccccccccccccccc Use calcFApprox function for approximate GA ccccccccccccccccccccccc
			Call calcFNewton(Fnew,t,tp,WF(1)%WC,Sav,WF(2)%FAmt)
!			Call calcFApprox(Fnew,t,tp,WF(1)%WC,Sav,WF(2)%FAmt)
			infilRate =  ksat + ((ksat * Sav*(WCS-
     +						WF(1)%WC))/WF(2)%FAmt)
			cummInfil = WF(1)%FAmt + WF(2)%FAmt
			deltaCummInfil = cummInfil - Fhold
	    !END // end redist = 20 or 21
	    ELSE IF (redistStatus.EQ.30.OR.redistStatus.EQ.31) THEN
	    !BEGIN
ccccccccccccccccccccccc  update values since WC1 is changing everytime step ccccccccccccccccccc
			Bfp = WF(3)%FAmt
			tp = t - deltaTerr
			CALL calcTpp(WF(2)%WC,Sav)
			tpp = tpp
			Fnew = WF(3)%FAmt + amtWater
cccccccccccccccccccccc Use calcFNewton function for regular GA cccccccccccccccccccccccccc
cccccccccccccccccccccc Use calcFApprox function for approximate GA ccccccccccccccccccccccc
			Call calcFNewton(Fnew,t,tp,WF(2)%WC,Sav,WF(3)%FAmt)
!			Call calcFApprox(Fnew,t,tp,WF(2)%WC,Sav,WF(3)%FAmt)
			infilRate =  ksat + ((ksat * Sav*(WCS-
     +						  WF(2)%WC))/WF(3)%FAmt)
			cummInfil = WF(1)%FAmt + WF(2)%FAmt +
     +								WF(3)%FAmt
			deltaCummInfil = cummInfil - Fhold
	    !END // end redist = 30 or 31
	    ELSE !// redist = 40
	    !BEGIN
ccccccccccccccccccccccc  update values since WC1 is changing everytime step ccccccccccccccccccc
			Bfp = WF(4)%FAmt
			tp = t - deltaTerr
			CALL calcTpp(WF(3)%WC,Sav)
			tpp = tpp
			Fnew = WF(4)%FAmt + amtWater
cccccccccccccccccccccc Use calcFNewton function for regular GA cccccccccccccccccccccccccc
cccccccccccccccccccccc Use calcFApprox function for approximate GA ccccccccccccccccccccccc
			Call calcFNewton(Fnew,t,tp,WF(3)%WC,Sav,WF(4)%FAmt)
!			Call calcFApprox(Fnew,t,tp,WF(3)%WC,Sav,WF(4)%FAmt)
			infilRate =  ksat + ((ksat * Sav*(WCS-
     +						WF(3)%WC))/WF(4)%FAmt)
			cummInfil = WF(1)%FAmt + WF(2)%FAmt +
     +								WF(3)%FAmt + WF(4)%FAmt
			deltaCummInfil = cummInfil - Fhold
	    ENDIF !// end else

	    fp = infilRate

	    IF (amtWater.GT.deltaCummInfil) THEN
	    !BEGIN
			pondingAmt = amtWater - deltaCummInfil
			IF (pondingAmt.GT.smax )THEN
			!BEGIN
				runoffAmt = runoffAmt + (pondingAmt - smax)
				pondingAmt = smax
			!END // end if pond > smax
			ELSE
				runoffAmt = runoffAmt
			ENDIF
	    !END //end if water > deltacumm
	    ELSE
	    !BEGIN
			runoffAmt = runoffAmt
			pondingAmt = 0.0
	    ENDIF !// end else

ccccccccccccccccccccccccccccccc redistribute water ccccccccccccccccccccccccccccccccccccc
	    !calc generic gamma coefficient need to adjust for intermediate time step
	    DO J = 1,3 
			 CALL adjFactorCalc(wf(j)%numRedist,(wf(j)%redistTime-
     +								deltaT+deltaTerr),wf(j)%adjFactor)
		ENDDO

	    IF (redistStatus.EQ.11) THEN
	    !BEGIN
			IF (pondingAmt.GE.0.0) THEN
			!BEGIN
				redistStatus = 10 !// do not redistribute wf 1
				CALL calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
											
			!END
			ELSE
				CALL Redist11(R,interTime-deltaTerr,intertime,
     +									deltaTerr)
			ENDIF
	    !END
	    ELSE IF (redistStatus.EQ.20) THEN
	    !BEGIN
			CALL Redist20(R,interTime-deltaTerr,intertime,deltaTerr)
	    !END
	    ELSE IF (redistStatus.EQ.21) THEN
	    !BEGIN
			IF (pondingAmt.GE.0.0) THEN
			!BEGIN
				redistStatus = 20
				CALL Redist20(R,interTime-deltaTerr,intertime,
     +									deltaTerr)
			!END
			ELSE
				CALL Redist21(R,interTime-deltaTerr,intertime,
     +									deltaTerr)
			ENDIF
	    !END
	    ELSE IF (redistStatus.EQ.30) THEN
	    !BEGIN
			CALL Redist30(R,interTime-deltaTerr,intertime,deltaTerr)
	    !END
	    ELSE IF (redistStatus.EQ.31) THEN
	    !BEGIN
			IF (pondingAmt.GE.0.0) THEN
			!BEGIN
				redistStatus = 30
				CALL Redist30(R,interTime-deltaTerr,intertime,
     +									deltaTerr)
			!END
			ELSE
				CALL Redist31(R,interTime-deltaTerr,intertime,
     +									deltaTerr)
			ENDIF
	    !END
	    ELSE IF (redistStatus.EQ.40) THEN
	    !BEGIN
			CALL Redist40(R,interTime-deltaTerr,intertime,deltaTerr)
	    !END
	    ELSE !//redistStatus = 10, do not redistribute , only move wf1 down
			CALL calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
								  
		ENDIF

	    IF (WF(1)%WC.GE.WCS) THEN
	    !BEGIN
			WF(1)%WC = WCS
			WF(1)%WCHold = WCS
			CALL calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
								  
	    ENDIF  !/end wc1 > wsc

	    CALL CalcWCO
	    CALL CalcRWC
	    CALL CalcLayerAvgWC
!	    Append(SomeTxtFile)
!	    Write(SomeTxtFile,intertime10,''1)
	    DO I = 1,4 
	    !BEGIN
!			Write(SomeTxtFile,WF(I)%FAmt36,''1)
!			Write(SomeTxtFile,WF(I)%z36,''1)
!			Write(SomeTxtFile,WF(I)%WC36,''1)
	    ENDDO
!	    Write(SomeTxtFile,WCO36,''1)
!	    Write(SomeTxtFile,RWC36,''1)
!	    Write(SomeTxtFile,cummInfil36,''1)

	    DO I = 1,(NumObsLayers-1) 
!			WRite(SomeTxtFile,LayerAvgWC(I)36,''1)
		ENDDO
!	    WriteLn(SomeTxtFile,LayerAvgWC(NumObsLayers)36,''1)
	  ENDIF  !//if deltaterr
	END





ccccccccccccccccc redistribution calculations cccccccccccccccccccccccccccccccccccccccc

cccccc calculates depth based on cumulative infiltration and water content cccccccccccc
	SUBROUTINE CalcRedistZf(F, WCA, WCB,Z) !//calcs z	
	USE ModuleDefs
	Implicit None
	Double Precision F, WCA, WCB,Z
	
	IF (WCA.EQ.WCI) THEN
      
		Z = F/(WCS - WCI)
      ELSE
		Z = F/(WCA - WCB)
      ENDIF
	END

ccccccccccc calculates SURFACE relative water content cccccccccccccccccccccccccccccccccccccccccccc
	Subroutine CalcRWC
	USE ModuleDefs
	Implicit None
  
	IF ((redistStatus.EQ.10).OR.(redistStatus.EQ.11)) THEN
		RWC = (WF(1)%WC - WCR)/(WCS - WCR)
	ELSE IF ((redistStatus.EQ.20).OR.(redistStatus.EQ.21)) THEN
		RWC = (WF(2)%WC - WCR)/(WCS - WCR)
	ELSE IF ((redistStatus.EQ.30).OR.(redistStatus.EQ.31)) THEN
		RWC = (WF(3)%WC - WCR)/(WCS - WCR)
	ELSE
		RWC = (WF(4)%WC - WCR)/(WCS - WCR)
	ENDIF
	
	END

ccccccc   calculates surface water content based on redistribution status ccccccccccccccccccccccccccccc
	Subroutine CalcWCO
	USE ModuleDefs
		
	IF (redistStatus.EQ.10.OR.redistStatus.EQ.11) THEN
	    WCO = WF(1)%WC
	ELSE IF (redistStatus.EQ.20.OR.redistStatus.EQ.21) THEN
	    WCO = WF(2)%WC
	ELSE IF (redistStatus.EQ.30.OR.redistStatus.EQ.31) THEN
	    WCO =  WF(3)%WC
	ELSE
	    WCO =  WF(4)%WC
	ENDIF
	  
	END

cccccccccc    calculates average water conent of observation layers ccccccccccccccccccccccccccccccccccc
	Subroutine CalcLayerAvgWC
	USE ModuleDefs
	Implicit None
	
	Integer L 
	
	  Do L = 1,NumObsLayers 
	  !BEGIN
	    IF (redistStatus.EQ.10.OR.redistStatus.EQ.11) THEN
	    !BEGIN
			IF (WF(1)%z.GE.LayerBottomZ(L)) THEN
			!BEGIN
				LayerAvgWC(L) = WF(1)%WC 
			!END
			ELSE IF ( WF(1)%z.LE.(LayerBottomZ(L)-
     +						ProfileDepth(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WCI 
			!END
			ELSE IF ((WF(1).z.GE.LayerBottomZ(L) - ProfileDepth(L)).
     +						AND.(WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(1)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(1)%WC +
     +									(LayerBottomZ(L)- WF(1)%z)*
     +									WCI)/ProfileDepth(L)
			ENDIF
	    !END  //end if redist = 10 or 11
	    ELSE IF (redistStatus.EQ.20.OR.redistStatus.EQ.21) THEN
	    !BEGIN
			IF (WF(2)%z.GE.LayerBottomZ(L)) THEN
			!BEGIN
				LayerAvgWC(L) = WF(2)%WC 
			!END
			ELSE IF (WF(1)%z.LE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WCI 
			!END
			ELSE IF ((WF(2)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					(WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(2)%z - LayerBottomZ(L) +
     +									ProfileDepth(L)) * WF(2)%WC +
     +									(LayerBottomZ(L) - WF(2)%z)
     +	 								* WF(1)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(2)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					(WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(2)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(2)%WC +
     +									(WF(1)%z-WF(2)%z)*
     +									WF(1)%WC+(LayerBottomZ(L) -
     +									WF(1)%z)* WCI)/ProfileDepth(L)
			!END
			ELSE IF ((WF(2)%z.LE.(LayerBottomZ(L) - ProfileDepth(L))).
     +				   AND.(WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WF(1)%WC 
			!END
			ELSE IF ((WF(1)%z.GE.(LayerBottomZ(L) - ProfileDepth(L))).
     +						AND.(WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(1)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(1)%WC +
     +									(LayerBottomZ(L) - WF(1)%z)
     +									* WCI)/ProfileDepth(L)
			ENDIF
	    !END //end if redist = 20 or 21
	    ELSE IF (redistStatus.EQ.30.OR.redistStatus.EQ.31) THEN
	    !BEGIN
			IF (WF(3)%z.GE.LayerBottomZ(L)) THEN
			!BEGIN
				LayerAvgWC(L) = WF(3)%WC 
			!END
			ELSE IF (WF(1)%z.LE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WCI 
			!END
			ELSE IF ((WF(3)%z.GE.(LayerBottomZ(L) -
     +						ProfileDepth(L))).AND.
     +						(WF(2)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(3)%z - LayerBottomZ(L) +
     +									  ProfileDepth(L))* WF(3)%WC +
     +									  (LayerBottomZ(L) - WF(3)%z)
     +									  * WF(2)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(3)%z.GE.(LayerBottomZ(L) -
     +						ProfileDepth(L))).AND.
     +						(WF(2)%z.LE.LayerBottomZ(L)).AND.
     +						(WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(3)%z - LayerBottomZ(L) +
     +									  ProfileDepth(L))* WF(3)%WC +
     +									  (WF(2)%z-WF(3)%z)*WF(2)%WC+
     +									  (LayerBottomZ(L) - WF(2)%z)*
     +									  WF(1)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(3)%z.GE.(LayerBottomZ(L) -
     +						  ProfileDepth(L))).AND.
     +						  (WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(3)%z - LayerBottomZ(L) +
     +									  ProfileDepth(L))* WF(3)%WC +
     +									  (WF(2)%z-WF(3)%z)*WF(2)%WC+
     +									  (WF(1)%z-WF(2)%z)*WF(1)%WC+
     +									  (LayerBottomZ(L) - WF(1)%z)*
     +									  WCI)/ProfileDepth(L)
			!END
			ELSE IF ((WF(3)%z.LE.(LayerBottomZ(L) -
     +					ProfileDepth(L))).AND.
     +					(WF(2)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WF(2)%WC 
			!END
			ELSE IF ((WF(2)%z.GE.(LayerBottomZ(L) -
     +					ProfileDepth(L))).AND.
     +					(WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(2)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(2)%WC +
     +									(LayerBottomZ(L) - WF(2)%z)*
     +									WF(1)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(2)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(2)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(2)%WC +
     +									(WF(1)%z-WF(2)%z)*WF(1)%WC +
     +									(LayerBottomZ(L) - WF(1)%z)*
     +									WCI)/ProfileDepth(L)
			!END
			ELSE IF ((WF(2)%z.LE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WF(1)%WC 
			!END
			ELSE IF ((WF(1)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(1)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(1)%WC +
     +									(LayerBottomZ(L) - WF(1)%z)*
     +									WCI)/ProfileDepth(L)
			ENDIF
	    !END //end if redist = 30 or 31
	    ELSE IF (redistStatus.EQ.40) THEN
	    !BEGIN
			IF (WF(4)%z.GE.LayerBottomZ(L)) THEN
			!BEGIN
				LayerAvgWC(L) = WF(4)%WC 
			!END
			ELSE IF (WF(1)%z.LE.(LayerBottomZ(L) -
     +					ProfileDepth(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WCI 
			!END
			ELSE IF ((WF(4)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(3)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(4)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(4)%WC +
     +									(LayerBottomZ(L) - WF(4)%z)*
     +									WF(3)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(4)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(3)%z.LE.LayerBottomZ(L)).AND.
     +					  (WF(2)%z.GE.LayerBottomZ(L)))THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(4)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(4)%WC +
     +									(WF(3)%z-WF(4)%z)*WF(3)%WC +
     +									(LayerBottomZ(L) - WF(3)%z)*
     +									WF(2)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(4)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(2)%z.LE.LayerBottomZ(L)).AND.
     +					  (WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(4)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(4)%WC +
     +									(WF(3)%z-WF(4)%z)*WF(3)%WC +
     +									(WF(2)%z-WF(3)%z)*WF(2)%WC +
     +									(LayerBottomZ(L) - WF(2)%z)*
     +									WF(1)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(4)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(4)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(4)%WC +
     +									(WF(3)%z-WF(4)%z)*WF(3)%WC +
     +									(WF(2)%z-WF(3)%z)*WF(2)%WC +
     +									(WF(1)%z-WF(2)%z)*WF(1)%WC +
     +									(LayerBottomZ(L) - WF(1)%z)*
     +									WCI)/ProfileDepth(L)
			!END
			ELSE IF ((WF(4)%z.LE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(3)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WF(3)%WC 
			!END
			ELSE IF ((WF(3)%z.GE.(LayerBottomZ(L) -
     +						ProfileDepth(L))).AND.
     +						(WF(2)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(3)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(3)%WC +
     +									(LayerBottomZ(L) - WF(3)%z)*
     +									WF(2)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(3)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(2)%z.LE.LayerBottomZ(L)).AND.
     +					  (WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(3)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(3)%WC +
     +									(WF(2)%z-WF(3)%z)*WF(2)%WC +
     +									(LayerBottomZ(L) - WF(2)%z)*
     +									WF(1)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(3)%z.GE.(LayerBottomZ(L) -
     +						ProfileDepth(L))).AND.
     +						(WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(3)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(3)%WC +
     +									(WF(2)%z-WF(3)%z)*WF(2)%WC +
     +									(WF(1)%z-WF(2)%z)*WF(1)%WC +
     +									(LayerBottomZ(L)- WF(1)%z)*
     +									WCI)/ProfileDepth(L)
			!END
			ELSE IF ((WF(3)%z.LE.(LayerBottomZ(L) -
     +						ProfileDepth(L))).AND.
     +						(WF(2)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WF(2)%WC 
			!END
			ELSE IF ((WF(2)%z.GE.(LayerBottomZ(L) -
     +						ProfileDepth(L))).AND.
     +						(WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(2)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(2)%WC +
     +									(LayerBottomZ(L) - WF(2)%z) *
     +									WF(1)%WC)/ProfileDepth(L)
			!END
			ELSE IF ((WF(2)%z.GE.(LayerBottomZ(L) -
     +						ProfileDepth(L))).AND.
     +						(WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) =  ((WF(2)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(2)%WC +
     +									(WF(1)%z-WF(2)%z)*WF(1)%WC +
     +									(LayerBottomZ(L) - WF(1)%z) *
     +									WCI)/ProfileDepth(L)
			!END
			ELSE IF ((WF(2)%z.LE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(1)%z.GE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = WF(1)%WC 
			!END
			ELSE IF ((WF(1)%z.GE.(LayerBottomZ(L) -
     +					  ProfileDepth(L))).AND.
     +					  (WF(1)%z.LE.LayerBottomZ(L))) THEN
			!BEGIN
				LayerAvgWC(L) = ((WF(1)%z - LayerBottomZ(L) +
     +									ProfileDepth(L))* WF(1)%WC +
     +									(LayerBottomZ(L) - WF(1)%z)*
     +									WCI)/ProfileDepth(L)
			ENDIF
	    ENDIF ! //end if redist = 40
	  ENDDO ! //end for l = 1 to numlayers
	
	END


cccccccccc calculates potential infitration rate ccccccccccccccccccccccccccccccccccccccccccccccccc
	Subroutine FpCalc
	USE ModuleDefs
	Implicit None

	IF (redistStatus.EQ.10.OR.redistStatus.EQ.11) THEN
		fp = ksat + ((ksat * Sav*(WCS-WCI))/WF(1)%FAmt)
      ELSE IF (redistStatus.EQ.20.OR.redistStatus.EQ.21) THEN
		fp = ksat + ((ksat * Sav*(WCS-WF(1)%WC))/WF(2)%FAmt)
      ELSE IF (redistStatus.EQ.30.OR.redistStatus.EQ.31) THEN
		fp = ksat + ((ksat * Sav*(WCS-WF(2)%WC))/WF(3)%FAmt)
      ELSE
		fp = ksat + ((ksat * Sav*(WCS-WF(3)%WC))/WF(4)%FAmt)
                      
	ENDIF
	END


ccccccccccccccc calculates the generic gamma coefficient ccccccccccccccccccccccccccccc
	Subroutine AdjFactorCalc(num,redistT,ADJ) 
	USE ModuleDefs
	Implicit None
	
	Integer num
	Double Precision redistT, ADJ
	Double Precision a,b,c,d,e,f
	Double Precision alpha,beta,delta
	
	IF (num.EQ.0) THEN
		ADJ = 0.0
	ELSE
	
		!//wave 86
		a = 4.2951615111 
		b = 154.6101111175 
		c = 0.0020393887 
		d = -0.0010402988
		e = -14.0032425382
		f = -61.5428564782

		alpha = ksat/(a* ksat + b)
		beta = c + d*(ksat**0.5)
		delta = ksat/(e* ksat + f)
		ADJ = alpha + beta*log(redistT) + delta/num
	ENDIF

	IF (ADJ.LT.0.0) THEN
		ADJ = 0.0
	ENDIF

	END


cccccccccccc calculates change in water content during redistribution ccccccccccccccccccccccc
	Subroutine Derivs(Rh, F, WC, y, dydx)
	Use ModuleDefs
	Implicit None
	
	Double Precision Rh, F, WC, y, dydx
	Double Precision RelSatI, RelSatO, expon1, expon2 
	Double Precision GWC, KWCI, KWCO, z 

	z = F/(y - wc)
	RelSatI = (wc - WCR)/(WCS - WCR)
	RelSatO = (y - WCR)/(WCS - WCR)
	expon1 = 3 + (1/lambda)
	GWC = Sav*((relSatO**expon1) - (relSatI**expon1))/
     +			(1-(relSatI**expon1))
	expon2 = 3 + (2/lambda)
	KWCI = ksat*(relSatI**expon2)
	KWCO = ksat*(relSatO**expon2)
	dydx =(Rh - KWCI - KWCO - (ksat*GWC/z))/(z) 
	END


cccccc Procedures RK4, RKQC and OdeInt perform 4th order Runge Kutta procedure
cccccc Adapted from:  Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1989.
cccccc Numerical Recipes in Pascal, Cambridge University Press: Port Chester, New York.

	Subroutine RK4(Rh,F,WC,y,dydx,h,yout)
c	inputs x = start time, y = y @ x, dydx = dydx @ x, h = time step
c	output yout
	Implicit None
	
	Double Precision hh,dym,dyt,yt,wc,y,dydx,h,yout
	Double Precision Rh,F
	
	hh = h/2.0
	yt = y + hh*dydx
	CALL derivs(Rh,F,WC,yt,dyt)
	yt = y + hh*dyt
  	CALL derivs(Rh,F,WC,yt,dym)
  	yt = y + h*dym
  	dym = dym + dyt
  	CALL derivs(Rh,F,WC, yt, dyt)
  	yout = y + h*(dydx + dyt + 2.0 * dym )/6.0

	END

	Subroutine RKQC(Rh,F,WC,y,dydx,x,htry,eps,yscal,hdid, hnext)
c	inputs : y and dydx @ x1; htry; eps; yscal
c	outputs : y and x replaced by new values; hdid is stepsize used; hnext is next stepsize
	Implicit None
	
	Double Precision pgrow,pshrnk,fcor,safety,errcon
	Double Precision xsav, hh, h, temp, errmax
	Double Precision dysav, ysav, ytemp, Rh,F
	Double Precision x, htry, eps,yscal,hdid,hnext,y,dydx,wc
  
	pgrow = -0.2
  	pshrnk = -0.25
  	fcor = 1.0/15.0
  	safety = 0.9
  	errcon = 6.0e-4
	xsav = x
  	ysav = y
  	dysav = dydx
  	h = htry
  	
	DO  
  	   	hh = h/2.0
    		CALL rk4(Rh,F,WC,ysav,dysav,hh,ytemp)
    		x = xsav + hh
    		CALL derivs(Rh,F,WC,ytemp,dydx)
    		CALL rk4(Rh,F,WC,ytemp,dydx,hh,y)
    		x = xsav + h
    	!test to see if x = xsav
    		CALL rk4(Rh,F,WC,ysav,dysav,h,ytemp)
    		errmax = 0.0
    		ytemp = y - ytemp
    		temp = abs(ytemp/yscal)
    		IF (errmax.LT.temp) THEN
      		errmax = temp
		ENDIF
    		errmax = errmax/eps
    		IF (errmax.LE.1.0) THEN
    			hdid = h
      		IF (errmax.GT.errcon) THEN
        			hnext = safety*h*exp(pgrow*log(errmax))
			ELSE
        			hnext = 4.0*h
			ENDIF
      		EXIT
    		ELSE
      		h = safety*h*exp(pshrnk*log(errmax))
		ENDIF
  	ENDDO
  	y = y + ytemp*fcor
	END

	Subroutine OdeInt(Rh,F,WC,ystart,x1, x2,eps,h1)
	Implicit None
	
	Double Precision x,hnext,tiny,hdid,h,yscal,y,dydx
	Double Precision Rh,F,wc,ystart,x1,x2,eps,h1
	Integer MaxStp,nstp
  	
	MaxStp = 100
  	tiny = 1.0e-3
  	x = x1
  	h = abs(h1)
  	y = ystart
  	DO nstp = 1 , maxstp 
  	  	CALL derivs(Rh,F,WC,y,dydx)
    		yscal = abs(y) + abs(dydx*h)+tiny
    		IF ((x+h-x2)*(x+h-x1).GT.0.0) THEN
      		h = x2-x
		ENDIF
    		CALL rkqc(Rh,F,WC,y,dydx,x,h,eps,yscal,hdid,hnext)
    		IF ((x-x2)*(x2-x1).GE.0.0) THEN
    			ystart = y
      		GOTO 99
    		ENDIF
     		h = hnext
  	ENDDO
99    MaxStp = 0
	END


cccccc Procedures Redist11, Redist20, Redist21, Redist30, Redist31 and Redist40 perform
cccccc the redistribution process depending on what redistribution status the system is in---}

	Subroutine Redist11(rh, ts,te,dt)
	USE ModuleDefs
	Implicit None
	Double Precision rh, ts, te,dt
	Double Precision y
	
	y = WF(1)%WCHold
	CALL OdeInt(rh,WF(1)%FAmt,WCI,y,ts,te,0.00001d0,dT)
	
	WF(1)%WCHold = y
	WF(1)%WC = y - wf(1)%adjFactor
	IF  (WF(1)%WC.LT.WCMin) THEN
		WF(1)%WC = WCMin
		IF (WF(1)%WCHold.LT.WCMin) THEN
			WF(1)%WCHold = WCMin
		ENDIF
	ENDIF
	
	CALL calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
                          
	CALL FpCalc !//update fp since wc has changed
	END
 
	Subroutine Redist20(rh,ts,te,dt)
c----- redist routine if two fronts & only one front redist -------------------------------------------------
	USE ModuleDefs
	Implicit None

	Double Precision y, rh
	Double Precision ts, te, dt
	Integer n, j

  	!redist wf 1
  	y = WF(1)%WChold
  	CALL OdeInt(0.d0,WF(1)%Famt,WCI,y,ts,te,0.00001d0,dT)
  	WF(1)%WChold = y
  	WF(1)%WC = y - wf(1)%adjFactor
	IF  (WF(1)%WC.LT.wcmin) THEN
    		WF(1)%WC = wcmin
       	IF (WF(1)%WC.LT.wcmin) THEN
         		WF(1)%WC = wcmin
		ENDIF
	ENDIF
  	
	Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
      
cccccc wf 2 does not redist since it is infiltrating
	Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
	
	!check conditions
  	IF (WF(2)%z.GE.WF(1)%z)THEN
  	!2 wf to 1 wf
    		WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    		WF(1)%WC = WF(2)%WC
    		WF(1)%WChold = WF(2)%WChold
    		Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
		WF(2)%FAmt = 0.0
		WF(2)%WC = 0.0
		WF(2)%WCHold = 0.0
		WF(2)%z = 0.0
		redistStatus = 10
  	ENDIF !if z2 > z1

	Call FpCalc
	
	END

	Subroutine Redist21(rh,ts,te,dt)
c----- redist routine if two fronts & both fronts redist --------------------------------------------
	USE ModuleDefs
	Implicit None

	Double Precision y,rh
	Double Precision ts, te, dt
	Integer n,j
  	
	!!decrease wc1 first
  	y = WF(1)%WCHold
  	CALL OdeInt(0.d0,WF(1)%FAmt,WCI,y,ts,te,0.00001d0,dT)
  	WF(1)%WCHold = y
  	WF(1)%WC = y  - wf(1)%adjFactor
  	IF  (WF(1)%WC.LT.wci) THEN
  		WF(1)%WC = wcmin
    		IF (WF(1)%WCHold.LT.wcmin)THEN
       		WF(1)%WCHold = wcmin
		ENDIF
  	ENDIF
  	Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
  	
	!!//decrease wc2 second
  	y = WF(2)%WCHold
  	CALL OdeInt(rh,WF(2)%FAmt,WF(1)%WCHold,y,ts,te,0.00001d0,dT)
  	WF(2)%WCHold = y
  	WF(2)%WC = y - wf(2)%adjFactor
  	Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
	
	IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC.OR.
     +	WF(2)%z.GE.WF(1)%z)THEN !! //merge 1 an 2
 	 	!begin two fronts to one
    		IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC)THEN
          	!BEGIN
      		WF(1)%WC = WF(1)%WC ! //wf one predominates
!      		IF (WF(2)%WCHold.LT.wcmin) THEN
!         			WF(2)%WCHold = wcmin
!			ENDIF
    		ELSE
    		!BEGIN
     			WF(1)%WC = WF(2)%WC !wf two predominates
      		WF(1)%WCHold = WF(2)%WCHold
    		ENDIF
    		WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    		Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
		WF(2)%wc = 0.0
		WF(2)%wchold = 0.0
		WF(2)%Famt = 0.0
		WF(2)%z = 0.0
		redistStatus = 11
		WF(1)%NumRedist = WF(2)%NumRedist
		WF(1)%redistTime = WF(2)%redistTime
		WF(1)%adjFactor = WF(2)%adjFactor
  	! ENDIF wc2 < wcmin
  	ENDIF
	CALL FpCalc	
	
	END


	Subroutine Redist30(rh,ts,te,dt)
c----- redist routine if three fronts & only two fronts redist ---------------------------------------------
	USE ModuleDefs
	Implicit None

	Double Precision y
	Double Precision ts, te, dt, rh
	Integer j, n

	!!decrease wc1 first
  	y = WF(1)%WCHold
  	CALL OdeInt(0.d0,WF(1)%FAmt,WCI,y,ts,te,0.00001d0,dT)
  	WF(1)%WCHold = y
  	WF(1)%WC = y  - wf(1)%adjFactor
  	IF  (WF(1)%WC.LT.wci) THEN
  		WF(1)%WC = wcmin
    		IF (WF(1)%WCHold.LT.wcmin)THEN
       		WF(1)%WCHold = wcmin
		ENDIF
  	ENDIF
  	Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
  	
	!!//decrease wc2 second
  	y = WF(2)%WCHold
  	CALL OdeInt(0.d0,WF(2)%FAmt,WF(1)%WCHold,y,ts,te,0.00001d0,dT)
  	WF(2)%WCHold = y
  	WF(2)%WC = y - wf(2)%adjFactor
  	Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)

  	!!checks
  	IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC.OR.
     +	WF(2)%z.GE.WF(1)%z)THEN !!merge 1 an 2
  	!! 3 fronts down to 2 fronts
    		IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC)THEN
        !BEGIN
      		WF(1)%WC = WF(1)%WC ! //wf one predominates
!      		IF (WF(2)%WCHold.LT.wcmin) THEN
!         			WF(2)%WCHold = wcmin
!			ENDIF
    		ELSE
    		!BEGIN
     			WF(1)%WC = WF(2)%WC !wf two predominates
      		WF(1)%WCHold = WF(2)%WCHold
    		ENDIF
    		WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    		Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)	
		redistStatus = 20
    		WF(2)%FAmt = WF(3)%FAmt
		WF(2)%WC = WF(3)%WC 
		WF(2)%WCHold = WF(3)%WCHold
		WF(2)%z = WF(3)%z
		WF(3)%WC = 0.0 
		WF(3)%WCHold = 0.0 
		WF(3)%FAmt = 0.0
		WF(3)%z = 0.0
		WF(1)%NumRedist = WF(2)%NumRedist
		WF(1)%redistTime = WF(2)%redistTime
		WF(1)%adjFactor = WF(2)%adjFactor
    	!!now move wf 2 down which was wf 3 and is only infil, not redist
    		Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
    	!!check
    		IF (WF(2)%z.GE.WF(1)%z) THEN
    	!!BEGIN //two fronts to one
      		WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    			WF(1)%WC = WF(2)%WC
			WF(1)%WCHold = WF(2)%WCHold
			Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
			WF(2)%WC = 0.0
			WF(2)%WCHold = 0.0
			WF(2)%FAmt = 0.0 
			WF(2)%z = 0.0 
      		redistStatus = 10
     		ENDIF
  	!!END//wc2 < wcmin,z2>z1
  	ELSE !!//still 3 fronts
  	!BEGIN
    			 !just move down since only infil
		Call calcRedistZf(WF(3)%FAmt,WF(3)%WC,WF(2)%WC,WF(3)%z) 
	    	!!checks
    		IF (WF(3)%z.GE.WF(2)%z) THEN
    	!!BEGIN //three fronts to two
      		WF(2)%FAmt = WF(2)%FAmt + WF(3)%FAmt
    			WF(2)%WC = WF(3)%WC
			WF(2)%WCHold = WF(3)%WCHold
			Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
    			WF(3)%wc = 0.0
			WF(3)%WCHold = 0.0
			WF(3)%FAmt = 0.0 
      		redistStatus = 20		
			
			IF (WF(2)%z.GE.WF(1)%z) THEN	!!recheck 
				WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    				WF(1)%WC = WF(2)%WC
				WF(1)%WCHold = WF(2)%WCHold
				Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
				WF(2)%WC = 0.0
				WF(2)%WCHold = 0.0
				WF(2)%FAmt = 0.0 
				WF(2)%z = 0.0
      			redistStatus = 10
			ENDIF
	!!endif 3 > 2
		ENDIF
	ENDIF !!else still 3 fronts

	CALL FpCalc

	END


	Subroutine Redist31(rh,ts,te,dt)
c----- redist routine if three fronts & all fronts redist ------------------------------------------------------
	USE ModuleDefs
	Implicit None

	Double Precision y
	Double Precision ts, te, dt, rh
	Integer j, n

	!!decrease wc1 first
  	y = WF(1)%WCHold
  	CALL OdeInt(0.d0,WF(1)%FAmt,WCI,y,ts,te,0.00001d0,dT)
  	WF(1)%WCHold = y
  	WF(1)%WC = y  - wf(1)%adjFactor
  	IF  (WF(1)%WC.LT.wci) THEN
  		WF(1)%WC = wcmin
    		IF (WF(1)%WCHold.LT.wcmin)THEN
       		WF(1)%WCHold = wcmin
		ENDIF
  	ENDIF
  	Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
  	
	!!//decrease wc2 second
  	y = WF(2)%WCHold
  	CALL OdeInt(0.d0,WF(2)%FAmt,WF(1)%WCHold,y,ts,te,0.00001d0,dT)
  	WF(2)%WCHold = y
  	WF(2)%WC = y - wf(2)%adjFactor
  	Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
     	!!checks
  	  	IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC.OR.
     +	WF(2)%z.GE.WF(1)%z)THEN !!merge 1 an 2
  	!! 3 fronts down to 2 fronts
    		IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC)THEN
        !BEGIN
      		WF(1)%WC = WF(1)%WC ! //wf one predominates
!      		IF (WF(2)%WCHold.LT.wcmin) THEN
!         			WF(2)%WCHold = wcmin
!			ENDIF
    		ELSE
    		!BEGIN
     			WF(1)%WC = WF(2)%WC !wf two predominates
      		WF(1)%WCHold = WF(2)%WCHold
    		ENDIF
    		WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    		Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)	
		redistStatus = 21
    		WF(2)%FAmt = WF(3)%FAmt
		WF(2)%WC = WF(3)%WC 
		WF(2)%WCHold = WF(3)%WCHold
		WF(2)%z = WF(3)%z
		WF(3)%WC = 0.0 
		WF(3)%WCHold = 0.0 
		WF(3)%FAmt = 0.0
		WF(3)%z = 0.0
		WF(1)%NumRedist = WF(2)%NumRedist
		WF(1)%redistTime = WF(2)%redistTime
		WF(1)%adjFactor = WF(2)%adjFactor
		WF(2)%NumRedist = WF(3)%NumRedist
		WF(2)%redistTime = WF(3)%redistTime
		WF(2)%adjFactor = WF(3)%adjFactor


	   	!!//need to decrease wc2 which was wc3
    		y = WF(2)%WCHold
  		CALL OdeInt(rh,WF(2)%FAmt,WF(1)%WCHold,y,ts,te,0.00001d0,dT)
  		WF(2)%WCHold = y
		WF(2)%WC = y - WF(2)%adjFactor
  	    Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
		
    		!!checks
  		IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC.OR.
     +		WF(2)%z.GE.WF(1)%z)THEN !!merge 1 an 2
  		!! 2 fronts down to 1 fronts
    			IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC)THEN
        		!BEGIN
      			WF(1)%WC = WF(1)%WC! //wf one predominates
!      			IF (WF(2)%WC.LT.wcmin) THEN
!         				WF(2)%WCHold = wcmin
!            		ENDIF
    			ELSE
    		!BEGIN
     				WF(1)%WCHold = WF(2)%WCHold !wf two predominates
      			WF(1)%WC = WF(2)%WC
    			ENDIF
    			WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    			Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
			WF(2)%WC = 0.0
			WF(2)%WCHold = 0.0
			WF(2)%FAmt = 0.0
			WF(2)%z = 0.0 
      		redistStatus = 11
			WF(1)%NumRedist = WF(2)%NumRedist
			WF(1)%redistTime = WF(2)%redistTime
			WF(1)%adjFactor = WF(2)%adjFactor
		!!END //wc2 < wcmin; z2 > z1 for 2 to 1fronts
	    ENDIF !! //2 fronts remain

  	!!END //wc2 < wcmin ; z2>z1 3 to 2 fronts
  	ELSE !!//still 3 wf for now
  	!!BEGIN
    	!!//decrease wc3 third
    		y = WF(3)%WCHold
  		CALL OdeInt(rh,WF(3)%FAmt,WF(2)%WCHold,y,ts,te,0.00001d0,dT)
  		WF(3)%WCHold = y
  		WF(3)%WC = y - WF(3)%adjFactor
  		Call calcRedistZf(WF(3)%FAmt,WF(3)%WC,WF(2)%WC,WF(3)%z)
    	 	IF  (WF(3)%WC.LE.WF(2)%WC.OR.WF(3)%z.GE.
     +		 WF(2)%z) THEN !!  //merge 3 and 2
    	!!BEGIN  // 3 fronts to 2
      		IF (WF(3)%WC.LE.WF(2)%WC) THEN 
		   		WF(2)%WC = WF(2)%WC
      		ELSE
      			!!BEGIN
        			WF(2)%WC = WF(3)%WC 
        			WF(2)%WCHold = WF(3)%WCHold 
      		ENDIF
      		WF(2)%FAmt = WF(2)%FAmt + WF(3)%FAmt
			Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
    			WF(3)%WCHold = 0.0
     			WF(3)%WC = 0.0
			WF(3)%FAmt = 0.0
			WF(3)%z = 0.0 
      		redistStatus = 21		    		
			WF(2)%NumRedist = WF(3)%NumRedist
			WF(2)%redistTime = WF(3)%redistTime
			WF(2)%adjFactor = WF(3)%adjFactor	
					
			IF (WF(2)%z.GE.WF(1)%z) THEN	!!recheck 
				WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    				WF(1)%WC = WF(2)%WC
				WF(1)%WCHold = WF(2)%WCHold
				Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
     				WF(2)%WC = 0.0
				WF(2)%WCHold = 0.0
				WF(2)%z = 0.0
				WF(2)%FAmt = 0.0
				WF(1)%NumRedist = WF(2)%NumRedist
				WF(1)%redistTime = WF(2)%redistTime
				WF(1)%adjFactor = WF(2)%adjFactor
      			redistStatus = 11
			ENDIF
			!!//no need to decrease or do anymore checks since taken care of above
		ENDIF !! //wc2 < wcmin; z2 > z1 for 3 to 2 fronts
  	ENDIF !! else, still 3 wf
	
	CALL FpCalc	
	
	END


	Subroutine Redist40(rh,ts,te,dt)
c----- redist routine if four fronts & only three fronts redist -----------------------------------------------------
	USE ModuleDefs
	Implicit None

	Double Precision y
	Double Precision ts, te, dt, rh
	Integer j, n

	!!decrease wc1 first
  	y = WF(1)%WCHold
  	CALL OdeInt(0.d0,WF(1)%FAmt,WCI,y,ts,te,0.00001d0,dT)
  	WF(1)%WCHold = y
  	WF(1)%WC = y  - wf(1)%adjFactor
  	IF  (WF(1)%WC.LT.wci) THEN
  		WF(1)%WC = wcmin
    		IF (WF(1)%WCHold.LT.wcmin)THEN
       		WF(1)%WCHold = wcmin
		ENDIF
  	ENDIF
  	Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
  	
	!!//decrease wc2 second
  	y = WF(2)%WCHold
  	CALL OdeInt(0.d0,WF(2)%FAmt,WF(1)%WCHold,y,ts,te,0.00001d0,dT)
  	WF(2)%WCHold = y
  	WF(2)%WC = y - wf(2)%adjFactor
  	Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
     	!!checks
  	  	IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC.OR.
     +	WF(2)%z.GE.WF(1)%z)THEN !!merge 1 an 2
  	!! 3 fronts down to 2 fronts
    		IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC)THEN
        !BEGIN
      		WF(1)%WC = WF(1)%WC ! //wf one predominates
!      		IF (WF(2)%WCHold.LT.wcmin) THEN
!         			WF(2)%WCHold = wcmin
!			ENDIF
    		ELSE
    		!BEGIN
     			WF(1)%WC = WF(2)%WC !wf two predominates
      		WF(1)%WCHold = WF(2)%WCHold
    		ENDIF
    		WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    		Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)	
		redistStatus = 30
    		WF(2)%FAmt = WF(3)%FAmt
		WF(2)%WC = WF(3)%WC 
		WF(2)%WCHold = WF(3)%WCHold
		WF(2)%z = WF(3)%z
		WF(3)%WC = WF(4)%WC 
		WF(3)%WCHold = WF(4)%WCHold 
		WF(3)%FAmt = WF(4)%FAmt
		WF(3)%z = WF(4)%z
		WF(4)%WC = 0.0 
		WF(4)%WCHold = 0.0 
		WF(4)%FAmt = 0.0
		WF(4)%z = 0.0
		WF(1)%NumRedist = WF(2)%NumRedist
		WF(1)%redistTime = WF(2)%redistTime
		WF(1)%adjFactor = WF(2)%adjFactor
		WF(2)%NumRedist = WF(3)%NumRedist
		WF(2)%redistTime = WF(3)%redistTime
		WF(2)%adjFactor = WF(3)%adjFactor
		!//no need to change new wf3 since infil only
 		
    	!!//need to decrease wc2 which was wc3
    		y = WF(2)%WCHold
  		CALL OdeInt(0.d0,WF(2)%FAmt,WF(1)%WCHold,y,ts,te,0.00001d0,dT)
  		WF(2)%WCHold = y
  		WF(2)%WC = y - WF(2)%adjFactor
  		Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
		
   		!!checks
  		IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC.OR.
     +		WF(2)%z.GE.WF(1)%z)THEN !!merge 1 an 2
  			!! 3 fronts down to 2 fronts
    			IF (WF(2)%WC.LT.wcmin.OR.WF(2)%WC.LE.WF(1)%WC)THEN
        		!BEGIN
      			WF(1)%WC = WF(1)%WC ! //wf one predominates
!      			IF (WF(2)%WC.LT.wcmin) THEN
!         				WF(2)%WC = wcmin
!            		ENDIF
    			ELSE
    			!BEGIN
     				WF(1)%WC = WF(2)%WC !wf two predominates
      			WF(1)%WCHold = WF(2)%WCHold
    			ENDIF
    			
			WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    			Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
			redistStatus = 20
    			WF(2)%FAmt = WF(3)%FAmt 
			WF(2)%WC = WF(3)%WC 
			WF(2)%WCHold = WF(3)%WCHold
			WF(3)%WC = 0.0 
			WF(3)%WCHold = 0.0 
			WF(3)%FAmt = 0.0
			WF(3)%z = 0.0
			WF(1)%NumRedist = WF(2)%NumRedist
			WF(1)%redistTime = WF(2)%redistTime
			WF(1)%adjFactor = WF(2)%adjFactor
		
		!!//move down z2; was originally 4th wf so do not decrease
      		Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)	
      		!!//check last time
      		IF (WF(2)%z.GE.WF(1)%z)THEN 
      			!!BEGIN
    				WF(1)%WC = WF(2)%WC
				WF(1)%WCHold = WF(2)%WCHold
				WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    				Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
    				WF(2)%WC = 0.0
				WF(2)%WCHold = 0.0
				WF(2)%FAmt = 0.0 
      			redistStatus = 10
      		ENDIF
    		!!END //wc2 < sd%wcmin; z2 > z1 for 3 to 2 fronts
    		ELSE !!//still 3 fronts for now
    		!!BEGIN
    	  	!!//move down z3; originally wf 4; do not redist, only infil
      		Call calcRedistZf(WF(3)%FAmt,WF(3)%WC,WF(2)%WC,WF(3)%z)
      		IF (WF(3)%z.GE.WF(2)%z) THEN 
      			!!BEGIN //merge 3 and 2; wf1 stays the same
    				WF(2)%FAmt = WF(2)%FAmt + WF(3)%FAmt
    				WF(2)%WC = WF(3)%WC
				WF(2)%WCHold = WF(3)%WCHold
				Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
    				WF(3)%WC = 0.0
				WF(3)%WCHold = 0.0
				WF(3)%FAmt = 0.0
				WF(3)%z = 0.0
				redistStatus = 20		 
				IF (WF(2)%z.GE.WF(1)%z) THEN	!!recheck 
					WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    					WF(1)%WC = WF(2)%WC
					WF(1)%WCHold = WF(2)%WCHold
					Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
    					WF(2)%WC = 0.0
					WF(2)%WCHold = 0.0
					WF(2)%FAmt = 0.0
					WF(2)%z = 0.0
					
      				redistStatus = 10
				ENDIF			      
      		ENDIF !!//z3>z2
    		ENDIF !!; //still  3 fronts for now
  	!!END //wc2 < wcmin ; z2>z1 4 to 3 fronts
  	ELSE !! //still 4 wf for now
  	!!BEGIN
    	!!//decrease wc3 third
    		y = WF(3)%WCHold
    		CALL OdeInt(0.d0,WF(3)%FAmt,WF(2)%WCHold,y,ts,te,0.00001d0,dT)
    		WF(3)%WCHold = y
    		WF(3)%WC = y - WF(3)%adjFactor
    		Call calcRedistZf(WF(3)%FAmt,WF(3)%WC,WF(2)%WC,WF(3)%z)
    		IF  (WF(3)%WC.LE.WF(2)%WC.OR.WF(3)%z.GE.
     +     		WF(2)%z) THEN !merge 3 and 2
    		!!BEGIN  // 4 fronts to 3
      		IF (WF(3)%WC.GE.WF(2)%WC) THEN 
        			WF(2)%WC = WF(2)%WC 
      		ELSE
      		!!BEGIN
        			WF(2)%WC = WF(3)%WC
        			WF(2)%WCHold = WF(3)%WCHold 
      		ENDIF
      		WF(2)%FAmt = WF(2)%FAmt + WF(3)%FAmt
    			Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)	
	     	redistStatus = 30
			WF(2)%NumRedist = WF(3)%NumRedist
			WF(2)%redistTime = WF(3)%redistTime
			WF(2)%adjFactor = WF(3)%adjFactor

      		IF (WF(2)%z.GE.WF(1)%z) THEN	!!recheck 
					!!merge 1 and 2 which was 3, move down 4 to 2
				WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    				WF(1)%WC = WF(2)%WC
				Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)	
				WF(1)%NumRedist = WF(2)%NumRedist
				WF(1)%redistTime = WF(2)%redistTime
				WF(1)%adjFactor = WF(2)%adjFactor
				
				WF(2)%FAmt = WF(4)%FAmt !!wf2 now wf4 since 3 was merged w/ 2 
        			WF(2)%WC = WF(4)%WC
				WF(2)%WCHold = WF(4)%WCHold

				Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,WF(2)%z)
      			redistStatus = 20
				WF(3)%WC = 0.0
				WF(3)%WCHold = 0.0
				WF(3)%FAmt = 0.0
				WF(3)%z = 0.0			
				WF(4)%WC = 0.0
				WF(4)%WCHold = 0.0
				WF(4)%FAmt = 0.0
				WF(4)%z = 0.0	
				
				IF (WF(2)%z.GE.WF(1)%z) THEN	!!recheck 
					WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
    					WF(1)%WC = WF(2)%WC
					WF(1)%WCHold = WF(2)%WC
					Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,WF(1)%z)
					WF(2)%WC = 0.0
					WF(2)%WCHold = 0.0
					WF(2)%FAmt = 0.0
					WF(2)%z = 0.0
					redistStatus = 10
				ENDIF			      
			!end recheck		
			ELSE !still 3 fronts, move down 4 front which is now 3
				!!//wf 1 stays the same; wf4 becomes wf3
        			WF(3)%WC = WF(4)%WC
				WF(3)%WCHold = WF(4)%WCHold
				WF(3)%FAmt = WF(4)%FAmt
				Call calcRedistZf(WF(3)%FAmt,WF(3)%WC,WF(2)%WC,
     +									WF(3)%z)	
      			WF(4)%WC = 0.0
				WF(4)%WCHold = 0.0
				WF(4)%FAmt = 0.0
				WF(4)%z = 0.0
				
      			redistStatus = 30      		
      			IF (WF(3)%z.GE.WF(2)%z) THEN 
      			!!BEGIN  //3 fronts to 2
      				WF(2)%WC = WF(3)%WC
					WF(2)%WCHold = WF(3)%WCHold
					WF(2)%FAmt = WF(2)%FAmt + WF(3)%FAmt
      				Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(1)%WC,
     +										WF(2)%z)	
      				WF(3)%WC = 0.0
					WF(3)%WCHold = 0.0
					WF(3)%FAmt = 0.0
					WF(3)%z = 0.0
					redistStatus = 20             		
        				IF (WF(2)%z.GE.WF(1)%z) THEN	!!recheck 
						WF(1)%WC = WF(2)%WC
						WF(1)%WC = WF(2)%WC
						WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt	
    						Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,
     +											WF(1)%z)
						WF(2)%WCHold = 0.0
						WF(2)%WC = 0.0
						WF(2)%FAmt = 0.0
						WF(2)%z = 0.0
						redistStatus = 10
					ENDIF		
				
				ENDIF
      		ENDIF
    		!!END// wc3 >wc2 or z3 > z2 ;4 to 3
    		ELSE !! //still 4 fronts
    		!!BEGIN
      		!!//calculate wf 4 depth
      		Call calcRedistZf(WF(4)%FAmt,WF(4)%WC,WF(3)%WC,WF(4)%z)
      		IF  (WF(4)%z.GE.WF(3)%z) THEN 
      		!!BEGIN  // 4 fronts to 3 ; 
      			!!//wf 1 stays the same; wf4 becomes wf3
        			WF(3)%WC = WF(4)%WC
				WF(3)%WCHold = WF(4)%WCHold
				WF(3)%FAmt = WF(3)%FAmt + WF(4)%FAmt
				Call calcRedistZf(WF(3)%FAmt,WF(3)%WC,WF(2)%WC,WF(3)%z)	
      			WF(4)%WC = 0.0
				WF(4)%WCHold = 0.0
				WF(4)%FAmt = 0.0
				WF(4)%z = 0.0
				redistStatus = 30      				
				IF (WF(3)%z.GE.WF(2)%z) THEN !check if 3 > 2
					WF(2)%WC = WF(3)%WC
					WF(2)%WCHold = WF(3)%WCHold
					WF(2)%FAmt = WF(2)%FAmt + WF(3)%FAmt
					Call calcRedistZf(WF(2)%FAmt,WF(2)%WC,WF(3)%WC,
     +										WF(2)%z)	
      				WF(3)%WC = 0.0
					WF(3)%WCHold = 0.0
					WF(3)%FAmt = 0.0
					WF(3)%z = 0.0
      				redistStatus = 20 
					IF (WF(2)%z.GE.WF(1)%z) THEN	!!recheck 
						WF(1)%WC = WF(2)%WC
						WF(1)%WCHold = WF(2)%WCHold
						WF(1)%FAmt = WF(1)%FAmt + WF(2)%FAmt
						Call calcRedistZf(WF(1)%FAmt,WF(1)%WC,WCI,
     +											WF(1)%z)	
      					WF(2)%WC = 0.0
						WF(2)%WCHold = 0.0
						WF(2)%FAmt = 0.0
						WF(2)%z = 0.0
      					redistStatus = 10
					ENDIF		
      			ENDIF
			ENDIF !!//z4 > z3 ; 
      	
    		ENDIF !! //else, still 4 fronts
  	ENDIF !!// else, still 4 fronts for now
	
	CALL FpCalc	
	
	END

cccccccc	performs redistribution process when precip < ks and no ponding ccccccccccccccccccc
	Subroutine Redistribution(tend, dT)
	USE ModuleDefs
	Implicit None
	Double Precision tend, dT
	Double Precision ts, te
	INTEGER J
	  ts = tend - dT
	  te = tend
	!calculate gamma coefficient for each wetting front
	  DO J = 1,3 
	     Call adjFactorCalc(wf(j)%numRedist,wf(j)%redistTime,
     +								wf(j)%adjFactor)
	  ENDDO
	  IF (redistStatus.EQ.11) THEN
	  !BEGIN
	    CALL Redist11(0.D0,ts,te,dT)
	  !END //end redist 11
	  ELSE IF (redistStatus.EQ.21) THEN
	  !BEGIN
	    CALL Redist21(0.D0,ts,te,dT)
	  !END  //end redist = 21
	  ELSE IF (redistStatus.EQ.31) THEN
	  !BEGIN
	    CALL Redist31(0.D0,ts,te,dT)
	  ENDIF !  //end redist = 31 or 40

	END



