c     Assemble and solve the solute transport equation
c     Finite differences
c
************************************************************************
      subroutine TrnsPort(N,NMat,volConv,t,dt,z,MatNum,chemPar,thOld,
     ;                    thNew,vOld,vNew,Dold,DNew,cTop,cBot,iRoot,
     ;                    chmBCTop,chmBCBot,Conc,TConc,CFix,NLossRat,
     ;                    B,D,E,F,cvTop,cvBot,Peclet,Courant,dtMaxC,
     ;                    RConc,Sink,Psi,MaxCh,tPotMax,tpotN,tPotNMax,
     ;                    tMassR,alphaG,Beta,sConvN,prodNSea,TVarCoef,
     ;                    lTCoef,dateInit,cStarN,sumDel,cropProp,
     ;                    crpIndex,nConc,cropNum,compCnt,lNitrate,ierr)
      implicit none
      include "defs.fi"
      integer       N,NMat,nConc,chmBCTop,chmBCBot,MaxCh,compCnt,
     ;              MatNum,iRoot,crpIndex,cropNum,ierr
      logical       lNitrate,lTCoef
      character*10  dateInit
      double precision dt,z,dz,chemPar,thOld,thNew,vOld,vNew,Dold,DNew,
     ;                 cTop,cBot,Conc,B,D,E,F,cvTop,cvBot,Peclet,
     ;                 Courant,dtMaxC,RConc,CFix,tPotNMax,
     ;                 TConc,tPotMax,tpotN,tMassR,volConv,prodNSea,
     ;                 alphaG,NLossRat,sConvN,t,Sink,cStarN,sumDel,
     ;                 cropProp,TVarCoef,Beta,Psi
      dimension z(N),MatNum(N),chemPar(nChemPar,NMat),
     ;          thOld(N),thNew(N),vOld(N),vNew(N),Dold(N),DNew(N),
     ;          Conc(nConcD,N),CFix(nFixC,N),TConc(nConcD,N),
     ;          cBot(nIrriC),cTop(nIrriC),cvTop(nIrriC),cvBot(nIrriC),
     ;          Beta(N),RConc(nConcD,N),B(N),D(N),E(N),F(N),Sink(N),
     ;          Psi(N),cropProp(nCropPar,*),TVarCoef(nMonths)
      logical   lArtD,trdone,lLagr
      integer   i,j,k,m
      double precision wc,alfa,CourMax,DD,Disp,Vel,Pec,vMax,dtMax,
     ;                 DDO,DDN,wc1,wcN,BN,FN,E1,F1,D1,a1,b1,Cour,thP
      dimension wc(NumNPD),thP(N)
      character*80 errStr
      data wc / NumNPD*0.0D00 /
      double precision zero,one,two,six,epsi,x,y,err,epserr,del,z1
      data    zero,one,two,six /0.0D+00,1.0D+00,2.0D+00,6.0D+00/
c     Use standard Crank-Nicholson
      data  epsi /0.5D+00/
c
      err(x,y)=dabs(1.-x/y)-1.0D-03
      epserr=1.0D-06
c
c     Eulerian-Lagrangian approach with reverse particle tracking
c     when lLagr is set to true. Useful for convection-dominated
c     problems
      lLagr=.FALSE.
c     No artificial dispersion
      lArtD=.FALSE.
      alfa=one-epsi
      Peclet=zero
      Courant=zero
      CourMax=one
      dtMaxC=1.0D+30
      do 12 i=N,1,-1
        M=MatNum(i)
c       chemPar(3,M) Dmolecular
c       chemPar(4,M) Dispersivity
        Dold(i)=chemPar(3,M)+chemPar(4,M)*dabs(vOld(i))/thOld(i)
        DNew(i)=chemPar(3,M)+chemPar(4,M)*dabs(vNew(i))/thNew(i)
        if(lArtD) then
          DD=vOld(i)*vOld(i)*dt/thOld(i)/thOld(i)/two
          if(DD.gt.Dold(i)) Dold(i)=DD
          DD=vNew(i)*vNew(i)*dt/thNew(i)/thNew(i)/two
          if(DD.gt.DNew(i)) DNew(i)=DD
        end if
        if(i.ne.N) then
          Disp=(Dold(i)+DDO+DNew(i)+DDN)/4.0D+00
          Vel=(vOld(i)/thOld(i)+vOld(i+1)/thOld(i+1) +
     ;         vNew(i)/thNew(i)+vNew(i+1)/thNew(i+1))/4.0D+00
          dz=z(i+1)-z(i)
c         Calculate Peclet and Courant numbers
          Pec=9.9999D+04
          dtMax=1.D+30
          vMax=dmax1(dabs(vNew(i))/thNew(i),dabs(vNew(i+1))/thNew(i+1))
          if(Disp.gt.zero) Pec=dabs(Vel)*dz/Disp
          Cour=vMax*dt/dz
          Peclet=dmax1(Peclet,Pec)
          Courant=dmax1(Courant,Cour)
          if(vMax.gt.zero) dtMax=CourMax*dz/vMax
          dtMaxC=dmin1(dtMaxC,dtMax)
        end if
        DDO=Dold(i)
        DDN=DNew(i)
        if(.not.lArtD) then
          Dold(i)=Dold(i)+vOld(i)**2*dt/six/thOld(i)**2
          DNew(i)=DNew(i)-vNew(i)**2*dt/six/thNew(i)**2
        end if
12    continue

c     Construct the left hand side of the matrix equation
      if(lLagr) then
        wc1=zero
        wcN=zero
      else
        wc1=wc(1)
        wcN=wc(N-1)
      end if
      if(lLagr) then
        k=1
        do 14 i=2,N-1
          z1=z(i)-vNew(i)/thNew(i)*dt
          if(z1.le.z(1)) then
            thP(i)=(alfa*thOld(1)+epsi*thNew(1))/two
          else if(z1.ge.z(n)) then
            thP(i)=(alfa*thOld(N)+epsi*thNew(N))/two
          else  
13          if(z1.lt.z(k+1)) then
              thP(i)=thOld(k)+(thOld(k+1)-thOld(k))*
     ;               (z1-z(k))/(z(k+1)-z(k))
              thP(i)=alfa*thP(i)+epsi*thNew(i)
            else
              k=k+1
              goto 13
            end if
          end if
14      continue
      end if
  
      b1=z(2)-z(1)
      D1=b1/two/dt*thNew(1)+epsi*(thNew(1)*DNew(1)+
     ;                   thNew(2)*DNew(2))/two/b1+
     ;                   epsi*((two+3.0D+00*wc1)*vNew(1)+
     ;                   vNew(2))/six
      E1=-epsi*(thNew(1)*DNew(1)+thNew(2)*DNew(2))/two/b1+
     ;    epsi*(vNew(1)+(two-3.0D+00*wc1)*vNew(2))/six

      if(chmBCBot.ge.1) then
        D(1)=one
        E(1)=zero
      else if(chmBCBot.le.-1) then
        E(1)=E1
        D(1)=D1
      else if(chmBCBot.eq.0) then
        D(1)=-one
        E(1)=one
      end if
      do i=2,N-1
        a1=b1
        b1=z(i+1)-z(i)
        dz=(z(i+1)-z(i-1))/two
        B(i)=-epsi*(thNew(i-1)*DNew(i-1)+thNew(i)*DNew(i))/two/a1
        D(i)= epsi*(thNew(i-1)*DNew(i-1)+thNew(i)*DNew(i))/two/a1+
     ;        epsi*(thNew(i+1)*DNew(i+1)+thNew(i)*DNew(i))/two/b1
        E(i)=-epsi*(thNew(i+1)*DNew(i+1)+thNew(i)*DNew(i))/two/b1
c
************************************************************************
        if(lLagr) then
          D(i)=D(i)+dz/dt*thP(i)-epsi*Sink(i)*dz
        else
          B(i)=B(i)-epsi*((two+3.0D+00*wc(i-1))*vNew(i-1)+vNew(i))/six
          D(i)=D(i)+epsi*(vNew(i+1)+
     ;         3.0D+00*(wc(i-1)+wc(i))*vNew(i)-vNew(i-1))/six+
     ;         dz/dt*thNew(i)
          E(i)=E(i)+epsi*(vNew(i)+(two-3.0D+00*wc(i))*vNew(i+1))/six
        end if
      enddo
      BN=-epsi*(thNew(N-1)*DNew(N-1)+thNew(N)*DNew(N))/b1/two-
     ;    epsi*((two+3.0D+00*wcN)*vNew(N-1)+vNew(N))/six
      DD=b1/two/dt*thNew(N)-
     ;   epsi*(vNew(N-1)+(two-3.0D+00*wcN)*vNew(N))/six+
     ;   epsi*(thNew(N-1)*DNew(N-1)+thNew(N)*DNew(N))/two/b1
      if(chmBCTop.gt.0) then
        B(N)=zero
        D(N)=one
      else
        B(N)=BN
        D(N)=DD
      end if
      do i=2,N
        D(i)=D(i)-B(i)*E(i-1)/D(i-1)
      enddo
c
c     Begin main transport calculation
c     Number of solute species is nConc
c
      j=1
      trdone=.FALSE.
c     top of the loop for species, construct right hand side of
c     the matrix equation 
 
301   if (trdone) then
        go to 23
      else
        b1=z(2)-z(1)
************************************************************************
        F1=Conc(j,1)*(alfa*(-(thOld(1)*Dold(1)+thOld(2)*Dold(2))/b1/two-
     ;                      ((two+3.0D+00*wc1)*vOld(1)+vOld(2))/six)+
     ;                b1/two/dt*thOld(1))+
     ;     Conc(j,2)*alfa*((thOld(1)*Dold(1)+thOld(2)*Dold(2))/b1/two-
     ;                    (vOld(1)+(two-3.0D+00*wc1)*vOld(2))/six)
        if(chmBCBot.ge.1)  F(1)=cBot(j)
        if(chmBCBot.le.-1) F(1)=F1+cBot(j)*(alfa*vOld(1)+epsi*vNew(1))
        if(chmBCBot.eq.0)  F(1)=zero
        do i=2,N-1
          a1=b1
          b1=z(i+1)-z(i)
          dz=(z(i+1)-z(i-1))/two
************************************************************************
          F(i)=Conc(j,i-1)*alfa*(thOld(i-1)*Dold(i-1)+
     ;                                    thOld(i)*Dold(i))/two/a1-
     ;         Conc(j,i) *(alfa*(thOld(i-1)*Dold(i-1)+
     ;                                    thOld(i)*Dold(i))/two/a1+
     ;           alfa*(thOld(i+1)*Dold(i+1)+thOld(i)*Dold(i))/two/b1)+
     ;         Conc(j,i+1)*alfa*(thOld(i+1)*Dold(i+1)+
     ;                                    thOld(i)*Dold(i))/two/b1
          if(lLagr) then
            F(i)=F(i)+wc(i)*thP(i)*dz/dt+Conc(j,i)*alfa*Sink(i)*dz
          else
            F(i)=F(i)+Conc(j,i-1)*alfa*((two+3.0D+00*wc(i-1))*vOld(i-1)+
     ;           vOld(i))/six+Conc(j,i)*(-alfa*(vOld(i+1)+
     ;           3.0D+00*(wc(i-1)+wc(i))*vNew(i)-vOld(i-1))/six+
     ;           dz/dt*thOld(i))-Conc(j,i+1)*alfa*(vOld(i)+
     ;           (two-3.0D+00*wc(i))*vOld(i+1))/six
          end if
        enddo
        FN=Conc(j,N-1)*alfa*((thOld(N-1)*Dold(N-1)+
     ;                     thOld(N)*Dold(N))/b1/two+
     ;                     ((two+3.0D+00*wcN)*vOld(N-1)+vOld(N))/six)+
     ;     Conc(j,N)*(alfa*((vOld(N-1)+(two-3.0D+00*wcN)*vOld(N))/six-
     ;                (thOld(N-1)*Dold(N-1)+thOld(N)*Dold(N))/two/b1)+
     ;                b1/two/dt*thOld(N))
        if(chmBCTop.gt.0) then
          F(N)=cTop(j)
        else
          F(N)=FN
          if(vOld(N).lt.zero) F(N)=F(N)-vOld(N)*cTop(j)*alfa
          if(vNew(N).lt.zero) F(N)=F(N)-vNew(N)*cTop(j)*epsi
        end if

c       Solve matrix equation
        do i=2,N
          F(i)=F(i)-B(i)*F(i-1)/D(i-1)
        enddo
        RConc(j,N)=F(N)/D(N)
        do i=2,N
          k=N-i+1
          RConc(j,k)=(F(k)-E(k)*RConc(j,k+1))/D(k)
        enddo
c
c       Set up mass flux
        if(chmBCTop.lt.0) then
          cvTop(j)=zero
          if(vOld(N).lt.zero) cvTop(j)=vOld(N)*alfa*cTop(j)
          if(vNew(N).lt.zero) cvTop(j)=cvTop(j)+vNew(N)*epsi*cTop(j)
        else
          cvTop(j)=FN-BN*RConc(j,N-1)-DD*RConc(j,N)
        end if

        if(chmBCBot.lt.0) then
          cvBot(j)=cBot(j)*(vOld(1)*alfa+vNew(1)*epsi)
        else
          cvBot(j)=D1*RConc(j,1)+E1*RConc(j,2)-F1
        endif
c
c       Increment counter for species and go to top of loop
c
        j=j+1
        if (.not.lNitrate.and.j.eq.3) then
          trDone=.TRUE.
        elseif (j.eq.4) then
          trDone=.TRUE.
        endif
      end if
      go to 301
23    continue
c     Before calling the Nitrogen sink copy the new RConc
c     matrix into TConc for storage and convergence testing
      call Mat1(N,RConc,nConcD,TConc)
************************************************************************
      if (lNitrate) then
        call setSrcN(N,NMat,volConv,z,MatNum,t,dt,CFix,nFixC,RConc,
     ;             nConcD,chemPar,nChemPar,Psi,thNew,NLossRat,sConvN,
     ;             TVarCoef,lTCoef,dateInit,prodNSea,nMonths)
        if (cropNum.gt.0) then
          call setSnkN(N,iRoot,volConv,t,dt,z,tpotN,tPotMax,tPotNMax,
     ;                 RConc,tMassR,thNew,Beta,alphaG,cStarN,sumDel,
     ;                 cropProp,crpIndex,ierr,compCnt)
c         Iterations in compensation routine exceeded 50 note in
c         error log and proceed
          if (ierr.gt.0) then
            call countErr(nErrLog,ierr,errStr,compCnt)
            ierr=0
          endif
        endif
      endif
      call Mat1(N,RConc,nConcD,Conc)
901   return
      end
c
c     Calculate cumulative mass flux for each chemical species
c     at the top and bottom of the column
      subroutine C_Cumul(N,NLayers,ldc,z,volConv,dt,vTop,vBot,cTop,
     ;                   Conc,theta,sumTopQ,sumBotQ,tMassR,tpotN,
     ;                   tPotNMax,sumRootN,sMassN,sPotN,initFlg,
     ;                   initNO3,layerDep,cvBot,leachN,sLeachN,
     ;                   tLeachN,stSeason,freeD,lPlant,lNitrate)
      implicit none
      integer N,NLayers,ldc
      logical initFlg,initNO3,stSeason,freeD,lPlant,lNitrate
      double precision z,dt,vTop,vBot,cTop,sumTopQ,sumBotQ,theta,
     ;                 tMassR,sumRootN,sMassN,sPotN,tpotN,volConv,
     ;                 cvBot,leachN,sLeachN,Conc,sMaxN,layerDep,
     ;                 tLeachN,tPotNMax
      dimension        z(N),cTop(ldc),sumTopQ(ldc),Conc(ldc,N),
     ;                 sumBotQ(ldc),cvBot(ldc),theta(N),
     ;                 layerDep(NLayers)
c
      integer j
      double precision zero,two,sinkSum,dz,MWN,amtLeach,massNLay
      dimension massNLay(NLayers)
      data zero,two,MWN /0.0D00, 2.0D00, 1.4D01/

      if (initFlg) then
        do j=1,ldc
          sumTopQ(j)=zero
          sumBotQ(j)=zero
        enddo
      endif
      if (initNO3) then
c       Daily initialization of leachN
        leachN=zero
        initNO3=.FALSE.
      else
        sMassN=sMassN+tMassR
        sPotN=sPotN+tPotNMax*dt !sum of maximum potential N uptake
      endif
      do j=1,ldc
        sumTopQ(j)=sumTopQ(j) + cTop(j)*volConv*vTop*dt !mmol/cm^2
        if (freeD) then
          sumBotQ(j)=sumBotQ(j) + cvBot(j)*volConv*dt !mmol/cm^2
        else
          sumBotQ(j)=sumBotQ(j) + Conc(j,1)*volConv*vBot*dt
        endif
      enddo
c     cvBot is (mmol/L)*(cm3/cm2)/d must convert the L to cm^3 to get
c     flux units of mmol/cm^2/d volConv=0.001 then convert to kg/ha/d
      if (lNitrate) then
        amtLeach=cvBot(3)*volConv*MWN*dt*1.0D02
        leachN=leachN + amtLeach
c       Total leached N
        tLeachN = tLeachN + amtLeach
        if (stSeason) then
c         Reset the seasonal sum at the start of the season
          sLeachN=amtLeach
          stSeason=.False.
        else
c         seasonal sum of leached N kg/ha
          sLeachN=sLeachN+amtLeach
        endif
      else
        leachN=zero
      endif
c
c     add to the cumulative root nitrogen flux needed for mass
c     balance in chemBal routine
      if (lPlant .and. lNitrate) then
        sumRootN=sumRootN + tMassR
      endif
      return
      end


************************************************************************
      subroutine Mat1(n,x1,ldc,x2)
      implicit none
      double precision x1,x2
      integer n,ldc
      dimension x1(ldc,n),x2(ldc,n)
c
      integer i,j
c
      do 12 i=1,ldc
        do 11 j=1,n
          x2(i,j)=x1(i,j)
11      continue
12    continue
      return
      end 
