chem/module_cam_mam_mz_aerosols_intr.F

   1 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   2 #define WRF_PORT
   3 #define MODAL_AERO
   4 module mz_aerosols_intr
   5
   6   use shr_kind_mod,only: r8 => shr_kind_r8
   7 #ifndef WRF_PORT
   8   use spmd_utils,  only: masterproc
   9   use pmgrid,      only: plev, plevp
  10   use ppgrid,      only: pcols, pver
  11   use constituents,only: pcnst, cnst_name, cnst_get_ind
  12   use phys_control,only: phys_getopts
  13   use abortutils,  only: endrun
  14 #else
  15   use module_cam_support, only: masterproc, pcnst =>pcnst_runtime, pcols, pver, endrun
  16   use constituents,only: cnst_name, cnst_get_ind
  17 #endif
  18 #ifdef MODAL_AERO
  19     use modal_aero_data, only: ntot_amode
  20 #endif
  21 #ifndef WRF_PORT
  22   use cam_logfile, only: iulog
  23 #else
  24   use module_cam_support, only: iulog
  25 #endif
  26
  27   implicit none
  28
  29   private          ! Make default type private to the module
  30
  31   save
  32
  33   !
  34   ! Public interfaces
  35   !
  36
  37   public :: mz_aero_initialize
  38   public :: mz_aero_wet_intr                             ! interface to wet deposition
  39 #ifndef WRF_PORT
  40   public :: mz_aero_setopts
  41   public :: mz_aero_defaultopts
  42 #endif
  43   public :: aer_wetdep_list
  44   public :: do_cam_sulfchem
  45 #ifdef MODAL_AERO
  46 #ifndef WRF_PORT
  47   public :: mz_aero_dry_intr                             ! interface to dry deposition
  48   public :: aer_drydep_list
  49 #endif
  50   public :: modal_aero_bcscavcoef_init  ! initialize impaction scavenging table
  51   public :: modal_aero_bcscavcoef_get   ! get impaction scavenging coefficients
  52 #endif
  53 #ifndef WRF_PORT
  54   public :: mz_prescribed_dust              ! return true if MOZART code provides prescribed dust
  55 #endif
  56
  57   integer :: num_mz_aerosols
  58   integer, allocatable :: mz_aerosol_ids(:) ! (/ id_cb2, id_oc2, id_so4, id_sa1, id_sa2, id_sa3, id_sa4 /)
  59 #ifndef WRF_PORT
  60   character(len=16), dimension(pcnst) :: aer_wetdep_list = ' '
  61 #else
  62   character(len=16), allocatable, dimension(:) :: aer_wetdep_list
  63 #endif
  64   logical :: use_cam_sulfchem = .false.
  65   logical :: do_cam_sulfchem
  66   logical :: inv_o3, inv_oh, inv_no3, inv_ho2
  67   integer, pointer :: id_so2, id_so4, id_dms, id_o3, id_h2o2, id_oh, id_no3, id_ho2
  68   integer, target  :: spc_ids(8)
  69
  70 #ifdef MODAL_AERO
  71 #ifndef WRF_PORT
  72   character(len=16), dimension(pcnst) :: aer_drydep_list = ''
  73 #else
  74   character(len=16), allocatable,dimension(:) :: aer_drydep_lis
  75 #endif
  76 ! variables for table lookup of aerosol impaction/interception scavenging rates
  77  integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
  78  real(r8) dlndg_nimptblgrow
  79  real(r8) scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
  80  real(r8) scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
  81 #endif
  82
  83 contains
  84 #ifndef WRF_PORT
  85   !===============================================================================
  86 #ifdef MODAL_AERO
  87   subroutine mz_aero_defaultopts( aer_wetdep_list_out, aer_drydep_list_out, use_cam_sulfchem_out )
  88 #else
  89   subroutine mz_aero_defaultopts( aer_wetdep_list_out, use_cam_sulfchem_out )
  90 #endif
  91     implicit none
  92     character(len=*), intent(out), optional :: aer_wetdep_list_out(:)
  93 #ifdef MODAL_AERO
  94     character(len=*), intent(out), optional :: aer_drydep_list_out(:)
  95 #endif
  96     logical,          intent(out), optional :: use_cam_sulfchem_out
  97
  98 #ifdef MODAL_AERO
  99     if ( present(aer_drydep_list_out) ) then
 100        aer_drydep_list_out(:) = aer_drydep_list(:)
 101     endif
 102 #endif
 103     if ( present(aer_wetdep_list_out) ) then
 104        aer_wetdep_list_out(:) = aer_wetdep_list(:)
 105     endif
 106     if ( present(use_cam_sulfchem_out) ) then
 107        use_cam_sulfchem_out = use_cam_sulfchem
 108     endif
 109
 110   end subroutine mz_aero_defaultopts
 111
 112   !===============================================================================
 113 #ifdef MODAL_AERO
 114   subroutine mz_aero_setopts( aer_wetdep_list_in, aer_drydep_list_in, use_cam_sulfchem_in )
 115 #else
 116   subroutine mz_aero_setopts( aer_wetdep_list_in, use_cam_sulfchem_in )
 117 #endif
 118     implicit none
 119     character(len=*), intent(in), optional :: aer_wetdep_list_in(:)
 120 #ifdef MODAL_AERO
 121     character(len=*), intent(in), optional :: aer_drydep_list_in(:)
 122 #endif
 123     logical,          intent(in), optional :: use_cam_sulfchem_in
 124
 125 #ifdef MODAL_AERO
 126     if ( present(aer_drydep_list_in) ) then
 127        aer_drydep_list(:) = aer_drydep_list_in(:)
 128     endif
 129 #endif
 130     if ( present(aer_wetdep_list_in) ) then
 131        aer_wetdep_list(:) = aer_wetdep_list_in(:)
 132     endif
 133     if ( present(use_cam_sulfchem_in) ) then
 134        use_cam_sulfchem = use_cam_sulfchem_in
 135     endif
 136
 137   end subroutine mz_aero_setopts
 138 #endif
 139   !===============================================================================
 140   subroutine mz_aero_initialize( )
 141 #ifndef WRF_PORT
 142     use cam_history,      only : addfld, add_default, phys_decomp
 143     use sulchem,          only : inisulchem
 144     use mo_chem_utls,     only : get_inv_ndx
 145     use dust_intr,        only : dst_wet_dep=>dust_has_wet_dep
 146     use progseasalts_intr,only : sst_wet_dep=>progseasalts_has_wet_dep
 147     use dust_intr,        only : dust_names
 148     use progseasalts_intr,only : progseasalts_names
 149     use gas_wetdep_opts,  only : gas_wetdep_list, gas_wetdep_cnt
 150 #else
 151     use module_cam_support, only: addfld, add_default, phys_decomp
 152 #endif
 153
 154     implicit none
 155
 156     integer :: m
 157     integer :: mm
 158     integer :: astat, id
 159
 160     logical :: history_aerosol      ! Output the MAM aerosol tendencies
 161
 162     !-----------------------------------------------------------------------
 163 #ifndef WRF_PORT
 164     call phys_getopts( history_aerosol_out        = history_aerosol   )
 165
 166     dst_wet_dep(:) = .false.
 167     sst_wet_dep(:) = .false.
 168 #else
 169     history_aerosol = .FALSE.
 170 #endif
 171
 172     num_mz_aerosols = 0
 173 #ifdef WRF_PORT
 174     !*NOTE* -only when invoked inside WRF
 175     !Set num_mz_aerosols to some value greater than 0 so that the code inside 'mz_aero_wet_intr'
 176     !subroutine executes
 177     num_mz_aerosols = pcnst
 178 #endif
 179 #ifndef WRF_PORT
 180     count_species: do m = 1,pcnst
 181
 182        if ( len_trim(aer_wetdep_list(m)) == 0 ) exit count_species
 183
 184        do mm=1,gas_wetdep_cnt
 185           if ( trim(gas_wetdep_list(mm)) == trim(aer_wetdep_list(m)) ) then
 186              call endrun('mz_aero_initialize: '//trim(aer_wetdep_list(m))//' is in both het_lst and aer_wetdep_list')
 187           endif
 188        enddo
 189
 190        call cnst_get_ind ( aer_wetdep_list(m), id, abort=.false. )
 191        if ( id < 1 ) then
 192           write(iulog,*) 'mz_aero_initialize: '//trim(aer_wetdep_list(m))//' does not exit in simulation'
 193 #ifdef WRF_PORT
 194           call wrf_message(iulog)
 195 #endif
 196           call endrun('mz_aero_initialize: invalid washout species')
 197        else
 198           where( aer_wetdep_list(m) == dust_names )
 199              dst_wet_dep = .true.
 200           endwhere
 201           where( aer_wetdep_list(m) == progseasalts_names )
 202              sst_wet_dep = .true.
 203           endwhere
 204
 205           if (masterproc) then
 206              write(iulog,*) 'mz_aero_initialize: '//aer_wetdep_list(m)//' will have wet removal'
 207 #ifdef WRF_PORT
 208              call wrf_message(iulog)
 209 #endif
 210           endif
 211           call addfld (trim(aer_wetdep_list(m))//'SFWET','kg/m2/s ',1,  'A','Wet deposition flux at surface',phys_decomp)
 212           call addfld (trim(aer_wetdep_list(m))//'SFSIC','kg/m2/s ', &
 213                 1,  'A','Wet deposition flux (incloud, convective) at surface',phys_decomp)
 214           call addfld (trim(aer_wetdep_list(m))//'SFSIS','kg/m2/s ', &
 215                 1,  'A','Wet deposition flux (incloud, stratiform) at surface',phys_decomp)
 216           call addfld (trim(aer_wetdep_list(m))//'SFSBC','kg/m2/s ', &
 217                 1,  'A','Wet deposition flux (belowcloud, convective) at surface',phys_decomp)
 218           call addfld (trim(aer_wetdep_list(m))//'SFSBS','kg/m2/s ', &
 219                 1,  'A','Wet deposition flux (belowcloud, stratiform) at surface',phys_decomp)
 220           call addfld (trim(aer_wetdep_list(m))//'WET','kg/kg/s ',pver, 'A','wet deposition tendency',phys_decomp)
 221           call addfld (trim(aer_wetdep_list(m))//'SIC','kg/kg/s ',pver, 'A', &
 222                trim(aer_wetdep_list(m))//' ic wet deposition',phys_decomp)
 223           call addfld (trim(aer_wetdep_list(m))//'SIS','kg/kg/s ',pver, 'A', &
 224                trim(aer_wetdep_list(m))//' is wet deposition',phys_decomp)
 225           call addfld (trim(aer_wetdep_list(m))//'SBC','kg/kg/s ',pver, 'A', &
 226                trim(aer_wetdep_list(m))//' bc wet deposition',phys_decomp)
 227           call addfld (trim(aer_wetdep_list(m))//'SBS','kg/kg/s ',pver, 'A', &
 228                trim(aer_wetdep_list(m))//' bs wet deposition',phys_decomp)
 229           if ( history_aerosol ) then
 230              call add_default (trim(aer_wetdep_list(m))//'SFWET', 1, ' ')
 231              call add_default (trim(aer_wetdep_list(m))//'SFSIC', 1, ' ')
 232              call add_default (trim(aer_wetdep_list(m))//'SFSIS', 1, ' ')
 233              call add_default (trim(aer_wetdep_list(m))//'SFSBC', 1, ' ')
 234              call add_default (trim(aer_wetdep_list(m))//'SFSBS', 1, ' ')
 235           endif
 236           num_mz_aerosols = num_mz_aerosols + 1
 237        endif
 238
 239     enddo count_species
 240
 241 #ifdef MODAL_AERO
 242     count_dry_species: do m = 1,pcnst
 243
 244        if ( len_trim(aer_drydep_list(m)) == 0 ) exit count_dry_species
 245
 246        call cnst_get_ind ( aer_drydep_list(m), id, abort=.false. )
 247        if ( id < 1 ) then
 248           write(*,*) 'mz_aero_initialize: '//trim(aer_drydep_list(m))//' does not exist in simulation'
 249           call endrun('mz_aero_initialize: invalid drydep species')
 250        else
 251           if (masterproc) then
 252              write(*,*) 'mz_aero_initialize: '//aer_drydep_list(m)//' will have dry deposition'
 253           endif
 254           call addfld (trim(aer_drydep_list(m))//'DDF','kg/m2/s ',   1, 'A', &
 255                        trim(aer_drydep_list(m))//' dry deposition flux at bottom (grav + turb)',phys_decomp)
 256           call addfld (trim(aer_drydep_list(m))//'TBF','kg/m2/s ',   1, 'A', &
 257                        trim(aer_drydep_list(m))//' turbulent dry deposition flux',phys_decomp)
 258           call addfld (trim(aer_drydep_list(m))//'GVF','kg/m2/s ',   1, 'A', &
 259                        trim(aer_drydep_list(m))//' gravitational dry deposition flux',phys_decomp)
 260           call addfld (trim(aer_drydep_list(m))//'DTQ','kg/kg/s ',pver, 'A', &
 261                        trim(aer_drydep_list(m))//' dry deposition',phys_decomp)
 262           call addfld (trim(aer_drydep_list(m))//'DDV','m/s     ',pver, 'A', &
 263                        trim(aer_drydep_list(m))//' deposition velocity',phys_decomp)
 264           if ( history_aerosol ) then
 265              call add_default (trim(aer_drydep_list(m))//'DDF', 1, ' ')
 266              call add_default (trim(aer_drydep_list(m))//'TBF', 1, ' ')
 267              call add_default (trim(aer_drydep_list(m))//'GVF', 1, ' ')
 268           endif
 269        endif
 270
 271     enddo count_dry_species
 272 #endif
 273
 274     allocate( mz_aerosol_ids(num_mz_aerosols) )
 275
 276     do m = 1, num_mz_aerosols
 277        call cnst_get_ind ( aer_wetdep_list(m), mz_aerosol_ids(m), abort=.false. )
 278     enddo
 279
 280     id_so2  => spc_ids(1)
 281     id_so4  => spc_ids(2)
 282     id_o3   => spc_ids(3)
 283     id_h2o2 => spc_ids(4)
 284     id_oh   => spc_ids(5)
 285     id_no3  => spc_ids(6)
 286     id_ho2  => spc_ids(7)
 287     id_dms  => spc_ids(8)
 288
 289     inv_o3   = get_inv_ndx('O3') > 0
 290     inv_oh   = get_inv_ndx('OH') > 0
 291     inv_no3  = get_inv_ndx('NO3') > 0
 292     inv_ho2  = get_inv_ndx('HO2') > 0
 293     call cnst_get_ind ( 'SO2',  id_so2,  abort=.false. )
 294     call cnst_get_ind ( 'SO4',  id_so4,  abort=.false. )
 295     call cnst_get_ind ( 'H2O2', id_h2o2, abort=.false. )
 296     call cnst_get_ind ( 'DMS',  id_dms,  abort=.false. )
 297
 298     if (inv_o3) then
 299        id_o3 = get_inv_ndx('O3')
 300     else
 301        call cnst_get_ind ( 'O3',   id_o3,   abort=.false. )
 302     endif
 303     if (inv_oh) then
 304        id_oh = get_inv_ndx('OH')
 305     else
 306        call cnst_get_ind ( 'OH',   id_oh,   abort=.false. )
 307     endif
 308     if (inv_no3) then
 309        id_no3 = get_inv_ndx('NO3')
 310     else
 311        call cnst_get_ind ( 'NO3',  id_no3,  abort=.false. )
 312     endif
 313     if (inv_ho2) then
 314        id_ho2 = get_inv_ndx('HO2')
 315     else
 316        call cnst_get_ind ( 'HO2',  id_ho2,  abort=.false. )
 317     endif
 318
 319     do_cam_sulfchem = use_cam_sulfchem .and. all(spc_ids > 0)
 320     if ( do_cam_sulfchem ) then
 321        call inisulchem()
 322     endif
 323
 324     if (masterproc) then
 325        write(iulog,*) 'mz_aero_initialize: do_cam_sulfchem = ',do_cam_sulfchem
 326 #ifdef WRF_PORT
 327        call wrf_message(iulog)
 328 #endif
 329     endif
 330 #endif
 331   end subroutine mz_aero_initialize
 332
 333   !===============================================================================
 334   subroutine mz_aero_wet_intr ( lchnk_in, ncol_in, state_q, state_pdel,          &
 335        state_pmid, state_t, ptend_name, ptend_lq, ptend_q, nstep, dt, cme, prain,    &
 336        evapr, cldv, cldvcu, cldvst, cldc, cldn, fracis, calday, cmfdqr, evapc, conicw, rainmr  &
 337 #if ( defined MODAL_AERO )
 338      , rate1ord_cw2pr_st                                                & ! rce 2010/05/01
 339      , dgncur_awet,  qqcw, qaerwat                                       &
 340 #endif
 341      , cam_out , dlf                                                         )
 342     !-----------------------------------------------------------------------
 343     !-----------------------------------------------------------------------
 344 #ifndef WRF_PORT
 345     use cam_history,   only: outfld
 346     use physics_types, only: physics_state, physics_ptend
 347     use camsrfexch_types, only: cam_out_t
 348 #else
 349     use module_cam_support, only: outfld
 350 #endif
 351     use wetdep,        only: wetdepa_v1, wetdepa_v2
 352 #ifndef WRF_PORT
 353     use sulchem,       only: sulf_chemdr, cldychmini, dicor
 354 #endif
 355     use physconst,     only: gravit
 356     use constituents,  only: cnst_mw
 357     use physconst,     only: amass => mwdry         ! molecular weight dry air ~ kg/kmole
 358     use physconst,     only: boltz                  ! J/K/molecule
 359 #ifndef WRF_PORT
 360     use dust_intr,        only: dust_names
 361     use progseasalts_intr,only: progseasalts_names
 362     use tracer_cnst,   only: get_cnst_data
 363     use chem_mods,     only: fix_mass
 364 #endif
 365 #ifdef MODAL_AERO
 366     use modal_aero_data
 367 #ifndef WRF_PORT
 368     use modal_aero_deposition, only: set_srf_wetdep
 369 #endif
 370 #endif
 371
 372     !-----------------------------------------------------------------------
 373     implicit none
 374     !-----------------------------------------------------------------------
 375     !
 376     ! Arguments:
 377     !
 378     real(r8),            intent(in)  :: dt             ! time step
 379 #ifndef WRF_PORT
 380     type(physics_state), intent(in ) :: state          ! Physics state variables
 381 #else
 382     integer,  intent(in) :: lchnk_in,ncol_in
 383     real(r8), intent(in) :: state_q(pcols,pver,pcnst)
 384     real(r8), intent(in) :: state_pdel(pcols,pver)
 385     real(r8), intent(in) :: state_t(pcols,pver)
 386     real(r8), intent(in) :: state_pmid(pcols,pver)
 387 #endif
 388     integer,  intent(in) :: nstep
 389     real(r8), intent(in) :: cme(pcols,pver)            ! local condensation of cloud water
 390     real(r8), intent(in) :: prain(pcols,pver)          ! production of rain
 391     real(r8), intent(in) :: evapr(pcols,pver)          ! evaporation of rain
 392     real(r8), intent(in) :: cldn(pcols,pver)           ! cloud fraction
 393     real(r8), intent(in) :: cldc(pcols,pver)           ! convective cloud fraction
 394     real(r8), intent(in) :: cldv(pcols,pver)           ! cloudy volume undergoing scavenging
 395     real(r8), intent(in) :: cldvcu(pcols,pver)         ! Convective precipitation area at the top interface of current layer
 396     real(r8), intent(in) :: cldvst(pcols,pver)         ! Stratiform precipitation area at the top interface of current layer
 397     real(r8), intent(in) :: conicw(pcols, pver)
 398     real(r8), intent(in) :: cmfdqr(pcols, pver)
 399     real(r8), intent(in) :: evapc(pcols, pver)         ! Evaporation rate of convective precipitation
 400     real(r8), intent(in) :: rainmr(pcols, pver)         ! rain mixing ratio
 401     real(r8), intent(in) :: dlf(pcols,pver)            ! Detrainment of convective condensate
 402 #ifndef WRF_PORT
 403     type(physics_ptend), intent(inout) :: ptend        ! indivdual parameterization tendencies
 404 #else
 405     character(*), intent(inout) :: ptend_name
 406     logical,      intent(inout) :: ptend_lq(pcnst)
 407     real(r8),     intent(inout) :: ptend_q(pcols,pver,pcnst)
 408     real(r8), target, intent(inout) :: qqcw(pcols,pver,pcnst) ! cloudborne tracer MR array
 409 #endif
 410     real(r8),            intent(inout) :: fracis(pcols,pver,pcnst) ! fraction of transported species that are insoluble
 411     real(r8), intent(in) :: calday        ! current calendar day
 412 #if ( defined MODAL_AERO )
 413     real(r8), intent(in) :: rate1ord_cw2pr_st(pcols,pver) ! 1st order rate for strat cw to precip (1/s) ! rce 2010/05/01
 414     real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode) ! wet/ambient geom. mean diameter (m)
 415                                                                ! for number distribution
 416     real(r8), intent(inout) :: qaerwat(pcols,pver,ntot_amode)
 417 #endif
 418 #ifndef WRF_PORT
 419     type(cam_out_t), intent(inout) :: cam_out
 420 #else
 421     !Balwinder.Singh@pnnl.gov: Temporarily defined as 'real' as it is not currently being used in any computations
 422     real(r8), intent(in) :: cam_out
 423 #endif
 424
 425     !
 426     ! Local variables
 427     !
 428     integer  :: m                                  ! tracer index
 429     integer  :: ixcldice, ixcldliq
 430     integer  :: lchnk                              ! chunk identifier
 431     integer  :: ncol                               ! number of atmospheric columns
 432     real(r8) :: obuf(1)
 433     real(r8) :: iscavt(pcols, pver)
 434     real(r8) :: totcond(pcols, pver) ! sum of cldice and cldliq
 435     integer  :: yr, mon, day, ncsec
 436     integer  :: ncdate
 437     integer  :: mm
 438     integer  :: nphob
 439     real(r8), dimension(pcols,pver) :: h2o23d,o3,oh,no3,ho2,so2,so4,dms,h2o2,ekso2,ekh2o2
 440
 441     real(r8), dimension(pcols,pver) ::&
 442          ph,          &! ph of drops
 443          hion          ! Hydrogen ion concentration in drops
 444
 445     integer ::&
 446          indcp(pcols,pver),    &! Indices of cloudy grid points
 447          ncldypts(pver)         ! Number of cloudy grid points
 448
 449     real(r8) :: icwmr2(pcols,pver) ! in cloud water mixing ratio for hack scheme
 450     real(r8) :: icscavt(pcols, pver)
 451     real(r8) :: isscavt(pcols, pver)
 452     real(r8) :: bcscavt(pcols, pver)
 453     real(r8) :: bsscavt(pcols, pver)
 454     real(r8) :: sol_factb, sol_facti
 455     real(r8) :: sol_factic(pcols,pver)
 456     real(r8) :: sol_factbi, sol_factii, sol_factiic
 457     real(r8) :: xhnm(pcols, pver)
 458     real(r8) :: sflx(pcols)            ! deposition flux
 459     real(r8) :: dicorfac(pcols)        ! factors to increase HO2 diurnal averages
 460     integer :: i,k
 461     real(r8) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1)
 462 #ifdef MODAL_AERO
 463     integer :: jnv                     ! index for scavcoefnv 3rd dimension
 464     integer :: lphase                  ! index for interstitial / cloudborne aerosol
 465     integer :: lspec                   ! index for aerosol number / chem-mass / water-mass
 466     integer :: lcoardust, lcoarnacl    ! indices for coarse mode dust and seasalt masses
 467     real(r8) :: dqdt_tmp(pcols,pver)   ! temporary array to hold tendency for 1 species
 468     real(r8) :: f_act_conv(pcols,pver) ! prescribed aerosol activation fraction for convective cloud   ! rce 2010/05/01
 469     real(r8) :: f_act_conv_coarse(pcols,pver) ! similar but for coarse mode                            ! rce 2010/05/02
 470     real(r8) :: f_act_conv_coarse_dust, f_act_conv_coarse_nacl                                         ! rce 2010/05/02
 471     real(r8) :: fracis_cw(pcols,pver)
 472     real(r8) :: hygro_sum_old(pcols,pver)  ! before removal    [sum of (mass*hydro/dens)]
 473     real(r8) :: hygro_sum_del(pcols,pver)  ! removal change to [sum of (mass*hydro/dens)]
 474     real(r8) :: hygro_sum_old_ik, hygro_sum_new_ik
 475     real(r8) :: prec(pcols)                ! precipitation rate
 476     real(r8) :: q_tmp(pcols,pver)          ! temporary array to hold "most current" mixing ratio for 1 species
 477     real(r8) :: qqcw_tmp(pcols,pver)       ! temporary array to hold qqcw   ! rce 2010/05/01
 478     real(r8) :: scavcoefnv(pcols,pver,0:2) ! Dana and Hales coefficient (/mm) for
 479                                            ! cloud-borne num & vol (0),
 480                                            ! interstitial num (1), interstitial vol (2)
 481     real(r8) :: tmpa, tmpb
 482     real(r8) :: tmpdust, tmpnacl
 483     real(r8) :: water_old, water_new   ! temporary old/new aerosol water mix-rat
 484     logical :: isprx(pcols,pver) ! true if precipation
 485     real(r8) :: aerdepwetis(pcols,pcnst)  ! aerosol wet deposition (interstitial)
 486     real(r8) :: aerdepwetcw(pcols,pcnst)  ! aerosol wet deposition (cloud water)
 487     real(r8), pointer :: fldcw(:,:)
 488 #endif
 489
 490     !-----------------------------------------------------------------------
 491     call cnst_get_ind('CLDICE', ixcldice)
 492     call cnst_get_ind('CLDLIQ', ixcldliq)
 493 #ifndef WRF_PORT
 494     lchnk = state%lchnk
 495     ncol  = state%ncol
 496 #else
 497     lchnk = lchnk_in
 498     ncol  = ncol_in
 499 #endif
 500
 501     totcond(:ncol, :) = state_q(:ncol,:,ixcldliq) + state_q(:ncol,:,ixcldice)
 502
 503     if ( num_mz_aerosols > 0 ) then
 504
 505        ! Wet deposition of mozart aerosol species.
 506        ptend_name  = ptend_name//'+mz_aero_wetdep'
 507
 508 #ifdef MODAL_AERO
 509      prec(:ncol)=0._r8
 510      do k=1,pver
 511         where (prec(:ncol) >= 1.e-7)
 512             isprx(:ncol,k) = .true.
 513         elsewhere
 514             isprx(:ncol,k) = .false.
 515         endwhere
 516         prec(:ncol) = prec(:ncol) + (prain(:ncol,k) + cmfdqr(:ncol,k) - evapr(:ncol,k)) &
 517                     *state_pdel(:ncol,k)/gravit
 518      end do
 519
 520
 521 ! calculate the mass-weighted sol_factic for coarse mode species
 522 !    sol_factic_coarse(:,:) = 0.30_r8   ! tuned 1/4
 523      f_act_conv_coarse(:,:) = 0.60_r8   ! rce 2010/05/02
 524      f_act_conv_coarse_dust = 0.40_r8   ! rce 2010/05/02
 525      f_act_conv_coarse_nacl = 0.80_r8   ! rce 2010/05/02
 526      if (modeptr_coarse > 0) then
 527         lcoardust = lptr_dust_a_amode(modeptr_coarse)
 528         lcoarnacl = lptr_nacl_a_amode(modeptr_coarse)
 529         if ((lcoardust > 0) .and. (lcoarnacl > 0)) then
 530            do k = 1, pver
 531            do i = 1, ncol
 532               tmpdust = max( 0.0_r8, state_q(i,k,lcoardust) + ptend_q(i,k,lcoardust)*dt )
 533               tmpnacl = max( 0.0_r8, state_q(i,k,lcoarnacl) + ptend_q(i,k,lcoarnacl)*dt )
 534               if ((tmpdust+tmpnacl) > 1.0e-30_r8) then
 535 !                sol_factic_coarse(i,k) = (0.2_r8*tmpdust + 0.4_r8*tmpnacl)/(tmpdust+tmpnacl)  ! tuned 1/6
 536                  f_act_conv_coarse(i,k) = (f_act_conv_coarse_dust*tmpdust &
 537                                          + f_act_conv_coarse_nacl*tmpnacl)/(tmpdust+tmpnacl)  ! rce 2010/05/02
 538               end if
 539            end do
 540            end do
 541         end if
 542      end if
 543
 544
 545     scavcoefnv(:,:,0) = 0.0_r8   ! below-cloud scavcoef = 0.0 for cloud-borne species
 546
 547     do m = 1, ntot_amode   ! main loop over aerosol modes
 548
 549        do lphase = 1, 2   ! loop over interstitial (1) and cloud-borne (2) forms
 550
 551 ! sol_factb and sol_facti values
 552 !    sol_factb - currently this is basically a tuning factor
 553 !    sol_facti & sol_factic - currently has a physical basis, and reflects activation fraction
 554 !
 555 ! 2008-mar-07 rce - sol_factb  (interstitial) changed from 0.3 to 0.1
 556 !                 - sol_factic (interstitial, dust modes) changed from 1.0 to 0.5
 557 !                 - sol_factic (cloud-borne, pcarb modes) no need to set it to 0.0
 558 !                       because the cloud-borne pcarbon == 0 (no activation)
 559 !
 560 ! rce 2010/05/02
 561 ! prior to this date, sol_factic was used for convective in-cloud wet removal,
 562 !    and its value reflected a combination of an activation fraction (which varied between modes)
 563 !    and a tuning factor
 564 ! from this date forward, two parameters are used for convective in-cloud wet removal
 565 !    f_act_conv is the activation fraction
 566 !       note that "non-activation" of aerosol in air entrained into updrafts should
 567 !          be included here
 568 !       eventually we might use the activate routine (with w ~= 1 m/s) to calculate
 569 !          this, but there is still the entrainment issue
 570 !    sol_factic is strictly a tuning factor
 571 !
 572           if (lphase == 1) then   ! interstial aerosol
 573              hygro_sum_old(:,:) = 0.0_r8
 574              hygro_sum_del(:,:) = 0.0_r8
 575              call modal_aero_bcscavcoef_get( m, ncol, isprx, dgncur_awet,   &
 576                                              scavcoefnv(:,:,1), scavcoefnv(:,:,2) )
 577
 578              sol_factb  = 0.1_r8   ! all below-cloud scav ON (0.1 "tuning factor")
 579 !            sol_factb  = 0.03_r8   ! all below-cloud scav ON (0.1 "tuning factor")  ! tuned 1/6
 580
 581              sol_facti  = 0.0_r8   ! strat  in-cloud scav totally OFF for institial
 582
 583              sol_factic = 0.4_r8      ! xl 2010/05/20
 584
 585              if (m == modeptr_pcarbon) then
 586 !               sol_factic = 0.0_r8   ! conv   in-cloud scav OFF (0.0 activation fraction)
 587                 f_act_conv = 0.0_r8   ! rce 2010/05/02
 588              else if ((m == modeptr_finedust) .or. (m == modeptr_coardust)) then
 589 !               sol_factic = 0.2_r8   ! conv   in-cloud scav ON  (0.5 activation fraction)  ! tuned 1/4
 590                 f_act_conv = 0.4_r8   ! rce 2010/05/02
 591              else
 592 !               sol_factic = 0.4_r8   ! conv   in-cloud scav ON  (1.0 activation fraction)  ! tuned 1/4
 593                 f_act_conv = 0.8_r8   ! rce 2010/05/02
 594              end if
 595
 596           else   ! cloud-borne aerosol (borne by stratiform cloud drops)
 597              sol_factb  = 0.0_r8   ! all below-cloud scav OFF (anything cloud-borne is located "in-cloud")
 598
 599              sol_facti  = 1.0_r8   ! strat  in-cloud scav totally ON for cloud-borne
 600 !            sol_facti  = 0.6_r8   ! strat  in-cloud scav totally ON for cloud-borne  ! tuned 1/6
 601              sol_factic = 0.0_r8   ! conv   in-cloud scav OFF (having this on would mean
 602                                    !        that conv precip collects strat droplets)
 603              f_act_conv = 0.0_r8   ! conv   in-cloud scav OFF (having this on would mean
 604
 605           end if
 606 !
 607 ! rce 2010/05/03
 608 ! wetdepa has 6 "sol_fact" parameters:
 609 !      sol_facti, sol_factic, sol_factb for liquid cloud
 610 !      sol_factii, sol_factiic, sol_factbi for ice cloud
 611 ! the ice cloud parameters are optional, and if not provided, they default to
 612 !    one of the other sol_fact parameters (see subr. wetdepa about this)
 613 ! for now, we set the ice cloud parameters equal
 614 !    to their liquid cloud counterparts
 615 ! currently the ice parameters are not used in wetdepa as
 616 !    wetdepa sets "weight" (the ice cloud fraction) to 0.0
 617 ! if this changes, we will have to give more thought to
 618 !    the ice cloud parameter values
 619 !
 620           sol_factbi  = sol_factb
 621           sol_factii  = sol_facti
 622           sol_factiic = sol_factic(1,1)
 623
 624
 625           do lspec = 0, nspec_amode(m)+1   ! loop over number + chem constituents + water
 626
 627              if (lspec == 0) then   ! number
 628                 if (lphase == 1) then
 629                    mm = numptr_amode(m)
 630                    jnv = 1
 631                 else
 632                    mm = numptrcw_amode(m)
 633                    jnv = 0
 634                 endif
 635              else if (lspec <= nspec_amode(m)) then   ! non-water mass
 636                 if (lphase == 1) then
 637                    mm = lmassptr_amode(lspec,m)
 638                    jnv = 2
 639                 else
 640                    mm = lmassptrcw_amode(lspec,m)
 641                    jnv = 0
 642                 endif
 643              else   ! water mass
 644 !   bypass wet removal of aerosol water
 645                 cycle
 646                 if (lphase == 1) then
 647                    mm = 0
 648 !                  mm = lwaterptr_amode(m)
 649                    jnv = 2
 650                 else
 651                    mm = 0
 652                    jnv = 0
 653                 endif
 654              endif
 655
 656           if (mm <= 0) cycle
 657
 658
 659 ! set f_act_conv for interstitial (lphase=1) coarse mode species
 660 ! for the convective in-cloud, we conceptually treat the coarse dust and seasalt
 661 !    as being externally mixed, and apply f_act_conv = f_act_conv_coarse_dust/nacl to dust/seasalt
 662 ! number and sulfate are conceptually partitioned to the dust and seasalt
 663 !    on a mass basis, so the f_act_conv for number and sulfate are
 664 !    mass-weighted averages of the values used for dust/seasalt
 665           if ((lphase == 1) .and. (m == modeptr_coarse)) then
 666 !            sol_factic = sol_factic_coarse
 667              f_act_conv = f_act_conv_coarse   ! rce 2010/05/02
 668              if (lspec > 0) then
 669                 if (lmassptr_amode(lspec,m) == lptr_dust_a_amode(m)) then
 670 !                  sol_factic = 0.2_r8   ! tuned 1/4
 671                    f_act_conv = f_act_conv_coarse_dust   ! rce 2010/05/02
 672                 else if (lmassptr_amode(lspec,m) == lptr_nacl_a_amode(m)) then
 673 !                  sol_factic = 0.4_r8   ! tuned 1/6
 674                    f_act_conv = f_act_conv_coarse_nacl   ! rce 2010/05/02
 675                 end if
 676              end if
 677           end if
 678           if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
 679              ptend_lq(mm) = .TRUE.
 680              dqdt_tmp(:,:) = 0.0_r8
 681 ! q_tmp reflects changes from modal_aero_calcsize and is the "most current" q
 682              q_tmp(1:ncol,:) = state_q(1:ncol,:,mm) + ptend_q(1:ncol,:,mm)*dt
 683 #ifndef WRF_PORT
 684              fldcw => qqcw_get_field(mm,lchnk)
 685 #else
 686              fldcw => qqcw(:,:,mm)
 687 #endif
 688              call wetdepa_v2( state_t, state_pmid, state_q, state_pdel,  &
 689                   cldn, cldc, cmfdqr, evapc, conicw, prain, cme,                     &
 690                   evapr, totcond, q_tmp(:,:), dt,            &
 691                   dqdt_tmp(:,:), iscavt, cldv, cldvcu, cldvst, dlf, fracis(:,:,mm), sol_factb, ncol, &
 692                   scavcoefnv(:,:,jnv), &
 693                   .false., rate1ord_cw2pr_st, fldcw, f_act_conv, &     ! rce 2010/05/01
 694                   icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
 695                   sol_facti_in=sol_facti, sol_factbi_in=sol_factbi, sol_factii_in=sol_factii, &   ! rce 2010/05/03
 696                   sol_factic_in=sol_factic, sol_factiic_in=sol_factiic )                          ! rce 2010/05/03
 697
 698              ptend_q(1:ncol,:,mm) = ptend_q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:)
 699
 700              call outfld( trim(cnst_name(mm))//'WET', dqdt_tmp(:,:), pcols, lchnk)
 701              call outfld( trim(cnst_name(mm))//'SIC', icscavt, pcols, lchnk)
 702              call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
 703              call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
 704              call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)
 705
 706              sflx(:)=0._r8
 707              do k=1,pver
 708                 do i=1,ncol
 709                    sflx(i)=sflx(i)+dqdt_tmp(i,k)*state_pdel(i,k)/gravit
 710                 enddo
 711              enddo
 712              call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)
 713              aerdepwetis(:ncol,mm) = sflx(:ncol)
 714
 715              sflx(:)=0._r8
 716              do k=1,pver
 717                 do i=1,ncol
 718                    sflx(i)=sflx(i)+icscavt(i,k)*state_pdel(i,k)/gravit
 719                 enddo
 720              enddo
 721              call outfld( trim(cnst_name(mm))//'SFSIC', sflx, pcols, lchnk)
 722              sflx(:)=0._r8
 723              do k=1,pver
 724                 do i=1,ncol
 725                    sflx(i)=sflx(i)+isscavt(i,k)*state_pdel(i,k)/gravit
 726                 enddo
 727              enddo
 728              call outfld( trim(cnst_name(mm))//'SFSIS', sflx, pcols, lchnk)
 729              sflx(:)=0._r8
 730              do k=1,pver
 731                 do i=1,ncol
 732                    sflx(i)=sflx(i)+bcscavt(i,k)*state_pdel(i,k)/gravit
 733                 enddo
 734              enddo
 735              call outfld( trim(cnst_name(mm))//'SFSBC', sflx, pcols, lchnk)
 736              sflx(:)=0._r8
 737              do k=1,pver
 738                 do i=1,ncol
 739                    sflx(i)=sflx(i)+bsscavt(i,k)*state_pdel(i,k)/gravit
 740                 enddo
 741              enddo
 742              call outfld( trim(cnst_name(mm))//'SFSBS', sflx, pcols, lchnk)
 743
 744              if (lspec > 0) then
 745                 tmpa = spechygro(lspectype_amode(lspec,m))/ &
 746                        specdens_amode(lspectype_amode(lspec,m))
 747                 tmpb = tmpa*dt
 748                 hygro_sum_old(1:ncol,:) = hygro_sum_old(1:ncol,:) &
 749                                          + tmpa*q_tmp(1:ncol,:)
 750                 hygro_sum_del(1:ncol,:) = hygro_sum_del(1:ncol,:) &
 751                                          + tmpb*dqdt_tmp(1:ncol,:)
 752              end if
 753
 754           else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then
 755 ! aerosol water -- because of how wetdepa treats evaporation of stratiform
 756 !    precip, it is not appropriate to apply wetdepa to aerosol water
 757 ! instead, "hygro_sum" = [sum of (mass*hygro/dens)] is calculated before and
 758 !    after wet removal, and new water is calculated using
 759 !    new_water = old_water*min(10,(hygro_sum_new/hygro_sum_old))
 760 ! the "min(10,...)" is to avoid potential problems when hygro_sum_old ~= 0
 761 ! also, individual wet removal terms (ic,is,bc,bs) are not output to history
 762 !            ptend%lq(mm) = .TRUE.
 763 !            dqdt_tmp(:,:) = 0.0_r8
 764              do k = 1, pver
 765              do i = 1, ncol
 766 !               water_old = max( 0.0_r8, state%q(i,k,mm)+ptend%q(i,k,mm)*dt )
 767                 water_old = max( 0.0_r8, qaerwat(i,k,mm) )
 768                 hygro_sum_old_ik = max( 0.0_r8, hygro_sum_old(i,k) )
 769                 hygro_sum_new_ik = max( 0.0_r8, hygro_sum_old_ik+hygro_sum_del(i,k) )
 770                 if (hygro_sum_new_ik >= 10.0_r8*hygro_sum_old_ik) then
 771                    water_new = 10.0_r8*water_old
 772                 else
 773                    water_new = water_old*(hygro_sum_new_ik/hygro_sum_old_ik)
 774                 end if
 775 !               dqdt_tmp(i,k) = (water_new - water_old)/dt
 776                 qaerwat(i,k,mm) = water_new
 777              end do
 778              end do
 779
 780 !            ptend%q(1:ncol,:,mm) = ptend%q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:)
 781
 782 !            call outfld( trim(cnst_name(mm))//'WET', dqdt_tmp(:,:), pcols, lchnk)
 783
 784 !            sflx(:)=0._r8
 785 !            do k=1,pver
 786 !               do i=1,ncol
 787 !                  sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
 788 !               enddo
 789 !            enddo
 790 !            call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)
 791
 792           else   ! lphase == 2
 793              dqdt_tmp(:,:) = 0.0_r8
 794              qqcw_tmp(:,:) = 0.0_r8    ! rce 2010/05/01
 795 #ifndef WRF_PORT
 796              fldcw => qqcw_get_field(mm,lchnk)
 797 #else
 798              fldcw => qqcw(:,:,mm)
 799 #endif
 800              call wetdepa_v2( state_t, state_pmid, state_q, state_pdel,  &
 801                   cldn, cldc, cmfdqr, evapc, conicw, prain, cme,                     &
 802                   evapr, totcond, fldcw, dt,            &
 803                   dqdt_tmp(:,:), iscavt, cldv, cldvcu, cldvst, dlf, fracis_cw(:,:), sol_factb, ncol, &
 804                   scavcoefnv(:,:,jnv), &
 805                   .true., rate1ord_cw2pr_st, qqcw_tmp, f_act_conv, &     ! rce 2010/05/01
 806                   icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
 807                   sol_facti_in=sol_facti, sol_factbi_in=sol_factbi, sol_factii_in=sol_factii, &   ! rce 2010/05/03
 808                   sol_factic_in=sol_factic, sol_factiic_in=sol_factiic )                          ! rce 2010/05/03
 809
 810              fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
 811
 812              sflx(:)=0._r8
 813              do k=1,pver
 814                 do i=1,ncol
 815                    sflx(i)=sflx(i)+dqdt_tmp(i,k)*state_pdel(i,k)/gravit
 816                 enddo
 817              enddo
 818              call outfld( trim(cnst_name_cw(mm))//'SFWET', sflx, pcols, lchnk)
 819              aerdepwetcw(:ncol,mm) = sflx(:ncol)
 820
 821              sflx(:)=0._r8
 822              do k=1,pver
 823                 do i=1,ncol
 824                    sflx(i)=sflx(i)+icscavt(i,k)*state_pdel(i,k)/gravit
 825                 enddo
 826              enddo
 827              call outfld( trim(cnst_name_cw(mm))//'SFSIC', sflx, pcols, lchnk)
 828              sflx(:)=0._r8
 829              do k=1,pver
 830                 do i=1,ncol
 831                    sflx(i)=sflx(i)+isscavt(i,k)*state_pdel(i,k)/gravit
 832                 enddo
 833              enddo
 834              call outfld( trim(cnst_name_cw(mm))//'SFSIS', sflx, pcols, lchnk)
 835              sflx(:)=0._r8
 836              do k=1,pver
 837                 do i=1,ncol
 838                    sflx(i)=sflx(i)+bcscavt(i,k)*state_pdel(i,k)/gravit
 839                 enddo
 840              enddo
 841              call outfld( trim(cnst_name_cw(mm))//'SFSBC', sflx, pcols, lchnk)
 842              sflx(:)=0._r8
 843              do k=1,pver
 844                 do i=1,ncol
 845                    sflx(i)=sflx(i)+bsscavt(i,k)*state_pdel(i,k)/gravit
 846                 enddo
 847              enddo
 848              call outfld( trim(cnst_name_cw(mm))//'SFSBS', sflx, pcols, lchnk)
 849
 850           endif
 851
 852           enddo   ! lspec = 0, nspec_amode(m)+1
 853        enddo   ! lphase = 1, 2
 854     enddo   ! m = 1, ntot_amode
 855 #ifndef WRF_PORT
 856     call set_srf_wetdep(aerdepwetis, aerdepwetcw, cam_out)
 857 #endif
 858
 859 #else
 860        scavcoef(:ncol,:)=0.1_r8
 861        do m = 1, num_mz_aerosols
 862
 863           mm = mz_aerosol_ids(m)
 864
 865           if (mm > 0 .and. (.not. any(dust_names==cnst_name(mm)))  .and. (.not. any(progseasalts_names==cnst_name(mm))) ) then
 866              ptend%lq(mm) = .TRUE.
 867
 868              sol_factb = 0.3_r8
 869              sol_facti = 0.3_r8
 870              call wetdepa_v1( state%t, state%pmid, state%q, state%pdel,  &
 871                   cldn, cldc, cmfdqr, conicw, prain, cme,                     &
 872                   evapr, totcond, state%q(:,:,mm), dt,            &
 873                   ptend%q(:,:,mm), iscavt, cldv, fracis(:,:,mm), sol_factb, ncol, &
 874                   scavcoef, icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
 875                   sol_facti_in=sol_facti)
 876
 877              call outfld( trim(cnst_name(mm))//'WET', ptend%q(:,:,mm), pcols, lchnk)
 878              call outfld( trim(cnst_name(mm))//'SIC', icscavt , pcols, lchnk)
 879              call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
 880              call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
 881              call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)
 882
 883              sflx(:)=0._r8
 884              do k=1,pver
 885                 do i=1,ncol
 886                    sflx(i)=sflx(i)+ptend%q(i,k,mm)*state%pdel(i,k)/gravit
 887                 enddo
 888              enddo
 889              call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)
 890
 891              ! export deposition fluxes to coupler
 892              select case (trim(cnst_name(mm)))
 893              case('CB2')
 894                 do i = 1, ncol
 895                    cam_out%bcphiwet(i) = max(-sflx(i), 0._r8)
 896                 end do
 897              case('OC2')
 898                 do i = 1, ncol
 899                    cam_out%ocphiwet(i) = max(-sflx(i), 0._r8)
 900                 end do
 901              end select
 902
 903           endif
 904
 905        end do
 906 #endif
 907     endif
 908 #ifndef WRF_PORT
 909     if ( do_cam_sulfchem ) then
 910
 911        so4(:ncol,:)  = state%q(:ncol,:,id_so4) + ptend%q(:ncol,:,id_so4)*dt
 912        so2(:ncol,:)  = state%q(:ncol,:,id_so2) + ptend%q(:ncol,:,id_so2)*dt
 913        dms(:ncol,:)  = state%q(:ncol,:,id_dms)
 914        h2o2(:ncol,:) = state%q(:ncol,:,id_h2o2)
 915
 916        ! need oh,ho2,no3 in molec/cm3, o3 in kg/kg
 917
 918        if ( inv_o3 ) then
 919          call get_cnst_data( 'O3', o3(:ncol,:),  ncol, lchnk ) ! vmr
 920        else
 921          o3(:ncol,:) = state%q(:ncol,:,id_o3)*amass/fix_mass(id_o3) ! vmr
 922        endif
 923
 924        ! vmr -> molec/cm3
 925        xhnm(:ncol,:) = 1.e-6_r8 * state%pmiddry(:ncol,:) / (boltz*state%t(:ncol,:))
 926
 927        if ( inv_oh ) then
 928          call get_cnst_data( 'OH', oh(:ncol,:),  ncol, lchnk ) ! vmr
 929          oh(:ncol,:) = oh(:ncol,:) * xhnm(:ncol,:)             ! molec/cm3
 930        else
 931          oh(:ncol,:) = state%q(:ncol,:,id_oh)* (amass/cnst_mw(id_oh )) * xhnm(:ncol,:)
 932        endif
 933        if ( inv_no3 ) then
 934          call get_cnst_data( 'NO3', no3(:ncol,:),  ncol, lchnk ) ! vmr
 935          no3(:ncol,:) = no3(:ncol,:) * xhnm(:ncol,:)             ! molec/cm3
 936        else
 937          no3(:ncol,:) = state%q(:ncol,:,id_no3)* (amass/cnst_mw(id_no3)) * xhnm(:ncol,:)
 938        endif
 939        if ( inv_ho2 ) then
 940          call get_cnst_data( 'HO2', ho2(:ncol,:),  ncol, lchnk ) ! vmr
 941          ho2(:ncol,:) = ho2(:ncol,:) * xhnm(:ncol,:)             ! molec/cm3
 942        else
 943          ho2(:ncol,:) = state%q(:ncol,:,id_ho2)* (amass/cnst_mw(id_ho2)) * xhnm(:ncol,:)
 944        endif
 945
 946        ! alter HO2 from diurnal average to instantaneous value
 947        call dicor( calday, state%lat, ncol, dicorfac)
 948        do k = 1, pver
 949           do i = 1, ncol
 950              ho2(i,k) = ho2(i,k)*dicorfac(i)
 951           end do
 952        end do
 953
 954        icwmr2(:,:) = 0._r8
 955        call cldychmini( state%t,  state%pmid, totcond, rainmr, cldn, &
 956             cldv, conicw, icwmr2, so2, so4, &
 957             ph, hion, ekso2, ekh2o2, indcp, &
 958             ncldypts, ncol )
 959        call outfld( 'PH', ph, pcols, lchnk )
 960
 961        call sulf_chemdr( state%t, state%pmid, totcond, rainmr, cldn, cldv, &
 962                          conicw, icwmr2, indcp, ncldypts, hion, &
 963                          so2, so4, dms,  h2o2, ekso2, ekh2o2, &
 964                          o3, h2o23d, oh, no3, ho2, state%q(:,:,1),  &
 965                          state%zm, dt, calday, state%lat, state%lon, ncol, lchnk  )
 966
 967        ptend%q(:ncol,:,id_so2) = (so2(:ncol,:)-state%q(:ncol,:,id_so2))/dt
 968        ptend%lq(id_so2) = .true.
 969
 970        ptend%q(:ncol,:,id_so4) = (so4(:ncol,:)-state%q(:ncol,:,id_so4))/dt
 971        ptend%lq(id_so4) = .true.
 972
 973        ptend%q(:ncol,:,id_dms) = (dms(:ncol,:)-state%q(:ncol,:,id_dms))/dt
 974        ptend%lq(id_dms) = .true.
 975
 976        ptend%q(:ncol,:,id_h2o2) = (h2o2(:ncol,:)-state%q(:ncol,:,id_h2o2))/dt
 977        ptend%lq(id_h2o2) = .true.
 978
 979     endif
 980 #endif
 981     return
 982
 983   end subroutine mz_aero_wet_intr
 984
 985
 986 #ifdef MODAL_AERO
 987 #ifndef WRF_PORT
 988   !===============================================================================
 989   subroutine mz_aero_dry_intr ( state, ptend, wvflx, dt, lat, clat, &
 990        fsds, obklen, ts, ustar, prect, snowh, pblh, hflx, month, landfrac, &
 991        icefrac, ocnfrac, fvin,ram1in, dgncur_awet, wetdens,  qaerwat, &
 992        cam_out)
 993     !-----------------------------------------------------------------------
 994     !
 995     ! Purpose:
 996     ! Interface to dry deposition and sedimentation of aerosols
 997     !
 998     ! Author: Xiaohong Liu and Dick Easter
 999     !
1000     !-----------------------------------------------------------------------
1001     use cam_history,       only: outfld
1002     use ppgrid,            only: pverp
1003     use physics_types,     only: physics_state, physics_ptend
1004     use camsrfexch_types,  only: cam_out_t
1005     use physconst,         only: gravit, rair, rhoh2o
1006     use drydep_mod,        only: setdvel,  d3ddflux, calcram
1007     use dust_sediment_mod, only: dust_sediment_tend, dust_sediment_vel
1008     use modal_aero_data
1009     use modal_aero_deposition, only: set_srf_drydep
1010     !-----------------------------------------------------------------------
1011     implicit none
1012     !-----------------------------------------------------------------------
1013     !
1014     ! Arguments:
1015     !
1016     real(r8),            intent(in)  :: dt             ! time step
1017     type(physics_state), intent(in ) :: state          ! Physics state variables
1018     integer, intent(in) :: lat(pcols)                  ! latitude index for S->N storage
1019     real(r8), intent(in) :: clat(pcols)                 ! latitude
1020     real(r8), intent(in) :: fsds(pcols)                 ! longwave down at sfc
1021     real(r8), intent(in) :: obklen(pcols)                 ! obklen
1022     real(r8), intent(in) :: ustar(pcols)                  ! sfc fric vel--used over oceans and sea ice.
1023     real(r8), intent(in) :: ts(pcols)                     ! sfc temp
1024     real(r8), intent(in) :: landfrac(pcols)               ! land fraction
1025     real(r8), intent(in) :: icefrac(pcols)                ! ice fraction
1026     real(r8), intent(in) :: ocnfrac(pcols)                ! ocean fraction
1027     real(r8), intent(in) :: hflx(pcols)                  ! sensible heat flux
1028     real(r8), intent(in) :: prect(pcols)                     ! prect
1029     real(r8), intent(in) :: snowh(pcols)                     ! snow depth
1030     real(r8), intent(in) :: pblh(pcols)                     ! pbl height
1031     integer, intent(in)  :: month
1032     real(r8), intent(in) :: wvflx(pcols)       ! water vapor flux
1033     real(r8), intent(in) :: fvin(pcols)        ! for dry dep velocities from land model
1034     real(r8), intent(in) :: ram1in(pcols)       ! for dry dep velocities from land model
1035     real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
1036     real(r8), intent(in) :: wetdens(pcols,pver,ntot_amode)
1037     real(r8), intent(inout) :: qaerwat(pcols,pver,ntot_amode)
1038
1039     type(physics_ptend), intent(inout) :: ptend          ! indivdual parameterization tendencies
1040     type(cam_out_t),     intent(inout) :: cam_out
1041     !
1042     ! Local variables
1043     !
1044     integer :: i,k
1045     integer :: jvlc                    ! index for last dimension of vlc_xxx arrays
1046     integer :: lphase                  ! index for interstitial / cloudborne aerosol
1047     integer :: lspec                   ! index for aerosol number / chem-mass / water-mass
1048     integer :: m                       ! aerosol mode index
1049     integer :: mm                      ! tracer index
1050     integer :: lchnk                   ! chunk identifier
1051     integer :: ncol                    ! number of atmospheric columns
1052     real(r8) :: tvs(pcols,pver)
1053     real(r8) :: sflx(pcols)            ! deposition flux
1054     real(r8)::  dep_trb(pcols)       !kg/m2/s
1055     real(r8)::  dep_dry(pcols)       !kg/m2/s (total of grav and trb)
1056     real(r8)::  dep_grv(pcols)       !kg/m2/s (total of grav and trb)
1057     real (r8) :: rho(pcols,pver)                    ! air density in kg/m3
1058     real(r8) :: pvmzaer(pcols,pverp)    ! sedimentation velocity in Pa
1059     real(r8) :: fv(pcols)         ! for dry dep velocities, from land modified over ocean & ice
1060     real(r8) :: ram1(pcols)       ! for dry dep velocities, from land modified over ocean & ice
1061     real(r8) :: dqdt_tmp(pcols,pver)   ! temporary array to hold tendency for 1 species
1062
1063     real(r8) :: rad_drop(pcols,pver)
1064     real(r8) :: dens_drop(pcols,pver)
1065     real(r8) :: sg_drop(pcols,pver)
1066     real(r8) :: rad_aer(pcols,pver)
1067     real(r8) :: dens_aer(pcols,pver)
1068     real(r8) :: sg_aer(pcols,pver)
1069
1070     real(r8) :: vlc_dry(pcols,pver,4)       ! dep velocity
1071     real(r8) :: vlc_grv(pcols,pver,4)       ! dep velocity
1072     real(r8)::  vlc_trb(pcols,4)            ! dep velocity
1073     real(r8) :: aerdepdryis(pcols,pcnst)  ! aerosol dry deposition (interstitial)
1074     real(r8) :: aerdepdrycw(pcols,pcnst)  ! aerosol dry deposition (cloud water)
1075     real(r8), pointer :: fldcw(:,:)
1076     !-----------------------------------------------------------------------
1077
1078     lchnk = state%lchnk
1079     ncol  = state%ncol
1080     tvs(:ncol,:) = state%t(:ncol,:)!*(1+state%q(:ncol,k)
1081     rho(:ncol,:)=  state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
1082
1083     ! Dry deposition of mozart aerosol species.
1084     ptend%name  = ptend%name//'+mz_aero_drydep'
1085
1086     !   #################################
1087     !    call setdvel( ncol, landfrac, icefrac, ocnfrac, .001_r8, .001_r8, .001_r8, dvel )
1088     !  we get the ram1,fv from the land model as ram1in,fvin, but need to calculate it over oceans and ice.
1089     !  better if we got thse from the ocean and ice model
1090     !  for friction velocity, we use ustar (from taux and tauy), except over land, where we use fv from the land model.
1091
1092     call calcram(ncol,landfrac,icefrac,ocnfrac,obklen,&
1093          ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
1094          state%pdel(:,pver),fvin,fv)
1095 !   call outfld( 'airFV', fv(:), pcols, lchnk )
1096 !   call outfld( 'RAM1', ram1(:), pcols, lchnk )
1097     !       call outfld( 'icefrc', icefrac(:), pcols, lchnk )
1098
1099
1100 !
1101 ! calc settling/deposition velocities for cloud droplets (and cloud-borne aerosols)
1102 !
1103 ! *** mean drop radius should eventually be computed from ndrop and qcldwtr
1104     rad_drop(:,:) = 5.0e-6_r8
1105     dens_drop(:,:) = rhoh2o
1106     sg_drop(:,:) = 1.46
1107     jvlc = 3
1108     call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv,  &
1109                      vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
1110                      rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 0, lchnk)
1111     jvlc = 4
1112     call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv,  &
1113                      vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
1114                      rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 3, lchnk)
1115
1116
1117
1118     do m = 1, ntot_amode   ! main loop over aerosol modes
1119
1120        do lphase = 1, 2   ! loop over interstitial / cloud-borne forms
1121
1122           if (lphase == 1) then   ! interstial aerosol - calc settling/dep velocities of mode
1123
1124 ! rad_aer = volume mean wet radius (m)
1125 ! dgncur_awet = geometric mean wet diameter for number distribution (m)
1126              rad_aer(1:ncol,:) = 0.5_r8*dgncur_awet(1:ncol,:,m)   &
1127                                  *exp(1.5*(alnsg_amode(m)**2))
1128 ! dens_aer(1:ncol,:) = wet density (kg/m3)
1129              dens_aer(1:ncol,:) = wetdens(1:ncol,:,m)
1130              sg_aer(1:ncol,:) = sigmag_amode(m)
1131
1132              jvlc = 1
1133              call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv,  &
1134                         vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
1135                         rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 0, lchnk)
1136              jvlc = 2
1137              call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv,  &
1138                         vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
1139                         rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 3, lchnk)
1140           end if
1141
1142           do lspec = 0, nspec_amode(m)+1   ! loop over number + constituents + water
1143
1144              if (lspec == 0) then   ! number
1145                 if (lphase == 1) then
1146                    mm = numptr_amode(m)
1147                    jvlc = 1
1148                 else
1149                    mm = numptrcw_amode(m)
1150                    jvlc = 3
1151                 endif
1152              else if (lspec <= nspec_amode(m)) then   ! non-water mass
1153                 if (lphase == 1) then
1154                    mm = lmassptr_amode(lspec,m)
1155                    jvlc = 2
1156                 else
1157                    mm = lmassptrcw_amode(lspec,m)
1158                    jvlc = 4
1159                 endif
1160              else   ! water mass
1161 !   bypass dry deposition of aerosol water
1162                 cycle
1163                 if (lphase == 1) then
1164                    mm = 0
1165 !                  mm = lwaterptr_amode(m)
1166                    jvlc = 2
1167                 else
1168                    mm = 0
1169                    jvlc = 4
1170                 endif
1171              endif
1172
1173
1174           if (mm <= 0) cycle
1175
1176 !         if (lphase == 1) then
1177           if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
1178              ptend%lq(mm) = .TRUE.
1179
1180              ! use pvprogseasalts instead (means making the top level 0)
1181              pvmzaer(:ncol,1)=0._r8
1182              pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
1183
1184              call outfld( trim(cnst_name(mm))//'DDV', pvmzaer(:,2:pverp), pcols, lchnk )
1185
1186              if(.true.) then ! use phil's method
1187              !      convert from meters/sec to pascals/sec
1188              !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
1189                 pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
1190
1191              !      calculate the tendencies and sfc fluxes from the above velocities
1192                 call dust_sediment_tend( &
1193                      ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
1194                      state%q(:,:,mm),  pvmzaer,  ptend%q(:,:,mm), sflx  )
1195              else   !use charlie's method
1196                 call d3ddflux( ncol, vlc_dry(:,:,jvlc), state%q(:,:,mm), state%pmid, &
1197                                state%pdel, tvs, sflx, ptend%q(:,:,mm), dt )
1198              endif
1199
1200              ! apportion dry deposition into turb and gravitational settling for tapes
1201              do i=1,ncol
1202                 dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
1203                 dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
1204              enddo
1205
1206              call outfld( trim(cnst_name(mm))//'DDF', sflx, pcols, lchnk)
1207              call outfld( trim(cnst_name(mm))//'TBF', dep_trb, pcols, lchnk )
1208              call outfld( trim(cnst_name(mm))//'GVF', dep_grv, pcols, lchnk )
1209              call outfld( trim(cnst_name(mm))//'DTQ', ptend%q(:,:,mm), pcols, lchnk)
1210              aerdepdryis(:ncol,mm) = sflx(:ncol)
1211
1212           else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then  ! aerosol water
1213              ! use pvprogseasalts instead (means making the top level 0)
1214              pvmzaer(:ncol,1)=0._r8
1215              pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
1216
1217              if(.true.) then ! use phil's method
1218              !      convert from meters/sec to pascals/sec
1219              !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
1220                 pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
1221
1222              !      calculate the tendencies and sfc fluxes from the above velocities
1223                 call dust_sediment_tend( &
1224                      ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
1225                      qaerwat(:,:,mm),  pvmzaer,  dqdt_tmp(:,:), sflx  )
1226              else   !use charlie's method
1227                 call d3ddflux( ncol, vlc_dry(:,:,jvlc), qaerwat(:,:,mm), state%pmid, &
1228                                state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
1229              endif
1230
1231              ! apportion dry deposition into turb and gravitational settling for tapes
1232              do i=1,ncol
1233                 dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
1234                 dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
1235              enddo
1236
1237              qaerwat(1:ncol,:,mm) = qaerwat(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) * dt
1238
1239           else  ! lphase == 2
1240              ! use pvprogseasalts instead (means making the top level 0)
1241              pvmzaer(:ncol,1)=0._r8
1242              pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
1243              fldcw => qqcw_get_field(mm,lchnk)
1244              if(.true.) then ! use phil's method
1245              !      convert from meters/sec to pascals/sec
1246              !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
1247                 pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
1248
1249              !      calculate the tendencies and sfc fluxes from the above velocities
1250                 call dust_sediment_tend( &
1251                      ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
1252                      fldcw(:,:),  pvmzaer,  dqdt_tmp(:,:), sflx  )
1253              else   !use charlie's method
1254                 call d3ddflux( ncol, vlc_dry(:,:,jvlc), fldcw(:,:), state%pmid, &
1255                                state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
1256              endif
1257
1258              ! apportion dry deposition into turb and gravitational settling for tapes
1259              do i=1,ncol
1260                 dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
1261                 dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
1262              enddo
1263
1264              fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
1265
1266              call outfld( trim(cnst_name_cw(mm))//'DDF', sflx, pcols, lchnk)
1267              call outfld( trim(cnst_name_cw(mm))//'TBF', dep_trb, pcols, lchnk )
1268              call outfld( trim(cnst_name_cw(mm))//'GVF', dep_grv, pcols, lchnk )
1269              aerdepdrycw(:ncol,mm) = sflx(:ncol)
1270
1271           endif
1272
1273           enddo   ! lspec = 0, nspec_amode(m)+1
1274        enddo   ! lphase = 1, 2
1275     enddo   ! m = 1, ntot_amode
1276
1277     call set_srf_drydep(aerdepdryis, aerdepdrycw, cam_out)
1278
1279     return
1280   end subroutine mz_aero_dry_intr
1281
1282
1283   !===============================================================================
1284   subroutine modal_aero_depvel_part( ncol, t, pmid, ram1, fv, vlc_dry, vlc_trb, vlc_grv,  &
1285                                      radius_part, density_part, sig_part, moment, lchnk )
1286
1287 !    calculates surface deposition velocity of particles
1288 !    L. Zhang, S. Gong, J. Padro, and L. Barrie
1289 !    A size-seggregated particle dry deposition scheme for an atmospheric aerosol module
1290 !    Atmospheric Environment, 35, 549-560, 2001.
1291 !
1292 !    Authors: X. Liu
1293
1294     !
1295     ! !USES
1296     !
1297     use physconst,     only: pi,boltz, gravit, rair
1298     use mo_drydep,     only: n_land_type, fraction_landuse
1299
1300     ! !ARGUMENTS:
1301     !
1302     implicit none
1303     !
1304     real(r8), intent(in) :: t(pcols,pver)       !atm temperature (K)
1305     real(r8), intent(in) :: pmid(pcols,pver)    !atm pressure (Pa)
1306     real(r8), intent(in) :: fv(pcols)           !friction velocity (m/s)
1307     real(r8), intent(in) :: ram1(pcols)         !aerodynamical resistance (s/m)
1308     real(r8), intent(in) :: radius_part(pcols,pver)    ! mean (volume/number) particle radius (m)
1309     real(r8), intent(in) :: density_part(pcols,pver)   ! density of particle material (kg/m3)
1310     real(r8), intent(in) :: sig_part(pcols,pver)       ! geometric standard deviation of particles
1311     integer,  intent(in) :: moment ! moment of size distribution (0 for number, 2 for surface area, 3 for volume)
1312     integer,  intent(in) :: ncol
1313     integer,  intent(in) :: lchnk
1314
1315     real(r8), intent(out) :: vlc_trb(pcols)       !Turbulent deposn velocity (m/s)
1316     real(r8), intent(out) :: vlc_grv(pcols,pver)       !grav deposn velocity (m/s)
1317     real(r8), intent(out) :: vlc_dry(pcols,pver)       !dry deposn velocity (m/s)
1318     !------------------------------------------------------------------------
1319
1320     !------------------------------------------------------------------------
1321     ! Local Variables
1322     integer  :: m,i,k,ix                !indices
1323     real(r8) :: rho     !atm density (kg/m**3)
1324     real(r8) :: vsc_dyn_atm(pcols,pver)   ![kg m-1 s-1] Dynamic viscosity of air
1325     real(r8) :: vsc_knm_atm(pcols,pver)   ![m2 s-1] Kinematic viscosity of atmosphere
1326     real(r8) :: shm_nbr       ![frc] Schmidt number
1327     real(r8) :: stk_nbr       ![frc] Stokes number
1328     real(r8) :: mfp_atm(pcols,pver)       ![m] Mean free path of air
1329     real(r8) :: dff_aer       ![m2 s-1] Brownian diffusivity of particle
1330     real(r8) :: slp_crc(pcols,pver) ![frc] Slip correction factor
1331     real(r8) :: rss_trb       ![s m-1] Resistance to turbulent deposition
1332     real(r8) :: rss_lmn       ![s m-1] Quasi-laminar layer resistance
1333     real(r8) :: brownian      ! collection efficiency for Browning diffusion
1334     real(r8) :: impaction     ! collection efficiency for impaction
1335     real(r8) :: interception  ! collection efficiency for interception
1336     real(r8) :: stickfrac     ! fraction of particles sticking to surface
1337     real(r8) :: radius_moment(pcols,pver) ! median radius (m) for moment
1338     real(r8) :: lnsig         ! ln(sig_part)
1339     real(r8) :: dispersion    ! accounts for influence of size dist dispersion on bulk settling velocity
1340                               ! assuming radius_part is number mode radius * exp(1.5 ln(sigma))
1341
1342     integer  :: lt
1343     real(r8) :: lnd_frc
1344     real(r8) :: wrk1, wrk2, wrk3
1345
1346     ! constants
1347     real(r8) gamma(11)      ! exponent of schmidt number
1348 !   data gamma/0.54d+00,  0.56d+00,  0.57d+00,  0.54d+00,  0.54d+00, &
1349 !              0.56d+00,  0.54d+00,  0.54d+00,  0.54d+00,  0.56d+00, &
1350 !              0.50d+00/
1351     data gamma/0.56d+00,  0.54d+00,  0.54d+00,  0.56d+00,  0.56d+00, &
1352                0.56d+00,  0.50d+00,  0.54d+00,  0.54d+00,  0.54d+00, &
1353                0.54d+00/
1354     save gamma
1355
1356     real(r8) alpha(11)      ! parameter for impaction
1357 !   data alpha/50.00d+00,  0.95d+00,  0.80d+00,  1.20d+00,  1.30d+00, &
1358 !               0.80d+00, 50.00d+00, 50.00d+00,  2.00d+00,  1.50d+00, &
1359 !             100.00d+00/
1360     data alpha/1.50d+00,   1.20d+00,  1.20d+00,  0.80d+00,  1.00d+00, &
1361                0.80d+00, 100.00d+00, 50.00d+00,  2.00d+00,  1.20d+00, &
1362               50.00d+00/
1363     save alpha
1364
1365     real(r8) radius_collector(11) ! radius (m) of surface collectors
1366 !   data radius_collector/-1.00d+00,  5.10d-03,  3.50d-03,  3.20d-03, 10.00d-03, &
1367 !                          5.00d-03, -1.00d+00, -1.00d+00, 10.00d-03, 10.00d-03, &
1368 !                         -1.00d+00/
1369     data radius_collector/10.00d-03,  3.50d-03,  3.50d-03,  5.10d-03,  2.00d-03, &
1370                            5.00d-03, -1.00d+00, -1.00d+00, 10.00d-03,  3.50d-03, &
1371                           -1.00d+00/
1372     save radius_collector
1373
1374     integer            :: iwet(11) ! flag for wet surface = 1, otherwise = -1
1375 !   data iwet/1,   -1,   -1,   -1,   -1,  &
1376 !            -1,   -1,   -1,    1,   -1,  &
1377 !             1/
1378     data iwet/-1,  -1,   -1,   -1,   -1,  &
1379               -1,   1,   -1,    1,   -1,  &
1380               -1/
1381     save iwet
1382
1383
1384     !------------------------------------------------------------------------
1385     do k=1,pver
1386        do i=1,ncol
1387
1388           lnsig = log(sig_part(i,k))
1389 ! use a maximum radius of 50 microns when calculating deposition velocity
1390           radius_moment(i,k) = min(50.0d-6,radius_part(i,k))*   &
1391                           exp((float(moment)-1.5)*lnsig*lnsig)
1392           dispersion = exp(2.*lnsig*lnsig)
1393
1394           rho=pmid(i,k)/rair/t(i,k)
1395
1396           ! Quasi-laminar layer resistance: call rss_lmn_get
1397           ! Size-independent thermokinetic properties
1398           vsc_dyn_atm(i,k) = 1.72e-5_r8 * ((t(i,k)/273.0_r8)**1.5_r8) * 393.0_r8 / &
1399                (t(i,k)+120.0_r8)      ![kg m-1 s-1] RoY94 p. 102
1400           mfp_atm(i,k) = 2.0_r8 * vsc_dyn_atm(i,k) / &   ![m] SeP97 p. 455
1401                (pmid(i,k)*sqrt(8.0_r8/(pi*rair*t(i,k))))
1402           vsc_knm_atm(i,k) = vsc_dyn_atm(i,k) / rho ![m2 s-1] Kinematic viscosity of air
1403
1404           slp_crc(i,k) = 1.0_r8 + mfp_atm(i,k) * &
1405                   (1.257_r8+0.4_r8*exp(-1.1_r8*radius_moment(i,k)/(mfp_atm(i,k)))) / &
1406                   radius_moment(i,k)   ![frc] Slip correction factor SeP97 p. 464
1407           vlc_grv(i,k) = (4.0_r8/18.0_r8) * radius_moment(i,k)*radius_moment(i,k)*density_part(i,k)* &
1408                   gravit*slp_crc(i,k) / vsc_dyn_atm(i,k) ![m s-1] Stokes' settling velocity SeP97 p. 466
1409           vlc_grv(i,k) = vlc_grv(i,k) * dispersion
1410
1411           vlc_dry(i,k)=vlc_grv(i,k)
1412        enddo
1413     enddo
1414     k=pver  ! only look at bottom level for next part
1415     do i=1,ncol
1416        dff_aer = boltz * t(i,k) * slp_crc(i,k) / &    ![m2 s-1]
1417                  (6.0_r8*pi*vsc_dyn_atm(i,k)*radius_moment(i,k)) !SeP97 p.474
1418        shm_nbr = vsc_knm_atm(i,k) / dff_aer                        ![frc] SeP97 p.972
1419
1420        wrk2 = 0._r8
1421        wrk3 = 0._r8
1422        do lt = 1,n_land_type
1423           lnd_frc = fraction_landuse(i,lt,lchnk)
1424           if ( lnd_frc /= 0._r8 ) then
1425              brownian = shm_nbr**(-gamma(lt))
1426              if (radius_collector(lt) > 0.0d0) then
1427 !       vegetated surface
1428                 stk_nbr = vlc_grv(i,k) * fv(i) / (gravit*radius_collector(lt))
1429                 interception = 2.0_r8*(radius_moment(i,k)/radius_collector(lt))**2.0_r8
1430              else
1431 !       non-vegetated surface
1432                 stk_nbr = vlc_grv(i,k) * fv(i) * fv(i) / (gravit*vsc_knm_atm(i,k))  ![frc] SeP97 p.965
1433                 interception = 0.0_r8
1434              endif
1435              impaction = (stk_nbr/(alpha(lt)+stk_nbr))**2.0_r8
1436
1437              if (iwet(lt) > 0) then
1438                 stickfrac = 1.0_r8
1439              else
1440                 stickfrac = exp(-sqrt(stk_nbr))
1441                 if (stickfrac < 1.0d-10) stickfrac = 1.0d-10
1442              endif
1443              rss_lmn = 1.0_r8 / (3.0_r8 * fv(i) * stickfrac * (brownian+interception+impaction))
1444              rss_trb = ram1(i) + rss_lmn + ram1(i)*rss_lmn*vlc_grv(i,k)
1445
1446              wrk1 = 1.0_r8 / rss_trb
1447              wrk2 = wrk2 + lnd_frc*( wrk1 )
1448              wrk3 = wrk3 + lnd_frc*( wrk1 + vlc_grv(i,k) )
1449           endif
1450        enddo  ! n_land_type
1451        vlc_trb(i) = wrk2
1452        vlc_dry(i,k) = wrk3
1453     enddo !ncol
1454
1455     return
1456   end subroutine modal_aero_depvel_part
1457
1458 #endif
1459 !(WRF_PORT end)
1460   !===============================================================================
1461   subroutine modal_aero_bcscavcoef_get( m, ncol, isprx, dgn_awet, scavcoefnum, scavcoefvol )
1462
1463     use modal_aero_data
1464     !-----------------------------------------------------------------------
1465     implicit none
1466
1467           integer,intent(in) :: m, ncol
1468           logical,intent(in):: isprx(pcols,pver)
1469           real(r8), intent(in) :: dgn_awet(pcols,pver,ntot_amode)
1470           real(r8), intent(out) :: scavcoefnum(pcols,pver), scavcoefvol(pcols,pver)
1471
1472           integer i, k, jgrow
1473           real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum
1474
1475
1476       do k = 1, pver
1477       do i = 1, ncol
1478
1479 !        do only if no precip
1480          if ( isprx(i,k)  ) then
1481 !
1482 !           interpolate table values using log of (actual-wet-size)/(base-dry-size)
1483
1484             dumdgratio = dgn_awet(i,k,m)/dgnum_amode(m)
1485
1486             if ((dumdgratio .ge. 0.99) .and. (dumdgratio .le. 1.01)) then
1487                 scavimpvol = scavimptblvol(0,m)
1488                 scavimpnum = scavimptblnum(0,m)
1489             else
1490                 xgrow = log( dumdgratio ) / dlndg_nimptblgrow
1491                 jgrow = int( xgrow )
1492                 if (xgrow .lt. 0.) jgrow = jgrow - 1
1493                 if (jgrow .lt. nimptblgrow_mind) then
1494                     jgrow = nimptblgrow_mind
1495                     xgrow = jgrow
1496                 else
1497                     jgrow = min( jgrow, nimptblgrow_maxd-1 )
1498                 end if
1499
1500                 dumfhi = xgrow - jgrow
1501                 dumflo = 1. - dumfhi
1502
1503                 scavimpvol = dumflo*scavimptblvol(jgrow,m) +   &
1504                           dumfhi*scavimptblvol(jgrow+1,m)
1505                 scavimpnum = dumflo*scavimptblnum(jgrow,m) +   &
1506                           dumfhi*scavimptblnum(jgrow+1,m)
1507
1508             end if
1509
1510 !           impaction scavenging removal amount for volume
1511             scavcoefvol(i,k) = exp( scavimpvol  )
1512 !           impaction scavenging removal amount to number
1513             scavcoefnum(i,k) = exp( scavimpnum  )
1514
1515 !           scavcoef = impaction scav rate (1/h) for precip = 1 mm/h
1516 !           scavcoef = impaction scav rate (1/s) for precip = pfx_inrain
1517 !           (scavcoef/3600) = impaction scav rate (1/s) for precip = 1 mm/h
1518 !           (pfx_inrain*3600) = in-rain-area precip rate (mm/h)
1519 !           impactrate = (scavcoef/3600) * (pfx_inrain*3600)
1520           else
1521             scavcoefvol(i,k) = 0._r8
1522             scavcoefnum(i,k) = 0._r8
1523           end if
1524
1525       end do
1526       end do
1527
1528         return
1529         end subroutine modal_aero_bcscavcoef_get
1530
1531
1532   !===============================================================================
1533   subroutine modal_aero_bcscavcoef_init
1534 !-----------------------------------------------------------------------
1535 !
1536 ! Purpose:
1537 ! Computes lookup table for aerosol impaction/interception scavenging rates
1538 !
1539 ! Authors: R. Easter
1540 !
1541 !-----------------------------------------------------------------------
1542
1543   use shr_kind_mod,only: r8 => shr_kind_r8
1544   use modal_aero_data
1545 #ifndef WRF_PORT
1546   use abortutils,  only: endrun
1547 #else
1548   use module_cam_support, only: endrun
1549 #endif
1550
1551   implicit none
1552
1553
1554 !   local variables
1555         integer nnfit_maxd
1556         parameter (nnfit_maxd=27)
1557
1558         integer i, jgrow, jdens, jpress, jtemp, ll, mode, nnfit
1559         integer lunerr
1560
1561         real(r8) dg0, dg0_cgs, press, &
1562         rhodryaero, rhowetaero, rhowetaero_cgs, rmserr, &
1563         scavratenum, scavratevol, sigmag,                &
1564         temp, wetdiaratio, wetvolratio
1565         real(r8) aafitnum(1), xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd)
1566         real(r8) aafitvol(1), xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd)
1567
1568         lunerr = 6
1569         dlndg_nimptblgrow = log( 1.25d00 )
1570
1571         do 7900 mode = 1, ntot_amode
1572
1573         sigmag = sigmag_amode(mode)
1574
1575         ll = lspectype_amode(1,mode)
1576         rhodryaero = specdens_amode(ll)
1577
1578         do 7800 jgrow = nimptblgrow_mind, nimptblgrow_maxd
1579
1580         wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
1581         dg0 = dgnum_amode(mode)*wetdiaratio
1582
1583         wetvolratio = exp( jgrow*dlndg_nimptblgrow*3. )
1584         rhowetaero = 1.0 + (rhodryaero-1.0)/wetvolratio
1585         rhowetaero = min( rhowetaero, rhodryaero )
1586
1587 !
1588 !   compute impaction scavenging rates at 1 temp-press pair and save
1589 !
1590         nnfit = 0
1591
1592         temp = 273.16
1593         press = 0.75e6   ! dynes/cm2
1594         rhowetaero = rhodryaero
1595
1596         dg0_cgs = dg0*1.0e2   ! m to cm
1597         rhowetaero_cgs = rhowetaero*1.0e-3   ! kg/m3 to g/cm3
1598         call calc_1_impact_rate( &
1599                 dg0_cgs, sigmag, rhowetaero_cgs, temp, press, &
1600                 scavratenum, scavratevol, lunerr )
1601
1602         nnfit = nnfit + 1
1603         if (nnfit .gt. nnfit_maxd) then
1604             write(lunerr,9110)
1605             call endrun()
1606         end if
1607 9110    format( '*** subr. modal_aero_bcscavcoef_init -- nnfit too big' )
1608
1609         xxfitnum(1,nnfit) = 1.
1610         yyfitnum(nnfit) = log( scavratenum )
1611
1612         xxfitvol(1,nnfit) = 1.
1613         yyfitvol(nnfit) = log( scavratevol )
1614
1615 5900    continue
1616
1617 !
1618 ! skip mlinfit stuff because scav table no longer has dependencies on
1619 !    air temp, air press, and particle wet density
1620 ! just load the log( scavrate--- ) values
1621 !
1622 !!
1623 !!   do linear regression
1624 !!      log(scavrate) = a1 + a2*log(wetdens)
1625 !!
1626 !       call mlinft( xxfitnum, yyfitnum, aafitnum, nnfit, 1, 1, rmserr )
1627 !       call mlinft( xxfitvol, yyfitvol, aafitvol, nnfit, 1, 1, rmserr )
1628 !
1629 !       scavimptblnum(jgrow,mode) = aafitnum(1)
1630 !       scavimptblvol(jgrow,mode) = aafitvol(1)
1631
1632         scavimptblnum(jgrow,mode) = yyfitnum(1)
1633         scavimptblvol(jgrow,mode) = yyfitvol(1)
1634
1635
1636
1637 7800    continue
1638 7900    continue
1639
1640         return
1641         end subroutine modal_aero_bcscavcoef_init
1642
1643
1644   !===============================================================================
1645         subroutine calc_1_impact_rate(             &
1646                 dg0, sigmag, rhoaero, temp, press, &
1647                 scavratenum, scavratevol, lunerr )
1648 !
1649 !   routine computes a single impaction scavenging rate
1650 !       for precipitation rate of 1 mm/h
1651 !
1652 !   dg0 = geometric mean diameter of aerosol number size distrib. (cm)
1653 !   sigmag = geometric standard deviation of size distrib.
1654 !   rhoaero = density of aerosol particles (g/cm^3)
1655 !   temp = temperature (K)
1656 !   press = pressure (dyne/cm^2)
1657 !   scavratenum = number scavenging rate (1/h)
1658 !   scavratevol = volume or mass scavenging rate (1/h)
1659 !   lunerr = logical unit for error message
1660 !
1661      use shr_kind_mod, only: r8 => shr_kind_r8
1662
1663         implicit none
1664
1665 !   subr. parameters
1666         integer lunerr
1667         real(r8) dg0, sigmag, rhoaero, temp, press, scavratenum, scavratevol
1668
1669 !   local variables
1670         integer nrainsvmax
1671         parameter (nrainsvmax=50)
1672         real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
1673                 vfallrainsv(nrainsvmax)
1674
1675         integer naerosvmax
1676         parameter (naerosvmax=51)
1677         real(r8) aaerosv(naerosvmax), &
1678         ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)
1679
1680         integer i, ja, jr, na, nr
1681         real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
1682         real(r8) anumsum, avolsum, boltzmann, cair, chi
1683         real(r8) d, dr, dum, dumfuchs, dx
1684         real(r8) ebrown, eimpact, eintercept, etotal, freepath, gravity
1685         real(r8) pi, precip, precipmmhr, precipsum
1686         real(r8) r, rainsweepout, reynolds, rhi, rhoair, rhowater, rlo, rnumsum
1687         real(r8) scavsumnum, scavsumnumbb
1688         real(r8) scavsumvol, scavsumvolbb
1689         real(r8) schmidt, sqrtreynolds, sstar, stokes, sx
1690         real(r8) taurelax, vfall, vfallstp
1691         real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair
1692
1693
1694         rlo = .005
1695         rhi = .250
1696         dr = 0.005
1697         nr = 1 + nint( (rhi-rlo)/dr )
1698         if (nr .gt. nrainsvmax) then
1699             write(lunerr,9110)
1700             call endrun()
1701         end if
1702 9110    format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' )
1703
1704         precipmmhr = 1.0
1705         precip = precipmmhr/36000.
1706
1707         ag0 = dg0/2.
1708         sx = log( sigmag )
1709         xg0 = log( ag0 )
1710         xg3 = xg0 + 3.*sx*sx
1711
1712         xlo = xg3 - 4.*sx
1713         xhi = xg3 + 4.*sx
1714         dx = 0.2*sx
1715
1716         dx = max( 0.2_r8*sx, 0.01_r8 )
1717         xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
1718         xhi = xg3 + max( 4._r8*sx, 2._r8*dx )
1719
1720         na = 1 + nint( (xhi-xlo)/dx )
1721         if (na .gt. naerosvmax) then
1722             write(lunerr,9120)
1723             call endrun()
1724         end if
1725 9120    format( '*** subr. calc_1_impact_rate -- na > naerosvmax' )
1726
1727 !   air molar density
1728         cair = press/(8.31436e7*temp)
1729 !   air mass density
1730         rhoair = 28.966*cair
1731 !   molecular freepath
1732         freepath = 2.8052e-10/cair
1733 !   boltzmann constant
1734         boltzmann = 1.3807e-16
1735 !   water density
1736         rhowater = 1.0
1737 !   gravity
1738         gravity = 980.616
1739 !   air dynamic viscosity
1740         airdynvisc = 1.8325e-4 * (416.16/(temp+120.)) *    &
1741                                         ((temp/296.16)**1.5)
1742 !   air kinemaic viscosity
1743         airkinvisc = airdynvisc/rhoair
1744 !   ratio of water viscosity to air viscosity (from Slinn)
1745         xmuwaterair = 60.0
1746
1747         pi = 3.1415926536
1748
1749 !
1750 !   compute rain drop number concentrations
1751 !       rrainsv = raindrop radius (cm)
1752 !       xnumrainsv = raindrop number concentration (#/cm^3)
1753 !               (number in the bin, not number density)
1754 !       vfallrainsv = fall velocity (cm/s)
1755 !
1756         precipsum = 0.
1757         do i = 1, nr
1758             r = rlo + (i-1)*dr
1759             rrainsv(i) = r
1760             xnumrainsv(i) = exp( -r/2.7e-2 )
1761
1762             d = 2.*r
1763             if (d .le. 0.007) then
1764                 vfallstp = 2.88e5 * d**2.
1765             else if (d .le. 0.025) then
1766                 vfallstp = 2.8008e4 * d**1.528
1767             else if (d .le. 0.1) then
1768                 vfallstp = 4104.9 * d**1.008
1769             else if (d .le. 0.25) then
1770                 vfallstp = 1812.1 * d**0.638
1771             else
1772                 vfallstp = 1069.8 * d**0.235
1773             end if
1774
1775             vfall = vfallstp * sqrt(1.204e-3/rhoair)
1776             vfallrainsv(i) = vfall
1777             precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
1778         end do
1779         precipsum = precipsum*pi*1.333333
1780
1781         rnumsum = 0.
1782         do i = 1, nr
1783             xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
1784             rnumsum = rnumsum + xnumrainsv(i)
1785         end do
1786
1787 !
1788 !   compute aerosol concentrations
1789 !       aaerosv = particle radius (cm)
1790 !       fnumaerosv = fraction of total number in the bin (--)
1791 !       fvolaerosv = fraction of total volume in the bin (--)
1792 !
1793         anumsum = 0.
1794         avolsum = 0.
1795         do i = 1, na
1796             x = xlo + (i-1)*dx
1797             a = exp( x )
1798             aaerosv(i) = a
1799             dum = (x - xg0)/sx
1800             ynumaerosv(i) = exp( -0.5*dum*dum )
1801             yvolaerosv(i) = ynumaerosv(i)*1.3333*pi*a*a*a
1802             anumsum = anumsum + ynumaerosv(i)
1803             avolsum = avolsum + yvolaerosv(i)
1804         end do
1805
1806         do i = 1, na
1807             ynumaerosv(i) = ynumaerosv(i)/anumsum
1808             yvolaerosv(i) = yvolaerosv(i)/avolsum
1809         end do
1810
1811
1812 !
1813 !   compute scavenging
1814 !
1815         scavsumnum = 0.
1816         scavsumvol = 0.
1817 !
1818 !   outer loop for rain drop radius
1819 !
1820         do 5900 jr = 1, nr
1821
1822         r = rrainsv(jr)
1823         vfall = vfallrainsv(jr)
1824
1825         reynolds = r * vfall / airkinvisc
1826         sqrtreynolds = sqrt( reynolds )
1827
1828 !
1829 !   inner loop for aerosol particle radius
1830 !
1831         scavsumnumbb = 0.
1832         scavsumvolbb = 0.
1833
1834         do 5500 ja = 1, na
1835
1836         a = aaerosv(ja)
1837
1838         chi = a/r
1839
1840         dum = freepath/a
1841         dumfuchs = 1. + 1.246*dum + 0.42*dum*exp(-0.87/dum)
1842         taurelax = 2.*rhoaero*a*a*dumfuchs/(9.*rhoair*airkinvisc)
1843
1844         aeromass = 4.*pi*a*a*a*rhoaero/3.
1845         aerodiffus = boltzmann*temp*taurelax/aeromass
1846
1847         schmidt = airkinvisc/aerodiffus
1848         stokes = vfall*taurelax/r
1849
1850         ebrown = 4.*(1. + 0.4*sqrtreynolds*(schmidt**0.3333333)) /  &
1851                         (reynolds*schmidt)
1852
1853         dum = (1. + 2.*xmuwaterair*chi) /         &
1854                         (1. + xmuwaterair/sqrtreynolds)
1855         eintercept = 4.*chi*(chi + dum)
1856
1857         dum = log( 1. + reynolds )
1858         sstar = (1.2 + dum/12.) / (1. + dum)
1859         eimpact = 0.
1860         if (stokes .gt. sstar) then
1861             dum = stokes - sstar
1862             eimpact = (dum/(dum+0.6666667)) ** 1.5
1863         end if
1864
1865         etotal = ebrown + eintercept + eimpact
1866         etotal = min( etotal, 1.0_r8 )
1867
1868         rainsweepout = xnumrainsv(jr)*4.*pi*r*r*vfall
1869
1870         scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
1871         scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
1872
1873 5500    continue
1874
1875         scavsumnum = scavsumnum + scavsumnumbb
1876         scavsumvol = scavsumvol + scavsumvolbb
1877 5900    continue
1878
1879         scavratenum = scavsumnum*3600.
1880         scavratevol = scavsumvol*3600.
1881
1882
1883
1884   return
1885   end subroutine calc_1_impact_rate
1886
1887 #endif
1888 ! (defined MODAL_AERO)
1889
1890 #ifndef WRF_PORT
1891 !-------------------------------------------------------------------
1892 !-------------------------------------------------------------------
1893   !===============================================================================
1894
1895   function mz_prescribed_dust()
1896
1897      ! Return true if MOZART is supplying prescribed dust.
1898
1899      use dust_intr,    only: dust_names
1900      use tracer_cnst,  only: num_tracer_cnst, tracer_cnst_flds
1901
1902      logical :: mz_prescribed_dust
1903
1904      integer :: i
1905      !-----------------------------------------------------------------------
1906
1907      mz_prescribed_dust = .false.
1908      do i = 1, num_tracer_cnst
1909         if ( trim(tracer_cnst_flds(i)) == dust_names(1) ) then
1910            mz_prescribed_dust = .true.
1911            exit
1912         end if
1913      end do
1914
1915   end function mz_prescribed_dust
1916
1917   !===============================================================================
1918 #endif
1919
1920 end module mz_aerosols_intr