	Subroutine ThreeDGAR

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: 3-dimensional Green-Ampt infiltration and redistribution method (3DGAR)c
c	Code is written in Fortran language											    c	
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	
	USE ModuleDefs3								
	Implicit None

	Integer j,redistStatus
	Integer k, i, nTimeSteps!, runFlag
	Double Precision dM !,Pi
	Double Precision deltat,ts,tstop,tend
	Double Precision EmitterRates(5000)
!	Parameter (RunFlag = 1)
!	Parameter (Pi = 3.14159265358979d0)

	Open(Unit = 500, file = 
     +	'Soil_In3D.txt',status = 'old')	
      Open(Unit = 200, file = '3DInfilOutput.txt')
	Open(unit = 300, file = '3DInfilOutput2.txt')
	Open(unit = 400, file = '3DInfilOutput3.txt')
	Read(500,*) SD%wcs,SD%wcr,SD%wci,SD%wcmin,SD%lambda,
     +			SD%ksat,SD%H,Q,nTimeSteps,tstop
	Read(500,*) (EmitterRates(j), j=1,nTimeSteps)	
	

	Do i = 1,ncols ! column
		Do j = 1,nrows ! row
			cell(j,i)%width = 5.0
			cell(j,i)%thick = 5.0
		EndDo
	EndDo

!cccccc just for testing purposes; root abstraction will be given by main calling program cccc	
	call RAGridSetup 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc	
	
	
	
	deltaT = tstop/nTimeSteps
	
	ts = 0.d0
	tend = ts + deltat
	
	j=1

	IF (runFlag.EQ.1) THEN

		Call ChuInitialize(redistStatus,dM)
	
		write (200,*) '90 Degree Line Results'
		write (200,*) 'Time(h) 5(cm) 10(cm) 15(cm) 20(cm) 25(cm) 30(cm)'
		write (400,*) 'tend Vol1 F1 WC1 RA1 Vol2 F2 WC2 RA2  '
!		write (300,*) '0 Degree Line Results'
!		write (300,*) 'Time(h) 5(cm) 10(cm) 15(cm) 20(cm) 25(cm) 30(cm)'
!		write (300,*) 'Time old new angle rate xrate yrate'

	EndIf

!	cell(1,1)%avgWCold = 0.38
!	cell(1,1)%avgWC = 0.3
!	cell(1,1)%area(1) = 1
!	cell(1,1)%area(2) = 3
!	cell(1,1)%centerX = 0.0
!	cell(1,1)%centerY = 1.0
!	wf3(1)%kint = 0.0
!	wf3(1)%a = 6.0
!	wf3(1)%b = 12.0
		
!	wf3(2)%kint = 0.0
!	wf3(2)%a = 2.0
!	wf3(2)%b = 4.0
!	redistStatus = 21
!	Call calcNormalAngle(redistStatus)
!	Call calcFlowRates(deltaT)
	
	Do While (tend.LE.tstop)
		Q = EmitterRates(j)
		Call ChuInfilTest(redistStatus,ts,tend,deltat)
!		cell(2,1)%normalangle = cell(2,1)%normalangle*180/pi
		write(200,7000) tend, (cell(k,1)%avgWC, k = 2,nrows)
!		write(300,7000) tend, (cell(1,k)%avgWC, k = 2,ncols)
!		write(300,8000) tend, cell(1,1)%AvgWcold, cell(1,1)%avgwc,
!     +		cell(1,1)%normalangle,cell(1,1)%rate%normalRate, 
!     +		cell(1,1)%rate%xdirrate, cell(1,1)%rate%ydirrate
        write(300,9000) tend, wf3(1)%R(1),wf3(1)%R(8),wf3(1)%vol,
     +                  wf3(1)%FAmt,wf3(1)%wc,wf3(1)%RootAbsActual
!        write(400,5000)   tend,(wf3(1)%R(k),k=1,8)
!        write(400,5000)   tend,(wf3(1)%x(k),k=1,8)
!        write(400,5000)   tend,(wf3(1)%y(k),k=1,8)
!     		write(400,*) tend, wf3(1)%kint,wf3(1)%a,wf3(1)%b
        write(400,5000) tend,(WF3(k)%Vol,WF3(k)%FAmt,WF3(k)%WC,
     +                      WF3(k)%RootAbsActual,k= 1,2)
       
     		ts = ts + deltat
		tend = tend + deltat
		j = j + 1
	ENDDO
	
	Write (200,*) ' '

7000	Format (F9.3, 15F6.3)
8000  Format (F9.3, 2F6.3,F8.2,3F9.3)
9000  Format (F9.3,2F7.1,F12.1,F10.1,F9.3,F10.5)
!5000  Format (F9.3,8F9.2)
5000  Format (F9.3,2F10.1,F8.3,3F10.1,F8.3,F10.1)



	Close (200)
	Close (300)
	Close (500)
	END
	
	Subroutine ChuInitialize(redistStatus,dM)
c-----initializes variables ---------------------------------------------------------------
	USE ModuleDefs3
	Implicit None
	
	Integer i,j,redistStatus
	Double Precision dM,ro
	
		
c-----Angles for points along wetting front ---------------------------------------------------------		
	angles(1) = (pi/8.d0)
	angles(2) = (3.d0*pi/16.d0)
	angles(3) = (pi/4.0)
	angles(4) = (5.d0*pi/16.d0)
	angles(5) = (3.d0*pi/8.d0)
	angles(6) = (7.d0*pi/16.d0)
	angles(7) = (pi/2.d0)
	
	
	dM = sd%wcs - sd%wci
	ro = Q*0.0265*(sd%H**(-0.6098))*(sd%Ksat**(-0.6555))
!	print*, ro
c	Simplified Philip Romin Calculation
c	ro = 0.9327*Q*(sd%H**-1.0)*(sd%Ksat**-1.0)/(2.0*PI) 
	
	Do i = 1,4
		WF3(i)%FAmt = 0.0
		WF3(i)%wc = 0.0
		Do j = 1, numA+1
			WF3(i)%R(j) = 0.0
			WF3(i)%x(j) = 0.0
			WF3(i)%y(j) = 0.0	
          Enddo
		WF3(i)%kint = 0.0
		WF3(i)%a = 0.0
		WF3(i)%b = 0.0
		WF3(i)%vol = 0.0
		WF3(i)%SA = 0.0
		WF3(i)%TInfil = 0.0
		WF3(i)%RootAbs = 0.0
	EndDo
      
	WF3(1)%WC = sd%wcs
	redistStatus = 10
	Do i = 1,numA+1
		WF3(1)%R(i) = ro
	EndDo


c-----Average water content H2D soils -----------------------------------------------------------
c     If obs. pts fall on x and/or y axis so that top row and left column are half cells
!	Do i = 1,ncols ! column
!		Do j = 1,nrows ! row
!			If(i.eq.1.and.j.eq.1) Then !corner cell
!			  cell(j,i)%y1 = 0	! top coordinate
!			  cell(j,i)%y2 = cell(j,i)%thick/2.0	!bottom coordinate
!			  cell(j,i)%x1 = 0	!left coordinate
!			  cell(j,i)%x2 = cell(j,i)%width/2.0 ! right coordinate
!			  cell(j,i)%CenterX = cell(j,i)%width/4.0
!			  cell(j,i)%CenterY = cell(j,i)%thick/4.0
!			ElseIF(i.eq.1.and.j.ne.1) THEN !!!far left column
!			  cell(j,i)%y1 = cell(j,i)%thick/2.0 + 
!     +							(j-2)*cell(j,i)%thick
!			  cell(j,i)%y2 = cell(j,i)%y1 + cell(j,i)%thick
!			  cell(j,i)%x1 = 0	
!			  cell(j,i)%x2 = cell(j,i)%width/2.0
!			  cell(j,i)%CenterX = cell(j,i)%width/4.0
!			  cell(j,i)%CenterY = (j-1)*cell(j,i)%thick
!			ElseIF(j.eq.1.and.i.ne.1) THEN !!!top row
!			  cell(j,i)%y1 = 0
!			  cell(j,i)%y2 = cell(j,i)%thick/2.0
!			  cell(j,i)%x1 = cell(j,i)%width/2.0 + 
!     +							(i-2)*cell(j,i)%width	
!			  cell(j,i)%x2 = cell(j,i)%x1 + cell(j,i)%width
!			  cell(j,i)%CenterX = (i-1)*cell(j,i)%width
!			  cell(j,i)%CenterY = cell(j,i)%thick/4.0
!			Else !! all other cells
!			  cell(j,i)%y1 = cell(j,i)%thick/2.0 + 
!     +							  (j-2)*cell(j,i)%thick
!			  cell(j,i)%y2 = cell(j,i)%y1 + cell(j,i)%thick
!			  cell(j,i)%x1 = cell(j,i)%width/2.0 + 
!     +							  (i-2)*cell(j,i)%width	
!			  cell(j,i)%x2 = cell(j,i)%x1 + cell(j,i)%width
!			  cell(j,i)%CenterX = (i-1)*cell(j,i)%width
!			  cell(j,i)%CenterY = (j-1)*cell(j,i)%thick
!			EndIF	
			
!			
!			cell(j,i)%avgWC = SD%wci
!			cell(j,i)%avgWCold = SD%wci
			
!		Enddo
!	Enddo


c-----standard grid -------------------------------------------------------------------------------------

	Do i = 1,ncols !!column
		Do j = 1,nrows !!!row
			cell(j,i)%x1 = (i-1)*cell(j,i)%width
			cell(j,i)%x2 = 	cell(j,i)%x1 + cell(j,i)%width
			cell(j,i)%y1 = (j-1)*cell(j,i)%thick	
			cell(j,i)%y2 = cell(j,i)%y1 + cell(j,i)%thick
			cell(j,i)%CenterX = cell(j,i)%width/2.0 + 
     +								(i-1)*cell(j,i)%width
			cell(j,i)%CenterY = cell(j,i)%width/2.0 + 
     +								(j-1)*cell(j,i)%thick
			cell(j,i)%avgWC = SD%wci
			cell(j,i)%avgWCold = SD%wci
			
		Enddo
	Enddo
	
	RunFlag = 0
	
	
	END
	
	Subroutine ChuInfilTest (redistStatus,ts,tend,deltat)
c-----updates redistribution status; calls infiltration and redistribution process -----------------------	
	USE ModuleDefs3
	Implicit None
	
	Double Precision Rold(7),Rnew(7),SA,Vol
	Double Precision Sorp,a,b,times(6) !,pi
	Double Precision cummInfil, deltat, ts,tstop,tend
	Double Precision t, R0new,rootHold
	Integer redistStatus,n
	Integer i,j,k,k1,k2,k3
!	Parameter(Pi = 3.14159265359)

cccc---------remove after testing ---------------------------
cccc------------this used for setting up root abs grid --------

c--------------------------- Update redistStatus --------------------------------------------------		
		IF (redistStatus.EQ.10.AND.Q.EQ.0.d0)THEN
        		redistStatus = 11
		ELSEIF (redistStatus.EQ.11.AND.Q.GT.0.d0)THEN
			redistStatus = 20
			WF3(2)%WC = sd%wcs
			Do i = 1,numA+1
				WF3(2)%R(i) = romin
			EndDo
		ELSEIF (redistStatus.EQ.20.AND.Q.EQ.0.d0)THEN
        		redistStatus = 21
		ELSEIF (redistStatus.EQ.21.AND.Q.GT.0.d0)THEN
			redistStatus = 30
			WF3(3)%WC = sd%wcs
			Do i = 1,numA+1
				WF3(3)%R(i) = romin
			EndDo
		ELSEIF (redistStatus.EQ.30.AND.Q.EQ.0.d0)THEN
        		redistStatus = 31
		ELSEIF (redistStatus.EQ.31.AND.Q.GT.0.d0)THEN
			redistStatus = 40
			WF3(4)%WC = sd%wcs
			Do i = 1,numA+1
				WF3(4)%R(i) = romin
			EndDo
		ELSEIF (redistStatus.EQ.40.AND.Q.EQ.0.d0) THEN
        		redistStatus = 31
			n = 1
			WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
			WF3(1)%WC = WF3(2)%WC
			Call calcR(n,WF3(1)%WC-sd%wci,WF3(1)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(1)%y(numA+1))
			Call calcSA(n,WF3(1)%y(numA+1))
c-------------------------------- reconfigure R,x,y,vol ---------------------------------------------
			WF3(2)%TInfil = WF3(3)%TInfil 
			WF3(2)%WC = WF3(3)%WC
			Do j = 1, numA+1
				WF3(2)%R(j) = WF3(3)%R(j)
				WF3(2)%x(j) = WF3(3)%x(j)
				WF3(2)%y(j) = WF3(3)%y(j)	
			Enddo
			WF3(2)%kint = WF3(3)%kint
			WF3(2)%a = WF3(3)%a
			WF3(2)%b = WF3(3)%b
			WF3(2)%vol = WF3(3)%vol
		
			WF3(3)%TInfil = WF3(4)%TInfil 
			WF3(3)%WC = WF3(4)%WC
			Do j = 1, numA+1
				WF3(3)%R(j) = WF3(4)%R(j)
				WF3(3)%x(j) = WF3(4)%x(j)
				WF3(3)%y(j) = WF3(4)%y(j)	
			Enddo
			WF3(3)%kint = WF3(4)%kint
			WF3(3)%a = WF3(4)%a
			WF3(3)%b = WF3(4)%b
			WF3(3)%vol = WF3(4)%vol

			WF3(4)%FAmt = 0.0
			WF3(4)%wc = 0.0
			Do j = 1, numA+1
				WF3(4)%R(j) = 0.0
				WF3(4)%x(j) = 0.0
				WF3(4)%y(j) = 0.0	
			Enddo
			WF3(4)%kint = 0.0
			WF3(4)%a = 0.0
			WF3(4)%b = 0.0
			WF3(4)%vol = 0.0		
			
		EndIf
		
		
		call calcCellRA
!		call calcWFRootAbs
		
		
		IF (RedistStatus.EQ.10) THEN
			n = 1
			WF3(1)%WC = sd%wcs
			romin = Q*0.0265*(sd%H**(-0.6098))*(sd%Ksat**(-0.6555))
!			romin = 0.9327*Q*(sd%H**-1.0)*(sd%Ksat**-1.0)/(2*pi)
			WF3(1)%TInfil = WF3(1)%TInfil + deltaT
			Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(1)%y(numA+1))
			Call calcSA(n,WF3(1)%y(numA+1))
!!!!!!!!!!!!!!calc root abs!!!!!!!!!!!!!!!!!!!!!!!1			
			Call calcAvgWC(11)
            Call calcWFRootAbs
            WF3(1)%thetaRA = WF3(1)%RootAbs/WF3(1)%Vol
	      WF3(1)%WC = WF3(1)%WC + WF3(1)%ThetaRA
	      RootHold = 0.0
	      If (WF3(1)%WC.LT.sd%wcmin) Then
	          RootHold = (WF3(1)%WC - sd%wcmin)*WF3(1)%Vol
	          WF3(1)%WC = sd%wcmin
	      Endif
            If(WF3(1)%WC.EQ.sd%wcmin) THEN
                WF3(1)%RootAbsTotal = WF3(1)%RootAbsTotal+ 
     +                                  WF3(1)%RootAbs - RootHold
	          WF3(1)%RootAbsActual = WF3(1)%RootAbs - RootHold
	      ELSE
	          WF3(1)%RootAbsTotal = WF3(1)%RootAbsTotal + WF3(1)%RootAbs
	          WF3(1)%RootAbsActual = WF3(1)%RootAbs
	      ENDIF
                       
			WF3(1)%FAmt = (WF3(1)%Vol)*(WF3(1)%WC-sd%wci)-WF3(1)%RootAbsTotal
!			WF3(1)%FAmt = (WF3(1)%Vol-2.0*pi*(romin**3.0)/3.0)*(WF3(1)%WC-sd%wci)
		ELSEIF (RedistStatus.EQ.20) THEN
			n = 2
			WF3(2)%TInfil = WF3(2)%TInfil + deltaT
!			Call calcR(n,WF3(2)%WC - WF3(1)%WC,
!     +			WF3(2)%TInfil)
!			Call calcXY(n)
!			Call calcEllipse(n)
!			Call calcVol(n,WF3(2)%y(numA+1))
!			Call calcSA(n,WF3(2)%y(numA+1))
!			WF3(1)%FAmt = WF3(1)%Vol*(WF3(1)%WC-sd%wci)
		ELSEIF (RedistStatus.EQ.30) THEN
			n = 3
			WF3(3)%TInfil = WF3(3)%TInfil + deltaT
			Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(3)%y(numA+1))
			Call calcSA(n,WF3(3)%y(numA+1))
		ELSEIF (RedistStatus.EQ.40) THEN
			n = 4
			WF3(4)%TInfil = WF3(4)%TInfil + deltaT
			Call calcR(n,WF3(4)%WC - WF3(3)%WC,WF3(4)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(4)%y(numA+1))
			Call calcSA(n,WF3(4)%y(numA+1))

		ENDIF		
		
		
		
		IF (redistStatus.EQ.11) THEN
        		CALL Redist11_3(ts,tend,deltaT)
 !       		WF3(1)%FAmt = WF3(1)%Vol*(WF3(1)%WC-sd%wci)
 !       		WF3(1)%FAmt = (WF3(1)%Vol-2.0*pi*(romin**3.0)/3.0)*
 !    +                                  (WF3(1)%WC-sd%wci)
		ELSEIF (redistStatus.EQ.20) THEN
        		CALL Redist20_3(ts,tend,deltaT,redistStatus)
		ELSEIF (redistStatus.EQ.21) THEN
        		CALL Redist21_3(ts,tend,deltaT,redistStatus)
		ELSEIF (redistStatus.EQ.30) THEN
        		CALL Redist30_3(ts,tend,deltaT,redistStatus)
		ELSEIF (redistStatus.EQ.31) THEN
        		CALL Redist31_3(ts,tend,deltaT,redistStatus)            
		ELSEIF (redistStatus.EQ.40) THEN
       		CALL Redist40_3(ts,tend,deltaT,redistStatus)		
		ENDIF

		call calcAvgwc(rediststatus)
		call calcNormalAngle
		call calcFlowRates(deltaT)
		
		cummF = cummF + Q*deltat
		
	END


	Subroutine CalcR(n,dM,t)
c-----calculates wetting front radius during infiltration -------------------------------------------
	USE ModuleDefs3
	Implicit None
	
	Integer i,redistStatus,n
	Double Precision  t, dgdr, g, rold1,dM
	Double Precision rold(numA), error, tolerance
	
	tolerance = 0.000006d0
		
	Do i = 1,numA
		error = 1.0
		Rold(i) = WF3(n)%R(i+1) + tolerance
		g = 0.5d0*(Rold(i)**2.d0 - romin**2.d0) - 
     +		(romin + (sd%H/sin(angles(i))))*(Rold(i) - romin) + 
     +		(sd%H/sin(angles(i)))*(romin + sd%H/sin(angles(i)))* 
     +		log((Rold(i)*sin(angles(i)) + sd%H)/(romin*sin(angles(i))  
     +        + sd%H))- (sd%Ksat*romin*sin(angles(i))*t)/dM

		DO WHILE (error.GT.tolerance)
			dgdr = Rold(i) - romin - (sd%H/sin(angles(i))) + 
     +			sd%H*(romin + (sd%H/sin(angles(i))))/
     +			(Rold(i)*sin(angles(i)) + sd%H)
			WF3(n)%R(i+1) = Rold(i) - g/dgdr
			Rold(i) = WF3(n)%R(i+1)
			g = 0.5d0*(Rold(i)**2.d0 - romin**2.d0) - 
     +			(romin + (sd%H/sin(angles(i))))*(Rold(i) - romin) + 
     +			(sd%H/sin(angles(i)))*(romin + sd%H/sin(angles(i)))* 
     +			log((Rold(i)*sin(angles(i)) + sd%H)/(romin* 
     +			sin(angles(i))+sd%H)) - (sd%Ksat*romin*
     +			sin(angles(i))*t)/dM
			error = abs(g)
		ENDDO
	EndDo

c----------------------------------- Calculate Ground Surface radius ----------------------------------
	error = 1.d0
	Rold1 = WF3(n)%R(1) + tolerance
	g = (Rold1**3.d0 - romin**3.d0)/3.d0 - romin*(Rold1**2.d0 -  
     +	 romin**2.d0)/2.d0- (sd%Ksat*romin*t*sd%H)/dM

	DO WHILE (error.GT.tolerance)
		dgdr = Rold1**2.d0 - romin*Rold1
		WF3(n)%R(1) = Rold1 - g/dgdr
		Rold1 = WF3(n)%R(1)
		g = (Rold1**3.d0 - romin**3.d0)/3.d0 - romin*(Rold1**2.d0 - 
     +		romin**2.d0)/2.d0 - (sd%Ksat*romin*t*sd%H)/dM 
		error = abs(g)
	ENDDO
	END

	Subroutine CalcXY(n)
c-----Calculates XY coordinates of radii ---------------------------------------------------------------
	USE ModuleDefs3
	Implicit None
	
	Integer i,try,n
	Double Precision Rhold(numA),R0

	WF3(n)%x(1) = WF3(n)%R(1)
	WF3(n)%y(1) = 0.d0
	Do i = 2,numA
	 WF3(n)%x(i) = WF3(n)%R(i)*cos(angles(i-1))
	 WF3(n)%y(i) = WF3(n)%R(i)*sin(angles(i-1))
	EndDo
	WF3(n)%x(numA+1) = 0.0
	WF3(n)%y(numA+1) = WF3(n)%R(numA+1)
	END
		
		
	Subroutine CalcEllipse(n)
c-----Calcs wetting front ellipse parameters --------------------------------------------------------
	USE ModuleDefs3
	Implicit none
	
	Integer i, n
	Double Precision  xmax,ymax
		
c--------------------------	find max x to locate center -----------------------------------------
	xmax = WF3(n)%x(1)
	ymax = WF3(n)%y(1)
	Do i = 2,NumA+1
		If (WF3(n)%x(i).GT.xmax) THEN
			xmax = WF3(n)%x(i)
			ymax = WF3(n)%y(i)
		EndIF
	ENDdo
	WF3(n)%kint = ymax
	WF3(n)%a = xmax
	WF3(n)%b = WF3(n)%y(numA+1) - ymax
	end	


	Subroutine CalcVol(n,R)
c-----Calcs vol of wet bulb --------------------------------------------------------------------------
	USE ModuleDefs3
	Implicit none
	
	Double Precision p1,p2,p3,R
	Integer n

	p1 = WF3(n)%kint*((R/WF3(n)%b)**2.d0)
	p2 = (WF3(n)%kint**2.d0)*R/(WF3(n)%b**2.d0)
	p3 = (R**3.d0)/(3.d0*(WF3(n)%b**2.d0))
	WF3(n)%vol = pi*(WF3(n)%a**2.d0)*(R+p1-p2-p3)
	END

	Subroutine CalcSA(n,yint)
c-----Calcs SA of wet bulb ------------------------------------------------------------------------------
	USE ModuleDefs3
	Implicit none
	
	Integer i, npoint,n
	Parameter (npoint = 5)
	double precision yint
	double precision x(npoint),w(npoint),SAfunc
	
	call gauleg(0.d0,yint,x,w,npoint)
	WF3(n)%SA = 0.d0
	do i = 1,npoint
		WF3(n)%SA = WF3(n)%SA + w(i)*SAfunc(n,x(i))
	enddo
	End

	Real Function SAFunc(n,x)
c-----used in CalcSA subroutine ----------------------------------------------------------------
	USE ModuleDefs3
	Implicit None

	double precision x, f, g
	integer n
	
	f = ((WF3(n)%a**2.d0)-((WF3(n)%a/WF3(n)%b)**2.d0)*
     +		(x-WF3(n)%kint)**2.d0)**(0.5d0)
	g = -((WF3(n)%a/WF3(n)%b)**2.d0)*(x-WF3(n)%kint)*
     +	((WF3(n)%a**2.d0)-((WF3(n)%a/WF3(n)%b)**2.d0)*
     +	(x-WF3(n)%kint)**2.d0)**(-0.5d0)

	SAfunc = 2.d0 * pi*f*(1.0+(g**2.d0))**(0.5d0)
	End

	Subroutine gauleg(x1,x2,x,w,n)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc	
c     used in CalcSA subroutine															c
c	Calculates n-point Gauss Legendre abscissas and weights                             c
c	Adapted from: Press, W.H., B.P. Flannery, S.A. Teukolsky, W.T. Vetterling           c 
c	Numerical Recipes in Fortran: The Art of Scientific Computing                       c
c	(Cambridge University Press, 1996)                                                  c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	Implicit None
	
	Integer n,m,j,i
	double precision eps, z1,z,xm,xl,pp,p3,p2,p1,x2,x1
	double precision x(n), w(n)
	Parameter (eps = 3.0e-11)

	m = (n+1)/2.d0
	xm = 0.5d0*(x2+x1)
	xl = 0.5d0*(x2-x1)
	do i = 1, m
		z = cos (3.141592654d0*(i -0.25d0)/(n + 0.5d0))
		do while (abs(z-z1).ge.eps)
			p1 = 1.d0
			p2 = 0.d0
			do j = 1, n
				p3 = p2
				p2 = p1
				p1 = ((2.d0*j -1.d0)*z*p2-(j-1.d0)*p3)/j
			enddo
			pp = n*(z*p1-p2)/(z*z-1.d0)
			z1 = z
			z = z1 - p1/pp
		enddo
		x(i) = xm-xl*z
		x(n+1-i) = xm+xl*z
		w(i) = 2.d0*xl/((1.d0-z*z)*pp*pp)
		w(n+1-i) = w(i)
	enddo
	End


	Subroutine Derivs_3 (vol,SA,wc,y,dydx)
c----- used to calc decrease in WC in redistXX subroutines -------------------------------------------------
	USE ModuleDefs3
	Implicit None
	
	Double Precision RelSatI, RelSatO, KWCI, KWCO,GWC
	Double Precision expon1, expon2 
	Double Precision RAvg,	 y,dydx,wc
	Double Precision SA,Vol
    
	RAvg = (3.d0*vol/(2.d0*pi))**(1.d0/3.d0)

	RelSatI = (wc - sd%wcr)/(sd%wcs - sd%wcr)
	RelSatO = (y - sd%wcr)/(sd%wcs - sd%wcr)
  	expon1 = 3.0 + (1.0/sd%lambda);
  	expon2 = 3.0 + (2.0/sd%lambda)
	
	IF (relSatO.LT.0.d0) THEN
		GWC = sd%H*((-(abs(relSatO))**expon1) + (relSatI**expon1))/
     +		(1-(relSatI**expon1))
  		KWCO = -sd%Ksat*((abs(relSatO))**expon2)
	ELSE
		GWC = sd%H*((relSatO**expon1) - (relSatI**expon1))/
     +		(1-(relSatI**expon1))
  		KWCO = sd%Ksat*(relSatO**expon2)
	EndIf
	KWCI = sd%Ksat*(relSatI**expon2)
  	dydx =(-SA*KWCI - SA*KWCO - (sd%Ksat*Ravg*GWC))
     +			/(vol) 
	END

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	Subroutines RK4_3, RKQC_3 and OdeInt_3 perform perform 4th order Runge Kutta procedure c
c	Adapted from: Press, W.H., B.P. Flannery, S.A. Teukolsky, W.T. Vetterling              c
c	Numerical Recipes in Fortran: The Art of Scientific Computing                          c
c	(Cambridge University Press, 1996)                                                     c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	
	Subroutine RK4_3(vol,SA,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 vol,SA
	
	hh = h/2.0
	yt = y + hh*dydx
	CALL Derivs_3(vol,SA,WC,yt,dyt)
	yt = y + hh*dyt
  	CALL Derivs_3(vol,SA,WC,yt,dym)
  	yt = y + h*dym
  	dym = dym + dyt
  	CALL Derivs_3(vol,SA,WC, yt, dyt)
  	yout = y + h*(dydx + dyt + 2.0 * dym )/6.0

	END

	Subroutine RKQC_3(vol,SA,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, vol, SA
	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_3(vol,SA,WC,ysav,dysav,hh,ytemp)
    		x = xsav + hh
    		CALL Derivs_3(vol,SA,WC,ytemp,dydx)
    		CALL RK4_3(vol,SA,WC,ytemp,dydx,hh,y)
    		x = xsav + h
    	!test to see if x = xsav
    		CALL RK4_3(vol,SA,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_3(vol,SA,WC,ystart,x1, x2,eps,h1)
	Implicit None
	
	Double Precision x,hnext,tiny,hdid,h,yscal,y,dydx
	Double Precision vol, SA,wc,ystart,x1,x2,eps,h1,RA
	Integer MaxStp,nstp
  	
	MaxStp = 10000
  	tiny = 1.0e-3
  	x = x1
  	h = abs(h1)
  	y = ystart
  	DO nstp = 1 , maxstp 
  	  	CALL Derivs_3(vol,SA,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_3(vol,SA,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

	
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	Procedures Redist11_3, Redist20_3, Redist21_3, Redist30_3, Redist31_3 and Redist40_3 perform            c
c	the redistribution process depending on what redistribution status the system is in         c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	
	Subroutine Redist11_3(ts,te,dt)
c----- redistribution routine if one front --------------------------------------------------------
	USE ModuleDefs3
	Implicit None 
	Double Precision ts, te, dt,volWi1,volF1,volW2
	Double Precision  y,dydx,wc1Hold,vol1Hold,thetao
      Double Precision VolWi, VolW,RootHold
	Integer n
!!!!first calculate the vol of water want to maintain for mass conservation
      volWi1 = WF3(1)%Vol*(sd%wci)
      vol1Hold = WF3(1)%Vol
      volF1 = WF3(1)%Vol*(WF3(1)%WC-sd%wci)
!!!!!Then redistribute!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	y = WF3(1)%WC
	CALL OdeInt_3(WF3(1)%vol,WF3(1)%SA,sd%wci,y,ts,te,
     +                      0.00001d0,dT)
	WF3(1)%WC = y !+ dydx
	WF3(1)%WC = y !+ dydx ! - QP^[Q]^.adjFactor
	IF  (WF3(1)%WC.LT.sd%wcmin) THEN	
		WF3(1)%WC = sd%wcmin
		IF (WF3(1)%WC.LT.sd%wcmin) THEN
			WF3(1)%WC = sd%wcmin
		ENDIF
	ENDIF
	
!!!!!!!!!calculate equivalent vol for mass balance/allow for expansion from redist!!!!!!!!!!
	WF3(1)%Vol = (volWi1+VolF1-sd%wci*vol1Hold)/(WF3(1)%WC-sd%wci)
!!!!!!!!calculate new Radii, a and b based on new Vol!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      volZ = WF3(1)%vol
      xk = WF3(1)%kint
      xe = WF3(1)%a/WF3(1)%b
      Call volZ_zbrent(WF3(1)%R(1),WF3(1)%R(numA+1))
      WF3(1)%a = xaZ
      WF3(1)%b = xbZ 
   
!!!!!!!	calculate new SA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      n=1
      Call calcXY(n)
      Call calcSA(n,WF3(1)%y(numA+1))

!	Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
!	Call calcXY(n)
!	Call calcEllipse(n)
!	Call calcVol(n,WF3(1)%y(numA+1))
!	Call calcSA(n,WF3(1)%y(numA+1))

!!!!!!!!!calculate change of water content from root abstraction!!!!!!!!!!!!
!!!!!!!!!vol and shape should remain the same and wc decrease only !!!!!!!!!!!
      Call calcAvgWC(11)
      Call calcWFRootAbs
      
      WF3(1)%thetaRA = WF3(1)%RootAbs/WF3(1)%Vol
	WF3(1)%WC = WF3(1)%WC + WF3(1)%ThetaRA
	RootHold = 0.0
	If (WF3(1)%WC.LT.sd%wcmin) Then
	  RootHold = (WF3(1)%WC - sd%wcmin)*WF3(1)%Vol
	  WF3(1)%WC = sd%wcmin
	Endif
      If(WF3(1)%WC.EQ.sd%wcmin) THEN
        WF3(1)%RootAbsTotal = WF3(1)%RootAbsTotal+ WF3(1)%RootAbs - 
     +                                  RootHold
	  WF3(1)%RootAbsActual = WF3(1)%RootAbs - RootHold
	ELSE
	  WF3(1)%RootAbsTotal = WF3(1)%RootAbsTotal + WF3(1)%RootAbs
	  WF3(1)%RootAbsActual = WF3(1)%RootAbs
	ENDIF
      
      volW = WF3(1)%WC*WF3(1)%Vol
      volWi = WF3(1)%Vol*sd%wci
      WF3(1)%FAmt = volw - volWi-WF3(1)%RootAbsTotal
	END
	
	Subroutine Redist20_3(ts,te,dt,redistStatus)
c----- redist routine if two fronts & only one front redist -------------------------------------------------
	USE ModuleDefs3
	Implicit None

	Double Precision y,volWi1,Vol1Hold,VolF1,rootHold,VolW,VolWi
	Double Precision ts, te, dt
	Integer n, j,redistStatus

!!!!first calculate the vol of water want to maintain for mass conservation
      volWi1 = WF3(1)%Vol*(sd%wci)
      vol1Hold = WF3(1)%Vol
      volF1 = WF3(1)%Vol*(WF3(1)%WC-sd%wci)
!!!!!Then redistribute!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	y = WF3(1)%WC
	CALL OdeInt_3(WF3(1)%vol,WF3(1)%SA,sd%wci,y,ts,te,
     +                      0.00001d0,dT)
	WF3(1)%WC = y !+ dydx
	WF3(1)%WC = y !+ dydx ! - QP^[Q]^.adjFactor
	IF  (WF3(1)%WC.LT.sd%wcmin) THEN	
		WF3(1)%WC = sd%wcmin
		IF (WF3(1)%WC.LT.sd%wcmin) THEN
			WF3(1)%WC = sd%wcmin
		ENDIF
	ENDIF
	
!!!!!!!!!calculate equivalent vol for mass balance/allow for expansion from redist!!!!!!!!!!
	WF3(1)%Vol = (volWi1+VolF1-sd%wci*vol1Hold)/(WF3(1)%WC-sd%wci)
!!!!!!!!calculate new Radii, a and b based on new Vol!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      volZ = WF3(1)%vol
      xk = WF3(1)%kint
      xe = WF3(1)%a/WF3(1)%b
      Call volZ_zbrent(WF3(1)%R(1),WF3(1)%R(numA+1))
      WF3(1)%a = xaZ
      WF3(1)%b = xbZ 
   
!!!!!!!	calculate new SA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      n=1
      Call calcXY(n)
      Call calcSA(n,WF3(1)%y(numA+1))

!!!!!!!!!calculate change of water content from root abstraction!!!!!!!!!!!!
!!!!!!!!!vol and shape should remain the same and wc decrease only !!!!!!!!!!!
      Call calcAvgWC(11)
      Call calcWFRootAbs
      
      WF3(1)%thetaRA = WF3(1)%RootAbs/WF3(1)%Vol
	WF3(1)%WC = WF3(1)%WC + WF3(1)%ThetaRA
	RootHold = 0.0
	If (WF3(1)%WC.LT.sd%wcmin) Then
	  RootHold = (WF3(1)%WC - sd%wcmin)*WF3(1)%Vol
	  WF3(1)%WC = sd%wcmin
	Endif
      If(WF3(1)%WC.EQ.sd%wcmin) THEN
        WF3(1)%RootAbsTotal = WF3(1)%RootAbsTotal+ WF3(1)%RootAbs - 
     +                                  RootHold
	  WF3(1)%RootAbsActual = WF3(1)%RootAbs - RootHold
	ELSE
	  WF3(1)%RootAbsTotal = WF3(1)%RootAbsTotal + WF3(1)%RootAbs
	  WF3(1)%RootAbsActual = WF3(1)%RootAbs
	ENDIF
      
      volW = WF3(1)%WC*WF3(1)%Vol
      volWi = WF3(1)%Vol*sd%wci
      WF3(1)%FAmt = volw - volWi-WF3(1)%RootAbsTotal
  	
  	
  	! WF3 2 does not redist since it is infiltrating
  	n = 2
	WF3(2)%WC = sd%wcs
	Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(2)%y(numA+1))
	Call calcSA(n,WF3(2)%y(numA+1))
	!!!!!!!!!!!!!!calc root abs!!!!!!!!!!!!!!!!!!!!!!!1			
			Call calcAvgWC(20)
            Call calcWFRootAbs
            WF3(2)%thetaRA = WF3(2)%RootAbs/WF3(2)%Vol
	      WF3(2)%WC = WF3(2)%WC + WF3(2)%ThetaRA
	      RootHold = 0.0
	      If (WF3(2)%WC.LT.sd%wcmin) Then
	          RootHold = (WF3(2)%WC - sd%wcmin)*WF3(2)%Vol
	          WF3(2)%WC = sd%wcmin
	      Endif
            If(WF3(2)%WC.EQ.sd%wcmin) THEN
                WF3(2)%RootAbsTotal = WF3(2)%RootAbsTotal+ 
     +                                  WF3(2)%RootAbs - RootHold
	          WF3(2)%RootAbsActual = WF3(2)%RootAbs - RootHold
	      ELSE
	          WF3(2)%RootAbsTotal = WF3(2)%RootAbsTotal + WF3(2)%RootAbs
	          WF3(2)%RootAbsActual = WF3(2)%RootAbs
	      ENDIF
                       
			WF3(2)%FAmt = (WF3(2)%Vol)*(WF3(2)%WC-sd%wci)-WF3(2)%RootAbsTotal
!!!!!!!!ToDo : Need to add updates if they merge considering the new RA!!!!!!!!!!!!!!!!!!!!!
	!check conditions
  	IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN
  	!2 WF3 to 1 WF3
    		WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    		WF3(1)%WC = WF3(2)%WC
    		WF3(1)%WC = WF3(2)%WC
    		WF3(1)%RootAbsTotal = WF3(1)%RootAbsTotal + WF3(2)%RootAbsTotal
    		WF3(1)%RootAbsActual = WF3(1)%RootAbsActual + WF3(2)%RootAbsActual
    		! no need to update root abs
    		n = 1
		Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(1)%y(numA+1))
		Call calcSA(n,WF3(1)%y(numA+1))	
		WF3(2)%wc = 0.0
		Do j = 1, numA+1
			WF3(2)%R(j) = 0.0
			WF3(2)%x(j) = 0.0
			WF3(2)%y(j) = 0.0	
          Enddo
		WF3(2)%kint = 0.0
		WF3(2)%a = 0.0
		WF3(2)%b = 0.0
		WF3(2)%vol = 0.0
		WF3(2)%SA = 0.0
		WF3(2)%TInfil = 0.0
		WF3(2)%FAmt = 0.0
		WF3(2)%RootAbsTotal = 0.0
		WF3(2)%RootAbsActual = 0.0
		redistStatus = 10
  	ENDIF !if z2 > z1

	
	
	
	
	END
	
	Subroutine Redist21_3(ts,te,dt,redistStatus)
c----- redist routine if two fronts & both fronts redist --------------------------------------------
	USE ModuleDefs3
	Implicit None

	Double Precision y
	Double Precision ts, te, dt
	Integer redistStatus,n,j
  	
	!!decrease wc1 first
  	y = WF3(1)%WC
  	CALL OdeInt_3(WF3(1)%vol,WF3(1)%SA,sd%wci,y,ts,te,
     +                      0.00001d0,dT)
  	WF3(1)%WC = y
  	WF3(1)%WC = y ! - QP^[Q]^.adjFactor;
  	IF  (WF3(1)%WC.LT.sd%wci) THEN
  		WF3(1)%WC = sd%wcmin
    		IF (WF3(1)%WC.LT.sd%wcmin)THEN
       		WF3(1)%WC = sd%wcmin
		ENDIF
  	ENDIF
  	n = 1
	Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(1)%y(numA+1))
	Call calcSA(n,WF3(1)%y(numA+1))
  	!!//decrease wc2 second
  	y = WF3(2)%WC
  	CALL OdeInt_3(WF3(2)%vol,WF3(2)%SA,WF3(1)%WC,y,ts,te,
     +                      0.00001d0,dT)
  	WF3(2)%WC = y
  	
  	n = 2
	Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(2)%y(numA+1))
	Call calcSA(n,WF3(2)%y(numA+1))
	IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC.OR.
     +	WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN !! //merge 1 an 2
 	 	!begin two fronts to one
    		IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC)THEN
          	!BEGIN
      		WF3(1)%WC = WF3(1)%WC ! //WF3 one predominates
      		IF (WF3(2)%WC.LT.sd%wcmin) THEN
         			WF3(2)%WC = sd%wcmin
			ENDIF
    		ELSE
    		!BEGIN
     			WF3(1)%WC = WF3(2)%WC !WF3 two predominates
      		WF3(1)%WC = WF3(2)%WC
    		ENDIF
    		WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    		! no need to update root abs
    		n = 1
		Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(1)%y(numA+1))
		Call calcSA(n,WF3(1)%y(numA+1))	
		WF3(2)%wc = 0.0
		Do j = 1, numA+1
			WF3(2)%R(j) = 0.0
			WF3(2)%x(j) = 0.0
			WF3(2)%y(j) = 0.0	
          Enddo
		WF3(2)%kint = 0.0
		WF3(2)%a = 0.0
		WF3(2)%b = 0.0
		WF3(2)%vol = 0.0
		WF3(2)%SA = 0.0
		WF3(2)%TInfil = 0.0
		redistStatus = 11
  	! ENDIF wc2 < wcmin
  	ENDIF
	END

	Subroutine Redist30_3(ts,te,dt,redistStatus)
c----- redist routine if three fronts & only two fronts redist ---------------------------------------------
	USE ModuleDefs3
	Implicit None

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

	!!decrease wc1 first
  	y = WF3(1)%WC 
  	CALL OdeInt_3(WF3(1)%vol,WF3(1)%SA,sd%wci,y,ts,te,
     +                      0.00001d0,dT)
  	WF3(1)%WC = y
  	WF3(1)%WC = y ! - QP^[Q]^.adjFactor;
  	IF  (WF3(1)%WC.LT.sd%wci) THEN
  	 	 WF3(1)%WC = sd%wcmin
    		 IF (WF3(1)%WC.LT.sd%wcmin)THEN
       		WF3(1)%WC = sd%wcmin
		ENDIF
  	ENDIF
  	n = 1
	Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(1)%y(numA+1))
	Call calcSA(n,WF3(1)%y(numA+1))
 	!decrease wc2 second
  	!!//decrease wc2 second
  	y = WF3(2)%WC
  	CALL OdeInt_3(WF3(2)%vol,WF3(2)%SA,WF3(1)%WC,y,ts,te,
     +                      0.00001d0,dT)
  	WF3(2)%WC = y
  	
  	n = 2
	Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(2)%y(numA+1))
	Call calcSA(n,WF3(2)%y(numA+1))
  	!!checks
  	IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC.OR.
     +	WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN !!merge 1 an 2
  	!! 3 fronts down to 2 fronts
    		IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC)THEN
        !BEGIN
      		WF3(1)%WC = WF3(1)%WC! //WF3 one predominates
      		IF (WF3(2)%WC.LT.sd%wcmin) THEN
         			WF3(2)%WC = sd%wcmin
			ENDIF
    		ELSE
    		!BEGIN
     			WF3(1)%WC = WF3(2)%WC !WF3 two predominates
      		WF3(1)%WC = WF3(2)%WC
    		ENDIF
    		WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    		! no need to update root abs
    		n = 1
		Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(1)%y(numA+1))
		Call calcSA(n,WF3(1)%y(numA+1))	
		redistStatus = 20
    		WF3(2)%TInfil = WF3(3)%TInfil 
		WF3(2)%WC = WF3(3)%WC 
		WF3(2)%WC = WF3(3)%WC
		n = 2
		Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(2)%y(numA+1))
		Call calcSA(n,WF3(2)%y(numA+1))	 
 		WF3(3)%WC = 0.0 
		WF3(3)%WC = 0.0 
		Do j = 1, numA+1
			WF3(3)%R(j) = 0.0
			WF3(3)%x(j) = 0.0
			WF3(3)%y(j) = 0.0	
          Enddo
		WF3(3)%kint = 0.0
		WF3(3)%a = 0.0
		WF3(3)%b = 0.0
		WF3(3)%vol = 0.0
		WF3(3)%SA = 0.0
		WF3(3)%TInfil = 0.0
    	!!now move WF3 2 down which was WF3 3 and is only infil, not redist
    		n = 2
		Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(2)%y(numA+1))
		Call calcSA(n,WF3(2)%y(numA+1))
    	!!check
    		IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN
    	!!BEGIN //two fronts to one
      		WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    			WF3(1)%WC = WF3(2)%WC
			n = 1
			Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(1)%y(numA+1))
			Call calcSA(n,WF3(1)%y(numA+1))	
    			WF3(2)%wc = 0.0
			Do j = 1, numA+1
				WF3(2)%R(j) = 0.0
				WF3(2)%x(j) = 0.0
				WF3(2)%y(j) = 0.0	
			Enddo
			WF3(2)%kint = 0.0
			WF3(2)%a = 0.0
			WF3(2)%b = 0.0
			WF3(2)%vol = 0.0
			WF3(2)%SA = 0.0
			WF3(2)%TInfil = 0.0 
      		redistStatus = 10
     		ENDIF
  	!!END//wc2 < wcmin,z2>z1
  	ELSE !!//still 3 fronts
  	!BEGIN
    			n = 3 !just move down since only infil
			Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(3)%y(numA+1))
			Call calcSA(n,WF3(3)%y(numA+1)) 
	    	!!checks
    		IF (WF3(3)%R(numA+1).GE.WF3(2)%R(numA+1)) THEN
    	!!BEGIN //three fronts to two
      		WF3(2)%TInfil = WF3(2)%TInfil + WF3(3)%TInfil
    			n = 2
			Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(2)%y(numA+1))
			Call calcSA(n,WF3(2)%y(numA+1))	
    			WF3(3)%wc = 0.0
			Do j = 1, numA+1
				WF3(3)%R(j) = 0.0
				WF3(3)%x(j) = 0.0
				WF3(3)%y(j) = 0.0	
			Enddo
			WF3(3)%kint = 0.0
			WF3(3)%a = 0.0
			WF3(3)%b = 0.0
			WF3(3)%vol = 0.0
			WF3(3)%SA = 0.0
			WF3(3)%TInfil = 0.0 
      		redistStatus = 20		
			IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN	!!recheck 
				WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    				n = 1
				Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(1)%y(numA+1))
				Call calcSA(n,WF3(1)%y(numA+1))	
    	    				WF3(2)%wc = 0.0
				Do j = 1, numA+1
					WF3(2)%R(j) = 0.0
					WF3(2)%x(j) = 0.0
					WF3(2)%y(j) = 0.0	
				Enddo
				WF3(2)%kint = 0.0
				WF3(2)%a = 0.0
				WF3(2)%b = 0.0
				WF3(2)%vol = 0.0
				WF3(2)%SA = 0.0
				WF3(2)%TInfil = 0.0 
      			redistStatus = 10
			ENDIF
	!!endif 3 > 2
		ENDIF
	ENDIF !!else still 3 fronts
	END

	Subroutine Redist31_3(ts,te,dt,redistStatus)
c----- redist routine if three fronts & all fronts redist ------------------------------------------------------
	USE ModuleDefs3
	Implicit None

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

	!!decrease wc1 first
  	y = WF3(1)%WC 
  	CALL OdeInt_3(WF3(1)%vol,WF3(1)%SA,sd%wci,y,ts,te,
     +                      0.00001d0,dT)
  	WF3(1)%WC = y
  	
  	IF  (WF3(1)%WC.LT.sd%wci) THEN
  	 	 WF3(1)%WC = sd%wcmin
    		 IF (WF3(1)%WC.LT.sd%wcmin)THEN
       		WF3(1)%WC = sd%wcmin
		 ENDIF
  	ENDIF
  	n = 1
	Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(1)%y(numA+1))
	Call calcSA(n,WF3(1)%y(numA+1))
  	!decrease wc2 second
  	!!//decrease wc2 second
  	y = WF3(2)%WC
  	CALL OdeInt_3(WF3(2)%vol,WF3(2)%SA,WF3(1)%WC,y,ts,te,
     +                      0.00001d0,dT)
  	WF3(2)%WC = y
  	
  	n = 2
	Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(2)%y(numA+1))
	Call calcSA(n,WF3(2)%y(numA+1))
     	!!checks
  	IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC.OR.
     +	WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN !!merge 1 an 2
  	!! 3 fronts down to 2 fronts
    		IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC)THEN
        !BEGIN
      		WF3(1)%WC = WF3(1)%WC ! //WF3 one predominates
      		IF (WF3(2)%WC.LT.sd%wcmin) THEN
         			WF3(2)%WC = sd%wcmin
			ENDIF
    		ELSE
    		!BEGIN
     			WF3(1)%WC = WF3(2)%WC !WF3 two predominates
      		WF3(1)%WC = WF3(2)%WC
      		
    		ENDIF
    		WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    		WF3(1)%RootAbs = WF3(2)%RootAbs
    		n = 1
		Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(1)%y(numA+1))
		Call calcSA(n,WF3(1)%y(numA+1))	
    		redistStatus = 21
    		WF3(2)%TInfil = WF3(3)%TInfil 
		WF3(2)%WC = WF3(3)%WC 
		WF3(2)%WC = WF3(3)%WC
		WF3(2)%RootAbs = WF3(3)%RootAbs
		n = 2
		Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(2)%y(numA+1))
		Call calcSA(n,WF3(2)%y(numA+1))	 
    		WF3(3)%WC = 0.0 
		WF3(3)%WC = 0.0 
		Do j = 1, numA+1
			WF3(3)%R(j) = 0.0
			WF3(3)%x(j) = 0.0
			WF3(3)%y(j) = 0.0	
          Enddo
		WF3(3)%kint = 0.0
		WF3(3)%a = 0.0
		WF3(3)%b = 0.0
		WF3(3)%vol = 0.0
		WF3(3)%SA = 0.0
		WF3(3)%TInfil = 0.0
		WF3(3)%RootAbs = 0.0
	   	!!//need to decrease wc2 which was wc3
    		y = WF3(2)%WC
  		CALL OdeInt_3(WF3(2)%vol,WF3(2)%SA,WF3(1)%WC,y,ts,te,
     +						0.00001d0,dT)
  		WF3(2)%WC = y
  	    
		n = 2
		Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(2)%y(numA+1))
		Call calcSA(n,WF3(2)%y(numA+1))
    		!!checks
  		IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC.OR.
     +		WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN !!merge 1 an 2
  		!! 2 fronts down to 1 fronts
    			IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC)THEN
        		!BEGIN
      			WF3(1)%WC = WF3(1)%WC! //WF3 one predominates
      			IF (WF3(2)%WC.LT.sd%wcmin) THEN
         				WF3(2)%WC = sd%wcmin
            		ENDIF
    			ELSE
    		!BEGIN
     				WF3(1)%WC = WF3(2)%WC !WF3 two predominates
      			WF3(1)%WC = WF3(2)%WC
    			ENDIF
    				WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    				! no need to update root abs
    				n = 1
				Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(1)%y(numA+1))
				Call calcSA(n,WF3(1)%y(numA+1))	
    				WF3(2)%wc = 0.0
				Do j = 1, numA+1
					WF3(2)%R(j) = 0.0
					WF3(2)%x(j) = 0.0
					WF3(2)%y(j) = 0.0	
				Enddo
				WF3(2)%kint = 0.0
				WF3(2)%a = 0.0
				WF3(2)%b = 0.0
				WF3(2)%vol = 0.0
				WF3(2)%SA = 0.0
				WF3(2)%TInfil = 0.0 
      			redistStatus = 11
		!!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 WF3 for now
  	!!BEGIN
    	!!//decrease wc3 third
    		y = WF3(3)%WC
  		CALL OdeInt_3(WF3(3)%vol,WF3(3)%SA,WF3(2)%WC,y,ts,te,
     +							0.00001d0,dT)
  		WF3(3)%WC = y
  	
  		n = 3
		Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(3)%y(numA+1))
		Call calcSA(n,WF3(3)%y(numA+1))
    	 		IF  (WF3(3)%WC.LE.WF3(2)%WC.OR.WF3(3)%R(numA+1).GE.
     +		 WF3(2)%R(numA+1)) THEN !!  //merge 3 and 2
    	!!BEGIN  // 3 fronts to 2
      		IF (WF3(3)%WC.LE.WF3(2)%WC) THEN 
		    	WF3(2)%WC = WF3(2)%WC
      		ELSE
      		!!BEGIN
        			WF3(2)%WC = WF3(3)%WC 
        			WF3(2)%WC = WF3(3)%WC 
      		ENDIF
      		WF3(2)%TInfil = WF3(2)%TInfil + WF3(3)%TInfil
    			! no need to update root abs
    			n = 2
			Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(2)%y(numA+1))
			Call calcSA(n,WF3(2)%y(numA+1))	
     			WF3(3)%wc = 0.0
			Do j = 1, numA+1
				WF3(3)%R(j) = 0.0
				WF3(3)%x(j) = 0.0
				WF3(3)%y(j) = 0.0	
			Enddo
			WF3(3)%kint = 0.0
			WF3(3)%a = 0.0
			WF3(3)%b = 0.0
			WF3(3)%vol = 0.0
			WF3(3)%SA = 0.0
			WF3(3)%TInfil = 0.0 
      		redistStatus = 21		    		
	    		IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN	!!recheck 
				WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    				n = 1
				Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(1)%y(numA+1))
				Call calcSA(n,WF3(1)%y(numA+1))	
     				WF3(2)%wc = 0.0
				Do j = 1, numA+1
					WF3(2)%R(j) = 0.0
					WF3(2)%x(j) = 0.0
					WF3(2)%y(j) = 0.0	
				Enddo
				WF3(2)%kint = 0.0
				WF3(2)%a = 0.0
				WF3(2)%b = 0.0
				WF3(2)%vol = 0.0
				WF3(2)%SA = 0.0
				WF3(2)%TInfil = 0.0 
      			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 WF3
	END

	Subroutine Redist40_3(ts,te,dt,redistStatus)
c----- redist routine if four fronts & only three fronts redist -----------------------------------------------------
	USE ModuleDefs3
	Implicit None

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

	!!decrease wc1 first
  	y = WF3(1)%WC 
  	CALL OdeInt_3(WF3(1)%vol,WF3(1)%SA,sd%wci,y,ts,te,
     +                              0.00001d0,dT)
  	WF3(1)%WC = y
  
  	IF  (WF3(1)%WC.LT.sd%wci) THEN
  	 	 WF3(1)%WC = sd%wcmin
    		 IF (WF3(1)%WC.LT.sd%wcmin)THEN
       		WF3(1)%WC = sd%wcmin
		 ENDIF
  	ENDIF
  	n = 1
	Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(1)%y(numA+1))
	Call calcSA(n,WF3(1)%y(numA+1))
  	
	!decrease wc2 second
  	!!//decrease wc2 second
  	y = WF3(2)%WC
  	CALL OdeInt_3(WF3(2)%vol,WF3(2)%SA,WF3(1)%WC,y,ts,te,
     +                          0.00001d0,dT)
  	WF3(2)%WC = y
  	
  	n = 2
	Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
	Call calcXY(n)
	Call calcEllipse(n)
	Call calcVol(n,WF3(2)%y(numA+1))
	Call calcSA(n,WF3(2)%y(numA+1))
  	!!checks
  	IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC.OR.
     +  	WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN !!merge 1 an 2
  	!! 4 fronts down to 3 fronts
    		IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC)THEN
		!BEGIN
      		WF3(1)%WC = WF3(1)%WC! //WF3 one predominates
!      		IF (WF3(2)%WC.LT.sd%wcmin) THEN
!        			WF3(2)%WC = sd%wcmin
!			ENDIF
    		ELSE
    		!BEGIN
     			WF3(1)%WC = WF3(2)%WC !WF3 two predominates
    		ENDIF
    		WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    		WF3(1)%RootAbs = WF3(2)%RootAbs
    		n = 1
		Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(1)%y(numA+1))
		Call calcSA(n,WF3(1)%y(numA+1))	
		redistStatus = 30
    		WF3(2)%TInfil = WF3(3)%TInfil 
		WF3(2)%WC = WF3(3)%WC 
		WF3(2)%WC = WF3(3)%WC
		WF3(2)%RootAbs = WF3(3)%RootAbs
		n = 2
		Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(2)%y(numA+1))
		Call calcSA(n,WF3(2)%y(numA+1))	 
 		WF3(3)%TInfil = WF3(4)%TInfil 
		WF3(3)%WC = WF3(4)%WC 
		WF3(3)%WC = WF3(4)%WC
		WF3(3)%rootAbs = 0.0
		!there is not root abstraction for WF4 since infil only
		n = 3
		Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(3)%y(numA+1))
		Call calcSA(n,WF3(3)%y(numA+1))
		WF3(4)%WC = 0.0 
		WF3(4)%WC = 0.0 
		Do j = 1, numA+1
			WF3(4)%R(j) = 0.0
			WF3(4)%x(j) = 0.0
			WF3(4)%y(j) = 0.0	
          Enddo
		WF3(4)%kint = 0.0
		WF3(4)%a = 0.0
		WF3(4)%b = 0.0
		WF3(4)%vol = 0.0
		WF3(4)%SA = 0.0
		WF3(4)%TInfil = 0.0
    	!!//need to decrease wc2 which was wc3
    		y = WF3(2)%WC
  		CALL OdeInt_3(WF3(2)%vol,WF3(2)%SA,WF3(1)%WC,y,ts,te,
     +								0.00001d0,dT)
  		WF3(2)%WC = y
  	
  		n = 2
		Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(2)%y(numA+1))
		Call calcSA(n,WF3(2)%y(numA+1))	 
   		!!checks
  		IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC.OR.
     +		WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN !!merge 1 an 2
  			!! 3 fronts down to 2 fronts
    			IF (WF3(2)%WC.LT.sd%wcmin.OR.WF3(2)%WC.LE.WF3(1)%WC)THEN
        		!BEGIN
      			WF3(1)%WC = WF3(1)%WC! //WF3 one predominates
!     			IF (WF3(2)%WC.LT.sd%wcmin) THEN
!         				WF3(2)%WC = sd%wcmin
!            		ENDIF
    			ELSE
    			!BEGIN
     				WF3(1)%WC = WF3(2)%WC !WF3 two predominates
      			WF3(1)%WC = WF3(2)%WC
    			ENDIF
    			
			WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    			WF3(1)%RootAbs = WF3(2)%RootAbs
    			n = 1
			Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(1)%y(numA+1))
			Call calcSA(n,WF3(1)%y(numA+1))	
			redistStatus = 20
    			WF3(2)%TInfil = WF3(3)%TInfil 
			WF3(2)%WC = WF3(3)%WC 
			WF3(2)%WC = WF3(3)%WC
			WF3(2)%RootAbs = WF3(3)%RootAbs
			n = 2
			Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(2)%y(numA+1))
			Call calcSA(n,WF3(2)%y(numA+1))	 
 			WF3(3)%WC = 0.0 
			WF3(3)%WC = 0.0 
			Do j = 1, numA+1
				WF3(3)%R(j) = 0.0
				WF3(3)%x(j) = 0.0
				WF3(3)%y(j) = 0.0	
			Enddo
			WF3(3)%kint = 0.0
			WF3(3)%a = 0.0
			WF3(3)%b = 0.0
			WF3(3)%vol = 0.0
			WF3(3)%SA = 0.0
			WF3(3)%TInfil = 0.0
		!!//move down z2; was originally 4th WF3 so do not decrease
      		n = 2
			Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(2)%y(numA+1))
			Call calcSA(n,WF3(2)%y(numA+1))	
      		!!//check last time
      		IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1))THEN 
      			!!BEGIN
    				WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    				WF3(1)%WC = WF3(2)%WC
    				!no need to update root abs
				n = 1
				Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(1)%y(numA+1))
				Call calcSA(n,WF3(1)%y(numA+1))	
    				WF3(2)%wc = 0.0
				Do j = 1, numA+1
					WF3(2)%R(j) = 0.0
					WF3(2)%x(j) = 0.0
					WF3(2)%y(j) = 0.0	
				Enddo
				WF3(2)%kint = 0.0
				WF3(2)%a = 0.0
				WF3(2)%b = 0.0
				WF3(2)%vol = 0.0
				WF3(2)%SA = 0.0
				WF3(2)%TInfil = 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 WF3 4; do not redist, only infil
      		n = 3
			Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(3)%y(numA+1))
			Call calcSA(n,WF3(3)%y(numA+1))
      		IF (WF3(3)%R(numA+1).GE.WF3(2)%R(numA+1)) THEN 
      			!!BEGIN //merge 3 and 2; WF31 stays the same
    				WF3(2)%WC = WF3(3)%WC
				WF3(2)%TInfil = WF3(2)%TInfil + WF3(3)%TInfil
    				!no need to update root abs
    				n = 2
    				Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(2)%y(numA+1))
				Call calcSA(n,WF3(2)%y(numA+1))	
    				WF3(3)%wc = 0.0
				Do j = 1, numA+1
					WF3(3)%R(j) = 0.0
					WF3(3)%x(j) = 0.0
					WF3(3)%y(j) = 0.0	
				Enddo
				WF3(3)%kint = 0.0
				WF3(3)%a = 0.0
				WF3(3)%b = 0.0
				WF3(3)%vol = 0.0
				WF3(3)%SA = 0.0
				WF3(3)%TInfil = 0.0 
      			redistStatus = 20		 
				IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN	!!recheck 
					WF3(1)%WC = WF3(2)%WC
					WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    					n = 1
					Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
					Call calcXY(n)
					Call calcEllipse(n)
					Call calcVol(n,WF3(1)%y(numA+1))
					Call calcSA(n,WF3(1)%y(numA+1))	
    		
    					WF3(2)%wc = 0.0
					Do j = 1, numA+1
						WF3(2)%R(j) = 0.0
						WF3(2)%x(j) = 0.0
						WF3(2)%y(j) = 0.0	
					Enddo
					WF3(2)%kint = 0.0
					WF3(2)%a = 0.0
					WF3(2)%b = 0.0
					WF3(2)%vol = 0.0
					WF3(2)%SA = 0.0
					WF3(2)%TInfil = 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 WF3 for now
  	!!BEGIN
    	!!//decrease wc3 third
    		y = WF3(3)%WC
    		CALL OdeInt_3(WF3(3)%vol,WF3(3)%SA,WF3(2)%WC,y,
     +							ts,te,0.00001d0,dT)
    		WF3(3)%WC = y
    	
    		n = 3
		Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
		Call calcXY(n)
		Call calcEllipse(n)
		Call calcVol(n,WF3(3)%y(numA+1))
		Call calcSA(n,WF3(3)%y(numA+1))
    		IF  (WF3(3)%WC.LE.WF3(2)%WC.OR.WF3(3)%R(numA+1).GE.
     +     		WF3(2)%R(numA+1)) THEN !merge 3 and 2
    		!!BEGIN  // 4 fronts to 3
      		IF (WF3(3)%WC.GE.WF3(2)%WC) THEN 
        			WF3(2)%WC = WF3(2)%WC 
      		ELSE
      		!!BEGIN
        			WF3(2)%WC = WF3(3)%WC
        		ENDIF
      			WF3(2)%TInfil = WF3(2)%TInfil + WF3(3)%TInfil
    				!no need to update root abs
    				n = 2
				Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(2)%y(numA+1))
				Call calcSA(n,WF3(2)%y(numA+1))	
	     			redistStatus = 30
      		IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN	!!recheck 
					!!merge 1 and 2 which was 3, move down 4 to 2
				WF3(1)%WC = WF3(2)%WC
				WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    				n = 1
				Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(1)%y(numA+1))
				Call calcSA(n,WF3(1)%y(numA+1))	
				WF3(2)%TInfil = WF3(4)%TInfil !!WF32 now WF34 since 3 was merged w/ 2 
        			n = 2
				Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil) ! move down WF34
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(2)%y(numA+1))
				Call calcSA(n,WF3(2)%y(numA+1)) 
      			redistStatus = 20
				WF3(3)%wc = 0.0
				WF3(4)%WC = 0.0
				Do j = 1, numA+1
					WF3(3)%R(j) = 0.0
					WF3(3)%x(j) = 0.0
					WF3(3)%y(j) = 0.0
					WF3(4)%R(j) = 0.0
					WF3(4)%x(j) = 0.0
					WF3(4)%y(j) = 0.0	
				Enddo
				WF3(3)%kint = 0.0
				WF3(3)%a = 0.0
				WF3(3)%b = 0.0
				WF3(3)%vol = 0.0
				WF3(3)%SA = 0.0
				WF3(3)%TInfil = 0.0 
      			WF3(4)%TInfil = 0.0
				WF3(4)%kint = 0.0
				WF3(4)%a = 0.0
				WF3(4)%b = 0.0
				WF3(4)%vol = 0.0
				WF3(4)%SA = 0.0
				IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN	!!recheck 
					WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
    					WF3(1)%WC = WF3(2)%WC
					n = 1
					Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
					Call calcXY(n)
					Call calcEllipse(n)
					Call calcVol(n,WF3(1)%y(numA+1))
					Call calcSA(n,WF3(1)%y(numA+1))	
    	    					WF3(2)%wc = 0.0
					Do j = 1, numA+1
						WF3(2)%R(j) = 0.0
						WF3(2)%x(j) = 0.0
						WF3(2)%y(j) = 0.0	
					Enddo
					WF3(2)%kint = 0.0
					WF3(2)%a = 0.0
					WF3(2)%b = 0.0
					WF3(2)%vol = 0.0
					WF3(2)%SA = 0.0
					WF3(2)%TInfil = 0.0 
      				redistStatus = 10
				ENDIF			      
			!end recheck		
			ELSE !still 3 fronts, move down 4 front which is now 3
				!!//WF3 1 stays the same; WF34 becomes WF33
        			WF3(3)%TInfil = WF3(4)%TInfil
				n = 3
				Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(3)%y(numA+1))
				Call calcSA(n,WF3(3)%y(numA+1))	
      			WF3(4)%wc = 0.0
				Do j = 1, numA+1
					WF3(4)%R(j) = 0.0
					WF3(4)%x(j) = 0.0
					WF3(4)%y(j) = 0.0	
				Enddo
				WF3(4)%kint = 0.0
				WF3(4)%a = 0.0
				WF3(4)%b = 0.0
				WF3(4)%vol = 0.0
				WF3(4)%SA = 0.0
				WF3(4)%TInfil = 0.0 
      			redistStatus = 30      		
      			IF (WF3(3)%R(numA+1).GE.WF3(2)%R(numA+1)) THEN 
      			!!BEGIN  //3 fronts to 2
      				WF3(2)%WC = WF3(3)%WC
					WF3(2)%TInfil = WF3(2)%TInfil + WF3(3)%TInfil
      				n = 2
					Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
					Call calcXY(n)
					Call calcEllipse(n)
					Call calcVol(n,WF3(2)%y(numA+1))
					Call calcSA(n,WF3(2)%y(numA+1))	
      				WF3(3)%wc = 0.0
					Do j = 1, numA+1
						WF3(3)%R(j) = 0.0
						WF3(3)%x(j) = 0.0
						WF3(3)%y(j) = 0.0	
					Enddo
					WF3(3)%kint = 0.0
					WF3(3)%a = 0.0
					WF3(3)%b = 0.0
					WF3(3)%vol = 0.0
					WF3(3)%SA = 0.0
					WF3(3)%TInfil = 0.0  
					redistStatus = 20             		
        				IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN	!!recheck 
						WF3(1)%WC = WF3(2)%WC
						WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
						n = 1
						Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
						Call calcXY(n)
						Call calcEllipse(n)
						Call calcVol(n,WF3(1)%y(numA+1))
						Call calcSA(n,WF3(1)%y(numA+1))	
    						WF3(2)%wc = 0.0
						Do j = 1, numA+1
							WF3(2)%R(j) = 0.0
							WF3(2)%x(j) = 0.0
							WF3(2)%y(j) = 0.0	
						Enddo
						WF3(2)%kint = 0.0
						WF3(2)%a = 0.0
						WF3(2)%b = 0.0
						WF3(2)%vol = 0.0
						WF3(2)%SA = 0.0
						WF3(2)%TInfil = 0.0 
      					redistStatus = 10
					ENDIF		
				
				ENDIF
      		ENDIF
    		!!END// wc3 >wc2 or z3 > z2 ;4 to 3
    		ELSE !! //still 4 fronts
    		!!BEGIN
      		!!//calculate WF3 4 depth
      		n = 4
			Call calcR(n,WF3(4)%WC - WF3(3)%WC,WF3(4)%TInfil)
			Call calcXY(n)
			Call calcEllipse(n)
			Call calcVol(n,WF3(4)%y(numA+1))
			Call calcSA(n,WF3(4)%y(numA+1))
      		IF  (WF3(4)%R(numA+1).GE.WF3(3)%R(numA+1)) THEN 
      		!!BEGIN  // 4 fronts to 3 ; 
      			!!//WF3 1 stays the same; WF34 becomes WF33
        			WF3(3)%WC = WF3(4)%WC
				WF3(3)%TInfil = WF3(3)%TInfil + WF3(4)%TInfil
				n = 3
				Call calcR(n,WF3(3)%WC - WF3(2)%WC,WF3(3)%TInfil)
				Call calcXY(n)
				Call calcEllipse(n)
				Call calcVol(n,WF3(3)%y(numA+1))
				Call calcSA(n,WF3(3)%y(numA+1))	
      			WF3(4)%wc = 0.0
				Do j = 1, numA+1
					WF3(4)%R(j) = 0.0
					WF3(4)%x(j) = 0.0
					WF3(4)%y(j) = 0.0	
				Enddo
				WF3(4)%kint = 0.0
				WF3(4)%a = 0.0
				WF3(4)%b = 0.0
				WF3(4)%vol = 0.0
				WF3(4)%SA = 0.0
				WF3(4)%TInfil = 0.0 
      			redistStatus = 30      				
				IF (WF3(3)%R(numA+1).GE.WF3(2)%R(numA+1)) THEN !check if 3 > 2
					WF3(2)%WC = WF3(3)%WC
					WF3(2)%TInfil = WF3(2)%TInfil + WF3(3)%TInfil
      				n = 2
					Call calcR(n,WF3(2)%WC - WF3(1)%WC,WF3(2)%TInfil)
					Call calcXY(n)
					Call calcEllipse(n)
					Call calcVol(n,WF3(2)%y(numA+1))
					Call calcSA(n,WF3(2)%y(numA+1))	
      				WF3(3)%wc = 0.0
					Do j = 1, numA+1
						WF3(3)%R(j) = 0.0
						WF3(3)%x(j) = 0.0
						WF3(3)%y(j) = 0.0	
					Enddo
					WF3(3)%kint = 0.0
					WF3(3)%a = 0.0
					WF3(3)%b = 0.0
					WF3(3)%vol = 0.0
					WF3(3)%SA = 0.0
					WF3(3)%TInfil = 0.0  
            			redistStatus = 20 
					IF (WF3(2)%R(numA+1).GE.WF3(1)%R(numA+1)) THEN	!!recheck 
						WF3(1)%WC = WF3(2)%WC
						WF3(1)%TInfil = WF3(1)%TInfil + WF3(2)%TInfil
						n = 1
						Call calcR(n,WF3(1)%WC - sd%wci,WF3(1)%TInfil)
						Call calcXY(n)
						Call calcEllipse(n)
						Call calcVol(n,WF3(1)%y(numA+1))
						Call calcSA(n,WF3(1)%y(numA+1))	
    						WF3(2)%wc = 0.0
						Do j = 1, numA+1
							WF3(2)%R(j) = 0.0
							WF3(2)%x(j) = 0.0
							WF3(2)%y(j) = 0.0	
						Enddo
						WF3(2)%kint = 0.0
						WF3(2)%a = 0.0
						WF3(2)%b = 0.0
						WF3(2)%vol = 0.0
						WF3(2)%SA = 0.0
						WF3(2)%TInfil = 0.0 
      					redistStatus = 10
					ENDIF		
      			ENDIF
			ENDIF !!//z4 > z3 ; 
      	
    		ENDIF !! //else, still 4 fronts
  	ENDIF !!// else, still 4 fronts for now


	END


	Subroutine CalcAvgWC(redistStatus)
c----- calcs Average WC of each cell in grid -------------------------------------------------------------
	USE ModuleDefs3
	Implicit None

	Integer i,j,k,l, redistStatus, sst
	Double Precision ErrorTest,TotalArea,CellPointCalc, areaTemp(4)
	
!	sst =2
c -----------store previous time step avg water content -----------------------
	Do i = 1,ncols
		Do j = 1,nrows
			cell(j,i)%avgWCold = cell(j,i)%avgWC
		EndDo
	EndDo
	
	Do i = 1,ncols
		Do j = 1,nrows
			Do k = 1,4
				cell(j,i)%d(k) = 0.d0
				cell(j,i)%e(k) = 0.d0	
				cell(j,i)%f(k) = 0.d0
				cell(j,i)%g(k) = 0.d0
				cell(j,i)%area(k) = 0.d0
				cell(j,i)%areaTemp(k) = 0.d0
			Enddo
		Enddo
	Enddo		
	
	If (redistStatus.EQ.10.OR.redistStatus.EQ.11) Then
		l = 1
	ElseIf (redistStatus.EQ.20.OR.redistStatus.EQ.21) Then
		l = 2
	ElseIf (redistStatus.EQ.30.OR.redistStatus.EQ.31) Then
		l = 3
	ELSE
		l = 4
	ENDIF

c-------------------------------- calculate WF3 intersections ---------------------------------------
	
	Do i = 1,ncols
		Do j= 1,nrows
			Do k = 1,l
			!!!!!!!!!!!!crosses top? !!!!!!!!!!!!!!!!!!
				ErrorTest = ((cell(j,i)%y1-WF3(k)%kint)/WF3(k)%b)**2.0 
				IF (ErrorTest.GT.1.d0) THEN
					Cell(j,i)%d(k) = 0.d0
				ELSE
					Cell(j,i)%d(k) =WF3(k)%a*(1-ErrorTest)**0.5
				ENDIF
				IF(Cell(j,i)%d(k).LT.cell(j,i)%x1.OR.cell(j,i)%d(k).GT.
     +			cell(j,i)%x2) THEN
					cell(j,i)%d(k) = 0.0
				ENDIF
			!!!!!!!!!!!!crosses bottom?!!!!!!!!!!!!!!!!!!!
				IF ((cell(j,i)%y2-WF3(k)%kint)/WF3(k)%b.GT.1.d0) THEN
					Cell(j,i)%e(k) = 0.d0
				ELSE			
					Cell(j,i)%e(k) =WF3(k)%a*(1-((cell(j,i)%y2-
     +				WF3(k)%kint)/WF3(k)%b)**2.d0)**0.5
				ENDIF
				IF(Cell(j,i)%e(k).LT.cell(j,i)%x1.OR.cell(j,i)%e(k).GT.
     +			cell(j,i)%x2) THEN
					cell(j,i)%e(k) = 0.0
				ENDIF
			!!!!!!!!!!!!!!!crosses left?!!!!!!!!!!!!!!!!!
				ErrorTest = (cell(j,i)%x1/WF3(k)%a)**2.0
							IF (ErrorTest.GT.1.d0) THEN
					Cell(j,i)%f(k) = 0.d0
				ELSE
					Cell(j,i)%f(k) =WF3(k)%b*(1-ErrorTest)**0.5
     +				+WF3(k)%kint
				ENDIF
				IF(Cell(j,i)%y1.LT.WF3(k)%kint)THEN 
					cell(j,i)%f(k)=-cell(j,i)%f(k)+2.0*WF3(k)%kint
				EndIF
				IF(Cell(j,i)%f(k).LT.cell(j,i)%y1.OR.cell(j,i)%f(k).GT.
     +			cell(j,i)%y2) THEN
					cell(j,i)%f(k) = 0.0
				ENDIF
			!!!!!!!!!!!!!!!crosses right?!!!!!!!!!!!!!!!!!
				ErrorTest = (cell(j,i)%x2/WF3(k)%a)**2.0
							IF (ErrorTest.GT.1.d0) THEN
					Cell(j,i)%g(k) = 0.d0
				ELSE
					Cell(j,i)%g(k) =WF3(k)%b*(1-ErrorTest)**0.5
     +				+WF3(k)%kint
				ENDIF
				IF(Cell(j,i)%y1.LT.WF3(k)%kint)THEN 
					cell(j,i)%g(k)=-cell(j,i)%g(k)+2.0*WF3(k)%kint
				EndIF
				IF(Cell(j,i)%g(k).LT.cell(j,i)%y1.OR.cell(j,i)%g(k).GT.
     +			cell(j,i)%y2) THEN
					cell(j,i)%g(k) = 0.0
				ENDIF
			Enddo
		Enddo
	Enddo
c------------------------------------calculate areas occupied by each front ---------------------------------
cccc Note: these areas are not necessarily the true area occupied by each front; these are temporary calculations of areas;
cccc Actual areas are calculated below and held in cell(j,i)%area
      Do i = 1,ncols
		Do j = 1,nrows
			Do k = 1,l
				IF (cell(j,i)%d(k).GT.0.d0.AND.
     +				cell(j,i)%f(k).GT.0.0) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%d(k) - 
     +					cell(j,i)%x1)*(cell(j,i)%f(k) - cell(j,i)%y1)
				ELSEIF(cell(j,i)%d(k).GT.0.d0.AND.
     +				cell(j,i)%e(k).GT.0.0.AND.cell(j,i)%d(k).
     +				GT.cell(j,i)%e(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%d(k) - 
     +					cell(j,i)%e(k))*(cell(j,i)%y2-cell(j,i)%y1)+
     +					(cell(j,i)%e(k)-cell(j,i)%x1)*(cell(j,i)%y2-
     +					cell(j,i)%y1)
				ELSEIF(cell(j,i)%g(k).GT.0.d0.AND.
     +				cell(j,i)%e(k).GT.0.0) THEN
					cell(j,i)%AreaTemp(k) = (cell(j,i)%x2-cell(j,i)%x1)* 
     +					(cell(j,i)%y2-cell(j,i)%y1)-0.5*(cell(j,i)%y2-
     +					cell(j,i)%g(k))*(cell(j,i)%x2-cell(j,i)%e(k))
      			ELSEIF(cell(j,i)%g(k).GT.0.d0.AND.
     +				cell(j,i)%f(k).GT.0.0.AND.cell(j,i)%f(k).
     +				GT.cell(j,i)%g(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%f(k) - 
     +					cell(j,i)%g(k))*(cell(j,i)%x2-cell(j,i)%x1)+
     +					(cell(j,i)%x2-cell(j,i)%x1)*(cell(j,i)%g(k)-
     +					cell(j,i)%y1)
				ELSEIF(cell(j,i)%d(k).GT.0.d0.AND.
     +				cell(j,i)%g(k).GT.0.0) THEN
					cell(j,i)%AreaTemp(k) = (cell(j,i)%x2-cell(j,i)%x1)* 
     +					(cell(j,i)%y2-cell(j,i)%y1)-0.5*(cell(j,i)%x2-
     +					cell(j,i)%d(k))*(cell(j,i)%g(k)-cell(j,i)%y1)			
				ELSEIF(cell(j,i)%f(k).GT.0.d0.AND.
     +				cell(j,i)%g(k).GT.0.0.AND.cell(j,i)%g(k).
     +				GT.cell(j,i)%f(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%g(k) - 
     +					cell(j,i)%f(k))*(cell(j,i)%x2-cell(j,i)%x1)+
     +					(cell(j,i)%x2-cell(j,i)%x1)*(cell(j,i)%y2-
     +					cell(j,i)%g(k))			
				ELSEIF(cell(j,i)%f(k).GT.0.d0.AND.
     +				cell(j,i)%e(k).GT.0.0) THEN
					cell(j,i)%AreaTemp = 0.5*(cell(j,i)%y2-
     +					cell(j,i)%f(k))*(cell(j,i)%e(k)-cell(j,i)%x1)				
				ELSEIF(cell(j,i)%e(k).GT.0.d0.AND.
     +				cell(j,i)%d(k).GT.0.0.AND.cell(j,i)%e(k).
     +				GT.cell(j,i)%d(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%e(k) - 
     +					cell(j,i)%d(k))*(cell(j,i)%y2-cell(j,i)%y1)+
     +					(cell(j,i)%y2-cell(j,i)%y1)*(cell(j,i)%d(k)-
     +					cell(j,i)%x1)
     				ELSE
					cell(j,i)%AreaTemp(k) = 0.d0																					
				ENDIF

			EndDo
		EndDo
	EndDo
  

c--------------------------------------- calculate average wc --------------------------------------------
	Do i = 1,ncols
		Do j = 1,nrows
			TotalArea = (cell(j,i)%x2-cell(j,i)%x1)*(cell(j,i)%y2-
     +			cell(j,i)%y1)
			If(cell(j,i)%AreaTemp(1).GT.0.0.AND.cell(j,i)%AreaTemp(2).
     +            LE.0.0.AND.cell(j,i)%AreaTemp(3).LE.0.0.AND.
     +            cell(j,i)%AreaTemp(4).LE.0.0) THEN
				cell(j,i)%area(1) = cell(j,i)%areaTemp(1)
				cell(j,i)%AvgWC = (cell(j,i)%Area(1)*WF3(1)%WC + 
     +			  (TotalArea - cell(j,i)%Area(1))*sd%wci)/TotalArea
			ElseIf(cell(j,i)%AreaTemp(1).GT.0.0.AND.cell(j,i)%AreaTemp(2).GT.
     +            0.0.AND.cell(j,i)%AreaTemp(3).LE.0.0.AND.
     +            cell(j,i)%AreaTemp(4).LE.0.0) THEN
				cell(j,i)%area(1) = cell(j,i)%areaTemp(1)-cell(j,i)%areaTemp(2)
				cell(j,i)%area(2) = cell(j,i)%areaTemp(2)
                cell(j,i)%AvgWC = (cell(j,i)%Area(2)*WF3(2)%WC + 
     +			  cell(j,i)%Area(1)*WF3(1)%WC + (TotalArea - 
     +			  cell(j,i)%AreaTemp(1))*sd%wci)/TotalArea   
			ElseIf(cell(j,i)%AreaTemp(1).LE.0.0.AND.cell(j,i)%AreaTemp(2).GT.
     +            0.0.AND.cell(j,i)%AreaTemp(3).LE.0.0.AND.
     +            cell(j,i)%AreaTemp(4).LE.0.0) THEN
				cell(j,i)%area(2) = cell(j,i)%areaTemp(2)
				cell(j,i)%Area(1) = TotalArea - cell(j,i)%Area(2)
				cell(j,i)%AvgWC = (cell(j,i)%Area(2)*WF3(2)%WC + 
     +			  cell(j,i)%Area(1)*WF3(1)%WC)/TotalArea
			ElseIf(cell(j,i)%AreaTemp(1).GT.0.0.AND.cell(j,i)%AreaTemp(2).GT.
     +            0.0.AND.cell(j,i)%AreaTemp(3).GT.0.0.AND.
     +            cell(j,i)%AreaTemp(4).LE.0.0) THEN
				cell(j,i)%area(3) = cell(j,i)%areaTemp(3)
				cell(j,i)%area(2) = cell(j,i)%areaTemp(2)-cell(j,i)%area(3)
				cell(j,i)%area(1) = cell(j,i)%areaTemp(1)-cell(j,i)%area(2)-
     +              cell(j,i)%area(3)
				cell(j,i)%AvgWC = (cell(j,i)%Area(3)*WF3(3)%WC + 
     +			  cell(j,i)%Area(2)*WF3(2)%WC + cell(j,i)%Area(1)*
     +			  WF3(1)%WC + (TotalArea -cell(j,i)%AreaTemp(1))*
     +			  sd%wci)/TotalArea
			ElseIf(cell(j,i)%AreaTemp(1).LE.0.0.AND.cell(j,i)%AreaTemp(2).GT.
     +            0.0.AND.cell(j,i)%AreaTemp(3).GT.0.0.AND.
     +            cell(j,i)%AreaTemp(4).LE.0.0) THEN
				cell(j,i)%Area(3) = cell(j,i)%areaTemp(3)
				cell(j,i)%Area(2) = cell(j,i)%areaTemp(2) - cell(j,i)%area(3)
				cell(j,i)%Area(1) = TotalArea - cell(j,i)%areaTemp(2)
				cell(j,i)%AvgWC = (cell(j,i)%Area(3)*WF3(3)%WC + 
     +			  cell(j,i)%Area(2)*WF3(2)%WC + cell(j,i)%Area(1)*
     +			  WF3(1)%WC)/TotalArea		
			ElseIf(cell(j,i)%AreaTemp(1).LE.0.0.AND.cell(j,i)%AreaTemp(2).LE.
     +            0.0.AND.cell(j,i)%AreaTemp(3).GT.0.0.AND.
     +            cell(j,i)%AreaTemp(4).LE.0.0) THEN
				cell(j,i)%Area(3) = cell(j,i)%areaTemp(3)
				cell(j,i)%Area(2) = TotalArea - cell(j,i)%Area(3)
				cell(j,i)%AvgWC = (cell(j,i)%Area(3)*WF3(3)%WC + 
     +			  cell(j,i)%Area(2)*WF3(2)%WC)/TotalArea			
			ElseIf(cell(j,i)%AreaTemp(1).GT.0.0.AND.cell(j,i)%AreaTemp(2).GT.
     +            0.0.AND.cell(j,i)%AreaTemp(3).GT.0.0.AND.
     +            cell(j,i)%AreaTemp(4).GT.0.0) THEN
				cell(j,i)%Area(4) = cell(j,i)%AreaTemp(4)
				cell(j,i)%Area(3) = cell(j,i)%AreaTemp(3)-cell(j,i)%Area(4)
				cell(j,i)%Area(2) = cell(j,i)%AreaTemp(2)-cell(j,i)%AreaTemp(3)
				cell(j,i)%Area(1) = cell(j,i)%AreaTemp(1)-cell(j,i)%AreaTemp(2)
				cell(j,i)%AvgWC = (cell(j,i)%Area(4)*WF3(4)%WC + 
     +			  cell(j,i)%Area(3)*WF3(3)%WC +cell(j,i)%Area(2)*WF3(2)%WC +
     +			  cell(j,i)%Area(1)*WF3(1)%WC +(TotalArea - 
     +			  cell(j,i)%AreaTemp(1))*sd%wci)/TotalArea
			ElseIf(cell(j,i)%AreaTemp(1).LE.0.0.AND.cell(j,i)%AreaTemp(2).GT.
     +            0.0.AND.cell(j,i)%AreaTemp(3).GT.0.0.AND.
     +            cell(j,i)%AreaTemp(4).GT.0.0) THEN
				cell(j,i)%Area(4) = cell(j,i)%AreaTemp(4)
				cell(j,i)%Area(3) = cell(j,i)%AreaTemp(3)-cell(j,i)%Area(4)
				cell(j,i)%Area(2) = cell(j,i)%AreaTemp(2)-cell(j,i)%AreaTemp(3)
				cell(j,i)%Area(1) = TotalArea-cell(j,i)%AreaTemp(2)
				cell(j,i)%AvgWC = (cell(j,i)%Area(4)*WF3(4)%WC + 
     +			  cell(j,i)%Area(3)*WF3(3)%WC +cell(j,i)%Area(2)*WF3(2)%WC +
     +			  cell(j,i)%Area(1)*WF3(1)%WC)/TotalArea		
			ElseIf(cell(j,i)%AreaTemp(1).LE.0.0.AND.cell(j,i)%AreaTemp(2).LE.
     +            0.0.AND.cell(j,i)%AreaTemp(3).GT.0.0.AND.
     +            cell(j,i)%AreaTemp(4).GT.0.0) THEN
				cell(j,i)%Area(4) = cell(j,i)%AreaTemp(4)
				cell(j,i)%Area(3) = cell(j,i)%AreaTemp(3)-cell(j,i)%Area(4)
				cell(j,i)%Area(2) = TotalArea-cell(j,i)%AreaTemp(3)
				cell(j,i)%AvgWC = (cell(j,i)%Area(4)*WF3(4)%WC + 
     +			  cell(j,i)%Area(3)*WF3(3)%WC +cell(j,i)%Area(2)*WF3(2)%WC)/ +
     +			  TotalArea		
			ElseIf(cell(j,i)%AreaTemp(1).LE.0.0.AND.cell(j,i)%AreaTemp(2).LE.
     +            0.0.AND.cell(j,i)%AreaTemp(3).LE.0.0.AND.
     +            cell(j,i)%AreaTemp(4).GT.0.0) THEN
				cell(j,i)%Area(4) = cell(j,i)%AreaTemp(4)
				cell(j,i)%Area(3) = TotalArea - cell(j,i)%Area(4)
				cell(j,i)%AvgWC = (cell(j,i)%Area(4)*WF3(4)%WC + 
     +			  cell(j,i)%Area(3)*WF3(3)%WC)/TotalArea
			Else
				Do k = l,1,-1
					CellPointCalc = (cell(j,i)%x1/WF3(k)%a)**2.0+
     +					  ((cell(j,i)%y1-WF3(k)%kint)/WF3(k)%b)**2.0
					If (CellPointCalc.LT.1.d0) THEN
						cell(j,i)%Area(k) = TotalArea
						cell(j,i)%AvgWC = WF3(k)%WC
						GOTO 101
					ENDIF
				EndDo
				cell(j,i)%AvgWC = sd%wci
			EndIf

101	    EndDo
	EndDo
	
	End	
	
      Subroutine CalcNormalAngle
	USE ModuleDefs3
	Implicit None
	Integer i, j, l, k
	Integer wfs, wfe
	Double Precision NA(4)
	Double Precision Ax(4),Ay(4)
	Double Precision theta(4), t(4), AxSum,AySum

	DO I = 1,ncols
		DO J = 1,nrows
			IF (cell(j,i)%Area(1).GT.0.d0.AND.cell(j,i)%Area(2).LE.0.0.AND.
     +			    cell(j,i)%Area(3).LE.0.0.AND.cell(j,i)%Area(4).LE.0.0) THEN
			   ! influenced by WF1 only
			   wfs = 1
			   wfe = 1
			ELSEIF (cell(j,i)%Area(1).LE.0.0.AND.cell(j,i)%Area(2).GT.0.d0.AND.
     +		    cell(j,i)%Area(3).LE.0.0.AND.cell(j,i)%Area(4).LE.0.0) THEN
			   ! influenced by WF2 only
			   wfs = 2
			   wfe = 2
			ELSEIF (cell(j,i)%Area(1).LE.0.0.AND.cell(j,i)%Area(2).LE.0.0.AND.
     +		    cell(j,i)%Area(3).GT.0.d0.AND.cell(j,i)%Area(4).LE.0.0) THEN
			   ! influenced by WF3 only
			   wfs = 3
			   wfe = 3
			ELSEIF (cell(j,i)%Area(1).LE.0.0.AND.cell(j,i)%Area(2).LE.0.0.AND.
     +		    cell(j,i)%Area(3).LE.0.0.AND.cell(j,i)%Area(4).GT.0.d0) THEN
			   ! influenced by WF4 only
			   wfs = 4
			   wfe = 4
		    ELSEIF (cell(j,i)%Area(1).GT.0.D0.AND.cell(j,i)%Area(2).
     +		        GT.0.D0.AND.cell(j,i)%Area(3).LE.0.0.AND.
     +              cell(j,i)%Area(4).LE.0.0) THEN
			   ! influenced by WF1 and WF2
			   wfs = 1
			   wfe = 2 
		    ELSEIF (cell(j,i)%Area(1).GT.0.D0.AND.cell(j,i)%Area(2).
     +		        GT.0.D0.AND.cell(j,i)%Area(3).GT.0.D0.AND.
     +              cell(j,i)%Area(4).LE.0.0) THEN
			   ! influenced by WF1 and WF2 and WF3
			   wfs = 1
			   wfe = 3
			ELSEIF (cell(j,i)%Area(1).GT.0.D0.AND.cell(j,i)%Area(2).
     +		        GT.0.D0.AND.cell(j,i)%Area(3).GT.0.D0.AND.
     +              cell(j,i)%Area(4).GT.0.D0) THEN
			   ! influenced by WF1 and WF2 and WF3 and WF4
			   wfs = 1
			   wfe = 4
			ELSEIF (cell(j,i)%Area(1).LE.0.0.AND.cell(j,i)%Area(2).
     +		        GT.0.D0.AND.cell(j,i)%Area(3).GT.0.D0.AND.
     +              cell(j,i)%Area(4).LE.0.0) THEN
			   ! influenced by WF2 and WF3 
			   wfs = 2
			   wfe = 3   
		    ELSEIF (cell(j,i)%Area(1).LE.0.0.AND.cell(j,i)%Area(2).
     +		        GT.0.D0.AND.cell(j,i)%Area(3).GT.0.D0.AND.
     +              cell(j,i)%Area(4).GT.0.D0) THEN
			   ! influenced by WF2 and WF3 and WF4
			   wfs = 2
			   wfe = 4
			ELSEIF (cell(j,i)%Area(1).LE.0.0.AND.cell(j,i)%Area(2).
     +		        LE.0.0.AND.cell(j,i)%Area(3).GT.0.D0.AND.
     +              cell(j,i)%Area(4).GT.0.D0) THEN
			   ! influenced by WF3 and WF4
			   wfs = 3
			   wfe = 4 
			ELSE
			   !outside all wetting fronts 
			   wfs = 5
			   wfe = 5
		    ENDIF       
				
			AxSum = 0.0
            AySum = 0.0
			Do k = wfs,wfe
			    IF (wfs.EQ.5) THEN
			        cell(j,i)%normalAngle = -9999.0
			        GOTO 999
			    EndIf     
			    theta(k) = ATan((cell(j,i)%CenterY - 
     +								   wf3(k)%kint)/cell(j,i)%CenterX) 
                IF (theta(k).LT.0.0) THEN ! adjust for correct branch of ATan function
                    theta(k) = theta(k) + 2*pi
                ENDIF
                t(k) = ATan(wf3(k)%a*Tan(theta(k))/wf3(k)%b)
                IF (t(k).LT.0.0) THEN ! adjust for correct branch of ATan function
                    t(k) = t(k) + 2*pi
                ENDIF
                NA(k) = ATan(wf3(k)%a*Tan(t(k))/wf3(k)%b)
                IF (NA(k).LT.0.0) THEN ! adjust for correct branch of ATan function
                    NA(k) = NA(k) + 2*pi
                ENDIF
                Ax(k) = cell(j,i)%Area(k)*cos(NA(k))   
			    Ay(k) = cell(j,i)%Area(k)*sin(NA(k))
			    AxSum = AxSum + Ax(k)
			    AySum = AySum + Ay(k)
            EndDo
            
            cell(j,i)%normalAngle = ATan(AySum/AxSum)
            IF (cell(j,i)%normalAngle.LT.0.0) THEN
                cell(j,i)%normalAngle = cell(j,i)%normalAngle + 2*pi
            EndIF
            
999         cell(j,i)%normalAngleDeg = cell(j,i)%normalAngle*180.0/pi
	
		ENDDO
	ENDDO

	

      END	
	

	Subroutine CalcFlowRates(deltaT)
c---------------------calculates flow rates from cells-----------------------------
	USE ModuleDefs3
	Implicit None

	Integer i,j
	Double Precision ya,yb,xa,xb,m
	Double Precision deltaT,deltaX !, pi
!	Parameter(pi = 3.14159265359d0)

	!flow rate calc q = deltaAvgWC*deltaX/deltaT
	Do i = 1,ncols
		Do j = 1,nrows
			
			m = Tan(cell(j,i)%normalAngle)
			If(cell(j,i)%normalAngle.GE.(7.0*pi/4.0).AND.
     +			cell(j,i)%normalAngle.LT.(pi/4.0))THEN
				
				ya = m*(cell(j,i)%x1-cell(j,i)%centerX) +
     +						cell(j,i)%centerY
				yb = m*(cell(j,i)%x2-cell(j,i)%centerX) +
     +						cell(j,i)%centerY
				deltaX = ((cell(j,i)%x2-cell(j,i)%x1)**2.0 +
     +					   (yb - ya)**2.0)**0.5																		
			
			ELSE
				xa = (cell(j,i)%y1-cell(j,i)%centerY-
     +					m*cell(j,i)%centerX)/m
				xb = (cell(j,i)%y2-cell(j,i)%centerY-
     +					m*cell(j,i)%centerX)/m
				deltaX = ((cell(j,i)%y2-cell(j,i)%y1)**2.0 +
     +					   (xb - xa)**2.0)**0.5

			ENDIF
				
			
			If (cell(j,i)%AvgWCold.LT.SD%wcs.AND.
     +					cell(j,i)%AvgWC.EQ.SD%wcs) THEN
				cell(j,i)%Rate%NormalRate = Q/2.d0
			Else IF (cell(j,i)%AvgWCold.EQ.SD%wcs.AND.
     +					cell(j,i)%AvgWC.EQ.SD%wcs) THEN
				cell(j,i)%Rate%NormalRate = Q/2.d0
			Else
				cell(j,i)%Rate%NormalRate = (cell(j,i)%AvgWCOld - 
     +					cell(j,i)%AvgWC)*deltaX/deltaT	
			EndIF
			
			cell(j,i)%Rate%XDirRate = cell(j,i)%Rate%NormalRate *
     +									cos(cell(j,i)%NormalAngle)
			cell(j,i)%Rate%YDirRate = cell(j,i)%Rate%NormalRate *
     +									sin(cell(j,i)%NormalAngle)
     										
		
		EndDo
	EndDo

	End
	
	Subroutine RAGridSetup
	Use ModuleDefs3
	Implicit None
	
	  Hour = 0
	  cellRA(1,1)%Percent = 0.1125
        cellRA(1,2)%Percent = 0.1125
        cellRA(1,3)%Percent = 0.027273
        cellRA(1,4)%Percent = 0.027273
        cellRA(1,5)%Percent = 0.011111
        cellRA(1,6)%Percent = 0.011111
        cellRA(1,7)%Percent = 0.003125
        cellRA(2,1)%Percent = 0.1125
        cellRA(2,2)%Percent = 0.1125
        cellRA(2,3)%Percent = 0.027273
        cellRA(2,4)%Percent = 0.027273
        cellRA(2,5)%Percent = 0.011111
        cellRA(2,6)%Percent = 0.011111
        cellRA(2,7)%Percent = 0.003125
        cellRA(3,1)%Percent = 0.027273
        cellRA(3,2)%Percent = 0.027273
        cellRA(3,3)%Percent = 0.027273
        cellRA(3,4)%Percent = 0.027273
        cellRA(3,5)%Percent = 0.011111
        cellRA(3,6)%Percent = 0.011111
        cellRA(3,7)%Percent = 0.003125
        cellRA(4,1)%Percent = 0.027273
        cellRA(4,2)%Percent = 0.027273
        cellRA(4,3)%Percent = 0.027273
        cellRA(4,4)%Percent = 0.011111
        cellRA(4,5)%Percent = 0.011111
        cellRA(4,6)%Percent = 0.011111
        cellRA(4,7)%Percent = 0.003125
        cellRA(5,1)%Percent = 0.011111
        cellRA(5,2)%Percent = 0.011111
        cellRA(5,3)%Percent = 0.011111
        cellRA(5,4)%Percent = 0.011111
        cellRA(5,5)%Percent = 0.011111
        cellRA(5,6)%Percent = 0.003125
        cellRA(5,7)%Percent = 0.003125
        cellRA(6,1)%Percent = 0.011111
        cellRA(6,2)%Percent = 0.011111
        cellRA(6,3)%Percent = 0.011111
        cellRA(6,4)%Percent = 0.011111
        cellRA(6,5)%Percent = 0.003125
        cellRA(6,6)%Percent = 0.003125
        cellRA(6,7)%Percent = 0.003125
        cellRA(7,1)%Percent = 0.003125
        cellRA(7,2)%Percent = 0.003125
        cellRA(7,3)%Percent = 0.003125
        cellRA(7,4)%Percent = 0.003125
        cellRA(7,5)%Percent = 0.003125
        cellRA(7,6)%Percent = 0.003125
        cellRA(7,7)%Percent = 0.003125

	End
	
	Subroutine calcCellRA
	Use ModuleDefs3
	Implicit None
	
	Double Precision HourlyRA,TotalRAPerDay
	Integer i, j
	
	TotalRAPerDay = 700.0
		
	
	
	If (Hour.EQ.24) Then
	  Hour = 0
	Endif
	
	IF ((Hour.GE.0.AND.Hour.LT.7).OR.(Hour.GE.18.AND.Hour.LT.24))THEN
	  HourlyRA = 0.0 !1.0/36.0
	ELSE IF (Hour.EQ.7.OR.Hour.EQ.17) THEN
	  HourlyRA = 1.0/36.0
	ELSE IF (Hour.EQ.8.OR.Hour.EQ.16) THEN
	  HourlyRA = 2.0/36.0
	ELSE IF (Hour.EQ.9.OR.Hour.EQ.15) THEN
	  HourlyRA = 3.0/36.0  
	ELSE IF (Hour.EQ.10.OR.Hour.EQ.14) THEN
	  HourlyRA = 4.0/36.0  
	ELSE IF (Hour.EQ.11.OR.Hour.EQ.13) THEN
	  HourlyRA = 5.0/36.0
	ELSE IF (Hour.EQ.12) THEN
	  HourlyRA = 6.0/36.0  
	ELSE
	  print *, 'RA Setup Error'
	ENDIF
	
!	write (300,*) 'j i RA'      
      Do i = 1, ncols
        Do j = 1, nrows
            cell(j,i)%RA = -TotalRAPerDay*HourlyRA*
     +                     cellRA(j,i)%Percent
!            write(300,*) j,i,cell(j,i)%RA
            
        ENDDO
      ENDDO
	
      Hour = Hour + 1
      
      END
      
      Subroutine CalcWFRootAbs
      Use ModuleDefs3
      Implicit None
      Double Precision TotalArea
      Integer j,i
      
      Do i = 1, 4
        WF3(i)%RootAbs = 0.d0
      EndDo
      
      Do i = 1, ncols
        Do j = 1, nrows
               TotalArea = (cell(j,i)%x2-cell(j,i)%x1)*(cell(j,i)%y2-
     +			            cell(j,i)%y1)
               WF3(1)%RootAbs = WF3(1)%RootAbs + (cell(j,i)%RA*
     +                          cell(j,i)%Area(1)/TotalArea) 
               WF3(2)%RootAbs = WF3(2)%RootAbs + (cell(j,i)%RA*
     +                          cell(j,i)%Area(2)/TotalArea)
               WF3(3)%RootAbs = WF3(3)%RootAbs + (cell(j,i)%RA*
     +                          cell(j,i)%Area(3)/TotalArea)
               WF3(4)%RootAbs = WF3(4)%RootAbs + (cell(j,i)%RA*
     +                          cell(j,i)%Area(4)/TotalArea)
                
        EndDo
      EndDo
      
      ENd 
      
      Subroutine calcWFAreas(redistStatus)
      Use ModuleDefs3
      Implicit None
      Integer i,j,l,k,redistStatus
      Double Precision ErrorTest
      
c ------------reinitialize WF3/cell intersections/areas ------------------------		
	Do i = 1,ncols
		Do j = 1,nrows
			Do k = 1,4
				cell(j,i)%d(k) = 0.d0
				cell(j,i)%e(k) = 0.d0	
				cell(j,i)%f(k) = 0.d0
				cell(j,i)%g(k) = 0.d0
				cell(j,i)%area(k) = 0.d0
				cell(j,i)%areaTemp(k) = 0.d0
			Enddo
		Enddo
	Enddo		
	
	If (redistStatus.EQ.10.OR.redistStatus.EQ.11) Then
		l = 1
	ElseIf (redistStatus.EQ.20.OR.redistStatus.EQ.21) Then
		l = 2
	ElseIf (redistStatus.EQ.30.OR.redistStatus.EQ.31) Then
		l = 3
	ELSE
		l = 4
	ENDIF

c-------------------------------- calculate WF3 intersections ---------------------------------------
	
	Do i = 1,ncols
		Do j= 1,nrows
			Do k = 1,l
			!!!!!!!!!!!!crosses top? !!!!!!!!!!!!!!!!!!
				ErrorTest = ((cell(j,i)%y1-WF3(k)%kint)/WF3(k)%b)**2.0 
				IF (ErrorTest.GT.1.d0) THEN
					Cell(j,i)%d(k) = 0.d0
				ELSE
					Cell(j,i)%d(k) =WF3(k)%a*(1-ErrorTest)**0.5
				ENDIF
				IF(Cell(j,i)%d(k).LT.cell(j,i)%x1.OR.cell(j,i)%d(k).GT.
     +			cell(j,i)%x2) THEN
					cell(j,i)%d(k) = 0.0
				ENDIF
			!!!!!!!!!!!!crosses bottom?!!!!!!!!!!!!!!!!!!!
				IF ((cell(j,i)%y2-WF3(k)%kint)/WF3(k)%b.GT.1.d0) THEN
					Cell(j,i)%e(k) = 0.d0
				ELSE			
					Cell(j,i)%e(k) =WF3(k)%a*(1-((cell(j,i)%y2-
     +				WF3(k)%kint)/WF3(k)%b)**2.d0)**0.5
				ENDIF
				IF(Cell(j,i)%e(k).LT.cell(j,i)%x1.OR.cell(j,i)%e(k).GT.
     +			cell(j,i)%x2) THEN
					cell(j,i)%e(k) = 0.0
				ENDIF
			!!!!!!!!!!!!!!!crosses left?!!!!!!!!!!!!!!!!!
				ErrorTest = (cell(j,i)%x1/WF3(k)%a)**2.0
							IF (ErrorTest.GT.1.d0) THEN
					Cell(j,i)%f(k) = 0.d0
				ELSE
					Cell(j,i)%f(k) =WF3(k)%b*(1-ErrorTest)**0.5
     +				+WF3(k)%kint
				ENDIF
				IF(Cell(j,i)%y1.LT.WF3(k)%kint)THEN 
					cell(j,i)%f(k)=-cell(j,i)%f(k)+2.0*WF3(k)%kint
				EndIF
				IF(Cell(j,i)%f(k).LT.cell(j,i)%y1.OR.cell(j,i)%f(k).GT.
     +			cell(j,i)%y2) THEN
					cell(j,i)%f(k) = 0.0
				ENDIF
			!!!!!!!!!!!!!!!crosses right?!!!!!!!!!!!!!!!!!
				ErrorTest = (cell(j,i)%x2/WF3(k)%a)**2.0
							IF (ErrorTest.GT.1.d0) THEN
					Cell(j,i)%g(k) = 0.d0
				ELSE
					Cell(j,i)%g(k) =WF3(k)%b*(1-ErrorTest)**0.5
     +				+WF3(k)%kint
				ENDIF
				IF(Cell(j,i)%y1.LT.WF3(k)%kint)THEN 
					cell(j,i)%g(k)=-cell(j,i)%g(k)+2.0*WF3(k)%kint
				EndIF
				IF(Cell(j,i)%g(k).LT.cell(j,i)%y1.OR.cell(j,i)%g(k).GT.
     +			cell(j,i)%y2) THEN
					cell(j,i)%g(k) = 0.0
				ENDIF
			Enddo
		Enddo
	Enddo
c------------------------------------calculate areas occupied by each front ---------------------------------
cccc Note: these areas are not necessarily the true area occupied by each front; these are temporary calculations of areas;
cccc Actual areas are calculated below and held in cell(j,i)%area
      Do i = 1,ncols
		Do j = 1,nrows
			Do k = 1,l
				IF (cell(j,i)%d(k).GT.0.d0.AND.
     +				cell(j,i)%f(k).GT.0.0) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%d(k) - 
     +					cell(j,i)%x1)*(cell(j,i)%f(k) - cell(j,i)%y1)
				ELSEIF(cell(j,i)%d(k).GT.0.d0.AND.
     +				cell(j,i)%e(k).GT.0.0.AND.cell(j,i)%d(k).
     +				GT.cell(j,i)%e(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%d(k) - 
     +					cell(j,i)%e(k))*(cell(j,i)%y2-cell(j,i)%y1)+
     +					(cell(j,i)%e(k)-cell(j,i)%x1)*(cell(j,i)%y2-
     +					cell(j,i)%y1)
				ELSEIF(cell(j,i)%g(k).GT.0.d0.AND.
     +				cell(j,i)%e(k).GT.0.0) THEN
					cell(j,i)%AreaTemp(k) = (cell(j,i)%x2-cell(j,i)%x1)* 
     +					(cell(j,i)%y2-cell(j,i)%y1)-0.5*(cell(j,i)%y2-
     +					cell(j,i)%g(k))*(cell(j,i)%x2-cell(j,i)%e(k))
      			ELSEIF(cell(j,i)%g(k).GT.0.d0.AND.
     +				cell(j,i)%f(k).GT.0.0.AND.cell(j,i)%f(k).
     +				GT.cell(j,i)%g(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%f(k) - 
     +					cell(j,i)%g(k))*(cell(j,i)%x2-cell(j,i)%x1)+
     +					(cell(j,i)%x2-cell(j,i)%x1)*(cell(j,i)%g(k)-
     +					cell(j,i)%y1)
				ELSEIF(cell(j,i)%d(k).GT.0.d0.AND.
     +				cell(j,i)%g(k).GT.0.0) THEN
					cell(j,i)%AreaTemp(k) = (cell(j,i)%x2-cell(j,i)%x1)* 
     +					(cell(j,i)%y2-cell(j,i)%y1)-0.5*(cell(j,i)%x2-
     +					cell(j,i)%d(k))*(cell(j,i)%g(k)-cell(j,i)%y1)			
				ELSEIF(cell(j,i)%f(k).GT.0.d0.AND.
     +				cell(j,i)%g(k).GT.0.0.AND.cell(j,i)%g(k).
     +				GT.cell(j,i)%f(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%g(k) - 
     +					cell(j,i)%f(k))*(cell(j,i)%x2-cell(j,i)%x1)+
     +					(cell(j,i)%x2-cell(j,i)%x1)*(cell(j,i)%y2-
     +					cell(j,i)%g(k))			
				ELSEIF(cell(j,i)%f(k).GT.0.d0.AND.
     +				cell(j,i)%e(k).GT.0.0) THEN
					cell(j,i)%AreaTemp = 0.5*(cell(j,i)%y2-
     +					cell(j,i)%f(k))*(cell(j,i)%e(k)-cell(j,i)%x1)				
				ELSEIF(cell(j,i)%e(k).GT.0.d0.AND.
     +				cell(j,i)%d(k).GT.0.0.AND.cell(j,i)%e(k).
     +				GT.cell(j,i)%d(k)) THEN
					cell(j,i)%AreaTemp(k) = 0.5*(cell(j,i)%e(k) - 
     +					cell(j,i)%d(k))*(cell(j,i)%y2-cell(j,i)%y1)+
     +					(cell(j,i)%y2-cell(j,i)%y1)*(cell(j,i)%d(k)-
     +					cell(j,i)%x1)
     				ELSE
					cell(j,i)%AreaTemp(k) = 0.d0																					
				ENDIF

			EndDo
		EndDo
	EndDo

      End
      
      Subroutine volZ_zbrent(R0,R1)
C	implicit solution of the ellipsoid volume equ. 5-16 to 5-19 in Leslie's Ph.D
C   driver for routine zbrent  --- R. Munoz-Carpena--
	Use ModuleDefs3
	implicit Double Precision (a-h,o-z)
	PARAMETER(N=100,NBMAX=20,X1=0.0,X2=100.0)
	EXTERNAL func
      
c --- inputs, need to be passed on to this subroutine
!	common/volZ/xk,xe,volZ,xaZ,xbZ
!	xe=.700
!	xk=2.0
!	volZ=1330.02466582
c --- end of inputs

	nb=NBMAX	
c --- automatically bracket the solution
	call zbrak(func,X1,X2,N,xbZ1,xbZ2,nb)
	tol=(1.0e-8)*(xbZ1+xbZ2)/2.0
!	Print*, volZ,xaz,xbz,xk,xe
c --- find optimized R1 and R0 for a given volume given k, e
	R1=zbrent(func,xbZ1,xbZ2,tol)
	R0=SQRT((1.0-(xk**2.d0)/(xbZ**2.d0))*xaZ**2.d0)
	
!	print *,'R0,R1=', R0,R1,xe

	return
	END


	function func(R1)
C	equ. 5-16 to 5-19 in Leslie's Ph.D
	Use ModuleDefs3
	implicit Double Precision (a-h,o-z)
!	common/volZ/xk,xe,volZ,xaZ,xbZ
	
	
	xbZ=R1-xk
	xaZ=xbZ*xe
	p1 =xk*(R1/xbZ)**2
	p2 =xk**2*R1/xbZ**2
	p3=R1**3/(3*xbZ**2)
	func=pi*xaZ**2*(R1+p1-p2-p3)-volZ
		
	return
	end



      SUBROUTINE zbrak(fx,x1,x2,n,xbZ1,xbZ2,nb)
C     (C) Copr. 1986-92 Numerical Recipes Software #>0K!.
	Use ModuleDefs3  
	  implicit Double Precision (a-h,o-z)
	  DIMENSION xbZ1(nb),xbZ2(nb)
      EXTERNAL fx
!	  common/volZ/xk,xe,volZ,xaZ,xbZ

      nbb=0
      x=x1
      dx=(x2-x1)/n
      fp=fx(x)
      do 11 i=1,n
        x=x+dx
        fc=fx(x)
        if(fc*fp.lt.0.) then
          nbb=nbb+1
          xbZ1(nbb)=x-dx
          xbZ2(nbb)=x
          if(nbb.eq.nb)goto 1
        endif
        fp=fc
11    continue
1     continue
      nb=nbb
      return
      END


      FUNCTION zbrent(func,x1,x2,tol)
C     (C) Copr. 1986-92 Numerical Recipes Software #>0K!.
	Use ModuleDefs3
	  implicit Double Precision (a-h,o-z)
      EXTERNAL func
      PARAMETER (ITMAX=500,EPS=3.e-8)
!	  common/volZ/xk,xe,volZ,xaZ,xbZ
	  
	a=x1
      b=x2
      fa=func(a)
      fb=func(b)
      if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.))pause
     *'root must be bracketed for zbrent'
      c=b
      fc=fb
      do 11 iter=1,ITMAX
        if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
          c=a
          fc=fa
          d=b-a
          e=d
        endif
        if(abs(fc).lt.abs(fb)) then
          a=b
          b=c
          c=a
          fa=fb
          fb=fc
          fc=fa
        endif
        tol1=2.*EPS*abs(b)+0.5*tol
        xm=.5*(c-b)
        if(abs(xm).le.tol1 .or. fb.eq.0.)then
          zbrent=b
          return
        endif
        if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
          s=fb/fa
          if(a.eq.c) then
            p=2.*xm*s
            q=1.-s
          else
            q=fa/fc
            r=fb/fc
            p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
            q=(q-1.)*(r-1.)*(s-1.)
          endif
          if(p.gt.0.) q=-q
          p=abs(p)
          if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
            e=d
            d=p/q
          else
            d=xm
            e=d
          endif
        else
          d=xm
          e=d
        endif
        a=b
        fa=fb
        if(abs(d) .gt. tol1) then
          b=b+d
        else
          b=b+sign(tol1,xm)
        endif
        fb=func(b)
11    continue
      pause 'zbrent exceeding maximum iterations'
      zbrent=b
      return
      END

	



