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