LCOV - code coverage report
Current view: | top level - physics/cam - ndrop.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 638 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 5 | 0.0 %

              Line data    Source code
           1             : module ndrop
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : ! Purpose:
           5             : !   CAM Interface for droplet activation by modal aerosols
           6             : !
           7             : ! ***N.B.*** This module is currently hardcoded to recognize only the modes that
           8             : !            affect the climate calculation.  This is implemented by using list
           9             : !            index 0 in all the calls to rad_constituent interfaces.
          10             : !---------------------------------------------------------------------------------
          11             :
          12             : use shr_kind_mod,     only: r8 => shr_kind_r8, shr_kind_cs
          13             : use ppgrid,           only: pcols, pver
          14             : use physconst,        only: pi, rhoh2o, mwh2o, r_universal, rh2o, &
          15             :                             gravit, latvap, cpair, rair
          16             : use constituents,     only: pcnst, cnst_get_ind, cnst_name, cnst_spec_class_gas, cnst_species_class
          17             : use physics_types,    only: physics_state, physics_ptend, physics_ptend_init
          18             : use physics_buffer,   only: physics_buffer_desc, pbuf_get_index, pbuf_get_field
          19             :
          20             : use wv_saturation,    only: qsat
          21             : use phys_control,     only: phys_getopts
          22             : use ref_pres,         only: top_lev => trop_cloud_top_lev
          23             : use shr_spfn_mod,     only: erf => shr_spfn_erf
          24             : use cam_history,      only: addfld, add_default, horiz_only, fieldname_len, outfld
          25             : use cam_abortutils,   only: endrun
          26             : use cam_logfile,      only: iulog
          27             :
          28             : use aerosol_properties_mod, only: aerosol_properties
          29             : use aerosol_state_mod, only: aerosol_state, ptr2d_t
          30             :
          31             : implicit none
          32             : private
          33             : save
          34             :
          35             : public ndrop_init, dropmixnuc, activate_aerosol
          36             :
          37             : ! mathematical constants
          38             : real(r8), parameter :: zero     = 0._r8
          39             : real(r8), parameter :: third    = 1._r8/3._r8
          40             : real(r8), parameter :: twothird = 2._r8*third
          41             : real(r8), parameter :: sixth    = 1._r8/6._r8
          42             : real(r8), parameter :: sq2      = sqrt(2._r8)
          43             : real(r8), parameter :: sq2pi    = sqrt(2._r8*pi)
          44             : real(r8), parameter :: sqpi     = sqrt(pi)
          45             : real(r8), parameter :: surften  = 0.076_r8
          46             : real(r8), parameter :: tmelt    = 273._r8
          47             :
          48             : real(r8) :: aten
          49             :
          50             : ! CCN diagnostic fields
          51             : integer,  parameter :: psat=6    ! number of supersaturations to calc ccn concentration
          52             : real(r8), parameter :: supersat(psat)= & ! supersaturation (%) to determine ccn concentration
          53             :                        (/ 0.02_r8, 0.05_r8, 0.1_r8, 0.2_r8, 0.5_r8, 1.0_r8 /)
          54             : character(len=8) :: ccn_name(psat)= &
          55             :                     (/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
          56             :
          57             : ! indices in state and pbuf structures
          58             : integer :: numliq_idx = -1
          59             : integer :: kvh_idx    = -1
          60             :
          61             : logical :: history_aerosol      ! Output the aerosol tendencies
          62             : character(len=fieldname_len), allocatable :: fieldname(:)    ! names for drop nuc tendency output fields
          63             : character(len=fieldname_len), allocatable :: fieldname_cw(:) ! names for drop nuc tendency output fields
          64             :
          65             : ! Indices for aerosol species in the ptend%q array.
          66             : integer, allocatable :: aer_cnst_idx(:,:)
          67             :
          68             : logical :: lq(pcnst) = .false. ! set flags true for constituents with non-zero tendencies
          69             :                                ! in the ptend object
          70             :
          71             : !===============================================================================
          72             : contains
          73             : !===============================================================================
          74             :
          75           0 : subroutine ndrop_init(aero_props)
          76             :
          77             :    class(aerosol_properties), intent(in) :: aero_props
          78             :
          79             :    integer :: l, m, mm
          80             :    integer :: idxtmp = -1
          81             :    character(len=32)   :: tmpname
          82             :    character(len=32)   :: tmpname_cw
          83             :    character(len=128)  :: long_name
          84             :    character(len=8)    :: unit
          85             :    logical :: history_amwg         ! output the variables used by the AMWG diag package
          86             :
          87             :    !-------------------------------------------------------------------------------
          88             :
          89             :    ! get indices into state%q and pbuf structures
          90           0 :    call cnst_get_ind('NUMLIQ', numliq_idx)
          91             :
          92           0 :    kvh_idx = pbuf_get_index('kvh')
          93             :
          94           0 :    aten = 2._r8*mwh2o*surften/(r_universal*tmelt*rhoh2o)
          95             :
          96             :    allocate( &
          97           0 :       aer_cnst_idx(aero_props%nbins(),0:maxval(aero_props%nmasses())), &
          98           0 :       fieldname(aero_props%ncnst_tot()),                 &
          99           0 :       fieldname_cw(aero_props%ncnst_tot())               )
         100             :
         101             :    ! Add dropmixnuc tendencies for all modal aerosol species
         102             :
         103             :    call phys_getopts(history_amwg_out = history_amwg, &
         104           0 :                      history_aerosol_out = history_aerosol)
         105             :
         106             :
         107           0 :    do m = 1, aero_props%nbins()
         108           0 :       do l = 0, aero_props%nmasses(m)
         109             :
         110           0 :          mm = aero_props%indexer(m,l)
         111             :
         112           0 :          unit = 'kg/m2/s'
         113           0 :          if (l == 0) then   ! number
         114           0 :             unit = '#/m2/s'
         115             :          end if
         116             :
         117           0 :          if (l == 0) then   ! number
         118           0 :             call aero_props%num_names( m, tmpname, tmpname_cw)
         119             :          else
         120           0 :             call aero_props%mmr_names( m,l, tmpname, tmpname_cw)
         121             :          end if
         122             :
         123           0 :          fieldname(mm)    = trim(tmpname) // '_mixnuc1'
         124           0 :          fieldname_cw(mm) = trim(tmpname_cw) // '_mixnuc1'
         125             :
         126             :          ! To set tendencies in the ptend object need to get the constituent indices
         127             :          ! for the prognostic species
         128             :
         129           0 :          call cnst_get_ind(tmpname, idxtmp, abort=.false.)
         130           0 :          aer_cnst_idx(m,l) = idxtmp
         131             :
         132           0 :          if (idxtmp>0) then
         133           0 :             lq(idxtmp) = .true.
         134             :          end if
         135             :
         136             :          ! Add tendency fields to the history only when prognostic MAM is enabled.
         137           0 :          long_name = trim(tmpname) // ' dropmixnuc mixnuc column tendency'
         138           0 :          call addfld(fieldname(mm),    horiz_only, 'A', unit, long_name)
         139             :
         140           0 :          long_name = trim(tmpname_cw) // ' dropmixnuc mixnuc column tendency'
         141           0 :          call addfld(fieldname_cw(mm), horiz_only, 'A', unit, long_name)
         142             :
         143           0 :          if (history_aerosol) then
         144           0 :             call add_default(fieldname(mm), 1, ' ')
         145           0 :             call add_default(fieldname_cw(mm), 1, ' ')
         146             :          end if
         147             :
         148             :       end do
         149             :    end do
         150             :
         151           0 :    call addfld('CCN1',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.02%')
         152           0 :    call addfld('CCN2',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.05%')
         153           0 :    call addfld('CCN3',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.1%')
         154           0 :    call addfld('CCN4',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.2%')
         155           0 :    call addfld('CCN5',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.5%')
         156           0 :    call addfld('CCN6',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=1.0%')
         157             :
         158             :
         159           0 :    call addfld('WTKE',     (/ 'lev' /), 'A', 'm/s', 'Standard deviation of updraft velocity')
         160           0 :    call addfld('NDROPMIX', (/ 'lev' /), 'A', '#/kg/s', 'Droplet number mixing')
         161           0 :    call addfld('NDROPSRC', (/ 'lev' /), 'A', '#/kg/s', 'Droplet number source')
         162           0 :    call addfld('NDROPSNK', (/ 'lev' /), 'A', '#/kg/s', 'Droplet number loss by microphysics')
         163           0 :    call addfld('NDROPCOL', horiz_only,  'A', '#/m2', 'Column droplet number')
         164             :
         165             :    ! set the add_default fields
         166           0 :    if (history_amwg) then
         167           0 :       call add_default('CCN3', 1, ' ')
         168             :    endif
         169             :
         170           0 : end subroutine ndrop_init
         171             :
         172             : !===============================================================================
         173             :
         174           0 : subroutine dropmixnuc( aero_props, aero_state, &
         175             :    state, ptend, dtmicro, pbuf, wsub, wmixmin, &
         176           0 :    cldn, cldo, cldliqf, tendnd, factnum, from_spcam)
         177             :
         178             :    ! vertical diffusion and nucleation of cloud droplets
         179             :    ! assume cloud presence controlled by cloud fraction
         180             :    ! doesn't distinguish between warm, cold clouds
         181             :
         182             :    ! arguments
         183             :    type(physics_state), target, intent(in)    :: state
         184             :    type(physics_ptend),         intent(out)   :: ptend
         185             :    real(r8),                    intent(in)    :: dtmicro     ! time step for microphysics (s)
         186             :    real(r8),                    intent(in)    :: wmixmin     ! minimum turbulence vertical velocity (m/s)
         187             :
         188             :    type(physics_buffer_desc), pointer :: pbuf(:)
         189             :
         190             :    class(aerosol_properties), intent(in) :: aero_props
         191             :    class(aerosol_state), intent(in) :: aero_state
         192             :
         193             :    ! arguments
         194             :    real(r8), intent(in) :: wsub(pcols,pver)    ! subgrid vertical velocity
         195             :    real(r8), intent(in) :: cldn(pcols,pver)    ! cloud fraction
         196             :    real(r8), intent(in) :: cldo(pcols,pver)    ! cloud fraction on previous time step
         197             :    real(r8), intent(in) :: cldliqf(pcols,pver) ! liquid cloud fraction (liquid / (liquid + ice))
         198             :    logical,  intent(in),optional :: from_spcam ! value insignificant - if variable present, is called from spcam
         199             :
         200             :    ! output arguments
         201             :    real(r8), intent(out) :: tendnd(pcols,pver) ! change in droplet number concentration (#/kg/s)
         202             :    real(r8), intent(out) :: factnum(:,:,:)     ! activation fraction for aerosol number
         203             :    !--------------------Local storage-------------------------------------
         204             :
         205             :    integer  :: lchnk               ! chunk identifier
         206             :    integer  :: ncol                ! number of columns
         207             :    integer  :: nbin                ! number of modes/bins
         208             :    integer  :: nele_tot            ! total number of aerosol elements
         209             :
         210           0 :    real(r8), pointer :: ncldwtr(:,:) ! droplet number concentration (#/kg)
         211           0 :    real(r8), pointer :: temp(:,:)    ! temperature (K)
         212           0 :    real(r8), pointer :: pmid(:,:)    ! mid-level pressure (Pa)
         213           0 :    real(r8), pointer :: pint(:,:)    ! pressure at layer interfaces (Pa)
         214           0 :    real(r8), pointer :: pdel(:,:)    ! pressure thickess of layer (Pa)
         215           0 :    real(r8), pointer :: rpdel(:,:)   ! inverse of pressure thickess of layer (/Pa)
         216           0 :    real(r8), pointer :: zm(:,:)      ! geopotential height of level (m)
         217             :
         218           0 :    real(r8), pointer :: kvh(:,:)     ! vertical diffusivity (m2/s)
         219             :
         220           0 :    type(ptr2d_t), allocatable :: raer(:)     ! aerosol mass, number mixing ratios
         221           0 :    type(ptr2d_t), allocatable :: qqcw(:)
         222             :    real(r8) :: raertend(pver)  ! tendency of aerosol mass, number mixing ratios
         223             :    real(r8) :: qqcwtend(pver)  ! tendency of cloudborne aerosol mass, number mixing ratios
         224             :
         225             :    real(r8), parameter :: zkmin = 0.01_r8, zkmax = 100._r8
         226             :    integer  :: i, k, l, m, mm, n
         227             :    integer  :: km1, kp1
         228             :    integer  :: nnew, nsav, ntemp
         229             :    integer  :: lptr
         230             :    integer  :: nsubmix, nsubmix_bnd
         231             :    integer, save :: count_submix(100)
         232             :    integer  :: phase ! phase of aerosol
         233             :
         234             :    real(r8) :: arg
         235             :    real(r8) :: dtinv
         236             :    real(r8) :: dtmin, tinv, dtt
         237             :    real(r8) :: lcldn(pcols,pver)
         238             :    real(r8) :: lcldo(pcols,pver)
         239             :
         240             :    real(r8) :: zs(pver) ! inverse of distance between levels (m)
         241             :    real(r8) :: qcld(pver) ! cloud droplet number mixing ratio (#/kg)
         242             :    real(r8) :: qncld(pver)     ! droplet number nucleated on cloud boundaries
         243             :    real(r8) :: srcn(pver)       ! droplet source rate (/s)
         244             :    real(r8) :: cs(pcols,pver)      ! air density (kg/m3)
         245             :    real(r8) :: csbot(pver)       ! air density at bottom (interface) of layer (kg/m3)
         246             :    real(r8) :: csbot_cscen(pver) ! csbot(i)/cs(i,k)
         247             :    real(r8) :: dz(pcols,pver)      ! geometric thickness of layers (m)
         248             :
         249             :    real(r8) :: wtke(pcols,pver)     ! turbulent vertical velocity at base of layer k (m/s)
         250             :    real(r8) :: wtke_cen(pcols,pver) ! turbulent vertical velocity at center of layer k (m/s)
         251             :    real(r8) :: wbar, wmix, wmin, wmax
         252             :
         253             :    real(r8) :: zn(pver)   ! g/pdel (m2/g) for layer
         254             :    real(r8) :: flxconv    ! convergence of flux into lowest layer
         255             :
         256             :    real(r8) :: wdiab           ! diabatic vertical velocity
         257             :    real(r8) :: ekd(pver)       ! diffusivity for droplets (m2/s)
         258             :    real(r8) :: ekk(0:pver)     ! density*diffusivity for droplets (kg/m3 m2/s)
         259             :    real(r8) :: ekkp(pver)      ! zn*zs*density*diffusivity
         260             :    real(r8) :: ekkm(pver)      ! zn*zs*density*diffusivity
         261             :
         262             :    real(r8) :: dum, dumc
         263             :    real(r8) :: tmpa
         264             :    real(r8) :: dact
         265             :    real(r8) :: fluxntot         ! (#/cm2/s)
         266             :    real(r8) :: dtmix
         267             :    real(r8) :: alogarg
         268             :    real(r8) :: overlapp(pver), overlapm(pver) ! cloud overlap
         269             :
         270             :    real(r8) :: nsource(pcols,pver)            ! droplet number source (#/kg/s)
         271             :    real(r8) :: ndropmix(pcols,pver)           ! droplet number mixing (#/kg/s)
         272             :    real(r8) :: ndropcol(pcols)               ! column droplet number (#/m2)
         273             :    real(r8) :: cldo_tmp, cldn_tmp
         274             :    real(r8) :: tau_cld_regenerate
         275             :    real(r8) :: taumix_internal_pver_inv ! 1/(internal mixing time scale for k=pver) (1/s)
         276             :
         277             :
         278           0 :    real(r8), allocatable :: nact(:,:)  ! fractional aero. number  activation rate (/s)
         279           0 :    real(r8), allocatable :: mact(:,:)  ! fractional aero. mass    activation rate (/s)
         280             :
         281           0 :    real(r8), allocatable :: raercol(:,:,:)    ! single column of aerosol mass, number mixing ratios
         282           0 :    real(r8), allocatable :: raercol_cw(:,:,:) ! same as raercol but for cloud-borne phase
         283             :
         284             :
         285             :    real(r8) :: na(pcols), va(pcols), hy(pcols)
         286           0 :    real(r8), allocatable :: naermod(:)  ! (1/m3)
         287           0 :    real(r8), allocatable :: hygro(:)    ! hygroscopicity of aerosol mode
         288           0 :    real(r8), allocatable :: vaerosol(:) ! interstit+activated aerosol volume conc (cm3/cm3)
         289             :
         290             :    real(r8) :: source(pver)
         291             :
         292           0 :    real(r8), allocatable :: fn(:)              ! activation fraction for aerosol number
         293           0 :    real(r8), allocatable :: fm(:)              ! activation fraction for aerosol mass
         294             :
         295           0 :    real(r8), allocatable :: fluxn(:)           ! number  activation fraction flux (cm/s)
         296           0 :    real(r8), allocatable :: fluxm(:)           ! mass    activation fraction flux (cm/s)
         297             :    real(r8)              :: flux_fullact(pver) ! 100%    activation fraction flux (cm/s)
         298             :    !     note:  activation fraction fluxes are defined as
         299             :    !     fluxn = [flux of activated aero. number into cloud (#/cm2/s)]
         300             :    !           / [aero. number conc. in updraft, just below cloudbase (#/cm3)]
         301             :
         302             :
         303           0 :    real(r8), allocatable :: coltend(:,:)       ! column tendency for diagnostic output
         304           0 :    real(r8), allocatable :: coltend_cw(:,:)    ! column tendency
         305             :    real(r8) :: ccn(pcols,pver,psat)    ! number conc of aerosols activated at supersat
         306             :
         307             :    !for gas species turbulent mixing
         308           0 :    real(r8), pointer :: rgas(:, :, :)
         309           0 :    real(r8), allocatable :: rgascol(:, :, :)
         310           0 :    real(r8), allocatable :: coltendgas(:)
         311             :    real(r8) :: zerogas(pver)
         312             :    character*200 fieldnamegas
         313             :
         314             :    logical  :: called_from_spcam
         315             :    integer :: errnum
         316             :    character(len=shr_kind_cs) :: errstr
         317             :    !-------------------------------------------------------------------------------
         318             :
         319           0 :    lchnk = state%lchnk
         320           0 :    ncol  = state%ncol
         321           0 :    nbin  = aero_props%nbins()
         322           0 :    nele_tot  = aero_props%ncnst_tot()
         323             :
         324           0 :    ncldwtr  => state%q(:,:,numliq_idx)
         325           0 :    temp     => state%t
         326           0 :    pmid     => state%pmid
         327           0 :    pint     => state%pint
         328           0 :    pdel     => state%pdel
         329           0 :    rpdel    => state%rpdel
         330           0 :    zm       => state%zm
         331             :
         332           0 :    call pbuf_get_field(pbuf, kvh_idx, kvh)
         333             :
         334             :    ! Create the liquid weighted cloud fractions that were passsed in
         335             :    ! before. This doesn't seem like the best variable, since the cloud could
         336             :    ! have liquid condensate, but the part of it that is changing could be the
         337             :    ! ice portion; however, this is what was done before.
         338           0 :    lcldo(:ncol,:)  = cldo(:ncol,:)  * cldliqf(:ncol,:)
         339           0 :    lcldn(:ncol,:) = cldn(:ncol,:) * cldliqf(:ncol,:)
         340             :
         341             :
         342           0 :    arg = 1.0_r8
         343           0 :    if (abs(0.8427_r8 - erf(arg))/0.8427_r8 > 0.001_r8) then
         344           0 :       write(iulog,*) 'erf(1.0) = ',ERF(arg)
         345           0 :       call endrun('dropmixnuc: Error function error')
         346             :    endif
         347           0 :    arg = 0.0_r8
         348           0 :    if (erf(arg) /= 0.0_r8) then
         349           0 :       write(iulog,*) 'erf(0.0) = ',erf(arg)
         350           0 :       write(iulog,*) 'dropmixnuc: Error function error'
         351           0 :       call endrun('dropmixnuc: Error function error')
         352             :    endif
         353             :
         354           0 :    dtinv = 1._r8/dtmicro
         355             :
         356             :    allocate( &
         357             :       nact(pver,nbin),             &
         358             :       mact(pver,nbin),             &
         359           0 :       raer(nele_tot),              &
         360           0 :       qqcw(nele_tot),              &
         361             :       raercol(pver,nele_tot,2),    &
         362             :       raercol_cw(pver,nele_tot,2), &
         363             :       coltend(pcols,nele_tot),     &
         364             :       coltend_cw(pcols,nele_tot),  &
         365           0 :       naermod(nbin),               &
         366           0 :       hygro(nbin),                 &
         367           0 :       vaerosol(nbin),              &
         368           0 :       fn(nbin),                    &
         369           0 :       fm(nbin),                    &
         370           0 :       fluxn(nbin),                 &
         371           0 :       fluxm(nbin)               )
         372             :
         373             :    ! Init pointers to mode number and specie mass mixing ratios in
         374             :    ! intersitial and cloud borne phases.
         375           0 :    call aero_state%get_states( aero_props, raer, qqcw )
         376             :
         377           0 :    called_from_spcam = (present(from_spcam))
         378             :
         379           0 :    if (called_from_spcam) then
         380           0 :       rgas  => state%q
         381           0 :       allocate(rgascol(pver, pcnst, 2))
         382           0 :       allocate(coltendgas(pcols))
         383             :    endif
         384             :
         385           0 :    factnum = 0._r8
         386           0 :    wtke = 0._r8
         387           0 :    nsource = 0._r8
         388           0 :    ndropmix = 0._r8
         389           0 :    ndropcol = 0._r8
         390           0 :    tendnd = 0._r8
         391             :
         392             :    ! initialize aerosol tendencies
         393           0 :    call physics_ptend_init(ptend, state%psetcols, 'ndrop', lq=lq)
         394             :
         395             :    ! overall_main_i_loop
         396           0 :    do i = 1, ncol
         397             :
         398           0 :       do k = top_lev, pver-1
         399           0 :          zs(k) = 1._r8/(zm(i,k) - zm(i,k+1))
         400             :       end do
         401           0 :       zs(pver) = zs(pver-1)
         402             :
         403             :       ! load number nucleated into qcld on cloud boundaries
         404             :
         405           0 :       do k = top_lev, pver
         406             :
         407           0 :          qcld(k)  = ncldwtr(i,k)
         408           0 :          qncld(k) = 0._r8
         409           0 :          srcn(k)  = 0._r8
         410           0 :          cs(i,k)  = pmid(i,k)/(rair*temp(i,k))        ! air density (kg/m3)
         411           0 :          dz(i,k)  = 1._r8/(cs(i,k)*gravit*rpdel(i,k)) ! layer thickness in m
         412             :
         413           0 :          do m = 1, nbin
         414           0 :             nact(k,m) = 0._r8
         415           0 :             mact(k,m) = 0._r8
         416             :          end do
         417             :
         418           0 :          zn(k) = gravit*rpdel(i,k)
         419             :
         420           0 :          if (k < pver) then
         421           0 :             ekd(k)   = kvh(i,k+1)
         422           0 :             ekd(k)   = max(ekd(k), zkmin)
         423           0 :             ekd(k)   = min(ekd(k), zkmax)
         424           0 :             csbot(k) = 2.0_r8*pint(i,k+1)/(rair*(temp(i,k) + temp(i,k+1)))
         425           0 :             csbot_cscen(k) = csbot(k)/cs(i,k)
         426             :          else
         427           0 :             ekd(k)   = 0._r8
         428           0 :             csbot(k) = cs(i,k)
         429           0 :             csbot_cscen(k) = 1.0_r8
         430             :          end if
         431             :
         432             :          ! rce-comment - define wtke at layer centers for new-cloud activation
         433             :          !    and at layer boundaries for old-cloud activation
         434           0 :          wtke_cen(i,k) = wsub(i,k)
         435           0 :          wtke(i,k)     = wsub(i,k)
         436           0 :          wtke_cen(i,k) = max(wtke_cen(i,k), wmixmin)
         437           0 :          wtke(i,k)     = max(wtke(i,k), wmixmin)
         438             :
         439           0 :          nsource(i,k) = 0._r8
         440             :
         441             :       end do
         442             :
         443           0 :       nsav = 1
         444           0 :       nnew = 2
         445             :
         446           0 :       do mm = 1,nele_tot
         447           0 :          raercol_cw(:,mm,nsav) = 0.0_r8
         448           0 :          raercol(:,mm,nsav)    = 0.0_r8
         449           0 :          raercol_cw(top_lev:pver,mm,nsav) = qqcw(mm)%fld(i,top_lev:pver)
         450           0 :          raercol(top_lev:pver,mm,nsav)    = raer(mm)%fld(i,top_lev:pver)
         451             :       end do
         452             :
         453           0 :       if (called_from_spcam) then
         454             :       !
         455             :       ! In the MMF model, turbulent mixing for tracer species are turned off.
         456             :       ! So the turbulent for gas species mixing are added here.
         457             :       ! (Previously, it had the turbulent mixing for aerosol species)
         458             :       !
         459           0 :          do m=1, pcnst
         460           0 :             if (cnst_species_class(m) == cnst_spec_class_gas) rgascol(:,m,nsav) = rgas(i,:,m)
         461             :          end do
         462             :
         463             :       endif
         464             :
         465             :       ! droplet nucleation/aerosol activation
         466             :
         467             :       ! tau_cld_regenerate = time scale for regeneration of cloudy air
         468             :       !    by (horizontal) exchange with clear air
         469             :       tau_cld_regenerate = 3600.0_r8 * 3.0_r8
         470             :
         471             :       if (called_from_spcam) then
         472             :       ! when this is called  in the MMF part, no cloud regeneration and decay.
         473             :       ! set the time scale be very long so that no cloud regeneration.
         474             :            tau_cld_regenerate = 3600.0_r8 * 24.0_r8 * 365.0_r8
         475             :       endif
         476             :
         477             :
         478             :       ! k-loop for growing/shrinking cloud calcs .............................
         479             :       ! grow_shrink_main_k_loop: &
         480           0 :       do k = top_lev, pver
         481             :
         482             :          ! This code was designed for liquid clouds, but the cloudbourne
         483             :          ! aerosol can be either from liquid or ice clouds. For the ice clouds,
         484             :          ! we do not do regeneration, but as cloud fraction decreases the
         485             :          ! aerosols should be returned interstitial. The lack of a liquid cloud
         486             :          ! should not mean that all of the aerosol is realease. Therefor a
         487             :          ! section has been added for shrinking ice clouds and checks were added
         488             :          ! to protect ice cloudbourne aerosols from being released when no
         489             :          ! liquid cloud is present.
         490             :
         491             :          ! shrinking ice cloud ......................................................
         492           0 :          cldo_tmp = cldo(i,k) * (1._r8 - cldliqf(i,k))
         493           0 :          cldn_tmp = cldn(i,k) * (1._r8 - cldliqf(i,k))
         494             :
         495           0 :          if (cldn_tmp < cldo_tmp) then
         496             :
         497             :             ! convert activated aerosol to interstitial in decaying cloud
         498             :
         499           0 :             dumc = (cldn_tmp - cldo_tmp)/cldo_tmp * (1._r8 - cldliqf(i,k))
         500           0 :             do mm = 1,nele_tot
         501           0 :                dact   = raercol_cw(k,mm,nsav)*dumc
         502           0 :                raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact   ! cloud-borne aerosol
         503           0 :                raercol(k,mm,nsav)    = raercol(k,mm,nsav) - dact
         504             :             end do
         505             :
         506             :          end if
         507             :
         508             :          ! shrinking liquid cloud ......................................................
         509             :          !    treat the reduction of cloud fraction from when cldn(i,k) < cldo(i,k)
         510             :          !    and also dissipate the portion of the cloud that will be regenerated
         511           0 :          cldo_tmp = lcldo(i,k)
         512           0 :          cldn_tmp = lcldn(i,k) * exp( -dtmicro/tau_cld_regenerate )
         513             :          !    alternate formulation
         514             :          !    cldn_tmp = cldn(i,k) * max( 0.0_r8, (1.0_r8-dtmicro/tau_cld_regenerate) )
         515             :
         516             :          ! fraction is also provided.
         517           0 :          if (cldn_tmp < cldo_tmp) then
         518             :             !  droplet loss in decaying cloud
         519             :             !++ sungsup
         520           0 :             nsource(i,k) = nsource(i,k) + qcld(k)*(cldn_tmp - cldo_tmp)/cldo_tmp*cldliqf(i,k)*dtinv
         521           0 :             qcld(k)      = qcld(k)*(1._r8 + (cldn_tmp - cldo_tmp)/cldo_tmp)
         522             :             !-- sungsup
         523             :
         524             :             ! convert activated aerosol to interstitial in decaying cloud
         525             :
         526           0 :             dumc = (cldn_tmp - cldo_tmp)/cldo_tmp * cldliqf(i,k)
         527           0 :             do mm = 1,nele_tot
         528           0 :                dact   = raercol_cw(k,mm,nsav)*dumc
         529           0 :                raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact   ! cloud-borne aerosol
         530           0 :                raercol(k,mm,nsav)    = raercol(k,mm,nsav) - dact
         531             :             end do
         532             :
         533             :          end if
         534             :
         535             :          ! growing liquid cloud ......................................................
         536             :          !    treat the increase of cloud fraction from when cldn(i,k) > cldo(i,k)
         537             :          !    and also regenerate part of the cloud
         538           0 :          cldo_tmp = cldn_tmp
         539           0 :          cldn_tmp = lcldn(i,k)
         540             :
         541           0 :          if (cldn_tmp-cldo_tmp > 0.01_r8) then
         542             :
         543             :             ! rce-comment - use wtke at layer centers for new-cloud activation
         544           0 :             wbar  = wtke_cen(i,k)
         545           0 :             wmix  = 0._r8
         546           0 :             wmin  = 0._r8
         547           0 :             wmax  = 10._r8
         548           0 :             wdiab = 0._r8
         549             :
         550             :             ! load aerosol properties, assuming external mixtures
         551             :
         552           0 :             phase = 1 ! interstitial
         553           0 :             do m = 1, nbin
         554             :                call aero_state%loadaer( aero_props, &
         555             :                   i, i, k, &
         556             :                   m, cs, phase, na, va, &
         557           0 :                   hy, errnum, errstr)
         558           0 :                if (errnum/=0) then
         559           0 :                   call endrun('dropmixnuc : '//trim(errstr))
         560             :                end if
         561           0 :                naermod(m)  = na(i)
         562           0 :                vaerosol(m) = va(i)
         563           0 :                hygro(m)    = hy(i)
         564             :             end do
         565             :
         566             :             call activate_aerosol( &
         567             :                wbar, wmix, wdiab, wmin, wmax,                       &
         568           0 :                temp(i,k), cs(i,k), naermod, nbin,             &
         569             :                vaerosol, hygro, aero_props, fn, fm, fluxn,     &
         570           0 :                fluxm,flux_fullact(k))
         571             :
         572           0 :             factnum(i,k,:) = fn
         573             :
         574           0 :             dumc = (cldn_tmp - cldo_tmp)
         575             :
         576           0 :             do m = 1, nbin
         577           0 :                mm = aero_props%indexer(m,0)
         578           0 :                dact   = dumc*fn(m)*raer(mm)%fld(i,k) ! interstitial only
         579           0 :                qcld(k) = qcld(k) + dact
         580           0 :                nsource(i,k) = nsource(i,k) + dact*dtinv
         581           0 :                raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact  ! cloud-borne aerosol
         582           0 :                raercol(k,mm,nsav)    = raercol(k,mm,nsav) - dact
         583           0 :                dum = dumc*fm(m)
         584           0 :                do l = 1,aero_props%nmasses(m)
         585           0 :                   mm = aero_props%indexer(m,l)
         586           0 :                   dact    = dum*raer(mm)%fld(i,k) ! interstitial only
         587           0 :                   raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact  ! cloud-borne aerosol
         588           0 :                   raercol(k,mm,nsav)    = raercol(k,mm,nsav) - dact
         589             :                enddo
         590             :             enddo
         591             :
         592             :          endif
         593             :
         594             :       enddo  ! grow_shrink_main_k_loop
         595             :       ! end of k-loop for growing/shrinking cloud calcs ......................
         596             :
         597             :       ! ......................................................................
         598             :       ! start of k-loop for calc of old cloud activation tendencies ..........
         599             :       !
         600             :       ! rce-comment
         601             :       !    changed this part of code to use current cloud fraction (cldn) exclusively
         602             :       !    consider case of cldo(:)=0, cldn(k)=1, cldn(k+1)=0
         603             :       !    previous code (which used cldo below here) would have no cloud-base activation
         604             :       !       into layer k.  however, activated particles in k mix out to k+1,
         605             :       !       so they are incorrectly depleted with no replacement
         606             :
         607             :       ! old_cloud_main_k_loop
         608           0 :       do k = top_lev, pver
         609           0 :          kp1 = min0(k+1, pver)
         610           0 :          taumix_internal_pver_inv = 0.0_r8
         611             :
         612           0 :          if (lcldn(i,k) > 0.01_r8) then
         613             :
         614           0 :             wdiab = 0._r8
         615           0 :             wmix  = 0._r8                       ! single updraft
         616           0 :             wbar  = wtke(i,k)                   ! single updraft
         617           0 :             if (k == pver) wbar = wtke_cen(i,k) ! single updraft
         618           0 :             wmax  = 10._r8
         619           0 :             wmin  = 0._r8
         620             :
         621           0 :             if (lcldn(i,k) - lcldn(i,kp1) > 0.01_r8 .or. k == pver) then
         622             :
         623             :                ! cloud base
         624             :
         625             :                ! ekd(k) = wtke(i,k)*dz(i,k)/sq2pi
         626             :                ! rce-comments
         627             :                !   first, should probably have 1/zs(k) here rather than dz(i,k) because
         628             :                !      the turbulent flux is proportional to ekd(k)*zs(k),
         629             :                !      while the dz(i,k) is used to get flux divergences
         630             :                !      and mixing ratio tendency/change
         631             :                !   second and more importantly, using a single updraft velocity here
         632             :                !      means having monodisperse turbulent updraft and downdrafts.
         633             :                !      The sq2pi factor assumes a normal draft spectrum.
         634             :                !      The fluxn/fluxm from activate must be consistent with the
         635             :                !      fluxes calculated in explmix.
         636           0 :                ekd(k) = wbar/zs(k)
         637             :
         638           0 :                alogarg = max(1.e-20_r8, 1._r8/lcldn(i,k) - 1._r8)
         639           0 :                wmin    = wbar + wmix*0.25_r8*sq2pi*log(alogarg)
         640           0 :                phase   = 1   ! interstitial
         641             :
         642           0 :                do m = 1, nbin
         643             :                   ! rce-comment - use kp1 here as old-cloud activation involves
         644             :                   !   aerosol from layer below
         645             :                   call aero_state%loadaer( aero_props, &
         646             :                      i, i, kp1,  &
         647             :                      m, cs, phase, na, va,   &
         648           0 :                      hy, errnum, errstr)
         649           0 :                   if (errnum/=0) then
         650           0 :                      call endrun('dropmixnuc : '//trim(errstr))
         651             :                   end if
         652           0 :                   naermod(m)  = na(i)
         653           0 :                   vaerosol(m) = va(i)
         654           0 :                   hygro(m)    = hy(i)
         655             :                end do
         656             :
         657             :                call activate_aerosol( &
         658             :                   wbar, wmix, wdiab, wmin, wmax,                       &
         659           0 :                   temp(i,k), cs(i,k), naermod, nbin,             &
         660             :                   vaerosol, hygro, aero_props, fn, fm, fluxn,     &
         661           0 :                   fluxm, flux_fullact(k))
         662             :
         663           0 :                factnum(i,k,:) = fn
         664             :
         665           0 :                if (k < pver) then
         666           0 :                   dumc = lcldn(i,k) - lcldn(i,kp1)
         667             :                else
         668           0 :                   dumc = lcldn(i,k)
         669             :                endif
         670             :
         671           0 :                fluxntot = 0
         672             :
         673             :                ! rce-comment 1
         674             :                !    flux of activated mass into layer k (in kg/m2/s)
         675             :                !       = "actmassflux" = dumc*fluxm*raercol(kp1,lmass)*csbot(k)
         676             :                !    source of activated mass (in kg/kg/s) = flux divergence
         677             :                !       = actmassflux/(cs(i,k)*dz(i,k))
         678             :                !    so need factor of csbot_cscen = csbot(k)/cs(i,k)
         679             :                !                   dum=1./(dz(i,k))
         680           0 :                dum=csbot_cscen(k)/(dz(i,k))
         681             :
         682             :                ! rce-comment 2
         683             :                !    code for k=pver was changed to use the following conceptual model
         684             :                !    in k=pver, there can be no cloud-base activation unless one considers
         685             :                !       a scenario such as the layer being partially cloudy,
         686             :                !       with clear air at bottom and cloudy air at top
         687             :                !    assume this scenario, and that the clear/cloudy portions mix with
         688             :                !       a timescale taumix_internal = dz(i,pver)/wtke_cen(i,pver)
         689             :                !    in the absence of other sources/sinks, qact (the activated particle
         690             :                !       mixratio) attains a steady state value given by
         691             :                !          qact_ss = fcloud*fact*qtot
         692             :                !       where fcloud is cloud fraction, fact is activation fraction,
         693             :                !       qtot=qact+qint, qint is interstitial particle mixratio
         694             :                !    the activation rate (from mixing within the layer) can now be
         695             :                !       written as
         696             :                !          d(qact)/dt = (qact_ss - qact)/taumix_internal
         697             :                !                     = qtot*(fcloud*fact*wtke/dz) - qact*(wtke/dz)
         698             :                !    note that (fcloud*fact*wtke/dz) is equal to the nact/mact
         699             :                !    also, d(qact)/dt can be negative.  in the code below
         700             :                !       it is forced to be >= 0
         701             :                !
         702             :                ! steve --
         703             :                !    you will likely want to change this.  i did not really understand
         704             :                !       what was previously being done in k=pver
         705             :                !    in the cam3_5_3 code, wtke(i,pver) appears to be equal to the
         706             :                !       droplet deposition velocity which is quite small
         707             :                !    in the cam3_5_37 version, wtke is done differently and is much
         708             :                !       larger in k=pver, so the activation is stronger there
         709             :                !
         710           0 :                if (k == pver) then
         711           0 :                   taumix_internal_pver_inv = flux_fullact(k)/dz(i,k)
         712             :                end if
         713             :
         714           0 :                do m = 1, nbin
         715           0 :                   mm = aero_props%indexer(m,0)
         716           0 :                   fluxn(m) = fluxn(m)*dumc
         717           0 :                   fluxm(m) = fluxm(m)*dumc
         718           0 :                   nact(k,m) = nact(k,m) + fluxn(m)*dum
         719           0 :                   mact(k,m) = mact(k,m) + fluxm(m)*dum
         720           0 :                   if (k < pver) then
         721             :                      ! note that kp1 is used here
         722             :                      fluxntot = fluxntot &
         723           0 :                         + fluxn(m)*raercol(kp1,mm,nsav)*cs(i,k)
         724             :                   else
         725           0 :                      tmpa = raercol(kp1,mm,nsav)*fluxn(m) &
         726             :                           + raercol_cw(kp1,mm,nsav)*(fluxn(m) &
         727           0 :                           - taumix_internal_pver_inv*dz(i,k))
         728           0 :                      fluxntot = fluxntot + max(0.0_r8, tmpa)*cs(i,k)
         729             :                   end if
         730             :                end do
         731           0 :                srcn(k)      = srcn(k) + fluxntot/(cs(i,k)*dz(i,k))
         732           0 :                nsource(i,k) = nsource(i,k) + fluxntot/(cs(i,k)*dz(i,k))
         733             :
         734             :             endif  ! (cldn(i,k) - cldn(i,kp1) > 0.01 .or. k == pver)
         735             :
         736             :          else
         737             :
         738             :             ! no liquid cloud
         739           0 :             nsource(i,k) = nsource(i,k) - qcld(k)*dtinv
         740           0 :             qcld(k)      = 0
         741             :
         742           0 :             if (cldn(i,k) < 0.01_r8) then
         743             :                ! no ice cloud either
         744             :
         745             :                ! convert activated aerosol to interstitial in decaying cloud
         746             :
         747           0 :                do mm = 1,nele_tot
         748           0 :                   raercol(k,mm,nsav)    = raercol(k,mm,nsav) + raercol_cw(k,mm,nsav)  ! cloud-borne aerosol
         749           0 :                   raercol_cw(k,mm,nsav) = 0._r8
         750             :                end do
         751             :
         752             :             end if
         753             :          end if
         754             :
         755             :       end do  ! old_cloud_main_k_loop
         756             :
         757             :       ! switch nsav, nnew so that nnew is the updated aerosol
         758           0 :       ntemp = nsav
         759           0 :       nsav  = nnew
         760           0 :       nnew  = ntemp
         761             :
         762             :       ! load new droplets in layers above, below clouds
         763             :
         764           0 :       dtmin     = dtmicro
         765           0 :       ekk(top_lev-1)    = 0.0_r8
         766           0 :       ekk(pver) = 0.0_r8
         767           0 :       do k = top_lev, pver-1
         768             :          ! rce-comment -- ekd(k) is eddy-diffusivity at k/k+1 interface
         769             :          !   want ekk(k) = ekd(k) * (density at k/k+1 interface)
         770             :          !   so use pint(i,k+1) as pint is 1:pverp
         771             :          !           ekk(k)=ekd(k)*2.*pint(i,k)/(rair*(temp(i,k)+temp(i,k+1)))
         772             :          !           ekk(k)=ekd(k)*2.*pint(i,k+1)/(rair*(temp(i,k)+temp(i,k+1)))
         773           0 :          ekk(k) = ekd(k)*csbot(k)
         774             :       end do
         775             :
         776           0 :       do k = top_lev, pver
         777           0 :          km1     = max0(k-1, top_lev)
         778           0 :          ekkp(k) = zn(k)*ekk(k)*zs(k)
         779           0 :          ekkm(k) = zn(k)*ekk(k-1)*zs(km1)
         780           0 :          tinv    = ekkp(k) + ekkm(k)
         781             :
         782             :          ! rce-comment -- tinv is the sum of all first-order-loss-rates
         783             :          !    for the layer.  for most layers, the activation loss rate
         784             :          !    (for interstitial particles) is accounted for by the loss by
         785             :          !    turb-transfer to the layer above.
         786             :          !    k=pver is special, and the loss rate for activation within
         787             :          !    the layer must be added to tinv.  if not, the time step
         788             :          !    can be too big, and explmix can produce negative values.
         789             :          !    the negative values are reset to zero, resulting in an
         790             :          !    artificial source.
         791           0 :          if (k == pver) tinv = tinv + taumix_internal_pver_inv
         792             :
         793           0 :          if (tinv .gt. 1.e-6_r8) then
         794           0 :             dtt   = 1._r8/tinv
         795           0 :             dtmin = min(dtmin, dtt)
         796             :          end if
         797             :       end do
         798             :
         799           0 :       dtmix   = 0.9_r8*dtmin
         800           0 :       nsubmix = int(dtmicro/dtmix) + 1
         801           0 :       if (nsubmix > 100) then
         802             :          nsubmix_bnd = 100
         803             :       else
         804             :          nsubmix_bnd = nsubmix
         805             :       end if
         806           0 :       count_submix(nsubmix_bnd) = count_submix(nsubmix_bnd) + 1
         807           0 :       dtmix = dtmicro/nsubmix
         808             :
         809           0 :       do k = top_lev, pver
         810           0 :          kp1 = min(k+1, pver)
         811           0 :          km1 = max(k-1, top_lev)
         812             :          ! maximum overlap assumption
         813           0 :          if (cldn(i,kp1) > 1.e-10_r8) then
         814           0 :             overlapp(k) = min(cldn(i,k)/cldn(i,kp1), 1._r8)
         815             :          else
         816           0 :             overlapp(k) = 1._r8
         817             :          end if
         818           0 :          if (cldn(i,km1) > 1.e-10_r8) then
         819           0 :             overlapm(k) = min(cldn(i,k)/cldn(i,km1), 1._r8)
         820             :          else
         821           0 :             overlapm(k) = 1._r8
         822             :          end if
         823             :       end do
         824             :
         825             :
         826             :       ! rce-comment
         827             :       !    the activation source(k) = mact(k,m)*raercol(kp1,lmass)
         828             :       !       should not exceed the rate of transfer of unactivated particles
         829             :       !       from kp1 to k which = ekkp(k)*raercol(kp1,lmass)
         830             :       !    however it might if things are not "just right" in subr activate
         831             :       !    the following is a safety measure to avoid negatives in explmix
         832           0 :       do k = top_lev, pver-1
         833           0 :          do m = 1, nbin
         834           0 :             nact(k,m) = min( nact(k,m), ekkp(k) )
         835           0 :             mact(k,m) = min( mact(k,m), ekkp(k) )
         836             :          end do
         837             :       end do
         838             :
         839             :
         840             :       ! old_cloud_nsubmix_loop
         841           0 :       do n = 1, nsubmix
         842           0 :          qncld(:) = qcld(:)
         843             :          ! switch nsav, nnew so that nsav is the updated aerosol
         844           0 :          ntemp   = nsav
         845           0 :          nsav    = nnew
         846           0 :          nnew    = ntemp
         847           0 :          srcn(:) = 0.0_r8
         848             :
         849           0 :          do m = 1, nbin
         850           0 :             mm = aero_props%indexer(m,0)
         851             :
         852             :             ! update droplet source
         853             :             ! rce-comment- activation source in layer k involves particles from k+1
         854             :             !          srcn(:)=srcn(:)+nact(:,m)*(raercol(:,mm,nsav))
         855           0 :             srcn(top_lev:pver-1) = srcn(top_lev:pver-1) + nact(top_lev:pver-1,m)*(raercol(top_lev+1:pver,mm,nsav))
         856             :
         857             :             ! rce-comment- new formulation for k=pver
         858             :             !              srcn(  pver  )=srcn(  pver  )+nact(  pver  ,m)*(raercol(  pver,mm,nsav))
         859             :             tmpa = raercol(pver,mm,nsav)*nact(pver,m) &
         860           0 :                  + raercol_cw(pver,mm,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
         861           0 :             srcn(pver) = srcn(pver) + max(0.0_r8,tmpa)
         862             :          end do
         863             :          call explmix(  &
         864             :             qcld, srcn, ekkp, ekkm, overlapp,  &
         865             :             overlapm, qncld, zero, zero, pver, &
         866           0 :             dtmix, .false.)
         867             :
         868             :          ! rce-comment
         869             :          !    the interstitial particle mixratio is different in clear/cloudy portions
         870             :          !    of a layer, and generally higher in the clear portion.  (we have/had
         871             :          !    a method for diagnosing the the clear/cloudy mixratios.)  the activation
         872             :          !    source terms involve clear air (from below) moving into cloudy air (above).
         873             :          !    in theory, the clear-portion mixratio should be used when calculating
         874             :          !    source terms
         875           0 :          do m = 1, nbin
         876           0 :             mm = aero_props%indexer(m,0)
         877             :             ! rce-comment -   activation source in layer k involves particles from k+1
         878             :             !                 source(:)= nact(:,m)*(raercol(:,mm,nsav))
         879           0 :             source(top_lev:pver-1) = nact(top_lev:pver-1,m)*(raercol(top_lev+1:pver,mm,nsav))
         880             :             ! rce-comment - new formulation for k=pver
         881             :             !               source(  pver  )= nact(  pver,  m)*(raercol(  pver,mm,nsav))
         882             :             tmpa = raercol(pver,mm,nsav)*nact(pver,m) &
         883           0 :                  + raercol_cw(pver,mm,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
         884           0 :             source(pver) = max(0.0_r8, tmpa)
         885           0 :             flxconv = 0._r8
         886             :
         887             :             call explmix( &
         888             :                  raercol_cw(:,mm,nnew), source, ekkp, ekkm, overlapp, &
         889             :                  overlapm, raercol_cw(:,mm,nsav), zero, zero, pver,   &
         890           0 :                  dtmix, .false.)
         891             :
         892             :             call explmix( &
         893             :                  raercol(:,mm,nnew), source, ekkp, ekkm, overlapp,  &
         894             :                  overlapm, raercol(:,mm,nsav), zero, flxconv, pver, &
         895           0 :                  dtmix, .true., raercol_cw(:,mm,nsav))
         896             :
         897           0 :             do l = 1,aero_props%nmasses(m)
         898           0 :                mm = aero_props%indexer(m,l)
         899             :                ! rce-comment -   activation source in layer k involves particles from k+1
         900             :                !                  source(:)= mact(:,m)*(raercol(:,mm,nsav))
         901           0 :                source(top_lev:pver-1) = mact(top_lev:pver-1,m)*(raercol(top_lev+1:pver,mm,nsav))
         902             :                ! rce-comment- new formulation for k=pver
         903             :                !                 source(  pver  )= mact(  pver  ,m)*(raercol(  pver,mm,nsav))
         904             :                tmpa = raercol(pver,mm,nsav)*mact(pver,m) &
         905           0 :                     + raercol_cw(pver,mm,nsav)*(mact(pver,m) - taumix_internal_pver_inv)
         906           0 :                source(pver) = max(0.0_r8, tmpa)
         907           0 :                flxconv = 0._r8
         908             :
         909             :                call explmix( &
         910             :                     raercol_cw(:,mm,nnew), source, ekkp, ekkm, overlapp, &
         911             :                     overlapm, raercol_cw(:,mm,nsav), zero, zero, pver,   &
         912           0 :                     dtmix, .false.)
         913             :
         914             :                call explmix( &
         915             :                     raercol(:,mm,nnew), source, ekkp, ekkm, overlapp,  &
         916             :                     overlapm, raercol(:,mm,nsav), zero, flxconv, pver, &
         917           0 :                     dtmix, .true., raercol_cw(:,mm,nsav))
         918             :
         919             :             end do
         920             :          end do
         921             :
         922           0 :          if (called_from_spcam) then
         923             :          !
         924             :          ! turbulent mixing for gas species .
         925             :          !
         926           0 :                do m=1, pcnst
         927           0 :                   if (cnst_species_class(m) == cnst_spec_class_gas) then
         928           0 :                     flxconv = 0.0_r8
         929           0 :                     zerogas(:) = 0.0_r8
         930           0 :                     call explmix(rgascol(1,m,nnew),zerogas,ekkp,ekkm,overlapp,overlapm,  &
         931           0 :                                  rgascol(1,m,nsav),zero, flxconv, pver,dtmix,&
         932           0 :                                    .true., zerogas)
         933             :                   end if
         934             :                end do
         935             :          endif
         936             :
         937             :       end do ! old_cloud_nsubmix_loop
         938             :
         939             :       ! evaporate particles again if no cloud (either ice or liquid)
         940             :
         941           0 :       do k = top_lev, pver
         942           0 :          if (cldn(i,k) == 0._r8) then
         943             :             ! no ice or liquid cloud
         944           0 :             qcld(k)=0._r8
         945             :
         946             :             ! convert activated aerosol to interstitial in decaying cloud
         947           0 :             do mm = 1,nele_tot
         948           0 :                raercol(k,mm,nnew)    = raercol(k,mm,nnew) + raercol_cw(k,mm,nnew)
         949           0 :                raercol_cw(k,mm,nnew) = 0._r8
         950             :             end do
         951             :
         952             :          end if
         953             :       end do
         954             :
         955             :       ! droplet number
         956             :
         957           0 :       ndropcol(i) = 0._r8
         958           0 :       do k = top_lev, pver
         959           0 :          ndropmix(i,k) = (qcld(k) - ncldwtr(i,k))*dtinv - nsource(i,k)
         960           0 :          tendnd(i,k)   = (max(qcld(k), 1.e-6_r8) - ncldwtr(i,k))*dtinv
         961           0 :          ndropcol(i)   = ndropcol(i) + ncldwtr(i,k)*pdel(i,k)
         962             :       end do
         963           0 :       ndropcol(i) = ndropcol(i)/gravit
         964             :
         965           0 :       raertend = 0._r8
         966           0 :       qqcwtend = 0._r8
         967             :
         968           0 :       do m = 1, nbin
         969           0 :          do l = 0, aero_props%nmasses(m)
         970             :
         971           0 :             mm = aero_props%indexer(m,l)
         972           0 :             lptr = aer_cnst_idx(m,l)
         973             :
         974           0 :             raertend(top_lev:pver) = (raercol(top_lev:pver,mm,nnew) - raer(mm)%fld(i,top_lev:pver))*dtinv
         975           0 :             qqcwtend(top_lev:pver) = (raercol_cw(top_lev:pver,mm,nnew) - qqcw(mm)%fld(i,top_lev:pver))*dtinv
         976             :
         977           0 :             coltend(i,mm)    = sum( pdel(i,:)*raertend )/gravit
         978           0 :             coltend_cw(i,mm) = sum( pdel(i,:)*qqcwtend )/gravit
         979             :
         980             :             ! check for advected aerosol constituents
         981           0 :             if (lptr>0) then ! advected aerosol parts
         982           0 :                ptend%q(i,:,lptr) = 0.0_r8
         983           0 :                ptend%q(i,top_lev:pver,lptr) = raertend(top_lev:pver)         ! set tendencies for interstitial aerosol
         984             :             else
         985           0 :                raer(mm)%fld(i,:) = 0.0_r8
         986           0 :                raer(mm)%fld(i,top_lev:pver)  = raercol(top_lev:pver,mm,nnew) ! update non-advected interstitial aerosol (pbuf)
         987             :             end if
         988             :
         989           0 :             qqcw(mm)%fld(i,:) = 0.0_r8
         990           0 :             qqcw(mm)%fld(i,top_lev:pver) = raercol_cw(top_lev:pver,mm,nnew)  ! update cloud-borne aerosol
         991             :
         992             :          end do
         993             :       end do
         994             :
         995           0 :       if (called_from_spcam) then
         996             :       !
         997             :       ! Gas tendency
         998             :       !
         999           0 :            do m=1, pcnst
        1000           0 :               if (cnst_species_class(m) == cnst_spec_class_gas) then
        1001           0 :                 ptend%lq(m) = .true.
        1002           0 :                 ptend%q(i, :, m) = (rgascol(:,m,nnew)-rgas(i,:,m)) * dtinv
        1003             :               end if
        1004             :            end do
        1005             :       endif
        1006             :
        1007             :    end do  ! overall_main_i_loop
        1008             :    ! end of main loop over i/longitude ....................................
        1009             :
        1010           0 :    call outfld('NDROPCOL', ndropcol, pcols, lchnk)
        1011           0 :    call outfld('NDROPSRC', nsource,  pcols, lchnk)
        1012           0 :    call outfld('NDROPMIX', ndropmix, pcols, lchnk)
        1013           0 :    call outfld('WTKE    ', wtke,     pcols, lchnk)
        1014             :
        1015           0 :    if(called_from_spcam) then
        1016           0 :         call outfld('SPLCLOUD  ', cldn    , pcols, lchnk   )
        1017           0 :         call outfld('SPKVH     ', kvh     , pcols, lchnk   )
        1018             :    endif
        1019             :
        1020           0 :    call ccncalc(aero_state, aero_props, state, cs, ccn)
        1021           0 :    do l = 1, psat
        1022           0 :       call outfld(ccn_name(l), ccn(1,1,l), pcols, lchnk)
        1023             :    enddo
        1024             :
        1025             :    ! do column tendencies
        1026           0 :    do m = 1, nbin
        1027           0 :       do l = 0,aero_props%nmasses(m)
        1028           0 :          mm = aero_props%indexer(m,l)
        1029           0 :          call outfld(fieldname(mm),    coltend(:,mm),    pcols, lchnk)
        1030           0 :          call outfld(fieldname_cw(mm), coltend_cw(:,mm), pcols, lchnk)
        1031             :       end do
        1032             :    end do
        1033             :
        1034           0 :    if(called_from_spcam) then
        1035             :    !
        1036             :    ! output column-integrated Gas tendency (this should be zero)
        1037             :    !
        1038           0 :         do m=1, pcnst
        1039           0 :             if (cnst_species_class(m) == cnst_spec_class_gas) then
        1040           0 :               do i=1, ncol
        1041           0 :                  coltendgas(i) = sum( pdel(i,:)*ptend%q(i,:,m) )/gravit
        1042             :               end do
        1043           0 :               fieldnamegas = trim(cnst_name(m)) // '_mixnuc1sp'
        1044           0 :               call outfld( trim(fieldnamegas), coltendgas, pcols, lchnk)
        1045             :             end if
        1046             :         end do
        1047           0 :         deallocate(rgascol, coltendgas)
        1048             :    end if
        1049             :
        1050           0 :    deallocate( &
        1051             :       nact,       &
        1052           0 :       mact,       &
        1053           0 :       raer,       &
        1054           0 :       qqcw,       &
        1055           0 :       raercol,    &
        1056           0 :       raercol_cw, &
        1057           0 :       coltend,    &
        1058           0 :       coltend_cw, &
        1059           0 :       naermod,    &
        1060           0 :       hygro,      &
        1061           0 :       vaerosol,   &
        1062           0 :       fn,         &
        1063           0 :       fm,         &
        1064           0 :       fluxn,      &
        1065           0 :       fluxm       )
        1066             :
        1067           0 : end subroutine dropmixnuc
        1068             :
        1069             : !===============================================================================
        1070             :
        1071           0 : subroutine explmix( q, src, ekkp, ekkm, overlapp, overlapm, &
        1072           0 :    qold, surfrate, flxconv, pver, dt, is_unact, qactold )
        1073             :
        1074             :    !  explicit integration of droplet/aerosol mixing
        1075             :    !     with source due to activation/nucleation
        1076             :
        1077             :
        1078             :    integer, intent(in) :: pver ! number of levels
        1079             :    real(r8), intent(out) :: q(pver) ! mixing ratio to be updated
        1080             :    real(r8), intent(in) :: qold(pver) ! mixing ratio from previous time step
        1081             :    real(r8), intent(in) :: src(pver) ! source due to activation/nucleation (/s)
        1082             :    real(r8), intent(in) :: ekkp(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
        1083             :    ! below layer k  (k,k+1 interface)
        1084             :    real(r8), intent(in) :: ekkm(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
        1085             :    ! above layer k  (k,k+1 interface)
        1086             :    real(r8), intent(in) :: overlapp(pver) ! cloud overlap below
        1087             :    real(r8), intent(in) :: overlapm(pver) ! cloud overlap above
        1088             :    real(r8), intent(in) :: surfrate ! surface exchange rate (/s)
        1089             :    real(r8), intent(in) :: flxconv ! convergence of flux from surface
        1090             :    real(r8), intent(in) :: dt ! time step (s)
        1091             :    logical, intent(in) :: is_unact ! true if this is an unactivated species
        1092             :    real(r8), intent(in),optional :: qactold(pver)
        1093             :    ! mixing ratio of ACTIVATED species from previous step
        1094             :    ! *** this should only be present
        1095             :    !     if the current species is unactivated number/sfc/mass
        1096             :
        1097             :    integer k,kp1,km1
        1098             :
        1099           0 :    if ( is_unact ) then
        1100             :       !     the qactold*(1-overlap) terms are resuspension of activated material
        1101           0 :       do k=top_lev,pver
        1102           0 :          kp1=min(k+1,pver)
        1103           0 :          km1=max(k-1,top_lev)
        1104           0 :          q(k) = qold(k) + dt*( - src(k) + ekkp(k)*(qold(kp1) - qold(k) +       &
        1105           0 :             qactold(kp1)*(1.0_r8-overlapp(k)))               &
        1106           0 :             + ekkm(k)*(qold(km1) - qold(k) +     &
        1107           0 :             qactold(km1)*(1.0_r8-overlapm(k))) )
        1108             :          !        force to non-negative
        1109             :          !        if(q(k)<-1.e-30)then
        1110             :          !           write(iulog,*)'q=',q(k),' in explmix'
        1111           0 :          q(k)=max(q(k),0._r8)
        1112             :          !        endif
        1113             :       end do
        1114             :
        1115             :       !     diffusion loss at base of lowest layer
        1116           0 :       q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
        1117             :       !        force to non-negative
        1118             :       !        if(q(pver)<-1.e-30)then
        1119             :       !           write(iulog,*)'q=',q(pver),' in explmix'
        1120           0 :       q(pver)=max(q(pver),0._r8)
        1121             :       !        endif
        1122             :    else
        1123           0 :       do k=top_lev,pver
        1124           0 :          kp1=min(k+1,pver)
        1125           0 :          km1=max(k-1,top_lev)
        1126           0 :          q(k) = qold(k) + dt*(src(k) + ekkp(k)*(overlapp(k)*qold(kp1)-qold(k)) +      &
        1127           0 :             ekkm(k)*(overlapm(k)*qold(km1)-qold(k)) )
        1128             :          !        force to non-negative
        1129             :          !        if(q(k)<-1.e-30)then
        1130             :          !           write(iulog,*)'q=',q(k),' in explmix'
        1131           0 :          q(k)=max(q(k),0._r8)
        1132             :          !        endif
        1133             :       end do
        1134             :       !     diffusion loss at base of lowest layer
        1135           0 :       q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
        1136             :       !        force to non-negative
        1137             :       !        if(q(pver)<-1.e-30)then
        1138             :       !           write(iulog,*)'q=',q(pver),' in explmix'
        1139           0 :       q(pver)=max(q(pver),0._r8)
        1140             :
        1141             :    end if
        1142             :
        1143           0 : end subroutine explmix
        1144             :
        1145             : !===============================================================================
        1146             :
        1147           0 : subroutine activate_aerosol(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,  &
        1148           0 :    na, nbins, volume, hygro, aero_props, &
        1149           0 :    fn, fm, fluxn, fluxm, flux_fullact, smax_prescribed, in_cloud_in, smax_f)
        1150             :
        1151             :    !      calculates number, surface, and mass fraction of aerosols activated as CCN
        1152             :    !      calculates flux of cloud droplets, surface area, and aerosol mass into cloud
        1153             :    !      assumes an internal mixture within each of up to nbin multiple aerosol bins
        1154             :    !      a gaussiam spectrum of updrafts can be treated.
        1155             :
        1156             :    !      mks units
        1157             :
        1158             :    !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
        1159             :    !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
        1160             :
        1161             :    !      input
        1162             :
        1163             :    real(r8), intent(in) :: wbar          ! grid cell mean vertical velocity (m/s)
        1164             :    real(r8), intent(in) :: sigw          ! subgrid standard deviation of vertical vel (m/s)
        1165             :    real(r8), intent(in) :: wdiab         ! diabatic vertical velocity (0 if adiabatic)
        1166             :    real(r8), intent(in) :: wminf         ! minimum updraft velocity for integration (m/s)
        1167             :    real(r8), intent(in) :: wmaxf         ! maximum updraft velocity for integration (m/s)
        1168             :    real(r8), intent(in) :: tair          ! air temperature (K)
        1169             :    real(r8), intent(in) :: rhoair        ! air density (kg/m3)
        1170             :    real(r8), intent(in) :: na(:)      ! aerosol number concentration (/m3)
        1171             :    integer,  intent(in) :: nbins      ! number of aerosol bins
        1172             :    real(r8), intent(in) :: volume(:)  ! aerosol volume concentration (m3/m3)
        1173             :    real(r8), intent(in) :: hygro(:)   ! hygroscopicity of aerosol mode
        1174             :
        1175             :    class(aerosol_properties), intent(in) :: aero_props
        1176             :
        1177             :    !      output
        1178             :
        1179             :    real(r8), intent(out) :: fn(:)      ! number fraction of aerosols activated
        1180             :    real(r8), intent(out) :: fm(:)      ! mass fraction of aerosols activated
        1181             :    real(r8), intent(out) :: fluxn(:)   ! flux of activated aerosol number fraction into cloud (cm/s)
        1182             :    real(r8), intent(out) :: fluxm(:)   ! flux of activated aerosol mass fraction into cloud (cm/s)
        1183             :    real(r8), intent(out) :: flux_fullact   ! flux of activated aerosol fraction assuming 100% activation (cm/s)
        1184             :    !    rce-comment
        1185             :    !    used for consistency check -- this should match (ekd(k)*zs(k))
        1186             :    !    also, fluxm/flux_fullact gives fraction of aerosol mass flux
        1187             :    !       that is activated
        1188             :
        1189             :    !      optional
        1190             :    real(r8), optional, intent(in) :: smax_prescribed  ! prescribed max. supersaturation for secondary activation
        1191             :    logical,  optional, intent(in) :: in_cloud_in      ! switch to modify calculations when above cloud base
        1192             :    real(r8), optional, intent(in) :: smax_f           ! droplet and rain size distr factor in the smax calculation
        1193             :                                                       ! used when in_cloud=.true.
        1194             :
        1195             :    !      local
        1196             :
        1197             :    integer, parameter:: nx=200
        1198             :    real(r8) integ,integf
        1199             :    real(r8), parameter :: p0 = 1013.25e2_r8    ! reference pressure (Pa)
        1200             :    real(r8) pres ! pressure (Pa)
        1201             :    real(r8) diff0,conduct0
        1202             :    real(r8) es ! saturation vapor pressure
        1203             :    real(r8) qs ! water vapor saturation mixing ratio
        1204             :    real(r8) dqsdt ! change in qs with temperature
        1205             :    real(r8) g ! thermodynamic function (m2/s)
        1206           0 :    real(r8) zeta(nbins), eta(nbins)
        1207             :    real(r8) alpha
        1208             :    real(r8) gamma
        1209             :    real(r8) beta
        1210             :    real(r8) sqrtg
        1211           0 :    real(r8) :: amcube(nbins) ! cube of dry bin radius (m)
        1212           0 :    real(r8) smc(nbins) ! critical supersaturation for number bin radius
        1213             :    real(r8) sumflx_fullact
        1214           0 :    real(r8) sumflxn(nbins)
        1215           0 :    real(r8) sumflxm(nbins)
        1216           0 :    real(r8) sumfn(nbins)
        1217           0 :    real(r8) sumfm(nbins)
        1218           0 :    real(r8) fnold(nbins)   ! number fraction activated
        1219           0 :    real(r8) fmold(nbins)   ! mass fraction activated
        1220             :    real(r8) wold,gold
        1221             :    real(r8) wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb
        1222             :    real(r8) dfmin,dfmax,fnew,fold,fnmin,fnbar,fmbar
        1223             :    real(r8) alw,sqrtalw
        1224             :    real(r8) smax
        1225             :    real(r8) z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
        1226           0 :    real(r8) etafactor1,etafactor2(nbins),etafactor2max
        1227             :    real(r8) grow
        1228             :    character(len=*), parameter :: subname='activate_aerosol'
        1229             :
        1230             :    logical :: in_cloud
        1231             :    integer m,n
        1232             :    !      numerical integration parameters
        1233             :    real(r8), parameter :: eps=0.3_r8,fmax=0.99_r8,sds=3._r8
        1234             :
        1235             :    real(r8), parameter :: namin=1.e6_r8   ! minimum aerosol number concentration (/m3)
        1236             :
        1237             :    integer ndist(nx)  ! accumulates frequency distribution of integration bins required
        1238             :    data ndist/nx*0/
        1239             :    save ndist
        1240             :
        1241           0 :    if (present(in_cloud_in)) then
        1242           0 :       if (.not. present(smax_f)) call endrun(subname//' error: smax_f must be supplied when in_cloud is used')
        1243           0 :       in_cloud = in_cloud_in
        1244             :    else
        1245             :       in_cloud = .false.
        1246             :    end if
        1247             :
        1248           0 :    fn(:)=0._r8
        1249           0 :    fm(:)=0._r8
        1250           0 :    fluxn(:)=0._r8
        1251           0 :    fluxm(:)=0._r8
        1252           0 :    flux_fullact=0._r8
        1253             :
        1254           0 :    if(nbins.eq.1.and.na(1).lt.1.e-20_r8)return
        1255             :
        1256           0 :    if(sigw.le.1.e-5_r8.and.wbar.le.0._r8)return
        1257             :
        1258           0 :    if ( present( smax_prescribed ) ) then
        1259           0 :       if (smax_prescribed <= 0.0_r8) return
        1260             :    end if
        1261             :
        1262           0 :    pres=rair*rhoair*tair
        1263           0 :    diff0=0.211e-4_r8*(p0/pres)*(tair/tmelt)**1.94_r8
        1264           0 :    conduct0=(5.69_r8+0.017_r8*(tair-tmelt))*4.186e2_r8*1.e-5_r8 ! convert to J/m/s/deg
        1265           0 :    call qsat(tair, pres, es, qs)
        1266           0 :    dqsdt=latvap/(rh2o*tair*tair)*qs
        1267           0 :    alpha=gravit*(latvap/(cpair*rh2o*tair*tair)-1._r8/(rair*tair))
        1268           0 :    gamma=(1.0_r8+latvap/cpair*dqsdt)/(rhoair*qs)
        1269           0 :    etafactor2max=1.e10_r8/(alpha*wmaxf)**1.5_r8 ! this should make eta big if na is very small.
        1270             :
        1271             :    grow  = 1._r8/(rhoh2o/(diff0*rhoair*qs)  &
        1272           0 :            + latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair) - 1._r8))
        1273           0 :    sqrtg = sqrt(grow)
        1274           0 :    beta  = 2._r8*pi*rhoh2o*grow*gamma
        1275             :
        1276           0 :    do m=1,nbins
        1277             :
        1278           0 :       if(volume(m).gt.1.e-39_r8.and.na(m).gt.1.e-39_r8)then
        1279             :          ! number mode radius (m)
        1280           0 :          amcube(m)=aero_props%amcube(m, volume(m),na(m))
        1281             :          ! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
        1282             :          ! should depend on mean radius of mode to account for gas kinetic effects
        1283             :          ! see Fountoukis and Nenes, JGR2005 and Meskhidze et al., JGR2006
        1284             :          ! for approriate size to use for effective diffusivity.
        1285           0 :          etafactor2(m)=1._r8/(na(m)*beta*sqrtg)
        1286           0 :          if(hygro(m).gt.1.e-10_r8)then
        1287           0 :             smc(m)=2._r8*aten*sqrt(aten/(27._r8*hygro(m)*amcube(m))) ! only if variable size dist
        1288             :          else
        1289           0 :             smc(m)=100._r8
        1290             :          endif
        1291             :       else
        1292           0 :          smc(m)=1._r8
        1293           0 :          etafactor2(m)=etafactor2max ! this should make eta big if na is very small.
        1294             :       endif
        1295             :
        1296             :    enddo
        1297             :
        1298           0 :    if(sigw.gt.1.e-5_r8)then ! spectrum of updrafts
        1299             :
        1300           0 :       wmax=min(wmaxf,wbar+sds*sigw)
        1301           0 :       wmin=max(wminf,-wdiab)
        1302           0 :       wmin=max(wmin,wbar-sds*sigw)
        1303           0 :       w=wmin
        1304           0 :       dwmax=eps*sigw
        1305           0 :       dw=dwmax
        1306           0 :       dfmax=0.2_r8
        1307           0 :       dfmin=0.1_r8
        1308           0 :       if (wmax <= w) return
        1309           0 :       do m=1,nbins
        1310           0 :          sumflxn(m)=0._r8
        1311           0 :          sumfn(m)=0._r8
        1312           0 :          fnold(m)=0._r8
        1313           0 :          sumflxm(m)=0._r8
        1314           0 :          sumfm(m)=0._r8
        1315           0 :          fmold(m)=0._r8
        1316             :       enddo
        1317           0 :       sumflx_fullact=0._r8
        1318             :
        1319           0 :       fold=0._r8
        1320           0 :       wold=0._r8
        1321           0 :       gold=0._r8
        1322             :
        1323           0 :       dwmin = min( dwmax, 0.01_r8 )
        1324           0 :       do n = 1, nx
        1325             :
        1326           0 : 100      wnuc=w+wdiab
        1327             :          !           write(iulog,*)'wnuc=',wnuc
        1328           0 :          alw=alpha*wnuc
        1329           0 :          sqrtalw=sqrt(alw)
        1330           0 :          etafactor1=alw*sqrtalw
        1331             :
        1332           0 :          do m=1,nbins
        1333           0 :             eta(m)=etafactor1*etafactor2(m)
        1334           0 :             zeta(m)=twothird*sqrtalw*aten/sqrtg
        1335             :          enddo
        1336             :
        1337           0 :          if ( present( smax_prescribed ) ) then
        1338           0 :             smax = smax_prescribed
        1339             :          else
        1340           0 :             smax = aero_props%maxsat(zeta,eta,smc)
        1341             :          endif
        1342             :
        1343           0 :          call aero_props%actfracs( nbins, smc(nbins), smax, fnew, fm(nbins) )
        1344             :
        1345           0 :          dwnew = dw
        1346           0 :          if(fnew-fold.gt.dfmax.and.n.gt.1)then
        1347             :             !              reduce updraft increment for greater accuracy in integration
        1348           0 :             if (dw .gt. 1.01_r8*dwmin) then
        1349           0 :                dw=0.7_r8*dw
        1350           0 :                dw=max(dw,dwmin)
        1351           0 :                w=wold+dw
        1352           0 :                go to 100
        1353             :             else
        1354             :                dwnew = dwmin
        1355             :             endif
        1356             :          endif
        1357             :
        1358           0 :          if(fnew-fold.lt.dfmin)then
        1359             :             !              increase updraft increment to accelerate integration
        1360           0 :             dwnew=min(1.5_r8*dw,dwmax)
        1361             :          endif
        1362           0 :          fold=fnew
        1363             :
        1364           0 :          z=(w-wbar)/(sigw*sq2)
        1365           0 :          g=exp(-z*z)
        1366           0 :          fnmin=1._r8
        1367             :
        1368           0 :          do m=1,nbins
        1369             :             !              modal
        1370           0 :             call aero_props%actfracs( m, smc(m), smax, fn(m), fm(m) )
        1371           0 :             fnmin=min(fn(m),fnmin)
        1372             :             !               integration is second order accurate
        1373             :             !               assumes linear variation of f*g with w
        1374           0 :             fnbar=(fn(m)*g+fnold(m)*gold)
        1375           0 :             fmbar=(fm(m)*g+fmold(m)*gold)
        1376           0 :             wb=(w+wold)
        1377           0 :             if(w.gt.0._r8)then
        1378             :                sumflxn(m)=sumflxn(m)+sixth*(wb*fnbar           &
        1379           0 :                   +(fn(m)*g*w+fnold(m)*gold*wold))*dw
        1380             :                sumflxm(m)=sumflxm(m)+sixth*(wb*fmbar           &
        1381           0 :                   +(fm(m)*g*w+fmold(m)*gold*wold))*dw
        1382             :             endif
        1383           0 :             sumfn(m)=sumfn(m)+0.5_r8*fnbar*dw
        1384           0 :             fnold(m)=fn(m)
        1385           0 :             sumfm(m)=sumfm(m)+0.5_r8*fmbar*dw
        1386           0 :             fmold(m)=fm(m)
        1387             :          enddo
        1388             :          !           same form as sumflxm but replace the fm with 1.0
        1389             :          sumflx_fullact = sumflx_fullact &
        1390           0 :             + sixth*(wb*(g+gold) + (g*w+gold*wold))*dw
        1391           0 :          gold=g
        1392           0 :          wold=w
        1393           0 :          dw=dwnew
        1394           0 :          if (n > 1 .and. (w > wmax .or. fnmin > fmax)) exit
        1395           0 :          w=w+dw
        1396           0 :          if (n == nx) then
        1397           0 :             write(iulog,*)'do loop is too short in activate'
        1398           0 :             write(iulog,*)'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw
        1399           0 :             write(iulog,*)'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab
        1400           0 :             write(iulog,*)'wnuc=',wnuc
        1401           0 :             write(iulog,*)'na=',(na(m),m=1,nbins)
        1402           0 :             write(iulog,*)'fn=',(fn(m),m=1,nbins)
        1403             :             !   dump all subr parameters to allow testing with standalone code
        1404             :             !   (build a driver that will read input and call activate)
        1405           0 :             write(iulog,*)'wbar,sigw,wdiab,tair,rhoair,nbins='
        1406           0 :             write(iulog,*) wbar,sigw,wdiab,tair,rhoair,nbins
        1407           0 :             write(iulog,*)'na=',na
        1408           0 :             write(iulog,*)'volume=', (volume(m),m=1,nbins)
        1409           0 :             write(iulog,*)'hygro='
        1410           0 :             write(iulog,*) hygro
        1411           0 :             call endrun(subname)
        1412             :          end if
        1413             :
        1414             :       enddo
        1415             :
        1416           0 :       ndist(n)=ndist(n)+1
        1417           0 :       if(w.lt.wmaxf)then
        1418             :
        1419             :          !            contribution from all updrafts stronger than wmax
        1420             :          !            assuming constant f (close to fmax)
        1421           0 :          wnuc=w+wdiab
        1422             :
        1423           0 :          z1=(w-wbar)/(sigw*sq2)
        1424           0 :          z2=(wmaxf-wbar)/(sigw*sq2)
        1425           0 :          g=exp(-z1*z1)
        1426           0 :          integ=sigw*0.5_r8*sq2*sqpi*(erf(z2)-erf(z1))
        1427             :          !            consider only upward flow into cloud base when estimating flux
        1428           0 :          wf1=max(w,zero)
        1429           0 :          zf1=(wf1-wbar)/(sigw*sq2)
        1430           0 :          gf1=exp(-zf1*zf1)
        1431           0 :          wf2=max(wmaxf,zero)
        1432           0 :          zf2=(wf2-wbar)/(sigw*sq2)
        1433           0 :          gf2=exp(-zf2*zf2)
        1434           0 :          gf=(gf1-gf2)
        1435           0 :          integf=wbar*sigw*0.5_r8*sq2*sqpi*(erf(zf2)-erf(zf1))+sigw*sigw*gf
        1436             :
        1437           0 :          do m=1,nbins
        1438           0 :             sumflxn(m)=sumflxn(m)+integf*fn(m)
        1439           0 :             sumfn(m)=sumfn(m)+fn(m)*integ
        1440           0 :             sumflxm(m)=sumflxm(m)+integf*fm(m)
        1441           0 :             sumfm(m)=sumfm(m)+fm(m)*integ
        1442             :          enddo
        1443             :          !           same form as sumflxm but replace the fm with 1.0
        1444           0 :          sumflx_fullact = sumflx_fullact + integf
        1445             :          !            sumg=sumg+integ
        1446             :       endif
        1447             :
        1448             :
        1449           0 :       do m=1,nbins
        1450           0 :          fn(m)=sumfn(m)/(sq2*sqpi*sigw)
        1451             :          !            fn(m)=sumfn(m)/(sumg)
        1452           0 :          if(fn(m).gt.1.01_r8)then
        1453           0 :             write(iulog,*)'fn=',fn(m),' > 1 in activate'
        1454           0 :             write(iulog,*)'w,m,na,amcube=',w,m,na(m),amcube(m)
        1455           0 :             write(iulog,*)'integ,sumfn,sigw=',integ,sumfn(m),sigw
        1456           0 :             call endrun('activate')
        1457             :          endif
        1458           0 :          fluxn(m)=sumflxn(m)/(sq2*sqpi*sigw)
        1459           0 :          fm(m)=sumfm(m)/(sq2*sqpi*sigw)
        1460             :          !            fm(m)=sumfm(m)/(sumg)
        1461           0 :          if(fm(m).gt.1.01_r8)then
        1462           0 :             write(iulog,*)'fm=',fm(m),' > 1 in activate'
        1463             :          endif
        1464           0 :          fluxm(m)=sumflxm(m)/(sq2*sqpi*sigw)
        1465             :       enddo
        1466             :       !        same form as fluxm
        1467           0 :       flux_fullact = sumflx_fullact/(sq2*sqpi*sigw)
        1468             :
        1469             :    else
        1470             :
        1471             :       !        single updraft
        1472           0 :       wnuc=wbar+wdiab
        1473             :
        1474           0 :       if(wnuc.gt.0._r8)then
        1475             :
        1476           0 :          w=wbar
        1477             :
        1478           0 :          if(in_cloud) then
        1479             :
        1480           0 :             if (smax_f > 0._r8) then
        1481           0 :                smax = alpha*w/(2.0_r8*pi*rhoh2o*grow*gamma*smax_f)
        1482             :             else
        1483           0 :                smax = 1.e-20_r8
        1484             :             end if
        1485             :
        1486             :          else ! at cloud base
        1487           0 :             alw        = alpha*wnuc
        1488           0 :             sqrtalw    = sqrt(alw)
        1489           0 :             etafactor1 = alw*sqrtalw
        1490             :
        1491           0 :             do m = 1, nbins
        1492           0 :                eta(m)  = etafactor1*etafactor2(m)
        1493           0 :                zeta(m) = twothird*sqrtalw*aten/sqrtg
        1494             :             end do
        1495           0 :             if ( present(smax_prescribed) ) then
        1496           0 :                smax = smax_prescribed
        1497             :             else
        1498           0 :                smax = aero_props%maxsat(zeta,eta,smc)
        1499             :             end if
        1500             :          end if
        1501             :
        1502           0 :          do m=1,nbins
        1503             :
        1504           0 :             call aero_props%actfracs( m, smc(m), smax, fn(m), fm(m) )
        1505             :
        1506           0 :             if(wbar.gt.0._r8)then
        1507           0 :                fluxn(m)=fn(m)*w
        1508           0 :                fluxm(m)=fm(m)*w
        1509             :             endif
        1510             :          enddo
        1511           0 :          flux_fullact = w
        1512             :       endif
        1513             :
        1514             :    endif
        1515             :
        1516             : end subroutine activate_aerosol
        1517             :
        1518             : !===============================================================================
        1519             :
        1520           0 : subroutine ccncalc(aero_state, aero_props, state, cs, ccn)
        1521             :
        1522             :    ! calculates number concentration of aerosols activated as CCN at
        1523             :    ! supersaturation supersat.
        1524             :    ! assumes an internal mixture of a multiple externally-mixed aerosol modes
        1525             :    ! cgs units
        1526             :
        1527             :    ! Ghan et al., Atmos. Res., 1993, 198-221.
        1528             :
        1529             :    ! arguments
        1530             :    class(aerosol_state), intent(in) :: aero_state
        1531             :    class(aerosol_properties), intent(in) :: aero_props
        1532             :
        1533             :    type(physics_state), target, intent(in)    :: state
        1534             :
        1535             :    real(r8), intent(in)  :: cs(pcols,pver)       ! air density (kg/m3)
        1536             :    real(r8), intent(out) :: ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat (#/m3)
        1537             :
        1538             :    ! local
        1539             :
        1540             :    integer :: ncol  ! number of columns
        1541             :    integer :: nbin  ! number of bins
        1542           0 :    real(r8), pointer :: tair(:,:)     ! air temperature (K)
        1543             :
        1544             :    real(r8) naerosol(pcols) ! interstit+activated aerosol number conc (/m3)
        1545             :    real(r8) vaerosol(pcols) ! interstit+activated aerosol volume conc (m3/m3)
        1546             :
        1547             :    real(r8) amcube(pcols)
        1548           0 :    real(r8), allocatable :: argfactor(:)
        1549             :    real(r8) surften_coef
        1550             :    real(r8) a(pcols) ! surface tension parameter
        1551             :    real(r8) hygro(pcols)  ! aerosol hygroscopicity
        1552             :    real(r8) sm(pcols)  ! critical supersaturation at mode radius
        1553             :    real(r8) arg(pcols)
        1554             :    integer l,m,i,k, astat
        1555             :    real(r8) smcoef(pcols)
        1556             :    integer phase ! phase of aerosol
        1557             :
        1558             :    integer :: errnum
        1559             :    character(len=shr_kind_cs) :: errstr
        1560             :
        1561             :    !     mathematical constants
        1562             :    real(r8), parameter :: super(psat) = supersat(:psat)*0.01_r8
        1563             :    real(r8), parameter :: smcoefcoef  = 2._r8/sqrt(27._r8)
        1564             :
        1565             :    !-------------------------------------------------------------------------------
        1566             :
        1567           0 :    nbin  = aero_props%nbins()
        1568           0 :    ncol  = state%ncol
        1569           0 :    tair  => state%t
        1570             :
        1571           0 :    allocate( argfactor(nbin), stat=astat )
        1572           0 :    if (astat/=0) then
        1573           0 :       call endrun('ndrop::ccncalc : not able to allocate argfactor')
        1574             :    end if
        1575             :
        1576           0 :    surften_coef=2._r8*mwh2o*surften/(r_universal*rhoh2o)
        1577             :
        1578           0 :    do m=1,nbin
        1579           0 :       argfactor(m)=twothird/(sq2*aero_props%alogsig(m))
        1580             :    end do
        1581             :
        1582           0 :    ccn = 0._r8
        1583           0 :    do k=top_lev,pver
        1584             :
        1585           0 :       do i=1,ncol
        1586           0 :          a(i)=surften_coef/tair(i,k)
        1587           0 :          smcoef(i)=smcoefcoef*a(i)*sqrt(a(i))
        1588             :       end do
        1589             :
        1590           0 :       do m=1,nbin
        1591             :
        1592           0 :          phase=3 ! interstitial+cloudborne
        1593             :
        1594             :          call aero_state%loadaer( aero_props, &
        1595             :             1, ncol, k, &
        1596             :             m, cs, phase, naerosol, vaerosol, &
        1597           0 :             hygro, errnum, errstr)
        1598           0 :          if (errnum/=0) then
        1599           0 :             call endrun('ccncalc : '//trim(errstr))
        1600             :          end if
        1601             :
        1602           0 :          where(naerosol(:ncol)>1.e-3_r8 .and. hygro(:ncol)>1.e-10_r8)
        1603           0 :             amcube(:ncol)=aero_props%amcube(m, vaerosol(:ncol), naerosol(:ncol) )
        1604           0 :             sm(:ncol)=smcoef(:ncol)/sqrt(hygro(:ncol)*amcube(:ncol)) ! critical supersaturation
        1605             :          elsewhere
        1606             :             sm(:ncol)=1._r8 ! value shouldn't matter much since naerosol is small
        1607             :          endwhere
        1608           0 :          do l=1,psat
        1609           0 :             do i=1,ncol
        1610           0 :                arg(i)=argfactor(m)*log(sm(i)/super(l))
        1611           0 :                ccn(i,k,l)=ccn(i,k,l)+naerosol(i)*0.5_r8*(1._r8-erf(arg(i)))
        1612             :             enddo
        1613             :          enddo
        1614             :       enddo
        1615             :    enddo
        1616           0 :    ccn(:ncol,:,:)=ccn(:ncol,:,:)*1.e-6_r8 ! convert from #/m3 to #/cm3
        1617             :
        1618           0 :    deallocate( argfactor )
        1619             :
        1620           0 : end subroutine ccncalc
        1621             :
        1622             : !===============================================================================
        1623             : end module ndrop

Generated by: LCOV version 1.14