LCOV - code coverage report
Current view: | top level - chemistry/modal_aero - aero_model.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 580 | 708 | 81.9 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 12 | 18 | 66.7 %

              Line data    Source code
           1             : !===============================================================================
           2             : ! Modal Aerosol Model
           3             : !===============================================================================
           4             : module aero_model
           5             :   use shr_kind_mod,   only: r8 => shr_kind_r8
           6             :   use constituents,   only: pcnst, cnst_name, cnst_get_ind
           7             :   use ppgrid,         only: pcols, pver, pverp
           8             :   use phys_control,   only: phys_getopts, cam_physpkg_is
           9             :   use cam_abortutils, only: endrun
          10             :   use cam_logfile,    only: iulog
          11             :   use perf_mod,       only: t_startf, t_stopf
          12             :   use camsrfexch,     only: cam_in_t, cam_out_t
          13             :   use aerodep_flx,    only: aerodep_flx_prescribed
          14             :   use physics_types,  only: physics_state, physics_ptend, physics_ptend_init
          15             :   use physics_buffer, only: physics_buffer_desc
          16             :   use physics_buffer, only: pbuf_get_field, pbuf_get_index, pbuf_set_field
          17             :   use physconst,      only: gravit, rair, rhoh2o
          18             :   use spmd_utils,     only: masterproc
          19             :   use infnan,         only: nan, assignment(=)
          20             :
          21             :   use cam_history,    only: outfld, fieldname_len
          22             :   use chem_mods,      only: gas_pcnst, adv_mass
          23             :   use mo_tracname,    only: solsym
          24             :
          25             :   use modal_aero_data,only: cnst_name_cw, lptr_so4_cw_amode
          26             :   use modal_aero_data,only: ntot_amode, modename_amode, nspec_max
          27             :
          28             :   use ref_pres,       only: top_lev => clim_modal_aero_top_lev
          29             :
          30             :   use modal_aero_wateruptake, only: modal_strat_sulfate
          31             :   use mo_setsox,              only: setsox, has_sox
          32             :   use modal_aerosol_properties_mod, only: modal_aerosol_properties
          33             :
          34             :   implicit none
          35             :   private
          36             :
          37             :   public :: aero_model_readnl
          38             :   public :: aero_model_register
          39             :   public :: aero_model_init
          40             :   public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols.
          41             :   public :: aero_model_drydep     ! aerosol dry deposition and sediment
          42             :   public :: aero_model_wetdep     ! aerosol wet removal
          43             :   public :: aero_model_emissions  ! aerosol emissions
          44             :   public :: aero_model_surfarea  ! tropopspheric aerosol wet surface area for chemistry
          45             :   public :: aero_model_strat_surfarea ! stratospheric aerosol wet surface area for chemistry
          46             :
          47             :   public :: calc_1_impact_rate
          48             :   public :: nimptblgrow_mind, nimptblgrow_maxd
          49             :
          50             :   ! Accessor functions
          51             :   public ::  get_scavimptblvol, get_scavimptblnum, get_dlndg_nimptblgrow
          52             :
          53             :   ! Misc private data
          54             :
          55             :   ! number of modes
          56             :   integer :: nmodes
          57             :   integer :: pblh_idx            = 0
          58             :   integer :: dgnum_idx           = 0
          59             :   integer :: dgnumwet_idx        = 0
          60             :   integer :: rate1_cw2pr_st_idx  = 0
          61             :
          62             :   integer :: wetdens_ap_idx      = 0
          63             :   integer :: qaerwat_idx         = 0
          64             :
          65             :   integer :: fracis_idx          = 0
          66             :   integer :: prain_idx           = 0
          67             :   integer :: rprddp_idx          = 0
          68             :   integer :: rprdsh_idx          = 0
          69             :   integer :: nevapr_shcu_idx     = 0
          70             :   integer :: nevapr_dpcu_idx     = 0
          71             :
          72             :   integer :: sulfeq_idx = -1
          73             :
          74             :   integer :: nh3_ndx    = 0
          75             :   integer :: nh4_ndx    = 0
          76             :
          77             :   ! variables for table lookup of aerosol impaction/interception scavenging rates
          78             :   integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
          79             :   real(r8) :: dlndg_nimptblgrow
          80             :   real(r8),allocatable :: scavimptblnum(:,:)
          81             :   real(r8),allocatable :: scavimptblvol(:,:)
          82             :
          83             :   ! for surf_area_dens
          84             :   integer,allocatable :: num_idx(:)
          85             :   integer,allocatable :: index_tot_mass(:,:)
          86             :   integer,allocatable :: index_chm_mass(:,:)
          87             :
          88             :   integer :: ndx_h2so4
          89             :   character(len=fieldname_len), allocatable :: dgnum_name(:), dgnumwet_name(:)
          90             :
          91             :   ! Namelist variables
          92             :   character(len=16) :: drydep_list(pcnst) = ' '
          93             :   real(r8)          :: seasalt_emis_scale
          94             :
          95             :   integer :: ndrydep = 0
          96             :   integer,allocatable :: drydep_indices(:)
          97             :   logical :: drydep_lq(pcnst)
          98             :
          99             :   logical :: modal_accum_coarse_exch = .false.
         100             :
         101             :   type(modal_aerosol_properties), pointer :: aero_props=>null()
         102             :
         103             : contains
         104             :
         105             :   !=============================================================================
         106             :   ! reads aerosol namelist options
         107             :   !=============================================================================
         108        1536 :   subroutine aero_model_readnl(nlfile)
         109             :
         110             :     use namelist_utils,  only: find_group_name
         111             :     use units,           only: getunit, freeunit
         112             :     use mpishorthand
         113             :     use aero_wetdep_cam, only: aero_wetdep_readnl
         114             :
         115             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         116             :
         117             :     ! Local variables
         118             :     integer :: unitn, ierr
         119             :     character(len=*), parameter :: subname = 'aero_model_readnl'
         120             :
         121             :     ! Namelist variables
         122             :     character(len=16) :: aer_drydep_list(pcnst) = ' '
         123             :
         124             :     namelist /aerosol_nl/ aer_drydep_list, modal_strat_sulfate, modal_accum_coarse_exch, seasalt_emis_scale
         125             :
         126             :     !-----------------------------------------------------------------------------
         127             :
         128             :     ! Read namelist
         129        1536 :     if (masterproc) then
         130           2 :        unitn = getunit()
         131           2 :        open( unitn, file=trim(nlfile), status='old' )
         132           2 :        call find_group_name(unitn, 'aerosol_nl', status=ierr)
         133           2 :        if (ierr == 0) then
         134           2 :           read(unitn, aerosol_nl, iostat=ierr)
         135           2 :           if (ierr /= 0) then
         136           0 :              call endrun(subname // ':: ERROR reading namelist')
         137             :           end if
         138             :        end if
         139           2 :        close(unitn)
         140           2 :        call freeunit(unitn)
         141             :     end if
         142             :
         143             : #ifdef SPMD
         144             :     ! Broadcast namelist variables
         145        1536 :     call mpibcast(aer_drydep_list,   len(aer_drydep_list(1))*pcnst, mpichar, 0, mpicom)
         146        1536 :     call mpibcast(modal_strat_sulfate,     1,                       mpilog,  0, mpicom)
         147        1536 :     call mpibcast(seasalt_emis_scale, 1,                            mpir8,   0, mpicom)
         148        1536 :     call mpibcast(modal_accum_coarse_exch, 1,                       mpilog,  0, mpicom)
         149             : #endif
         150             :
         151       64512 :     drydep_list = aer_drydep_list
         152             :
         153        1536 :     call aero_wetdep_readnl(nlfile)
         154             :
         155        1536 :   end subroutine aero_model_readnl
         156             :
         157             :   !=============================================================================
         158             :   !=============================================================================
         159        1536 :   subroutine aero_model_register()
         160        1536 :     use modal_aero_data, only: modal_aero_data_reg
         161             :
         162        1536 :     call modal_aero_data_reg()
         163             :
         164        1536 :   end subroutine aero_model_register
         165             :
         166             :   !=============================================================================
         167             :   !=============================================================================
         168        1536 :   subroutine aero_model_init( pbuf2d )
         169             :
         170             :     use mo_chem_utls,    only: get_inv_ndx
         171             :     use cam_history,     only: addfld, add_default, horiz_only
         172        1536 :     use mo_chem_utls,    only: get_spc_ndx
         173             :     use modal_aero_data, only: cnst_name_cw
         174             :     use modal_aero_data, only: modal_aero_data_init
         175             :     use rad_constituents,only: rad_cnst_get_info
         176             :     use dust_model,      only: dust_init, dust_names, dust_active, dust_nbin, dust_nnum
         177             :     use seasalt_model,   only: seasalt_init, seasalt_names, seasalt_active,seasalt_nbin
         178             :     use aer_drydep_mod,  only: inidrydep
         179             :     use aero_wetdep_cam, only: aero_wetdep_init
         180             :
         181             :     use modal_aero_calcsize,   only: modal_aero_calcsize_init
         182             :     use modal_aero_coag,       only: modal_aero_coag_init
         183             :     use aero_deposition_cam, only: aero_deposition_cam_init
         184             :     use modal_aero_gasaerexch, only: modal_aero_gasaerexch_init
         185             :     use modal_aero_newnuc,     only: modal_aero_newnuc_init
         186             :     use modal_aero_rename,     only: modal_aero_rename_init
         187             :
         188             :     ! args
         189             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         190             :
         191             :     ! local vars
         192             :     character(len=*), parameter :: subrname = 'aero_model_init'
         193             :     integer :: m, n, id
         194             :     character(len=20) :: dummy
         195             :
         196             :     logical  :: history_aerosol ! Output MAM or SECT aerosol tendencies
         197             :     logical  :: history_chemistry, history_cesm_forcing, history_dust
         198             :
         199             :     integer :: l
         200             :     character(len=6) :: test_name
         201             :     character(len=64) :: errmes
         202             :
         203             :     character(len=2)  :: unit_basename  ! Units 'kg' or '1'
         204             :     integer :: errcode
         205             :     character(len=fieldname_len) :: field_name
         206             :
         207             :     character(len=32) :: spec_name
         208             :     character(len=32) :: spec_type
         209             :     character(len=32) :: mode_type
         210             :     integer :: nspec
         211             :
         212        1536 :     dgnum_idx       = pbuf_get_index('DGNUM')
         213        1536 :     dgnumwet_idx    = pbuf_get_index('DGNUMWET')
         214        1536 :     fracis_idx      = pbuf_get_index('FRACIS')
         215        1536 :     prain_idx       = pbuf_get_index('PRAIN')
         216        1536 :     rprddp_idx      = pbuf_get_index('RPRDDP')
         217        1536 :     rprdsh_idx      = pbuf_get_index('RPRDSH')
         218        1536 :     nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
         219        1536 :     nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
         220        1536 :     sulfeq_idx      = pbuf_get_index('MAMH2SO4EQ',errcode)
         221             :
         222             :     call phys_getopts(history_aerosol_out = history_aerosol, &
         223             :                       history_chemistry_out=history_chemistry, &
         224             :                       history_cesm_forcing_out=history_cesm_forcing, &
         225        1536 :                       history_dust_out=history_dust)
         226             :
         227        1536 :     call rad_cnst_get_info(0, nmodes=nmodes)
         228             :
         229        1536 :     call modal_aero_data_init(pbuf2d)
         230        1536 :     call modal_aero_bcscavcoef_init()
         231             :
         232        1536 :     call modal_aero_rename_init( modal_accum_coarse_exch )
         233             :     !   calcsize call must follow rename call
         234        1536 :     call modal_aero_calcsize_init( pbuf2d )
         235        1536 :     call modal_aero_gasaerexch_init
         236             :     !   coag call must follow gasaerexch call
         237        1536 :     call modal_aero_coag_init
         238        1536 :     call modal_aero_newnuc_init
         239             :
         240             :     ! call aero_deposition_cam_init only if the user has not specified
         241             :     ! prescribed aerosol deposition fluxes
         242        1536 :     if (.not.aerodep_flx_prescribed()) then
         243        1536 :        aero_props => modal_aerosol_properties()
         244        1536 :        call aero_deposition_cam_init(aero_props)
         245             :     endif
         246             :
         247        1536 :     call dust_init()
         248        1536 :     call seasalt_init(seasalt_emis_scale)
         249             :
         250        1536 :     ndrydep = 0
         251             :
         252       64512 :     count_species: do m = 1,pcnst
         253       64512 :        if ( len_trim(drydep_list(m)) /= 0 ) then
         254       29184 :           ndrydep = ndrydep+1
         255             :        endif
         256             :     enddo count_species
         257             :
         258        1536 :     if (ndrydep>0) &
         259        4608 :          allocate(drydep_indices(ndrydep))
         260             :
         261       30720 :     do m = 1,ndrydep
         262       29184 :        call cnst_get_ind ( drydep_list(m), id, abort=.false. )
         263       29184 :        if (id>0) then
         264       29184 :           drydep_indices(m) = id
         265             :        else
         266           0 :           call endrun(subrname//': invalid drydep species: '//trim(drydep_list(m)) )
         267             :        endif
         268             :
         269       30720 :        if (masterproc) then
         270          38 :           write(iulog,*) subrname//': '//drydep_list(m)//' will have drydep applied'
         271             :        endif
         272             :     enddo
         273             :
         274        1536 :     if (ndrydep>0) then
         275             :
         276        1536 :        call inidrydep(rair, gravit)
         277             :
         278        1536 :        dummy = 'RAM1'
         279        1536 :        call addfld (dummy,horiz_only, 'A','frac','RAM1')
         280        1536 :        if ( history_aerosol ) then
         281           0 :           call add_default (dummy, 1, ' ')
         282             :        endif
         283        1536 :        dummy = 'airFV'
         284        1536 :        call addfld (dummy,horiz_only, 'A','frac','FV')
         285        1536 :        if ( history_aerosol ) then
         286           0 :           call add_default (dummy, 1, ' ')
         287             :        endif
         288             :
         289             :     endif
         290             :
         291        1536 :     if (dust_active) then
         292             :        ! emissions diagnostics ....
         293             :
         294       10752 :        do m = 1, dust_nbin+dust_nnum
         295        9216 :           dummy = trim(dust_names(m)) // 'SF'
         296        9216 :           call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(dust_names(m))//' dust surface emission')
         297       10752 :           if (history_aerosol.or.history_chemistry) then
         298        9216 :              call add_default (dummy, 1, ' ')
         299             :           endif
         300             :        enddo
         301             :
         302        1536 :        dummy = 'DSTSFMBL'
         303        1536 :        call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
         304        1536 :        if (history_aerosol .or. history_dust) then
         305           0 :           call add_default (dummy, 1, ' ')
         306             :        endif
         307             :
         308        1536 :        dummy = 'LND_MBL'
         309        1536 :        call addfld (dummy,horiz_only, 'A','frac','Soil erodibility factor')
         310        1536 :        if (history_aerosol) then
         311           0 :           call add_default (dummy, 1, ' ')
         312             :        endif
         313             :
         314             :     endif
         315             :
         316        1536 :     if (seasalt_active) then
         317             :
         318        1536 :        dummy = 'SSTSFMBL'
         319        1536 :        call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
         320        1536 :        if (history_aerosol) then
         321           0 :           call add_default (dummy, 1, ' ')
         322             :        endif
         323             :
         324        6144 :        do m = 1, seasalt_nbin
         325        4608 :           dummy = trim(seasalt_names(m)) // 'SF'
         326        4608 :           call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(seasalt_names(m))//' seasalt surface emission')
         327        6144 :           if (history_aerosol.or.history_chemistry) then
         328        4608 :              call add_default (dummy, 1, ' ')
         329             :           endif
         330             :        enddo
         331             :
         332             :     endif
         333             :
         334             :
         335             :     ! set flags for drydep tendencies
         336        1536 :     drydep_lq(:) = .false.
         337       30720 :     do m=1,ndrydep
         338       29184 :        id = drydep_indices(m)
         339       30720 :        drydep_lq(id) =  .true.
         340             :     enddo
         341             :
         342        1536 :     wetdens_ap_idx = pbuf_get_index('WETDENS_AP')
         343        1536 :     qaerwat_idx    = pbuf_get_index('QAERWAT')
         344        1536 :     pblh_idx       = pbuf_get_index('pblh')
         345             :
         346        1536 :     rate1_cw2pr_st_idx  = pbuf_get_index('RATE1_CW2PR_ST')
         347        1536 :     call pbuf_set_field(pbuf2d, rate1_cw2pr_st_idx, 0.0_r8)
         348             :
         349       30720 :     do m = 1,ndrydep
         350             :
         351             :        ! units
         352       29184 :        if (drydep_list(m)(1:3) == 'num') then
         353        6144 :           unit_basename = ' 1'
         354             :        else
         355       23040 :           unit_basename = 'kg'
         356             :        endif
         357             :
         358             :        call addfld (trim(drydep_list(m))//'DDF', horiz_only,  'A',unit_basename//'/m2/s ', &
         359       29184 :             trim(drydep_list(m))//' dry deposition flux at bottom (grav + turb)')
         360             :        call addfld (trim(drydep_list(m))//'TBF', horiz_only,  'A',unit_basename//'/m2/s',  &
         361       29184 :             trim(drydep_list(m))//' turbulent dry deposition flux')
         362             :        call addfld (trim(drydep_list(m))//'GVF', horiz_only,  'A',unit_basename//'/m2/s ', &
         363       29184 :             trim(drydep_list(m))//' gravitational dry deposition flux')
         364             :        call addfld (trim(drydep_list(m))//'DTQ', (/ 'lev' /), 'A',unit_basename//'/kg/s ', &
         365       58368 :             trim(drydep_list(m))//' dry deposition')
         366             :        call addfld (trim(drydep_list(m))//'DDV', (/ 'lev' /), 'A','m/s',                   &
         367       58368 :             trim(drydep_list(m))//' deposition velocity')
         368             :
         369       29184 :        if ( history_aerosol.or.history_chemistry ) then
         370       29184 :           call add_default (trim(drydep_list(m))//'DDF', 1, ' ')
         371             :        endif
         372       30720 :        if ( history_aerosol ) then
         373           0 :           call add_default (trim(drydep_list(m))//'TBF', 1, ' ')
         374           0 :           call add_default (trim(drydep_list(m))//'GVF', 1, ' ')
         375             :        endif
         376             :
         377             :     enddo
         378             :
         379       49152 :     do m = 1,gas_pcnst
         380             :
         381       47616 :        if  ( solsym(m)(1:3) == 'num') then
         382        6144 :           unit_basename = ' 1'  ! Units 'kg' or '1'
         383             :        else
         384       41472 :           unit_basename = 'kg'  ! Units 'kg' or '1'
         385             :        end if
         386             :
         387             :        call addfld( 'GS_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
         388       47616 :                     trim(solsym(m))//' gas chemistry/wet removal (for gas species)')
         389             :        call addfld( 'AQ_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
         390       47616 :                     trim(solsym(m))//' aqueous chemistry (for gas species)')
         391       49152 :        if ( history_aerosol ) then
         392           0 :           call add_default( 'AQ_'//trim(solsym(m)), 1, ' ')
         393             :        endif
         394             :
         395             :     enddo
         396       64512 :     do n = 1,pcnst
         397       64512 :        if( .not. (cnst_name_cw(n) == ' ') ) then
         398             :
         399       29184 :           if (cnst_name_cw(n)(1:3) == 'num') then
         400        6144 :              unit_basename = ' 1'
         401             :           else
         402       23040 :              unit_basename = 'kg'
         403             :           endif
         404             :
         405             :           call addfld( cnst_name_cw(n),                (/ 'lev' /), 'A', unit_basename//'/kg ',   &
         406       58368 :                trim(cnst_name_cw(n))//' in cloud water')
         407       29184 :           call addfld (trim(cnst_name_cw(n))//'DDF',   horiz_only,  'A', unit_basename//'/m2/s ', &
         408       58368 :                trim(cnst_name_cw(n))//' dry deposition flux at bottom (grav + turb)')
         409       29184 :           call addfld (trim(cnst_name_cw(n))//'TBF',   horiz_only,  'A', unit_basename//'/m2/s ', &
         410       58368 :                trim(cnst_name_cw(n))//' turbulent dry deposition flux')
         411       29184 :           call addfld (trim(cnst_name_cw(n))//'GVF',   horiz_only,  'A', unit_basename//'/m2/s ', &
         412       58368 :                trim(cnst_name_cw(n))//' gravitational dry deposition flux')
         413             :
         414       29184 :           if ( history_aerosol.or. history_chemistry ) then
         415       29184 :              call add_default( cnst_name_cw(n), 1, ' ' )
         416             :           endif
         417       29184 :           if ( history_aerosol ) then
         418           0 :              call add_default (trim(cnst_name_cw(n))//'GVF', 1, ' ')
         419           0 :              call add_default (trim(cnst_name_cw(n))//'TBF', 1, ' ')
         420           0 :              call add_default (trim(cnst_name_cw(n))//'DDF', 1, ' ')
         421             :           endif
         422             :        endif
         423             :     enddo
         424             :
         425        6144 :     allocate(dgnum_name(ntot_amode), dgnumwet_name(ntot_amode))
         426        7680 :     do n=1,ntot_amode
         427        6144 :        dgnum_name(n) = ' '
         428        6144 :        dgnumwet_name(n) = ' '
         429        6144 :        write(dgnum_name(n),fmt='(a,i1)') 'dgnum',n
         430        6144 :        write(dgnumwet_name(n),fmt='(a,i1)') 'dgnumwet',n
         431       12288 :        call addfld( dgnum_name(n),    (/ 'lev' /), 'I', 'm', 'Aerosol mode dry diameter' )
         432       12288 :        call addfld( dgnumwet_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode wet diameter' )
         433        6144 :        if ( history_aerosol ) then
         434           0 :           call add_default( dgnum_name(n), 1, ' ' )
         435           0 :           call add_default( dgnumwet_name(n), 1, ' ' )
         436             :        endif
         437        6144 :        if ( history_cesm_forcing .and. n<4 ) then
         438           0 :           call add_default( dgnumwet_name(n), 8, ' ' )
         439             :        endif
         440             :
         441        7680 :        if (modal_strat_sulfate) then
         442           0 :           field_name = ' '
         443           0 :           write(field_name,fmt='(a,i1)') 'wtpct_a',n
         444           0 :           call addfld( field_name, (/ 'lev' /), 'I', '%', 'Aerosol mode weight percent H2SO4' )
         445           0 :           if ( history_aerosol ) then
         446           0 :              call add_default (field_name, 0, 'I')
         447             :           endif
         448             :
         449           0 :           field_name = ' '
         450           0 :           write(field_name,fmt='(a,i1)') 'sulfeq_a',n
         451           0 :           call addfld( field_name, (/ 'lev' /), 'I', 'kg/kg', 'H2SO4 equilibrium mixing ratio' )
         452           0 :           if ( history_aerosol ) then
         453           0 :              call add_default (field_name, 0, 'I')
         454             :           endif
         455             :
         456           0 :           field_name = ' '
         457           0 :           write(field_name,fmt='(a,i1)') 'sulden_a',n
         458           0 :           call addfld( field_name, (/ 'lev' /), 'I', 'g/cm3', 'Sulfate aerosol particle mass density' )
         459           0 :           if ( history_aerosol ) then
         460           0 :              call add_default (field_name, 0, 'I')
         461             :           endif
         462             :
         463             :        end if
         464             :     end do
         465             :
         466        1536 :     ndx_h2so4 = get_spc_ndx('H2SO4')
         467        1536 :     nh3_ndx = get_spc_ndx('NH3')
         468        1536 :     nh4_ndx = get_spc_ndx('NH4')
         469             :
         470        4608 :     allocate(num_idx(ntot_amode))
         471        7680 :     num_idx = -1
         472             :
         473             :     ! for aero_model_surfarea called from mo_usrrxt
         474        7680 :     do l=1,ntot_amode
         475        6144 :        test_name = ' '
         476        6144 :        write(test_name,fmt='(a5,i1)') 'num_a',l
         477        6144 :        num_idx(l) = get_spc_ndx( trim(test_name) )
         478        7680 :        if (num_idx(l) < 0) then
         479           0 :           write(errmes,fmt='(a,i1)') 'usrrxt_inti: cannot find MAM num_idx ',l
         480           0 :           write(iulog,*) errmes
         481           0 :           call endrun(errmes)
         482             :        endif
         483             :     end do
         484             :
         485        6144 :     allocate(index_tot_mass(nmodes,nspec_max))
         486        4608 :     allocate(index_chm_mass(nmodes,nspec_max))
         487       47616 :     index_tot_mass = -1
         488       47616 :     index_chm_mass = -1
         489             :
         490             :     ! for surf_area_dens
         491             :     ! define indices associated with the various aerosol types
         492        7680 :     do n = 1,nmodes
         493        6144 :        call rad_cnst_get_info(0, n, mode_type=mode_type, nspec=nspec)
         494        7680 :        if ( trim(mode_type) /= 'primary_carbon') then ! ignore the primary_carbon mode
         495       24576 :           do l = 1, nspec
         496       19968 :              call rad_cnst_get_info(0, n, l, spec_type=spec_type, spec_name=spec_name)
         497       19968 :              index_tot_mass(n,l) = get_spc_ndx(spec_name)
         498             :              if ( trim(spec_type) == 'sulfate'   .or. &
         499             :                   trim(spec_type) == 's-organic' .or. &
         500             :                   trim(spec_type) == 'p-organic' .or. &
         501       19968 :                   trim(spec_type) == 'black-c'   .or. &
         502        4608 :                   trim(spec_type) == 'ammonium') then
         503       10752 :                 index_chm_mass(n,l) = get_spc_ndx(spec_name)
         504             :              endif
         505             :           enddo
         506             :        endif
         507             :     enddo
         508             :
         509        1536 :     if (has_sox) then
         510        7680 :        do m = 1, ntot_amode
         511             :
         512        6144 :           l = lptr_so4_cw_amode(m)
         513        7680 :           if (l > 0) then
         514             :              call addfld (&
         515        4608 :                   trim(cnst_name_cw(l))//'AQSO4',horiz_only,  'A','kg/m2/s', &
         516        9216 :                   trim(cnst_name_cw(l))//' aqueous phase chemistry')
         517             :              call addfld (&
         518        4608 :                   trim(cnst_name_cw(l))//'AQH2SO4',horiz_only,  'A','kg/m2/s', &
         519        9216 :                   trim(cnst_name_cw(l))//' aqueous phase chemistry')
         520        4608 :              if ( history_aerosol ) then
         521           0 :                 call add_default (trim(cnst_name_cw(l))//'AQSO4', 1, ' ')
         522           0 :                 call add_default (trim(cnst_name_cw(l))//'AQH2SO4', 1, ' ')
         523             :              endif
         524             :           end if
         525             :
         526             :        end do
         527             :
         528        3072 :        call addfld( 'XPH_LWC',    (/ 'lev' /), 'A','kg/kg',   'pH value multiplied by lwc')
         529        1536 :        call addfld ('AQSO4_H2O2', horiz_only,  'A','kg/m2/s', 'SO4 aqueous phase chemistry due to H2O2')
         530        1536 :        call addfld ('AQSO4_O3',   horiz_only,  'A','kg/m2/s', 'SO4 aqueous phase chemistry due to O3')
         531             :
         532        1536 :        if ( history_aerosol ) then
         533           0 :           call add_default ('XPH_LWC', 1, ' ')
         534           0 :           call add_default ('AQSO4_H2O2', 1, ' ')
         535           0 :           call add_default ('AQSO4_O3', 1, ' ')
         536             :        endif
         537             :     endif
         538             :
         539        1536 :     call aero_wetdep_init()
         540             :
         541        3072 :   end subroutine aero_model_init
         542             :
         543             :   !=============================================================================
         544             :   !=============================================================================
         545    64034568 :   subroutine aero_model_drydep  ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend )
         546             :
         547        1536 :     use dust_sediment_mod, only: dust_sediment_tend
         548             :     use aer_drydep_mod,    only: d3ddflux, calcram
         549             :     use modal_aero_data,   only: qqcw_get_field
         550             :     use modal_aero_data,   only: cnst_name_cw
         551             :     use modal_aero_data,   only: alnsg_amode
         552             :     use modal_aero_data,   only: sigmag_amode
         553             :     use modal_aero_data,   only: nspec_amode
         554             :     use modal_aero_data,   only: numptr_amode
         555             :     use modal_aero_data,   only: numptrcw_amode
         556             :     use modal_aero_data,   only: lmassptr_amode
         557             :     use modal_aero_data,   only: lmassptrcw_amode
         558             :     use aero_deposition_cam,only: aero_deposition_cam_setdry
         559             :
         560             :   ! args
         561             :     type(physics_state),    intent(in)    :: state     ! Physics state variables
         562             :     real(r8),               intent(in)    :: obklen(:)
         563             :     real(r8),               intent(in)    :: ustar(:)  ! sfc fric vel
         564             :     type(cam_in_t), target, intent(in)    :: cam_in    ! import state
         565             :     real(r8),               intent(in)    :: dt             ! time step
         566             :     type(cam_out_t),        intent(inout) :: cam_out   ! export state
         567             :     type(physics_ptend),    intent(out)   :: ptend     ! indivdual parameterization tendencies
         568             :     type(physics_buffer_desc),    pointer :: pbuf(:)
         569             :
         570             :   ! local vars
         571     1489176 :     real(r8), pointer :: landfrac(:) ! land fraction
         572     1489176 :     real(r8), pointer :: icefrac(:)  ! ice fraction
         573     1489176 :     real(r8), pointer :: ocnfrac(:)  ! ocean fraction
         574     1489176 :     real(r8), pointer :: fvin(:)     !
         575     1489176 :     real(r8), pointer :: ram1in(:)   ! for dry dep velocities from land model for progseasalts
         576             :
         577             :     real(r8) :: fv(pcols)            ! for dry dep velocities, from land modified over ocean & ice
         578             :     real(r8) :: ram1(pcols)          ! for dry dep velocities, from land modified over ocean & ice
         579             :
         580             :     integer :: lchnk                   ! chunk identifier
         581             :     integer :: ncol                    ! number of atmospheric columns
         582             :     integer :: jvlc                    ! index for last dimension of vlc_xxx arrays
         583             :     integer :: lphase                  ! index for interstitial / cloudborne aerosol
         584             :     integer :: lspec                   ! index for aerosol number / chem-mass / water-mass
         585             :     integer :: m                       ! aerosol mode index
         586             :     integer :: mm                      ! tracer index
         587             :     integer :: i
         588             :
         589             :     real(r8) :: tvs(pcols,pver)
         590             :     real(r8) :: rho(pcols,pver)      ! air density in kg/m3
         591             :     real(r8) :: sflx(pcols)          ! deposition flux
         592             :     real(r8) :: dep_trb(pcols)       !kg/m2/s
         593             :     real(r8) :: dep_grv(pcols)       !kg/m2/s (total of grav and trb)
         594             :     real(r8) :: pvmzaer(pcols,pverp) ! sedimentation velocity in Pa
         595             :     real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species
         596             :
         597             :     real(r8) :: rad_drop(pcols,pver)
         598             :     real(r8) :: dens_drop(pcols,pver)
         599             :     real(r8) :: sg_drop(pcols,pver)
         600             :     real(r8) :: rad_aer(pcols,pver)
         601             :     real(r8) :: dens_aer(pcols,pver)
         602             :     real(r8) :: sg_aer(pcols,pver)
         603             :
         604             :     real(r8) :: vlc_dry(pcols,pver,4)     ! dep velocity
         605             :     real(r8) :: vlc_grv(pcols,pver,4)     ! dep velocity
         606             :     real(r8)::  vlc_trb(pcols,4)          ! dep velocity
         607             :     real(r8) :: aerdepdryis(pcols,pcnst)  ! aerosol dry deposition (interstitial)
         608             :     real(r8) :: aerdepdrycw(pcols,pcnst)  ! aerosol dry deposition (cloud water)
         609     1489176 :     real(r8), pointer :: fldcw(:,:)
         610     1489176 :     real(r8), pointer :: dgncur_awet(:,:,:)
         611     1489176 :     real(r8), pointer :: wetdens(:,:,:)
         612     1489176 :     real(r8), pointer :: qaerwat(:,:,:)
         613             :
         614     1489176 :     landfrac => cam_in%landfrac(:)
         615     1489176 :     icefrac  => cam_in%icefrac(:)
         616     1489176 :     ocnfrac  => cam_in%ocnfrac(:)
         617     1489176 :     fvin     => cam_in%fv(:)
         618     1489176 :     ram1in   => cam_in%ram1(:)
         619             :
         620     1489176 :     lchnk = state%lchnk
         621     1489176 :     ncol  = state%ncol
         622             :
         623             :     ! calc ram and fv over ocean and sea ice ...
         624             :     call calcram( ncol,landfrac,icefrac,ocnfrac,obklen,&
         625             :                   ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
         626     1489176 :                   state%pdel(:,pver),fvin,fv)
         627             :
         628     1489176 :     call outfld( 'airFV', fv(:), pcols, lchnk )
         629     1489176 :     call outfld( 'RAM1', ram1(:), pcols, lchnk )
         630             :
         631             :     ! note that tendencies are not only in sfc layer (because of sedimentation)
         632             :     ! and that ptend is updated within each subroutine for different species
         633             :
         634     1489176 :     call physics_ptend_init(ptend, state%psetcols, 'aero_model_drydep', lq=drydep_lq)
         635             :
         636     5956704 :     call pbuf_get_field(pbuf, dgnumwet_idx,   dgncur_awet, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
         637     5956704 :     call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens,     start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
         638     5956704 :     call pbuf_get_field(pbuf, qaerwat_idx,    qaerwat,     start=(/1,1,1/), kount=(/pcols,pver,nmodes/) )
         639             :
         640  2314006344 :     tvs(:ncol,:) = state%t(:ncol,:)!*(1+state%q(:ncol,k)
         641  2314006344 :     rho(:ncol,:)=  state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
         642             :
         643             : !
         644             : ! calc settling/deposition velocities for cloud droplets (and cloud-borne aerosols)
         645             : !
         646             : ! *** mean drop radius should eventually be computed from ndrop and qcldwtr
         647  2355876432 :     rad_drop(:,:) = 5.0e-6_r8
         648  2355876432 :     dens_drop(:,:) = rhoh2o
         649  2355876432 :     sg_drop(:,:) = 1.46_r8
         650     1489176 :     jvlc = 3
         651             :     call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv,  &
         652             :                      vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
         653     1489176 :                      rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 0, lchnk)
         654     1489176 :     jvlc = 4
         655             :     call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv,  &
         656             :                      vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
         657     1489176 :                      rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 3, lchnk)
         658             :
         659             :
         660             :
         661     7445880 :     do m = 1, ntot_amode   ! main loop over aerosol modes
         662             :
         663    19359288 :        do lphase = 1, 2   ! loop over interstitial / cloud-borne forms
         664             :
         665    11913408 :           if (lphase == 1) then   ! interstial aerosol - calc settling/dep velocities of mode
         666             :
         667             : ! rad_aer = volume mean wet radius (m)
         668             : ! dgncur_awet = geometric mean wet diameter for number distribution (m)
         669    11913408 :              rad_aer(1:ncol,:) = 0.5_r8*dgncur_awet(1:ncol,:,m)   &
         670  9267938784 :                                  *exp(1.5_r8*(alnsg_amode(m)**2))
         671             : ! dens_aer(1:ncol,:) = wet density (kg/m3)
         672  9256025376 :              dens_aer(1:ncol,:) = wetdens(1:ncol,:,m)
         673  9256025376 :              sg_aer(1:ncol,:) = sigmag_amode(m)
         674             :
         675     5956704 :              jvlc = 1
         676             :              call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv,  &
         677             :                         vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
         678     5956704 :                         rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 0, lchnk)
         679     5956704 :              jvlc = 2
         680             :              call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv,  &
         681             :                         vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
         682     5956704 :                         rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 3, lchnk)
         683             :           end if
         684             :
         685    86372208 :           do lspec = 0, nspec_amode(m)+1   ! loop over number + constituents + water
         686             :
         687    68502096 :              if (lspec == 0) then   ! number
         688    11913408 :                 if (lphase == 1) then
         689     5956704 :                    mm = numptr_amode(m)
         690     5956704 :                    jvlc = 1
         691             :                 else
         692     5956704 :                    mm = numptrcw_amode(m)
         693     5956704 :                    jvlc = 3
         694             :                 endif
         695    56588688 :              else if (lspec <= nspec_amode(m)) then   ! non-water mass
         696    44675280 :                 if (lphase == 1) then
         697    22337640 :                    mm = lmassptr_amode(lspec,m)
         698    22337640 :                    jvlc = 2
         699             :                 else
         700    22337640 :                    mm = lmassptrcw_amode(lspec,m)
         701    22337640 :                    jvlc = 4
         702             :                 endif
         703             :              else   ! water mass
         704             : !   bypass dry deposition of aerosol water
         705             :                 cycle
         706             :                 if (lphase == 1) then
         707             :                    mm = 0
         708             : !                  mm = lwaterptr_amode(m)
         709             :                    jvlc = 2
         710             :                 else
         711             :                    mm = 0
         712             :                    jvlc = 4
         713             :                 endif
         714             :              endif
         715             :
         716             :
         717    56588688 :           if (mm <= 0) cycle
         718             :
         719             : !         if (lphase == 1) then
         720    68502096 :           if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
         721    28294344 :              ptend%lq(mm) = .TRUE.
         722             :
         723             :              ! use pvprogseasalts instead (means making the top level 0)
         724   472449744 :              pvmzaer(:ncol,1)=0._r8
         725 43966120536 :              pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
         726             :
         727    28294344 :              call outfld( trim(cnst_name(mm))//'DDV', pvmzaer(:,2:pverp), pcols, lchnk )
         728             :
         729             :              if(.true.) then ! use phil's method
         730             :              !      convert from meters/sec to pascals/sec
         731             :              !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
         732 43966120536 :                 pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
         733             :
         734             :              !      calculate the tendencies and sfc fluxes from the above velocities
         735             :                 call dust_sediment_tend( &
         736             :                      ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
         737    28294344 :                      state%q(:,:,mm),  pvmzaer,  ptend%q(:,:,mm), sflx  )
         738             :              else   !use charlie's method
         739             :                 call d3ddflux( ncol, vlc_dry(:,:,jvlc), state%q(:,:,mm), state%pmid, &
         740             :                                state%pdel, tvs, sflx, ptend%q(:,:,mm), dt )
         741             :              endif
         742             :
         743             :              ! apportion dry deposition into turb and gravitational settling for tapes
         744    28294344 :              dep_trb = 0._r8
         745    28294344 :              dep_grv = 0._r8
         746   472449744 :              do i=1,ncol
         747   472449744 :                 if (vlc_dry(i,pver,jvlc) /= 0._r8) then
         748   444155400 :                    dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
         749   444155400 :                    dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
         750             :                 end if
         751             :              enddo
         752             :
         753    28294344 :              call outfld( trim(cnst_name(mm))//'DDF', sflx, pcols, lchnk)
         754    28294344 :              call outfld( trim(cnst_name(mm))//'TBF', dep_trb, pcols, lchnk )
         755    28294344 :              call outfld( trim(cnst_name(mm))//'GVF', dep_grv, pcols, lchnk )
         756    28294344 :              call outfld( trim(cnst_name(mm))//'DTQ', ptend%q(:,:,mm), pcols, lchnk)
         757   472449744 :              aerdepdryis(:ncol,mm) = sflx(:ncol)
         758             :
         759    28294344 :           else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then  ! aerosol water
         760             :              ! use pvprogseasalts instead (means making the top level 0)
         761           0 :              pvmzaer(:ncol,1)=0._r8
         762           0 :              pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
         763             :
         764             :              if(.true.) then ! use phil's method
         765             :              !      convert from meters/sec to pascals/sec
         766             :              !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
         767           0 :                 pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
         768             :
         769             :              !      calculate the tendencies and sfc fluxes from the above velocities
         770             :                 call dust_sediment_tend( &
         771             :                      ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
         772           0 :                      qaerwat(:,:,mm),  pvmzaer,  dqdt_tmp(:,:), sflx  )
         773             :              else   !use charlie's method
         774             :                 call d3ddflux( ncol, vlc_dry(:,:,jvlc), qaerwat(:,:,mm), state%pmid, &
         775             :                                state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
         776             :              endif
         777             :
         778             :              ! apportion dry deposition into turb and gravitational settling for tapes
         779           0 :              dep_trb = 0._r8
         780           0 :              dep_grv = 0._r8
         781           0 :              do i=1,ncol
         782           0 :                 if (vlc_dry(i,pver,jvlc) /= 0._r8) then
         783           0 :                    dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
         784           0 :                    dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
         785             :                 end if
         786             :              enddo
         787             :
         788           0 :              qaerwat(1:ncol,:,mm) = qaerwat(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) * dt
         789             :
         790             :           else  ! lphase == 2
         791             :              ! use pvprogseasalts instead (means making the top level 0)
         792   472449744 :              pvmzaer(:ncol,1)=0._r8
         793 43966120536 :              pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
         794    28294344 :              fldcw => qqcw_get_field(pbuf, mm,lchnk)
         795             :
         796             :              if(.true.) then ! use phil's method
         797             :              !      convert from meters/sec to pascals/sec
         798             :              !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
         799 43966120536 :                 pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit
         800             :
         801             :              !      calculate the tendencies and sfc fluxes from the above velocities
         802             :                 call dust_sediment_tend( &
         803             :                      ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
         804    28294344 :                      fldcw(:,:),  pvmzaer,  dqdt_tmp(:,:), sflx  )
         805             :              else   !use charlie's method
         806             :                 call d3ddflux( ncol, vlc_dry(:,:,jvlc), fldcw(:,:), state%pmid, &
         807             :                                state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
         808             :              endif
         809             :
         810             :              ! apportion dry deposition into turb and gravitational settling for tapes
         811    28294344 :              dep_trb = 0._r8
         812    28294344 :              dep_grv = 0._r8
         813   472449744 :              do i=1,ncol
         814   472449744 :                 if (vlc_dry(i,pver,jvlc) /= 0._r8) then
         815   444155400 :                    dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
         816   444155400 :                    dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
         817             :                 end if
         818             :              enddo
         819             :
         820 43966120536 :              fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
         821             :
         822    28294344 :              call outfld( trim(cnst_name_cw(mm))//'DDF', sflx, pcols, lchnk)
         823    28294344 :              call outfld( trim(cnst_name_cw(mm))//'TBF', dep_trb, pcols, lchnk )
         824    28294344 :              call outfld( trim(cnst_name_cw(mm))//'GVF', dep_grv, pcols, lchnk )
         825   472449744 :              aerdepdrycw(:ncol,mm) = sflx(:ncol)
         826             :
         827             :           endif
         828             :
         829             :           enddo   ! lspec = 0, nspec_amode(m)+1
         830             :        enddo   ! lphase = 1, 2
         831             :     enddo   ! m = 1, ntot_amode
         832             :
         833             :     ! if the user has specified prescribed aerosol dep fluxes then
         834             :     ! do not set cam_out dep fluxes according to the prognostic aerosols
         835     1489176 :     if (.not.aerodep_flx_prescribed()) then
         836     1489176 :        call aero_deposition_cam_setdry(aerdepdryis, aerdepdrycw, cam_out)
         837             :     endif
         838             :
         839     2978352 :   endsubroutine aero_model_drydep
         840             :
         841             :   !=============================================================================
         842             :   !=============================================================================
         843    62545392 :   subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf)
         844             :
         845     1489176 :     use aero_wetdep_cam, only: aero_wetdep_tend
         846             :
         847             :     ! args
         848             :
         849             :     type(physics_state), intent(in)    :: state       ! Physics state variables
         850             :     real(r8),            intent(in)    :: dt          ! time step
         851             :     real(r8),            intent(in)    :: dlf(:,:)    ! shallow+deep convective detrainment [kg/kg/s]
         852             :     type(cam_out_t),     intent(inout) :: cam_out     ! export state
         853             :     type(physics_ptend), intent(out)   :: ptend       ! indivdual parameterization tendencies
         854             :     type(physics_buffer_desc), pointer :: pbuf(:)
         855             :
         856     1489176 :     call aero_wetdep_tend(state, dt, dlf, cam_out, ptend, pbuf)
         857             :
         858     1489176 :   end subroutine aero_model_wetdep
         859             :
         860             :   !-------------------------------------------------------------------------
         861             :   ! provides wet tropospheric aerosol surface area info for modal aerosols
         862             :   ! called from mo_usrrxt
         863             :   !-------------------------------------------------------------------------
         864           0 :   subroutine aero_model_surfarea( &
         865           0 :                   mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, rho, ltrop, &
         866           0 :                   dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_trop, reff_trop )
         867             :
         868             :     ! dummy args
         869             :     real(r8), intent(in)    :: pmid(:,:)
         870             :     real(r8), intent(in)    :: temp(:,:)
         871             :     real(r8), intent(in)    :: mmr(:,:,:)
         872             :     real(r8), intent(in)    :: radmean      ! mean radii in cm
         873             :     real(r8), intent(in)    :: strato_sad(:,:)
         874             :     integer,  intent(in)    :: ncol
         875             :     integer,  intent(in)    :: ltrop(:)
         876             :     real(r8), intent(in)    :: dlat(:)                    ! degrees latitude
         877             :     integer,  intent(in)    :: het1_ndx
         878             :     real(r8), intent(in)    :: relhum(:,:)
         879             :     real(r8), intent(in)    :: rho(:,:) ! total atm density (/cm^3)
         880             :     real(r8), intent(in)    :: sulfate(:,:)
         881             :     type(physics_buffer_desc), pointer :: pbuf(:)
         882             :
         883             :     real(r8), intent(inout) :: sfc(:,:,:)
         884             :     real(r8), intent(inout) :: dm_aer(:,:,:)
         885             :     real(r8), intent(inout) :: sad_trop(:,:)
         886             :     real(r8), intent(out)   :: reff_trop(:,:)
         887             :
         888             :     ! local vars
         889           0 :     real(r8), pointer, dimension(:,:,:) :: dgnumwet
         890           0 :     integer :: beglev(ncol)
         891           0 :     integer :: endlev(ncol)
         892             :     integer :: i,k
         893             :
         894           0 :     call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
         895             :
         896           0 :     beglev(:ncol)=ltrop(:ncol)+1
         897           0 :     endlev(:ncol)=pver
         898           0 :     call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, sad_trop, reff_trop, sfc=sfc )
         899             :
         900           0 :     do i = 1,ncol
         901           0 :        do k = ltrop(i)+1,pver
         902           0 :           dm_aer(i,k,:) = dgnumwet(i,k,:) * 1.e2_r8 ! convert m to cm
         903             :        enddo
         904             :     enddo
         905             :
         906     1489176 :   end subroutine aero_model_surfarea
         907             :
         908             :   !-------------------------------------------------------------------------
         909             :   ! provides WET stratospheric aerosol surface area info for modal aerosols
         910             :   ! if modal_strat_sulfate = TRUE -- called from mo_gas_phase_chemdr
         911             :   !-------------------------------------------------------------------------
         912           0 :   subroutine aero_model_strat_surfarea( ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat )
         913             :
         914             :     ! dummy args
         915             :     integer,  intent(in)    :: ncol
         916             :     real(r8), intent(in)    :: mmr(:,:,:)
         917             :     real(r8), intent(in)    :: pmid(:,:)
         918             :     real(r8), intent(in)    :: temp(:,:)
         919             :     integer,  intent(in)    :: ltrop(:) ! tropopause level indices
         920             :     type(physics_buffer_desc), pointer :: pbuf(:)
         921             :     real(r8), intent(out)   :: strato_sad(:,:)
         922             :     real(r8), intent(out)   :: reff_strat(:,:)
         923             :
         924             :     ! local vars
         925           0 :     real(r8), pointer, dimension(:,:,:) :: dgnumwet
         926           0 :     integer :: beglev(ncol)
         927           0 :     integer :: endlev(ncol)
         928             :
         929           0 :     reff_strat = 0._r8
         930           0 :     strato_sad = 0._r8
         931             :
         932           0 :     if (.not.modal_strat_sulfate) return
         933             :
         934           0 :     call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet )
         935             :
         936           0 :     beglev(:ncol)=top_lev
         937           0 :     endlev(:ncol)=ltrop(:ncol)
         938           0 :     call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, strato_sad, reff_strat )
         939             :
         940           0 :   end subroutine aero_model_strat_surfarea
         941             :
         942             :   !=============================================================================
         943             :   !=============================================================================
         944     1489176 :   subroutine aero_model_gasaerexch( loffset, ncol, lchnk, troplev, delt, reaction_rates, &
         945     1489176 :                                     tfld, pmid, pdel, mbar, relhum, &
         946     2978352 :                                     zm,  qh2o, cwat, cldfr, cldnum, &
         947     2978352 :                                     airdens, invariants, del_h2so4_gasprod,  &
         948     1489176 :                                     vmr0, vmr, pbuf )
         949             :
         950             :     use time_manager,          only : get_nstep
         951             :     use modal_aero_coag,       only : modal_aero_coag_sub
         952             :     use modal_aero_gasaerexch, only : modal_aero_gasaerexch_sub
         953             :     use modal_aero_newnuc,     only : modal_aero_newnuc_sub
         954             :     use modal_aero_data,       only : cnst_name_cw, qqcw_get_field
         955             :
         956             :     !-----------------------------------------------------------------------
         957             :     !      ... dummy arguments
         958             :     !-----------------------------------------------------------------------
         959             :     integer,  intent(in) :: loffset                ! offset applied to modal aero "pointers"
         960             :     integer,  intent(in) :: ncol                   ! number columns in chunk
         961             :     integer,  intent(in) :: lchnk                  ! chunk index
         962             :     integer,  intent(in) :: troplev(pcols)
         963             :     real(r8), intent(in) :: delt                   ! time step size (sec)
         964             :     real(r8), intent(in) :: reaction_rates(:,:,:)  ! reaction rates
         965             :     real(r8), intent(in) :: tfld(:,:)              ! temperature (K)
         966             :     real(r8), intent(in) :: pmid(:,:)              ! pressure at model levels (Pa)
         967             :     real(r8), intent(in) :: pdel(:,:)              ! pressure thickness of levels (Pa)
         968             :     real(r8), intent(in) :: mbar(:,:)              ! mean wet atmospheric mass ( amu )
         969             :     real(r8), intent(in) :: relhum(:,:)            ! relative humidity
         970             :     real(r8), intent(in) :: airdens(:,:)           ! total atms density (molec/cm**3)
         971             :     real(r8), intent(in) :: invariants(:,:,:)
         972             :     real(r8), intent(in) :: del_h2so4_gasprod(:,:)
         973             :     real(r8), intent(in) :: zm(:,:)
         974             :     real(r8), intent(in) :: qh2o(:,:)
         975             :     real(r8), intent(in) :: cwat(:,:)          ! cloud liquid water content (kg/kg)
         976             :     real(r8), intent(in) :: cldfr(:,:)
         977             :     real(r8), intent(in) :: cldnum(:,:)       ! droplet number concentration (#/kg)
         978             :     real(r8), intent(in) :: vmr0(:,:,:)       ! initial mixing ratios (before gas-phase chem changes)
         979             :     real(r8), intent(inout) :: vmr(:,:,:)         ! mixing ratios ( vmr )
         980             :
         981             :     type(physics_buffer_desc), pointer :: pbuf(:)
         982             :
         983             :     ! local vars
         984             :
         985             :     integer :: n, m
         986             :     integer :: i,k,l
         987             :     integer :: nstep
         988             :
         989     2978352 :     real(r8) :: del_h2so4_aeruptk(ncol,pver)
         990             :
         991     1489176 :     real(r8), pointer :: dgnum(:,:,:), dgnumwet(:,:,:), wetdens(:,:,:)
         992     1489176 :     real(r8), pointer :: pblh(:)                    ! pbl height (m)
         993             :
         994     2978352 :     real(r8), dimension(ncol) :: wrk
         995             :     character(len=32)         :: name
         996     2978352 :     real(r8) :: dvmrcwdt(ncol,pver,gas_pcnst)
         997     2978352 :     real(r8) :: dvmrdt(ncol,pver,gas_pcnst)
         998     2978352 :     real(r8) :: vmrcw(ncol,pver,gas_pcnst)            ! cloud-borne aerosol (vmr)
         999             :
        1000     2978352 :     real(r8) ::  aqso4(ncol,ntot_amode)               ! aqueous phase chemistry
        1001     2978352 :     real(r8) ::  aqh2so4(ncol,ntot_amode)             ! aqueous phase chemistry
        1002     2978352 :     real(r8) ::  aqso4_h2o2(ncol)                     ! SO4 aqueous phase chemistry due to H2O2
        1003     2978352 :     real(r8) ::  aqso4_o3(ncol)                       ! SO4 aqueous phase chemistry due to O3
        1004     2978352 :     real(r8) ::  xphlwc(ncol,pver)                    ! pH value multiplied by lwc
        1005     2978352 :     real(r8) ::  nh3_beg(ncol,pver)
        1006     1489176 :     real(r8), pointer :: fldcw(:,:)
        1007     1489176 :     real(r8), pointer :: sulfeq(:,:,:)
        1008             :
        1009             :     logical :: is_spcam_m2005
        1010             : !
        1011             : ! ... initialize nh3
        1012             : !
        1013     1489176 :     if ( nh3_ndx > 0 ) then
        1014           0 :       nh3_beg = vmr(1:ncol,:,nh3_ndx)
        1015             :     end if
        1016             : !
        1017     1489176 :     is_spcam_m2005   = cam_physpkg_is('spcam_m2005')
        1018             :
        1019     1489176 :     call pbuf_get_field(pbuf, dgnum_idx,      dgnum)
        1020     1489176 :     call pbuf_get_field(pbuf, dgnumwet_idx,   dgnumwet )
        1021     1489176 :     call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens )
        1022     1489176 :     call pbuf_get_field(pbuf, pblh_idx,       pblh)
        1023             :
        1024     7445880 :     do n=1,ntot_amode
        1025     5956704 :        call outfld(dgnum_name(n), dgnum(1:ncol,1:pver,n), ncol, lchnk )
        1026     7445880 :        call outfld(dgnumwet_name(n), dgnumwet(1:ncol,1:pver,n), ncol, lchnk )
        1027             :     end do
        1028             :
        1029             : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
        1030             :
        1031     1489176 :     nstep = get_nstep()
        1032             :
        1033             :     ! calculate tendency due to gas phase chemistry and processes
        1034 71735685840 :     dvmrdt(:ncol,:,:) = (vmr(:ncol,:,:) - vmr0(:ncol,:,:)) / delt
        1035    47653632 :     do m = 1, gas_pcnst
        1036   770839056 :       wrk(:) = 0._r8
        1037  4339458864 :       do k = 1,pver
        1038 71734196664 :         wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m)*adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
        1039             :       end do
        1040    46164456 :       name = 'GS_'//trim(solsym(m))
        1041    47653632 :       call outfld( name, wrk(:ncol), ncol, lchnk )
        1042             :     enddo
        1043             :
        1044             : !
        1045             : ! Aerosol processes ...
        1046             : !
        1047     1489176 :     call qqcw2vmr( lchnk, vmrcw, mbar, ncol, loffset, pbuf )
        1048             :
        1049     1489176 :     if (.not. is_spcam_m2005) then  ! regular CAM
        1050 71735685840 :       dvmrdt(:ncol,:,:) = vmr(:ncol,:,:)
        1051 71735685840 :       dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:)
        1052             :
        1053             :     ! aqueous chemistry ...
        1054             :
        1055     1489176 :       if( has_sox ) then
        1056             :          call setsox(   &
        1057             :               ncol,     &
        1058             :               lchnk,    &
        1059             :               loffset,  &
        1060             :               delt,     &
        1061             :               pmid,     &
        1062             :               pdel,     &
        1063             :               tfld,     &
        1064             :               mbar,     &
        1065             :               cwat,     &
        1066             :               cldfr,    &
        1067             :               cldnum,   &
        1068             :               airdens,  &
        1069             :               invariants, &
        1070             :               vmrcw,    &
        1071             :               vmr,      &
        1072             :               xphlwc,   &
        1073             :               aqso4,    &
        1074             :               aqh2so4,  &
        1075             :               aqso4_h2o2, &
        1076             :               aqso4_o3  &
        1077     1489176 :               )
        1078             :
        1079     7445880 :          do n = 1, ntot_amode
        1080     5956704 :             l = lptr_so4_cw_amode(n)
        1081     7445880 :             if (l > 0) then
        1082     4467528 :                call outfld( trim(cnst_name_cw(l))//'AQSO4',   aqso4(:ncol,n),   ncol, lchnk)
        1083     4467528 :                call outfld( trim(cnst_name_cw(l))//'AQH2SO4', aqh2so4(:ncol,n), ncol, lchnk)
        1084             :             end if
        1085             :          end do
        1086             :
        1087     1489176 :          call outfld( 'AQSO4_H2O2', aqso4_h2o2(:ncol), ncol, lchnk)
        1088     1489176 :          call outfld( 'AQSO4_O3',   aqso4_o3(:ncol),   ncol, lchnk)
        1089     1489176 :          call outfld( 'XPH_LWC',    xphlwc(:ncol,:),   ncol, lchnk )
        1090             :
        1091             :       endif
        1092             :
        1093             : !   Tendency due to aqueous chemistry
        1094 71737175016 :     dvmrdt = (vmr - dvmrdt) / delt
        1095 71735685840 :     dvmrcwdt = (vmrcw - dvmrcwdt) / delt
        1096    47653632 :     do m = 1, gas_pcnst
        1097   770839056 :       wrk(:) = 0._r8
        1098  4339458864 :       do k = 1,pver
        1099 71734196664 :         wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m) * adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
        1100             :       end do
        1101    46164456 :       name = 'AQ_'//trim(solsym(m))
        1102    47653632 :       call outfld( name, wrk(:ncol), ncol, lchnk )
        1103             :     enddo
        1104             :
        1105             :    else if (is_spcam_m2005) then  ! SPCAM ECPP
        1106             : ! when ECPP is used, aqueous chemistry is done in ECPP,
        1107             : ! and not updated here.
        1108             : ! Minghuai Wang, 2010-02 (Minghuai.Wang@pnl.gov)
        1109             : !
        1110           0 :       dvmrdt = 0.0_r8
        1111           0 :       dvmrcwdt = 0.0_r8
        1112             :    endif
        1113             :
        1114             : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
        1115             :
        1116     1489176 :     if (ndx_h2so4 > 0) then
        1117  2314006344 :        del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4)
        1118             :     else
        1119           0 :        del_h2so4_aeruptk(:,:) = 0.0_r8
        1120             :     endif
        1121             :
        1122     1489176 :     call t_startf('modal_gas-aer_exchng')
        1123             :
        1124     1489176 :     if ( sulfeq_idx>0 ) then
        1125           0 :        call pbuf_get_field( pbuf, sulfeq_idx, sulfeq )
        1126             :     else
        1127     1489176 :        nullify( sulfeq )
        1128             :     endif
        1129             :
        1130             :     call modal_aero_gasaerexch_sub(            &
        1131             :          lchnk,    ncol,     nstep,            &
        1132             :          loffset,            delt,             &
        1133             :          tfld,     pmid,     pdel,             &
        1134             :          qh2o,               troplev,          &
        1135             :          vmr,                vmrcw,            &
        1136             :          dvmrdt,             dvmrcwdt,         &
        1137             :          dgnum,              dgnumwet,         &
        1138     1489176 :          sulfeq     )
        1139             :
        1140     1489176 :     if (ndx_h2so4 > 0) then
        1141  2314006344 :        del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) - del_h2so4_aeruptk(1:ncol,:)
        1142             :     endif
        1143             :
        1144     1489176 :     call t_stopf('modal_gas-aer_exchng')
        1145             :
        1146     1489176 :     call t_startf('modal_nucl')
        1147             :
        1148             :     ! do aerosol nucleation (new particle formation)
        1149             :     call modal_aero_newnuc_sub(                             &
        1150             :          lchnk,    ncol,     nstep,            &
        1151             :          loffset,            delt,             &
        1152             :          tfld,     pmid,     pdel,             &
        1153             :          zm,       pblh,                       &
        1154             :          qh2o,     cldfr,                      &
        1155             :          vmr,                                  &
        1156   543653136 :          del_h2so4_gasprod,  del_h2so4_aeruptk )
        1157             :
        1158     1489176 :     call t_stopf('modal_nucl')
        1159             :
        1160     1489176 :     call t_startf('modal_coag')
        1161             :
        1162             :     ! do aerosol coagulation
        1163             :     call modal_aero_coag_sub(                               &
        1164             :          lchnk,    ncol,     nstep,            &
        1165             :          loffset,            delt,             &
        1166             :          tfld,     pmid,     pdel,             &
        1167             :          vmr,                                  &
        1168             :          dgnum,              dgnumwet,         &
        1169     1489176 :          wetdens                          )
        1170             :
        1171     1489176 :     call t_stopf('modal_coag')
        1172             :
        1173     1489176 :     call vmr2qqcw( lchnk, vmrcw, mbar, ncol, loffset, pbuf )
        1174             :
        1175             :     ! diagnostics for cloud-borne aerosols...
        1176    62545392 :     do n = 1,pcnst
        1177    61056216 :        fldcw => qqcw_get_field(pbuf,n,lchnk,errorhandle=.true.)
        1178    62545392 :        if(associated(fldcw)) then
        1179    28294344 :           call outfld( cnst_name_cw(n), fldcw(:,:), pcols, lchnk )
        1180             :        endif
        1181             :     end do
        1182             : !
        1183             : ! ... put missing NH3 into NH4
        1184             : !
        1185     1489176 :     if ( nh3_ndx > 0 .and. nh4_ndx > 0 ) then
        1186           0 :       vmr(1:ncol,:,nh4_ndx) = vmr(1:ncol,:,nh4_ndx) + (nh3_beg-vmr(1:ncol,:,nh3_ndx))
        1187           0 :       vmr(1:ncol,:,nh4_ndx) = max(0._r8,vmr(1:ncol,:,nh4_ndx))
        1188             :     end if
        1189             :
        1190     2978352 :   end subroutine aero_model_gasaerexch
        1191             :
        1192             :   !=============================================================================
        1193             :   !=============================================================================
        1194     1489176 :   subroutine aero_model_emissions( state, cam_in )
        1195     1489176 :     use seasalt_model, only: seasalt_emis, seasalt_names, seasalt_indices, seasalt_active,seasalt_nbin
        1196             :     use dust_model,    only: dust_emis, dust_names, dust_indices, dust_active,dust_nbin, dust_nnum
        1197             :     use physics_types, only: physics_state
        1198             :
        1199             :     ! Arguments:
        1200             :
        1201             :     type(physics_state),    intent(in)    :: state   ! Physics state variables
        1202             :     type(cam_in_t),         intent(inout) :: cam_in  ! import state
        1203             :
        1204             :     ! local vars
        1205             :
        1206             :     integer :: lchnk, ncol
        1207             :     integer :: m, mm
        1208             :     real(r8) :: soil_erod_tmp(pcols)
        1209             :     real(r8) :: sflx(pcols)   ! accumulate over all bins for output
        1210             :     real(r8) :: u10cubed(pcols)
        1211             :     real (r8), parameter :: z0=0.0001_r8  ! m roughness length over oceans--from ocean model
        1212             :
        1213     1489176 :     lchnk = state%lchnk
        1214     1489176 :     ncol = state%ncol
        1215             :
        1216     1489176 :     if (dust_active) then
        1217             :
        1218     1489176 :        call dust_emis( ncol, lchnk, cam_in%dstflx, cam_in%cflx, soil_erod_tmp )
        1219             :
        1220             :        ! some dust emis diagnostics ...
        1221     1489176 :        sflx(:)=0._r8
        1222    10424232 :        do m=1,dust_nbin+dust_nnum
        1223     8935056 :           mm = dust_indices(m)
        1224    79064856 :           if (m<=dust_nbin) sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
        1225    10424232 :           call outfld(trim(dust_names(m))//'SF',cam_in%cflx(:,mm),pcols, lchnk)
        1226             :        enddo
        1227     1489176 :        call outfld('DSTSFMBL',sflx(:),pcols,lchnk)
        1228     1489176 :        call outfld('LND_MBL',soil_erod_tmp(:),pcols, lchnk )
        1229             :     endif
        1230             :
        1231     1489176 :     if (seasalt_active) then
        1232    24865776 :        u10cubed(:ncol)=sqrt(state%u(:ncol,pver)**2+state%v(:ncol,pver)**2)
        1233             :        ! move the winds to 10m high from the midpoint of the gridbox:
        1234             :        ! follows Tie and Seinfeld and Pandis, p.859 with math.
        1235             :
        1236    24865776 :        u10cubed(:ncol)=u10cubed(:ncol)*log(10._r8/z0)/log(state%zm(:ncol,pver)/z0)
        1237             :
        1238             :        ! we need them to the 3.41 power, according to Gong et al., 1997:
        1239    24865776 :        u10cubed(:ncol)=u10cubed(:ncol)**3.41_r8
        1240             :
        1241     1489176 :        sflx(:)=0._r8
        1242             :
        1243     1489176 :        call seasalt_emis( u10cubed, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx )
        1244             :
        1245     5956704 :        do m=1,seasalt_nbin
        1246     4467528 :           mm = seasalt_indices(m)
        1247    74597328 :           sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
        1248     5956704 :           call outfld(trim(seasalt_names(m))//'SF',cam_in%cflx(:,mm),pcols,lchnk)
        1249             :        enddo
        1250     1489176 :        call outfld('SSTSFMBL',sflx(:),pcols,lchnk)
        1251             :     endif
        1252             :
        1253     1489176 :   end subroutine aero_model_emissions
        1254             :
        1255             :   !===============================================================================
        1256             :   ! private methods
        1257             :
        1258             :
        1259             :   !=============================================================================
        1260             :   !=============================================================================
        1261           0 :   subroutine surf_area_dens( ncol, mmr, pmid, temp, diam, beglev, endlev, sad, reff, sfc )
        1262     1489176 :     use mo_constants,    only : pi
        1263             :     use modal_aero_data, only : nspec_amode, alnsg_amode
        1264             :
        1265             :     ! dummy args
        1266             :     integer,  intent(in)  :: ncol
        1267             :     real(r8), intent(in)  :: mmr(:,:,:)
        1268             :     real(r8), intent(in)  :: pmid(:,:)
        1269             :     real(r8), intent(in)  :: temp(:,:)
        1270             :     real(r8), intent(in)  :: diam(:,:,:)
        1271             :     integer,  intent(in)  :: beglev(:)
        1272             :     integer,  intent(in)  :: endlev(:)
        1273             :     real(r8), intent(out) :: sad(:,:)
        1274             :     real(r8), intent(out) :: reff(:,:)
        1275             :     real(r8),optional, intent(out) :: sfc(:,:,:)
        1276             :
        1277             :     ! local vars
        1278           0 :     real(r8) :: sad_mode(pcols,pver,ntot_amode),radeff(pcols,pver)
        1279           0 :     real(r8) :: vol(pcols,pver),vol_mode(pcols,pver,ntot_amode)
        1280             :     real(r8) :: rho_air
        1281             :     integer  :: i,k,l,m
        1282             :     real(r8) :: chm_mass, tot_mass
        1283             :
        1284             :     !
        1285             :     ! Compute surface aero for each mode.
        1286             :     ! Total over all modes as the surface area for chemical reactions.
        1287             :     !
        1288             :
        1289           0 :     sad = 0._r8
        1290           0 :     sad_mode = 0._r8
        1291           0 :     vol = 0._r8
        1292           0 :     vol_mode = 0._r8
        1293           0 :     reff = 0._r8
        1294             :
        1295           0 :     do i = 1,ncol
        1296           0 :        do k = beglev(i),endlev(i)
        1297           0 :           rho_air = pmid(i,k)/(temp(i,k)*287.04_r8)
        1298           0 :           do l=1,ntot_amode
        1299             :              !
        1300             :              ! compute a mass weighting of the number
        1301             :              !
        1302           0 :              tot_mass = 0._r8
        1303           0 :              chm_mass = 0._r8
        1304           0 :              do m=1,nspec_amode(l)
        1305           0 :                if ( index_tot_mass(l,m) > 0 ) &
        1306           0 :                     tot_mass = tot_mass + mmr(i,k,index_tot_mass(l,m))
        1307           0 :                if ( index_chm_mass(l,m) > 0 ) &
        1308           0 :                     chm_mass = chm_mass + mmr(i,k,index_chm_mass(l,m))
        1309             :              end do
        1310           0 :              if ( tot_mass > 0._r8 ) then
        1311             :               ! surface area density
        1312           0 :                sad_mode(i,k,l) = chm_mass/tot_mass &
        1313           0 :                                * mmr(i,k,num_idx(l))*rho_air*pi*diam(i,k,l)**2._r8 &
        1314           0 :                                * exp(2._r8*alnsg_amode(l)**2._r8)  ! m^2/m^3
        1315           0 :                sad_mode(i,k,l) = 1.e-2_r8 * sad_mode(i,k,l) ! cm^2/cm^3
        1316             :
        1317             :               ! volume calculation, for use in effective radius calculation
        1318             :                vol_mode(i,k,l) = chm_mass/tot_mass &
        1319           0 :                                * mmr(i,k,num_idx(l))*rho_air*pi/6._r8*diam(i,k,l)**3._r8  &
        1320           0 :                                * exp(4.5_r8*alnsg_amode(l)**2._r8)  ! m^3/m^3 = cm^3/cm^3
        1321             :              else
        1322           0 :                sad_mode(i,k,l) = 0._r8
        1323           0 :                vol_mode(i,k,l) = 0._r8
        1324             :              end if
        1325             :           end do
        1326           0 :           sad(i,k) = sum(sad_mode(i,k,:))
        1327           0 :           vol(i,k) = sum(vol_mode(i,k,:))
        1328           0 :           reff(i,k) = 3._r8*vol(i,k)/sad(i,k)
        1329             :
        1330             :        enddo
        1331             :     enddo
        1332             :
        1333           0 :     if (present(sfc)) then
        1334           0 :        sfc(:,:,:) = sad_mode(:,:,:)
        1335             :     endif
        1336             :
        1337           0 :   end subroutine surf_area_dens
        1338             :
        1339             :   !===============================================================================
        1340             :   !===============================================================================
        1341        1536 :   subroutine modal_aero_bcscavcoef_init
        1342             :     !-----------------------------------------------------------------------
        1343             :     !
        1344             :     ! Purpose:
        1345             :     ! Computes lookup table for aerosol impaction/interception scavenging rates
        1346             :     !
        1347             :     ! Authors: R. Easter
        1348             :     !
        1349             :     !-----------------------------------------------------------------------
        1350             :
        1351           0 :     use shr_kind_mod,    only: r8 => shr_kind_r8
        1352             :     use modal_aero_data
        1353             :     use cam_abortutils,  only: endrun
        1354             :
        1355             :     implicit none
        1356             :
        1357             :
        1358             :     !   local variables
        1359             :     integer nnfit_maxd
        1360             :     parameter (nnfit_maxd=27)
        1361             :
        1362             :     integer i, jgrow, jdens, jpress, jtemp, mode, nnfit
        1363             :     integer lunerr
        1364             :
        1365             :     real(r8) dg0, dg0_cgs, press, &
        1366             :          rhodryaero, rhowetaero, rhowetaero_cgs, rmserr, &
        1367             :          scavratenum, scavratevol, sigmag,                &
        1368             :          temp, wetdiaratio, wetvolratio
        1369             :     real(r8) aafitnum(1), xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd)
        1370             :     real(r8) aafitvol(1), xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd)
        1371             :
        1372        4608 :     allocate(scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode))
        1373        3072 :     allocate(scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode))
        1374             :
        1375        1536 :     lunerr = 6
        1376        1536 :     dlndg_nimptblgrow = log( 1.25_r8 )
        1377             :
        1378        7680 :     modeloop: do mode = 1, ntot_amode
        1379             :
        1380        6144 :        sigmag = sigmag_amode(mode)
        1381             :
        1382        6144 :        rhodryaero = specdens_amode(1,mode)
        1383             :
        1384      130560 :        growloop: do jgrow = nimptblgrow_mind, nimptblgrow_maxd
        1385             :
        1386      122880 :           wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
        1387      122880 :           dg0 = dgnum_amode(mode)*wetdiaratio
        1388             :
        1389             :           wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 )
        1390      122880 :           rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio
        1391      122880 :           rhowetaero = min( rhowetaero, rhodryaero )
        1392             :
        1393             :           !
        1394             :           !   compute impaction scavenging rates at 1 temp-press pair and save
        1395             :           !
        1396      122880 :           nnfit = 0
        1397             :
        1398      122880 :           temp = 273.16_r8
        1399      122880 :           press = 0.75e6_r8   ! dynes/cm2
        1400      122880 :           rhowetaero = rhodryaero
        1401             :
        1402      122880 :           dg0_cgs = dg0*1.0e2_r8   ! m to cm
        1403      122880 :           rhowetaero_cgs = rhowetaero*1.0e-3_r8   ! kg/m3 to g/cm3
        1404             :           call calc_1_impact_rate( &
        1405             :                dg0_cgs, sigmag, rhowetaero_cgs, temp, press, &
        1406      122880 :                scavratenum, scavratevol, lunerr )
        1407             :
        1408      122880 :           nnfit = nnfit + 1
        1409             :           if (nnfit > nnfit_maxd) then
        1410             :              write(lunerr,9110)
        1411             :              call endrun()
        1412             :           end if
        1413             : 9110      format( '*** subr. modal_aero_bcscavcoef_init -- nnfit too big' )
        1414             :
        1415      122880 :           xxfitnum(1,nnfit) = 1._r8
        1416      122880 :           yyfitnum(nnfit) = log( scavratenum )
        1417             :
        1418      122880 :           xxfitvol(1,nnfit) = 1._r8
        1419      122880 :           yyfitvol(nnfit) = log( scavratevol )
        1420             :
        1421             : 5900      continue
        1422             :
        1423             :           !
        1424             :           ! skip mlinfit stuff because scav table no longer has dependencies on
        1425             :           !    air temp, air press, and particle wet density
        1426             :           ! just load the log( scavrate--- ) values
        1427             :           !
        1428             :           !!
        1429             :           !!   do linear regression
        1430             :           !!    log(scavrate) = a1 + a2*log(wetdens)
        1431             :           !!
        1432             :           !     call mlinft( xxfitnum, yyfitnum, aafitnum, nnfit, 1, 1, rmserr )
        1433             :           !     call mlinft( xxfitvol, yyfitvol, aafitvol, nnfit, 1, 1, rmserr )
        1434             :           !
        1435             :           !     scavimptblnum(jgrow,mode) = aafitnum(1)
        1436             :           !     scavimptblvol(jgrow,mode) = aafitvol(1)
        1437             :
        1438      122880 :           scavimptblnum(jgrow,mode) = yyfitnum(1)
        1439      129024 :           scavimptblvol(jgrow,mode) = yyfitvol(1)
        1440             :
        1441             :        enddo growloop
        1442             :     enddo modeloop
        1443        1536 :     return
        1444        1536 :   end subroutine modal_aero_bcscavcoef_init
        1445             :
        1446             :   !===============================================================================
        1447             :   !===============================================================================
        1448    14891760 :   subroutine modal_aero_depvel_part( ncol, t, pmid, ram1, fv, vlc_dry, vlc_trb, vlc_grv,  &
        1449             :                                      radius_part, density_part, sig_part, moment, lchnk )
        1450             :
        1451             : !    calculates surface deposition velocity of particles
        1452             : !    L. Zhang, S. Gong, J. Padro, and L. Barrie
        1453             : !    A size-seggregated particle dry deposition scheme for an atmospheric aerosol module
        1454             : !    Atmospheric Environment, 35, 549-560, 2001.
        1455             : !
        1456             : !    Authors: X. Liu
        1457             :
        1458             :     !
        1459             :     ! !USES
        1460             :     !
        1461        1536 :     use physconst,     only: pi,boltz, gravit, rair
        1462             :     use mo_drydep,     only: n_land_type, fraction_landuse
        1463             :
        1464             :     ! !ARGUMENTS:
        1465             :     !
        1466             :     implicit none
        1467             :     !
        1468             :     real(r8), intent(in) :: t(pcols,pver)       !atm temperature (K)
        1469             :     real(r8), intent(in) :: pmid(pcols,pver)    !atm pressure (Pa)
        1470             :     real(r8), intent(in) :: fv(pcols)           !friction velocity (m/s)
        1471             :     real(r8), intent(in) :: ram1(pcols)         !aerodynamical resistance (s/m)
        1472             :     real(r8), intent(in) :: radius_part(pcols,pver)    ! mean (volume/number) particle radius (m)
        1473             :     real(r8), intent(in) :: density_part(pcols,pver)   ! density of particle material (kg/m3)
        1474             :     real(r8), intent(in) :: sig_part(pcols,pver)       ! geometric standard deviation of particles
        1475             :     integer,  intent(in) :: moment ! moment of size distribution (0 for number, 2 for surface area, 3 for volume)
        1476             :     integer,  intent(in) :: ncol
        1477             :     integer,  intent(in) :: lchnk
        1478             :
        1479             :     real(r8), intent(out) :: vlc_trb(pcols)       !Turbulent deposn velocity (m/s)
        1480             :     real(r8), intent(out) :: vlc_grv(pcols,pver)       !grav deposn velocity (m/s)
        1481             :     real(r8), intent(out) :: vlc_dry(pcols,pver)       !dry deposn velocity (m/s)
        1482             :     !------------------------------------------------------------------------
        1483             :
        1484             :     !------------------------------------------------------------------------
        1485             :     ! Local Variables
        1486             :     integer  :: m,i,k,ix                !indices
        1487             :     real(r8) :: rho     !atm density (kg/m**3)
        1488             :     real(r8) :: vsc_dyn_atm(pcols,pver)   ![kg m-1 s-1] Dynamic viscosity of air
        1489             :     real(r8) :: vsc_knm_atm(pcols,pver)   ![m2 s-1] Kinematic viscosity of atmosphere
        1490             :     real(r8) :: shm_nbr       ![frc] Schmidt number
        1491             :     real(r8) :: stk_nbr       ![frc] Stokes number
        1492             :     real(r8) :: mfp_atm(pcols,pver)       ![m] Mean free path of air
        1493             :     real(r8) :: dff_aer       ![m2 s-1] Brownian diffusivity of particle
        1494             :     real(r8) :: slp_crc(pcols,pver) ![frc] Slip correction factor
        1495             :     real(r8) :: rss_trb       ![s m-1] Resistance to turbulent deposition
        1496             :     real(r8) :: rss_lmn       ![s m-1] Quasi-laminar layer resistance
        1497             :     real(r8) :: brownian      ! collection efficiency for Browning diffusion
        1498             :     real(r8) :: impaction     ! collection efficiency for impaction
        1499             :     real(r8) :: interception  ! collection efficiency for interception
        1500             :     real(r8) :: stickfrac     ! fraction of particles sticking to surface
        1501             :     real(r8) :: radius_moment(pcols,pver) ! median radius (m) for moment
        1502             :     real(r8) :: lnsig         ! ln(sig_part)
        1503             :     real(r8) :: dispersion    ! accounts for influence of size dist dispersion on bulk settling velocity
        1504             :                               ! assuming radius_part is number mode radius * exp(1.5 ln(sigma))
        1505             :
        1506             :     integer  :: lt
        1507             :     real(r8) :: lnd_frc
        1508             :     real(r8) :: wrk1, wrk2, wrk3
        1509             :
        1510             :     ! constants
        1511             :     real(r8) gamma(11)      ! exponent of schmidt number
        1512             : !   data gamma/0.54d+00,  0.56d+00,  0.57d+00,  0.54d+00,  0.54d+00, &
        1513             : !              0.56d+00,  0.54d+00,  0.54d+00,  0.54d+00,  0.56d+00, &
        1514             : !              0.50d+00/
        1515             :     data gamma/0.56e+00_r8,  0.54e+00_r8,  0.54e+00_r8,  0.56e+00_r8,  0.56e+00_r8, &
        1516             :                0.56e+00_r8,  0.50e+00_r8,  0.54e+00_r8,  0.54e+00_r8,  0.54e+00_r8, &
        1517             :                0.54e+00_r8/
        1518             :     save gamma
        1519             :
        1520             :     real(r8) alpha(11)      ! parameter for impaction
        1521             : !   data alpha/50.00d+00,  0.95d+00,  0.80d+00,  1.20d+00,  1.30d+00, &
        1522             : !               0.80d+00, 50.00d+00, 50.00d+00,  2.00d+00,  1.50d+00, &
        1523             : !             100.00d+00/
        1524             :     data alpha/1.50e+00_r8,   1.20e+00_r8,  1.20e+00_r8,  0.80e+00_r8,  1.00e+00_r8, &
        1525             :                0.80e+00_r8, 100.00e+00_r8, 50.00e+00_r8,  2.00e+00_r8,  1.20e+00_r8, &
        1526             :               50.00e+00_r8/
        1527             :     save alpha
        1528             :
        1529             :     real(r8) radius_collector(11) ! radius (m) of surface collectors
        1530             : !   data radius_collector/-1.00d+00,  5.10d-03,  3.50d-03,  3.20d-03, 10.00d-03, &
        1531             : !                          5.00d-03, -1.00d+00, -1.00d+00, 10.00d-03, 10.00d-03, &
        1532             : !                         -1.00d+00/
        1533             :     data radius_collector/10.00e-03_r8,  3.50e-03_r8,  3.50e-03_r8,  5.10e-03_r8,  2.00e-03_r8, &
        1534             :                            5.00e-03_r8, -1.00e+00_r8, -1.00e+00_r8, 10.00e-03_r8,  3.50e-03_r8, &
        1535             :                           -1.00e+00_r8/
        1536             :     save radius_collector
        1537             :
        1538             :     integer            :: iwet(11) ! flag for wet surface = 1, otherwise = -1
        1539             : !   data iwet/1,   -1,   -1,   -1,   -1,  &
        1540             : !            -1,   -1,   -1,    1,   -1,  &
        1541             : !             1/
        1542             :     data iwet/-1,  -1,   -1,   -1,   -1,  &
        1543             :               -1,   1,   -1,    1,   -1,  &
        1544             :               -1/
        1545             :     save iwet
        1546             :
        1547             :
        1548    14891760 :     vlc_trb = 0._r8
        1549    14891760 :     vlc_grv = 0._r8
        1550    14891760 :     vlc_dry = 0._r8
        1551             :
        1552             :     !------------------------------------------------------------------------
        1553  1399825440 :     do k=top_lev,pver ! radius_part is not defined above top_lev
        1554 23140063440 :        do i=1,ncol
        1555             :
        1556 21740238000 :           lnsig = log(sig_part(i,k))
        1557             : ! use a maximum radius of 50 microns when calculating deposition velocity
        1558             :           radius_moment(i,k) = min(50.0e-6_r8,radius_part(i,k))*   &
        1559 21740238000 :                           exp((float(moment)-1.5_r8)*lnsig*lnsig)
        1560 21740238000 :           dispersion = exp(2._r8*lnsig*lnsig)
        1561             :
        1562 21740238000 :           rho=pmid(i,k)/rair/t(i,k)
        1563             :
        1564             :           ! Quasi-laminar layer resistance: call rss_lmn_get
        1565             :           ! Size-independent thermokinetic properties
        1566             :           vsc_dyn_atm(i,k) = 1.72e-5_r8 * ((t(i,k)/273.0_r8)**1.5_r8) * 393.0_r8 / &
        1567 21740238000 :                (t(i,k)+120.0_r8)      ![kg m-1 s-1] RoY94 p. 102
        1568             :           mfp_atm(i,k) = 2.0_r8 * vsc_dyn_atm(i,k) / &   ![m] SeP97 p. 455
        1569 21740238000 :                (pmid(i,k)*sqrt(8.0_r8/(pi*rair*t(i,k))))
        1570 21740238000 :           vsc_knm_atm(i,k) = vsc_dyn_atm(i,k) / rho ![m2 s-1] Kinematic viscosity of air
        1571             :
        1572             :           slp_crc(i,k) = 1.0_r8 + mfp_atm(i,k) * &
        1573             :                   (1.257_r8+0.4_r8*exp(-1.1_r8*radius_moment(i,k)/(mfp_atm(i,k)))) / &
        1574 21740238000 :                   radius_moment(i,k)   ![frc] Slip correction factor SeP97 p. 464
        1575             :           vlc_grv(i,k) = (4.0_r8/18.0_r8) * radius_moment(i,k)*radius_moment(i,k)*density_part(i,k)* &
        1576 21740238000 :                   gravit*slp_crc(i,k) / vsc_dyn_atm(i,k) ![m s-1] Stokes' settling velocity SeP97 p. 466
        1577 21740238000 :           vlc_grv(i,k) = vlc_grv(i,k) * dispersion
        1578             :
        1579 23125171680 :           vlc_dry(i,k)=vlc_grv(i,k)
        1580             :        enddo
        1581             :     enddo
        1582    14891760 :     k=pver  ! only look at bottom level for next part
        1583   248657760 :     do i=1,ncol
        1584   233766000 :        dff_aer = boltz * t(i,k) * slp_crc(i,k) / &    ![m2 s-1]
        1585   233766000 :                  (6.0_r8*pi*vsc_dyn_atm(i,k)*radius_moment(i,k)) !SeP97 p.474
        1586   233766000 :        shm_nbr = vsc_knm_atm(i,k) / dff_aer                        ![frc] SeP97 p.972
        1587             :
        1588   233766000 :        wrk2 = 0._r8
        1589   233766000 :        wrk3 = 0._r8
        1590  2805192000 :        do lt = 1,n_land_type
        1591  2571426000 :           lnd_frc = fraction_landuse(i,lt,lchnk)
        1592  2805192000 :           if ( lnd_frc /= 0._r8 ) then
        1593   570692070 :              brownian = shm_nbr**(-gamma(lt))
        1594   570692070 :              if (radius_collector(lt) > 0.0_r8) then
        1595             : !       vegetated surface
        1596   259898730 :                 stk_nbr = vlc_grv(i,k) * fv(i) / (gravit*radius_collector(lt))
        1597   259898730 :                 interception = 2.0_r8*(radius_moment(i,k)/radius_collector(lt))**2.0_r8
        1598             :              else
        1599             : !       non-vegetated surface
        1600   310793340 :                 stk_nbr = vlc_grv(i,k) * fv(i) * fv(i) / (gravit*vsc_knm_atm(i,k))  ![frc] SeP97 p.965
        1601   310793340 :                 interception = 0.0_r8
        1602             :              endif
        1603   570692070 :              impaction = (stk_nbr/(alpha(lt)+stk_nbr))**2.0_r8
        1604             :
        1605   570692070 :              if (iwet(lt) > 0) then
        1606             :                 stickfrac = 1.0_r8
        1607             :              else
        1608   353636010 :                 stickfrac = exp(-sqrt(stk_nbr))
        1609   353636010 :                 if (stickfrac < 1.0e-10_r8) stickfrac = 1.0e-10_r8
        1610             :              endif
        1611   570692070 :              rss_lmn = 1.0_r8 / (3.0_r8 * fv(i) * stickfrac * (brownian+interception+impaction))
        1612   570692070 :              rss_trb = ram1(i) + rss_lmn + ram1(i)*rss_lmn*vlc_grv(i,k)
        1613             :
        1614   570692070 :              wrk1 = 1.0_r8 / rss_trb
        1615   570692070 :              wrk2 = wrk2 + lnd_frc*( wrk1 )
        1616   570692070 :              wrk3 = wrk3 + lnd_frc*( wrk1 + vlc_grv(i,k) )
        1617             :           endif
        1618             :        enddo  ! n_land_type
        1619   233766000 :        vlc_trb(i) = wrk2
        1620   248657760 :        vlc_dry(i,k) = wrk3
        1621             :     enddo !ncol
        1622             :
        1623    14891760 :     return
        1624    14891760 :   end subroutine modal_aero_depvel_part
        1625             :
        1626             :   !===============================================================================
        1627             :   subroutine modal_aero_bcscavcoef_get( m, ncol, isprx, dgn_awet, scavcoefnum, scavcoefvol )
        1628             :
        1629    14891760 :     use modal_aero_data
        1630             :     !-----------------------------------------------------------------------
        1631             :     implicit none
        1632             :
        1633             :     integer,intent(in) :: m, ncol
        1634             :     logical,intent(in):: isprx(pcols,pver)
        1635             :     real(r8), intent(in) :: dgn_awet(pcols,pver,ntot_amode)
        1636             :     real(r8), intent(out) :: scavcoefnum(pcols,pver), scavcoefvol(pcols,pver)
        1637             :
        1638             :     integer i, k, jgrow
        1639             :     real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum
        1640             :
        1641             :
        1642             :     do k = 1, pver
        1643             :        do i = 1, ncol
        1644             :
        1645             :           ! do only if no precip
        1646             :           if ( isprx(i,k) ) then
        1647             :              !
        1648             :              ! interpolate table values using log of (actual-wet-size)/(base-dry-size)
        1649             :
        1650             :              dumdgratio = dgn_awet(i,k,m)/dgnum_amode(m)
        1651             :
        1652             :              if ((dumdgratio >= 0.99_r8) .and. (dumdgratio <= 1.01_r8)) then
        1653             :                 scavimpvol = scavimptblvol(0,m)
        1654             :                 scavimpnum = scavimptblnum(0,m)
        1655             :              else
        1656             :                 xgrow = log( dumdgratio ) / dlndg_nimptblgrow
        1657             :                 jgrow = int( xgrow )
        1658             :                 if (xgrow < 0._r8) jgrow = jgrow - 1
        1659             :                 if (jgrow < nimptblgrow_mind) then
        1660             :                    jgrow = nimptblgrow_mind
        1661             :                    xgrow = jgrow
        1662             :                 else
        1663             :                    jgrow = min( jgrow, nimptblgrow_maxd-1 )
        1664             :                 end if
        1665             :
        1666             :                 dumfhi = xgrow - jgrow
        1667             :                 dumflo = 1._r8 - dumfhi
        1668             :
        1669             :                 scavimpvol = dumflo*scavimptblvol(jgrow,m) + &
        1670             :                      dumfhi*scavimptblvol(jgrow+1,m)
        1671             :                 scavimpnum = dumflo*scavimptblnum(jgrow,m) + &
        1672             :                      dumfhi*scavimptblnum(jgrow+1,m)
        1673             :
        1674             :              end if
        1675             :
        1676             :              ! impaction scavenging removal amount for volume
        1677             :              scavcoefvol(i,k) = exp( scavimpvol )
        1678             :              ! impaction scavenging removal amount to number
        1679             :              scavcoefnum(i,k) = exp( scavimpnum )
        1680             :
        1681             :              ! scavcoef = impaction scav rate (1/h) for precip = 1 mm/h
        1682             :              ! scavcoef = impaction scav rate (1/s) for precip = pfx_inrain
        1683             :              ! (scavcoef/3600) = impaction scav rate (1/s) for precip = 1 mm/h
        1684             :              ! (pfx_inrain*3600) = in-rain-area precip rate (mm/h)
        1685             :              ! impactrate = (scavcoef/3600) * (pfx_inrain*3600)
        1686             :           else
        1687             :              scavcoefvol(i,k) = 0._r8
        1688             :              scavcoefnum(i,k) = 0._r8
        1689             :           end if
        1690             :
        1691             :        end do
        1692             :     end do
        1693             :
        1694             :     return
        1695             :   end subroutine modal_aero_bcscavcoef_get
        1696             :
        1697             :   !===============================================================================
        1698      122880 :         subroutine calc_1_impact_rate(             &
        1699             :                 dg0, sigmag, rhoaero, temp, press, &
        1700             :                 scavratenum, scavratevol, lunerr )
        1701             :    !
        1702             :    !   routine computes a single impaction scavenging rate
        1703             :    !    for precipitation rate of 1 mm/h
        1704             :    !
        1705             :    !   dg0 = geometric mean diameter of aerosol number size distrib. (cm)
        1706             :    !   sigmag = geometric standard deviation of size distrib.
        1707             :    !   rhoaero = density of aerosol particles (g/cm^3)
        1708             :    !   temp = temperature (K)
        1709             :    !   press = pressure (dyne/cm^2)
        1710             :    !   scavratenum = number scavenging rate (1/h)
        1711             :    !   scavratevol = volume or mass scavenging rate (1/h)
        1712             :    !   lunerr = logical unit for error message
        1713             :    !
        1714             :    use shr_kind_mod, only: r8 => shr_kind_r8
        1715             :    use mo_constants, only: boltz_cgs, pi, rhowater => rhoh2o_cgs, &
        1716             :                            gravity => gravity_cgs, rgas => rgas_cgs
        1717             :
        1718             :    implicit none
        1719             :
        1720             :    !   subr. parameters
        1721             :    integer lunerr
        1722             :    real(r8) dg0, sigmag, rhoaero, temp, press, scavratenum, scavratevol
        1723             :
        1724             :    !   local variables
        1725             :    integer nrainsvmax
        1726             :    parameter (nrainsvmax=50)
        1727             :    real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
        1728             :         vfallrainsv(nrainsvmax)
        1729             :
        1730             :    integer naerosvmax
        1731             :    parameter (naerosvmax=51)
        1732             :    real(r8) aaerosv(naerosvmax), &
        1733             :         ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)
        1734             :
        1735             :    integer i, ja, jr, na, nr
        1736             :    real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
        1737             :    real(r8) anumsum, avolsum, cair, chi
        1738             :    real(r8) d, dr, dum, dumfuchs, dx
        1739             :    real(r8) ebrown, eimpact, eintercept, etotal, freepath
        1740             :    real(r8) precip, precipmmhr, precipsum
        1741             :    real(r8) r, rainsweepout, reynolds, rhi, rhoair, rlo, rnumsum
        1742             :    real(r8) scavsumnum, scavsumnumbb
        1743             :    real(r8) scavsumvol, scavsumvolbb
        1744             :    real(r8) schmidt, sqrtreynolds, sstar, stokes, sx
        1745             :    real(r8) taurelax, vfall, vfallstp
        1746             :    real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair
        1747             :
        1748             :
        1749      122880 :    rlo = .005_r8
        1750      122880 :    rhi = .250_r8
        1751      122880 :    dr = 0.005_r8
        1752      122880 :    nr = 1 + nint( (rhi-rlo)/dr )
        1753             :    if (nr > nrainsvmax) then
        1754             :       write(lunerr,9110)
        1755             :       call endrun()
        1756             :    end if
        1757             :
        1758             : 9110 format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' )
        1759             :
        1760      122880 :    precipmmhr = 1.0_r8
        1761      122880 :    precip = precipmmhr/36000._r8
        1762             :
        1763      122880 :    ag0 = dg0/2._r8
        1764      122880 :    sx = log( sigmag )
        1765      122880 :    xg0 = log( ag0 )
        1766      122880 :    xg3 = xg0 + 3._r8*sx*sx
        1767             :
        1768      122880 :    xlo = xg3 - 4._r8*sx
        1769      122880 :    xhi = xg3 + 4._r8*sx
        1770      122880 :    dx = 0.2_r8*sx
        1771             :
        1772      122880 :    dx = max( 0.2_r8*sx, 0.01_r8 )
        1773      122880 :    xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
        1774      122880 :    xhi = xg3 + max( 4._r8*sx, 2._r8*dx )
        1775             :
        1776      122880 :    na = 1 + nint( (xhi-xlo)/dx )
        1777      122880 :    if (na > naerosvmax) then
        1778           0 :       write(lunerr,9120)
        1779           0 :       call endrun()
        1780             :    end if
        1781             :
        1782             : 9120 format( '*** subr. calc_1_impact_rate -- na > naerosvmax' )
        1783             :
        1784             :    !   air molar density
        1785      122880 :    cair = press/(rgas*temp)
        1786             :    !   air mass density
        1787      122880 :    rhoair = 28.966_r8*cair
        1788             :    !   molecular freepath
        1789      122880 :    freepath = 2.8052e-10_r8/cair
        1790             :    !   air dynamic viscosity
        1791             :    airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) *    &
        1792      122880 :         ((temp/296.16_r8)**1.5_r8)
        1793             :    !   air kinemaic viscosity
        1794      122880 :    airkinvisc = airdynvisc/rhoair
        1795             :    !   ratio of water viscosity to air viscosity (from Slinn)
        1796      122880 :    xmuwaterair = 60.0_r8
        1797             :
        1798             :    !
        1799             :    !   compute rain drop number concentrations
        1800             :    !    rrainsv = raindrop radius (cm)
        1801             :    !    xnumrainsv = raindrop number concentration (#/cm^3)
        1802             :    !            (number in the bin, not number density)
        1803             :    !    vfallrainsv = fall velocity (cm/s)
        1804             :    !
        1805      122880 :    precipsum = 0._r8
        1806     6266880 :    do i = 1, nr
        1807     6144000 :       r = rlo + (i-1)*dr
        1808     6144000 :       rrainsv(i) = r
        1809     6144000 :       xnumrainsv(i) = exp( -r/2.7e-2_r8 )
        1810             :
        1811     6144000 :       d = 2._r8*r
        1812     6144000 :       if (d <= 0.007_r8) then
        1813           0 :          vfallstp = 2.88e5_r8 * d**2._r8
        1814     6144000 :       else if (d <= 0.025_r8) then
        1815      245760 :          vfallstp = 2.8008e4_r8 * d**1.528_r8
        1816     5898240 :       else if (d <= 0.1_r8) then
        1817      983040 :          vfallstp = 4104.9_r8 * d**1.008_r8
        1818     4915200 :       else if (d <= 0.25_r8) then
        1819     1843200 :          vfallstp = 1812.1_r8 * d**0.638_r8
        1820             :       else
        1821     3072000 :          vfallstp = 1069.8_r8 * d**0.235_r8
        1822             :       end if
        1823             :
        1824     6144000 :       vfall = vfallstp * sqrt(1.204e-3_r8/rhoair)
        1825     6144000 :       vfallrainsv(i) = vfall
        1826     6266880 :       precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
        1827             :    end do
        1828      122880 :    precipsum = precipsum*pi*1.333333_r8
        1829             :
        1830      122880 :    rnumsum = 0._r8
        1831     6266880 :    do i = 1, nr
        1832     6144000 :       xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
        1833     6266880 :       rnumsum = rnumsum + xnumrainsv(i)
        1834             :    end do
        1835             :
        1836             :    !
        1837             :    !   compute aerosol concentrations
        1838             :    !    aaerosv = particle radius (cm)
        1839             :    !    fnumaerosv = fraction of total number in the bin (--)
        1840             :    !    fvolaerosv = fraction of total volume in the bin (--)
        1841             :    !
        1842             :    anumsum = 0._r8
        1843             :    avolsum = 0._r8
        1844     5160960 :    do i = 1, na
        1845     5038080 :       x = xlo + (i-1)*dx
        1846     5038080 :       a = exp( x )
        1847     5038080 :       aaerosv(i) = a
        1848     5038080 :       dum = (x - xg0)/sx
        1849     5038080 :       ynumaerosv(i) = exp( -0.5_r8*dum*dum )
        1850     5038080 :       yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a
        1851     5038080 :       anumsum = anumsum + ynumaerosv(i)
        1852     5160960 :       avolsum = avolsum + yvolaerosv(i)
        1853             :    end do
        1854             :
        1855     5160960 :    do i = 1, na
        1856     5038080 :       ynumaerosv(i) = ynumaerosv(i)/anumsum
        1857     5160960 :       yvolaerosv(i) = yvolaerosv(i)/avolsum
        1858             :    end do
        1859             :
        1860             :
        1861             :    !
        1862             :    !   compute scavenging
        1863             :    !
        1864             :    scavsumnum = 0._r8
        1865             :    scavsumvol = 0._r8
        1866             :    !
        1867             :    !   outer loop for rain drop radius
        1868             :    !
        1869     6266880 :    jr_loop: do jr = 1, nr
        1870             :
        1871     6144000 :       r = rrainsv(jr)
        1872     6144000 :       vfall = vfallrainsv(jr)
        1873             :
        1874     6144000 :       reynolds = r * vfall / airkinvisc
        1875     6144000 :       sqrtreynolds = sqrt( reynolds )
        1876             :
        1877             :       !
        1878             :       !   inner loop for aerosol particle radius
        1879             :       !
        1880     6144000 :       scavsumnumbb = 0._r8
        1881     6144000 :       scavsumvolbb = 0._r8
        1882             :
        1883   258048000 :       ja_loop: do ja = 1, na
        1884             :
        1885   251904000 :          a = aaerosv(ja)
        1886             :
        1887   251904000 :          chi = a/r
        1888             :
        1889   251904000 :          dum = freepath/a
        1890   251904000 :          dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum)
        1891   251904000 :          taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc)
        1892             :
        1893   251904000 :          aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8
        1894   251904000 :          aerodiffus = boltz_cgs*temp*taurelax/aeromass
        1895             :
        1896   251904000 :          schmidt = airkinvisc/aerodiffus
        1897   251904000 :          stokes = vfall*taurelax/r
        1898             :
        1899             :          ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) /  &
        1900   251904000 :               (reynolds*schmidt)
        1901             :
        1902             :          dum = (1._r8 + 2._r8*xmuwaterair*chi) /         &
        1903   251904000 :               (1._r8 + xmuwaterair/sqrtreynolds)
        1904   251904000 :          eintercept = 4._r8*chi*(chi + dum)
        1905             :
        1906   251904000 :          dum = log( 1._r8 + reynolds )
        1907   251904000 :          sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum)
        1908   251904000 :          eimpact = 0._r8
        1909   251904000 :          if (stokes > sstar) then
        1910    61905408 :             dum = stokes - sstar
        1911    61905408 :             eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8
        1912             :          end if
        1913             :
        1914   251904000 :          etotal = ebrown + eintercept + eimpact
        1915   251904000 :          etotal = min( etotal, 1.0_r8 )
        1916             :
        1917   251904000 :          rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall
        1918             :
        1919   251904000 :          scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
        1920   258048000 :          scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
        1921             :
        1922             :       enddo ja_loop
        1923             :
        1924     6144000 :       scavsumnum = scavsumnum + scavsumnumbb
        1925     6266880 :       scavsumvol = scavsumvol + scavsumvolbb
        1926             :
        1927             :    enddo jr_loop
        1928             :
        1929      122880 :    scavratenum = scavsumnum*3600._r8
        1930      122880 :    scavratevol = scavsumvol*3600._r8
        1931             :
        1932      122880 :    return
        1933             :  end subroutine calc_1_impact_rate
        1934             :
        1935             :   !=============================================================================
        1936             :   !=============================================================================
        1937     2978352 :   subroutine qqcw2vmr(lchnk, vmr, mbar, ncol, im, pbuf)
        1938             :     use modal_aero_data, only : qqcw_get_field
        1939             :     use physics_buffer, only : physics_buffer_desc
        1940             :     !-----------------------------------------------------------------
        1941             :     !   ... Xfrom from mass to volume mixing ratio
        1942             :     !-----------------------------------------------------------------
        1943             :
        1944             :     use chem_mods, only : adv_mass, gas_pcnst
        1945             :
        1946             :     implicit none
        1947             :
        1948             :     !-----------------------------------------------------------------
        1949             :     !   ... Dummy args
        1950             :     !-----------------------------------------------------------------
        1951             :     integer, intent(in)     :: lchnk, ncol, im
        1952             :     real(r8), intent(in)    :: mbar(ncol,pver)
        1953             :     real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst)
        1954             :     type(physics_buffer_desc), pointer :: pbuf(:)
        1955             :
        1956             :     !-----------------------------------------------------------------
        1957             :     !   ... Local variables
        1958             :     !-----------------------------------------------------------------
        1959             :     integer :: k, m
        1960     1489176 :     real(r8), pointer :: fldcw(:,:)
        1961             :
        1962    47653632 :     do m=1,gas_pcnst
        1963    47653632 :        if( adv_mass(m) /= 0._r8 ) then
        1964    46164456 :           fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.)
        1965    46164456 :           if(associated(fldcw)) then
        1966  2659668336 :              do k=1,pver
        1967 43966120536 :                 vmr(:ncol,k,m) = mbar(:ncol,k) * fldcw(:ncol,k) / adv_mass(m)
        1968             :              end do
        1969             :           else
        1970 27768076128 :              vmr(:,:,m) = 0.0_r8
        1971             :           end if
        1972             :        end if
        1973             :     end do
        1974     2978352 :   end subroutine qqcw2vmr
        1975             :
        1976             :
        1977             :   !=============================================================================
        1978             :   !=============================================================================
        1979     2978352 :   subroutine vmr2qqcw( lchnk, vmr, mbar, ncol, im, pbuf )
        1980             :     !-----------------------------------------------------------------
        1981             :     !   ... Xfrom from volume to mass mixing ratio
        1982             :     !-----------------------------------------------------------------
        1983             :
        1984     1489176 :     use m_spc_id
        1985             :     use chem_mods,       only : adv_mass, gas_pcnst
        1986             :     use modal_aero_data, only : qqcw_get_field
        1987             :     use physics_buffer,  only : physics_buffer_desc
        1988             :
        1989             :     implicit none
        1990             :
        1991             :     !-----------------------------------------------------------------
        1992             :     !   ... Dummy args
        1993             :     !-----------------------------------------------------------------
        1994             :     integer, intent(in)     :: lchnk, ncol, im
        1995             :     real(r8), intent(in)    :: mbar(ncol,pver)
        1996             :     real(r8), intent(in)    :: vmr(ncol,pver,gas_pcnst)
        1997             :     type(physics_buffer_desc), pointer :: pbuf(:)
        1998             :
        1999             :     !-----------------------------------------------------------------
        2000             :     !   ... Local variables
        2001             :     !-----------------------------------------------------------------
        2002             :     integer :: k, m
        2003     1489176 :     real(r8), pointer :: fldcw(:,:)
        2004             :     !-----------------------------------------------------------------
        2005             :     !   ... The non-group species
        2006             :     !-----------------------------------------------------------------
        2007    47653632 :     do m = 1,gas_pcnst
        2008    46164456 :        fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.)
        2009    47653632 :        if( adv_mass(m) /= 0._r8 .and. associated(fldcw)) then
        2010  2659668336 :           do k = 1,pver
        2011 43966120536 :              fldcw(:ncol,k) = adv_mass(m) * vmr(:ncol,k,m) / mbar(:ncol,k)
        2012             :           end do
        2013             :        end if
        2014             :     end do
        2015             :
        2016     2978352 :   end subroutine vmr2qqcw
        2017             :
        2018           0 :   function get_dlndg_nimptblgrow() result (dlndg_nimptblgrow_ret)
        2019             :     real(r8) ::  dlndg_nimptblgrow_ret
        2020           0 :     dlndg_nimptblgrow_ret =  dlndg_nimptblgrow
        2021     1489176 :   end function get_dlndg_nimptblgrow
        2022             :
        2023           0 :   function get_scavimptblvol() result (scavimptblvol_ret)
        2024             :     real(r8) :: scavimptblvol_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
        2025           0 :     scavimptblvol_ret = scavimptblvol
        2026           0 :   end function get_scavimptblvol
        2027             :
        2028           0 :   function get_scavimptblnum() result (scavimptblnum_ret)
        2029             :     real(r8) :: scavimptblnum_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
        2030           0 :     scavimptblnum_ret = scavimptblnum
        2031           0 :   end function get_scavimptblnum
        2032             :
        2033             : end module aero_model

Generated by: LCOV version 1.14