chem/module_cam_mam_wetdep.F

   1 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   2 #define WRF_PORT
   3 #define MODAL_AERO
   4 module wetdep
   5
   6 !-----------------------------------------------------------------------
   7 !
   8 ! Wet deposition routines for both aerosols and gas phase constituents.
   9 !
  10 !-----------------------------------------------------------------------
  11
  12 use shr_kind_mod, only: r8 => shr_kind_r8
  13 #ifndef WRF_PORT
  14 use ppgrid,       only: pcols, pver
  15 #endif
  16 use physconst,    only: gravit, rair, tmelt
  17 #ifndef WRF_PORT
  18 use phys_control, only: cam_physpkg_is
  19 use cam_logfile,  only: iulog
  20 #else
  21 use module_cam_support, only: pcols, pver, iulog
  22 #endif
  23
  24 implicit none
  25 save
  26 private
  27
  28 public :: wetdepa_v1  ! scavenging codes for very soluble aerosols -- CAM4 version
  29 public :: wetdepa_v2  ! scavenging codes for very soluble aerosols -- CAM5 version
  30 #ifndef WRF_PORT
  31 public :: wetdepg     ! scavenging of gas phase constituents by henry's law
  32 #endif
  33 public :: clddiag     ! calc of cloudy volume and rain mixing ratio
  34
  35 real(r8), parameter :: cmftau = 3600._r8
  36 real(r8), parameter :: rhoh2o = 1000._r8            ! density of water
  37 real(r8), parameter :: molwta = 28.97_r8            ! molecular weight dry air gm/mole
  38
  39 !==============================================================================
  40 contains
  41 !==============================================================================
  42
  43 subroutine clddiag(t, pmid, pdel, cmfdqr, evapc, &
  44                    cldt, cldcu, cldst, cme, evapr, &
  45                    prain, cldv, cldvcu, cldvst, rain, &
  46                    ncol)
  47
  48    ! ------------------------------------------------------------------------------------
  49    ! Estimate the cloudy volume which is occupied by rain or cloud water as
  50    ! the max between the local cloud amount or the
  51    ! sum above of (cloud*positive precip production)      sum total precip from above
  52    !              ----------------------------------   x ------------------------
  53    ! sum above of     (positive precip           )        sum positive precip from above
  54    ! Author: P. Rasch
  55    !         Sungsu Park. Mar.2010
  56    ! ------------------------------------------------------------------------------------
  57
  58    ! Input arguments:
  59    real(r8), intent(in) :: t(pcols,pver)        ! temperature (K)
  60    real(r8), intent(in) :: pmid(pcols,pver)     ! pressure at layer midpoints
  61    real(r8), intent(in) :: pdel(pcols,pver)     ! pressure difference across layers
  62    real(r8), intent(in) :: cmfdqr(pcols,pver)   ! dq/dt due to convective rainout
  63    real(r8), intent(in) :: evapc(pcols,pver)    ! Evaporation rate of convective precipitation ( >= 0 )
  64    real(r8), intent(in) :: cldt(pcols,pver)    ! total cloud fraction
  65    real(r8), intent(in) :: cldcu(pcols,pver)    ! Cumulus cloud fraction
  66    real(r8), intent(in) :: cldst(pcols,pver)    ! Stratus cloud fraction
  67    real(r8), intent(in) :: cme(pcols,pver)      ! rate of cond-evap within the cloud
  68    real(r8), intent(in) :: evapr(pcols,pver)    ! rate of evaporation of falling precipitation (kg/kg/s)
  69    real(r8), intent(in) :: prain(pcols,pver)    ! rate of conversion of condensate to precipitation (kg/kg/s)
  70    integer, intent(in) :: ncol
  71
  72    ! Output arguments:
  73    real(r8), intent(out) :: cldv(pcols,pver)     ! fraction occupied by rain or cloud water
  74    real(r8), intent(out) :: cldvcu(pcols,pver)   ! Convective precipitation volume
  75    real(r8), intent(out) :: cldvst(pcols,pver)   ! Stratiform precipitation volume
  76    real(r8), intent(out) :: rain(pcols,pver)     ! mixing ratio of rain (kg/kg)
  77
  78    ! Local variables:
  79    integer  i, k
  80    real(r8) convfw         ! used in fallspeed calculation; taken from findmcnew
  81    real(r8) sumppr(pcols)        ! precipitation rate (kg/m2-s)
  82    real(r8) sumpppr(pcols)       ! sum of positive precips from above
  83    real(r8) cldv1(pcols)         ! precip weighted cloud fraction from above
  84    real(r8) lprec                ! local production rate of precip (kg/m2/s)
  85    real(r8) lprecp               ! local production rate of precip (kg/m2/s) if positive
  86    real(r8) rho                  ! air density
  87    real(r8) vfall
  88    real(r8) sumppr_cu(pcols)     ! Convective precipitation rate (kg/m2-s)
  89    real(r8) sumpppr_cu(pcols)    ! Sum of positive convective precips from above
  90    real(r8) cldv1_cu(pcols)      ! Convective precip weighted convective cloud fraction from above
  91    real(r8) lprec_cu             ! Local production rate of convective precip (kg/m2/s)
  92    real(r8) lprecp_cu            ! Local production rate of convective precip (kg/m2/s) if positive
  93    real(r8) sumppr_st(pcols)     ! Stratiform precipitation rate (kg/m2-s)
  94    real(r8) sumpppr_st(pcols)    ! Sum of positive stratiform precips from above
  95    real(r8) cldv1_st(pcols)      ! Stratiform precip weighted stratiform cloud fraction from above
  96    real(r8) lprec_st             ! Local production rate of stratiform precip (kg/m2/s)
  97    real(r8) lprecp_st            ! Local production rate of stratiform precip (kg/m2/s) if positive
  98    ! -----------------------------------------------------------------------
  99
 100    convfw = 1.94_r8*2.13_r8*sqrt(rhoh2o*gravit*2.7e-4_r8)
 101    do i=1,ncol
 102       sumppr(i) = 0._r8
 103       cldv1(i) = 0._r8
 104       sumpppr(i) = 1.e-36_r8
 105       sumppr_cu(i)  = 0._r8
 106       cldv1_cu(i)   = 0._r8
 107       sumpppr_cu(i) = 1.e-36_r8
 108       sumppr_st(i)  = 0._r8
 109       cldv1_st(i)   = 0._r8
 110       sumpppr_st(i) = 1.e-36_r8
 111    end do
 112
 113    do k = 1,pver
 114       do i = 1,ncol
 115          cldv(i,k) = &
 116             max(min(1._r8, &
 117             cldv1(i)/sumpppr(i) &
 118             )*sumppr(i)/sumpppr(i), &
 119             cldt(i,k) &
 120             )
 121          lprec = pdel(i,k)/gravit &
 122             *(prain(i,k)+cmfdqr(i,k)-evapr(i,k))
 123          lprecp = max(lprec,1.e-30_r8)
 124          cldv1(i) = cldv1(i)  + cldt(i,k)*lprecp
 125          sumppr(i) = sumppr(i) + lprec
 126          sumpppr(i) = sumpppr(i) + lprecp
 127
 128          ! For convective precipitation volume at the top interface of each layer. Neglect the current layer.
 129          cldvcu(i,k)   = max(min(1._r8,cldv1_cu(i)/sumpppr_cu(i))*(sumppr_cu(i)/sumpppr_cu(i)),0._r8)
 130          lprec_cu      = (pdel(i,k)/gravit)*(cmfdqr(i,k)-evapc(i,k))
 131          lprecp_cu     = max(lprec_cu,1.e-30_r8)
 132          cldv1_cu(i)   = cldv1_cu(i) + cldcu(i,k)*lprecp_cu
 133          sumppr_cu(i)  = sumppr_cu(i) + lprec_cu
 134          sumpppr_cu(i) = sumpppr_cu(i) + lprecp_cu
 135
 136          ! For stratiform precipitation volume at the top interface of each layer. Neglect the current layer.
 137          cldvst(i,k)   = max(min(1._r8,cldv1_st(i)/sumpppr_st(i))*(sumppr_st(i)/sumpppr_st(i)),0._r8)
 138          lprec_st      = (pdel(i,k)/gravit)*(prain(i,k)-evapr(i,k))
 139          lprecp_st     = max(lprec_st,1.e-30_r8)
 140          cldv1_st(i)   = cldv1_st(i) + cldst(i,k)*lprecp_st
 141          sumppr_st(i)  = sumppr_st(i) + lprec_st
 142          sumpppr_st(i) = sumpppr_st(i) + lprecp_st
 143
 144          rain(i,k) = 0._r8
 145          if(t(i,k) .gt. tmelt) then
 146             rho = pmid(i,k)/(rair*t(i,k))
 147             vfall = convfw/sqrt(rho)
 148             rain(i,k) = sumppr(i)/(rho*vfall)
 149             if (rain(i,k).lt.1.e-14_r8) rain(i,k) = 0._r8
 150          endif
 151       end do
 152    end do
 153
 154 end subroutine clddiag
 155
 156 !==============================================================================
 157
 158 ! This is the CAM5 version of wetdepa.
 159
 160 subroutine wetdepa_v2(t, p, q, pdel, &
 161                    cldt, cldc, cmfdqr, evapc, conicw, precs, conds, &
 162                        evaps, cwat, tracer, deltat, &
 163                        scavt, iscavt, cldv, cldvcu, cldvst, dlf, fracis, sol_fact, ncol, &
 164                        scavcoef,  &
 165 #ifdef MODAL_AERO
 166                        is_strat_cloudborne, rate1ord_cw2pr_st, qqcw, f_act_conv, &     ! rce 2010/05/01
 167 #endif
 168                        icscavt, isscavt, bcscavt, bsscavt, &
 169                        sol_facti_in, sol_factbi_in, sol_factii_in, &
 170                        sol_factic_in, sol_factiic_in )
 171
 172       !-----------------------------------------------------------------------
 173       ! Purpose:
 174       ! scavenging code for very soluble aerosols
 175       !
 176       ! Author: P. Rasch
 177       ! Modified by T. Bond 3/2003 to track different removals
 178       ! Sungsu Park. Mar.2010 : Impose consistencies with a few changes in physics.
 179       !-----------------------------------------------------------------------
 180
 181
 182
 183       implicit none
 184
 185       real(r8), intent(in) ::&
 186          t(pcols,pver),        &! temperature
 187          p(pcols,pver),        &! pressure
 188          q(pcols,pver),        &! moisture
 189          pdel(pcols,pver),     &! pressure thikness
 190          cldt(pcols,pver),    &! total cloud fraction
 191          cldc(pcols,pver),     &! convective cloud fraction
 192          cmfdqr(pcols,pver),   &! rate of production of convective precip
 193 ! Sungsu
 194          evapc(pcols,pver),    &! Evaporation rate of convective precipitation
 195 ! Sungsu
 196          conicw(pcols,pver),   &! convective cloud water
 197          cwat(pcols,pver),     &! cloud water amount
 198          precs(pcols,pver),    &! rate of production of stratiform precip
 199          conds(pcols,pver),    &! rate of production of condensate
 200          evaps(pcols,pver),    &! rate of evaporation of precip
 201          cldv(pcols,pver),     &! total cloud fraction
 202 ! Sungsu
 203          cldvcu(pcols,pver),   &! Convective precipitation area at the top interface of each layer
 204          cldvst(pcols,pver),   &! Stratiform precipitation area at the top interface of each layer
 205          dlf(pcols,pver),      &! Detrainment of convective condensate [kg/kg/s]
 206 ! Sungsu
 207          deltat,               &! time step
 208          tracer(pcols,pver)     ! trace species
 209       ! If subroutine is called with just sol_fact:
 210             ! sol_fact is used for both in- and below-cloud scavenging
 211       ! If subroutine is called with optional argument sol_facti_in:
 212             ! sol_fact  is used for below cloud scavenging
 213             ! sol_facti is used for in cloud scavenging
 214          real(r8), intent(in) :: sol_fact ! solubility factor (fraction of aer scavenged below & in, or just below or sol_facti_in is provided)
 215          integer, intent(in) :: ncol
 216          real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
 217 #ifdef MODAL_AERO
 218          ! rce 2010/05/01
 219          ! is_strat_cloudborne = .true. if tracer is stratiform-cloudborne aerosol; else .false.
 220          logical, intent(in), optional :: is_strat_cloudborne
 221          ! rate1ord_cw2pr_st = 1st order rate for strat cw to precip (1/s)
 222          real(r8), intent(in), optional  :: rate1ord_cw2pr_st(pcols,pver)
 223          ! qqcw = strat-cloudborne aerosol corresponding to tracer when is_strat_cloudborne==.false.; else 0.0
 224          real(r8), intent(in), optional  :: qqcw(pcols,pver)
 225          ! f_act_conv = conv-cloud activation fraction when is_strat_cloudborne==.false.; else 0.0
 226          real(r8), intent(in), optional  :: f_act_conv(pcols,pver)
 227          ! end rce 2010/05/01
 228 #endif
 229
 230          real(r8), intent(in), optional :: sol_facti_in   ! solubility factor (frac of aerosol scavenged in cloud)
 231          real(r8), intent(in), optional :: sol_factbi_in  ! solubility factor (frac of aerosol scavenged below cloud by ice)
 232          real(r8), intent(in), optional :: sol_factii_in  ! solubility factor (frac of aerosol scavenged in cloud by ice)
 233          real(r8), intent(in), optional :: sol_factic_in(pcols,pver)  ! sol_facti_in for convective clouds
 234          real(r8), intent(in), optional :: sol_factiic_in ! sol_factii_in for convective clouds
 235
 236       real(r8), intent(out) ::&
 237          scavt(pcols,pver),    &! scavenging tend
 238          iscavt(pcols,pver),   &! incloud scavenging tends
 239          fracis(pcols,pver)     ! fraction of species not scavenged
 240
 241       real(r8), intent(out), optional ::    icscavt(pcols,pver)     ! incloud, convective
 242       real(r8), intent(out), optional ::    isscavt(pcols,pver)     ! incloud, stratiform
 243       real(r8), intent(out), optional ::    bcscavt(pcols,pver)     ! below cloud, convective
 244       real(r8), intent(out), optional ::    bsscavt(pcols,pver)     ! below cloud, stratiform
 245
 246       ! local variables
 247
 248       integer i                 ! x index
 249       integer k                 ! z index
 250
 251       real(r8) adjfac               ! factor stolen from cmfmca
 252       real(r8) aqfrac               ! fraction of tracer in aqueous phase
 253       real(r8) cwatc                ! local convective total water amount
 254       real(r8) cwats                ! local stratiform total water amount
 255       real(r8) cwatp                ! local water amount falling from above precip
 256       real(r8) fracev(pcols)        ! fraction of precip from above that is evaporating
 257 ! Sungsu
 258       real(r8) fracev_cu(pcols)     ! Fraction of convective precip from above that is evaporating
 259 ! Sungsu
 260       real(r8) fracp                ! fraction of cloud water converted to precip
 261       real(r8) gafrac               ! fraction of tracer in gas phasea
 262       real(r8) hconst               ! henry's law solubility constant when equation is expressed
 263                                 ! in terms of mixing ratios
 264       real(r8) mpla                 ! moles / liter H2O entering the layer from above
 265       real(r8) mplb                 ! moles / liter H2O leaving the layer below
 266       real(r8) omsm                 ! 1 - (a small number)
 267       real(r8) part                 !  partial pressure of tracer in atmospheres
 268       real(r8) patm                 ! total pressure in atmospheres
 269       real(r8) pdog                 ! work variable (pdel/gravit)
 270       real(r8) precabc(pcols)       ! conv precip from above (work array)
 271       real(r8) precabs(pcols)       ! strat precip from above (work array)
 272       real(r8) precbl               ! precip falling out of level (work array)
 273       real(r8) precmin              ! minimum convective precip causing scavenging
 274       real(r8) rat(pcols)           ! ratio of amount available to amount removed
 275       real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
 276       real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
 277       real(r8) srcc                 ! tend for convective rain
 278       real(r8) srcs                 ! tend for stratiform rain
 279       real(r8) srct(pcols)          ! work variable
 280       real(r8) tracab(pcols)        ! column integrated tracer amount
 281 !      real(r8) vfall                ! fall speed of precip
 282       real(r8) fins                 ! fraction of rem. rate by strat rain
 283       real(r8) finc                 ! fraction of rem. rate by conv. rain
 284       real(r8) srcs1                ! work variable
 285       real(r8) srcs2                ! work variable
 286       real(r8) tc                   ! temp in celcius
 287       real(r8) weight               ! fraction of condensate which is ice
 288       real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
 289       real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
 290       real(r8) odds                 ! limit on removal rate (proportional to prec)
 291       real(r8) dblchek(pcols)
 292       logical :: found
 293
 294     ! Jan.16.2009. Sungsu for wet scavenging below clouds.
 295     ! real(r8) cldovr_cu(pcols)     ! Convective precipitation area at the base of each layer
 296     ! real(r8) cldovr_st(pcols)     ! Stratiform precipitation area at the base of each layer
 297
 298       real(r8) tracer_incu
 299       real(r8) tracer_mean
 300
 301     ! End by Sungsu
 302
 303       real(r8) sol_facti,  sol_factb  ! in cloud and below cloud fraction of aerosol scavenged
 304       real(r8) sol_factii, sol_factbi ! in cloud and below cloud fraction of aerosol scavenged by ice
 305       real(r8) sol_factic(pcols,pver)             ! sol_facti for convective clouds
 306       real(r8) sol_factiic            ! sol_factii for convective clouds
 307       ! sol_factic & solfact_iic added for MODAL_AERO.
 308       ! For stratiform cloud, cloudborne aerosol is treated explicitly,
 309       !    and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
 310       ! For convective cloud, cloudborne aerosol is not treated explicitly,
 311       !    and sol_factic is 1.0 for both cloudborne and interstitial.
 312
 313       ! ------------------------------------------------------------------------
 314 !      omsm = 1.-1.e-10          ! used to prevent roundoff errors below zero
 315       omsm = 1._r8-2*epsilon(1._r8) ! used to prevent roundoff errors below zero
 316       precmin =  0.1_r8/8.64e4_r8      ! set critical value to 0.1 mm/day in kg/m2/s
 317
 318       adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
 319
 320       ! assume 4 m/s fall speed currently (should be improved)
 321 !      vfall = 4.
 322
 323       ! default (if other sol_facts aren't in call, set all to required sol_fact
 324       sol_facti = sol_fact
 325       sol_factb = sol_fact
 326       sol_factii = sol_fact
 327       sol_factbi = sol_fact
 328
 329       if ( present(sol_facti_in) )  sol_facti = sol_facti_in
 330       if ( present(sol_factii_in) )  sol_factii = sol_factii_in
 331       if ( present(sol_factbi_in) )  sol_factbi = sol_factbi_in
 332
 333       sol_factic  = sol_facti
 334       sol_factiic = sol_factii
 335       if ( present(sol_factic_in ) )  sol_factic  = sol_factic_in
 336       if ( present(sol_factiic_in) )  sol_factiic = sol_factiic_in
 337
 338       ! this section of code is for highly soluble aerosols,
 339       ! the assumption is that within the cloud that
 340       ! all the tracer is in the cloud water
 341       !
 342       ! for both convective and stratiform clouds,
 343       ! the fraction of cloud water converted to precip defines
 344       ! the amount of tracer which is pulled out.
 345       !
 346
 347       do i = 1,pcols
 348          precabs(i) = 0
 349          precabc(i) = 0
 350          scavab(i) = 0
 351          scavabc(i) = 0
 352          tracab(i) = 0
 353          cldmabs(i) = 0
 354          cldmabc(i) = 0
 355
 356        ! Jan.16. Sungsu
 357        ! I added below to compute vertically projected cumulus and stratus fractions from the top to the
 358        ! current model layer by assuming a simple independent maximum overlapping assumption for
 359        ! each cloud.
 360        ! cldovr_cu(i) = 0._r8
 361        ! cldovr_st(i) = 0._r8
 362        ! End by Sungsu
 363
 364       end do
 365
 366       do k = 1,pver
 367          do i = 1,ncol
 368             tc     = t(i,k) - tmelt
 369             weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
 370             weight = 0._r8                                 ! assume no ice
 371
 372             pdog = pdel(i,k)/gravit
 373
 374             ! ****************** Evaporation **************************
 375             ! calculate the fraction of strat precip from above
 376             !                 which evaporates within this layer
 377             fracev(i) = evaps(i,k)*pdel(i,k)/gravit &
 378                      /max(1.e-12_r8,precabs(i))
 379
 380             ! trap to ensure reasonable ratio bounds
 381             fracev(i) = max(0._r8,min(1._r8,fracev(i)))
 382
 383 ! Sungsu : Same as above but convective precipitation part
 384             fracev_cu(i) = evapc(i,k)*pdel(i,k)/gravit/max(1.e-12_r8,precabc(i))
 385             fracev_cu(i) = max(0._r8,min(1._r8,fracev_cu(i)))
 386 ! Sungsu
 387             ! ****************** Convection ***************************
 388             ! now do the convective scavenging
 389
 390             ! set odds proportional to fraction of the grid box that is swept by the
 391             ! precipitation =precabc/rhoh20*(area of sphere projected on plane
 392             !                                /volume of sphere)*deltat
 393             ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
 394             ! unless the fraction of the area that is cloud is less than odds, in which
 395             ! case use the cloud fraction (assumes precabs is in kg/m2/s)
 396             ! is really: precabs*3/4/1000./1e-3*deltat
 397             ! here I use .1 from Balkanski
 398             !
 399             ! use a local rate of convective rain production for incloud scav
 400             !odds=max(min(1._r8, &
 401             !     cmfdqr(i,k)*pdel(i,k)/gravit*0.1_r8*deltat),0._r8)
 402             !++mcb -- change cldc to cldt; change cldt to cldv (9/17/96)
 403             !            srcs1 =  cldt(i,k)*odds*tracer(i,k)*(1.-weight) &
 404             ! srcs1 =  cldv(i,k)*odds*tracer(i,k)*(1.-weight) &
 405             !srcs1 =  cldc(i,k)*odds*tracer(i,k)*(1.-weight) &
 406             !         /deltat
 407
 408             ! fraction of convective cloud water converted to rain
 409             ! Dec.29.2009 : Sungsu multiplied cldc(i,k) to conicw(i,k) below
 410             ! fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,conicw(i,k))
 411             ! fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,cldc(i,k)*conicw(i,k))
 412             ! Sungsu: Below new formula of 'fracp' is necessary since 'conicw' is a LWC/IWC
 413             !         that has already precipitated out, that is, 'conicw' does not contain
 414             !         precipitation at all !
 415               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.
 416             ! Dec.29.2009
 417             ! Note cmfdrq can be negative from evap of rain, so constrain it <-- This is wrong. cmfdqr does not
 418             ! contain evaporation of precipitation.
 419             fracp = max(min(1._r8,fracp),0._r8)
 420             ! remove that amount from within the convective area
 421 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat ! liquid only
 422 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)/deltat             ! any condensation
 423 !           srcs1 = 0.
 424 !           Jan.02.2010. Sungsu : cldt --> cldc below.
 425 #ifdef MODAL_AERO
 426             ! rce 2010/05/01
 427             if ( is_strat_cloudborne ) then
 428                ! only strat in-cloud removal affects strat-cloudborne aerosol
 429                srcs1 = 0._r8
 430             else
 431                tracer_incu = f_act_conv(i,k)*(tracer(i,k)+&
 432                              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)))))))
 433                srcs1 = sol_factic(i,k)*cldc(i,k)*fracp*tracer_incu*(1._r8-weight)/deltat &  ! Liquid
 434                      + sol_factiic    *cldc(i,k)*fracp*tracer_incu*(weight)/deltat          ! Ice
 435             end if
 436 #else
 437             srcs1 = sol_factic(i,k)*cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat &  ! liquid
 438                  +  sol_factiic*cldc(i,k)*fracp*tracer(i,k)*(weight)/deltat      ! ice
 439 #endif
 440
 441
 442             !--mcb
 443
 444             ! scavenge below cloud
 445
 446             !            cldmabc(i) = max(cldc(i,k),cldmabc(i))
 447             !            cldmabc(i) = max(cldt(i,k),cldmabc(i))
 448             ! cldmabc(i) = max(cldv(i,k),cldmabc(i))
 449             ! cldmabc(i) = cldv(i,k)
 450             cldmabc(i) = cldvcu(i,k)
 451
 452             ! Jan. 16. 2010. Sungsu
 453             ! cldmabc(i) = cldmabc(i) * cldovr_cu(i) / max( 0.01_r8, cldovr_cu(i) + cldovr_st(i) )
 454             ! End by Sungsu
 455
 456 #ifdef MODAL_AERO
 457             ! rce 2010/05/01
 458             if ( is_strat_cloudborne ) then
 459                ! only strat in-cloud removal affects strat-cloudborne aerosol
 460                srcs2 = 0._r8
 461             else
 462                tracer_mean = tracer(i,k)*(1._r8-cldc(i,k)*f_act_conv(i,k))-cldc(i,k)*f_act_conv(i,k)*&
 463                              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))))))
 464                tracer_mean = max(0._r8,tracer_mean)
 465                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)
 466                srcs2 = sol_factb *cldmabc(i)*odds*tracer_mean*(1._r8-weight)/deltat & ! Liquid
 467                      + sol_factbi*cldmabc(i)*odds*tracer_mean*(weight)/deltat         ! Ice
 468             end if
 469 #else
 470             odds=max( &
 471                  min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8) &
 472                  *scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
 473             srcs2 = sol_factb*cldmabc(i)*odds*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
 474                  +  sol_factbi*cldmabc(i)*odds*tracer(i,k)*(weight)/deltat    !ice
 475 #endif
 476
 477
 478             !Note that using the temperature-determined weight doesn't make much sense here
 479
 480
 481             srcc = srcs1 + srcs2  ! convective tend by both processes
 482             finc = srcs1/(srcc + 1.e-36_r8) ! fraction in-cloud
 483
 484             ! ****************** Stratiform ***********************
 485             ! now do the stratiform scavenging
 486
 487             ! incloud scavenging
 488
 489 #ifdef MODAL_AERO
 490             ! rce 2010/05/01
 491             if ( is_strat_cloudborne ) then
 492                ! new code for stratiform incloud scav of cloudborne (modal) aerosol
 493                ! >> use the 1st order cw to precip rate calculated in microphysics routine
 494                ! >> cloudborne aerosol resides in cloudy portion of grid cell, so do not apply "cldt" factor
 495              ! fracp = rate1ord_cw2pr_st(i,k)*deltat
 496              ! fracp = max(0._r8,min(1._r8,fracp))
 497                fracp = precs(i,k)*deltat/max(cwat(i,k)+precs(i,k)*deltat,1.e-12_r8) ! Sungsu. Mar.19.2010.
 498                fracp = max(0._r8,min(1._r8,fracp))
 499                srcs1 = sol_facti *fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! Liquid
 500                      + sol_factii*fracp*tracer(i,k)/deltat*(weight)          ! Ice
 501             else
 502                ! strat in-cloud removal only affects strat-cloudborne aerosol
 503                srcs1 = 0._r8
 504             end if
 505             ! end rce 2010/05/01
 506 #else
 507
 508             ! fracp is the fraction of cloud water converted to precip
 509             ! Sungsu modified fracp as the convectiv case.
 510             !        Below new formula by Sungsu of 'fracp' is necessary since 'cwat' is a LWC/IWC
 511             !        that has already precipitated out, that is, 'cwat' does not contain
 512             !        precipitation at all !
 513 !           fracp =  precs(i,k)*deltat/max(cwat(i,k),1.e-12_r8)
 514             fracp =  precs(i,k)*deltat/max(cwat(i,k)+precs(i,k)*deltat,1.e-12_r8) ! Sungsu. Mar.19.2010.
 515             fracp = max(0._r8,min(1._r8,fracp))
 516 !           fracp = 0.     ! for debug
 517
 518             ! assume the corresponding amnt of tracer is removed
 519             !++mcb -- remove cldc; change cldt to cldv
 520             !            srcs1 = (cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat
 521             !            srcs1 = cldv(i,k)*fracp*tracer(i,k)/deltat &
 522 !           srcs1 = cldt(i,k)*fracp*tracer(i,k)/deltat            ! all condensate
 523 !           Jan.02.2010. Sungsu : cldt --> cldt - cldc below.
 524             srcs1 = sol_facti*(cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! liquid
 525                  + sol_factii*(cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat*(weight)       ! ice
 526 #endif
 527
 528
 529             ! below cloud scavenging
 530
 531 !           volume undergoing below cloud scavenging
 532 !           cldmabs(i) = cldv(i,k)   ! precipitating volume
 533 !           cldmabs(i) = cldt(i,k)   ! local cloud volume
 534             cldmabs(i) = cldvst(i,k) ! Stratiform precipitation area at the top interface of current layer
 535
 536             ! Jan. 16. 2010. Sungsu
 537             ! cldmabs(i) = cldmabs(i) * cldovr_st(i) / max( 0.01_r8, cldovr_cu(i) + cldovr_st(i) )
 538             ! End by Sungsu
 539
 540 #ifdef MODAL_AERO
 541             ! rce 2010/05/01
 542             if ( is_strat_cloudborne ) then
 543                ! only strat in-cloud removal affects strat-cloudborne aerosol
 544                srcs2 = 0._r8
 545             else
 546                odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
 547                odds = max(min(1._r8,odds),0._r8)
 548                srcs2 = sol_factb *cldmabs(i)*odds*tracer_mean*(1._r8-weight)/deltat & ! Liquid
 549                      + sol_factbi*cldmabs(i)*odds*tracer_mean*(weight)/deltat         ! Ice
 550             end if
 551 #else
 552             odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
 553             odds = max(min(1._r8,odds),0._r8)
 554             srcs2 =sol_factb*(cldmabs(i)*odds) *tracer(i,k)*(1._r8-weight)/deltat & ! liquid
 555                  + sol_factbi*(cldmabs(i)*odds) *tracer(i,k)*(weight)/deltat       ! ice
 556
 557
 558             !Note that using the temperature-determined weight doesn't make much sense here
 559 #endif
 560
 561
 562             srcs = srcs1 + srcs2             ! total stratiform scavenging
 563             fins=srcs1/(srcs + 1.e-36_r8)    ! fraction taken by incloud processes
 564
 565             ! make sure we dont take out more than is there
 566             ! ratio of amount available to amount removed
 567             rat(i) = tracer(i,k)/max(deltat*(srcc+srcs),1.e-36_r8)
 568             if (rat(i).lt.1._r8) then
 569                srcs = srcs*rat(i)
 570                srcc = srcc*rat(i)
 571             endif
 572             srct(i) = (srcc+srcs)*omsm
 573
 574
 575             ! fraction that is not removed within the cloud
 576             ! (assumed to be interstitial, and subject to convective transport)
 577             fracp = deltat*srct(i)/max(cldmabs(i)*tracer(i,k),1.e-36_r8)  ! amount removed
 578             fracp = max(0._r8,min(1._r8,fracp))
 579             fracis(i,k) = 1._r8 - fracp
 580
 581             ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
 582             ! Sungsu added cumulus contribution in the below 3 blocks
 583
 584             scavt(i,k) = -srct(i) + (fracev(i)*scavab(i)+fracev_cu(i)*scavabc(i))*gravit/pdel(i,k)
 585             iscavt(i,k) = -(srcc*finc + srcs*fins)*omsm
 586
 587             if ( present(icscavt) ) icscavt(i,k) = -(srcc*finc) * omsm
 588             if ( present(isscavt) ) isscavt(i,k) = -(srcs*fins) * omsm
 589             if ( present(bcscavt) ) bcscavt(i,k) = -(srcc * (1-finc)) * omsm +  &
 590                  fracev_cu(i)*scavabc(i)*gravit/pdel(i,k)
 591             if ( present(bsscavt) ) bsscavt(i,k) = -(srcs * (1-fins)) * omsm +  &
 592                  fracev(i)*scavab(i)*gravit/pdel(i,k)
 593
 594             dblchek(i) = tracer(i,k) + deltat*scavt(i,k)
 595
 596             ! now keep track of scavenged mass and precip
 597             scavab(i) = scavab(i)*(1-fracev(i)) + srcs*pdel(i,k)/gravit
 598             precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdel(i,k)/gravit
 599             scavabc(i) = scavabc(i)*(1-fracev_cu(i)) + srcc*pdel(i,k)/gravit
 600             precabc(i) = precabc(i) + (cmfdqr(i,k) - evapc(i,k))*pdel(i,k)/gravit
 601             tracab(i) = tracab(i) + tracer(i,k)*pdel(i,k)/gravit
 602
 603        ! Jan.16.2010. Sungsu
 604        ! Compute convective and stratiform precipitation areas at the base interface
 605        ! of current layer. These are for computing 'below cloud scavenging' in the
 606        ! next layer below.
 607
 608        ! cldovr_cu(i) = max( cldovr_cu(i), cldc(i,k) )
 609        ! cldovr_st(i) = max( cldovr_st(i), max( 0._r8, cldt(i,k) - cldc(i,k) ) )
 610
 611        ! cldovr_cu(i) = max( 0._r8, min ( 1._r8, cldovr_cu(i) ) )
 612        ! cldovr_st(i) = max( 0._r8, min ( 1._r8, cldovr_st(i) ) )
 613
 614        ! End by Sungsu
 615
 616          end do ! End of i = 1, ncol
 617
 618          found = .false.
 619          do i = 1,ncol
 620             if ( dblchek(i) < 0._r8 ) then
 621                found = .true.
 622                exit
 623             end if
 624          end do
 625
 626          if ( found ) then
 627             do i = 1,ncol
 628                if (dblchek(i) .lt. 0._r8) then
 629                   write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
 630                        dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
 631 #ifdef WRF_PORT
 632                   call wrf_message(iulog)
 633 #endif
 634                endif
 635             end do
 636          endif
 637
 638       end do ! End of k = 1, pver
 639
 640    end subroutine wetdepa_v2
 641
 642
 643 !==============================================================================
 644
 645 ! This is the frozen CAM4 version of wetdepa.
 646
 647
 648    subroutine wetdepa_v1( t, p, q, pdel, &
 649                        cldt, cldc, cmfdqr, conicw, precs, conds, &
 650                        evaps, cwat, tracer, deltat, &
 651                        scavt, iscavt, cldv, fracis, sol_fact, ncol, &
 652                        scavcoef,icscavt, isscavt, bcscavt, bsscavt, &
 653                        sol_facti_in, sol_factbi_in, sol_factii_in, &
 654                        sol_factic_in, sol_factiic_in )
 655
 656       !-----------------------------------------------------------------------
 657       ! Purpose:
 658       ! scavenging code for very soluble aerosols
 659       !
 660       ! Author: P. Rasch
 661       ! Modified by T. Bond 3/2003 to track different removals
 662       !-----------------------------------------------------------------------
 663
 664       implicit none
 665
 666       real(r8), intent(in) ::&
 667          t(pcols,pver),        &! temperature
 668          p(pcols,pver),        &! pressure
 669          q(pcols,pver),        &! moisture
 670          pdel(pcols,pver),     &! pressure thikness
 671          cldt(pcols,pver),    &! total cloud fraction
 672          cldc(pcols,pver),     &! convective cloud fraction
 673          cmfdqr(pcols,pver),   &! rate of production of convective precip
 674          conicw(pcols,pver),   &! convective cloud water
 675          cwat(pcols,pver),     &! cloud water amount
 676          precs(pcols,pver),    &! rate of production of stratiform precip
 677          conds(pcols,pver),    &! rate of production of condensate
 678          evaps(pcols,pver),    &! rate of evaporation of precip
 679          cldv(pcols,pver),     &! total cloud fraction
 680          deltat,               &! time step
 681          tracer(pcols,pver)     ! trace species
 682       ! If subroutine is called with just sol_fact:
 683             ! sol_fact is used for both in- and below-cloud scavenging
 684       ! If subroutine is called with optional argument sol_facti_in:
 685             ! sol_fact  is used for below cloud scavenging
 686             ! sol_facti is used for in cloud scavenging
 687          real(r8), intent(in) :: sol_fact ! solubility factor (fraction of aer scavenged below & in, or just below or sol_facti_in is provided)
 688          real(r8), intent(in), optional :: sol_facti_in   ! solubility factor (frac of aerosol scavenged in cloud)
 689          real(r8), intent(in), optional :: sol_factbi_in  ! solubility factor (frac of aerosol scavenged below cloud by ice)
 690          real(r8), intent(in), optional :: sol_factii_in  ! solubility factor (frac of aerosol scavenged in cloud by ice)
 691          real(r8), intent(in), optional :: sol_factic_in(pcols,pver)  ! sol_facti_in for convective clouds
 692          real(r8), intent(in), optional :: sol_factiic_in ! sol_factii_in for convective clouds
 693          real(r8), intent(in) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1 if not MODAL_AERO)
 694
 695       integer, intent(in) :: ncol
 696
 697       real(r8), intent(out) ::&
 698          scavt(pcols,pver),    &! scavenging tend
 699          iscavt(pcols,pver),   &! incloud scavenging tends
 700          fracis(pcols,pver)     ! fraction of species not scavenged
 701
 702       real(r8), intent(out), optional ::    icscavt(pcols,pver)     ! incloud, convective
 703       real(r8), intent(out), optional ::    isscavt(pcols,pver)     ! incloud, stratiform
 704       real(r8), intent(out), optional ::    bcscavt(pcols,pver)     ! below cloud, convective
 705       real(r8), intent(out), optional ::    bsscavt(pcols,pver)     ! below cloud, stratiform
 706
 707       ! local variables
 708
 709       integer i                 ! x index
 710       integer k                 ! z index
 711
 712       real(r8) adjfac               ! factor stolen from cmfmca
 713       real(r8) aqfrac               ! fraction of tracer in aqueous phase
 714       real(r8) cwatc                ! local convective total water amount
 715       real(r8) cwats                ! local stratiform total water amount
 716       real(r8) cwatp                ! local water amount falling from above precip
 717       real(r8) fracev(pcols)        ! fraction of precip from above that is evaporating
 718       real(r8) fracp                ! fraction of cloud water converted to precip
 719       real(r8) gafrac               ! fraction of tracer in gas phasea
 720       real(r8) hconst               ! henry's law solubility constant when equation is expressed
 721                                 ! in terms of mixing ratios
 722       real(r8) mpla                 ! moles / liter H2O entering the layer from above
 723       real(r8) mplb                 ! moles / liter H2O leaving the layer below
 724       real(r8) omsm                 ! 1 - (a small number)
 725       real(r8) part                 !  partial pressure of tracer in atmospheres
 726       real(r8) patm                 ! total pressure in atmospheres
 727       real(r8) pdog                 ! work variable (pdel/gravit)
 728       real(r8) precabc(pcols)       ! conv precip from above (work array)
 729       real(r8) precabs(pcols)       ! strat precip from above (work array)
 730       real(r8) precbl               ! precip falling out of level (work array)
 731       real(r8) precmin              ! minimum convective precip causing scavenging
 732       real(r8) rat(pcols)           ! ratio of amount available to amount removed
 733       real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
 734       real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
 735       real(r8) srcc                 ! tend for convective rain
 736       real(r8) srcs                 ! tend for stratiform rain
 737       real(r8) srct(pcols)          ! work variable
 738       real(r8) tracab(pcols)        ! column integrated tracer amount
 739 !      real(r8) vfall                ! fall speed of precip
 740       real(r8) fins                 ! fraction of rem. rate by strat rain
 741       real(r8) finc                 ! fraction of rem. rate by conv. rain
 742       real(r8) srcs1                ! work variable
 743       real(r8) srcs2                ! work variable
 744       real(r8) tc                   ! temp in celcius
 745       real(r8) weight               ! fraction of condensate which is ice
 746       real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
 747       real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
 748       real(r8) odds                 ! limit on removal rate (proportional to prec)
 749       real(r8) dblchek(pcols)
 750       logical :: found
 751
 752       real(r8) sol_facti,  sol_factb  ! in cloud and below cloud fraction of aerosol scavenged
 753       real(r8) sol_factii, sol_factbi ! in cloud and below cloud fraction of aerosol scavenged by ice
 754       real(r8) sol_factic(pcols,pver)             ! sol_facti for convective clouds
 755       real(r8) sol_factiic            ! sol_factii for convective clouds
 756       ! sol_factic & solfact_iic added for MODAL_AERO.
 757       ! For stratiform cloud, cloudborne aerosol is treated explicitly,
 758       !    and sol_facti is 1.0 for cloudborne, 0.0 for interstitial.
 759       ! For convective cloud, cloudborne aerosol is not treated explicitly,
 760       !    and sol_factic is 1.0 for both cloudborne and interstitial.
 761
 762       ! ------------------------------------------------------------------------
 763 !      omsm = 1.-1.e-10          ! used to prevent roundoff errors below zero
 764       omsm = 1._r8-2*epsilon(1._r8) ! used to prevent roundoff errors below zero
 765       precmin =  0.1_r8/8.64e4_r8      ! set critical value to 0.1 mm/day in kg/m2/s
 766
 767       adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
 768
 769       ! assume 4 m/s fall speed currently (should be improved)
 770 !      vfall = 4.
 771
 772       ! default (if other sol_facts aren't in call, set all to required sol_fact
 773       sol_facti = sol_fact
 774       sol_factb = sol_fact
 775       sol_factii = sol_fact
 776       sol_factbi = sol_fact
 777
 778       if ( present(sol_facti_in) )  sol_facti = sol_facti_in
 779       if ( present(sol_factii_in) )  sol_factii = sol_factii_in
 780       if ( present(sol_factbi_in) )  sol_factbi = sol_factbi_in
 781
 782       sol_factic  = sol_facti
 783       sol_factiic = sol_factii
 784       if ( present(sol_factic_in ) )  sol_factic  = sol_factic_in
 785       if ( present(sol_factiic_in) )  sol_factiic = sol_factiic_in
 786
 787       ! this section of code is for highly soluble aerosols,
 788       ! the assumption is that within the cloud that
 789       ! all the tracer is in the cloud water
 790       !
 791       ! for both convective and stratiform clouds,
 792       ! the fraction of cloud water converted to precip defines
 793       ! the amount of tracer which is pulled out.
 794       !
 795
 796       do i = 1,pcols
 797          precabs(i) = 0
 798          precabc(i) = 0
 799          scavab(i) = 0
 800          scavabc(i) = 0
 801          tracab(i) = 0
 802          cldmabs(i) = 0
 803          cldmabc(i) = 0
 804       end do
 805
 806       do k = 1,pver
 807          do i = 1,ncol
 808             tc     = t(i,k) - tmelt
 809             weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
 810             weight = 0._r8                                 ! assume no ice
 811
 812             pdog = pdel(i,k)/gravit
 813
 814             ! ****************** Evaporation **************************
 815             ! calculate the fraction of strat precip from above
 816             !                 which evaporates within this layer
 817             fracev(i) = evaps(i,k)*pdel(i,k)/gravit &
 818                      /max(1.e-12_r8,precabs(i))
 819
 820             ! trap to ensure reasonable ratio bounds
 821             fracev(i) = max(0._r8,min(1._r8,fracev(i)))
 822
 823             ! ****************** Convection ***************************
 824             ! now do the convective scavenging
 825
 826             ! set odds proportional to fraction of the grid box that is swept by the
 827             ! precipitation =precabc/rhoh20*(area of sphere projected on plane
 828             !                                /volume of sphere)*deltat
 829             ! assume the radius of a raindrop is 1 e-3 m from Rogers and Yau,
 830             ! unless the fraction of the area that is cloud is less than odds, in which
 831             ! case use the cloud fraction (assumes precabs is in kg/m2/s)
 832             ! is really: precabs*3/4/1000./1e-3*deltat
 833             ! here I use .1 from Balkanski
 834             !
 835             ! use a local rate of convective rain production for incloud scav
 836             !odds=max(min(1._r8, &
 837             !     cmfdqr(i,k)*pdel(i,k)/gravit*0.1_r8*deltat),0._r8)
 838             !++mcb -- change cldc to cldt; change cldt to cldv (9/17/96)
 839             !            srcs1 =  cldt(i,k)*odds*tracer(i,k)*(1.-weight) &
 840             ! srcs1 =  cldv(i,k)*odds*tracer(i,k)*(1.-weight) &
 841             !srcs1 =  cldc(i,k)*odds*tracer(i,k)*(1.-weight) &
 842             !         /deltat
 843
 844             ! fraction of convective cloud water converted to rain
 845             fracp = cmfdqr(i,k)*deltat/max(1.e-8_r8,conicw(i,k))
 846             ! note cmfdrq can be negative from evap of rain, so constrain it
 847             fracp = max(min(1._r8,fracp),0._r8)
 848             ! remove that amount from within the convective area
 849 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat ! liquid only
 850 !           srcs1 = cldc(i,k)*fracp*tracer(i,k)/deltat             ! any condensation
 851 !           srcs1 = 0.
 852             srcs1 = sol_factic(i,k)*cldt(i,k)*fracp*tracer(i,k)*(1._r8-weight)/deltat &  ! liquid
 853                  +  sol_factiic*cldt(i,k)*fracp*tracer(i,k)*(weight)/deltat      ! ice
 854
 855
 856             !--mcb
 857
 858             ! scavenge below cloud
 859
 860             !            cldmabc(i) = max(cldc(i,k),cldmabc(i))
 861             !            cldmabc(i) = max(cldt(i,k),cldmabc(i))
 862             cldmabc(i) = max(cldv(i,k),cldmabc(i))
 863             cldmabc(i) = cldv(i,k)
 864
 865             odds=max( &
 866                  min(1._r8,precabc(i)/max(cldmabc(i),1.e-5_r8) &
 867                  *scavcoef(i,k)*deltat),0._r8) ! Dana and Hales coefficient (/mm)
 868             srcs2 = sol_factb*cldmabc(i)*odds*tracer(i,k)*(1._r8-weight)/deltat & ! liquid
 869                  +  sol_factbi*cldmabc(i)*odds*tracer(i,k)*(weight)/deltat    !ice
 870             !Note that using the temperature-determined weight doesn't make much sense here
 871
 872
 873             srcc = srcs1 + srcs2  ! convective tend by both processes
 874             finc = srcs1/(srcc + 1.e-36_r8) ! fraction in-cloud
 875
 876             ! ****************** Stratiform ***********************
 877             ! now do the stratiform scavenging
 878
 879             ! incloud scavenging
 880
 881             ! fracp is the fraction of cloud water converted to precip
 882             fracp =  precs(i,k)*deltat/max(cwat(i,k),1.e-12_r8)
 883             fracp = max(0._r8,min(1._r8,fracp))
 884 !           fracp = 0.     ! for debug
 885
 886             ! assume the corresponding amnt of tracer is removed
 887             !++mcb -- remove cldc; change cldt to cldv
 888             !            srcs1 = (cldt(i,k)-cldc(i,k))*fracp*tracer(i,k)/deltat
 889             !            srcs1 = cldv(i,k)*fracp*tracer(i,k)/deltat &
 890 !           srcs1 = cldt(i,k)*fracp*tracer(i,k)/deltat            ! all condensate
 891             srcs1 = sol_facti*cldt(i,k)*fracp*tracer(i,k)/deltat*(1._r8-weight) &  ! liquid
 892                  + sol_factii*cldt(i,k)*fracp*tracer(i,k)/deltat*(weight)       ! ice
 893
 894
 895             ! below cloud scavenging
 896
 897 !           volume undergoing below cloud scavenging
 898             cldmabs(i) = cldv(i,k)   ! precipitating volume
 899 !           cldmabs(i) = cldt(i,k)   ! local cloud volume
 900
 901             odds = precabs(i)/max(cldmabs(i),1.e-5_r8)*scavcoef(i,k)*deltat
 902             odds = max(min(1._r8,odds),0._r8)
 903             srcs2 =sol_factb*(cldmabs(i)*odds) *tracer(i,k)*(1._r8-weight)/deltat & ! liquid
 904                  + sol_factbi*(cldmabs(i)*odds) *tracer(i,k)*(weight)/deltat       ! ice
 905             !Note that using the temperature-determined weight doesn't make much sense here
 906
 907
 908             srcs = srcs1 + srcs2             ! total stratiform scavenging
 909             fins=srcs1/(srcs + 1.e-36_r8)    ! fraction taken by incloud processes
 910
 911             ! make sure we dont take out more than is there
 912             ! ratio of amount available to amount removed
 913             rat(i) = tracer(i,k)/max(deltat*(srcc+srcs),1.e-36_r8)
 914             if (rat(i).lt.1._r8) then
 915                srcs = srcs*rat(i)
 916                srcc = srcc*rat(i)
 917             endif
 918             srct(i) = (srcc+srcs)*omsm
 919
 920
 921             ! fraction that is not removed within the cloud
 922             ! (assumed to be interstitial, and subject to convective transport)
 923             fracp = deltat*srct(i)/max(cldmabs(i)*tracer(i,k),1.e-36_r8)  ! amount removed
 924             fracp = max(0._r8,min(1._r8,fracp))
 925             fracis(i,k) = 1._r8 - fracp
 926
 927             ! tend is all tracer removed by scavenging, plus all re-appearing from evaporation above
 928             scavt(i,k) = -srct(i) + fracev(i)*scavab(i)*gravit/pdel(i,k)
 929             iscavt(i,k) = -(srcc*finc + srcs*fins)*omsm
 930
 931             if ( present(icscavt) ) icscavt(i,k) = -(srcc*finc) * omsm
 932             if ( present(isscavt) ) isscavt(i,k) = -(srcs*fins) * omsm
 933             if ( present(bcscavt) ) bcscavt(i,k) = -(srcc * (1-finc)) * omsm
 934             if ( present(bsscavt) ) bsscavt(i,k) = -(srcs * (1-fins)) * omsm +  &
 935                  fracev(i)*scavab(i)*gravit/pdel(i,k)
 936
 937             dblchek(i) = tracer(i,k) + deltat*scavt(i,k)
 938
 939             ! now keep track of scavenged mass and precip
 940             scavab(i) = scavab(i)*(1-fracev(i)) + srcs*pdel(i,k)/gravit
 941             precabs(i) = precabs(i) + (precs(i,k) - evaps(i,k))*pdel(i,k)/gravit
 942             scavabc(i) = scavabc(i) + srcc*pdel(i,k)/gravit
 943             precabc(i) = precabc(i) + (cmfdqr(i,k))*pdel(i,k)/gravit
 944             tracab(i) = tracab(i) + tracer(i,k)*pdel(i,k)/gravit
 945
 946          end do
 947
 948          found = .false.
 949          do i = 1,ncol
 950             if ( dblchek(i) < 0._r8 ) then
 951                found = .true.
 952                exit
 953             end if
 954          end do
 955
 956          if ( found ) then
 957             do i = 1,ncol
 958                if (dblchek(i) .lt. 0._r8) then
 959                   write(iulog,*) ' wetdapa: negative value ', i, k, tracer(i,k), &
 960                        dblchek(i), scavt(i,k), srct(i), rat(i), fracev(i)
 961 #ifdef WRF_PORT
 962                   call wrf_message(iulog)
 963 #endif
 964                endif
 965             end do
 966          endif
 967
 968       end do
 969
 970    end subroutine wetdepa_v1
 971 #ifndef WRF_PORT
 972 !==============================================================================
 973
 974 ! wetdepg is currently being used for both CAM4 and CAM5 by making use of the
 975 ! cam_physpkg_is method.
 976
 977    subroutine wetdepg( t, p, q, pdel, &
 978                        cldt, cldc, cmfdqr, evapc, precs, evaps, &
 979                        rain, cwat, tracer, deltat, molwt, &
 980                        solconst, scavt, iscavt, cldv, icwmr1, &
 981                        icwmr2, fracis, ncol )
 982
 983       !-----------------------------------------------------------------------
 984       ! Purpose:
 985       ! scavenging of gas phase constituents by henry's law
 986       !
 987       ! Author: P. Rasch
 988       !-----------------------------------------------------------------------
 989
 990       real(r8), intent(in) ::&
 991          t(pcols,pver),        &! temperature
 992          p(pcols,pver),        &! pressure
 993          q(pcols,pver),        &! moisture
 994          pdel(pcols,pver),     &! pressure thikness
 995          cldt(pcols,pver),     &! total cloud fraction
 996          cldc(pcols,pver),     &! convective cloud fraction
 997          cmfdqr(pcols,pver),   &! rate of production of convective precip
 998          rain (pcols,pver),    &! total rainwater mixing ratio
 999          cwat(pcols,pver),     &! cloud water amount
1000          precs(pcols,pver),    &! rate of production of stratiform precip
1001          evaps(pcols,pver),    &! rate of evaporation of precip
1002 ! Sungsu
1003          evapc(pcols,pver),    &! Rate of evaporation of convective precipitation
1004 ! Sungsu
1005          cldv(pcols,pver),     &! estimate of local volume occupied by clouds
1006          icwmr1 (pcols,pver),  &! in cloud water mixing ration for zhang scheme
1007          icwmr2 (pcols,pver),  &! in cloud water mixing ration for hack  scheme
1008          deltat,               &! time step
1009          tracer(pcols,pver),   &! trace species
1010          molwt                  ! molecular weights
1011
1012       integer, intent(in) :: ncol
1013
1014       real(r8) &
1015          solconst(pcols,pver)   ! Henry's law coefficient
1016
1017       real(r8), intent(out) ::&
1018          scavt(pcols,pver),    &! scavenging tend
1019          iscavt(pcols,pver),   &! incloud scavenging tends
1020          fracis(pcols, pver)    ! fraction of constituent that is insoluble
1021
1022       ! local variables
1023
1024       integer i                 ! x index
1025       integer k                 ! z index
1026
1027       real(r8) adjfac               ! factor stolen from cmfmca
1028       real(r8) aqfrac               ! fraction of tracer in aqueous phase
1029       real(r8) cwatc                ! local convective total water amount
1030       real(r8) cwats                ! local stratiform total water amount
1031       real(r8) cwatl                ! local cloud liq water amount
1032       real(r8) cwatp                ! local water amount falling from above precip
1033       real(r8) cwatpl               ! local water amount falling from above precip (liq)
1034       real(r8) cwatt                ! local sum of strat + conv total water amount
1035       real(r8) cwatti               ! cwatt/cldv = cloudy grid volume mixing ratio
1036       real(r8) fracev               ! fraction of precip from above that is evaporating
1037       real(r8) fracp                ! fraction of cloud water converted to precip
1038       real(r8) gafrac               ! fraction of tracer in gas phasea
1039       real(r8) hconst               ! henry's law solubility constant when equation is expressed
1040                                 ! in terms of mixing ratios
1041       real(r8) mpla                 ! moles / liter H2O entering the layer from above
1042       real(r8) mplb                 ! moles / liter H2O leaving the layer below
1043       real(r8) omsm                 ! 1 - (a small number)
1044       real(r8) part                 !  partial pressure of tracer in atmospheres
1045       real(r8) patm                 ! total pressure in atmospheres
1046       real(r8) pdog                 ! work variable (pdel/gravit)
1047       real(r8) precab(pcols)        ! precip from above (work array)
1048       real(r8) precbl               ! precip work variable
1049       real(r8) precxx               ! precip work variable
1050       real(r8) precxx2               !
1051       real(r8) precic               ! precip work variable
1052       real(r8) rat                  ! ratio of amount available to amount removed
1053       real(r8) scavab(pcols)        ! scavenged tracer flux from above (work array)
1054       real(r8) scavabc(pcols)       ! scavenged tracer flux from above (work array)
1055       !      real(r8) vfall                ! fall speed of precip
1056       real(r8) scavmax              ! an estimate of the max tracer avail for removal
1057       real(r8) scavbl               ! flux removed at bottom of layer
1058       real(r8) fins                 ! in cloud fraction removed by strat rain
1059       real(r8) finc                 ! in cloud fraction removed by conv rain
1060       real(r8) rate                 ! max removal rate estimate
1061       real(r8) scavlimt             ! limiting value 1
1062       real(r8) scavt1               ! limiting value 2
1063       real(r8) scavin               ! scavenging by incloud processes
1064       real(r8) scavbc               ! scavenging by below cloud processes
1065       real(r8) tc
1066       real(r8) weight               ! ice fraction
1067       real(r8) wtpl                 ! work variable
1068       real(r8) cldmabs(pcols)       ! maximum cloud at or above this level
1069       real(r8) cldmabc(pcols)       ! maximum cloud at or above this level
1070       !-----------------------------------------------------------
1071
1072       omsm = 1._r8-2*epsilon(1._r8)   ! used to prevent roundoff errors below zero
1073
1074       adjfac = deltat/(max(deltat,cmftau)) ! adjustment factor from hack scheme
1075
1076       ! assume 4 m/s fall speed currently (should be improved)
1077       !      vfall = 4.
1078
1079       ! zero accumulators
1080       do i = 1,pcols
1081          precab(i) = 1.e-36_r8
1082          scavab(i) = 0._r8
1083          cldmabs(i) = 0._r8
1084       end do
1085
1086       do k = 1,pver
1087          do i = 1,ncol
1088
1089             tc     = t(i,k) - tmelt
1090             weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
1091
1092             cldmabs(i) = max(cldmabs(i),cldt(i,k))
1093
1094             ! partitioning coefs for gas and aqueous phase
1095             !              take as a cloud water amount, the sum of the stratiform amount
1096             !              plus the convective rain water amount
1097
1098             ! convective amnt is just the local precip rate from the hack scheme
1099             !              since there is no storage of water, this ignores that falling from above
1100             !            cwatc = cmfdqr(i,k)*deltat/adjfac
1101             !++mcb -- test cwatc
1102             cwatc = (icwmr1(i,k) + icwmr2(i,k)) * (1._r8-weight)
1103             !--mcb
1104
1105             ! strat cloud water amount and also ignore the part falling from above
1106             cwats = cwat(i,k)
1107
1108             ! cloud water as liq
1109             !++mcb -- add cwatc later (in cwatti)
1110             !            cwatl = (1.-weight)*(cwatc+cwats)
1111             cwatl = (1._r8-weight)*cwats
1112             ! cloud water as ice
1113             !*not used        cwati = weight*(cwatc+cwats)
1114
1115             ! total suspended condensate as liquid
1116             cwatt = cwatl + rain(i,k)
1117
1118             ! incloud version
1119             !++mcb -- add cwatc here
1120             cwatti = cwatt/max(cldv(i,k), 0.00001_r8) + cwatc
1121
1122             ! partitioning terms
1123             patm = p(i,k)/1.013e5_r8 ! pressure in atmospheres
1124             hconst = molwta*patm*solconst(i,k)*cwatti/rhoh2o
1125             aqfrac = hconst/(1._r8+hconst)
1126             gafrac = 1/(1._r8+hconst)
1127             fracis(i,k) = gafrac
1128
1129
1130             ! partial pressure of the tracer in the gridbox in atmospheres
1131             part = patm*gafrac*tracer(i,k)*molwta/molwt
1132
1133             ! use henrys law to give moles tracer /liter of water
1134             ! in this volume
1135             ! then convert to kg tracer /liter of water (kg tracer / kg water)
1136             mplb = solconst(i,k)*part*molwt/1000._r8
1137
1138
1139             pdog = pdel(i,k)/gravit
1140
1141             ! this part of precip will be carried downward but at a new molarity of mpl
1142             precic = pdog*(precs(i,k) + cmfdqr(i,k))
1143
1144             ! we cant take out more than entered, plus that available in the cloud
1145             !                  scavmax = scavab(i)+tracer(i,k)*cldt(i,k)/deltat*pdog
1146             scavmax = scavab(i)+tracer(i,k)*cldv(i,k)/deltat*pdog
1147
1148             ! flux of tracer by incloud processes
1149             scavin = precic*(1._r8-weight)*mplb
1150
1151             ! fraction of precip which entered above that leaves below
1152             if (cam_physpkg_is('cam5')) then
1153                ! Sungsu added evaporation of convective precipitation below.
1154                precxx = precab(i)-pdog*(evaps(i,k)+evapc(i,k))
1155             else
1156                precxx = precab(i)-pdog*evaps(i,k)
1157             end if
1158             precxx = max (precxx,0.0_r8)
1159
1160             ! flux of tracer by below cloud processes
1161             !++mcb -- removed wtpl because it is now not assigned and previously
1162             !          when it was assigned it was unnecessary:  if(tc.gt.0)wtpl=1
1163             if (tc.gt.0) then
1164                !               scavbc = precxx*wtpl*mplb ! if liquid
1165                scavbc = precxx*mplb ! if liquid
1166             else
1167                precxx2=max(precxx,1.e-36_r8)
1168                scavbc = scavab(i)*precxx2/(precab(i)) ! if ice
1169             endif
1170
1171             scavbl = min(scavbc + scavin, scavmax)
1172
1173             ! first guess assuming that henries law works
1174             scavt1 = (scavab(i)-scavbl)/pdog*omsm
1175
1176             ! pjr this should not be required, but we put it in to make sure we cant remove too much
1177             ! remember, scavt1 is generally negative (indicating removal)
1178             scavt1 = max(scavt1,-tracer(i,k)*cldv(i,k)/deltat)
1179
1180             !++mcb -- remove this limitation for gas species
1181             !c use the dana and hales or balkanski limit on scavenging
1182             !c            rate = precab(i)*0.1
1183             !            rate = (precic + precxx)*0.1
1184             !            scavlimt = -tracer(i,k)*cldv(i,k)
1185             !     $           *rate/(1.+rate*deltat)
1186
1187             !            scavt(i,k) = max(scavt1, scavlimt)
1188
1189             ! instead just set scavt to scavt1
1190             scavt(i,k) = scavt1
1191             !--mcb
1192
1193             ! now update the amount leaving the layer
1194             scavbl = scavab(i) - scavt(i,k)*pdog
1195
1196             ! in cloud amount is that formed locally over the total flux out bottom
1197             fins = scavin/(scavin + scavbc + 1.e-36_r8)
1198             iscavt(i,k) = scavt(i,k)*fins
1199
1200             scavab(i) = scavbl
1201             precab(i) = max(precxx + precic,1.e-36_r8)
1202
1203
1204
1205          end do
1206       end do
1207
1208    end subroutine wetdepg
1209 #endif
1210 !##############################################################################
1211
1212 end module wetdep