! Subroutine to determine the average infiltration rate fave during the time interval
! Using Green-Ampt infiltration model
! By P. Shen, July 2016, HKUST
!------------------------------------------
! Input:
!   ci - the cumulative infiltartion at the beginning of the time interval dt.
!   inflx - water flow in during dt.
! Output:
!   fave (ir) - the average infiltration rate during dt.
!   ci - the cumulative infiltration at the end of dt.

Subroutine infr(imx1,fave,ci,tempci,inflx,dt)
use input_vars; use grids
implicit none
integer imx1,i
double precision:: fave(imx1),ci(imx1),tempci(imx1),inflx(imx1),dt
double precision:: cinext,fnow,ftemp,cinow,tempcinext,dci,cip,dtp
double precision:: psiti,ksti,delth,tolerr

tolerr=1e-8 ! iteration error

do i=1,imx1

psiti=psit(zo(i))
ksti=kst(zo(i))
delth=thsatt(zo(i))-thinit(zo(i))
cinow=ci(i)
cinext=0.

! case 1. exfiltration at cells where water table is initially at at the ground surface
! 'depth' and 'rizero' are defined in 'grids' module
if (depth(i)==0. .and. rizero(i)<0.) then
    fave(i)=0.
    ! ci does not change
    
else 

! case 2. infiltration
! calculate infiltration rate fnow, at the beginning of dt  
if (cinow==0) then
    fnow=100. ! use a large enough number to replace infinity
else
    fnow=ksti*(psiti*delth+cinow)/cinow
end if

if (fnow<=inflx(i)) then 
! runoff occurs at the beginning of the time interval

    tempcinext=cinow ! assign an initial value to start iteration
    dci=1.
    do while (abs(dci)>=tolerr)
        cinext=cinow+psiti*delth*log((tempcinext+psiti*delth)/(cinow+psiti*delth))+ksti*dt
        dci=cinext-tempcinext
        tempcinext=cinext
    end do
    tempci(i)=cinext
    fave(i)=(cinext-cinow)/dt
    
else 
! No runoff at the beginning of the interval
    
    ! Calculate tentative values
    tempcinext=cinow+inflx(i)*dt
    ftemp=ksti*(psiti*delth+tempcinext)/tempcinext
    
    if (ftemp<=inflx(i)) then 
        ! runoff occurs during the interval
        cip=ksti*psiti*delth/(inflx(i)-ksti)
        dtp=(cip-cinow)/inflx(i)
        dci=1.
        tempcinext=ksti
        do while (abs(dci)>=tolerr)
            cinext=cip+psiti*delth*log((tempcinext+psiti*delth)/(cip+psiti*delth))+ksti*(dt-dtp)
            dci=cinext-tempcinext
            tempcinext=cinext
        end do
        tempci(i)=cinext
        fave(i)=(cinext-cinow)/dt
     else 
        ! no runoff throughout the interval
        tempci(i)=tempcinext
        fave(i)=inflx(i)
     end if
     
! if (fnow<=inflx(i)) then            
end if
! if (depth(i)==0. .and. rizero(i)<0.) then
end if 
  
end do

end
