c Water flow subroutines calculate the pressure head (hNew) as a 
c function of depth and the water content (ThNew). 
      subroutine WatFlow(NumNP,NTab,NMat,hTab,ConTab,CapTab,hNew,
     ;                   hOld,MatNum,ParD,Con,Cap,ConSat,Ah,AK,ATh,
     ;                   hSat,hTemp,KodTop,KodBot,rTop,rBot,dt,
     ;                   z,Sink,P,R,S,FreeD,hTop,hBot,WLayer,Iter,
     ;                   ItCum,TopInf,KTOld,KBOld,dtMin,t,tOld,
     ;                   ConvgF,TheTab,ThNew,ThOld,thr,ths,vTop,
     ;                   AThS,ThRR,ConO,ConR,AKS,hCritA,iModel,
     ;                   lTable,ierr)
      implicit none
      include "defs.fi"
      logical ConvgF,ItCrit,FreeD,TopInf,WLayer,lTable
      integer iModel,MatNum,KodTop,KodBot,KTOld,KBOld,ItCum,Iter,
     ;        NMat,NTab,NumNP,ierr
      double precision hTab,ConTab,CapTab,hNew,hOld,ParD,Con,Cap,
     ;                 Ah,Ak,ATh,hSat,hTemp,rTop,rBot,dt,z,Sink,hTop,
     ;                 hBot,dtMin,TheTab,ThNew,ThOld,thr,ths,AThs,vTop,
     ;                 ThRR,ConO,ConR,AKS,ConSat,t,tOld,hCritA
      double precision P,R,S,PB,RB,SB,PT,RT,ST
      dimension z(NumNP),hNew(NumNP),hOld(NumNP),hTemp(NumNP),
     ;          MatNum(NumNP),ParD(NTPar,*),Sink(NumNP),hTab(NTabD),
     ;          ConTab(NTabD,NMat),CapTab(NTabD,NMat),Con(NumNP),
     ;          Cap(NumNP),ConSat(NMat),P(NumNP),R(NumNP),S(NumNP),
     ;          hSat(NMat),Ah(NumNP),AK(NumNP),ATh(NumNP),ths(NMat),
     ;          TheTab(NTabD,NMat),ThNew(NumNP),ThOld(NumNP),
     ;          ConO(NumNP),AThS(NumNP),ThRR(NumNP),
     ;          ConR(NumNP),AKS(NumNP),thr(NMat)
c
c Local variables
      double precision rMax,rMin,EpsTh,EpsH,Th
      integer m,i
      double precision zero,one
      data zero,one /0.0D00,1.0D00/
c
      rMax=1.D+10
      rMin=1.D-100
      ierr=0
11    continue

      Iter=0
      ConvgF=.true.
12    continue

c     Generate terms of matrix equation and solve by Gauss elimination

      call SetMat(NumNP,NTab,NMat,Iter,hTab,ConTab,CapTab,hNew,MatNum,
     ;            ParD,Con,Cap,ConSat,Ah,AK,ATh,hSat,hTemp,TheTab,ThNew,
     ;            thr,ths,ConO,AThS,ThRR,ConR,AKS,iModel,lTable)
      call Reset (NumNP,rTop,rBot,dt,z,hOld,Con,Cap,WLayer,hNew,
     ;            Sink,P,R,S,PB,RB,SB,PT,RT,ST,FreeD,ThNew,ThOld,vTop)
      call Shift (KodTop,rTop,rBot,hTop,hBot,WLayer,hCritA,
     ;            Con(NumNP),Con(NumNP-1),hNew(NumNP),hNew(NumNP-1),
     ;            z(NumNP),z(NumNP-1),TopInf,KodBot,Con(1),Con(2),
     ;            hNew(1),hNew(2),z(1),z(2),ThNew(NumNP),
     ;            ThOld(NumNP),Sink(NumNP),ThNew(1),ThOld(1),Sink(1),dt,
     ;            hOld(NumNP))
      call veccpr(1,NumNP,hNew,hTemp,ierr)
      call Gauss (NumNP,KodTop,KodBot,hTop,hBot,hNew,P,R,S,PB,RB,SB,PT,
     ;            RT,ST,rMin,rMax,hCritA)
      Iter =Iter+1
      ItCum=ItCum+1

c     Test for convergence
      ItCrit=.true.
      do 14 i=1,NumNP
        m=MatNum(i)
        EpsTh=zero
        EpsH=zero
        if(hTemp(i).lt.hSat(m).and.hNew(i).lt.hSat(m)) then
          Th=ThNew(i)+Cap(i)*(hNew(i)-hTemp(i))/(ths(m)-thr(m))/ATh(i)
          EpsTh=dabs(ThNew(i)-Th)
        else
          EpsH=dabs(hNew(i)-hTemp(i))
        end if
        if(EpsTh.gt.5.0D-04.or.EpsH.gt.2.0D-01.or.
     ;     dabs(hNew(i)).gt.rMax*0.999D00) then
          ItCrit=.false.
          if(abs(hNew(i)).gt.rMax*0.999D00) Iter=20
          goto 15
        end if
14    continue
15    continue
      if(ItCrit) then
        do 26 i=1,NumNP
          ThNew(i)=ThNew(i)+Cap(i)*(hNew(i)-hTemp(i))
26      continue
      end if
      if(.not.ItCrit.or.Iter.le.one) then
        if(Iter.lt.20) then
          goto 12
        else if(dt.le.dtMin) then
          ierr=71
          return
        else
          do 16 i=1,NumNP
            hNew(i) =hOld(i)
            hTemp(i)=hOld(i)
16        continue
          KodTop=KTOld
          KodBot=KBOld
c         Reduce dt and recalculate the time(t)
          dt=max(dt/3.0D00,dtMin)
          t=tOld+dt
          goto 11
        end if
      end if
      if(Wlayer) then
        if(hNew(NumNP).gt.zero) then
          KodTop=4
          hTop=zero
        end if
      end if
      return
      end

c     Reset the boundary elements using the iterated values of pressure
c     head and recalculate the three vectors of the tridiagonal system
      subroutine Reset(N,rTop,rBot,dt,z,hOld,Con,Cap,WLayer,hNew,
     ;                 Sink,P,R,S,PB,RB,SB,PT,RT,ST,FreeD,ThNew,ThOld,
     ;                 vTop)
      implicit none
      logical WLayer,FreeD
      double precision z,hOld,P,R,S,PB,RB,SB,PT,RT,ST,Con,Cap,
     ;                 hNew,Sink,thNew,thOld,rTop,rBot,dt,vTop
      dimension z(N),hOld(N),hNew(N),P(N),R(N),S(N),Con(N),Cap(N),
     ;          Sink(N),ThNew(N),ThOld(N)
c
c     Local variables
      double precision dzA,dzB,dz,ConA,ConB,A2,A3,B,F2
      double precision zero,one,two
      data zero,one,two /0.0D00, 1.0D00, 2.0D00/
      integer i,N
c     Finite differences
      dzB=z(2)-z(1)
      dz=dzB/two
      ConB=(Con(1)+Con(2))/two
      S(1)=-ConB/dzB
      if(FreeD) then
        RB=ConB/dzB
        SB=-RB
        PB=zero
      else
        RB=ConB/dzB+Cap(1)*dz/dt
        SB=-ConB/dzB             
        PB=ConB-Sink(1)*dz+Cap(1)*dz*hNew(1)/dt-
     ;     (ThNew(1)-ThOld(1))*dz/dt+rBot
      end if
      do 11 i=2,N-1
        dzA=z(i)-z(i-1)
        dzB=z(i+1)-z(i)
        dz=(dzA+dzB)/two
        ConA=(Con(i)+Con(i-1))/two
        ConB=(Con(i)+Con(i+1))/two
        A2=ConA/dzA+ConB/dzB
        A3=-ConB/dzB
        B =(ConA-ConB)
        F2=Cap(i)*dz
        R(i)=A2+F2/dt
        P(i)=F2*hNew(i)/dt-(ThNew(i)-ThOld(i))*dz/dt-B-Sink(i)*dz
        S(i)=A3
11    continue
      dzA=z(N)-z(N-1)
      dz=dzA/two
      ConA=(Con(N)+Con(N-1))/two
      RT=ConA/dzA+Cap(N)*dz/dt
      ST=-ConA/dzA
      PT=Cap(N)*dz*hNew(N)/dt-(ThNew(N)-ThOld(N))*dz/dt-
     ;   Sink(N-1)*dz-ConA
      vTop=-dble(ST)*hNew(N-1)-dble(RT)*hNew(N)+dble(PT)
      PT=PT-rTop
      if(WLayer) then
        if(hNew(N).gt.zero) then
          RT=RT+one/dt
          PT=PT+max(hOld(N),zero)/dt
        else
          PT=PT+max(hOld(N),zero)/dt
        end if
      end if
      return
      end

************************************************************************
      subroutine Gauss(N,KodTop,KodBot,hTop,hBot,hNew,P,R,S,PB,RB,SB,PT,
     ;                 RT,ST,rMin,rMax,hCritA)
      implicit none
      integer  N,KodTop,KodBot
      double precision P,R,S,PB,RB,SB,PT,RT,ST,rMin,rMax,hTop,hBot,hNew,
     ;                 hCritA
      dimension hNew(N),P(N),R(N),S(N)
c
c Local variables
      integer i
c     Forward
      if(KodBot.ge.0) then
        P(2)=P(2)-S(1)*hBot
      else
        if(dabs(RB).lt.rMin) RB=rMin
        P(2)=P(2)-PB*S(1)/RB
        R(2)=R(2)-SB*S(1)/RB
      end if
      do i=3,N-1
        if(dabs(R(i-1)).lt.rMin) R(i-1)=rMin
        P(i)=P(i)-P(i-1)*S(i-1)/R(i-1)
        R(i)=R(i)-S(i-1)*S(i-1)/R(i-1)
      enddo
      if(KodTop.gt.0) then
        P(N-1)=P(N-1)-S(N-1)*hTop
      else
        if(dabs(R(N-1)).lt.rMin) R(N-1)=rMin
        P(N)=PT-P(N-1)*ST/R(N-1)
        R(N)=RT-S(N-1)*ST/R(N-1)
      end if

c     Back
      if(dabs(R(N-1)).lt.rMin) R(N-1)=rMin
      if(KodTop.ge.0) then
        hNew(N)=hTop
        hNew(N-1)=dble(P(N-1)/R(N-1))
      else
        hNew(N)=max(hCritA,P(N)/R(N))
        hNew(N-1)=dble((P(N-1)-S(N-1)*hNew(N))/R(N-1))
      end if
      do i=N-2,2,-1
        if(dabs(R(i)).lt.rMin) R(i)=rMin
        hNew(i)=dble((P(i)-S(i)*hNew(i+1))/R(i))
      enddo
      if(KodBot.ge.0) then
        hNew(1)=hBot
      else
        if(dabs(RB).lt.rMin) RB=rMin
        hNew(1)=dble((PB-SB*hNew(2))/RB)
      end if
      do i=1,N
        if(abs(hNew(i)).gt.rMax) hNew(i)=sign(rMax,hNew(i))
      enddo      
      return
      end

************************************************************************
      subroutine Shift(KodTop,rTop,rBot,hTop,hBot,WLayer,hCritA,
     ;                 ConTop,ConBlw,hNT,hNB,zTop,zBlw,TopInf,KodBot,
     ;                 ConBot,ConAbv,hN1,hN2,zBot,zAbv,ThNewT,
     ;                 ThOldT,SinkTop,ThNewB,ThOldB,SinkBot,dt,hOT)
      implicit none
      logical WLayer,TopInf
      integer KodTop,KodBot
      double precision rTop,rBot,hTop,hBot,hCritA,ConTop,ConBlw,hNt,
     ;                 hNb,zTop,zBlw,ConBot,ConAbv,hN1,hN2,zBot,
     ;                 zAbv,ThNewT,ThOldT,SinkTop,ThNewB,ThOldB,
     ;                 SinkBot,dt,hOT
c
      double precision vTop
      double precision zero,one,two,eps
      data zero,one,two,eps /0.0D00, 1.0D00, 2.0D00, 1.0D-03/
    
c     Atmospheric boundary condition
      if(TopInf.and.abs(KodTop).eq.4) then
        if(KodTop.gt.0) then
          vTop=-(ConTop+ConBlw)/two*((hNT-hNB)/(zTop-zBlw)+one)-
     ;       (ThNewT-ThOldT)*(zTop-zBlw)/two/dt-SinkTop*(zTop-zBlw)/two
          if(abs(vTop).gt.abs(rTop).or.vTop*rTop.le.zero .or.
     ;      (rTop.lt.zero .and.(abs(vTop)/abs(rTop).lt.eps))) then
            KodTop=-4
          endif
        else
          if(.not.WLayer) then
            if(hNT.gt.zero) then
              KodTop=4
              hTop=zero
            end if
          end if
          if(hNT.le.hCritA) then
            KodTop=4
            hTop=hCritA
          end if
        end if
      end if
      return
      end
c
c     Build interpolation tables for hydraulic conductivity, capacitance
c     and water content
      subroutine setTab(NTab,NMat,NumNP,hTab,iModel,ParD,hSat,thr,ths,
     ;                  conSat,conTab,capTab,theTab,ATh)
      implicit none
      include "defs.fi"
      integer NTab,NMat,NumNP,iModel
      double precision hTab,hSat,thr,ths,conSat,conTab,capTab,theTab,
     ;                 ParD,ATh
      dimension hTab(NTab),hSat(NMat),thr(NMat),ths(NMat),conSat(NMat),
     ;          conTab(NTabD,NMat),capTab(NTabD,NMat),
     ;          theTab(NTabD,NMat),ParD(NTPar,*),ATh(NumNP)
c
      integer i,iMat,k,j
      double precision Par,dlh,alh,hR,fH,fQ,fK,fC,zero,one
      external fH,fQ,fK,fC
      dimension Par(NTPar)
      data zero,one /0.0D00, 1.0D00/
c
      hTab(1) = -1.0D-06
      hTab(NTab) = -1.0D04
      dlh = (log10(-hTab(NTab)) - log10(-hTab(1)))/dble(NTab-1)
      do i=2,NTab
        alh = log10(-hTab(1)) + dble(i-1)*dlh
        hTab(i) = -(1.D01**alh)
      enddo
      do iMat=1,NMat
        do k=1,NTPar
          Par(k)=ParD(k,iMat)
        enddo
c       FH gets theta instead of scaled water content
        hSat(i) = FH(iModel,Par(2),ParD,iMat)
        hR = FH(iModel,Par(1),ParD,iMat)
        thr(iMat)   =FQ(iModel, Hr ,Par)
        ths(iMat)   =FQ(iModel,zero,Par)
        conSat(iMat)=FK(iModel,zero,Par)
        do j=1,NTab
          conTab(j,iMat)=FK(iModel,hTab(j),Par)
          capTab(j,iMat)=FC(iModel,hTab(j),Par)
          theTab(j,iMat)=FQ(iModel,hTab(j),Par)
        enddo
      enddo
c
c     ATh was defined in the earlier version of E-G but only took
c     on values of unity
      do i=1,NumNP
        ATh(i)=one
      enddo
      return
      end

************************************************************************
c
c     Calculate hydraulic conductivity, capacitance and theta for all nodes
      subroutine SetMat(NumNP,NTab,NMat,Iter,hTab,ConTab,CapTab,hNew,
     ;                  MatNum,ParD,Con,Cap,ConSat,Ah,AK,ATh,hSat,
     ;                  hTemp,TheTab,theta,thr,ths,ConO,AThS,ThRR,ConR,
     ;                  AKS,iModel,lTable)
      implicit none
      include "defs.fi"
      logical   lTable
      integer   NumNP,NTab,NMat,MatNum,iModel,Iter
      double precision hTab,conTab,capTab,hNew,ParD,Con,Cap,conSat,
     ;          Ah,AK,ATh,hSat,hTemp,theTab,theta,thr,ths,conO,
     ;          AThS,ThRR,conR,AKS,tempPar
      dimension hTab(NTab),ConTab(NTabD,NMat),CapTab(NTabD,NMat),
     ;          hNew(NumNP),MatNum(NumNP),ParD(NTPar,*),Con(NumNP),
     ;          Cap(NumNP),ConSat(NMat),Ah(NumNP),AK(NumNP),
     ;          ATh(NumNP),hSat(NMat),hTemp(NumNP),TheTab(NTabD,NMat),
     ;          theta(NumNP),thr(NMat),ths(NMat),ConO(NumNP),
     ;          AThS(NumNP),ThRR(NUmNP), ConR(NumNP),AKS(NumNP),
     ;          tempPar(NTPar)
c
      integer i,j,M,iT
      double precision alh1,dlh,hi1,hi2,hiM,Coni,dh,Capi,Thei,zero
      double precision FC,FQ,FK,small
      external FC,FQ,FK
c
      data zero,small /0.0D00, 1.D-37/
      alh1=log10(-hTab(1))
      dlh =(log10(-hTab(NTab))-alh1)/(NTab-1)
      do 11 i=1,NumNP
        M=MatNum(i)
        do j=1,NTPar
          tempPar(j)=ParD(j,M)
        enddo
c
        hi1=min(hSat(M),hTemp(i)/Ah(i))
        hi2=min(hSat(M), hNew(i)/Ah(i))
        hiM=0.1D00*hi1+0.9D00*hi2
        if(hi1.ge.hSat(M).and.hi2.ge.hSat(M)) then
          Coni=ConSat(M)
        else if(hiM.gt.hTab(NTab).and.hiM.le.hTab(1).and.lTable) then
          iT=int((log10(-hiM)-alh1)/dlh)+1
          dh=(hiM-hTab(iT))/(hTab(iT+1)-hTab(iT))
          Coni=ConTab(iT,M)+(ConTab(iT+1,M)-ConTab(iT,M))*dh
        else
          Coni=FK(iModel,hiM,tempPar)
        end if
        if(hiM.ge.hSat(M)) then
          Capi=zero
          Thei=ths(M)
        else if(hiM.ge.hTab(NTab).and.hiM.le.hTab(1).and.lTable) then
          iT=int((log10(-hiM)-alh1)/dlh)+1
          dh=(hiM-hTab(iT))/(hTab(iT+1)-hTab(iT))
          Capi=CapTab(iT,M) + (CapTab(iT+1,M)-CapTab(iT,M))*dh
          Thei=TheTab(iT,M) + (TheTab(iT+1,M)-TheTab(iT,M))*dh
        else
          Capi=FC(iModel,hiM,tempPar)
          Thei=FQ(iModel,hiM,tempPar)
        end if
        Con(i)=max(Coni*AK(i)*AKS(i),small)
        Cap(i)=Capi*ATh(i)*AThS(i)/Ah(i)
        theta(i)=thr(M)+(Thei-thr(M))*ATh(i)*AThS(i)
        if(Iter.eq.0) then
          ConO(i)=Con(i)
        end if
11    continue
      return
      end

      subroutine W_Flux(NumNP,dt,vTop,vBot,vNew,ThNew,ThOld,z,Sink)
      implicit none
      integer  NumNP
      double precision z,vNew,ThNew,ThOld,Sink,vTop,vBot,dt
      dimension z(NumNP),vNew(NumNP),ThNew(NumNP),ThOld(NumNP),
     ;          Sink(NumNP)
c
c     Local variables
      integer i
      double precision dV,SS,dz,zero,two
      data zero,two / 0.0D+00,2.0D+00 /
      dV=zero
      SS=zero
      vNew(NumNP)=vTop
      do 10 i=NumNP-1,1,-1
        dz=z(i+1)-z(i)
        dV=dV+dz*(ThNew(i)+ThNew(i+1)-ThOld(i)-ThOld(i+1))/two
        SS=SS+dz*(Sink(i)+Sink(i+1))/two
        vNew(i)=vTop+dV/dt+SS
10    continue
      vBot=vNew(1)
      return
      end

      subroutine initBal(accVTop,accVBot,accSink,totInfil)
      implicit none
      double precision accVTop,accVBot,accSink,totInfil
      double precision zero
      data zero /0.0D00/
      accVTop=zero
      accVBot=zero
      accSink=zero
      totInfil=zero
      return
      end
c
c     Calculation of cumulative variables related to water flow
      subroutine W_Cumul(N,NLayers,cropNum,dt,dz,vTop,vBot,vNew,vRoot,
     ;                   accVTop,totInfil,accVBot,accSink,alfaEff,
     ;                   alphaG,sumEff,sumAlfG,TPotMax,accETp,sNoComp,
     ;                   alfaGHS,sLayTopQ,sLayBotQ,sLaySink,nodeZ,
     ;                   Beta,cropProp,crpIndex,Sink,layerDep,counter,
     ;                   lPlant,ierr)
      implicit none
      include "defs.fi"
      integer N,NLayers,cropNum,crpIndex,counter,ierr
      logical lPlant
      double precision dz,dt,vTop,vNew,vRoot,accVTop,accVBot,accSink,
     ;                 alfaEff,alphaG,sumEff,sumAlfG,TPotMax,
     ;                 accETp,sNoComp,alfaGHS,vBot,sLayTopQ,sLayBotQ,
     ;                 sLaySink,nodeZ,Sink,layerDep,Beta,
     ;                 cropProp,totInfil
      dimension        vNew(N),Sink(N),slayTopQ(NLayers),nodeZ(N),
     ;                 sLayBotQ(NLayers),sLaySink(NLayers),Beta(N),
     ;                 layerDep(NLayers),
     ;                 cropProp(nCropPar,*)
c
      integer i,nodeTop,nodeNum,nodOfZ
      double precision zero,two,sinkSum,sumSinkL,currDep
      external sumSinkL,nodOfZ
      data zero,two /0.0D00, 2.0D00/
c
      counter=counter+1
      accVTop=accVTop + vTop*dt
      accVBot=accVBot + vBot*dt
      if (vTop.lt.zero) then
        totInfil=totInfil + vTop*dt
      endif
      if (lPlant) then
        accETp=accETp + TPotMax*dt
c
c       Add the total root water uptake for the time step
        accSink=accSink + vRoot*dt
c       Sum the reduction coefficients
        if (cropNum.ne.0) then
          sumEff=sumEff + alfaEff
          sumAlfG=sumAlfG + alphaG
          sNoComp=sNoComp + alfaGHS
        endif
      endif
c
c     Layer-specific accumulation
      nodeTop=N
      do i=NLayers,1,-1
        currDep=layerDep(i)
        nodeNum=nodOfZ(N,nodeZ,currDep,ierr)
        if (ierr.gt.0) return
        if (i.eq.NLayers) then
c         water flow into top layer at soil surface
          sLayTopQ(i)=sLayTopQ(i)+vTop*dt
          sLayBotQ(i)=sLayBotQ(i)+vNew(nodeNum)*dt
        elseif (i.eq.1) then
          sLayTopQ(1)=sLayBotQ(2)
c         water flow from bottom layer
          sLayBotQ(1)=sLayBotQ(1)+vBot*dt
        else
          sLayTopQ(i)=sLayTopQ(i)+sLayBotQ(i+1)
          sLayBotQ(i)=sLayBotQ(i)+vNew(nodeNum)*dt
        endif
        if (lPlant) then
          sLaySink(i)=sLaySink(i)+sumSinkL(N,NLayers,Sink,
     ;                                     nodeTop,nodeNum,dt)
        else
          sLaySink(i)=zero
        endif
        nodeTop=nodeNum-1
      enddo
      return
      end !W_Cumul
c
c     Hydraulic conductivity function: calculate hydraulic
c     conductivity for this pressure head. 
      double precision function FK(iModel,h,Par)
      implicit none
      include "defs.fi"
      double precision h,Par(NTPAR)
      integer iModel
c
      double precision Qr,Qs,Alfa,n,m,Ks,Kr,Kk,Lambda,BPar,
     ;                 HMin,HH,Qees,Qeek,Hs,Hk,Qee,Qek,FFQ,
     ;                 FFQk,t,Qm,Qk,Qa,Qe
      integer PPar
      double precision qnorm
      double precision zero,one
      data zero,one /0.0D00,1.0D00/
c     ENVIRO-GRO 2013 uses iModel=1 for hydraulic properties
c     calculated by the standard van Genuchten-Mualem model.
c     Code for other model is included for future development.
      Qr=Par(1)
      Qs=Par(2)
      Alfa=Par(3)
      n=Par(4)
      Ks=Par(5)
      BPar=Par(6)
      if(iModel.le.1) then
        PPar=2
        Qm=Par(7)
        Qa=Par(8)
        Qk=Par(9)
        Kk=Par(10)
        m=one - one/n
        HMin=-1.D300**(one/n)/max(Alfa,one)
        HH=max(h,HMin)
        Qees=min((Qs-Qa)/(Qm-Qa),.999999999999999D0)
        Qeek=min((Qk-Qa)/(Qm-Qa),Qees)
        Hs=-one/Alfa*(Qees**(-one/m)-one)**(one/n)
        Hk=-one/Alfa*(Qeek**(-one/m)-one)**(one/n)
        if(h.lt.Hk) then
          Qee=(one+(-Alfa*HH)**n)**(-m)
          Qe =(Qm-Qa)/(Qs-Qa)*Qee
          Qek=(Qm-Qa)/(Qs-Qa)*Qeek
          FFQ=one - (one-Qee**(one/m))**m
          FFQk=one - (one-Qeek**(one/m))**m
          if(FFQ.le.zero) FFQ=m*Qee**(one/m)
          Kr=(Qe/Qek)**Bpar*(FFQ/FFQk)**PPar*Kk/Ks
          FK=max(Ks*Kr,1.D-37)
          return
        end if
        if(h.ge.Hk .and. h.lt.Hs) then
          Kr=(one-Kk/Ks)/(Hs-Hk)*(h-Hs)+one
          FK=Ks*Kr
        end if
        if(h.ge.Hs) FK=Ks
      else if(iModel.eq.2) then
        Lambda=2.D0  !  !=2 for Mualem Model, =1.5 for Burdine model
        Hs=-one/Alfa
        if(h.lt.Hs) then
          Kr=one/(-Alfa*h)**(n*(BPar+Lambda)+2.D0)
          FK=max(Ks*Kr,1.D-37)
        else
          FK=Ks
        end if
      else if(iModel.eq.3) then
        Hs=zero
        if(h.lt.Hs) then
          Qee=qnorm(dlog(-h/Alfa)/n)
          t=qnorm(dlog(-h/Alfa)/n+n)
          Kr=Qee**Bpar*t*t
          FK=max(Ks*Kr,1.D-37)
        else
          FK=Ks
        end if
      end if
      return
      end !FK

c     Capacitance function
      double precision function FC(iModel,h,Par)
      implicit none
      include "defs.fi"
      double precision h,Par(NTPAR)
      integer iModel
c
      double precision Qr,Qs,Alfa,n,Qm,Qa,m,HMin,HH,Qees,
     ;                 Hs,C1,C2,t
      double precision zero,one,two,pi
      data zero,one,two,pi /0.0D00,1.0D00,2.0D00,
     ;                      3.1415926535897932D00/
      Qr=Par(1)
      Qs=Par(2)
      Alfa=Par(3)
      n=Par(4)
      if(iModel.le.1) then
        Qm=Par(7)
        Qa=Par(8)
        m=one-one/n
        HMin=-1.D300**(one/n)/max(Alfa,one)
        HH=max(h,HMin)
        Qees=min((Qs-Qa)/(Qm-Qa),.999999999999999D00)
        Hs=-one/Alfa*(Qees**(-one/m)-one)**(one/n)
        if(h.lt.Hs) then
          C1=(one + (-Alfa*HH)**n)**(-m-one)
          C2=(Qm-Qa)*m*n*(Alfa**n)*(-HH)**(n-one)*C1
          FC=max(C2,1.D-37)
          return
        else
          FC=zero
        end if
      else if(iModel.eq.2) then
        Hs=-one/Alfa
        if(h.lt.Hs) then
          C2=(Qs-Qr)*n*Alfa**(-n)*(-h)**(-n-one)
          FC=max(C2,1.D-37)
        else
          FC=zero
        end if
      else if(iModel.eq.3) then
        Hs=zero
        if(h.lt.Hs) then
          t=exp(-one*(dlog(-h/Alfa))**two/(two*n**two))
          C2=(Qs-Qr)/(two*pi)**0.5d0/n/(-h)*t
          FC=max(C2,1.D-37)
        else
          FC=zero
        end if
      end if
      return
      end ! FC
c
c     Water content from pressure head
      double precision function FQ(iModel,h,Par)
      implicit none
      include "defs.fi"
      double precision h,Par(NTPar)
      integer iModel
c
      double precision Qr,Qs,Alfa,n,Qm,Qa,m,HMin,HH,Qees,Hs,
     ;                 Qee,zero,one
      data zero,one /0.0D00,1.0D00/
      double precision qnorm
      Qr=Par(1)
      Qs=Par(2)
      Alfa=Par(3)
      n=Par(4)
      if(iModel.le.1) then
        Qm=Par(7)
        Qa=Par(8)
        m=one-one/n
        HMin=-1.d300**(one/n)/max(Alfa,one)
        HH=max(h,HMin)
        Qees=min((Qs-Qa)/(Qm-Qa),.999999999999999D00)
        Hs=-one/Alfa*(Qees**(-one/m)-one)**(one/n)
        if(h.lt.Hs) then
          Qee=(one+(-Alfa*HH)**n)**(-m)
          FQ=max(Qa+(Qm-Qa)*Qee,1.d-37)
          return
        else
          FQ=Qs
        end if
      else if(iModel.eq.2) then
        Hs=-one/Alfa
        if(h.lt.Hs) then
          Qee=(-Alfa*h)**(-n)
          FQ=max(Qr+(Qs-Qr)*Qee,1.D-37)
        else
          FQ=Qs
        end if
      else if(iModel.eq.3) then
        Hs=zero
        if(h.lt.Hs) then
          Qee=qnorm(dlog(-h/Alfa)/n)
          FQ=max(Qr+(Qs-Qr)*Qee,1.d-37)
        else
          FQ=Qs
        end if
      end if
      return
      end !FQ

c     Pressure head from theta
      double precision function FH(iModel,theta,ParD,iMat)
      implicit none
      include "defs.fi"
      integer iModel,iMat
      double precision theta,ParD(NTPar,*)
c
      double precision Par(nTPar),Qr,Qs,Alfa,n,Qe,Qm,Qa,
     ;                 m,HMin,QeeM,Qee,y,x,p,zero,one,two
      data zero,one,two /0.0D00,1.0D00,2.0D00/
c
      call layerPar(ParD,Par,nTPar,iMat)
      Qr=Par(1)
      Qs=Par(2)
      Alfa=Par(3)
      n=Par(4)
c     scaled water content
      Qe = max((theta-Qr)/(Qs-Qr),zero)
      Qe = min(Qe,one)
      if(iModel.le.1) then
        Qm=Par(7)
        Qa=Par(8)
        m=one-one/n
        HMin=-1.d300**(one/n)/max(Alfa,one)
        QeeM=(one+(-Alfa*HMin)**n)**(-m)
        Qee=min(max(Qe*(Qs-Qa)/(Qm-Qa),QeeM),one)
        FH=max(-one/Alfa*(Qee**(-one/m)-one)**(one/n),-1.D37)
      else if(iModel.eq.2) then
        FH=dble(max(-one/Alfa*Qe**(-one/n),-1.d37))
      else if(iModel.eq.3) then !not implemented
        if(Qe.gt.0.9999d0) then
          FH=zero
        elseif(Qe.lt.1.0D-05) then
          FH=-1.D+06
        else
          y=Qe*two
          if(y.lt.one) p=sqrt(-dlog(y/two))
          if(y.ge.one) p=sqrt(-dlog(one-y/two))
          x=p-(1.881796+0.9425908*p+0.0546028*p**3)/
     ;       (1.+2.356868*p+0.3087091*p**2+0.0937563*p**3+0.021914*p**4)
          if(y.ge.one) x=-x
          FH=-Alfa*exp(sqrt(two)*n*x)
        endif
      endif
      return
      end !FH

      double precision function qnorm(x)
      implicit none
      double precision x
c
      double precision z,t,erfc,zero,one,two,half
      data zero,one,two,half /0.0D00,1.0D00,2.0D00,5.0D-01/
      z=abs(x/two**half)
      t=one/(one+half*z)
      erfc=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+
     ;     t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+
     ;     t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))))
      if(x.lt.zero) erfc=two-erfc
      qnorm=erfc/two
      return
      end
c
c     Copy variables for water flow from "new" to "old" vectors
      subroutine copyWat(N,hOld,hNew,hTemp,thOld,thNew,vOld,vNew,
     ;                   dt,dtOld)
      implicit none
      integer N
      double precision hOld,hNew,hTemp,thOld,thNew,vOld,vNew,dt,
     ;                 dtOld
      dimension hOld(N),hNew(N),thOld(N),thNew(N),vOld(N),vNew(N),
     ;          hTemp(N)
c
      integer i
      do i=1,N
        hTemp(i)=hNew(i)+(hNew(i)-hOld(i))*dt/dtOld
        hOld(i) =hNew(i)
        hNew(i) =hTemp(i)
        thOld(i)=thNew(i)
        vOld(i) =vNew(i)
      enddo
      return
      end