module wetdep 3,5

!----------------------------------------------------------------------- 
!
! Wet deposition routines for both aerosols and gas phase constituents.
! 
!-----------------------------------------------------------------------

use shr_kind_mod, only: r8 => shr_kind_r8
use ppgrid,       only: pcols, pver
use physconst,    only: gravit, rair, tmelt
use phys_control, only: cam_physpkg_is
use cam_logfile,  only: iulog

implicit none
save
private

public :: wetdepa_v1  ! scavenging codes for very soluble aerosols -- CAM4 version
public :: wetdepa_v2  ! scavenging codes for very soluble aerosols -- CAM5 version
public :: wetdepg     ! scavenging of gas phase constituents by henry's law
public :: clddiag     ! calc of cloudy volume and rain mixing ratio

real(r8), parameter :: cmftau = 3600._r8
real(r8), parameter :: rhoh2o = 1000._r8            ! density of water
real(r8), parameter :: molwta = 28.97_r8            ! molecular weight dry air gm/mole

!==============================================================================
contains
!==============================================================================


subroutine clddiag(t, pmid, pdel, cmfdqr, evapc, & 1
                   cldt, cldcu, cldst, cme, evapr, &
                   prain, cldv, cldvcu, cldvst, rain, &
                   ncol)

   ! ------------------------------------------------------------------------------------ 
   ! Estimate the cloudy volume which is occupied by rain or cloud water as
   ! the max between the local cloud amount or the
   ! sum above of (cloud*positive precip production)      sum total precip from above
   !              ----------------------------------   x ------------------------
   ! sum above of     (positive precip           )        sum positive precip from above
   ! Author: P. Rasch
   !         Sungsu Park. Mar.2010 
   ! ------------------------------------------------------------------------------------

   ! Input arguments:
   real(r8), intent(in) :: t(pcols,pver)        ! temperature (K)
   real(r8), intent(in) :: pmid(pcols,pver)     ! pressure at layer midpoints
   real(r8), intent(in) :: pdel(pcols,pver)     ! pressure difference across layers
   real(r8), intent(in) :: cmfdqr(pcols,pver)   ! dq/dt due to convective rainout 
   real(r8), intent(in) :: evapc(pcols,pver)    ! Evaporation rate of convective precipitation ( >= 0 ) 
   real(r8), intent(in) :: cldt(pcols,pver)    ! total cloud fraction
   real(r8), intent(in) :: cldcu(pcols,pver)    ! Cumulus cloud fraction
   real(r8), intent(in) :: cldst(pcols,pver)    ! Stratus cloud fraction
   real(r8), intent(in) :: cme(pcols,pver)      ! rate of cond-evap within the cloud
   real(r8), intent(in) :: evapr(pcols,pver)    ! rate of evaporation of falling precipitation (kg/kg/s)
   real(r8), intent(in) :: prain(pcols,pver)    ! rate of conversion of condensate to precipitation (kg/kg/s)
   integer, intent(in) :: ncol

   ! Output arguments:
   real(r8), intent(out) :: cldv(pcols,pver)     ! fraction occupied by rain or cloud water 
   real(r8), intent(out) :: cldvcu(pcols,pver)   ! Convective precipitation volume
   real(r8), intent(out) :: cldvst(pcols,pver)   ! Stratiform precipitation volume
   real(r8), intent(out) :: rain(pcols,pver)     ! mixing ratio of rain (kg/kg)

   ! Local variables:
   integer  i, k
   real(r8) convfw         ! used in fallspeed calculation; taken from findmcnew
   real(r8) sumppr(pcols)        ! precipitation rate (kg/m2-s)
   real(r8) sumpppr(pcols)       ! sum of positive precips from above
   real(r8) cldv1(pcols)         ! precip weighted cloud fraction from above
   real(r8) lprec                ! local production rate of precip (kg/m2/s)
   real(r8) lprecp               ! local production rate of precip (kg/m2/s) if positive
   real(r8) rho                  ! air density
   real(r8) vfall
   real(r8) sumppr_cu(pcols)     ! Convective precipitation rate (kg/m2-s)
   real(r8) sumpppr_cu(pcols)    ! Sum of positive convective precips from above
   real(r8) cldv1_cu(pcols)      ! Convective precip weighted convective cloud fraction from above
   real(r8) lprec_cu             ! Local production rate of convective precip (kg/m2/s)
   real(r8) lprecp_cu            ! Local production rate of convective precip (kg/m2/s) if positive
   real(r8) sumppr_st(pcols)     ! Stratiform precipitation rate (kg/m2-s)
   real(r8) sumpppr_st(pcols)    ! Sum of positive stratiform precips from above
   real(r8) cldv1_st(pcols)      ! Stratiform precip weighted stratiform cloud fraction from above
   real(r8) lprec_st             ! Local production rate of stratiform precip (kg/m2/s)
   real(r8) lprecp_st            ! Local production rate of stratiform precip (kg/m2/s) if positive
   ! -----------------------------------------------------------------------

   convfw = 1.94_r8*2.13_r8*sqrt(rhoh2o*gravit*2.7e-4_r8)
   do i=1,ncol
      sumppr(i) = 0._r8
      cldv1(i) = 0._r8
      sumpppr(i) = 1.e-36_r8
      sumppr_cu(i)  = 0._r8
      cldv1_cu(i)   = 0._r8
      sumpppr_cu(i) = 1.e-36_r8
      sumppr_st(i)  = 0._r8
      cldv1_st(i)   = 0._r8
      sumpppr_st(i) = 1.e-36_r8
   end do

   do k = 1,pver
      do i = 1,ncol
         cldv(i,k) = &
            max(min(1._r8, &
            cldv1(i)/sumpppr(i) &
            )*sumppr(i)/sumpppr(i), &
            cldt(i,k) &
            )
         lprec = pdel(i,k)/gravit &
            *(prain(i,k)+cmfdqr(i,k)-evapr(i,k))
         lprecp = max(lprec,1.e-30_r8)
         cldv1(i) = cldv1(i)  + cldt(i,k)*lprecp
         sumppr(i) = sumppr(i) + lprec
         sumpppr(i) = sumpppr(i) + lprecp

         ! For convective precipitation volume at the top interface of each layer. Neglect the current layer.
         cldvcu(i,k)   = max(min(1._r8,cldv1_cu(i)/sumpppr_cu(i))*(sumppr_cu(i)/sumpppr_cu(i)),0._r8)
         lprec_cu      = (pdel(i,k)/gravit)*(cmfdqr(i,k)-evapc(i,k))
         lprecp_cu     = max(lprec_cu,1.e-30_r8)
         cldv1_cu(i)   = cldv1_cu(i) + cldcu(i,k)*lprecp_cu
         sumppr_cu(i)  = sumppr_cu(i) + lprec_cu
         sumpppr_cu(i) = sumpppr_cu(i) + lprecp_cu

         ! For stratiform precipitation volume at the top interface of each layer. Neglect the current layer.
         cldvst(i,k)   = max(min(1._r8,cldv1_st(i)/sumpppr_st(i))*(sumppr_st(i)/sumpppr_st(i)),0._r8)
         lprec_st      = (pdel(i,k)/gravit)*(prain(i,k)-evapr(i,k))
         lprecp_st     = max(lprec_st,1.e-30_r8)
         cldv1_st(i)   = cldv1_st(i) + cldst(i,k)*lprecp_st
         sumppr_st(i)  = sumppr_st(i) + lprec_st
         sumpppr_st(i) = sumpppr_st(i) + lprecp_st

         rain(i,k) = 0._r8
         if(t(i,k) .gt. tmelt) then
            rho = pmid(i,k)/(rair*t(i,k))
            vfall = convfw/sqrt(rho)
            rain(i,k) = sumppr(i)/(rho*vfall)
            if (rain(i,k).lt.1.e-14_r8) rain(i,k) = 0._r8
         endif
      end do
   end do

end subroutine clddiag

!==============================================================================

! This is the CAM5 version of wetdepa.


subroutine wetdepa_v2(t, p, q, pdel, &
                   cldt, cldc, cmfdqr, evapc, conicw, precs, conds, &
                       evaps, cwat, tracer, deltat, &
                       scavt, iscavt, cldv, cldvcu, cldvst, dlf, fracis, sol_fact, ncol, &
                       scavcoef,  &
#ifdef MODAL_AERO
                       is_strat_cloudborne, rate1ord_cw2pr_st, qqcw, f_act_conv, &     ! rce 2010/05/01
#endif
                       icscavt, isscavt, bcscavt, bsscavt, &
                       sol_facti_in, sol_factbi_in, sol_factii_in, &
                       sol_factic_in, sol_factiic_in )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! scavenging code for very soluble aerosols
      ! 
      ! Author: P. Rasch
      ! Modified by T. Bond 3/2003 to track different removals
      ! Sungsu Park. Mar.2010 : Impose consistencies with a few changes in physics.
      !-----------------------------------------------------------------------



      implicit none

      real(r8), intent(in) ::&
         t(pcols,pver),        &! temperature
         p(pcols,pver),        &! pressure
         q(pcols,pver),        &! moisture
         pdel(pcols,pver),     &! pressure thikness
         cldt(pcols,pver),    &! total cloud fraction
         cldc(pcols,pver),     &! convective cloud fraction
         cmfdqr(pcols,pver),   &! rate of production of convective precip
! Sungsu
         evapc(pcols,pver),    &! Evaporation rate of convective precipitation
! Sungsu
         conicw(pcols,pver),   &! convective cloud water
         cwat(pcols,pver),     &! cloud water amount 
         precs(pcols,pver),    &! rate of production of stratiform precip
         conds(pcols,pver),    &! rate of production of condensate
         evaps(pcols,pver),    &! rate of evaporation of precip
         cldv(pcols,pver),     &! total cloud fraction
! Sungsu
         cldvcu(pcols,pver),   &! Convective precipitation area at the top interface of each layer
         cldvst(pcols,pver),   &! Stratiform precipitation area at the top interface of each layer
         dlf(pcols,pver),      &! Detrainment of convective condensate [kg/kg/s]
! Sungsu
         deltat,               &! time step
         tracer(pcols,pver)     ! trace species
      ! If subroutine is called with just sol_fact:
            ! sol_fact is used for both in- and below-cloud scavenging
      ! If subroutine is called with optional argument sol_facti_in:
            ! sol_fact  is used for below cloud scavenging
            ! sol_facti is used for in cloud scavenging
         real(r8), intent(in) :: sol_fact ! solubility factor (fraction of aer scavenged below & in, or just below or sol_facti_in is provided)
         integer, intent(in) :: ncol
         real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
#ifdef MODAL_AERO
         ! rce 2010/05/01
         ! is_strat_cloudborne = .true. if tracer is stratiform-cloudborne aerosol; else .false. 
         logical, intent(in), optional :: is_strat_cloudborne   
         ! rate1ord_cw2pr_st = 1st order rate for strat cw to precip (1/s) 
         real(r8), intent(in), optional  :: rate1ord_cw2pr_st(pcols,pver)
         ! qqcw = strat-cloudborne aerosol corresponding to tracer when is_strat_cloudborne==.false.; else 0.0 
         real(r8), intent(in), optional  :: qqcw(pcols,pver)
         ! f_act_conv = conv-cloud activation fraction when is_strat_cloudborne==.false.; else 0.0 
         real(r8), intent(in), optional  :: f_act_conv(pcols,pver)
         ! end rce 2010/05/01
#endif

         real(r8), intent(in), optional :: sol_facti_in   ! solubility factor (frac of aerosol scavenged in cloud)
         real(r8), intent(in), optional :: sol_factbi_in  ! solubility factor (frac of aerosol scavenged below cloud by ice)
         real(r8), intent(in), optional :: sol_factii_in  ! solubility factor (frac of aerosol scavenged in cloud by ice)
         real(r8), intent(in), optional :: sol_factic_in(pcols,pver)  ! sol_facti_in for convective clouds
         real(r8), intent(in), optional :: sol_factiic_in ! sol_factii_in for convective clouds
         
      real(r8), intent(out) ::&
         scavt(pcols,pver),    &! scavenging tend 
         iscavt(pcols,pver),   &! incloud scavenging tends
         fracis(pcols,pver)     ! fraction of species not scavenged

      real(r8), intent(out), optional ::    icscavt(pcols,pver)     ! incloud, convective
      real(r8), intent(out), optional ::    isscavt(pcols,pver)     ! incloud, stratiform
      real(r8), intent(out), optional ::    bcscavt(pcols,pver)     ! below cloud, convective
      real(r8), intent(out), optional ::    bsscavt(pcols,pver)     ! below cloud, stratiform

      ! local variables

      integer i                 ! x index
      integer k                 ! z index

      real(r8) adjfac               ! factor stolen from cmfmca
      real(r8) aqfrac               ! fraction of tracer in aqueous phase
      real(r8) cwatc                ! local convective total water amount 
      real(r8) cwats                ! local stratiform total water amount 
      real(r8) cwatp                ! local water amount falling from above precip
      real(r8) fracev(pcols)        ! fraction of precip from above that is evaporating
! Sungsu
      real(r8) fracev_cu(pcols)     ! Fraction of convective precip from above that is evaporating
! Sungsu
      real(r8) fracp                ! fraction of cloud water converted to precip
      real(r8) gafrac               ! fraction of tracer in gas phasea
      real(r8) hconst               ! henry's law solubility constant when equation is expressed
                                ! in terms of mixing ratios
      real(r8) mpla                 ! moles / liter H2O entering the layer from above
      real(r8) mplb                 ! moles / liter H2O leaving the layer below
      real(r8) omsm                 ! 1 - (a small number)
      real(r8) part                 !  partial pressure of tracer in atmospheres
      real(r8) patm                 ! total pressure in atmospheres
      real(r8) pdog                 ! work variable (pdel/gravit)
      real(r8) precabc(pcols)       ! conv precip from above (work array)
      real(r8) precabs(pcols)       ! strat precip from above (work array)
      real(r8) precbl               ! precip falling out of level (work array)
      real(r8) precmin              ! minimum convective precip causing scavenging
      real(r8) rat(pcols)           ! ratio of amount available to amount removed
      real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
      real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
      real(r8) srcc                 ! tend for convective rain
      real(r8) srcs                 ! tend for stratiform rain
      real(r8) srct(pcols)          ! work variable
      real(r8) tracab(pcols)        ! column integrated tracer amount
!      real(r8) vfall                ! fall speed of precip
      real(r8) fins                 ! fraction of rem. rate by strat rain
      real(r8) finc                 ! fraction of rem. rate by conv. rain
      real(r8) srcs1                ! work variable
      real(r8) srcs2                ! work variable
      real(r8) tc                   ! temp in celcius
      real(r8) weight               ! fraction of condensate which is ice
      real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
      real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
      real(r8) odds                 ! limit on removal rate (proportional to prec)
      real(r8) dblchek(pcols)
      logical :: found

    ! Jan.16.2009. Sungsu for wet scavenging below clouds.
    ! real(r8) cldovr_cu(pcols)     ! Convective precipitation area at the base of each layer
    ! real(r8) cldovr_st(pcols)     ! Stratiform precipitation area at the base of each layer

      real(r8) tracer_incu
      real(r8) tracer_mean

    ! End by Sungsu

      real(r8) sol_facti,  sol_factb  ! in cloud and below cloud fraction of aerosol scavenged
      real(r8) sol_factii, sol_factbi ! in cloud and below cloud fraction of aerosol scavenged by ice
      real(r8) sol_factic(pcols,pver)             ! sol_facti for convective clouds
      real(r8) sol_factiic            ! sol_factii for convective clouds
      ! sol_factic & solfact_iic added for MODAL_AERO.  
      ! For stratiform cloud, cloudborne aerosol is treated explicitly,
      !    and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
      ! For convective cloud, cloudborne aerosol is not treated explicitly,
      !    and sol_factic is 1.0 for both cloudborne and interstitial.

      ! ------------------------------------------------------------------------
!      omsm = 1.-1.e-10          ! used to prevent roundoff errors below zero
      omsm = 1._r8-2*epsilon(1._r8) ! used to prevent roundoff errors below zero
      precmin =  0.1_r8/8.64e4_r8      ! set critical value to 0.1 mm/day in kg/m2/s

      adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme

      ! assume 4 m/s fall speed currently (should be improved)
!      vfall = 4.
	
      ! default (if other sol_facts aren't in call, set all to required sol_fact
      sol_facti = sol_fact
      sol_factb = sol_fact
      sol_factii = sol_fact
      sol_factbi = sol_fact

      if ( present(sol_facti_in) )  sol_facti = sol_facti_in
      if ( present(sol_factii_in) )  sol_factii = sol_factii_in
      if ( present(sol_factbi_in) )  sol_factbi = sol_factbi_in

      sol_factic  = sol_facti
      sol_factiic = sol_factii
      if ( present(sol_factic_in ) )  sol_factic  = sol_factic_in
      if ( present(sol_factiic_in) )  sol_factiic = sol_factiic_in

      ! this section of code is for highly soluble aerosols,
      ! the assumption is that within the cloud that
      ! all the tracer is in the cloud water
      !
      ! for both convective and stratiform clouds, 
      ! the fraction of cloud water converted to precip defines
      ! the amount of tracer which is pulled out.
      !

      do i = 1,pcols
         precabs(i) = 0
         precabc(i) = 0
         scavab(i) = 0
         scavabc(i) = 0
         tracab(i) = 0
         cldmabs(i) = 0
         cldmabc(i) = 0

       ! Jan.16. Sungsu 
       ! I added below to compute vertically projected cumulus and stratus fractions from the top to the
       ! current model layer by assuming a simple independent maximum overlapping assumption for 
       ! each cloud.
       ! cldovr_cu(i) = 0._r8
       ! cldovr_st(i) = 0._r8
       ! End by Sungsu

      end do

      do k = 1,pver
         do i = 1,ncol
            tc     = t(i,k) - tmelt
            weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
            weight = 0._r8                                 ! assume no ice

            pdog = pdel(i,k)/gravit

            ! ****************** Evaporation **************************
            ! calculate the fraction of strat precip from above 
            !                 which evaporates within this layer
            fracev(i) = evaps(i,k)*pdel(i,k)/gravit &
                     /max(1.e-12_r8,precabs(i))

            ! trap to ensure reasonable ratio bounds
            fracev(i) = max(0._r8,min(1._r8,fracev(i)))

! Sungsu : Same as above but convective precipitation part
            fracev_cu(i) = evapc(i,k)*pdel(i,k)/gravit/max(1.e-12_r8,precabc(i))
            fracev_cu(i) = max(0._r8,min(1._r8,fracev_cu(i)))
! Sungsu
            ! ****************** Convection ***************************
            ! now do the convective scavenging

            ! set odds proportional to fraction of the grid box that is swept by the 
            ! precipitation =precabc/rhoh20*(area of sphere projected on plane
            !                                /volume of sphere)*deltat
            ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
            ! unless the fraction of the area that is cloud is less than odds, in which
            ! case use the cloud fraction (assumes precabs is in kg/m2/s)
            ! is really: precabs*3/4/1000./1e-3*deltat
            ! here I use .1 from Balkanski
            !
            ! use a local rate of convective rain production for incloud scav
            !odds=max(min(1._r8, &
            !     cmfdqr(i,k)*pdel(i,k)/gravit*0.1_r8*deltat),0._r8)
            !++mcb -- change cldc to cldt; change cldt to cldv (9/17/96)
            !            srcs1 =  cldt(i,k)*odds*tracer(i,k)*(1.-weight) &
            ! srcs1 =  cldv(i,k)*odds*tracer(i,k)*(1.-weight) &
            !srcs1 =  cldc(i,k)*odds*tracer(i,k)*(1.-weight) &
            !         /deltat 

            ! fraction of convective cloud water converted to rain
            ! Dec.29.2009 : Sungsu multiplied cldc(i,k) to conicw(i,k) below
            ! fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,conicw(i,k))
            ! fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,cldc(i,k)*conicw(i,k))
            ! Sungsu: Below new formula of 'fracp' is necessary since 'conicw' is a LWC/IWC
            !         that has already precipitated out, that is, 'conicw' does not contain
            !         precipitation at all ! 
              fracp = cmfdqr(i,k)*deltat/max(1.e-12_r8,cldc(i,k)*conicw(i,k)+(cmfdqr(i,k)+dlf(i,k))*deltat) ! Sungsu.Mar.19.2010.
            ! Dec.29.2009
            ! Note cmfdrq can be negative from evap of rain, so constrain it <-- This is wrong. cmfdqr does not
            ! contain evaporation of precipitation.
            fracp = max(min(1._r8,fracp),0._r8)
            ! remove that amount from within the convective area
!           srcs1 = cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat ! liquid only
!           srcs1 = cldc(i,k)*fracp*tracer(i,k)/deltat             ! any condensation
!           srcs1 = 0.
!           Jan.02.2010. Sungsu : cldt --> cldc below.
#ifdef MODAL_AERO
            ! rce 2010/05/01
            if ( is_strat_cloudborne ) then
               ! only strat in-cloud removal affects strat-cloudborne aerosol
               srcs1 = 0._r8
            else
               tracer_incu = f_act_conv(i,k)*(tracer(i,k)+& 
                             min(qqcw(i,k),tracer(i,k)*((cldt(i,k)-cldc(i,k))/max(0.01_r8,(1._r8-(cldt(i,k)-cldc(i,k)))))))              
               srcs1 = sol_factic(i,k)*cldc(i,k)*fracp*tracer_incu*(1._r8-weight)/deltat &  ! Liquid
                     + sol_factiic    *cldc(i,k)*fracp*tracer_incu*(weight)/deltat          ! Ice
            end if
#else
            srcs1 = sol_factic(i,k)*cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat &  ! liquid
                 +  sol_factiic*cldc(i,k)*fracp*tracer(i,k)*(weight)/deltat      ! ice
#endif


            !--mcb

            ! scavenge below cloud

            !            cldmabc(i) = max(cldc(i,k),cldmabc(i))
            !            cldmabc(i) = max(cldt(i,k),cldmabc(i))
            ! cldmabc(i) = max(cldv(i,k),cldmabc(i))
            ! cldmabc(i) = cldv(i,k)
            cldmabc(i) = cldvcu(i,k)

            ! Jan. 16. 2010. Sungsu
            ! cldmabc(i) = cldmabc(i) * cldovr_cu(i) / max( 0.01_r8, cldovr_cu(i) + cldovr_st(i) )
            ! End by Sungsu

#ifdef MODAL_AERO
            ! rce 2010/05/01
            if ( is_strat_cloudborne ) then
               ! only strat in-cloud removal affects strat-cloudborne aerosol
               srcs2 = 0._r8
            else
               tracer_mean = tracer(i,k)*(1._r8-cldc(i,k)*f_act_conv(i,k))-cldc(i,k)*f_act_conv(i,k)*&
                             min(qqcw(i,k),tracer(i,k)*((cldt(i,k)-cldc(i,k))/max(0.01_r8,(1._r8-(cldt(i,k)-cldc(i,k))))))
               tracer_mean = max(0._r8,tracer_mean) 
               odds  = max(min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8)*scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
               srcs2 = sol_factb *cldmabc(i)*odds*tracer_mean*(1._r8-weight)/deltat & ! Liquid
                     + sol_factbi*cldmabc(i)*odds*tracer_mean*(weight)/deltat         ! Ice
            end if
#else
            odds=max( &
                 min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8) &
                 *scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
            srcs2 = sol_factb*cldmabc(i)*odds*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
                 +  sol_factbi*cldmabc(i)*odds*tracer(i,k)*(weight)/deltat    !ice
#endif


            !Note that using the temperature-determined weight doesn't make much sense here


            srcc = srcs1 + srcs2  ! convective tend by both processes
            finc = srcs1/(srcc + 1.e-36_r8) ! fraction in-cloud

            ! ****************** Stratiform ***********************
            ! now do the stratiform scavenging

            ! incloud scavenging

#ifdef MODAL_AERO
            ! rce 2010/05/01
            if ( is_strat_cloudborne ) then
               ! new code for stratiform incloud scav of cloudborne (modal) aerosol 
               ! >> use the 1st order cw to precip rate calculated in microphysics routine
               ! >> cloudborne aerosol resides in cloudy portion of grid cell, so do not apply "cldt" factor
             ! fracp = rate1ord_cw2pr_st(i,k)*deltat
             ! fracp = max(0._r8,min(1._r8,fracp))
               fracp = precs(i,k)*deltat/max(cwat(i,k)+precs(i,k)*deltat,1.e-12_r8) ! Sungsu. Mar.19.2010.
               fracp = max(0._r8,min(1._r8,fracp))
               srcs1 = sol_facti *fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! Liquid
                     + sol_factii*fracp*tracer(i,k)/deltat*(weight)          ! Ice
            else
               ! strat in-cloud removal only affects strat-cloudborne aerosol
               srcs1 = 0._r8
            end if
            ! end rce 2010/05/01
#else

            ! fracp is the fraction of cloud water converted to precip
            ! Sungsu modified fracp as the convectiv case.
            !        Below new formula by Sungsu of 'fracp' is necessary since 'cwat' is a LWC/IWC
            !        that has already precipitated out, that is, 'cwat' does not contain
            !        precipitation at all ! 
!           fracp =  precs(i,k)*deltat/max(cwat(i,k),1.e-12_r8)
            fracp =  precs(i,k)*deltat/max(cwat(i,k)+precs(i,k)*deltat,1.e-12_r8) ! Sungsu. Mar.19.2010.
            fracp = max(0._r8,min(1._r8,fracp))
!           fracp = 0.     ! for debug

            ! assume the corresponding amnt of tracer is removed
            !++mcb -- remove cldc; change cldt to cldv 
            !            srcs1 = (cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat
            !            srcs1 = cldv(i,k)*fracp*tracer(i,k)/deltat &
!           srcs1 = cldt(i,k)*fracp*tracer(i,k)/deltat            ! all condensate
!           Jan.02.2010. Sungsu : cldt --> cldt - cldc below.
            srcs1 = sol_facti*(cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! liquid
                 + sol_factii*(cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat*(weight)       ! ice
#endif


            ! below cloud scavenging

!           volume undergoing below cloud scavenging
!           cldmabs(i) = cldv(i,k)   ! precipitating volume
!           cldmabs(i) = cldt(i,k)   ! local cloud volume
            cldmabs(i) = cldvst(i,k) ! Stratiform precipitation area at the top interface of current layer

            ! Jan. 16. 2010. Sungsu
            ! cldmabs(i) = cldmabs(i) * cldovr_st(i) / max( 0.01_r8, cldovr_cu(i) + cldovr_st(i) )
            ! End by Sungsu

#ifdef MODAL_AERO
            ! rce 2010/05/01
            if ( is_strat_cloudborne ) then
               ! only strat in-cloud removal affects strat-cloudborne aerosol
               srcs2 = 0._r8
            else
               odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
               odds = max(min(1._r8,odds),0._r8)
               srcs2 = sol_factb *cldmabs(i)*odds*tracer_mean*(1._r8-weight)/deltat & ! Liquid
                     + sol_factbi*cldmabs(i)*odds*tracer_mean*(weight)/deltat         ! Ice
            end if
#else
            odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
            odds = max(min(1._r8,odds),0._r8)
            srcs2 =sol_factb*(cldmabs(i)*odds) *tracer(i,k)*(1._r8-weight)/deltat & ! liquid
                 + sol_factbi*(cldmabs(i)*odds) *tracer(i,k)*(weight)/deltat       ! ice


            !Note that using the temperature-determined weight doesn't make much sense here
#endif


            srcs = srcs1 + srcs2             ! total stratiform scavenging
            fins=srcs1/(srcs + 1.e-36_r8)    ! fraction taken by incloud processes

            ! make sure we dont take out more than is there
            ! ratio of amount available to amount removed
            rat(i) = tracer(i,k)/max(deltat*(srcc+srcs),1.e-36_r8)
            if (rat(i).lt.1._r8) then
               srcs = srcs*rat(i)
               srcc = srcc*rat(i)
            endif
            srct(i) = (srcc+srcs)*omsm

            
            ! fraction that is not removed within the cloud
            ! (assumed to be interstitial, and subject to convective transport)
            fracp = deltat*srct(i)/max(cldmabs(i)*tracer(i,k),1.e-36_r8)  ! amount removed
            fracp = max(0._r8,min(1._r8,fracp))
            fracis(i,k) = 1._r8 - fracp

            ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
            ! Sungsu added cumulus contribution in the below 3 blocks
         
            scavt(i,k) = -srct(i) + (fracev(i)*scavab(i)+fracev_cu(i)*scavabc(i))*gravit/pdel(i,k)
            iscavt(i,k) = -(srcc*finc + srcs*fins)*omsm

            if ( present(icscavt) ) icscavt(i,k) = -(srcc*finc) * omsm
            if ( present(isscavt) ) isscavt(i,k) = -(srcs*fins) * omsm
            if ( present(bcscavt) ) bcscavt(i,k) = -(srcc * (1-finc)) * omsm +  &
                 fracev_cu(i)*scavabc(i)*gravit/pdel(i,k)
            if ( present(bsscavt) ) bsscavt(i,k) = -(srcs * (1-fins)) * omsm +  &
                 fracev(i)*scavab(i)*gravit/pdel(i,k)

            dblchek(i) = tracer(i,k) + deltat*scavt(i,k)

            ! now keep track of scavenged mass and precip
            scavab(i) = scavab(i)*(1-fracev(i)) + srcs*pdel(i,k)/gravit
            precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdel(i,k)/gravit
            scavabc(i) = scavabc(i)*(1-fracev_cu(i)) + srcc*pdel(i,k)/gravit
            precabc(i) = precabc(i) + (cmfdqr(i,k) - evapc(i,k))*pdel(i,k)/gravit
            tracab(i) = tracab(i) + tracer(i,k)*pdel(i,k)/gravit

       ! Jan.16.2010. Sungsu
       ! Compute convective and stratiform precipitation areas at the base interface
       ! of current layer. These are for computing 'below cloud scavenging' in the 
       ! next layer below.

       ! cldovr_cu(i) = max( cldovr_cu(i), cldc(i,k) )
       ! cldovr_st(i) = max( cldovr_st(i), max( 0._r8, cldt(i,k) - cldc(i,k) ) )

       ! cldovr_cu(i) = max( 0._r8, min ( 1._r8, cldovr_cu(i) ) )
       ! cldovr_st(i) = max( 0._r8, min ( 1._r8, cldovr_st(i) ) )

       ! End by Sungsu

         end do ! End of i = 1, ncol

         found = .false.
         do i = 1,ncol
            if ( dblchek(i) < 0._r8 ) then
               found = .true.
               exit
            end if
         end do

         if ( found ) then
            do i = 1,ncol
               if (dblchek(i) .lt. 0._r8) then
                  write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
                       dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
               endif
            end do
         endif

      end do ! End of k = 1, pver

   end subroutine wetdepa_v2


!==============================================================================

! This is the frozen CAM4 version of wetdepa.



   subroutine wetdepa_v1( t, p, q, pdel, & 2
                       cldt, cldc, cmfdqr, conicw, precs, conds, &
                       evaps, cwat, tracer, deltat, &
                       scavt, iscavt, cldv, fracis, sol_fact, ncol, &
                       scavcoef,icscavt, isscavt, bcscavt, bsscavt, &
                       sol_facti_in, sol_factbi_in, sol_factii_in, &
                       sol_factic_in, sol_factiic_in )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! scavenging code for very soluble aerosols
      ! 
      ! Author: P. Rasch
      ! Modified by T. Bond 3/2003 to track different removals
      !-----------------------------------------------------------------------

      implicit none

      real(r8), intent(in) ::&
         t(pcols,pver),        &! temperature
         p(pcols,pver),        &! pressure
         q(pcols,pver),        &! moisture
         pdel(pcols,pver),     &! pressure thikness
         cldt(pcols,pver),    &! total cloud fraction
         cldc(pcols,pver),     &! convective cloud fraction
         cmfdqr(pcols,pver),   &! rate of production of convective precip
         conicw(pcols,pver),   &! convective cloud water
         cwat(pcols,pver),     &! cloud water amount 
         precs(pcols,pver),    &! rate of production of stratiform precip
         conds(pcols,pver),    &! rate of production of condensate
         evaps(pcols,pver),    &! rate of evaporation of precip
         cldv(pcols,pver),     &! total cloud fraction
         deltat,               &! time step
         tracer(pcols,pver)     ! trace species
      ! If subroutine is called with just sol_fact:
            ! sol_fact is used for both in- and below-cloud scavenging
      ! If subroutine is called with optional argument sol_facti_in:
            ! sol_fact  is used for below cloud scavenging
            ! sol_facti is used for in cloud scavenging
         real(r8), intent(in) :: sol_fact ! solubility factor (fraction of aer scavenged below & in, or just below or sol_facti_in is provided)
         real(r8), intent(in), optional :: sol_facti_in   ! solubility factor (frac of aerosol scavenged in cloud)
         real(r8), intent(in), optional :: sol_factbi_in  ! solubility factor (frac of aerosol scavenged below cloud by ice)
         real(r8), intent(in), optional :: sol_factii_in  ! solubility factor (frac of aerosol scavenged in cloud by ice)
         real(r8), intent(in), optional :: sol_factic_in(pcols,pver)  ! sol_facti_in for convective clouds
         real(r8), intent(in), optional :: sol_factiic_in ! sol_factii_in for convective clouds
         real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
         
      integer, intent(in) :: ncol

      real(r8), intent(out) ::&
         scavt(pcols,pver),    &! scavenging tend 
         iscavt(pcols,pver),   &! incloud scavenging tends
         fracis(pcols,pver)     ! fraction of species not scavenged

      real(r8), intent(out), optional ::    icscavt(pcols,pver)     ! incloud, convective
      real(r8), intent(out), optional ::    isscavt(pcols,pver)     ! incloud, stratiform
      real(r8), intent(out), optional ::    bcscavt(pcols,pver)     ! below cloud, convective
      real(r8), intent(out), optional ::    bsscavt(pcols,pver)     ! below cloud, stratiform

      ! local variables

      integer i                 ! x index
      integer k                 ! z index

      real(r8) adjfac               ! factor stolen from cmfmca
      real(r8) aqfrac               ! fraction of tracer in aqueous phase
      real(r8) cwatc                ! local convective total water amount 
      real(r8) cwats                ! local stratiform total water amount 
      real(r8) cwatp                ! local water amount falling from above precip
      real(r8) fracev(pcols)        ! fraction of precip from above that is evaporating
      real(r8) fracp                ! fraction of cloud water converted to precip
      real(r8) gafrac               ! fraction of tracer in gas phasea
      real(r8) hconst               ! henry's law solubility constant when equation is expressed
                                ! in terms of mixing ratios
      real(r8) mpla                 ! moles / liter H2O entering the layer from above
      real(r8) mplb                 ! moles / liter H2O leaving the layer below
      real(r8) omsm                 ! 1 - (a small number)
      real(r8) part                 !  partial pressure of tracer in atmospheres
      real(r8) patm                 ! total pressure in atmospheres
      real(r8) pdog                 ! work variable (pdel/gravit)
      real(r8) precabc(pcols)       ! conv precip from above (work array)
      real(r8) precabs(pcols)       ! strat precip from above (work array)
      real(r8) precbl               ! precip falling out of level (work array)
      real(r8) precmin              ! minimum convective precip causing scavenging
      real(r8) rat(pcols)           ! ratio of amount available to amount removed
      real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
      real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
      real(r8) srcc                 ! tend for convective rain
      real(r8) srcs                 ! tend for stratiform rain
      real(r8) srct(pcols)          ! work variable
      real(r8) tracab(pcols)        ! column integrated tracer amount
!      real(r8) vfall                ! fall speed of precip
      real(r8) fins                 ! fraction of rem. rate by strat rain
      real(r8) finc                 ! fraction of rem. rate by conv. rain
      real(r8) srcs1                ! work variable
      real(r8) srcs2                ! work variable
      real(r8) tc                   ! temp in celcius
      real(r8) weight               ! fraction of condensate which is ice
      real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
      real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
      real(r8) odds                 ! limit on removal rate (proportional to prec)
      real(r8) dblchek(pcols)
      logical :: found

      real(r8) sol_facti,  sol_factb  ! in cloud and below cloud fraction of aerosol scavenged
      real(r8) sol_factii, sol_factbi ! in cloud and below cloud fraction of aerosol scavenged by ice
      real(r8) sol_factic(pcols,pver)             ! sol_facti for convective clouds
      real(r8) sol_factiic            ! sol_factii for convective clouds
      ! sol_factic & solfact_iic added for MODAL_AERO.  
      ! For stratiform cloud, cloudborne aerosol is treated explicitly,
      !    and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
      ! For convective cloud, cloudborne aerosol is not treated explicitly,
      !    and sol_factic is 1.0 for both cloudborne and interstitial.

      ! ------------------------------------------------------------------------
!      omsm = 1.-1.e-10          ! used to prevent roundoff errors below zero
      omsm = 1._r8-2*epsilon(1._r8) ! used to prevent roundoff errors below zero
      precmin =  0.1_r8/8.64e4_r8      ! set critical value to 0.1 mm/day in kg/m2/s

      adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme

      ! assume 4 m/s fall speed currently (should be improved)
!      vfall = 4.
	
      ! default (if other sol_facts aren't in call, set all to required sol_fact
      sol_facti = sol_fact
      sol_factb = sol_fact
      sol_factii = sol_fact
      sol_factbi = sol_fact

      if ( present(sol_facti_in) )  sol_facti = sol_facti_in
      if ( present(sol_factii_in) )  sol_factii = sol_factii_in
      if ( present(sol_factbi_in) )  sol_factbi = sol_factbi_in

      sol_factic  = sol_facti
      sol_factiic = sol_factii
      if ( present(sol_factic_in ) )  sol_factic  = sol_factic_in
      if ( present(sol_factiic_in) )  sol_factiic = sol_factiic_in

      ! this section of code is for highly soluble aerosols,
      ! the assumption is that within the cloud that
      ! all the tracer is in the cloud water
      !
      ! for both convective and stratiform clouds, 
      ! the fraction of cloud water converted to precip defines
      ! the amount of tracer which is pulled out.
      !

      do i = 1,pcols
         precabs(i) = 0
         precabc(i) = 0
         scavab(i) = 0
         scavabc(i) = 0
         tracab(i) = 0
         cldmabs(i) = 0
         cldmabc(i) = 0
      end do

      do k = 1,pver
         do i = 1,ncol
            tc     = t(i,k) - tmelt
            weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
            weight = 0._r8                                 ! assume no ice

            pdog = pdel(i,k)/gravit

            ! ****************** Evaporation **************************
            ! calculate the fraction of strat precip from above 
            !                 which evaporates within this layer
            fracev(i) = evaps(i,k)*pdel(i,k)/gravit &
                     /max(1.e-12_r8,precabs(i))

            ! trap to ensure reasonable ratio bounds
            fracev(i) = max(0._r8,min(1._r8,fracev(i)))

            ! ****************** Convection ***************************
            ! now do the convective scavenging

            ! set odds proportional to fraction of the grid box that is swept by the 
            ! precipitation =precabc/rhoh20*(area of sphere projected on plane
            !                                /volume of sphere)*deltat
            ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
            ! unless the fraction of the area that is cloud is less than odds, in which
            ! case use the cloud fraction (assumes precabs is in kg/m2/s)
            ! is really: precabs*3/4/1000./1e-3*deltat
            ! here I use .1 from Balkanski
            !
            ! use a local rate of convective rain production for incloud scav
            !odds=max(min(1._r8, &
            !     cmfdqr(i,k)*pdel(i,k)/gravit*0.1_r8*deltat),0._r8)
            !++mcb -- change cldc to cldt; change cldt to cldv (9/17/96)
            !            srcs1 =  cldt(i,k)*odds*tracer(i,k)*(1.-weight) &
            ! srcs1 =  cldv(i,k)*odds*tracer(i,k)*(1.-weight) &
            !srcs1 =  cldc(i,k)*odds*tracer(i,k)*(1.-weight) &
            !         /deltat 

            ! fraction of convective cloud water converted to rain
            fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,conicw(i,k))
            ! note cmfdrq can be negative from evap of rain, so constrain it
            fracp = max(min(1._r8,fracp),0._r8)
            ! remove that amount from within the convective area
!           srcs1 = cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat ! liquid only
!           srcs1 = cldc(i,k)*fracp*tracer(i,k)/deltat             ! any condensation
!           srcs1 = 0.
            srcs1 = sol_factic(i,k)*cldt(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat &  ! liquid
                 +  sol_factiic*cldt(i,k)*fracp*tracer(i,k)*(weight)/deltat      ! ice


            !--mcb

            ! scavenge below cloud

            !            cldmabc(i) = max(cldc(i,k),cldmabc(i))
            !            cldmabc(i) = max(cldt(i,k),cldmabc(i))
            cldmabc(i) = max(cldv(i,k),cldmabc(i))
            cldmabc(i) = cldv(i,k)

            odds=max( &
                 min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8) &
                 *scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
            srcs2 = sol_factb*cldmabc(i)*odds*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
                 +  sol_factbi*cldmabc(i)*odds*tracer(i,k)*(weight)/deltat    !ice
            !Note that using the temperature-determined weight doesn't make much sense here


            srcc = srcs1 + srcs2  ! convective tend by both processes
            finc = srcs1/(srcc + 1.e-36_r8) ! fraction in-cloud

            ! ****************** Stratiform ***********************
            ! now do the stratiform scavenging

            ! incloud scavenging

            ! fracp is the fraction of cloud water converted to precip
            fracp =  precs(i,k)*deltat/max(cwat(i,k),1.e-12_r8)
            fracp = max(0._r8,min(1._r8,fracp))
!           fracp = 0.     ! for debug

            ! assume the corresponding amnt of tracer is removed
            !++mcb -- remove cldc; change cldt to cldv 
            !            srcs1 = (cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat
            !            srcs1 = cldv(i,k)*fracp*tracer(i,k)/deltat &
!           srcs1 = cldt(i,k)*fracp*tracer(i,k)/deltat            ! all condensate
            srcs1 = sol_facti*cldt(i,k)*fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! liquid
                 + sol_factii*cldt(i,k)*fracp*tracer(i,k)/deltat*(weight)       ! ice


            ! below cloud scavenging

!           volume undergoing below cloud scavenging
            cldmabs(i) = cldv(i,k)   ! precipitating volume
!           cldmabs(i) = cldt(i,k)   ! local cloud volume

            odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
            odds = max(min(1._r8,odds),0._r8)
            srcs2 =sol_factb*(cldmabs(i)*odds) *tracer(i,k)*(1._r8-weight)/deltat & ! liquid
                 + sol_factbi*(cldmabs(i)*odds) *tracer(i,k)*(weight)/deltat       ! ice
            !Note that using the temperature-determined weight doesn't make much sense here


            srcs = srcs1 + srcs2             ! total stratiform scavenging
            fins=srcs1/(srcs + 1.e-36_r8)    ! fraction taken by incloud processes

            ! make sure we dont take out more than is there
            ! ratio of amount available to amount removed
            rat(i) = tracer(i,k)/max(deltat*(srcc+srcs),1.e-36_r8)
            if (rat(i).lt.1._r8) then
               srcs = srcs*rat(i)
               srcc = srcc*rat(i)
            endif
            srct(i) = (srcc+srcs)*omsm

            
            ! fraction that is not removed within the cloud
            ! (assumed to be interstitial, and subject to convective transport)
            fracp = deltat*srct(i)/max(cldmabs(i)*tracer(i,k),1.e-36_r8)  ! amount removed
            fracp = max(0._r8,min(1._r8,fracp))
            fracis(i,k) = 1._r8 - fracp

            ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
            scavt(i,k) = -srct(i) + fracev(i)*scavab(i)*gravit/pdel(i,k)
            iscavt(i,k) = -(srcc*finc + srcs*fins)*omsm

            if ( present(icscavt) ) icscavt(i,k) = -(srcc*finc) * omsm
            if ( present(isscavt) ) isscavt(i,k) = -(srcs*fins) * omsm
            if ( present(bcscavt) ) bcscavt(i,k) = -(srcc * (1-finc)) * omsm
            if ( present(bsscavt) ) bsscavt(i,k) = -(srcs * (1-fins)) * omsm +  &
                 fracev(i)*scavab(i)*gravit/pdel(i,k)

            dblchek(i) = tracer(i,k) + deltat*scavt(i,k)

            ! now keep track of scavenged mass and precip
            scavab(i) = scavab(i)*(1-fracev(i)) + srcs*pdel(i,k)/gravit
            precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdel(i,k)/gravit
            scavabc(i) = scavabc(i) + srcc*pdel(i,k)/gravit
            precabc(i) = precabc(i) + (cmfdqr(i,k))*pdel(i,k)/gravit
            tracab(i) = tracab(i) + tracer(i,k)*pdel(i,k)/gravit

         end do

         found = .false.
         do i = 1,ncol
            if ( dblchek(i) < 0._r8 ) then
               found = .true.
               exit
            end if
         end do

         if ( found ) then
            do i = 1,ncol
               if (dblchek(i) .lt. 0._r8) then
                  write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
                       dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
               endif
            end do
         endif

      end do

   end subroutine wetdepa_v1

!==============================================================================

! wetdepg is currently being used for both CAM4 and CAM5 by making use of the
! cam_physpkg_is method.


   subroutine wetdepg( t, p, q, pdel, &,1
                       cldt, cldc, cmfdqr, evapc, precs, evaps, &
                       rain, cwat, tracer, deltat, molwt, &
                       solconst, scavt, iscavt, cldv, icwmr1, &
                       icwmr2, fracis, ncol )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! scavenging of gas phase constituents by henry's law
      ! 
      ! Author: P. Rasch
      !-----------------------------------------------------------------------

      real(r8), intent(in) ::&
         t(pcols,pver),        &! temperature
         p(pcols,pver),        &! pressure
         q(pcols,pver),        &! moisture
         pdel(pcols,pver),     &! pressure thikness
         cldt(pcols,pver),     &! total cloud fraction
         cldc(pcols,pver),     &! convective cloud fraction
         cmfdqr(pcols,pver),   &! rate of production of convective precip
         rain (pcols,pver),    &! total rainwater mixing ratio
         cwat(pcols,pver),     &! cloud water amount 
         precs(pcols,pver),    &! rate of production of stratiform precip
         evaps(pcols,pver),    &! rate of evaporation of precip
! Sungsu
         evapc(pcols,pver),    &! Rate of evaporation of convective precipitation
! Sungsu 
         cldv(pcols,pver),     &! estimate of local volume occupied by clouds
         icwmr1 (pcols,pver),  &! in cloud water mixing ration for zhang scheme
         icwmr2 (pcols,pver),  &! in cloud water mixing ration for hack  scheme
         deltat,               &! time step
         tracer(pcols,pver),   &! trace species
         molwt                  ! molecular weights

      integer, intent(in) :: ncol

      real(r8) &
         solconst(pcols,pver)   ! Henry's law coefficient

      real(r8), intent(out) ::&
         scavt(pcols,pver),    &! scavenging tend 
         iscavt(pcols,pver),   &! incloud scavenging tends
         fracis(pcols, pver)    ! fraction of constituent that is insoluble

      ! local variables

      integer i                 ! x index
      integer k                 ! z index

      real(r8) adjfac               ! factor stolen from cmfmca
      real(r8) aqfrac               ! fraction of tracer in aqueous phase
      real(r8) cwatc                ! local convective total water amount 
      real(r8) cwats                ! local stratiform total water amount 
      real(r8) cwatl                ! local cloud liq water amount 
      real(r8) cwatp                ! local water amount falling from above precip
      real(r8) cwatpl               ! local water amount falling from above precip (liq)
      real(r8) cwatt                ! local sum of strat + conv total water amount 
      real(r8) cwatti               ! cwatt/cldv = cloudy grid volume mixing ratio
      real(r8) fracev               ! fraction of precip from above that is evaporating
      real(r8) fracp                ! fraction of cloud water converted to precip
      real(r8) gafrac               ! fraction of tracer in gas phasea
      real(r8) hconst               ! henry's law solubility constant when equation is expressed
                                ! in terms of mixing ratios
      real(r8) mpla                 ! moles / liter H2O entering the layer from above
      real(r8) mplb                 ! moles / liter H2O leaving the layer below
      real(r8) omsm                 ! 1 - (a small number)
      real(r8) part                 !  partial pressure of tracer in atmospheres
      real(r8) patm                 ! total pressure in atmospheres
      real(r8) pdog                 ! work variable (pdel/gravit)
      real(r8) precab(pcols)        ! precip from above (work array)
      real(r8) precbl               ! precip work variable
      real(r8) precxx               ! precip work variable
      real(r8) precxx2               !
      real(r8) precic               ! precip work variable
      real(r8) rat                  ! ratio of amount available to amount removed
      real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
      real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
      !      real(r8) vfall                ! fall speed of precip
      real(r8) scavmax              ! an estimate of the max tracer avail for removal
      real(r8) scavbl               ! flux removed at bottom of layer
      real(r8) fins                 ! in cloud fraction removed by strat rain
      real(r8) finc                 ! in cloud fraction removed by conv rain
      real(r8) rate                 ! max removal rate estimate
      real(r8) scavlimt             ! limiting value 1
      real(r8) scavt1               ! limiting value 2
      real(r8) scavin               ! scavenging by incloud processes
      real(r8) scavbc               ! scavenging by below cloud processes
      real(r8) tc
      real(r8) weight               ! ice fraction
      real(r8) wtpl                 ! work variable
      real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
      real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
      !-----------------------------------------------------------

      omsm = 1._r8-2*epsilon(1._r8)   ! used to prevent roundoff errors below zero

      adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme

      ! assume 4 m/s fall speed currently (should be improved)
      !      vfall = 4.

      ! zero accumulators
      do i = 1,pcols
         precab(i) = 1.e-36_r8
         scavab(i) = 0._r8
         cldmabs(i) = 0._r8
      end do

      do k = 1,pver
         do i = 1,ncol

            tc     = t(i,k) - tmelt
            weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice

            cldmabs(i) = max(cldmabs(i),cldt(i,k))

            ! partitioning coefs for gas and aqueous phase
            !              take as a cloud water amount, the sum of the stratiform amount
            !              plus the convective rain water amount 

            ! convective amnt is just the local precip rate from the hack scheme
            !              since there is no storage of water, this ignores that falling from above
            !            cwatc = cmfdqr(i,k)*deltat/adjfac
            !++mcb -- test cwatc
            cwatc = (icwmr1(i,k) + icwmr2(i,k)) * (1._r8-weight)
            !--mcb 

            ! strat cloud water amount and also ignore the part falling from above
            cwats = cwat(i,k) 

            ! cloud water as liq
            !++mcb -- add cwatc later (in cwatti)
            !            cwatl = (1.-weight)*(cwatc+cwats)
            cwatl = (1._r8-weight)*cwats
            ! cloud water as ice
            !*not used        cwati = weight*(cwatc+cwats)

            ! total suspended condensate as liquid
            cwatt = cwatl + rain(i,k)

            ! incloud version 
            !++mcb -- add cwatc here
            cwatti = cwatt/max(cldv(i,k), 0.00001_r8) + cwatc

            ! partitioning terms
            patm = p(i,k)/1.013e5_r8 ! pressure in atmospheres
            hconst = molwta*patm*solconst(i,k)*cwatti/rhoh2o
            aqfrac = hconst/(1._r8+hconst)
            gafrac = 1/(1._r8+hconst)
            fracis(i,k) = gafrac


            ! partial pressure of the tracer in the gridbox in atmospheres
            part = patm*gafrac*tracer(i,k)*molwta/molwt

            ! use henrys law to give moles tracer /liter of water
            ! in this volume 
            ! then convert to kg tracer /liter of water (kg tracer / kg water)
            mplb = solconst(i,k)*part*molwt/1000._r8


            pdog = pdel(i,k)/gravit

            ! this part of precip will be carried downward but at a new molarity of mpl 
            precic = pdog*(precs(i,k) + cmfdqr(i,k))

            ! we cant take out more than entered, plus that available in the cloud
            !                  scavmax = scavab(i)+tracer(i,k)*cldt(i,k)/deltat*pdog
            scavmax = scavab(i)+tracer(i,k)*cldv(i,k)/deltat*pdog

            ! flux of tracer by incloud processes
            scavin = precic*(1._r8-weight)*mplb

            ! fraction of precip which entered above that leaves below
            if (cam_physpkg_is('cam5')) then
               ! Sungsu added evaporation of convective precipitation below.
               precxx = precab(i)-pdog*(evaps(i,k)+evapc(i,k))
            else
               precxx = precab(i)-pdog*evaps(i,k)
            end if
            precxx = max (precxx,0.0_r8)

            ! flux of tracer by below cloud processes
            !++mcb -- removed wtpl because it is now not assigned and previously
            !          when it was assigned it was unnecessary:  if(tc.gt.0)wtpl=1
            if (tc.gt.0) then
               !               scavbc = precxx*wtpl*mplb ! if liquid
               scavbc = precxx*mplb ! if liquid
            else
               precxx2=max(precxx,1.e-36_r8)
               scavbc = scavab(i)*precxx2/(precab(i)) ! if ice
            endif

            scavbl = min(scavbc + scavin, scavmax)

            ! first guess assuming that henries law works
            scavt1 = (scavab(i)-scavbl)/pdog*omsm

            ! pjr this should not be required, but we put it in to make sure we cant remove too much
            ! remember, scavt1 is generally negative (indicating removal)
            scavt1 = max(scavt1,-tracer(i,k)*cldv(i,k)/deltat)

            !++mcb -- remove this limitation for gas species
            !c use the dana and hales or balkanski limit on scavenging
            !c            rate = precab(i)*0.1
            !            rate = (precic + precxx)*0.1
            !            scavlimt = -tracer(i,k)*cldv(i,k)
            !     $           *rate/(1.+rate*deltat)

            !            scavt(i,k) = max(scavt1, scavlimt)

            ! instead just set scavt to scavt1
            scavt(i,k) = scavt1
            !--mcb

            ! now update the amount leaving the layer
            scavbl = scavab(i) - scavt(i,k)*pdog 

            ! in cloud amount is that formed locally over the total flux out bottom
            fins = scavin/(scavin + scavbc + 1.e-36_r8)
            iscavt(i,k) = scavt(i,k)*fins

            scavab(i) = scavbl
            precab(i) = max(precxx + precic,1.e-36_r8)

        
            
         end do
      end do
      
   end subroutine wetdepg

!##############################################################################

end module wetdep