      subroutine vrtslv(nlay,ii,jj,igrd,dtin,entrn,dilut,depth,conc,
     &                   lvupsolv,fluxtop,fluxlay,sens,spec,fc1,fc2,
     &                                                    fc3,ldoipts)
      use filunit
      use tracer
c
c----CAMx v7Beta6 190902
c 
c     VRTSLV performs column mass adjustments due to vertical advection and
c     changes in layer structure that result from the time- and space-varying
c     vertical coordinate system.  This version performs the mass adjustments
c     for both the concentrations and the sensitivities.
c     The system is solved using a fully implicit backward Euler approach in
c     time and a hybrid centered/upstream-donor approach in space.  The latter
c     (upstream) is used if the vertical CFL>1, or if centered difference results in
c     a local negative concentrations locally.
c
c     The difference equations are generalized to allow for three time 
c     differencing options, depending on the value of the parameter "mu":
c        1) Crank-Nicholson semi-implicit (mu=0.5)
c        2) Fully explicit forward Euler (mu=0)
c        3) Fully implicit backward Euler (mu=1) !Current implementation!
c     Upper boundary conditions:
c        This routine accomodates non-zero vertical velocities at the top of
c        the column; the top boundary condition is held in the conc vector
c        at location (nlay+1).
c     Lower boundary conditions:
c        Zero flux is specified at the ground
c
c     Copyright 1996 - 2018
c     Ramboll
c
c     Modifications: 
c        12/3/99   This routine is now fully implicit, and calculation
c                  of sub-timesteps has been removed.
c        3/06/06   Revised top BC approach; removed top dummy layer
c       07/16/07 -bkoo-     Added check for HDDM
c       07/16/08 -bkoo-     Added DDM turn-off flag
c       11/04/09   Revised for hybrid centered/upstream-donor technique
c
c     Input arguments: 
c        nlay                number of layers
c        ii                  column index
c        jj                  row index
c        dtin                timestep (s)
c        entrn               entrainment rate (m/s) 
c        dilut               dilution rate (m/s)
c        depth               layer depth (m)
c        conc                species concentration (umol/m3)
c        lvupsolv            flag to use upstream donor solver
c        sens                sensitivity coefficients(umol/m3/parameter unit)
c        spec                species name
c        ldoipts             flag to calculate Process Analysis data
c 
c     Output arguments: 
c        conc                species concentration (umol/m3)
c        fluxtop             flux across the top boundary (umol/(m2*s))
c        fluxlay             flux across top of each layer (umol/(m2*s))
c        sens                sensitivity coefficients(umol/m3/parameter unit)
c        fc1                 conc entrained through bottom of layer (umol/m3)
c        fc2                 conc entrained through top of layer (umol/m3)
c        fc3                 conc diluted by layer expansion/contraction
c                            (umol/m3)
c 
c     Routines Called: 
c        TRIDIAG
c 
c     Called by: 
c        ZADVEC 
c
      implicit none
      include "camx.prm"
c
      real mu
      parameter (mu = 1.0)

      character*10 spec
      integer nlay
      integer ii,jj,igrd
      real conc(nlay+1),entrn(nlay),dilut(nlay),depth(nlay)
c
      real dtin,deltac
      real fluxtop
c
      integer nsteps,k,kk,ineg
      real dt
      logical lvupsolv(nlay)
c
      real pn(MXLAYER)
      real pnp1(MXLAYER)
      real aa(MXLAYER)
      real bb(MXLAYER)
      real cc(MXLAYER)
      real rr(MXLAYER*(MXTRSP+1))
      real zrp(MXLAYER)
      real zrm(MXLAYER)
c
c===========================DDM Begin=================================
c
      integer ls,lrr,isen
      real sens((nlay+1)*nddmsp)
c
c     Note:  nddmsp is the total number of DDM parameters for which
c            sensitivities are calculated.
c
c===========================DDM End===================================
c
c
c================== Source Apportionmment Begin ======================
c
      real fluxlay(MXLAYER)
c
c=================== Source Apportionment End ========================
c
c
c================== Process Analysis Begin ===========================
c
      logical ldoipts
      real fc1(*)
      real fc2(*)
      real fc3(*)
c
c=================== Process Analysis End ============================
c
c
c-----Entry point
c
      dt = dtin
      nsteps = 1
c
c-----Calculate some constants
c
      fluxtop = 0.
      do k = 1,nlay
        pn(k) = (1. - mu)*dt/depth(k)
        pnp1(k) = mu*dt/depth(k)
      enddo
      do k = 1,nlay-2
        zrp(k) = depth(k+1)/(depth(k+1) + depth(k))
        zrm(k) = depth(k)/(depth(k+1) + depth(k))
      enddo
      do k = nlay-1,nlay
        zrp(k) = 0.
        zrm(k) = 1.
        if (entrn(k).lt.0.) then
          zrp(k) = 1.
          zrm(k) = 0.
        endif
      enddo
c
c-----Check upstream donor flag for hybrid case
c
      do k = 1,nlay-2
        if (lvupsolv(k)) then
          zrp(k) = 0.
          zrm(k) = 1.
          if (entrn(k).lt.0.) then
            zrp(k) = 1.
            zrm(k) = 0.
          endif
        endif
      enddo
c
c-----Calculate matrix coefficients A,B & C, and array R for the equation
c        n+1     n+1     n+1   n
c     A*c   + B*c   + C*c   = R
c        k-1     k       k+1
c
 100  continue
c
      do k = 2,nlay-1
        aa(k) = pnp1(k)*zrp(k-1)*entrn(k-1)
        bb(k) = 1. + pnp1(k)*
     &          (zrm(k-1)*entrn(k-1) - zrp(k)*entrn(k) + dilut(k))
        cc(k) = -pnp1(k)*zrm(k)*entrn(k)
        rr(k) = -pn(k)*zrp(k-1)*entrn(k-1)*conc(k-1) +
     &          conc(k)*(1. - pn(k)*
     &          (zrm(k-1)*entrn(k-1) - zrp(k)*entrn(k) + dilut(k))) +
     &          pn(k)*zrm(k)*entrn(k)*conc(k+1)
      enddo
c
c-----Lower boundary conditions
c
      k = 1
      aa(k) = 0.
      bb(k) = 1. + pnp1(k)*(-zrp(k)*entrn(k) + dilut(k))
      cc(k) = -pnp1(k)*zrm(k)*entrn(k)
      rr(k) = conc(k)*(1. - pn(k)*(-zrp(k)*entrn(k) + dilut(k))) +
     &        pn(k)*zrm(k)*entrn(k)*conc(k+1)
c
c-----Upper boundary conditions
c
      k = nlay
      aa(k) = pnp1(k)*zrp(k-1)*entrn(k-1)
      bb(k) = 1. + pnp1(k)*
     &        (zrm(k-1)*entrn(k-1) - zrp(k)*entrn(k) + dilut(k))
      cc(k) = 0.
      rr(k) = -pn(k)*zrp(k-1)*entrn(k-1)*conc(k-1) +
     &        conc(k)*(1. - pn(k)*
     &        (zrm(k-1)*entrn(k-1) - zrp(k)*entrn(k) + dilut(k))) +
     &        zrm(k)*entrn(k)*conc(k+1)*dt/depth(k)
c
c============================DDM Begin================================
c
c-----If sensitivities are requested, load additional arrays into rr,
c     one array for each sensitivity coefficient.  Additional right-hand
c     side arrays rr are calculated for the equation
c        n+1     n+1     n+1    n
c     A*s   + B*s   + C*s   = r
c        k-1     k       k+1
c     The right-hand side arrays for both the concentration and
c     sensitivity equations are loaded into the same array rr, and
c     the equations are solved together.
c
c     Load the lower boundary condition, then layers
c     2 through nlay-1, finally the upper boundary condition.
c
        if( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then
          ls = 0
          lrr = nlay
          do isen = 1,nddmsp
            ls = ls + 1
            lrr = lrr + 1
            k = 1
            rr(lrr) = sens(ls)*
     &                (1. - pn(k)*(-zrp(k)*entrn(k) + dilut(k))) +
     &                pn(k)*zrm(k)*entrn(k)*sens(ls+1)
            do k = 2,nlay-1
              ls = ls + 1
              lrr = lrr + 1
              rr(lrr) = -pn(k)*zrp(k-1)*entrn(k-1)*sens(ls-1) +
     &                  sens(ls)*(1. - pn(k)*
     &                           (zrm(k-1)*entrn(k-1) -
     &                            zrp(k)*entrn(k) +
     &                            dilut(k))) +
     &                  pn(k)*zrm(k)*entrn(k)*sens(ls+1)
            enddo
 
            ls = ls + 1
            lrr = lrr + 1
            k = nlay
            rr(lrr) = -pn(k)*zrp(k-1)*entrn(k-1)*sens(ls-1) +
     &                sens(ls)*(1. - pn(k)*
     &                         (zrm(k-1)*entrn(k-1) - 
     &                          zrp(k)*entrn(k) + 
     &                          dilut(k))) +
     &                zrm(k)*entrn(k)*sens(ls+1)*dt/depth(k)
            ls = ls + 1
          enddo
        endif
c
c=============================DDM End=================================
c
c-----Solve the equations
c
      if( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then
         call tridiag(aa,bb,cc,rr,nlay,nddmsp+1)
      else
         call tridiag(aa,bb,cc,rr,nlay,1)
      endif
c
      if (nsteps.lt.nlay) then
        ineg = 0
        do k = 1,nlay
          deltac = (rr(k) - conc(k))/conc(k)
          if (deltac.le.-0.99) then
            if( k.gt.1 ) then
              if( entrn(k-1).gt.0.) then
                zrp(k-1) = 0.
                zrm(k-1) = 1.
              endif
            endif
            if (entrn(k).lt.0.) then
              zrp(k) = 1.
              zrm(k) = 0.
            endif
            ineg = 1
          endif
        enddo
        if (ineg.eq.1) then
          nsteps = nsteps + 1
          goto 100
        endif
      endif
c
c===================== Source Apportionment Begin =====================
c
      fluxlay(1) = -(rr(1) - conc(1))*depth(1)/dt
      do k = 2,nlay
        fluxlay(k) = -(rr(k) - conc(k))*depth(k)/dt + fluxlay(k-1)
      enddo
c
c====================== Source Apportionment End ======================
c
c
c===================== Process Analysis Begin =========================
c
      if (ldoipts) then
        do k = 1,nlay
          fc3(k) = -dt/depth(k)*dilut(k)*conc(k)
        enddo
        fc1(1) = 0.
        do k = 2,nlay
          fc1(k) = fluxlay(k-1)/depth(k)*dt
        enddo
        do k = 1,nlay
          fc2(k) = -fluxlay(k)/depth(k)*dt
        enddo
      endif
c
c====================== Process Analysis End ==========================
c
      fluxtop =  fluxlay(nlay)
      do k = 1,nlay
        if (rr(k).le.0.) then 
          write(iout,'(//,a)') 'ERROR in VRTSLV:'
          write(iout,*) 'Negative concentration ',
     &                  'when doing advection in z-direction'
          write(iout,*) 'Grid: ',igrd
          write(iout,*) 'Location (I,J,K): ',ii,jj,k
          write(iout,*) 'Species: ',spec 
          write(iout,'(a,i3,a,i3)')'Number of steps: ',nsteps
          do kk = 1,nlay 
            write(iout,'(i3,6e12.4,f5.1,6f7.3)') kk,rr(kk), 
     &            conc(kk),depth(kk),entrn(kk),dilut(kk)
          enddo 
          call camxerr()
        endif
        conc(k) = rr(k)
      enddo
c
c==============================DDM Begin===============================
c
      if( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then
        ls = -1
        lrr = nlay
        do isen = 1,nddmsp
          ls = ls + 1
          do k = 1, nlay
            ls = ls + 1
            lrr = lrr + 1
            sens(ls) = rr(lrr)
          enddo
        enddo
      endif
c
c===============================DDM End================================
c
      return
      end
