LCOV - code coverage report
Current view: | top level - chemistry/aerosol - aero_convproc.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 643 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 8 | 0.0 %

              Line data    Source code
           1             : module aero_convproc
           2             : !---------------------------------------------------------------------------------
           3             : ! Purpose:
           4             : !
           5             : ! CAM interface to aerosol/trace-gas convective cloud processing scheme
           6             : !
           7             : ! currently these routines assume stratiform and convective clouds only interact
           8             : ! through the detrainment of convective cloudborne material into stratiform clouds
           9             : !
          10             : ! thus the stratiform-cloudborne aerosols (in the qqcw array) are not processed
          11             : ! by the convective up/downdrafts, but are affected by the detrainment
          12             : !
          13             : ! Author: R. C. Easter
          14             : !
          15             : !---------------------------------------------------------------------------------
          16             :
          17             : use shr_kind_mod,    only: r8=>shr_kind_r8
          18             : use shr_kind_mod,    only: shr_kind_cs
          19             :
          20             : use spmd_utils,      only: masterproc
          21             : use physconst,       only: gravit, rair
          22             : use ppgrid,          only: pver, pcols, pverp
          23             : use constituents,    only: pcnst, cnst_get_ind
          24             : use constituents,    only: cnst_species_class, cnst_spec_class_aerosol
          25             : use phys_control,    only: phys_getopts
          26             :
          27             : use physics_types,   only: physics_state, physics_ptend
          28             : use physics_buffer,  only: physics_buffer_desc, pbuf_get_index, pbuf_get_field
          29             : use time_manager,    only: get_nstep
          30             : use cam_history,     only: outfld, addfld, add_default, horiz_only
          31             : use cam_logfile,     only: iulog
          32             : use cam_abortutils,  only: endrun
          33             :
          34             : use aerosol_properties_mod, only: aerosol_properties
          35             : use aerosol_state_mod, only: aerosol_state, ptr2d_t
          36             :
          37             : implicit none
          38             : private
          39             :
          40             : public :: aero_convproc_readnl
          41             : public :: aero_convproc_init
          42             : public :: aero_convproc_intr
          43             :
          44             : ! namelist options
          45             : ! NOTE: These are the defaults for CAM6.
          46             : logical, protected, public :: deepconv_wetdep_history = .true.
          47             : logical, protected, public :: convproc_do_deep = .true.
          48             : ! NOTE: These are the defaults for the Eaton/Wang parameterization.
          49             : logical, protected, public :: convproc_do_evaprain_atonce = .false.
          50             : real(r8), protected, public    :: convproc_pom_spechygro = -1._r8
          51             : real(r8), protected, public    :: convproc_wup_max       = 4.0_r8
          52             :
          53             : logical, parameter :: use_cwaer_for_activate_maxsat = .false.
          54             : logical, parameter :: apply_convproc_tend_to_ptend = .true.
          55             :
          56             : real(r8) :: hund_ovr_g ! = 100.0_r8/gravit
          57             : !  used with zm_conv mass fluxes and delta-p
          58             : !     for mu = [mbar/s],   mu*hund_ovr_g = [kg/m2/s]
          59             : !     for dp = [mbar] and q = [kg/kg],   q*dp*hund_ovr_g = [kg/m2]
          60             :
          61             : !  method1_activate_nlayers = number of layers (including cloud base) where activation is applied
          62             : integer, parameter  :: method1_activate_nlayers = 2
          63             : !  method2_activate_smaxmax = the uniform or peak supersat value (as 0-1 fraction = percent*0.01)
          64             : real(r8), parameter :: method2_activate_smaxmax = 0.003_r8
          65             :
          66             : !  method_reduce_actfrac = 1 -- multiply activation fractions by factor_reduce_actfrac
          67             : !                               (this works ok with convproc_method_activate = 1 but not for ... = 2)
          68             : !                        = 2 -- do 2 iterations to get an overall reduction by factor_reduce_actfrac
          69             : !                               (this works ok with convproc_method_activate = 1 or 2)
          70             : !                        = other -- do nothing involving reduce_actfrac
          71             : integer, parameter  :: method_reduce_actfrac = 0
          72             : real(r8), parameter :: factor_reduce_actfrac = 0.5_r8
          73             :
          74             : !  convproc_method_activate - 1=apply abdulrazzak-ghan to entrained aerosols for lowest nlayers
          75             : !                             2=do secondary activation with prescribed supersat
          76             : integer, parameter :: convproc_method_activate = 2
          77             :
          78             : logical :: convproc_do_aer
          79             :
          80             : ! physics buffer indices
          81             : integer :: fracis_idx         = 0
          82             :
          83             : integer :: rprddp_idx         = 0
          84             : integer :: rprdsh_idx         = 0
          85             : integer :: nevapr_shcu_idx    = 0
          86             : integer :: nevapr_dpcu_idx    = 0
          87             :
          88             : integer :: icwmrdp_idx        = 0
          89             : integer :: icwmrsh_idx        = 0
          90             : integer :: sh_frac_idx        = 0
          91             : integer :: dp_frac_idx        = 0
          92             :
          93             : integer :: zm_eu_idx          = 0
          94             : integer :: zm_du_idx          = 0
          95             : integer :: zm_ed_idx          = 0
          96             : integer :: zm_dp_idx          = 0
          97             : integer :: zm_jt_idx          = 0
          98             : integer :: zm_maxg_idx        = 0
          99             : integer :: zm_ideep_idx       = 0
         100             :
         101             : integer :: cmfmc_sh_idx       = 0
         102             : integer :: sh_e_ed_ratio_idx  = 0
         103             :
         104             : integer :: istat
         105             :
         106             : integer :: nbins = 0
         107             : integer :: ncnstaer = 0
         108             :
         109             : integer, allocatable :: aer_cnst_ndx(:)
         110             :
         111             : character(len=32), allocatable :: cnst_name_extd(:,:) ! (2,ncnstaer)
         112             :
         113             : contains
         114             :
         115             : !=========================================================================================
         116           0 : subroutine aero_convproc_readnl(nlfile)
         117             :
         118             :   use namelist_utils, only: find_group_name
         119             :   use spmd_utils,     only: mpicom, masterprocid, mpi_real8, mpi_logical
         120             :
         121             :   character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         122             :
         123             :   ! Local variables
         124             :   integer :: unitn, ierr
         125             :   character(len=*), parameter :: subname = 'aero_convproc_readnl'
         126             :
         127             :   namelist /aerosol_convproc_opts/ deepconv_wetdep_history, convproc_do_deep, &
         128             :        convproc_do_evaprain_atonce, convproc_pom_spechygro, convproc_wup_max
         129             :
         130             :   ! Read namelist
         131           0 :   if (masterproc) then
         132           0 :      open( newunit=unitn, file=trim(nlfile), status='old' )
         133           0 :      call find_group_name(unitn, 'aerosol_convproc_opts', status=ierr)
         134           0 :      if (ierr == 0) then
         135           0 :         read(unitn, aerosol_convproc_opts, iostat=ierr)
         136           0 :         if (ierr /= 0) then
         137           0 :            call endrun(subname // ':: ERROR reading namelist')
         138             :         end if
         139             :      end if
         140           0 :      close(unitn)
         141             :   end if
         142             :
         143             :   ! Broadcast namelist variables
         144           0 :   call mpi_bcast( deepconv_wetdep_history,  1, mpi_logical, masterprocid, mpicom, ierr)
         145           0 :   call mpi_bcast( convproc_do_deep,  1, mpi_logical, masterprocid, mpicom, ierr)
         146           0 :   call mpi_bcast( convproc_do_evaprain_atonce,  1, mpi_logical, masterprocid, mpicom, ierr)
         147           0 :   call mpi_bcast( convproc_pom_spechygro,  1, mpi_real8, masterprocid, mpicom, ierr)
         148           0 :   call mpi_bcast( convproc_wup_max,  1, mpi_real8, masterprocid, mpicom, ierr)
         149             :
         150           0 :   if (masterproc) then
         151           0 :      write(iulog,*) subname//': deepconv_wetdep_history = ',deepconv_wetdep_history
         152           0 :      write(iulog,*) subname//': convproc_do_deep = ',convproc_do_deep
         153           0 :      write(iulog,*) subname//': convproc_do_evaprain_atonce = ',convproc_do_evaprain_atonce
         154           0 :      write(iulog,*) subname//': convproc_pom_spechygro = ',convproc_pom_spechygro
         155           0 :      write(iulog,*) subname//': convproc_wup_max = ', convproc_wup_max
         156             :   end if
         157             :
         158           0 : end subroutine aero_convproc_readnl
         159             :
         160             : !=========================================================================================
         161             :
         162           0 : subroutine aero_convproc_init(aero_props)
         163             :
         164             :   class(aerosol_properties), intent(in) :: aero_props
         165             :
         166             :    integer :: m, mm, l, ndx, astat
         167             :    integer :: npass_calc_updraft
         168             :    logical :: history_aerosol
         169             :    character(len=32) :: name_a, name_c
         170             :
         171             :    character(len=*), parameter :: prefix = 'aero_convproc_init: '
         172             :
         173           0 :    hund_ovr_g = 100.0_r8/gravit
         174             :    !  used with zm_conv mass fluxes and delta-p
         175             :    !     for mu = [mbar/s],   mu*hund_ovr_g = [kg/m2/s]
         176             :    !     for dp = [mbar] and q = [kg/kg],   q*dp*hund_ovr_g = [kg/m2]
         177             :
         178           0 :    nbins = aero_props%nbins()
         179           0 :    ncnstaer = aero_props%ncnst_tot()
         180             :
         181           0 :    allocate(aer_cnst_ndx(ncnstaer),stat=astat)
         182           0 :    if (astat/=0) then
         183           0 :       call endrun(prefix//'aer_cnst_ndx allocation error')
         184             :    end if
         185           0 :    allocate(cnst_name_extd(2,ncnstaer),stat=astat)
         186           0 :    if (astat/=0) then
         187           0 :       call endrun(prefix//'cnst_name_extd allocation error')
         188             :    end if
         189             :
         190           0 :    aer_cnst_ndx(:) = -1
         191             :
         192           0 :    do m = 1, aero_props%nbins()
         193           0 :       do l = 0, aero_props%nmasses(m)
         194           0 :          mm = aero_props%indexer(m,l)
         195           0 :          if (l==0) then
         196           0 :             call aero_props%num_names(m, name_a, name_c)
         197             :          else
         198           0 :             call aero_props%mmr_names(m,l, name_a, name_c)
         199             :          endif
         200           0 :          cnst_name_extd(1,mm) = name_a
         201           0 :          cnst_name_extd(2,mm) = name_c
         202             :
         203           0 :          call cnst_get_ind(trim(name_a), ndx, abort=.false.)
         204           0 :          aer_cnst_ndx(mm) = ndx
         205             :       end do
         206             :    end do
         207             :
         208             :    call phys_getopts( history_aerosol_out=history_aerosol, &
         209           0 :         convproc_do_aer_out = convproc_do_aer )
         210             :
         211             :    call addfld('DP_MFUP_MAX', horiz_only, 'A', 'kg/m2', &
         212           0 :                'Deep conv. column-max updraft mass flux' )
         213             :    call addfld('DP_WCLDBASE', horiz_only, 'A', 'm/s', &
         214           0 :                'Deep conv. cloudbase vertical velocity' )
         215             :    call addfld('DP_KCLDBASE', horiz_only, 'A', '1', &
         216           0 :         'Deep conv. cloudbase level index' )
         217             :
         218             :    ! output wet deposition fields to history
         219             :    !    I = in-cloud removal;     E = precip-evap resuspension
         220             :    !    C = convective (total);   D = deep convective
         221             :    ! note that the precip-evap resuspension includes that resulting from
         222             :    !    below-cloud removal, calculated in mz_aero_wet_intr
         223           0 :    if (convproc_do_aer .and. apply_convproc_tend_to_ptend ) then
         224             :
         225           0 :       do m = 1, aero_props%nbins()
         226           0 :          do l = 0, aero_props%nmasses(m)
         227           0 :             mm = aero_props%indexer(m,l)
         228             :
         229           0 :             ndx = aer_cnst_ndx(mm)
         230             :
         231           0 :             if ( deepconv_wetdep_history ) then
         232           0 :                call addfld (trim(cnst_name_extd(1,mm))//'SFSID', &
         233           0 :                     horiz_only,  'A','kg/m2/s','Wet deposition flux (incloud, deep convective) at surface')
         234           0 :                call addfld (trim(cnst_name_extd(1,mm))//'SFSED', &
         235           0 :                     horiz_only,  'A','kg/m2/s','Wet deposition flux (precip evap, deep convective) at surface')
         236           0 :                if (history_aerosol) then
         237           0 :                   call add_default(trim(cnst_name_extd(1,mm))//'SFSID', 1, ' ')
         238           0 :                   call add_default(trim(cnst_name_extd(1,mm))//'SFSED', 1, ' ')
         239             :                end if
         240             :             end if
         241             :
         242             :          end do
         243             :       end do
         244             :    end if
         245             :
         246           0 :    if ( history_aerosol .and. convproc_do_aer ) then
         247           0 :       call add_default( 'DP_MFUP_MAX', 1, ' ' )
         248           0 :       call add_default( 'DP_WCLDBASE', 1, ' ' )
         249           0 :       call add_default( 'DP_KCLDBASE', 1, ' ' )
         250             :    end if
         251             :
         252           0 :    fracis_idx      = pbuf_get_index('FRACIS')
         253             :
         254           0 :    rprddp_idx      = pbuf_get_index('RPRDDP')
         255           0 :    rprdsh_idx      = pbuf_get_index('RPRDSH')
         256           0 :    nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
         257           0 :    nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
         258             :
         259           0 :    icwmrdp_idx     = pbuf_get_index('ICWMRDP')
         260           0 :    icwmrsh_idx     = pbuf_get_index('ICWMRSH')
         261           0 :    dp_frac_idx     = pbuf_get_index('DP_FRAC')
         262           0 :    sh_frac_idx     = pbuf_get_index('SH_FRAC')
         263             :
         264           0 :    zm_eu_idx       = pbuf_get_index('ZM_EU')
         265           0 :    zm_du_idx       = pbuf_get_index('ZM_DU')
         266           0 :    zm_ed_idx       = pbuf_get_index('ZM_ED')
         267           0 :    zm_dp_idx       = pbuf_get_index('ZM_DP')
         268           0 :    zm_jt_idx       = pbuf_get_index('ZM_JT')
         269           0 :    zm_maxg_idx     = pbuf_get_index('ZM_MAXG')
         270           0 :    zm_ideep_idx    = pbuf_get_index('ZM_IDEEP')
         271             :
         272           0 :    cmfmc_sh_idx    = pbuf_get_index('CMFMC_SH')
         273           0 :    sh_e_ed_ratio_idx = pbuf_get_index('SH_E_ED_RATIO', istat)
         274             :
         275           0 :    if (masterproc ) then
         276             :
         277           0 :       write(iulog,'(a,l12)')     'aero_convproc_init - convproc_do_aer               = ', &
         278           0 :          convproc_do_aer
         279           0 :       write(iulog,'(a,l12)')     'aero_convproc_init - use_cwaer_for_activate_maxsat = ', &
         280           0 :          use_cwaer_for_activate_maxsat
         281           0 :       write(iulog,'(a,l12)')     'aero_convproc_init - apply_convproc_tend_to_ptend  = ', &
         282           0 :          apply_convproc_tend_to_ptend
         283           0 :       write(iulog,'(a,i12)')     'aero_convproc_init - convproc_method_activate      = ', &
         284           0 :          convproc_method_activate
         285           0 :       write(iulog,'(a,i12)')     'aero_convproc_init - method1_activate_nlayers      = ', &
         286           0 :          method1_activate_nlayers
         287           0 :       write(iulog,'(a,1pe12.4)') 'aero_convproc_init - method2_activate_smaxmax      = ', &
         288           0 :          method2_activate_smaxmax
         289           0 :       write(iulog,'(a,i12)')     'aero_convproc_init - method_reduce_actfrac         = ', &
         290           0 :          method_reduce_actfrac
         291           0 :       write(iulog,'(a,1pe12.4)') 'aero_convproc_init - factor_reduce_actfrac         = ', &
         292           0 :          factor_reduce_actfrac
         293             :
         294           0 :       npass_calc_updraft = 1
         295             :       if ( (method_reduce_actfrac == 2)      .and. &
         296             :          (factor_reduce_actfrac >= 0.0_r8) .and. &
         297             :          (factor_reduce_actfrac <= 1.0_r8) ) npass_calc_updraft = 2
         298           0 :       write(iulog,'(a,i12)')     'aero_convproc_init - npass_calc_updraft            = ', &
         299           0 :          npass_calc_updraft
         300             :
         301             :    end if
         302             :
         303           0 : end subroutine aero_convproc_init
         304             :
         305             : !=========================================================================================
         306             :
         307           0 : subroutine aero_convproc_intr( aero_props, aero_state, state, ptend, pbuf, ztodt,             &
         308           0 :                            nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr,  &
         309           0 :                            aerdepwetis, dcondt_resusp3d )
         310             : !-----------------------------------------------------------------------
         311             : !
         312             : ! Convective cloud processing (transport, activation/resuspension,
         313             : !    wet removal) of aerosols and trace gases.
         314             : !    (Currently no aqueous chemistry and no trace-gas wet removal)
         315             : ! Does aerosols    when convproc_do_aer is .true.
         316             : !
         317             : ! Does deep convection
         318             : ! Uses mass fluxes, cloud water, precip production from the
         319             : !    convective cloud routines
         320             : !
         321             : ! Author: R. Easter
         322             : !
         323             : !-----------------------------------------------------------------------
         324             :
         325             :
         326             :   ! Arguments
         327             :    class(aerosol_properties), intent(in) :: aero_props
         328             :    class(aerosol_state), intent(in) :: aero_state
         329             :
         330             :    type(physics_state),target,intent(in )   :: state      ! Physics state variables
         331             :    type(physics_ptend),       intent(inout) :: ptend      ! %lq set in aero_model_wetdep
         332             :    type(physics_buffer_desc), pointer       :: pbuf(:)
         333             :    real(r8), intent(in) :: ztodt                          ! 2 delta t (model time increment)
         334             :
         335             :    integer,  intent(in)    :: nsrflx_mzaer2cnvpr
         336             :    real(r8), intent(in)    :: qsrflx_mzaer2cnvpr(pcols,ncnstaer,nsrflx_mzaer2cnvpr)
         337             :    real(r8), intent(inout) :: aerdepwetis(pcols,pcnst)  ! aerosol wet deposition (interstitial)
         338             :    real(r8), intent(inout) :: dcondt_resusp3d(ncnstaer,pcols,pver)
         339             :
         340             :    ! Local variables
         341             :    integer, parameter :: nsrflx = 5        ! last dimension of qsrflx
         342             :    integer  :: l, m, mm, ndx, lchnk
         343             :    integer  :: ncol
         344             :
         345           0 :    real(r8) :: dqdt(pcols,pver,ncnstaer)
         346             :    real(r8) :: dt
         347             :
         348             :
         349             :
         350           0 :    real(r8) :: q(pcols,pver,ncnstaer)
         351           0 :    real(r8) :: qsrflx(pcols,ncnstaer,nsrflx)
         352           0 :    real(r8), pointer :: qptr(:,:)
         353             :
         354           0 :    real(r8) :: sflxic(pcols,ncnstaer)
         355           0 :    real(r8) :: sflxid(pcols,ncnstaer)
         356           0 :    real(r8) :: sflxec(pcols,ncnstaer)
         357           0 :    real(r8) :: sflxed(pcols,ncnstaer)
         358             :
         359           0 :    type(ptr2d_t) :: raer(ncnstaer)     ! aerosol mass, number mixing ratios
         360           0 :    type(ptr2d_t) :: qqcw(ncnstaer)
         361             :
         362             :    logical  :: dotend(pcnst)
         363             :    logical  :: applytend
         364             :
         365             :    !-------------------------------------------------------------------------------------------------
         366             :
         367           0 :    dotend = .false.
         368             :
         369             :    ! Initialize
         370           0 :    lchnk = state%lchnk
         371           0 :    ncol  = state%ncol
         372           0 :    dt = ztodt
         373             :
         374           0 :    sflxic(:,:) = 0.0_r8
         375           0 :    sflxid(:,:) = 0.0_r8
         376           0 :    sflxec(:,:) = 0.0_r8
         377           0 :    sflxed(:,:) = 0.0_r8
         378             :
         379           0 :    call aero_state%get_states( aero_props, raer, qqcw )
         380             :
         381             :    ! prepare for deep conv processing
         382           0 :    do m = 1, aero_props%nbins()
         383           0 :       do l = 0, aero_props%nmasses(m)
         384             :
         385           0 :          mm = aero_props%indexer(m,l)
         386           0 :          ndx = aer_cnst_ndx(mm)
         387             :
         388           0 :          sflxec(1:ncol,mm) = qsrflx_mzaer2cnvpr(1:ncol,mm,1)
         389           0 :          sflxed(1:ncol,mm) = qsrflx_mzaer2cnvpr(1:ncol,mm,2)
         390             :
         391           0 :          applytend = .false.
         392           0 :          if ( ndx > 0 ) then
         393           0 :             applytend = ptend%lq(ndx)
         394           0 :             dotend(ndx) = applytend
         395             :          endif
         396             :
         397           0 :          qptr => raer(mm)%fld
         398             :
         399           0 :          if ( applytend ) then
         400             :             ! calc new q (after calcaersize and mz_aero_wet_intr)
         401           0 :             q(1:ncol,:,mm) = max( 0.0_r8, qptr(1:ncol,:) + dt*ptend%q(1:ncol,:,ndx) )
         402             :          else
         403             :             ! use old q
         404           0 :             q(1:ncol,:,mm) = qptr(1:ncol,:)
         405             :          end if
         406             :
         407             :       end do
         408             :    end do
         409             :
         410           0 :    dqdt(:,:,:) = 0.0_r8
         411           0 :    qsrflx(:,:,:) = 0.0_r8
         412             :
         413           0 :    if (convproc_do_aer) then
         414             :
         415             :       ! do deep conv processing
         416           0 :       if (convproc_do_deep) then
         417             :          call aero_convproc_dp_intr( aero_props, &
         418             :             state, pbuf, dt,                          &
         419           0 :             q, dqdt, nsrflx, qsrflx, dcondt_resusp3d )
         420             :
         421             :          ! apply deep conv processing tendency
         422             :
         423           0 :          do m = 1, aero_props%nbins()
         424           0 :             do l = 0, aero_props%nmasses(m)
         425           0 :                mm = aero_props%indexer(m,l)
         426           0 :                ndx = aer_cnst_ndx(mm)
         427             :
         428             :                if ( apply_convproc_tend_to_ptend ) then
         429             :                   ! add dqdt onto ptend%q and set ptend%lq
         430           0 :                   if (ndx>0) then ! advected species
         431           0 :                      ptend%q(1:ncol,:,ndx) = ptend%q(1:ncol,:,ndx) + dqdt(1:ncol,:,mm)
         432             :                   else
         433           0 :                      raer(mm)%fld(1:ncol,:) = max( 0.0_r8, raer(mm)%fld(1:ncol,:) + dqdt(1:ncol,:,mm) * dt )
         434             :                   end if
         435             :                end if
         436             :
         437             :                ! these used for history file wetdep diagnostics
         438           0 :                sflxic(1:ncol,mm) = sflxic(1:ncol,mm) + qsrflx(1:ncol,mm,4)
         439           0 :                sflxid(1:ncol,mm) = sflxid(1:ncol,mm) + qsrflx(1:ncol,mm,4)
         440           0 :                sflxec(1:ncol,mm) = sflxec(1:ncol,mm) + qsrflx(1:ncol,mm,5)
         441           0 :                sflxed(1:ncol,mm) = sflxed(1:ncol,mm) + qsrflx(1:ncol,mm,5)
         442             :
         443             :                ! this used for surface coupling
         444           0 :                if (ndx>0) then
         445           0 :                   aerdepwetis(1:ncol,ndx) = aerdepwetis(1:ncol,ndx) &
         446           0 :                        + qsrflx(1:ncol,mm,4) + qsrflx(1:ncol,mm,5)
         447             :                end if
         448             :             end do
         449             :          end do
         450             :
         451             :       end if
         452             :
         453             :    end if ! (convproc_do_aer) then
         454             :
         455           0 :    if (convproc_do_aer .and. apply_convproc_tend_to_ptend ) then
         456             :
         457           0 :       do m = 1, aero_props%nbins()
         458           0 :          do l = 0, aero_props%nmasses(m)
         459           0 :             mm = aero_props%indexer(m,l)
         460             :
         461           0 :             ndx = aer_cnst_ndx(mm)
         462             :
         463           0 :             if (ndx>0) call outfld( trim(cnst_name_extd(1,mm))//'SFWET', aerdepwetis(:,ndx), pcols, lchnk )
         464           0 :             call outfld( trim(cnst_name_extd(1,mm))//'SFSIC', sflxic(:,mm), pcols, lchnk )
         465           0 :             call outfld( trim(cnst_name_extd(1,mm))//'SFSEC', sflxec(:,mm), pcols, lchnk )
         466             :
         467           0 :             if ( deepconv_wetdep_history ) then
         468           0 :                call outfld( trim(cnst_name_extd(1,mm))//'SFSID', sflxid(:,mm), pcols, lchnk )
         469           0 :                call outfld( trim(cnst_name_extd(1,mm))//'SFSED', sflxed(:,mm), pcols, lchnk )
         470             :             end if
         471             :          end do
         472             :       end do
         473             :
         474             :    end if
         475             :
         476           0 : end subroutine aero_convproc_intr
         477             :
         478             : !=========================================================================================
         479             :
         480           0 : subroutine aero_convproc_dp_intr( aero_props,  &
         481             :      state, pbuf, dt,                          &
         482           0 :      q, dqdt, nsrflx, qsrflx,  dcondt_resusp3d)
         483             : !-----------------------------------------------------------------------
         484             : !
         485             : ! Convective cloud processing (transport, activation/resuspension,
         486             : !    wet removal) of aerosols and trace gases.
         487             : !    (Currently no aqueous chemistry and no trace-gas wet removal)
         488             : ! Does aerosols    when convproc_do_aer is .true.
         489             : !
         490             : ! This routine does deep convection
         491             : ! Uses mass fluxes, cloud water, precip production from the
         492             : !    convective cloud routines
         493             : !
         494             : ! Author: R. Easter
         495             : !
         496             : !-----------------------------------------------------------------------
         497             :
         498             :    ! Arguments
         499             :    class(aerosol_properties), intent(in) :: aero_props
         500             :
         501             :    type(physics_state),       intent(in ) :: state          ! Physics state variables
         502             :    type(physics_buffer_desc), pointer     :: pbuf(:)
         503             :
         504             :    real(r8), intent(in) :: dt                         ! delta t (model time increment)
         505             :
         506             :    real(r8), intent(in)    :: q(pcols,pver,ncnstaer)
         507             :    real(r8), intent(inout) :: dqdt(pcols,pver,ncnstaer)
         508             :    integer,  intent(in)    :: nsrflx
         509             :    real(r8), intent(inout) :: qsrflx(pcols,ncnstaer,nsrflx)
         510             :    real(r8), intent(inout) :: dcondt_resusp3d(ncnstaer,pcols,pver)
         511             :
         512             :    integer :: i
         513             :    integer :: lchnk
         514             :    integer :: nstep
         515             :
         516             :    real(r8) :: dpdry(pcols,pver)     ! layer delta-p-dry (mb)
         517             :    real(r8) :: fracice(pcols,pver)   ! Ice fraction of cloud droplets
         518             :    real(r8) :: xx_mfup_max(pcols), xx_wcldbase(pcols), xx_kcldbase(pcols)
         519             :
         520             :    ! physics buffer fields
         521           0 :    real(r8), pointer :: fracis(:,:,:)  ! fraction of transported species that are insoluble
         522           0 :    real(r8), pointer :: rprddp(:,:)    ! Deep conv precip production (kg/kg/s - grid avg)
         523           0 :    real(r8), pointer :: evapcdp(:,:)   ! Deep conv precip evaporation (kg/kg/s - grid avg)
         524           0 :    real(r8), pointer :: icwmrdp(:,:)   ! Deep conv cloud condensate (kg/kg - in cloud)
         525           0 :    real(r8), pointer :: dp_frac(:,:)   ! Deep conv cloud frac (0-1)
         526             :
         527             :    ! deep conv variables
         528           0 :    real(r8), pointer :: du(:,:)        ! Mass detrain rate from updraft (pcols,pver)
         529           0 :    real(r8), pointer :: eu(:,:)        ! Mass entrain rate into updraft (pcols,pver)
         530           0 :    real(r8), pointer :: ed(:,:)        ! Mass entrain rate into downdraft (pcols,pver)
         531             :                                        ! eu, ed, du are "d(massflux)/dp" and are all positive
         532           0 :    real(r8), pointer :: dp(:,:)        ! Delta pressure between interfaces (pcols,pver)
         533           0 :    integer,  pointer :: jt(:)          ! Index of cloud top for each column (pcols)
         534           0 :    integer,  pointer :: maxg(:)        ! Index of cloud bottom for each column (pcols)
         535           0 :    integer,  pointer :: ideep(:)       ! Gathering array (pcols)
         536             :    integer           :: lengath        ! Gathered min lon indices over which to operate
         537             :
         538             :    ! Initialize
         539             :
         540           0 :    lchnk = state%lchnk
         541           0 :    nstep = get_nstep()
         542             :
         543             :    ! Associate pointers with physics buffer fields
         544           0 :    call pbuf_get_field(pbuf, rprddp_idx,      rprddp)
         545           0 :    call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp)
         546           0 :    call pbuf_get_field(pbuf, icwmrdp_idx,     icwmrdp)
         547           0 :    call pbuf_get_field(pbuf, dp_frac_idx,     dp_frac)
         548           0 :    call pbuf_get_field(pbuf, fracis_idx,      fracis)
         549           0 :    call pbuf_get_field(pbuf, zm_eu_idx,       eu)
         550           0 :    call pbuf_get_field(pbuf, zm_du_idx,       du)
         551           0 :    call pbuf_get_field(pbuf, zm_ed_idx,       ed)
         552           0 :    call pbuf_get_field(pbuf, zm_dp_idx,       dp)
         553           0 :    call pbuf_get_field(pbuf, zm_jt_idx,       jt)
         554           0 :    call pbuf_get_field(pbuf, zm_maxg_idx,     maxg)
         555           0 :    call pbuf_get_field(pbuf, zm_ideep_idx,    ideep)
         556             :
         557           0 :    lengath = count(ideep > 0)
         558             :
         559           0 :    fracice(:,:) = 0.0_r8
         560             :
         561             :    ! initialize dpdry (units=mb), which is used for tracers of dry mixing ratio type
         562           0 :    dpdry = 0._r8
         563           0 :    do i = 1, lengath
         564           0 :       dpdry(i,:) = state%pdeldry(ideep(i),:)/100._r8
         565             :    end do
         566             :
         567             :    call aero_convproc_tend( aero_props, 'deep', lchnk,   dt,      &
         568             :                      state%t,    state%pmid, q, du,      eu,      &
         569             :                      ed,         dp,         dpdry,      jt,      &
         570             :                      maxg,       ideep,      1,          lengath, &
         571             :                      dp_frac,    icwmrdp,    rprddp,     evapcdp, &
         572             :                      fracice,     dqdt,      nsrflx,     qsrflx,  &
         573             :                      xx_mfup_max, xx_wcldbase, xx_kcldbase,       &
         574           0 :                      dcondt_resusp3d  )
         575             :
         576           0 :    call outfld( 'DP_MFUP_MAX', xx_mfup_max, pcols, lchnk )
         577           0 :    call outfld( 'DP_WCLDBASE', xx_wcldbase, pcols, lchnk )
         578           0 :    call outfld( 'DP_KCLDBASE', xx_kcldbase, pcols, lchnk )
         579             :
         580           0 : end subroutine aero_convproc_dp_intr
         581             :
         582             : !=========================================================================================
         583             :
         584           0 : subroutine aero_convproc_tend( aero_props, convtype, lchnk, dt,  &
         585           0 :                      t,          pmid,       q,  du,     eu,         &
         586             :                      ed,         dp,         dpdry,      jt,         &
         587             :                      mx,         ideep,      il1g,       il2g,       &
         588             :                      cldfrac,    icwmr,      rprd,       evapc,      &
         589           0 :                      fracice,    dqdt,       nsrflx,     qsrflx,     &
         590             :                      xx_mfup_max, xx_wcldbase, xx_kcldbase,          &
         591           0 :                      dcondt_resusp3d )
         592             :
         593             : !-----------------------------------------------------------------------
         594             : !
         595             : ! Purpose:
         596             : ! Convective transport of trace species.
         597             : ! The trace species need not be conservative, and source/sink terms for
         598             : !    activation, resuspension, aqueous chemistry and gas uptake, and
         599             : !    wet removal are all applied.
         600             : ! Currently this works with the ZM deep convection, but we should be able
         601             : !    to adapt it for both Hack and McCaa shallow convection
         602             : !
         603             : ! Compare to subr convproc which does conservative trace species.
         604             : !
         605             : ! Method:
         606             : ! Computes tracer mixing ratios in updraft and downdraft "cells" in a
         607             : ! Lagrangian manner, with source/sinks applied in the updraft other.
         608             : ! Then computes grid-cell-mean tendencies by considering
         609             : !    updraft and downdraft fluxes across layer boundaries
         610             : !    environment subsidence/lifting fluxes across layer boundaries
         611             : !    sources and sinks in the updraft
         612             : !    resuspension of activated species in the grid-cell as a whole
         613             : !
         614             : ! Note1:  A better estimate or calculation of either the updraft velocity
         615             : !         or fractional area is needed.
         616             : ! Note2:  If updraft area is a small fraction of over cloud area,
         617             : !         then aqueous chemistry is underestimated.  These are both
         618             : !         research areas.
         619             : !
         620             : ! Authors: O. Seland and R. Easter, based on convtran by P. Rasch
         621             : !
         622             : !-----------------------------------------------------------------------
         623             :
         624             : !-----------------------------------------------------------------------
         625             : !
         626             : ! Input arguments
         627             : !
         628             :    class(aerosol_properties), intent(in) :: aero_props
         629             :
         630             :    character(len=*), intent(in) :: convtype  ! identifies the type of
         631             :                                              ! convection ("deep", "shcu")
         632             :    integer,  intent(in) :: lchnk             ! chunk identifier
         633             :    real(r8), intent(in) :: dt                ! Model timestep
         634             :    real(r8), intent(in) :: t(pcols,pver)     ! Temperature
         635             :    real(r8), intent(in) :: pmid(pcols,pver)  ! Pressure at model levels
         636             :    real(r8), intent(in) :: q(pcols,pver,ncnstaer) ! Tracer array including moisture
         637             :
         638             :    real(r8), intent(in) :: du(pcols,pver)    ! Mass detrain rate from updraft
         639             :    real(r8), intent(in) :: eu(pcols,pver)    ! Mass entrain rate into updraft
         640             :    real(r8), intent(in) :: ed(pcols,pver)    ! Mass entrain rate into downdraft
         641             : ! *** note1 - mu, md, eu, ed, du, dp, dpdry are GATHERED ARRAYS ***
         642             : ! *** note2 - mu and md units are (mb/s), which is used in the zm_conv code
         643             : !           - eventually these should be changed to (kg/m2/s)
         644             : ! *** note3 - eu, ed, du are "d(massflux)/dp" (with dp units = mb), and are all >= 0
         645             :
         646             :    real(r8), intent(in) :: dp(pcols,pver)    ! Delta pressure between interfaces (mb)
         647             :    real(r8), intent(in) :: dpdry(pcols,pver) ! Delta dry-pressure (mb)
         648             :    integer,  intent(in) :: jt(pcols)         ! Index of cloud top for each column
         649             :    integer,  intent(in) :: mx(pcols)         ! Index of cloud bottom for each column
         650             :    integer,  intent(in) :: ideep(pcols)      ! Gathering array indices
         651             :    integer,  intent(in) :: il1g              ! Gathered min lon indices over which to operate
         652             :    integer,  intent(in) :: il2g              ! Gathered max lon indices over which to operate
         653             : ! *** note4 -- for il1g <= i <= il2g,  icol = ideep(i) is the "normal" chunk column index
         654             :
         655             :    real(r8), intent(in) :: cldfrac(pcols,pver)  ! Convective cloud fractional area
         656             :    real(r8), intent(in) :: icwmr(pcols,pver)    ! Convective cloud water from zhang
         657             :    real(r8), intent(in) :: rprd(pcols,pver)     ! Convective precipitation formation rate
         658             :    real(r8), intent(in) :: evapc(pcols,pver)    ! Convective precipitation evaporation rate
         659             :    real(r8), intent(in) :: fracice(pcols,pver)  ! Ice fraction of cloud droplets
         660             :
         661             :    real(r8), intent(out):: dqdt(pcols,pver,ncnstaer)  ! Tracer tendency array
         662             :    integer,  intent(in) :: nsrflx            ! last dimension of qsrflx
         663             :    real(r8), intent(out):: qsrflx(pcols,ncnstaer,nsrflx)
         664             :                               ! process-specific column tracer tendencies
         665             :                               ! (1=activation,  2=resuspension, 3=aqueous rxn,
         666             :                               !  4=wet removal, 5=renaming)
         667             :    real(r8), intent(out) :: xx_mfup_max(pcols)
         668             :    real(r8), intent(out) :: xx_wcldbase(pcols)
         669             :    real(r8), intent(out) :: xx_kcldbase(pcols)
         670             :    real(r8), intent(inout) :: dcondt_resusp3d(ncnstaer,pcols,pver)
         671             :
         672             : !--------------------------Local Variables------------------------------
         673             :
         674             : ! cloudborne aerosol, so the arrays are dimensioned with pcnst_extd = pcnst*2
         675             :
         676             :    integer :: i, icol         ! Work index
         677             :    integer :: iconvtype       ! 1=deep, 2=uw shallow
         678             :    integer :: iflux_method    ! 1=as in convtran (deep), 2=simpler
         679             :    integer :: ipass_calc_updraft
         680             :    integer :: jtsub           ! Work index
         681             :    integer :: k               ! Work index
         682             :    integer :: kactcnt         ! Counter for no. of levels having activation
         683             :    integer :: kactcntb        ! Counter for activation diagnostic output
         684             :    integer :: kactfirst       ! Lowest layer with activation (= cloudbase)
         685             :    integer :: kbot            ! Cloud-flux bottom layer for current i (=mx(i))
         686             :    integer :: kbot_prevap     ! Lowest layer for doing resuspension from evaporating precip
         687             :    integer :: ktop            ! Cloud-flux top    layer for current i (=jt(i))
         688             :                               ! Layers between kbot,ktop have mass fluxes
         689             :                               !    but not all have cloud water, because the
         690             :                               !    updraft starts below the cloud base
         691             :    integer :: km1, km1x       ! Work index
         692             :    integer :: kp1, kp1x       ! Work index
         693             :    integer :: l, mm           ! Work index
         694             :    integer :: m, n, ndx       ! Work index
         695             :    integer :: nerr            ! number of errors for entire run
         696             :    integer :: nerrmax         ! maximum number of errors to report
         697             :    integer :: npass_calc_updraft
         698             :    integer :: ntsub           !
         699             :
         700             :    logical  do_act_this_lev             ! flag for doing activation at current level
         701             :
         702           0 :    real(r8) aqfrac(2,ncnstaer)       ! aqueous fraction of constituent in updraft
         703             :    real(r8) cldfrac_i(pver)          ! cldfrac at current i (with adjustments)
         704             :
         705           0 :    real(r8) chat(2,ncnstaer,pverp)   ! mix ratio in env at interfaces
         706           0 :    real(r8) cond(2,ncnstaer,pverp)   ! mix ratio in downdraft at interfaces
         707           0 :    real(r8) const(2,ncnstaer,pver)   ! gathered tracer array
         708           0 :    real(r8) conu(2,ncnstaer,pverp)   ! mix ratio in updraft at interfaces
         709             :
         710           0 :    real(r8) dcondt(2,ncnstaer,pver)  ! grid-average TMR tendency for current column
         711           0 :    real(r8) dcondt_prevap(2,ncnstaer,pver) ! portion of dcondt from precip evaporation
         712           0 :    real(r8) dcondt_resusp(2,ncnstaer,pver) ! portion of dcondt from resuspension
         713             :
         714           0 :    real(r8) dcondt_wetdep(2,ncnstaer,pver) ! portion of dcondt from wet deposition
         715           0 :    real(r8) dconudt_activa(2,ncnstaer,pverp) ! d(conu)/dt by activation
         716           0 :    real(r8) dconudt_aqchem(2,ncnstaer,pverp) ! d(conu)/dt by aqueous chem
         717           0 :    real(r8) dconudt_wetdep(2,ncnstaer,pverp) ! d(conu)/dt by wet removal
         718             :
         719           0 :    real(r8) maxflux(2,ncnstaer)      ! maximum (over layers) of fluxin and fluxout
         720           0 :    real(r8) maxflux2(2,ncnstaer)     ! ditto but computed using method-2 fluxes
         721           0 :    real(r8) maxprevap(2,ncnstaer)    ! maximum (over layers) of dcondt_prevap*dp
         722           0 :    real(r8) maxresusp(2,ncnstaer)    ! maximum (over layers) of dcondt_resusp*dp
         723           0 :    real(r8) maxsrce(2,ncnstaer)      ! maximum (over layers) of netsrce
         724             :
         725           0 :    real(r8) sumflux(2,ncnstaer)      ! sum (over layers) of netflux
         726           0 :    real(r8) sumflux2(2,ncnstaer)     ! ditto but computed using method-2 fluxes
         727           0 :    real(r8) sumsrce(2,ncnstaer)      ! sum (over layers) of dp*netsrce
         728           0 :    real(r8) sumchng(2,ncnstaer)      ! sum (over layers) of dp*dcondt
         729           0 :    real(r8) sumchng3(2,ncnstaer)     ! ditto but after call to resusp_conv
         730           0 :    real(r8) sumprevap(2,ncnstaer)    ! sum (over layers) of dp*dcondt_prevap
         731           0 :    real(r8) sumwetdep(2,ncnstaer)    ! sum (over layers) of dp*dconudt_wetdep
         732             :
         733             :    real(r8) cabv                 ! mix ratio of constituent above
         734             :    real(r8) cbel                 ! mix ratio of constituent below
         735             :    real(r8) cdifr                ! normalized diff between cabv and cbel
         736             :    real(r8) cdt(pver)            ! (in-updraft first order wet removal rate) * dt
         737             :    real(r8) clw_cut              ! threshold clw value for doing updraft
         738             :                                  ! transformation and removal
         739             :    real(r8) courantmax           ! maximum courant no.
         740             :    real(r8) dddp(pver)           ! dd(i,k)*dp(i,k) at current i
         741             :    real(r8) dp_i(pver)           ! dp(i,k) at current i
         742             :    real(r8) dt_u(pver)           ! lagrangian transport time in the updraft
         743             :    real(r8) dudp(pver)           ! du(i,k)*dp(i,k) at current i
         744           0 :    real(r8) dqdt_i(pver,ncnstaer)   ! dqdt(i,k,m) at current i
         745             :    real(r8) dtsub                ! dt/ntsub
         746             :    real(r8) dz                   ! working layer thickness (m)
         747             :    real(r8) eddp(pver)           ! ed(i,k)*dp(i,k) at current i
         748             :    real(r8) eudp(pver)           ! eu(i,k)*dp(i,k) at current i
         749             :    real(r8) expcdtm1             ! a work variable
         750             :    real(r8) fa_u(pver)           ! fractional area of in the updraft
         751             :    real(r8) fa_u_dp              ! current fa_u(k)*dp_i(k)
         752             :    real(r8) f_ent                ! fraction of the "before-detrainment" updraft
         753             :                                  ! massflux at k/k-1 interface resulting from
         754             :                                  ! entrainment of level k air
         755             :    real(r8) fluxin               ! a work variable
         756             :    real(r8) fluxout              ! a work variable
         757             :    real(r8) maxc                 ! a work variable
         758             :    real(r8) mbsth                ! Threshold for mass fluxes
         759             :    real(r8) minc                 ! a work variable
         760             :    real(r8) md_m_eddp            ! a work variable
         761             :    real(r8) md_i(pverp)          ! md(i,k) at current i (note pverp dimension)
         762             :    real(r8) md_x(pverp)          ! md(i,k) at current i (note pverp dimension)
         763             :    real(r8) mu_i(pverp)          ! mu(i,k) at current i (note pverp dimension)
         764             :    real(r8) mu_x(pverp)          ! mu(i,k) at current i (note pverp dimension)
         765             :    ! md_i, md_x, mu_i, mu_x are all "dry" mass fluxes
         766             :    ! the mu_x/md_x are initially calculated from the incoming mu/md by applying dp/dpdry
         767             :    ! the mu_i/md_i are next calculated by applying the mbsth threshold
         768             :    real(r8) mu_p_eudp(pver)      ! = mu_i(kp1) + eudp(k)
         769             :    real(r8) netflux              ! a work variable
         770             :    real(r8) netsrce              ! a work variable
         771           0 :    real(r8) q_i(pver,ncnstaer)      ! q(i,k,m) at current i
         772           0 :    real(r8) qsrflx_i(ncnstaer,nsrflx) ! qsrflx(i,m,n) at current i
         773             :    real(r8) rhoair_i(pver)       ! air density at current i
         774             :    real(r8) small                ! a small number
         775             :    real(r8) tmpa                 ! work variables
         776             :    real(r8) tmpf                 ! work variables
         777             :    real(r8) xinv_ntsub           ! 1.0/ntsub
         778             :    real(r8) wup(pver)            ! working updraft velocity (m/s)
         779           0 :    real(r8) conu2(pcols,pver,2,ncnstaer)
         780           0 :    real(r8) dcondt2(pcols,pver,2,ncnstaer)
         781             :
         782             :    !Fractional area of ensemble mean updrafts in ZM scheme set to 0.01
         783             :    !Chosen to reproduce vertical velocities in GATEIII GIGALES (Khairoutdinov etal 2009, JAMES)
         784             :    real(r8), parameter :: zm_areafrac = 0.01_r8
         785             :
         786             : !-----------------------------------------------------------------------
         787             : !
         788           0 :    iconvtype = -1
         789           0 :    iflux_method = -1
         790             :
         791           0 :    if (convtype == 'deep') then
         792             :       iconvtype = 1
         793             :       iflux_method = 1
         794           0 :    else if (convtype == 'uwsh') then
         795             :       iconvtype = 2
         796             :       iflux_method = 2
         797             :    else
         798           0 :       call endrun( '*** aero_convproc_tend -- convtype is not |deep| or |uwsh|' )
         799             :    end if
         800             :
         801           0 :    nerr = 0
         802           0 :    nerrmax = 99
         803             :
         804           0 :    dcondt_resusp3d(:,:,:) = 0._r8
         805             :
         806             :    small = 1.e-36_r8
         807             : ! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
         808             :    mbsth = 1.e-15_r8
         809             :
         810           0 :    qsrflx(:,:,:) = 0.0_r8
         811           0 :    dqdt(:,:,:) = 0.0_r8
         812           0 :    xx_mfup_max(:) = 0.0_r8
         813           0 :    xx_wcldbase(:) = 0.0_r8
         814           0 :    xx_kcldbase(:) = 0.0_r8
         815             :
         816           0 :    wup(:) = 0.0_r8
         817             :
         818           0 :    dcondt2 = 0.0_r8
         819           0 :    conu2 = 0.0_r8
         820           0 :    aqfrac = 0.0_r8
         821             :
         822             : ! inititialize aqfrac to 1.0 for activated aerosol species, 0.0 otherwise
         823           0 :    do m = 1, aero_props%nbins()
         824           0 :       do l = 0, aero_props%nmasses(m)
         825           0 :          mm = aero_props%indexer(m,l)
         826           0 :          aqfrac(2,mm) = 1.0_r8
         827             :       enddo
         828             :    enddo
         829             :
         830             : ! Loop ever each column that has convection
         831             : ! *** i is index to gathered arrays; ideep(i) is index to "normal" chunk arrays
         832             : i_loop_main_aa: &
         833           0 :    do i = il1g, il2g
         834           0 :    icol = ideep(i)
         835             :
         836             :
         837           0 :    if ( (jt(i) <= 0) .and. (mx(i) <= 0) .and. (iconvtype /= 1) ) then
         838             : ! shallow conv case with jt,mx <= 0, which means there is no shallow conv
         839             : ! in this column -- skip this column
         840             :       cycle i_loop_main_aa
         841             :
         842           0 :    else if ( (jt(i) < 1) .or. (mx(i) > pver) .or. (jt(i) > mx(i)) ) then
         843             : ! invalid cloudtop and cloudbase indices -- skip this column
         844           0 :       write(*,9010) 'illegal jt, mx', convtype, lchnk, icol, i,    &
         845           0 :                                       jt(i), mx(i)
         846             : 9010  format( '*** aero_convproc_tend error -- ', a, 5x, 'convtype = ', a /   &
         847             :               '*** lchnk, icol, il, jt, mx = ', 5(1x,i10) )
         848           0 :       cycle i_loop_main_aa
         849             :
         850           0 :    else if (jt(i) == mx(i)) then
         851             : ! cloudtop = cloudbase (1 layer cloud) -- skip this column
         852           0 :       write(*,9010) 'jt == mx', convtype, lchnk, icol, i, jt(i), mx(i)
         853           0 :       cycle i_loop_main_aa
         854             :
         855             :    end if
         856             :
         857             :
         858             : !
         859             : ! cloudtop and cloudbase indices are valid so proceed with calculations
         860             : !
         861             :
         862             : ! Load dp_i and cldfrac_i, and calc rhoair_i
         863           0 :       do k = 1, pver
         864           0 :          dp_i(k) = dpdry(i,k)
         865           0 :          cldfrac_i(k) = cldfrac(icol,k)
         866           0 :          rhoair_i(k) = pmid(icol,k)/(rair*t(icol,k))
         867             :       end do
         868             :
         869             : ! Calc dry mass fluxes
         870             : !    This is approximate because the updraft air is has different temp and qv than
         871             : !    the grid mean, but the whole convective parameterization is highly approximate
         872           0 :       mu_x(:) = 0.0_r8
         873           0 :       md_x(:) = 0.0_r8
         874             : ! (eu-du) = d(mu)/dp -- integrate upwards, multiplying by dpdry
         875           0 :       do k = pver, 1, -1
         876           0 :          mu_x(k) = mu_x(k+1) + (eu(i,k)-du(i,k))*dp_i(k)
         877           0 :          xx_mfup_max(icol) = max( xx_mfup_max(icol), mu_x(k) )
         878             :       end do
         879             : ! (ed) = d(md)/dp -- integrate downwards, multiplying by dpdry
         880           0 :       do k = 2, pver
         881           0 :          md_x(k) = md_x(k-1) - ed(i,k-1)*dp_i(k-1)
         882             :       end do
         883             :
         884             : ! Load mass fluxes over cloud layers
         885             : ! (Note - use of arrays dimensioned k=1,pver+1 simplifies later coding)
         886             : ! Zero out values below threshold
         887             : ! Zero out values at "top of cloudtop", "base of cloudbase"
         888           0 :       ktop = jt(i)
         889           0 :       kbot = mx(i)
         890             : ! usually the updraft ( & downdraft) start ( & end ) at kbot=pver, but sometimes kbot < pver
         891             : ! transport, activation, resuspension, and wet removal only occur between kbot >= k >= ktop
         892             : ! resuspension from evaporating precip can occur at k > kbot when kbot < pver
         893           0 :       kbot_prevap = pver
         894           0 :       mu_i(:) = 0.0_r8
         895           0 :       md_i(:) = 0.0_r8
         896           0 :       do k = ktop+1, kbot
         897           0 :          mu_i(k) = mu_x(k)
         898           0 :          if (mu_i(k) <= mbsth) mu_i(k) = 0.0_r8
         899           0 :          md_i(k) = md_x(k)
         900           0 :          if (md_i(k) >= -mbsth) md_i(k) = 0.0_r8
         901             :       end do
         902           0 :       mu_i(ktop) = 0.0_r8
         903           0 :       md_i(ktop) = 0.0_r8
         904           0 :       mu_i(kbot+1) = 0.0_r8
         905           0 :       md_i(kbot+1) = 0.0_r8
         906             :
         907             : !  Compute updraft and downdraft "entrainment*dp" from eu and ed
         908             : !  Compute "detrainment*dp" from mass conservation
         909           0 :       eudp(:) = 0.0_r8
         910           0 :       dudp(:) = 0.0_r8
         911           0 :       eddp(:) = 0.0_r8
         912           0 :       dddp(:) = 0.0_r8
         913           0 :       courantmax = 0.0_r8
         914           0 :       do k = ktop, kbot
         915           0 :          if ((mu_i(k) > 0) .or. (mu_i(k+1) > 0)) then
         916           0 :             if (du(i,k) <= 0.0_r8) then
         917           0 :                eudp(k) = mu_i(k) - mu_i(k+1)
         918             :             else
         919           0 :                eudp(k) = max( eu(i,k)*dp_i(k), 0.0_r8 )
         920           0 :                dudp(k) = (mu_i(k+1) + eudp(k)) - mu_i(k)
         921           0 :                if (dudp(k) < 1.0e-12_r8*eudp(k)) then
         922           0 :                   eudp(k) = mu_i(k) - mu_i(k+1)
         923           0 :                   dudp(k) = 0.0_r8
         924             :                end if
         925             :             end if
         926             :          end if
         927           0 :          if ((md_i(k) < 0) .or. (md_i(k+1) < 0)) then
         928           0 :             eddp(k) = max( ed(i,k)*dp_i(k), 0.0_r8 )
         929           0 :             dddp(k) = (md_i(k+1) + eddp(k)) - md_i(k)
         930           0 :             if (dddp(k) < 1.0e-12_r8*eddp(k)) then
         931           0 :                eddp(k) = md_i(k) - md_i(k+1)
         932           0 :                dddp(k) = 0.0_r8
         933             :             end if
         934             :          end if
         935           0 :          courantmax = max( courantmax, ( mu_i(k+1)+eudp(k)-md_i(k)+eddp(k) )*dt/dp_i(k) )
         936             :       end do ! k
         937             :
         938             : ! number of time substeps needed to maintain "courant number" <= 1
         939           0 :       ntsub = 1
         940           0 :       if (courantmax > (1.0_r8 + 1.0e-6_r8)) then
         941           0 :          ntsub = 1 + int( courantmax )
         942             :       end if
         943           0 :       xinv_ntsub = 1.0_r8/ntsub
         944           0 :       dtsub = dt*xinv_ntsub
         945           0 :       courantmax = courantmax*xinv_ntsub
         946             :
         947             : !  load tracer mixing ratio array, which will be updated at the end of each jtsub interation
         948           0 :       do m = 1, aero_props%nbins()
         949           0 :          do l = 0, aero_props%nmasses(m)
         950           0 :             mm = aero_props%indexer(m,l)
         951           0 :             q_i(1:pver,mm) = q(icol,1:pver,mm)
         952           0 :             conu2(icol,1:pver,1,mm) = q(icol,1:pver,mm)
         953             :          end do
         954             :       end do
         955             :
         956             : !
         957             : !   when method_reduce_actfrac = 2, need to do the updraft calc twice
         958             : !   (1st to get non-adjusted activation amount, 2nd to apply reduction factor)
         959             :       npass_calc_updraft = 1
         960             :       if ( (method_reduce_actfrac == 2)      .and. &
         961             :            (factor_reduce_actfrac >= 0.0_r8) .and. &
         962             :            (factor_reduce_actfrac <= 1.0_r8) ) npass_calc_updraft = 2
         963             :
         964             :
         965             : jtsub_loop_main_aa: &
         966           0 :       do jtsub = 1, ntsub
         967             :
         968             :
         969             : ipass_calc_updraft_loop: &
         970           0 :       do ipass_calc_updraft = 1, npass_calc_updraft
         971             :
         972           0 :       qsrflx_i(:,:) = 0.0_r8
         973           0 :       dqdt_i(:,:) = 0.0_r8
         974             :
         975           0 :       const = 0.0_r8 ! zero cloud-phase species
         976           0 :       chat = 0.0_r8 ! zero cloud-phase species
         977           0 :       conu = 0.0_r8
         978           0 :       cond = 0.0_r8
         979             :
         980           0 :       dcondt = 0.0_r8
         981           0 :       dcondt_resusp = 0.0_r8
         982           0 :       dcondt_wetdep = 0.0_r8
         983           0 :       dcondt_prevap = 0.0_r8
         984           0 :       dconudt_aqchem = 0.0_r8
         985           0 :       dconudt_wetdep = 0.0_r8
         986             :
         987             : ! only initialize the activation tendency on ipass=1
         988           0 :       if (ipass_calc_updraft == 1) dconudt_activa = 0.0_r8
         989             :
         990             :       ! initialize mixing ratio arrays (chat, const, conu, cond)
         991           0 :       do m = 1, aero_props%nbins()
         992           0 :          do l = 0, aero_props%nmasses(m)
         993           0 :             mm = aero_props%indexer(m,l)
         994             :
         995           0 :             const(1,mm,:) = q_i(:,mm)
         996             :
         997             :             ! From now on work only with gathered data
         998             :             ! Interpolate environment tracer values to interfaces
         999           0 :             do k = 1,pver
        1000           0 :                km1 = max(1,k-1)
        1001           0 :                minc = min(const(1,mm,km1),const(1,mm,k))
        1002           0 :                maxc = max(const(1,mm,km1),const(1,mm,k))
        1003           0 :                if (minc < 0) then
        1004             :                   cdifr = 0._r8
        1005             :                else
        1006           0 :                   cdifr = abs(const(1,mm,k)-const(1,mm,km1))/max(maxc,small)
        1007             :                endif
        1008             :
        1009             :                ! If the two layers differ significantly use a geometric averaging procedure
        1010             :                ! But only do that for deep convection.  For shallow, use the simple
        1011             :                ! averaging which is used in subr cmfmca
        1012           0 :                if (iconvtype /= 1) then
        1013           0 :                   chat(1,mm,k) = 0.5_r8* (const(1,mm,k)+const(1,mm,km1))
        1014           0 :                else if (cdifr > 1.E-6_r8) then
        1015           0 :                   cabv = max(const(1,mm,km1),maxc*1.e-12_r8)
        1016           0 :                   cbel = max(const(1,mm,k),maxc*1.e-12_r8)
        1017           0 :                   chat(1,mm,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
        1018             :                else             ! Small diff, so just arithmetic mean
        1019           0 :                   chat(1,mm,k) = 0.5_r8* (const(1,mm,k)+const(1,mm,km1))
        1020             :                end if
        1021             :
        1022             :                ! Set provisional up and down draft values, and tendencies
        1023           0 :                conu(1,mm,k) = chat(1,mm,k)
        1024           0 :                cond(1,mm,k) = chat(1,mm,k)
        1025             :             end do ! k
        1026             :
        1027             :             ! Values at surface inferface == values in lowest layer
        1028           0 :             chat(1,mm,pver+1) = const(1,mm,pver)
        1029           0 :             conu(1,mm,pver+1) = const(1,mm,pver)
        1030           0 :             cond(1,mm,pver+1) = const(1,mm,pver)
        1031             :          end do ! l
        1032             :       end do ! m
        1033             :
        1034             :
        1035             :
        1036             : ! Compute updraft mixing ratios from cloudbase to cloudtop
        1037             : ! No special treatment is needed at k=pver because arrays
        1038             : !    are dimensioned 1:pver+1
        1039             : ! A time-split approach is used.  First, entrainment is applied to produce
        1040             : !    an initial conu(m,k) from conu(m,k+1).  Next, chemistry/physics are
        1041             : !    applied to the initial conu(m,k) to produce a final conu(m,k).
        1042             : !    Detrainment from the updraft uses this final conu(m,k).
        1043             : ! Note that different time-split approaches would give somewhat different
        1044             : !    results
        1045           0 :       kactcnt = 0 ; kactcntb = 0 ; kactfirst = 1
        1046             : k_loop_main_bb: &
        1047           0 :       do k = kbot, ktop, -1
        1048           0 :          kp1 = k+1
        1049             :
        1050             : ! cldfrac = conv cloud fractional area.  This could represent anvil cirrus area,
        1051             : !    and may not useful for aqueous chem and wet removal calculations
        1052           0 :          cldfrac_i(k) = max( cldfrac_i(k), 0.005_r8 )
        1053             : ! mu_p_eudp(k) = updraft massflux at k, without detrainment between kp1,k
        1054           0 :          mu_p_eudp(k) = mu_i(kp1) + eudp(k)
        1055             :
        1056           0 :          fa_u(k) = 0.0_r8 !BSINGH(10/15/2014): Initialized so that it has a value if the following "if" check yeilds .false.
        1057           0 :          if (mu_p_eudp(k) > mbsth) then
        1058             : ! if (mu_p_eudp(k) <= mbsth) the updraft mass flux is negligible at base and top
        1059             : !    of current layer,
        1060             : ! so current layer is a "gap" between two unconnected updrafts,
        1061             : ! so essentially skip all the updraft calculations for this layer
        1062             :
        1063             : ! First apply changes from entrainment
        1064           0 :             f_ent = eudp(k)/mu_p_eudp(k)
        1065           0 :             f_ent = max( 0.0_r8, min( 1.0_r8, f_ent ) )
        1066           0 :             tmpa = 1.0_r8 - f_ent
        1067           0 :             do n = 1,2 ! phase
        1068           0 :                do m = 1, aero_props%nbins()
        1069           0 :                   do l = 0, aero_props%nmasses(m)
        1070           0 :                      mm = aero_props%indexer(m,l)
        1071           0 :                      conu(n,mm,k) = tmpa*conu(n,mm,kp1) + f_ent*const(n,mm,k)
        1072             :                   end do
        1073             :                end do
        1074             :             end do
        1075             :
        1076             : ! estimate updraft velocity (wup)
        1077           0 :             if (iconvtype /= 1) then
        1078             : ! shallow - wup = (mup in kg/m2/s) / [rhoair * (updraft area)]
        1079           0 :                wup(k) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
        1080           0 :                       / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
        1081             :             else
        1082             : ! deep - as in shallow, but assumed constant updraft_area with height zm_areafrac
        1083           0 :                wup(k) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
        1084           0 :                       / (rhoair_i(k) * zm_areafrac)
        1085             :             end if
        1086             :
        1087             : ! compute lagrangian transport time (dt_u) and updraft fractional area (fa_u)
        1088             : ! *** these must obey    dt_u(k)*mu_p_eudp(k) = dp_i(k)*fa_u(k)
        1089           0 :             dz = dp_i(k)*hund_ovr_g/rhoair_i(k)
        1090           0 :             dt_u(k) = dz/wup(k)
        1091           0 :             dt_u(k) = min( dt_u(k), dt )
        1092           0 :             fa_u(k) = dt_u(k)*(mu_p_eudp(k)/dp_i(k))
        1093             :
        1094             :
        1095             : ! Now apply transformation and removal changes
        1096             : !    Skip levels where icwmr(icol,k) <= clw_cut (= 1.0e-6) to eliminate
        1097             : !    occasional very small icwmr values from the ZM module
        1098           0 :             clw_cut = 1.0e-6_r8
        1099             :
        1100             :
        1101             :             if (convproc_method_activate <= 1) then
        1102             : ! aerosol activation - method 1
        1103             : !    skip levels that are completely glaciated (fracice(icol,k) == 1.0)
        1104             : !    when kactcnt=1 (first/lowest layer with cloud water) apply
        1105             : !       activatation to the entire updraft
        1106             : !    when kactcnt>1 apply activatation to the amount entrained at this level
        1107             :                if ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0_r8)) then
        1108             :                   kactcnt = kactcnt + 1
        1109             :
        1110             :                   if ((kactcnt == 1) .or. (f_ent > 0.0_r8)) then
        1111             :                      kactcntb = kactcntb + 1
        1112             :                   end if
        1113             :
        1114             :                   if (kactcnt == 1) then
        1115             :                      ! diagnostic fields
        1116             :                      ! xx_wcldbase = w at first cloudy layer, estimated from mu and cldfrac
        1117             :                      xx_wcldbase(icol) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
        1118             :                          / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
        1119             :                      xx_kcldbase(icol) = k
        1120             :
        1121             :                      kactfirst = k
        1122             :                      tmpa = 1.0_r8
        1123             :                      call activate_convproc( aero_props, &
        1124             :                         conu(:,:,k), dconudt_activa(:,:,k), conu(:,:,k),  &
        1125             :                         tmpa,       dt_u(k),            wup(k),           &
        1126             :                         t(icol,k),  rhoair_i(k), ipass_calc_updraft )
        1127             :                   else if (f_ent > 0.0_r8) then
        1128             :                      ! current layer is above cloud base (=first layer with activation)
        1129             :                      !    only allow activation at k = kactfirst thru kactfirst-(method1_activate_nlayers-1)
        1130             :                      if (k >= kactfirst-(method1_activate_nlayers-1)) then
        1131             :                         call activate_convproc( aero_props, &
        1132             :                            conu(:,:,k),  dconudt_activa(:,:,k), const(:,:,k), &
        1133             :                            f_ent,      dt_u(k),             wup(k),           &
        1134             :                            t(icol,k),  rhoair_i(k), ipass_calc_updraft  )
        1135             :                      end if
        1136             :                   end if
        1137             : ! the following was for cam2 shallow convection (hack),
        1138             : ! but is not appropriate for cam5 (uwshcu)
        1139             : !                 else if ((kactcnt > 0) .and. (iconvtype /= 1)) then
        1140             : ! !    for shallow conv, when you move from activation occuring to
        1141             : ! !       not occuring, reset kactcnt=0, because the hack scheme can
        1142             : ! !       produce multiple "1.5 layer clouds" separated by clear air
        1143             : !                    kactcnt = 0
        1144             : !                 end if
        1145             :                end if ! ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0)) then
        1146             :
        1147             :             else ! (convproc_method_activate >= 2)
        1148             : ! aerosol activation - method 2
        1149             : !    skip levels that are completely glaciated (fracice(icol,k) == 1.0)
        1150             : !    when kactcnt=1 (first/lowest layer with cloud water)
        1151             : !       apply "primary" activatation to the entire updraft
        1152             : !    when kactcnt>1
        1153             : !       apply secondary activatation to the entire updraft
        1154             : !       do this for all levels above cloud base (even if completely glaciated)
        1155             : !          (this is something for sensitivity testing)
        1156           0 :                do_act_this_lev = .false.
        1157           0 :                if (kactcnt <= 0) then
        1158           0 :                   if (icwmr(icol,k) > clw_cut) then
        1159           0 :                      do_act_this_lev = .true.
        1160           0 :                      kactcnt = 1
        1161           0 :                      kactfirst = k
        1162             :                      ! diagnostic fields
        1163             :                      ! xx_wcldbase = w at first cloudy layer, estimated from mu and cldfrac
        1164           0 :                      xx_wcldbase(icol) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
        1165           0 :                          / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
        1166           0 :                      xx_kcldbase(icol) = k
        1167             :                   end if
        1168             :                else
        1169             : !                 if ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0)) then
        1170           0 :                      do_act_this_lev = .true.
        1171           0 :                      kactcnt = kactcnt + 1
        1172             : !                 end if
        1173             :                end if
        1174             :
        1175             :                if ( do_act_this_lev ) then
        1176           0 :                   kactcntb = kactcntb + 1
        1177             :
        1178             :                   call activate_convproc_method2( aero_props, &
        1179             :                      conu(:,:,k),  dconudt_activa(:,:,k),     &
        1180             :                      f_ent,      dt_u(k),   wup(k),           &
        1181           0 :                      t(icol,k),  rhoair_i(k), k,              &
        1182           0 :                      kactfirst,  ipass_calc_updraft  )
        1183             :
        1184             :                end if
        1185           0 :                conu2(icol,k,:,:) = conu(:,:,k)
        1186             :
        1187             :             end if ! (convproc_method_activate <= 1)
        1188             :
        1189             : ! aqueous chemistry
        1190             : !    do glaciated levels as aqchem_conv will eventually do acid vapor uptake
        1191             : !    to ice, and aqchem_conv module checks fracice before doing liquid wtr stuff
        1192             : !            if (icwmr(icol,k) > clw_cut) then
        1193             : !              call aqchem_conv( conu(1,k), dconudt_aqchem(1,k), aqfrac,  &
        1194             : !                 t(icol,k), fracice(icol,k), icwmr(icol,k), rhoair_i(k), &
        1195             : !                 lh2o2(icol,k), lo3(icol,k), dt_u(k)                     )
        1196             : !            end if
        1197             :
        1198             : ! wet removal
        1199             : !
        1200             : ! mirage2
        1201             : !    rprd               = precip formation as a grid-cell average (kgW/kgA/s)
        1202             : !    icwmr              = cloud water MR within updraft area (kgW/kgA)
        1203             : !    fupdr              = updraft fractional area (--)
        1204             : !    A = rprd/fupdr     = precip formation rate within updraft area (kgW/kgA/s)
        1205             : !    B = A/icwmr = rprd/(icwmr*fupdr)
        1206             : !                       = first-order removal rate (1/s)
        1207             : !    C = dp/(mup/fupdr) = updraft air residence time in the layer (s)
        1208             : !
        1209             : !    fraction removed = (1.0 - exp(-cdt)) where
        1210             : !                 cdt = B*C = (dp/mup)*rprd/icwmr
        1211             : !
        1212             : !    Note1:  fupdr cancels out in cdt, so need not be specified
        1213             : !    Note2:  dp & mup units need only be consistent (e.g., mb & mb/s)
        1214             : !    Note3:  for shallow conv, cdt = 1-beta (beta defined in Hack scheme)
        1215             : !    Note4:  the "dp" in C above and code below should be the moist dp
        1216             : !
        1217             : ! cam5
        1218             : !    clw_preloss = cloud water MR before loss to precip
        1219             : !                = icwmr + dt*(rprd/fupdr)
        1220             : !    B = A/clw_preloss  = (rprd/fupdr)/(icwmr + dt*rprd/fupdr)
        1221             : !                       = rprd/(fupdr*icwmr + dt*rprd)
        1222             : !                       = first-order removal rate (1/s)
        1223             : !
        1224             : !    fraction removed = (1.0 - exp(-cdt)) where
        1225             : !                 cdt = B*C = (fupdr*dp/mup)*[rprd/(fupdr*icwmr + dt*rprd)]
        1226             : !
        1227             : !    Note1:  *** cdt is now sensitive to fupdr, which we do not really know,
        1228             : !                and is not the same as the convective cloud fraction
        1229             : !    Note2:  dt is appropriate in the above cdt expression, not dtsub
        1230             : !
        1231             : !    Apply wet removal at levels where
        1232             : !       icwmr(icol,k) > clw_cut  AND  rprd(icol,k) > 0.0
        1233             : !    as wet removal occurs in both liquid and ice clouds
        1234             : !
        1235           0 :             cdt(k) = 0.0_r8
        1236           0 :             if ((icwmr(icol,k) > clw_cut) .and. (rprd(icol,k) > 0.0_r8)) then
        1237             : !              if (iconvtype == 1) then
        1238           0 :                   tmpf = 0.5_r8*cldfrac_i(k)
        1239           0 :                   cdt(k) = (tmpf*dp(i,k)/mu_p_eudp(k)) * rprd(icol,k) / &
        1240           0 :                         (tmpf*icwmr(icol,k) + dt*rprd(icol,k))
        1241             : !              else if (k < pver) then
        1242             : !                 if (eudp(k+1) > 0) cdt(k) =   &
        1243             : !                       rprd(icol,k)*dp(i,k)/(icwmr(icol,k)*eudp(k+1))
        1244             : !              end if
        1245             :             end if
        1246           0 :             if (cdt(k) > 0.0_r8) then
        1247           0 :                expcdtm1 = exp(-cdt(k)) - 1.0_r8
        1248             :
        1249           0 :                do m = 1, aero_props%nbins()
        1250           0 :                   do l = 0, aero_props%nmasses(m)
        1251           0 :                      mm = aero_props%indexer(m,l)
        1252           0 :                      do n = 1,2
        1253           0 :                         dconudt_wetdep(n,mm,k) = conu(n,mm,k)*aqfrac(n,mm)*expcdtm1
        1254           0 :                         conu(n,mm,k) = conu(n,mm,k) + dconudt_wetdep(n,mm,k)
        1255           0 :                         dconudt_wetdep(n,mm,k) = dconudt_wetdep(n,mm,k) / dt_u(k)
        1256           0 :                         conu2(icol,k,n,mm) = conu(n,mm,k)
        1257             :                      enddo
        1258             :                   enddo
        1259             :                enddo
        1260             :
        1261             :             end if
        1262             :
        1263             :          end if    ! "(mu_p_eudp(k) > mbsth)"
        1264             :       end do k_loop_main_bb ! "k = kbot, ktop, -1"
        1265             :
        1266             : ! when doing updraft calcs twice, only need to go this far on the first pass
        1267             :       if ( (ipass_calc_updraft == 1) .and. &
        1268             :            (npass_calc_updraft == 2) ) cycle ipass_calc_updraft_loop
        1269             :
        1270             :
        1271             : ! Compute downdraft mixing ratios from cloudtop to cloudbase
        1272             : ! No special treatment is needed at k=2
        1273             : ! No transformation or removal is applied in the downdraft
        1274           0 :       do k = ktop, kbot
        1275           0 :          kp1 = k + 1
        1276             : ! md_m_eddp = downdraft massflux at kp1, without detrainment between k,kp1
        1277           0 :          md_m_eddp = md_i(k) - eddp(k)
        1278           0 :          if (md_m_eddp < -mbsth) then
        1279             :
        1280           0 :             do m = 1, aero_props%nbins()
        1281           0 :                do l = 0, aero_props%nmasses(m)
        1282           0 :                   mm = aero_props%indexer(m,l)
        1283           0 :                   do n = 1,2
        1284           0 :                      cond(n,mm,kp1) = ( md_i(k)*cond(n,mm,k) &
        1285           0 :                                       - eddp(k)*const(n,mm,k) ) / md_m_eddp
        1286             :                   end do
        1287             :                end do
        1288             :             end do
        1289             :          end if
        1290             :       end do ! k
        1291             :
        1292             :
        1293             : ! Now computes fluxes and tendencies
        1294             : ! NOTE:  The approach used in convtran applies to inert tracers and
        1295             : !        must be modified to include source and sink terms
        1296           0 :       sumflux = 0.0_r8
        1297           0 :       sumflux2 = 0.0_r8
        1298           0 :       sumsrce = 0.0_r8
        1299           0 :       sumchng = 0.0_r8
        1300           0 :       sumchng3 = 0.0_r8
        1301           0 :       sumwetdep = 0.0_r8
        1302           0 :       sumprevap = 0.0_r8
        1303             :
        1304           0 :       maxflux = 0.0_r8
        1305           0 :       maxflux2 = 0.0_r8
        1306           0 :       maxresusp = 0.0_r8
        1307           0 :       maxsrce = 0.0_r8
        1308           0 :       maxprevap = 0.0_r8
        1309             :
        1310             : k_loop_main_cc: &
        1311           0 :       do k = ktop, kbot
        1312           0 :          kp1 = k+1
        1313           0 :          km1 = k-1
        1314           0 :          kp1x = min( kp1, pver )
        1315           0 :          km1x = max( km1, 1 )
        1316           0 :          fa_u_dp = fa_u(k)*dp_i(k)
        1317             :
        1318           0 :          do m = 1, aero_props%nbins()
        1319           0 :             do l = 0, aero_props%nmasses(m)
        1320           0 :                mm = aero_props%indexer(m,l)
        1321           0 :                do n = 1,2
        1322             :
        1323             :                   ! First compute fluxes using environment subsidence/lifting and
        1324             :                   ! entrainment/detrainment into up/downdrafts,
        1325             :                   ! to provide an additional mass balance check
        1326             :                   ! (this could be deleted after the code is well tested)
        1327           0 :                   fluxin  = mu_i(k)*min(chat(n,mm,k),const(n,mm,km1x))       &
        1328           0 :                           - md_i(kp1)*min(chat(n,mm,kp1),const(n,mm,kp1x))   &
        1329           0 :                           + dudp(k)*conu(n,mm,k) + dddp(k)*cond(n,mm,kp1)
        1330             :                   fluxout = mu_i(kp1)*min(chat(n,mm,kp1),const(n,mm,k))      &
        1331             :                           - md_i(k)*min(chat(n,mm,k),const(n,mm,k))          &
        1332           0 :                           + (eudp(k) + eddp(k))*const(n,mm,k)
        1333             :
        1334           0 :                   netflux = fluxin - fluxout
        1335             :
        1336           0 :                   sumflux2(n,mm) = sumflux2(n,mm) + netflux
        1337           0 :                   maxflux2(n,mm) = max( maxflux2(n,mm), abs(fluxin), abs(fluxout) )
        1338             :
        1339             :                   ! Now compute fluxes as in convtran, and also source/sink terms
        1340             :                   ! (version 3 limit fluxes outside convection to mass in appropriate layer
        1341             :                   ! (these limiters are probably only safe for positive definite quantitities
        1342             :                   ! (it assumes that mu and md already satify a courant number limit of 1)
        1343           0 :                   if (iflux_method /= 2) then
        1344             :                      fluxin  =     mu_i(kp1)*conu(n,mm,kp1)                     &
        1345             :                                  + mu_i(k  )*min(chat(n,mm,k  ),const(n,mm,km1x))  &
        1346             :                                - ( md_i(k  )*cond(n,mm,k)                       &
        1347           0 :                                  + md_i(kp1)*min(chat(n,mm,kp1),const(n,mm,kp1x)) )
        1348             :                      fluxout =     mu_i(k  )*conu(n,mm,k)                       &
        1349             :                                  + mu_i(kp1)*min(chat(n,mm,kp1),const(n,mm,k   ))  &
        1350             :                                - ( md_i(kp1)*cond(n,mm,kp1)                     &
        1351           0 :                                  + md_i(k  )*min(chat(n,mm,k  ),const(n,mm,k   )) )
        1352             :                   else
        1353             :                      fluxin  =     mu_i(kp1)*conu(n,mm,kp1)                     &
        1354           0 :                                - ( md_i(k  )*cond(n,mm,k) )
        1355             :                      fluxout =     mu_i(k  )*conu(n,mm,k)                       &
        1356           0 :                                - ( md_i(kp1)*cond(n,mm,kp1) )
        1357             :
        1358             :                      ! new method -- simple upstream method for the env subsidence
        1359             :                      ! tmpa = net env mass flux (positive up) at top of layer k
        1360           0 :                      tmpa = -( mu_i(k  ) + md_i(k  ) )
        1361           0 :                      if (tmpa <= 0.0_r8) then
        1362           0 :                         fluxin  = fluxin  - tmpa*const(n,mm,km1x)
        1363             :                      else
        1364           0 :                         fluxout = fluxout + tmpa*const(n,mm,k   )
        1365             :                      end if
        1366             :                      ! tmpa = net env mass flux (positive up) at base of layer k
        1367           0 :                      tmpa = -( mu_i(kp1) + md_i(kp1) )
        1368           0 :                      if (tmpa >= 0.0_r8) then
        1369           0 :                         fluxin  = fluxin  + tmpa*const(n,mm,kp1x)
        1370             :                      else
        1371           0 :                         fluxout = fluxout - tmpa*const(n,mm,k   )
        1372             :                      end if
        1373             :                   end if
        1374             :
        1375           0 :                   netflux = fluxin - fluxout
        1376             :                   netsrce = fa_u_dp*(dconudt_aqchem(n,mm,k) + &
        1377           0 :                        dconudt_activa(n,mm,k) + dconudt_wetdep(n,mm,k))
        1378           0 :                   dcondt(n,mm,k) = (netflux+netsrce)/dp_i(k)
        1379             :
        1380           0 :                   dcondt_wetdep(n,mm,k) = fa_u_dp*dconudt_wetdep(n,mm,k)/dp_i(k)
        1381           0 :                   sumwetdep(n,mm) = sumwetdep(n,mm) + fa_u_dp*dconudt_wetdep(n,mm,k)
        1382             :
        1383           0 :                   dcondt2(icol,k,n,mm) = dcondt(n,mm,k)
        1384             :
        1385             :                end do
        1386             :             end do
        1387             :
        1388             :          end do
        1389             :       end do k_loop_main_cc ! "k = ktop, kbot"
        1390             :
        1391             : ! calculate effects of precipitation evaporation
        1392             :       call precpevap_convproc( aero_props, dcondt, dcondt_wetdep,  dcondt_prevap,   &
        1393             :                                   rprd,   evapc,          dp_i,            &
        1394           0 :                                   icol,   ktop            )
        1395             :
        1396             : ! make adjustments to dcondt for activated & unactivated aerosol species
        1397             : !    pairs to account any (or total) resuspension of convective-cloudborne aerosol
        1398           0 :       call resuspend_convproc( aero_props, dcondt, dcondt_resusp, ktop, kbot_prevap )
        1399             :
        1400             :       ! Do resuspension of aerosols from rain only when the rain has
        1401             :       ! totally evaporated.
        1402           0 :       if (convproc_do_evaprain_atonce) then
        1403             :
        1404           0 :          do m = 1, aero_props%nbins()
        1405           0 :             do l = 0, aero_props%nmasses(m)
        1406           0 :                mm = aero_props%indexer(m,l)
        1407           0 :                dcondt_resusp3d(mm,icol,:) = dcondt_resusp(2,mm,:)
        1408             :             end do
        1409             :          end do
        1410             :
        1411           0 :          dcondt_resusp(2,:,:) = 0._r8
        1412             :       end if
        1413             :
        1414             : ! calculate new column-tendency variables
        1415           0 :       do m = 1, aero_props%nbins()
        1416           0 :          do l = 0, aero_props%nmasses(m)
        1417           0 :             mm = aero_props%indexer(m,l)
        1418           0 :             do n = 1,2
        1419           0 :                do k = ktop, kbot_prevap
        1420           0 :                   sumprevap(n,mm) = sumprevap(n,mm) + dcondt_prevap(n,mm,k)*dp_i(k)
        1421             :                end do
        1422             :             end do
        1423             :          end do
        1424             :       end do
        1425             :
        1426             : !
        1427             : ! note again the aero_convproc_tend does not apply convective cloud processing
        1428             : !    to the stratiform-cloudborne aerosol
        1429             : ! within this routine, cloudborne aerosols are convective-cloudborne
        1430             : !
        1431             : ! before tendencies (dcondt, which is loaded into dqdt) are returned,
        1432             : !    the convective-cloudborne aerosol tendencies must be combined
        1433             : !    with the interstitial tendencies
        1434             : ! resuspend_convproc has already done this for the dcondt
        1435             : !
        1436             : ! the individual process column tendencies (sumwetdep, sumprevap, ...)
        1437             : !    are just diagnostic fields that can be written to history
        1438             : ! tendencies for interstitial and convective-cloudborne aerosol could
        1439             : !    both be passed back and output, if desired
        1440             : ! currently, however, the interstitial and convective-cloudborne tendencies
        1441             : !    are combined (in the next code block) before being passed back (in qsrflx)
        1442             : !
        1443             :
        1444           0 :       do m = 1, aero_props%nbins()
        1445           0 :          do l = 0, aero_props%nmasses(m)
        1446           0 :             mm = aero_props%indexer(m,l)
        1447           0 :             sumwetdep(1,mm) = sumwetdep(1,mm) + sumwetdep(2,mm)
        1448           0 :             sumprevap(1,mm) = sumprevap(1,mm) + sumprevap(2,mm)
        1449             :          enddo
        1450             :       enddo
        1451             :
        1452             : !
        1453             : ! scatter overall tendency back to full array
        1454             : !
        1455           0 :       do m = 1, aero_props%nbins()
        1456           0 :          do l = 0, aero_props%nmasses(m)
        1457           0 :             mm = aero_props%indexer(m,l)
        1458           0 :             ndx = aer_cnst_ndx(mm)
        1459           0 :             do k = ktop, kbot_prevap
        1460           0 :                dqdt_i(k,mm) = dcondt(1,mm,k)
        1461           0 :                dqdt(icol,k,mm) = dqdt(icol,k,mm) + dqdt_i(k,mm)*xinv_ntsub
        1462             :             end do
        1463             :
        1464             :          end do
        1465             :       end do ! m
        1466             :
        1467             : ! scatter column burden tendencies for various processes to qsrflx
        1468           0 :       do m = 1, aero_props%nbins()
        1469           0 :          do l = 0, aero_props%nmasses(m)
        1470           0 :             mm = aero_props%indexer(m,l)
        1471           0 :             qsrflx_i(mm,4) = sumwetdep(1,mm)*hund_ovr_g
        1472           0 :             qsrflx_i(mm,5) = sumprevap(1,mm)*hund_ovr_g
        1473           0 :             qsrflx(icol,mm,1:5) = qsrflx(icol,mm,1:5) + qsrflx_i(mm,1:5)*xinv_ntsub
        1474             :          end do
        1475             :       end do
        1476             :
        1477           0 :       if (jtsub < ntsub) then
        1478             :          ! update the q_i for the next interation of the jtsub loop
        1479           0 :          do m = 1, aero_props%nbins()
        1480           0 :             do l = 0, aero_props%nmasses(m)
        1481           0 :                mm = aero_props%indexer(m,l)
        1482           0 :                ndx = aer_cnst_ndx(mm)
        1483           0 :                do k = ktop, kbot_prevap
        1484           0 :                   q_i(k,mm) = max( (q_i(k,mm) + dqdt_i(k,mm)*dtsub), 0.0_r8 )
        1485             :                end do
        1486             :             end do
        1487             :          end do
        1488             :       end if
        1489             :
        1490             :       end do ipass_calc_updraft_loop
        1491             :
        1492             :       end do jtsub_loop_main_aa  ! of the main "do jtsub = 1, ntsub" loop
        1493             :
        1494             :
        1495             :    end do i_loop_main_aa  ! of the main "do i = il1g, il2g" loop
        1496             :
        1497           0 :    do m = 1, aero_props%nbins()
        1498           0 :       do l = 0, aero_props%nmasses(m)
        1499           0 :          mm = aero_props%indexer(m,l)
        1500             :
        1501           0 :          call outfld( trim(cnst_name_extd(1,mm))//'WETC', dcondt2(:,:,1,mm), pcols, lchnk )
        1502           0 :          call outfld( trim(cnst_name_extd(1,mm))//'CONU', conu2(:,:,1,mm), pcols, lchnk )
        1503           0 :          call outfld( trim(cnst_name_extd(2,mm))//'WETC', dcondt2(:,:,2,mm), pcols, lchnk )
        1504           0 :          call outfld( trim(cnst_name_extd(2,mm))//'CONU', conu2(:,:,2,mm), pcols, lchnk )
        1505             :
        1506             :       end do
        1507             :    end do
        1508             :
        1509           0 : end subroutine aero_convproc_tend
        1510             :
        1511             : !=========================================================================================
        1512           0 :    subroutine precpevap_convproc(  aero_props,      &
        1513           0 :               dcondt,  dcondt_wetdep, dcondt_prevap,           &
        1514             :               rprd,    evapc,         dp_i,                    &
        1515             :               icol,    ktop           )
        1516             : !-----------------------------------------------------------------------
        1517             : !
        1518             : ! Purpose:
        1519             : ! Calculate resuspension of wet-removed aerosol species resulting
        1520             : !    from precip evaporation
        1521             : !
        1522             : ! Author: R. Easter
        1523             : !
        1524             : !-----------------------------------------------------------------------
        1525             :
        1526             : !-----------------------------------------------------------------------
        1527             : ! arguments
        1528             : ! (note:  TMR = tracer mixing ratio)
        1529             :
        1530             :    class(aerosol_properties), intent(in) :: aero_props
        1531             :    real(r8), intent(inout) :: dcondt(2,ncnstaer,pver)
        1532             :                               ! overall TMR tendency from convection
        1533             :    real(r8), intent(in)    :: dcondt_wetdep(2,ncnstaer,pver)
        1534             :                               ! portion of TMR tendency due to wet removal
        1535             :    real(r8), intent(inout) :: dcondt_prevap(2,ncnstaer,pver)
        1536             :                               ! portion of TMR tendency due to precip evaporation
        1537             :                               ! (actually, due to the adjustments made here)
        1538             :                               ! (on entry, this is 0.0)
        1539             :
        1540             :    real(r8), intent(in)    :: rprd(pcols,pver)  ! conv precip production  rate (gathered)
        1541             :    real(r8), intent(in)    :: evapc(pcols,pver)  ! conv precip evaporation rate (gathered)
        1542             :    real(r8), intent(in)    :: dp_i(pver) ! pressure thickness of level (in mb)
        1543             :
        1544             :    integer,  intent(in)    :: icol  ! normal (ungathered) i index for current column
        1545             :    integer,  intent(in)    :: ktop  ! index of top cloud level for current column
        1546             :
        1547             : !-----------------------------------------------------------------------
        1548             : ! local variables
        1549             :    integer  :: k, l, m, mm, n
        1550             :    real(r8) :: del_pr_flux_prod      ! change to precip flux from production  [(kg/kg/s)*mb]
        1551             :    real(r8) :: del_pr_flux_evap      ! change to precip flux from evaporation [(kg/kg/s)*mb]
        1552             :    real(r8) :: del_wd_flux_evap      ! change to wet deposition flux from evaporation [(kg/kg/s)*mb]
        1553             :    real(r8) :: fdel_pr_flux_evap     ! fractional change to precip flux from evaporation
        1554             :    real(r8) :: pr_flux               ! precip flux at base of current layer [(kg/kg/s)*mb]
        1555             :    real(r8) :: pr_flux_old
        1556             :    real(r8) :: tmpdp                 ! delta-pressure (mb)
        1557           0 :    real(r8) :: wd_flux(2,ncnstaer)   ! tracer wet deposition flux at base of current layer [(kg/kg/s)*mb]
        1558             : !-----------------------------------------------------------------------
        1559             :
        1560           0 :    pr_flux = 0.0_r8
        1561           0 :    wd_flux = 0.0_r8
        1562             :
        1563           0 :    do k = ktop, pver
        1564           0 :       tmpdp = dp_i(k)
        1565             :
        1566           0 :       pr_flux_old = pr_flux
        1567           0 :       del_pr_flux_prod = tmpdp*max(0.0_r8, rprd(icol,k))
        1568           0 :       pr_flux = pr_flux_old + del_pr_flux_prod
        1569             :
        1570           0 :       del_pr_flux_evap = min( pr_flux, tmpdp*max(0.0_r8, evapc(icol,k)) )
        1571             :
        1572             :       ! Do resuspension of aerosols from rain only when the rain has
        1573             :       ! totally evaporated in one layer.
        1574           0 :       if (convproc_do_evaprain_atonce .and. &
        1575           0 :           (del_pr_flux_evap.ne.pr_flux)) del_pr_flux_evap = 0._r8
        1576             :
        1577           0 :       fdel_pr_flux_evap = del_pr_flux_evap / max(pr_flux, 1.0e-35_r8)
        1578             :
        1579           0 :       do m = 1, aero_props%nbins()
        1580           0 :          do l = 0, aero_props%nmasses(m)
        1581           0 :             mm = aero_props%indexer(m,l)
        1582           0 :             do n = 1,2
        1583             :
        1584             :                ! use -dcondt_wetdep(m,k) as it is negative (or zero)
        1585           0 :                wd_flux(n,mm) = wd_flux(n,mm) + tmpdp*max(0.0_r8, -dcondt_wetdep(n,mm,k))
        1586           0 :                del_wd_flux_evap = wd_flux(n,mm)*fdel_pr_flux_evap
        1587             :
        1588           0 :                dcondt_prevap(n,mm,k) = del_wd_flux_evap/tmpdp
        1589             :
        1590             :             end do
        1591             :          end do
        1592             :       end do
        1593             :
        1594             :       ! resuspension --> create larger aerosols
        1595           0 :       if (convproc_do_evaprain_atonce) then
        1596           0 :          call aero_props%resuspension_resize( dcondt_prevap(1,:,k) )
        1597             :       endif
        1598             :
        1599           0 :       do m = 1, aero_props%nbins()
        1600           0 :          do l = 0, aero_props%nmasses(m)
        1601           0 :             mm = aero_props%indexer(m,l)
        1602           0 :             do n = 1,2
        1603           0 :                dcondt(n,mm,k) = dcondt(n,mm,k) + dcondt_prevap(n,mm,k)
        1604             :             end do
        1605             :          end do
        1606             :       end do
        1607             :
        1608           0 :       pr_flux = max( 0.0_r8, pr_flux-del_pr_flux_evap )
        1609             :
        1610             :    end do ! k
        1611             :
        1612           0 :    end subroutine precpevap_convproc
        1613             :
        1614             : !=========================================================================================
        1615             :    subroutine activate_convproc( aero_props,    &
        1616             :               conu,       dconudt,   conent,    &
        1617             :               f_ent,      dt_u,      wup,       &
        1618             :               tair,       rhoair, ipass_calc_updraft )
        1619             : !-----------------------------------------------------------------------
        1620             : !
        1621             : ! Purpose:
        1622             : ! Calculate activation of aerosol species in convective updraft
        1623             : ! for a single column and level
        1624             : !
        1625             : ! Method:
        1626             : ! conu(l)    = Updraft TMR (tracer mixing ratio) at k/k-1 interface
        1627             : ! conent(l)  = TMR of air that is entrained into the updraft from level k
        1628             : ! f_ent      = Fraction of the "before-detrainment" updraft massflux at
        1629             : !              k/k-1 interface" resulting from entrainment of level k air
        1630             : !              (where k is the current level in subr aero_convproc_tend)
        1631             : !
        1632             : ! On entry to this routine, the conu(l) represents the updraft TMR
        1633             : ! after entrainment, but before chemistry/physics and detrainment,
        1634             : ! and is equal to
        1635             : !    conu(l) = f_ent*conent(l) + (1.0-f_ent)*conu_below(l)
        1636             : ! where
        1637             : !    conu_below(l) = updraft TMR at the k+1/k interface, and
        1638             : !    f_ent   = (eudp/mu_p_eudp) is the fraction of the updraft massflux
        1639             : !              from level k entrainment
        1640             : !
        1641             : ! This routine applies aerosol activation to the entrained tracer,
        1642             : ! then adjusts the conu so that on exit,
        1643             : !   conu(la) = conu_incoming(la) - f_ent*conent(la)*f_act(la)
        1644             : !   conu(lc) = conu_incoming(lc) + f_ent*conent(la)*f_act(la)
        1645             : ! where
        1646             : !   la, lc   = indices for an unactivated/activated aerosol component pair
        1647             : !   f_act    = fraction of conent(la) that is activated.  The f_act are
        1648             : !              calculated with the Razzak-Ghan activation parameterization.
        1649             : !              The f_act differ for each mode, and for number/surface/mass.
        1650             : !
        1651             : ! Note:  At the lowest layer with cloud water, subr convproc calls this
        1652             : ! routine with conent==conu and f_ent==1.0, with the result that
        1653             : ! activation is applied to the entire updraft tracer flux
        1654             : !
        1655             : ! *** The updraft velocity used for activation calculations is rather
        1656             : !     uncertain and needs more work.  However, an updraft of 1-3 m/s
        1657             : !     will activate essentially all of accumulation and coarse mode particles.
        1658             : !
        1659             : ! Author: R. Easter
        1660             : !
        1661             : !-----------------------------------------------------------------------
        1662             :
        1663             :    use ndrop, only: activate_aerosol
        1664             :
        1665             : !-----------------------------------------------------------------------
        1666             : ! arguments  (note:  TMR = tracer mixing ratio)
        1667             :
        1668             :    class(aerosol_properties), intent(in) :: aero_props
        1669             :
        1670             :    ! conu = tracer mixing ratios in updraft at top of this (current) level
        1671             :    !        The conu are changed by activation
        1672             :    real(r8), intent(inout) :: conu(2,ncnstaer)
        1673             :    ! conent = TMRs in the entrained air at this level
        1674             :    real(r8), intent(in)    :: conent(2,ncnstaer)
        1675             :    real(r8), intent(inout) :: dconudt(2,ncnstaer) ! TMR tendencies due to activation
        1676             :
        1677             :    real(r8), intent(in)    :: f_ent  ! fraction of updraft massflux that was
        1678             :                                      ! entrained across this layer == eudp/mu_p_eudp
        1679             :    real(r8), intent(in)    :: dt_u   ! lagrangian transport time (s) in the
        1680             :                                      ! updraft at current level
        1681             :    real(r8), intent(in)    :: wup    ! mean updraft vertical velocity (m/s)
        1682             :                                      ! at current level updraft
        1683             :
        1684             :    real(r8), intent(in)    :: tair   ! Temperature in Kelvin
        1685             :    real(r8), intent(in)    :: rhoair ! air density (kg/m3)
        1686             :
        1687             :    integer,  intent(in)    :: ipass_calc_updraft
        1688             :
        1689             : !-----------------------------------------------------------------------
        1690             : ! local variables
        1691             :    integer  :: l, m, mm
        1692             :
        1693             :    real(r8) :: delact           ! working variable
        1694             :    real(r8) :: dt_u_inv         ! 1.0/dt_u
        1695             :    real(r8) :: fluxm(nbins)     ! to understand this, see subr activate_aerosol
        1696             :    real(r8) :: fluxn(nbins)     ! to understand this, see subr activate_aerosol
        1697             :    real(r8) :: flux_fullact     ! to understand this, see subr activate_aerosol
        1698             :    real(r8) :: fm(nbins)        ! mass fraction of aerosols activated
        1699             :    real(r8) :: fn(nbins)        ! number fraction of aerosols activated
        1700             :    real(r8) :: hygro(nbins)     ! current hygroscopicity for int+act
        1701             :    real(r8) :: naerosol(nbins)  ! interstitial+activated number conc (#/m3)
        1702             :    real(r8) :: sigw             ! standard deviation of updraft velocity (cm/s)
        1703             :    real(r8) :: tmp_fact         ! working variable
        1704             :    real(r8) :: vaerosol(nbins)  ! int+act volume (m3/m3)
        1705             :    real(r8) :: wbar             ! mean updraft velocity (cm/s)
        1706             :    real(r8) :: wdiab            ! diabatic vertical velocity (cm/s)
        1707             :    real(r8) :: wminf, wmaxf     ! limits for integration over updraft spectrum (cm/s)
        1708             :
        1709             :    real(r8) :: spec_hygro
        1710             :    real(r8) :: spec_dens
        1711             :    character(len=32) :: spec_type
        1712             :
        1713             :    real(r8) :: tmpa, tmpb, tmpc ! working variable
        1714             :    real(r8) :: naerosol_a(1)    ! number conc (1/m3)
        1715             :    real(r8) :: vaerosol_a(1)    ! volume conc (m3/m3)
        1716             :
        1717             : !-----------------------------------------------------------------------
        1718             :
        1719             : ! when ipass_calc_updraft == 2, apply the activation tendencies
        1720             : !    from pass 1, but multiplied by factor_reduce_actfrac
        1721             : ! (can only have ipass_calc_updraft == 2 when method_reduce_actfrac = 2)
        1722             :    if (ipass_calc_updraft == 2) then
        1723             :
        1724             :       dt_u_inv = 1.0_r8/dt_u
        1725             :
        1726             :       do m = 1, aero_props%nbins()
        1727             :          do l = 0, aero_props%nmasses(m)
        1728             :             mm = aero_props%indexer(m,l)
        1729             :
        1730             :             delact = dconudt(2,mm)*dt_u * factor_reduce_actfrac
        1731             :             delact = min( delact, conu(1,mm) )
        1732             :             delact = max( delact, 0.0_r8 )
        1733             :             conu(1,mm) = conu(1,mm) - delact
        1734             :             conu(2,mm) = conu(2,mm) + delact
        1735             :             dconudt(1,mm) = -delact*dt_u_inv
        1736             :             dconudt(2,mm) =  delact*dt_u_inv
        1737             :
        1738             :          end do
        1739             :       end do
        1740             :
        1741             :       return
        1742             :
        1743             :    end if ! (ipass_calc_updraft == 2)
        1744             :
        1745             : ! check f_ent > 0
        1746             :    if (f_ent <= 0.0_r8) return
        1747             :
        1748             :    hygro = 0.0_r8
        1749             :    vaerosol = 0.0_r8
        1750             :    naerosol = 0.0_r8
        1751             :
        1752             :    do m = 1, nbins
        1753             : ! compute a (or a+cw) volume and hygroscopicity
        1754             :       tmpa = 0.0_r8
        1755             :       tmpb = 0.0_r8
        1756             :       do l = 1, aero_props%nmasses(m)
        1757             :
        1758             :          mm = aero_props%indexer(m,l)
        1759             :
        1760             :          call aero_props%get(m, l, spectype=spec_type, density=spec_dens, hygro=spec_hygro)
        1761             :
        1762             :          tmpc = max( conent(1,mm), 0.0_r8 )
        1763             :          if ( use_cwaer_for_activate_maxsat ) &
        1764             :          tmpc = tmpc + max( conent(2,mm), 0.0_r8 )
        1765             :          tmpc = tmpc / spec_dens
        1766             :          tmpa = tmpa + tmpc
        1767             :          tmpb = tmpb + tmpc * spec_hygro
        1768             :       end do
        1769             :       vaerosol(m) = tmpa * rhoair
        1770             :       if (tmpa < 1.0e-35_r8) then
        1771             :          hygro(m) = 0.2_r8
        1772             :       else
        1773             :          hygro(m) = tmpb/tmpa
        1774             :       end if
        1775             :
        1776             : ! load a (or a+cw) number and bound it
        1777             :       tmpa = max( conent(1,mm), 0.0_r8 )
        1778             :       if ( use_cwaer_for_activate_maxsat ) &
        1779             :       tmpa = tmpa + max( conent(2,mm), 0.0_r8 )
        1780             :       naerosol(m) = tmpa * rhoair
        1781             :
        1782             :       naerosol_a(1) = naerosol(m)
        1783             :       vaerosol_a(1) = vaerosol(m)
        1784             :
        1785             :       call aero_props%apply_number_limits( naerosol_a, vaerosol_a, 1, 1, m )
        1786             :
        1787             :       naerosol(m) = naerosol_a(1)
        1788             :    end do
        1789             :
        1790             : ! call Razzak-Ghan activation routine with single updraft
        1791             :    wbar = max( wup, 0.5_r8 )  ! force wbar >= 0.5 m/s for now
        1792             :    sigw = 0.0_r8
        1793             :    wdiab = 0.0_r8
        1794             :    wminf = wbar
        1795             :    wmaxf = wbar
        1796             :
        1797             :    call activate_aerosol(                                                    &
        1798             :          wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,                    &
        1799             :          naerosol, nbins, vaerosol, hygro, aero_props,            &
        1800             :          fn, fm, fluxn, fluxm, flux_fullact                                )
        1801             :
        1802             : ! apply the activation fractions to the updraft aerosol mixing ratios
        1803             :    dt_u_inv = 1.0_r8/dt_u
        1804             :
        1805             :    do m = 1, aero_props%nbins()
        1806             :       do l = 0, aero_props%nmasses(m)
        1807             :          mm = aero_props%indexer(m,l)
        1808             :
        1809             :          if ( (method_reduce_actfrac == 1)      .and. &
        1810             :               (factor_reduce_actfrac >= 0.0_r8) .and. &
        1811             :               (factor_reduce_actfrac <  1.0_r8) )     &
        1812             :               tmp_fact = tmp_fact * factor_reduce_actfrac
        1813             :
        1814             :          delact = min( conent(1,mm)*tmp_fact*f_ent, conu(1,mm) )
        1815             :          delact = max( delact, 0.0_r8 )
        1816             :          conu(1,mm) = conu(1,mm) - delact
        1817             :          conu(2,mm) = conu(2,mm) + delact
        1818             :          dconudt(1,mm) = -delact*dt_u_inv
        1819             :          dconudt(2,mm) =  delact*dt_u_inv
        1820             :       end do
        1821             :    end do
        1822             :
        1823             :    end subroutine activate_convproc
        1824             :
        1825             : !=========================================================================================
        1826           0 :    subroutine activate_convproc_method2( aero_props, &
        1827           0 :               conu,       dconudt,             &
        1828             :               f_ent,      dt_u,     wup,       &
        1829             :               tair,       rhoair,   k,         &
        1830             :               kactfirst,  ipass_calc_updraft    )
        1831             : !-----------------------------------------------------------------------
        1832             : !
        1833             : ! Purpose:
        1834             : ! Calculate activation of aerosol species in convective updraft
        1835             : ! for a single column and level
        1836             : !
        1837             : ! Method:
        1838             : ! conu(l)    = Updraft TMR (tracer mixing ratio) at k/k-1 interface
        1839             : ! f_ent      = Fraction of the "before-detrainment" updraft massflux at
        1840             : !              k/k-1 interface" resulting from entrainment of level k air
        1841             : !              (where k is the current level in subr aero_convproc_tend)
        1842             : !
        1843             : ! On entry to this routine, the conu(l) represents the updraft TMR
        1844             : ! after entrainment, but before chemistry/physics and detrainment.
        1845             : !
        1846             : ! This routine applies aerosol activation to the conu tracer mixing ratios,
        1847             : ! then adjusts the conu so that on exit,
        1848             : !   conu(la) = conu_incoming(la) - conu(la)*f_act(la)
        1849             : !   conu(lc) = conu_incoming(lc) + conu(la)*f_act(la)
        1850             : ! where
        1851             : !   la, lc   = indices for an unactivated/activated aerosol component pair
        1852             : !   f_act    = fraction of conu(la) that is activated.  The f_act are
        1853             : !              calculated with the Razzak-Ghan activation parameterization.
        1854             : !              The f_act differ for each mode, and for number/surface/mass.
        1855             : !
        1856             : ! At cloud base (k==kactfirst), primary activation is done using the
        1857             : ! "standard" code in subr activate do diagnose maximum supersaturation.
        1858             : ! Above cloud base, secondary activation is done using a
        1859             : ! prescribed supersaturation.
        1860             : !
        1861             : ! *** The updraft velocity used for activation calculations is rather
        1862             : !     uncertain and needs more work.  However, an updraft of 1-3 m/s
        1863             : !     will activate essentially all of accumulation and coarse mode particles.
        1864             : !
        1865             : ! Author: R. Easter
        1866             : !
        1867             : !-----------------------------------------------------------------------
        1868             :
        1869             :    use ndrop, only: activate_aerosol
        1870             :
        1871             : !-----------------------------------------------------------------------
        1872             : ! arguments  (note:  TMR = tracer mixing ratio)
        1873             :
        1874             :    class(aerosol_properties), intent(in) :: aero_props
        1875             :
        1876             :    ! conu = tracer mixing ratios in updraft at top of this (current) level
        1877             :    !        The conu are changed by activation
        1878             :    real(r8), intent(inout) :: conu(2,ncnstaer)
        1879             :    real(r8), intent(inout) :: dconudt(2,ncnstaer) ! TMR tendencies due to activation
        1880             :
        1881             :    real(r8), intent(in)    :: f_ent  ! fraction of updraft massflux that was
        1882             :                                      ! entrained across this layer == eudp/mu_p_eudp
        1883             :    real(r8), intent(in)    :: dt_u   ! lagrangian transport time (s) in the
        1884             :                                      ! updraft at current level
        1885             :    real(r8), intent(in)    :: wup    ! mean updraft vertical velocity (m/s)
        1886             :                                      ! at current level updraft
        1887             :
        1888             :    real(r8), intent(in)    :: tair   ! Temperature in Kelvin
        1889             :    real(r8), intent(in)    :: rhoair ! air density (kg/m3)
        1890             :                                      ! used as in-cloud wet removal rate
        1891             :    integer,  intent(in)    :: k      ! level index
        1892             :    integer,  intent(in)    :: kactfirst ! k at cloud base
        1893             :    integer,  intent(in)    :: ipass_calc_updraft
        1894             :
        1895             : !-----------------------------------------------------------------------
        1896             : ! local variables
        1897             :    integer  :: l, m, mm
        1898             :
        1899             :    real(r8) :: delact           ! working variable
        1900             :    real(r8) :: dt_u_inv         ! 1.0/dt_u
        1901           0 :    real(r8) :: fluxm(nbins)     ! to understand this, see subr activate_aerosol
        1902           0 :    real(r8) :: fluxn(nbins)     ! to understand this, see subr activate_aerosol
        1903             :    real(r8) :: flux_fullact     ! to understand this, see subr activate_aerosol
        1904           0 :    real(r8) :: fm(nbins)        ! mass fraction of aerosols activated
        1905           0 :    real(r8) :: fn(nbins)        ! number fraction of aerosols activated
        1906           0 :    real(r8) :: hygro(nbins)     ! current hygroscopicity for int+act
        1907           0 :    real(r8) :: naerosol(nbins)  ! interstitial+activated number conc (#/m3)
        1908             :    real(r8) :: sigw             ! standard deviation of updraft velocity (cm/s)
        1909             :    real(r8) :: smax_prescribed  ! prescribed supersaturation for secondary activation (0-1 fraction)
        1910             :    real(r8) :: tmp_fact         ! working variable
        1911           0 :    real(r8) :: vaerosol(nbins)  ! int+act volume (m3/m3)
        1912             :    real(r8) :: wbar             ! mean updraft velocity (cm/s)
        1913             :    real(r8) :: wdiab            ! diabatic vertical velocity (cm/s)
        1914             :    real(r8) :: wminf, wmaxf     ! limits for integration over updraft spectrum (cm/s)
        1915             :
        1916             :    real(r8) :: spec_hygro
        1917             :    real(r8) :: spec_dens
        1918             :    character(len=32) :: spec_type
        1919             :
        1920             :    real(r8) :: tmpa, tmpb, tmpc ! working variable
        1921             :    real(r8) :: naerosol_a(1)    ! number conc (1/m3)
        1922             :    real(r8) :: vaerosol_a(1)    ! volume conc (m3/m3)
        1923             :
        1924             : !-----------------------------------------------------------------------
        1925             :
        1926             : ! when ipass_calc_updraft == 2, apply the activation tendencies
        1927             : !    from pass 1, but multiplied by factor_reduce_actfrac
        1928             : ! (can only have ipass_calc_updraft == 2 when method_reduce_actfrac = 2)
        1929             :
        1930           0 :    if (ipass_calc_updraft == 2) then
        1931             :
        1932           0 :       dt_u_inv = 1.0_r8/dt_u
        1933             :
        1934           0 :       do m = 1, aero_props%nbins()
        1935           0 :          do l = 0, aero_props%nmasses(m)
        1936           0 :             mm = aero_props%indexer(m,l)
        1937             :
        1938           0 :             delact = dconudt(2,mm)*dt_u * factor_reduce_actfrac
        1939           0 :             delact = min( delact, conu(1,mm) )
        1940           0 :             delact = max( delact, 0.0_r8 )
        1941           0 :             conu(1,mm) = conu(1,mm) - delact
        1942           0 :             conu(2,mm) = conu(2,mm) + delact
        1943           0 :             dconudt(1,mm) = -delact*dt_u_inv
        1944           0 :             dconudt(2,mm) =  delact*dt_u_inv
        1945             :          end do
        1946             :       end do   ! "n = 1, ntot_amode"
        1947             :       return
        1948             :
        1949             :    end if ! (ipass_calc_updraft == 2)
        1950             :
        1951             : ! check f_ent > 0
        1952           0 :    if (f_ent <= 0.0_r8) return
        1953             :
        1954           0 :    hygro = 0.0_r8
        1955           0 :    vaerosol = 0.0_r8
        1956           0 :    naerosol = 0.0_r8
        1957             :
        1958           0 :    do m = 1, nbins
        1959             : ! compute a (or a+cw) volume and hygroscopicity
        1960           0 :       tmpa = 0.0_r8
        1961           0 :       tmpb = 0.0_r8
        1962           0 :       do l = 1, aero_props%nspecies(m)
        1963             :
        1964           0 :          mm = aero_props%indexer(m,l)
        1965             :
        1966           0 :          call aero_props%get(m, l, spectype=spec_type, density=spec_dens, hygro=spec_hygro)
        1967             :
        1968           0 :          tmpc = max( conu(1,mm), 0.0_r8 )
        1969             :          if ( use_cwaer_for_activate_maxsat ) &
        1970             :          tmpc = tmpc + max( conu(2,mm), 0.0_r8 )
        1971           0 :          tmpc = tmpc / spec_dens
        1972           0 :          tmpa = tmpa + tmpc
        1973             :
        1974             :          ! Change the hygroscopicity of POM based on the discussion with Prof.
        1975             :          ! Xiaohong Liu. Some observational studies found that the primary organic
        1976             :          ! material from biomass burning emission shows very high hygroscopicity.
        1977             :          ! Also, found that BC mass will be overestimated if all the aerosols in
        1978             :          ! the primary mode are free to be removed. Therefore, set the hygroscopicity
        1979             :          ! of POM here as 0.2 to enhance the wet scavenge of primary BC and POM.
        1980             :
        1981           0 :          if (spec_type=='p-organic' .and. convproc_pom_spechygro>0._r8) then
        1982           0 :             tmpb = tmpb + tmpc * convproc_pom_spechygro
        1983             :          else
        1984           0 :             tmpb = tmpb + tmpc * spec_hygro
        1985             :          end if
        1986             :       end do
        1987           0 :       vaerosol(m) = tmpa * rhoair
        1988           0 :       if (tmpa < 1.0e-35_r8) then
        1989           0 :          hygro(m) = 0.2_r8
        1990             :       else
        1991           0 :          hygro(m) = tmpb/tmpa
        1992             :       end if
        1993             :
        1994           0 :       mm = aero_props%indexer(m,0)
        1995             :
        1996             : ! load a (or a+cw) number and bound it
        1997           0 :       tmpa = max( conu(1,mm), 0.0_r8 )
        1998             :       if ( use_cwaer_for_activate_maxsat ) &
        1999             :       tmpa = tmpa + max( conu(2,mm), 0.0_r8 )
        2000           0 :       naerosol(m) = tmpa * rhoair
        2001             :
        2002           0 :       naerosol_a(1) = naerosol(m)
        2003           0 :       vaerosol_a(1) = vaerosol(m)
        2004             :
        2005           0 :       call aero_props%apply_number_limits( naerosol_a, vaerosol_a, 1, 1, m )
        2006             :
        2007           0 :       naerosol(m) = naerosol_a(1)
        2008             :
        2009             :    end do
        2010             :
        2011             : ! call Razzak-Ghan activation routine with single updraft
        2012           0 :    wbar = max( wup, 0.5_r8 )  ! force wbar >= 0.5 m/s for now
        2013           0 :    sigw = 0.0_r8
        2014           0 :    wdiab = 0.0_r8
        2015           0 :    wminf = wbar
        2016           0 :    wmaxf = wbar
        2017             :
        2018           0 :    if (k == kactfirst) then
        2019             :
        2020             :       call activate_aerosol(                                                 &
        2021             :          wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,                    &
        2022             :          naerosol, nbins, vaerosol, hygro, aero_props,            &
        2023           0 :          fn, fm, fluxn, fluxm, flux_fullact                                )
        2024             :
        2025             :
        2026             :    else
        2027             : ! above cloud base - do secondary activation with prescribed supersat
        2028             : ! that is constant with height
        2029           0 :       smax_prescribed = method2_activate_smaxmax
        2030             :       call activate_aerosol(                                                 &
        2031             :          wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,                    &
        2032             :          naerosol, nbins, vaerosol, hygro, aero_props,            &
        2033           0 :          fn, fm, fluxn, fluxm, flux_fullact, smax_prescribed               )
        2034             :    end if
        2035             :
        2036             : ! apply the activation fractions to the updraft aerosol mixing ratios
        2037           0 :    dt_u_inv = 1.0_r8/dt_u
        2038             :
        2039           0 :    do m = 1, aero_props%nbins()
        2040           0 :       do l = 0, aero_props%nmasses(m)
        2041           0 :          mm = aero_props%indexer(m,l)
        2042           0 :          if (l==0) then
        2043           0 :             tmp_fact = fn(m)
        2044             :          else
        2045           0 :             tmp_fact = fm(m)
        2046             :          end if
        2047             :
        2048             :          if ( (method_reduce_actfrac == 1)      .and. &
        2049             :               (factor_reduce_actfrac >= 0.0_r8) .and. &
        2050             :               (factor_reduce_actfrac <  1.0_r8) )     &
        2051             :               tmp_fact = tmp_fact * factor_reduce_actfrac
        2052             :
        2053           0 :          delact = min( conu(1,mm)*tmp_fact, conu(1,mm) )
        2054           0 :          delact = max( delact, 0.0_r8 )
        2055           0 :          conu(1,mm) = conu(1,mm) - delact
        2056           0 :          conu(2,mm) = conu(2,mm) + delact
        2057           0 :          dconudt(1,mm) = -delact*dt_u_inv
        2058           0 :          dconudt(2,mm) =  delact*dt_u_inv
        2059             :       end do
        2060             :    end do
        2061             :
        2062           0 :    end subroutine activate_convproc_method2
        2063             :
        2064             : !=========================================================================================
        2065           0 :    subroutine resuspend_convproc( aero_props, &
        2066           0 :               dcondt,  dcondt_resusp, ktop,  kbot_prevap )
        2067             : !-----------------------------------------------------------------------
        2068             : !
        2069             : ! Purpose:
        2070             : ! Calculate resuspension of activated aerosol species resulting from both
        2071             : !    detrainment from updraft and downdraft into environment
        2072             : !    subsidence and lifting of environment, which may move air from
        2073             : !       levels with large-scale cloud to levels with no large-scale cloud
        2074             : !
        2075             : ! Method:
        2076             : ! Three possible approaches were considered:
        2077             : !
        2078             : ! 1. Ad-hoc #1 approach.  At each level, adjust dcondt for the activated
        2079             : !    and unactivated portions of a particular aerosol species so that the
        2080             : !    ratio of dcondt (activated/unactivate) is equal to the ratio of the
        2081             : !    mixing ratios before convection.
        2082             : !    THIS WAS IMPLEMENTED IN MIRAGE2
        2083             : !
        2084             : ! 2. Ad-hoc #2 approach.  At each level, adjust dcondt for the activated
        2085             : !    and unactivated portions of a particular aerosol species so that the
        2086             : !    change to the activated portion is minimized (zero if possible).  The
        2087             : !    would minimize effects of convection on the large-scale cloud.
        2088             : !    THIS IS CURRENTLY IMPLEMENTED IN CAM5 where we assume that convective
        2089             : !    clouds have no impact on the stratiform-cloudborne aerosol
        2090             : !
        2091             : ! 3. Mechanistic approach that treats the details of interactions between
        2092             : !    the large-scale and convective clouds.  (Something for the future.)
        2093             : !
        2094             : ! Author: R. Easter
        2095             : !
        2096             : !-----------------------------------------------------------------------
        2097             :
        2098             : !-----------------------------------------------------------------------
        2099             : ! arguments
        2100             : ! (note:  TMR = tracer mixing ratio)
        2101             :
        2102             :    class(aerosol_properties), intent(in) :: aero_props
        2103             :    real(r8), intent(inout) :: dcondt(2,ncnstaer,pver)
        2104             :                               ! overall TMR tendency from convection
        2105             :    real(r8), intent(inout) :: dcondt_resusp(2,ncnstaer,pver)
        2106             :                               ! portion of TMR tendency due to resuspension
        2107             :                               ! (actually, due to the adjustments made here)
        2108             :    integer,  intent(in)    :: ktop, kbot_prevap ! indices of top and bottom cloud levels
        2109             :
        2110             : !-----------------------------------------------------------------------
        2111             : ! local variables
        2112             :    integer  :: k, l, m, mm
        2113             :    real(r8) :: qdota, qdotc, qdotac  ! working variables (MR tendencies)
        2114             :    !-----------------------------------------------------------------------
        2115             :
        2116             :    ! apply adjustments to dcondt for pairs of unactivated and
        2117             :    ! activated aerosol species
        2118           0 :    do m = 1, aero_props%nbins()
        2119           0 :       do l = 0, aero_props%nmasses(m)
        2120           0 :          mm = aero_props%indexer(m,l)
        2121             :
        2122           0 :          do k = ktop, kbot_prevap
        2123           0 :             if (convproc_do_evaprain_atonce) then
        2124           0 :                dcondt_resusp(1,mm,k) = dcondt(1,mm,k)
        2125           0 :                dcondt_resusp(2,mm,k) = dcondt(2,mm,k)
        2126             :             else
        2127           0 :                qdota = dcondt(1,mm,k)
        2128           0 :                qdotc = dcondt(2,mm,k)
        2129           0 :                qdotac = qdota + qdotc
        2130             :
        2131           0 :                dcondt(1,mm,k) = qdotac
        2132           0 :                dcondt(2,mm,k) = 0.0_r8
        2133             :
        2134           0 :                dcondt_resusp(1,mm,k) = (dcondt(1,mm,k) - qdota)
        2135           0 :                dcondt_resusp(2,mm,k) = (dcondt(2,mm,k) - qdotc)
        2136             :             end if
        2137             :          end do
        2138             :
        2139             :       end do
        2140             :    end do
        2141             :
        2142           0 :    end subroutine resuspend_convproc
        2143             :
        2144             : !=========================================================================================
        2145             :
        2146             : end module aero_convproc

Generated by: LCOV version 1.14