LCOV - code coverage report
Current view: | top level - chemistry/aerosol - wetdep.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 18 | 322 | 5.6 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 7 | 14.3 %

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

Generated by: LCOV version 1.14