subroutine samfdeepcnv_aerosols(im, ix, km, itc, ntc, ntr, delt,
     &  xlamde, xlamdd, cnvflg, jmin, kb, kmax, kbcon, ktcon, fscav,
     &  edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp,
     &  qtr, qaero)

      use machine , only : kind_phys
      use physcons, only : g => con_g, qamin

      implicit none

c     -- input arguments
      integer,                                      intent(in) :: im,
     &  ix, km, itc, ntc, ntr
      real(kind=kind_phys),                         intent(in) :: delt,
     &  xlamde, xlamdd
      logical,              dimension(im),          intent(in) :: cnvflg
      integer,              dimension(im),          intent(in) :: jmin,
     &  kb, kmax, kbcon, ktcon
      real(kind=kind_phys), dimension(im),          intent(in) :: edto,
     &  xlamd, xmb
      real(kind=kind_phys), dimension(ntc),         intent(in) :: fscav
      real(kind=kind_phys), dimension(im,km),       intent(in) :: c0t,
     &  eta, etad, zi, xlamue, xlamud
      real(kind=kind_phys), dimension(ix,km),       intent(in) :: delp
      real(kind=kind_phys), dimension(ix,km,ntr+2), intent(in) :: qtr
c     -- output arguments
      real(kind=kind_phys), dimension(im,km,ntc),   intent(out) :: qaero

c     -- local variables
c     -- general variables
      integer :: i, indx, it, k, kk, km1, kp1, n
      real(kind=kind_phys) :: adw, aup, dtime_max, dv1q, dv2q, dv3q,
     &  dtovdz, dz, factor, ptem, ptem1, qamax, tem, tem1
      real(kind=kind_phys), dimension(ix,km) :: xmbp
c     -- chemical transport variables
      real(kind=kind_phys), dimension(im,km,ntc) :: ctro2, ecko2, ecdo2,
     &  dellae2
c     -- additional variables for tracers for wet deposition,
      real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw,
     &  wet_dep
c     -- if reevaporation is enabled, uncomment lines below
c     real(kind=kind_phys), dimension(im,ntc) :: pwav
c     real(kind=kind_phys), dimension(im,km) :: pwdper
c     real(kind=kind_phys), dimension(im,km,ntr) :: chem_pwd
c     -- additional variables for fct
      real(kind=kind_phys), dimension(im,km) :: flx_lo, totlout, clipout

      real(kind=kind_phys), parameter :: one     = 1.0_kind_phys
      real(kind=kind_phys), parameter :: half    = 0.5_kind_phys
      real(kind=kind_phys), parameter :: quarter = 0.25_kind_phys
      real(kind=kind_phys), parameter :: zero    = 0.0_kind_phys
      real(kind=kind_phys), parameter :: epsil   = 1.e-22_kind_phys    ! prevent division by zero

c     -- begin

c     -- check if aerosols are present
      if ( ntc <= 0 .or. itc <= 0 .or. ntr <= 0 ) return
      if ( ntr < itc + ntc - 3 ) return

c     -- initialize work variables
      km1 = km - 1

      chem_c    = zero
      chem_pw   = zero
      ctro2     = zero
      dellae2   = zero
      ecdo2     = zero
      ecko2     = zero
      qaero     = zero

c     -- set work arrays

      do n = 1, ntc
        it = n + itc - 1
        do k = 1, km
          do i = 1, im
            if (k <= kmax(i)) qaero(i,k,n) = max(qamin, qtr(i,k,it))
          enddo
        enddo
      enddo

      do k = 1, km
        do i = 1, im
          xmbp(i,k) = g * xmb(i) / delp(i,k)
        enddo
      enddo

      do n = 1, ntc
c       -- interface level
        do k = 1, km1
          kp1 = k + 1
          do i = 1, im
            if (kp1 <= kmax(i)) ctro2(i,k,n) =
     &        half * (qaero(i,k,n) + qaero(i,kp1,n))
          enddo
        enddo
c       -- top level
        do i = 1, im
          ctro2(i,kmax(i),n) = qaero(i,kmax(i),n)
        enddo
      enddo

      do n = 1, ntc
        do k = 1, km
          do i = 1, im
            if (cnvflg(i) .and. (k <= kb(i)))
     &        ecko2(i,k,n) = ctro2(i,k,n)
          enddo
        enddo
      enddo

      do n = 1, ntc
        do i = 1, im
          if (cnvflg(i)) ecdo2(i,jmin(i),n) = ctro2(i,jmin(i),n)
        enddo
      enddo

c     do chemical tracers, first need to know how much reevaporates

c     aerosol re-evaporation is set to zero for now
c     uncomment and edit the following code to enable re-evaporation
c     chem_pwd  = zero
c     pwdper    = zero
c     pwav      = zero
c     do i = 1, im
c        do k=1,jmin(i)
c           pwdper(i,k)= -edto(i)*pwdo(i,k)/pwavo(i)
c        enddo
c     enddo
c
c     calculate include mixing ratio (ecko2), how much goes into
c     rainwater to be rained out (chem_pw), and total scavenged,
c     if not reevaporated (pwav)

      do n = 1, ntc
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i)) then
              if ((k > kb(i)) .and. (k < ktcon(i))) then
                dz   = zi(i,k) - zi(i,kk)
                tem  = half    * (xlamue(i,k)+xlamue(i,kk)) * dz
                tem1 = quarter * (xlamud(i,k)+xlamud(i,kk)) * dz
                factor = one + tem - tem1

c               if conserved (not scavenging) then
                ecko2(i,k,n) = ((one-tem1)*ecko2(i,kk,n)
     &            + half*tem*(ctro2(i,k,n)+ctro2(i,kk,n)))/factor

c               how much will be scavenged
c
c               this choice was used in GF, and is also described in a
c               successful implementation into CESM in GRL (Yu et al. 2019),
c               it uses dimesnsionless scavenging coefficients (fscav),
c               but includes henry coeffs with gas phase chemistry

c               fraction fscav is going into liquid
                chem_c(i,k,n)=fscav(n)*ecko2(i,k,n)

c               of that part is going into rain out (chem_pw)
                tem=chem_c(i,k,n)/(one+c0t(i,k)*dz)
                chem_pw(i,k,n)=c0t(i,k)*dz*tem*eta(i,kk) !etah
                ecko2(i,k,n)=tem+ecko2(i,k,n)-chem_c(i,k,n)

c               pwav needed fo reevaporation in downdraft
c               if including reevaporation, please uncomment code below
c               pwav(i,n)=pwav(i,n)+chem_pw(i,k,n)
              endif
            endif
          enddo
        enddo
        do k = 1, km1
          do i = 1, im
            if (k >= ktcon(i)) ecko2(i,k,n)=ctro2(i,k,n)
          enddo
        enddo
      enddo

c     reevaporation of some, pw and pwd terms needed later for dellae2

      do n = 1, ntc
        do k = km1, 1, -1
          kp1 = k + 1
          do i = 1, im
            if (cnvflg(i) .and. (k < jmin(i))) then
              dz = zi(i,kp1) - zi(i,k)
              if (k >= kbcon(i)) then
                tem  = xlamde * dz
                tem1 = half * xlamdd * dz
              else
                tem  = xlamde * dz
                tem1 = half * (xlamd(i)+xlamdd) * dz
              endif
              factor = one + tem - tem1
              ecdo2(i,k,n) = ((one-tem1)*ecdo2(i,kp1,n)
     &          +half*tem*(ctro2(i,k,n)+ctro2(i,kp1,n)))/factor
c             if including reevaporation, please uncomment code below
c             ecdo2(i,k,n)=ecdo2(i,k,n)+pwdper(i,kp1)*pwav(i,n)
c             chem_pwd(i,k,n)=max(zero,pwdper(i,kp1)*pwav(i,n))
            endif
          enddo
        enddo
      enddo

      do n = 1, ntc
        do i = 1, im
          if (cnvflg(i)) then
c           subsidence term treated in fct routine
            dellae2(i,1,n) = edto(i)*etad(i,1)*ecdo2(i,1,n)*xmbp(i,1)
          endif
        enddo
      enddo

      do n = 1, ntc
        do i = 1, im
          if (cnvflg(i)) then
            k = ktcon(i)
            kk = k - 1
c           for the subsidence term already is considered
            dellae2(i,k,n) = eta(i,kk) * ecko2(i,kk,n) * xmbp(i,k)
          endif
        enddo
      enddo

c     --- for updraft & downdraft vertical transport
c
c     initialize maximum allowed timestep for upstream difference approach
c
      dtime_max=delt
      do k=2,km1
        kk = k - 1
        do i = 1, im
          if (kk < ktcon(i)) dtime_max = min(dtime_max,half*delp(i,kk))
        enddo
      enddo

c     now for every chemistry tracer
      do n = 1, ntc
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i) .and. (k < ktcon(i))) then
              dz = zi(i,k) - zi(i,kk)
              aup = one
              if (k <= kb(i)) aup = zero
              adw = one
              if (k > jmin(i)) adw = zero

              dv1q = half * (ecko2(i,k,n) + ecko2(i,kk,n))
              dv2q = half * (ctro2(i,k,n) + ctro2(i,kk,n))
              dv3q = half * (ecdo2(i,k,n) + ecdo2(i,kk,n))

              tem  = half * (xlamue(i,k) + xlamue(i,kk))
              tem1 = half * (xlamud(i,k) + xlamud(i,kk))

              if (k <= kbcon(i)) then
                ptem  = xlamde
                ptem1 = xlamd(i) + xlamdd
              else
                ptem  = xlamde
                ptem1 = xlamdd
              endif
              dellae2(i,k,n) = dellae2(i,k,n) +
c                 detrainment from updraft
     &          ( aup*tem1*eta(i,kk)*dv1q
c                 entrainement into up and downdraft
     &          - (aup*tem*eta(i,kk)+adw*edto(i)*ptem*etad(i,k))*dv2q
c                 detrainment from downdraft
     &          + (adw*edto(i)*ptem1*etad(i,k)*dv3q) ) * dz * xmbp(i,k)

              wet_dep(i,k,n)=chem_pw(i,k,n)*g/delp(i,k)

c             sinks from where updraft and downdraft start
              if (k == jmin(i)+1) then
                dellae2(i,k,n) = dellae2(i,k,n)
     &            -edto(i)*etad(i,kk)*ctro2(i,kk,n)*xmbp(i,k)
              endif
              if (k == kb(i))then
                dellae2(i,k,n) = dellae2(i,k,n)
     &            -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
              endif
            endif
          enddo
        enddo

        do i = 1, im
          if (cnvflg(i)) then
            if (kb(i) == 1) then
              k=kb(i)
              dellae2(i,k,n) = dellae2(i,k,n)
     &          -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
            endif
          endif
        enddo

      enddo

c     for every tracer...

      do n = 1, ntc
        flx_lo  = zero
        totlout = zero
        clipout = zero
c       compute low-order mass flux, upstream
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i) .and. (kk < ktcon(i))) then
              tem = zero
              if (kk >= kb(i)  ) tem = eta(i,kk)
              if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
c             low-order flux,upstream
              if (tem > zero) then
                flx_lo(i,k) = -xmb(i) * tem * qaero(i,k,n)
              elseif (tem < zero) then
                flx_lo(i,k) = -xmb(i) * tem * qaero(i,kk,n)
              endif
            endif
          enddo
        enddo

c       --- make sure low-ord fluxes don't violate positive-definiteness
        do k=1,km1
          kp1 = k + 1
          do i=1,im
            if (cnvflg(i) .and. (k <= ktcon(i))) then
c             time step / grid spacing
              dtovdz = g * dtime_max / abs(delp(i,k))
c             total flux out
              totlout(i,k)=max(zero,flx_lo(i,kp1))-min(zero,flx_lo(i,k))
              clipout(i,k)=min(one ,qaero(i,k,n)/max(epsil,totlout(i,k))
     &           / (1.0001_kind_phys*dtovdz))
            endif
          enddo
        enddo

c       recompute upstream mass fluxes
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i) .and. (kk < ktcon(i))) then
              tem = zero
              if (kk >= kb(i)  ) tem = eta(i,kk)
              if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
              if (tem > zero) then
                flx_lo(i,k) = flx_lo(i,k) * clipout(i,k)
              elseif (tem < zero) then
                flx_lo(i,k) = flx_lo(i,k) * clipout(i,kk)
              endif
            endif
          enddo
        enddo

c      --- a positive-definite low-order (diffusive) solution for the subsidnce fluxes
       do k=1,km1
         kp1 = k + 1
         do i=1,im
           if (cnvflg(i) .and. (k <= ktcon(i))) then
             dtovdz = g * dtime_max / abs(delp(i,k)) ! time step /grid spacing
             dellae2(i,k,n) = dellae2(i,k,n)
     &         -(flx_lo(i,kp1)-flx_lo(i,k))*dtovdz/dtime_max
           endif
          enddo
        enddo

      enddo ! ctr

c     convert wet deposition to total mass deposited over dt and dp

      do n = 1, ntc
        do k = 1, km
          do i = 1, im
            if (cnvflg(i) .and. (k < ktcon(i)))
     &        wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*delt*delp(i,k)
          enddo
        enddo
      enddo

c     compute final aerosol concentrations

      do n = 1, ntc
        do k = 1, km
          do i = 1, im
            if (cnvflg(i) .and. (k <= min(kmax(i),ktcon(i)))) then
              qaero(i,k,n) = qaero(i,k,n) + dellae2(i,k,n) * delt
              if (qaero(i,k,n) < zero) then
c               add negative mass to wet deposition
                wet_dep(i,k,n) = wet_dep(i,k,n)-qaero(i,k,n)*delp(i,k)
                qaero(i,k,n) = qamin
              endif
            endif
          enddo
        enddo
      enddo

      return
      end


      subroutine samfshalcnv_aerosols(im, ix, km, itc, ntc, ntr, delt,
     &  cnvflg, kb, kmax, kbcon, ktcon, fscav,
     &  xmb, c0t, eta, zi, xlamue, xlamud, delp,
     &  qtr, qaero)

      use machine , only : kind_phys
      use physcons, only : g => con_g, qamin

      implicit none

c     -- input arguments
      integer,                                      intent(in) :: im,
     &  ix, km, itc, ntc, ntr
      real(kind=kind_phys),                         intent(in) :: delt
!     &  xlamde, xlamdd
      logical,              dimension(im),          intent(in) :: cnvflg
!     integer,              dimension(im),          intent(in) :: jmin,
      integer,              dimension(im),          intent(in) ::
     &  kb, kmax, kbcon, ktcon
      real(kind=kind_phys), dimension(im),          intent(in) ::
     &  xmb, xlamud
      real(kind=kind_phys), dimension(ntc),         intent(in) :: fscav
      real(kind=kind_phys), dimension(im,km),       intent(in) :: c0t,
     &  eta, zi, xlamue  !, xlamud
      real(kind=kind_phys), dimension(ix,km),       intent(in) :: delp
      real(kind=kind_phys), dimension(ix,km,ntr+2), intent(in) :: qtr
c     -- output arguments
      real(kind=kind_phys), dimension(im,km,ntc),   intent(out) :: qaero

c     -- local variables
c     -- general variables
      integer :: i, indx, it, k, kk, km1, kp1, n
!     real(kind=kind_phys) :: adw, aup, dtime_max, dv1q, dv2q, dv3q,
      real(kind=kind_phys) :: aup, dtime_max, dv1q, dv2q, dv3q,
     &  dtovdz, dz, factor, ptem, ptem1, qamax, tem, tem1
      real(kind=kind_phys), dimension(ix,km) :: xmbp
c     -- chemical transport variables
      real(kind=kind_phys), dimension(im,km,ntc) :: ctro2,ecko2,dellae2
c     -- additional variables for tracers for wet deposition,
      real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw,
     &  wet_dep
c     -- if reevaporation is enabled, uncomment lines below
c     real(kind=kind_phys), dimension(im,ntc) :: pwav
c     real(kind=kind_phys), dimension(im,km) :: pwdper
c     real(kind=kind_phys), dimension(im,km,ntr) :: chem_pwd
c     -- additional variables for fct
      real(kind=kind_phys), dimension(im,km) :: flx_lo, totlout, clipout

      real(kind=kind_phys), parameter :: one     = 1.0_kind_phys
      real(kind=kind_phys), parameter :: half    = 0.5_kind_phys
      real(kind=kind_phys), parameter :: quarter = 0.25_kind_phys
      real(kind=kind_phys), parameter :: zero    = 0.0_kind_phys
      real(kind=kind_phys), parameter :: epsil   = 1.e-22_kind_phys    ! prevent division by zero
      real(kind=kind_phys), parameter :: escav   = 0.8_kind_phys       ! wet scavenging efficiency

c     -- begin

c     -- check if aerosols are present
      if ( ntc <= 0 .or. itc <= 0 .or. ntr <= 0 ) return
      if ( ntr  < itc + ntc - 3 ) return

c     -- initialize work variables
      km1 = km - 1

      chem_c    = zero
      chem_pw   = zero
      ctro2     = zero
      dellae2   = zero
      !ecdo2     = zero
      ecko2     = zero
      qaero     = zero

c     -- set work arrays

      do n = 1, ntc
        it = n + itc - 1
        do k = 1, km
          do i = 1, im
            if (k <= kmax(i)) qaero(i,k,n) = max(qamin, qtr(i,k,it))
          enddo
        enddo
      enddo

      do k = 1, km
        do i = 1, im
          xmbp(i,k) = g * xmb(i) / delp(i,k)
        enddo
      enddo

      do n = 1, ntc
c       -- interface level
        do k = 1, km1
          kp1 = k + 1
          do i = 1, im
            if (kp1 <= kmax(i)) ctro2(i,k,n) =
     &        half * (qaero(i,k,n) + qaero(i,kp1,n))
          enddo
        enddo
c       -- top level
        do i = 1, im
          ctro2(i,kmax(i),n) = qaero(i,kmax(i),n)
        enddo
      enddo

      do n = 1, ntc
        do k = 1, km
          do i = 1, im
            if (cnvflg(i) .and. (k <= kb(i)))
     &        ecko2(i,k,n) = ctro2(i,k,n)
          enddo
        enddo
      enddo

      !do n = 1, ntc
      !  do i = 1, im
      !    if (cnvflg(i)) ecdo2(i,jmin(i),n) = ctro2(i,jmin(i),n)
      !  enddo
      !enddo

c     do chemical tracers, first need to know how much reevaporates

c     aerosol re-evaporation is set to zero for now
c     uncomment and edit the following code to enable re-evaporation
c     chem_pwd  = zero
c     pwdper    = zero
c     pwav      = zero
c     do i = 1, im
c        do k=1,jmin(i)
c           pwdper(i,k)= -edto(i)*pwdo(i,k)/pwavo(i)
c        enddo
c     enddo
c
c     calculate include mixing ratio (ecko2), how much goes into
c     rainwater to be rained out (chem_pw), and total scavenged,
c     if not reevaporated (pwav)

      do n = 1, ntc
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i)) then
              if ((k > kb(i)) .and. (k < ktcon(i))) then
                dz   = zi(i,k) - zi(i,kk)
                tem  = half    * (xlamue(i,k)+xlamue(i,kk)) * dz
!               tem1 = quarter * (xlamud(i,k)+xlamud(i,kk)) * dz
                tem1 = quarter * (xlamud(i  )+xlamud(i   )) * dz
                factor = one + tem - tem1

c               if conserved (not scavenging) then
                ecko2(i,k,n) = ((one-tem1)*ecko2(i,kk,n)
     &            + half*tem*(ctro2(i,k,n)+ctro2(i,kk,n)))/factor

c               how much will be scavenged
c
c               this choice was used in GF, and is also described in a
c               successful implementation into CESM in GRL (Yu et al. 2019),
c               it uses dimesnsionless scavenging coefficients (fscav),
c               but includes henry coeffs with gas phase chemistry

c               fraction fscav is going into liquid
                chem_c(i,k,n)=escav*fscav(n)*ecko2(i,k,n)

c               of that part is going into rain out (chem_pw)
                tem=chem_c(i,k,n)/(one+c0t(i,k)*dz)
                chem_pw(i,k,n)=c0t(i,k)*dz*tem*eta(i,kk) !etah
                ecko2(i,k,n)=tem+ecko2(i,k,n)-chem_c(i,k,n)

c               pwav needed fo reevaporation in downdraft
c               if including reevaporation, please uncomment code below
c               pwav(i,n)=pwav(i,n)+chem_pw(i,k,n)
              endif
            endif
          enddo
        enddo
        do k = 1, km1
          do i = 1, im
            if (k >= ktcon(i)) ecko2(i,k,n)=ctro2(i,k,n)
          enddo
        enddo
      enddo

c     reevaporation of some, pw and pwd terms needed later for dellae2

!      do n = 1, ntc
!        do k = km1, 1, -1
!          kp1 = k + 1
!          do i = 1, im
!            if (cnvflg(i) .and. (k < jmin(i))) then
!              dz = zi(i,kp1) - zi(i,k)
!              if (k >= kbcon(i)) then
!                tem  = xlamde * dz
!                tem1 = half * xlamdd * dz
!              else
!                tem  = xlamde * dz
!                tem1 = half * (xlamd(i)+xlamdd) * dz
!              endif
!              factor = one + tem - tem1
!              ecdo2(i,k,n) = ((one-tem1)*ecdo2(i,kp1,n)
!     &          +half*tem*(ctro2(i,k,n)+ctro2(i,kp1,n)))/factor
c             if including reevaporation, please uncomment code below
c             ecdo2(i,k,n)=ecdo2(i,k,n)+pwdper(i,kp1)*pwav(i,n)
c             chem_pwd(i,k,n)=max(zero,pwdper(i,kp1)*pwav(i,n))
!            endif
!          enddo
!        enddo
!      enddo

!      do n = 1, ntc
!        do i = 1, im
!          if (cnvflg(i)) then
c           subsidence term treated in fct routine
!            dellae2(i,1,n) = edto(i)*etad(i,1)*ecdo2(i,1,n)*xmbp(i,1)
!          endif
!        enddo
!      enddo

      do n = 1, ntc
        do i = 1, im
          if (cnvflg(i)) then
            k = ktcon(i)
            kk = k - 1
c           for the subsidence term already is considered
            dellae2(i,k,n) = eta(i,kk) * ecko2(i,kk,n) * xmbp(i,k)
          endif
        enddo
      enddo

c     --- for updraft & downdraft vertical transport
c
c     initialize maximum allowed timestep for upstream difference approach
c
      dtime_max=delt
      do k=2,km1
        kk = k - 1
        do i = 1, im
          if (kk < ktcon(i)) dtime_max = min(dtime_max,half*delp(i,kk))
        enddo
      enddo

c     now for every chemistry tracer
      do n = 1, ntc
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i) .and. (k < ktcon(i))) then
              dz = zi(i,k) - zi(i,kk)
              aup = one
              if (k <= kb(i)) aup = zero
!              adw = one
!              if (k > jmin(i)) adw = zero

              dv1q = half * (ecko2(i,k,n) + ecko2(i,kk,n))
              dv2q = half * (ctro2(i,k,n) + ctro2(i,kk,n))
c             dv3q = half * (ecdo2(i,k,n) + ecdo2(i,kk,n))

              tem  = half * (xlamue(i,k) + xlamue(i,kk))
             !tem1 = half * (xlamud(i,k) + xlamud(i,kk))
              tem1 = half * (xlamud(i  ) + xlamud(i   ))

!              if (k <= kbcon(i)) then
!                ptem  = xlamde
!                ptem1 = xlamd(i) + xlamdd
!              else
!                ptem  = xlamde
!                ptem1 = xlamdd
!              endif
              dellae2(i,k,n) = dellae2(i,k,n) +
c                 detrainment from updraft
     &          ( aup*tem1*eta(i,kk)*dv1q
c                 entrainement into up and downdraft
!    &          - (aup*tem*eta(i,kk)+adw*edto(i)*ptem*etad(i,k))*dv2q
     &          - (aup*tem*eta(i,kk))*dv2q
c                 detrainment from downdraft
!    &          + (adw*edto(i)*ptem1*etad(i,k)*dv3q)
     &            ) * dz * xmbp(i,k)

              wet_dep(i,k,n)=chem_pw(i,k,n)*g/delp(i,k)

c             sinks from where updraft and downdraft start
!              if (k == jmin(i)+1) then
!                dellae2(i,k,n) = dellae2(i,k,n)
!     &            -edto(i)*etad(i,kk)*ctro2(i,kk,n)*xmbp(i,k)
!              endif
              if (k == kb(i))then
                dellae2(i,k,n) = dellae2(i,k,n)
     &            -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
              endif
            endif
          enddo
        enddo

        do i = 1, im
          if (cnvflg(i)) then
            if (kb(i) == 1) then
              k=kb(i)
              dellae2(i,k,n) = dellae2(i,k,n)
     &          -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
            endif
          endif
        enddo

      enddo

c     for every tracer...

      do n = 1, ntc
        flx_lo  = zero
        totlout = zero
        clipout = zero
c       compute low-order mass flux, upstream
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i) .and. (kk < ktcon(i))) then
              tem = zero
              if (kk >= kb(i)  ) tem = eta(i,kk)
!              if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
c             low-order flux,upstream
              if (tem > zero) then
                flx_lo(i,k) = -xmb(i) * tem * qaero(i,k,n)
              elseif (tem < zero) then
                flx_lo(i,k) = -xmb(i) * tem * qaero(i,kk,n)
              endif
            endif
          enddo
        enddo

c       --- make sure low-ord fluxes don't violate positive-definiteness
        do k=1,km1
          kp1 = k + 1
          do i=1,im
            if (cnvflg(i) .and. (k <= ktcon(i))) then
c             time step / grid spacing
              dtovdz = g * dtime_max / abs(delp(i,k))
c             total flux out
              totlout(i,k)=max(zero,flx_lo(i,kp1))-min(zero,flx_lo(i,k))
              clipout(i,k)=min(one ,qaero(i,k,n)/max(epsil,totlout(i,k))
     &           / (1.0001_kind_phys*dtovdz))
            endif
          enddo
        enddo

c       recompute upstream mass fluxes
        do k = 2, km1
          kk = k - 1
          do i = 1, im
            if (cnvflg(i) .and. (kk < ktcon(i))) then
              tem = zero
              if (kk >= kb(i)  ) tem = eta(i,kk)
!             if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
              if (tem > zero) then
                flx_lo(i,k) = flx_lo(i,k) * clipout(i,k)
              elseif (tem < zero) then
                flx_lo(i,k) = flx_lo(i,k) * clipout(i,kk)
              endif
            endif
          enddo
        enddo

c      --- a positive-definite low-order (diffusive) solution for the subsidnce fluxes
       do k=1,km1
         kp1 = k + 1
         do i=1,im
           if (cnvflg(i) .and. (k <= ktcon(i))) then
             dtovdz = g * dtime_max / abs(delp(i,k)) ! time step /grid spacing
             dellae2(i,k,n) = dellae2(i,k,n)
     &         -(flx_lo(i,kp1)-flx_lo(i,k))*dtovdz/dtime_max
           endif
          enddo
        enddo

      enddo ! ctr

c     convert wet deposition to total mass deposited over dt and dp

      do n = 1, ntc
        do k = 1, km
          do i = 1, im
            if (cnvflg(i) .and. (k < ktcon(i)))
     &        wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*delt*delp(i,k)
          enddo
        enddo
      enddo

c     compute final aerosol concentrations

      do n = 1, ntc
        do k = 1, km
          do i = 1, im
            if (cnvflg(i) .and. (k <= min(kmax(i),ktcon(i)))) then
              qaero(i,k,n) = qaero(i,k,n) + dellae2(i,k,n) * delt
              if (qaero(i,k,n) < zero) then
c               add negative mass to wet deposition
                wet_dep(i,k,n) = wet_dep(i,k,n)-qaero(i,k,n)*delp(i,k)
                qaero(i,k,n) = qamin
              endif
            endif
          enddo
        enddo
      enddo

      return
      end