module mz_aerosols_intr

  use shr_kind_mod,only: r8 => shr_kind_r8
  use spmd_utils,  only: masterproc
  use pmgrid,      only: plev, plevp
  use ppgrid,      only: pcols, pver
  use ref_pres,    only: top_lev => clim_modal_aero_top_lev
  use constituents,only: pcnst, cnst_name, cnst_get_ind
  use phys_control,only: phys_getopts
  use abortutils,  only: endrun
#ifdef MODAL_AERO
    use modal_aero_data, only: ntot_amode
#endif
  use cam_logfile, only: iulog
  use aerodep_flx, only: aerodep_flx_prescribed

  implicit none

  private          ! Make default type private to the module

  save

  !
  ! Public interfaces
  !

  public :: mz_aero_initialize                           
  public :: mz_aero_wet_intr                             ! interface to wet deposition
  public :: mz_aero_setopts
  public :: mz_aero_defaultopts
  public :: sol_facti_cloud_borne
  public :: aer_wetdep_list
  public :: do_cam_sulfchem
#ifdef MODAL_AERO
  public :: mz_aero_dry_intr                             ! interface to dry deposition
  public :: aer_drydep_list
  public :: modal_aero_bcscavcoef_init  ! initialize impaction scavenging table
  public :: modal_aero_bcscavcoef_get   ! get impaction scavenging coefficients
#endif
  public :: mz_prescribed_dust              ! return true if MOZART code provides prescribed dust

  integer :: num_mz_aerosols
  integer, allocatable :: mz_aerosol_ids(:) ! (/ id_cb2, id_oc2, id_so4, id_sa1, id_sa2, id_sa3, id_sa4 /)
  character(len=16), dimension(pcnst) :: aer_wetdep_list = ' '
  logical :: use_cam_sulfchem = .false.
  logical :: do_cam_sulfchem
  logical :: inv_o3, inv_oh, inv_no3, inv_ho2
  integer, pointer :: id_so2, id_so4, id_dms, id_o3, id_h2o2, id_oh, id_no3, id_ho2
  integer, target  :: spc_ids(8)

  real(r8) :: sol_facti_cloud_borne

#ifdef MODAL_AERO
  character(len=16), dimension(pcnst) :: aer_drydep_list = ''
! variables for table lookup of aerosol impaction/interception scavenging rates
 integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
 real(r8) dlndg_nimptblgrow
 real(r8) scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
 real(r8) scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)
#endif

contains

  !===============================================================================
#ifdef MODAL_AERO
  subroutine mz_aero_defaultopts( sol_facti_cloud_borne_out, aer_wetdep_list_out, aer_drydep_list_out, use_cam_sulfchem_out)
#else
  subroutine mz_aero_defaultopts( sol_facti_cloud_borne_out, aer_wetdep_list_out, use_cam_sulfchem_out)
#endif
    implicit none

    real(r8),         intent(out)           :: sol_facti_cloud_borne_out

    character(len=*), intent(out), optional :: aer_wetdep_list_out(:)
#ifdef MODAL_AERO
    character(len=*), intent(out), optional :: aer_drydep_list_out(:)
#endif
    logical,          intent(out), optional :: use_cam_sulfchem_out


    sol_facti_cloud_borne_out = 1.0_r8

#ifdef MODAL_AERO
    if ( present(aer_drydep_list_out) ) then
       aer_drydep_list_out(:) = aer_drydep_list(:)
    endif
#endif
    if ( present(aer_wetdep_list_out) ) then
       aer_wetdep_list_out(:) = aer_wetdep_list(:)
    endif
    if ( present(use_cam_sulfchem_out) ) then
       use_cam_sulfchem_out = use_cam_sulfchem
    endif

  end subroutine mz_aero_defaultopts

  !===============================================================================
#ifdef MODAL_AERO
  subroutine mz_aero_setopts( sol_facti_cloud_borne_in, aer_wetdep_list_in, aer_drydep_list_in, use_cam_sulfchem_in )
#else
  subroutine mz_aero_setopts( sol_facti_cloud_borne_in, aer_wetdep_list_in, use_cam_sulfchem_in )
#endif
    implicit none

    real(r8),         intent(in)           :: sol_facti_cloud_borne_in

    character(len=*), intent(in), optional :: aer_wetdep_list_in(:)
#ifdef MODAL_AERO
    character(len=*), intent(in), optional :: aer_drydep_list_in(:)
#endif
    logical,          intent(in), optional :: use_cam_sulfchem_in

    sol_facti_cloud_borne = sol_facti_cloud_borne_in

#ifdef MODAL_AERO
    if ( present(aer_drydep_list_in) ) then
       aer_drydep_list(:) = aer_drydep_list_in(:)
    endif
#endif
    if ( present(aer_wetdep_list_in) ) then
       aer_wetdep_list(:) = aer_wetdep_list_in(:)
    endif
    if ( present(use_cam_sulfchem_in) ) then
       use_cam_sulfchem = use_cam_sulfchem_in
    endif
    
  end subroutine mz_aero_setopts

  !===============================================================================
  subroutine mz_aero_initialize( )
    use cam_history,      only : addfld, add_default, phys_decomp
    use sulchem,          only : inisulchem
    use mo_chem_utls,     only : get_inv_ndx
    use dust_intr,        only : dst_wet_dep=>dust_has_wet_dep
    use progseasalts_intr,only : sst_wet_dep=>progseasalts_has_wet_dep
    use dust_intr,        only : dust_names
    use progseasalts_intr,only : progseasalts_names
    use gas_wetdep_opts,  only : gas_wetdep_list, gas_wetdep_cnt

    implicit none

    integer :: m                                
    integer :: mm                               
    integer :: astat, id

    logical :: history_aerosol      ! Output the MAM aerosol tendencies
    character(len=2)  :: unit_basename  ! Units 'kg' or '1' 
    !-----------------------------------------------------------------------

    call phys_getopts( history_aerosol_out        = history_aerosol   )

    dst_wet_dep(:) = .false.
    sst_wet_dep(:) = .false.

    num_mz_aerosols = 0

    count_species: do m = 1,pcnst

       if ( len_trim(aer_wetdep_list(m)) == 0 ) exit count_species

       do mm=1,gas_wetdep_cnt
          if ( trim(gas_wetdep_list(mm)) == trim(aer_wetdep_list(m)) ) then
             call endrun('mz_aero_initialize: '//trim(aer_wetdep_list(m))//' is in both het_lst and aer_wetdep_list')
          endif
       enddo

       call cnst_get_ind ( aer_wetdep_list(m), id, abort=.false. )
       if ( id < 1 ) then 
          write(iulog,*) 'mz_aero_initialize: '//trim(aer_wetdep_list(m))//' does not exit in simulation'
          call endrun('mz_aero_initialize: invalid washout species')
       else
          where( aer_wetdep_list(m) == dust_names )
             dst_wet_dep = .true.
          endwhere
          where( aer_wetdep_list(m) == progseasalts_names )
             sst_wet_dep = .true.
          endwhere

          if (masterproc) then
             write(iulog,*) 'mz_aero_initialize: '//aer_wetdep_list(m)//' will have wet removal'
          endif

          ! units 
          if (aer_wetdep_list(m)(1:3) == 'num') then
             unit_basename = ' 1'
          else
             unit_basename = 'kg'  
          endif

          call addfld (trim(aer_wetdep_list(m))//'SFWET',unit_basename//'/m2/s ',1,  'A','Wet deposition flux at surface',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SFSIC',unit_basename//'/m2/s ', &
               1,  'A','Wet deposition flux (incloud, convective) at surface',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SFSIS',unit_basename//'/m2/s ', &
               1,  'A','Wet deposition flux (incloud, stratiform) at surface',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SFSBC',unit_basename//'/m2/s ', &
               1,  'A','Wet deposition flux (belowcloud, convective) at surface',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SFSBS',unit_basename//'/m2/s ', &
               1,  'A','Wet deposition flux (belowcloud, stratiform) at surface',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'WET',unit_basename//'/kg/s ',pver, 'A','wet deposition tendency',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SIC',unit_basename//'/kg/s ',pver, 'A', &
               trim(aer_wetdep_list(m))//' ic wet deposition',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SIS',unit_basename//'/kg/s ',pver, 'A', &
               trim(aer_wetdep_list(m))//' is wet deposition',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SBC',unit_basename//'/kg/s ',pver, 'A', &
               trim(aer_wetdep_list(m))//' bc wet deposition',phys_decomp)
          call addfld (trim(aer_wetdep_list(m))//'SBS',unit_basename//'/kg/s ',pver, 'A', &
               trim(aer_wetdep_list(m))//' bs wet deposition',phys_decomp)
       

          if ( history_aerosol ) then          
             call add_default (trim(aer_wetdep_list(m))//'SFWET', 1, ' ')
             call add_default (trim(aer_wetdep_list(m))//'SFSIC', 1, ' ')
             call add_default (trim(aer_wetdep_list(m))//'SFSIS', 1, ' ')
             call add_default (trim(aer_wetdep_list(m))//'SFSBC', 1, ' ')
             call add_default (trim(aer_wetdep_list(m))//'SFSBS', 1, ' ')
          endif
          num_mz_aerosols = num_mz_aerosols + 1
       endif

    enddo count_species

#ifdef MODAL_AERO
    count_dry_species: do m = 1,pcnst

       if ( len_trim(aer_drydep_list(m)) == 0 ) exit count_dry_species

       call cnst_get_ind ( aer_drydep_list(m), id, abort=.false. )
       if ( id < 1 ) then
          write(*,*) 'mz_aero_initialize: '//trim(aer_drydep_list(m))//' does not exist in simulation'
          call endrun('mz_aero_initialize: invalid drydep species')
       else
          if (masterproc) then
             write(*,*) 'mz_aero_initialize: '//aer_drydep_list(m)//' will have dry deposition'
          endif

          ! units 
          if (aer_drydep_list(m)(1:3) == 'num') then
             unit_basename = ' 1'
          else
             unit_basename = 'kg'  
          endif
 
          call addfld (trim(aer_drydep_list(m))//'DDF',unit_basename//'/m2/s ',   1, 'A', &
               trim(aer_drydep_list(m))//' dry deposition flux at bottom (grav + turb)',phys_decomp)
          call addfld (trim(aer_drydep_list(m))//'TBF',unit_basename//'/m2/s',   1, 'A', &
               trim(aer_drydep_list(m))//' turbulent dry deposition flux',phys_decomp)
          call addfld (trim(aer_drydep_list(m))//'GVF',unit_basename//'/m2/s ',   1, 'A', &
               trim(aer_drydep_list(m))//' gravitational dry deposition flux',phys_decomp)
          call addfld (trim(aer_drydep_list(m))//'DTQ',unit_basename//'/kg/s ',pver, 'A', &
               trim(aer_drydep_list(m))//' dry deposition',phys_decomp)
          call addfld (trim(aer_drydep_list(m))//'DDV','m/s     ',pver, 'A', &
               trim(aer_drydep_list(m))//' deposition velocity',phys_decomp)

          if ( history_aerosol ) then 
             call add_default (trim(aer_drydep_list(m))//'DDF', 1, ' ')
             call add_default (trim(aer_drydep_list(m))//'TBF', 1, ' ')
             call add_default (trim(aer_drydep_list(m))//'GVF', 1, ' ')
          endif
       endif

    enddo count_dry_species
#endif

    allocate( mz_aerosol_ids(num_mz_aerosols) )

    do m = 1, num_mz_aerosols
       call cnst_get_ind ( aer_wetdep_list(m), mz_aerosol_ids(m), abort=.false. )
    enddo

    id_so2  => spc_ids(1)
    id_so4  => spc_ids(2)
    id_o3   => spc_ids(3)
    id_h2o2 => spc_ids(4)
    id_oh   => spc_ids(5)
    id_no3  => spc_ids(6)
    id_ho2  => spc_ids(7)
    id_dms  => spc_ids(8)

    inv_o3   = get_inv_ndx('O3') > 0
    inv_oh   = get_inv_ndx('OH') > 0
    inv_no3  = get_inv_ndx('NO3') > 0
    inv_ho2  = get_inv_ndx('HO2') > 0
    call cnst_get_ind ( 'SO2',  id_so2,  abort=.false. )
    call cnst_get_ind ( 'SO4',  id_so4,  abort=.false. )
    call cnst_get_ind ( 'H2O2', id_h2o2, abort=.false. )
    call cnst_get_ind ( 'DMS',  id_dms,  abort=.false. )

    if (inv_o3) then
       id_o3 = get_inv_ndx('O3')
    else
       call cnst_get_ind ( 'O3',   id_o3,   abort=.false. )
    endif
    if (inv_oh) then
       id_oh = get_inv_ndx('OH')
    else
       call cnst_get_ind ( 'OH',   id_oh,   abort=.false. )
    endif
    if (inv_no3) then
       id_no3 = get_inv_ndx('NO3')
    else
       call cnst_get_ind ( 'NO3',  id_no3,  abort=.false. )
    endif
    if (inv_ho2) then
       id_ho2 = get_inv_ndx('HO2')
    else
       call cnst_get_ind ( 'HO2',  id_ho2,  abort=.false. )
    endif

    do_cam_sulfchem = use_cam_sulfchem .and. all(spc_ids > 0)
    if ( do_cam_sulfchem ) then
       call inisulchem()
    endif

    if (masterproc) then
       write(iulog,*) 'mz_aero_initialize: do_cam_sulfchem = ',do_cam_sulfchem
    endif

  end subroutine mz_aero_initialize

  !===============================================================================
  subroutine mz_aero_wet_intr ( state, ptend, nstep, dt, cme, prain,    &
       evapr, cldv, cldvcu, cldvst, cldc, cldn, fracis, calday, cmfdqr, evapc, conicw, rainmr  &
#if ( defined MODAL_AERO )
     , rate1ord_cw2pr_st                                                & ! rce 2010/05/01
     , dgncur_awet,  qaerwat                                       &
#endif
     , cam_out , dlf, pbuf                                         )
    !----------------------------------------------------------------------- 
    !-----------------------------------------------------------------------
    use cam_history,   only: outfld
    use physics_types, only: physics_state, physics_ptend
    use camsrfexch,    only: cam_out_t     
    use wetdep,        only: wetdepa_v1, wetdepa_v2
    use sulchem,       only: sulf_chemdr, cldychmini, dicor
    use physconst,     only: gravit
    use constituents,  only: cnst_mw
    use physconst,     only: mwdry    ! molecular weight dry air ~ kg/kmole
    use physconst,     only: boltz    ! J/K/molecule
    use dust_intr,        only: dust_names
    use progseasalts_intr,only: progseasalts_names
    use tracer_cnst,   only: get_cnst_data
    use chem_mods,     only: fix_mass
#ifdef MODAL_AERO
    use modal_aero_data
    use modal_aero_deposition, only: set_srf_wetdep
#endif
    use physics_buffer, only : physics_buffer_desc

    !-----------------------------------------------------------------------
    implicit none
    !-----------------------------------------------------------------------
    !
    ! Arguments:
    !
    real(r8),            intent(in)  :: dt             ! time step
    type(physics_state), intent(in ) :: state          ! Physics state variables
    integer,  intent(in) :: nstep
    real(r8), intent(in) :: cme(pcols,pver)            ! local condensation of cloud water
    real(r8), intent(in) :: prain(pcols,pver)          ! production of rain
    real(r8), intent(in) :: evapr(pcols,pver)          ! evaporation of rain
    real(r8), intent(in) :: cldn(pcols,pver)           ! cloud fraction
    real(r8), intent(in) :: cldc(pcols,pver)           ! convective cloud fraction
    real(r8), intent(in) :: cldv(pcols,pver)           ! cloudy volume undergoing scavenging
    real(r8), intent(in) :: cldvcu(pcols,pver)         ! Convective precipitation area at the top interface of current layer
    real(r8), intent(in) :: cldvst(pcols,pver)         ! Stratiform precipitation area at the top interface of current layer
    real(r8), intent(in) :: conicw(pcols, pver)
    real(r8), intent(in) :: cmfdqr(pcols, pver)
    real(r8), intent(in) :: evapc(pcols, pver)         ! Evaporation rate of convective precipitation  
    real(r8), intent(in) :: rainmr(pcols, pver)         ! rain mixing ratio
    real(r8), intent(in) :: dlf(pcols,pver)            ! Detrainment of convective condensate
    type(physics_ptend), intent(inout) :: ptend        ! indivdual parameterization tendencies
    real(r8),            intent(inout) :: fracis(pcols,pver,pcnst) ! fraction of transported species that are insoluble
    real(r8), intent(in) :: calday        ! current calendar day
#if ( defined MODAL_AERO )
    real(r8), intent(in) :: rate1ord_cw2pr_st(pcols,pver) ! 1st order rate for strat cw to precip (1/s) ! rce 2010/05/01
    real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode) ! wet/ambient geom. mean diameter (m)
                                                               ! for number distribution
    real(r8), intent(inout) :: qaerwat(pcols,pver,ntot_amode)
#endif
    type(cam_out_t), intent(inout) :: cam_out
    type(physics_buffer_desc), pointer :: pbuf(:)


    !
    ! Local variables
    !
    integer  :: m                                  ! tracer index
    integer  :: ixcldice, ixcldliq
    integer  :: lchnk                              ! chunk identifier
    integer  :: ncol                               ! number of atmospheric columns
    real(r8) :: obuf(1)
    real(r8) :: iscavt(pcols, pver)
    real(r8) :: totcond(pcols, pver) ! sum of cldice and cldliq
    integer  :: yr, mon, day, ncsec
    integer  :: ncdate
    integer  :: mm
    integer  :: nphob
    real(r8), dimension(pcols,pver) :: h2o23d,o3,oh,no3,ho2,so2,so4,dms,h2o2,ekso2,ekh2o2

    real(r8), dimension(pcols,pver) ::&
         ph,          &! ph of drops
         hion          ! Hydrogen ion concentration in drops

    integer ::&
         indcp(pcols,pver),    &! Indices of cloudy grid points
         ncldypts(pver)         ! Number of cloudy grid points

    real(r8) :: icwmr2(pcols,pver) ! in cloud water mixing ratio for hack scheme
    real(r8) :: icscavt(pcols, pver)
    real(r8) :: isscavt(pcols, pver)
    real(r8) :: bcscavt(pcols, pver)
    real(r8) :: bsscavt(pcols, pver)
    real(r8) :: sol_factb, sol_facti
    real(r8) :: sol_factic(pcols,pver)
    real(r8) :: sol_factbi, sol_factii, sol_factiic
    real(r8) :: xhnm(pcols, pver)
    real(r8) :: sflx(pcols)            ! deposition flux
    real(r8) :: dicorfac(pcols)        ! factors to increase HO2 diurnal averages
    integer :: i,k
    real(r8) :: scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm) (0.1)
#ifdef MODAL_AERO
    integer :: jnv                     ! index for scavcoefnv 3rd dimension
    integer :: lphase                  ! index for interstitial / cloudborne aerosol
    integer :: lspec                   ! index for aerosol number / chem-mass / water-mass
    integer :: lcoardust, lcoarnacl    ! indices for coarse mode dust and seasalt masses
    real(r8) :: dqdt_tmp(pcols,pver)   ! temporary array to hold tendency for 1 species
    real(r8) :: f_act_conv(pcols,pver) ! prescribed aerosol activation fraction for convective cloud   ! rce 2010/05/01
    real(r8) :: f_act_conv_coarse(pcols,pver) ! similar but for coarse mode                            ! rce 2010/05/02
    real(r8) :: f_act_conv_coarse_dust, f_act_conv_coarse_nacl                                         ! rce 2010/05/02
    real(r8) :: fracis_cw(pcols,pver)
    real(r8) :: hygro_sum_old(pcols,pver)  ! before removal    [sum of (mass*hydro/dens)]
    real(r8) :: hygro_sum_del(pcols,pver)  ! removal change to [sum of (mass*hydro/dens)]
    real(r8) :: hygro_sum_old_ik, hygro_sum_new_ik
    real(r8) :: prec(pcols)                ! precipitation rate
    real(r8) :: q_tmp(pcols,pver)          ! temporary array to hold "most current" mixing ratio for 1 species
    real(r8) :: qqcw_tmp(pcols,pver)       ! temporary array to hold qqcw   ! rce 2010/05/01
    real(r8) :: scavcoefnv(pcols,pver,0:2) ! Dana and Hales coefficient (/mm) for
                                           ! cloud-borne num & vol (0), 
                                           ! interstitial num (1), interstitial vol (2)
    real(r8) :: tmpa, tmpb
    real(r8) :: tmpdust, tmpnacl
    real(r8) :: water_old, water_new   ! temporary old/new aerosol water mix-rat
    logical :: isprx(pcols,pver) ! true if precipation
    real(r8) :: aerdepwetis(pcols,pcnst)  ! aerosol wet deposition (interstitial)
    real(r8) :: aerdepwetcw(pcols,pcnst)  ! aerosol wet deposition (cloud water)
    real(r8), pointer :: fldcw(:,:)
#endif

    !-----------------------------------------------------------------------
    call cnst_get_ind('CLDICE', ixcldice)
    call cnst_get_ind('CLDLIQ', ixcldliq)
    lchnk = state%lchnk
    ncol  = state%ncol

    totcond(:ncol, :) = state%q(:ncol,:,ixcldliq) + state%q(:ncol,:,ixcldice)

    if ( num_mz_aerosols > 0 ) then

       ! Wet deposition of mozart aerosol species.
       ptend%name  = ptend%name//'+mz_aero_wetdep'

#ifdef MODAL_AERO
     prec(:ncol)=0._r8
     do k=1,pver
        where (prec(:ncol) >= 1.e-7_r8)
            isprx(:ncol,k) = .true.
        elsewhere
            isprx(:ncol,k) = .false.
        endwhere
        prec(:ncol) = prec(:ncol) + (prain(:ncol,k) + cmfdqr(:ncol,k) - evapr(:ncol,k)) &
                    *state%pdel(:ncol,k)/gravit
     end do


! calculate the mass-weighted sol_factic for coarse mode species
!    sol_factic_coarse(:,:) = 0.30_r8   ! tuned 1/4
     f_act_conv_coarse(:,:) = 0.60_r8   ! rce 2010/05/02
     f_act_conv_coarse_dust = 0.40_r8   ! rce 2010/05/02
     f_act_conv_coarse_nacl = 0.80_r8   ! rce 2010/05/02
     if (modeptr_coarse > 0) then
        lcoardust = lptr_dust_a_amode(modeptr_coarse)
        lcoarnacl = lptr_nacl_a_amode(modeptr_coarse)
        if ((lcoardust > 0) .and. (lcoarnacl > 0)) then
           do k = 1, pver
           do i = 1, ncol
              tmpdust = max( 0.0_r8, state%q(i,k,lcoardust) + ptend%q(i,k,lcoardust)*dt )
              tmpnacl = max( 0.0_r8, state%q(i,k,lcoarnacl) + ptend%q(i,k,lcoarnacl)*dt )
              if ((tmpdust+tmpnacl) > 1.0e-30_r8) then
!                sol_factic_coarse(i,k) = (0.2_r8*tmpdust + 0.4_r8*tmpnacl)/(tmpdust+tmpnacl)  ! tuned 1/6
                 f_act_conv_coarse(i,k) = (f_act_conv_coarse_dust*tmpdust &
                                         + f_act_conv_coarse_nacl*tmpnacl)/(tmpdust+tmpnacl)  ! rce 2010/05/02
              end if
           end do
           end do
        end if
     end if


    scavcoefnv(:,:,0) = 0.0_r8   ! below-cloud scavcoef = 0.0 for cloud-borne species

    do m = 1, ntot_amode   ! main loop over aerosol modes

       do lphase = 1, 2   ! loop over interstitial (1) and cloud-borne (2) forms

! sol_factb and sol_facti values
!    sol_factb - currently this is basically a tuning factor
!    sol_facti & sol_factic - currently has a physical basis, and reflects activation fraction
!
! 2008-mar-07 rce - sol_factb  (interstitial) changed from 0.3 to 0.1
!                 - sol_factic (interstitial, dust modes) changed from 1.0 to 0.5
!                 - sol_factic (cloud-borne, pcarb modes) no need to set it to 0.0
!                       because the cloud-borne pcarbon == 0 (no activation)
!
! rce 2010/05/02
! prior to this date, sol_factic was used for convective in-cloud wet removal,
!    and its value reflected a combination of an activation fraction (which varied between modes)
!    and a tuning factor
! from this date forward, two parameters are used for convective in-cloud wet removal
!    f_act_conv is the activation fraction
!       note that "non-activation" of aerosol in air entrained into updrafts should
!          be included here
!       eventually we might use the activate routine (with w ~= 1 m/s) to calculate 
!          this, but there is still the entrainment issue
!    sol_factic is strictly a tuning factor
!
          if (lphase == 1) then   ! interstial aerosol
             hygro_sum_old(:,:) = 0.0_r8
             hygro_sum_del(:,:) = 0.0_r8
             call modal_aero_bcscavcoef_get( m, ncol, isprx, dgncur_awet,   &
                                             scavcoefnv(:,:,1), scavcoefnv(:,:,2) )

             sol_factb  = 0.1_r8   ! all below-cloud scav ON (0.1 "tuning factor")
!            sol_factb  = 0.03_r8   ! all below-cloud scav ON (0.1 "tuning factor")  ! tuned 1/6

             sol_facti  = 0.0_r8   ! strat  in-cloud scav totally OFF for institial

             sol_factic = 0.4_r8      ! xl 2010/05/20

             if (m == modeptr_pcarbon) then
!               sol_factic = 0.0_r8   ! conv   in-cloud scav OFF (0.0 activation fraction)
                f_act_conv = 0.0_r8   ! rce 2010/05/02
             else if ((m == modeptr_finedust) .or. (m == modeptr_coardust)) then
!               sol_factic = 0.2_r8   ! conv   in-cloud scav ON  (0.5 activation fraction)  ! tuned 1/4
                f_act_conv = 0.4_r8   ! rce 2010/05/02
             else
!               sol_factic = 0.4_r8   ! conv   in-cloud scav ON  (1.0 activation fraction)  ! tuned 1/4
                f_act_conv = 0.8_r8   ! rce 2010/05/02
             end if

          else   ! cloud-borne aerosol (borne by stratiform cloud drops)
             sol_factb  = 0.0_r8   ! all below-cloud scav OFF (anything cloud-borne is located "in-cloud")
             sol_facti  = sol_facti_cloud_borne   ! strat  in-cloud scav cloud-borne tuning factor
             sol_factic = 0.0_r8   ! conv   in-cloud scav OFF (having this on would mean
                                   !        that conv precip collects strat droplets)
             f_act_conv = 0.0_r8   ! conv   in-cloud scav OFF (having this on would mean

          end if
!
! rce 2010/05/03
! wetdepa has 6 "sol_fact" parameters: 
!      sol_facti, sol_factic, sol_factb for liquid cloud
!      sol_factii, sol_factiic, sol_factbi for ice cloud
! the ice cloud parameters are optional, and if not provided, they default to
!    one of the other sol_fact parameters (see subr. wetdepa about this)
! for now, we set the ice cloud parameters equal 
!    to their liquid cloud counterparts
! currently the ice parameters are not used in wetdepa as
!    wetdepa sets "weight" (the ice cloud fraction) to 0.0
! if this changes, we will have to give more thought to 
!    the ice cloud parameter values
!
          sol_factbi  = sol_factb
          sol_factii  = sol_facti
          sol_factiic = sol_factic(1,1)


          do lspec = 0, nspec_amode(m)+1   ! loop over number + chem constituents + water

             if (lspec == 0) then   ! number
                if (lphase == 1) then
                   mm = numptr_amode(m)
                   jnv = 1
                else
                   mm = numptrcw_amode(m)
                   jnv = 0
                endif
             else if (lspec <= nspec_amode(m)) then   ! non-water mass
                if (lphase == 1) then
                   mm = lmassptr_amode(lspec,m)
                   jnv = 2
                else
                   mm = lmassptrcw_amode(lspec,m)
                   jnv = 0
                endif
             else   ! water mass
!   bypass wet removal of aerosol water
                cycle
                if (lphase == 1) then
                   mm = 0
!                  mm = lwaterptr_amode(m)
                   jnv = 2
                else
                   mm = 0
                   jnv = 0
                endif
             endif

          if (mm <= 0) cycle


! set f_act_conv for interstitial (lphase=1) coarse mode species
! for the convective in-cloud, we conceptually treat the coarse dust and seasalt
!    as being externally mixed, and apply f_act_conv = f_act_conv_coarse_dust/nacl to dust/seasalt
! number and sulfate are conceptually partitioned to the dust and seasalt
!    on a mass basis, so the f_act_conv for number and sulfate are
!    mass-weighted averages of the values used for dust/seasalt
          if ((lphase == 1) .and. (m == modeptr_coarse)) then
!            sol_factic = sol_factic_coarse
             f_act_conv = f_act_conv_coarse   ! rce 2010/05/02
             if (lspec > 0) then
                if (lmassptr_amode(lspec,m) == lptr_dust_a_amode(m)) then
!                  sol_factic = 0.2_r8   ! tuned 1/4
                   f_act_conv = f_act_conv_coarse_dust   ! rce 2010/05/02
                else if (lmassptr_amode(lspec,m) == lptr_nacl_a_amode(m)) then
!                  sol_factic = 0.4_r8   ! tuned 1/6
                   f_act_conv = f_act_conv_coarse_nacl   ! rce 2010/05/02
                end if
             end if
          end if
          if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
             ptend%lq(mm) = .TRUE.
             dqdt_tmp(:,:) = 0.0_r8
! q_tmp reflects changes from modal_aero_calcsize and is the "most current" q
             q_tmp(1:ncol,:) = state%q(1:ncol,:,mm) + ptend%q(1:ncol,:,mm)*dt
             fldcw => qqcw_get_field(pbuf, mm,lchnk)
             call wetdepa_v2( state%t, state%pmid, state%q, state%pdel,  &
                  cldn, cldc, cmfdqr, evapc, conicw, prain, cme,                     &
                  evapr, totcond, q_tmp(:,:), dt,            &
                  dqdt_tmp(:,:), iscavt, cldv, cldvcu, cldvst, dlf, fracis(:,:,mm), sol_factb, ncol, &
                  scavcoefnv(:,:,jnv), &
                  .false., rate1ord_cw2pr_st, fldcw, f_act_conv, &     ! rce 2010/05/01
                  icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
                  sol_facti_in=sol_facti, sol_factbi_in=sol_factbi, sol_factii_in=sol_factii, &   ! rce 2010/05/03
                  sol_factic_in=sol_factic, sol_factiic_in=sol_factiic )                          ! rce 2010/05/03

             ptend%q(1:ncol,:,mm) = ptend%q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:)

             call outfld( trim(cnst_name(mm))//'WET', dqdt_tmp(:,:), pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SIC', icscavt, pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)

             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)
             aerdepwetis(:ncol,mm) = sflx(:ncol)

             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name(mm))//'SFSIC', sflx, pcols, lchnk)
             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name(mm))//'SFSIS', sflx, pcols, lchnk)
             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name(mm))//'SFSBC', sflx, pcols, lchnk)
             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name(mm))//'SFSBS', sflx, pcols, lchnk)

             if (lspec > 0) then
!                tmpa = spechygro(lspectype_amode(lspec,m))/ &
!                       specdens_amode(lspectype_amode(lspec,m))
                tmpa = spechygro(lspec,m)/ &
                       specdens_amode(lspec,m)
                tmpb = tmpa*dt
                hygro_sum_old(1:ncol,:) = hygro_sum_old(1:ncol,:) &
                                         + tmpa*q_tmp(1:ncol,:)
                hygro_sum_del(1:ncol,:) = hygro_sum_del(1:ncol,:) &
                                         + tmpb*dqdt_tmp(1:ncol,:)
             end if

          else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then
! aerosol water -- because of how wetdepa treats evaporation of stratiform
!    precip, it is not appropriate to apply wetdepa to aerosol water
! instead, "hygro_sum" = [sum of (mass*hygro/dens)] is calculated before and 
!    after wet removal, and new water is calculated using
!    new_water = old_water*min(10,(hygro_sum_new/hygro_sum_old))
! the "min(10,...)" is to avoid potential problems when hygro_sum_old ~= 0
! also, individual wet removal terms (ic,is,bc,bs) are not output to history
!            ptend%lq(mm) = .TRUE.
!            dqdt_tmp(:,:) = 0.0_r8
             do k = top_lev, pver
             do i = 1, ncol
!               water_old = max( 0.0_r8, state%q(i,k,mm)+ptend%q(i,k,mm)*dt )
                water_old = max( 0.0_r8, qaerwat(i,k,mm) )
                hygro_sum_old_ik = max( 0.0_r8, hygro_sum_old(i,k) )
                hygro_sum_new_ik = max( 0.0_r8, hygro_sum_old_ik+hygro_sum_del(i,k) )
                if (hygro_sum_new_ik >= 10.0_r8*hygro_sum_old_ik) then
                   water_new = 10.0_r8*water_old
                else
                   water_new = water_old*(hygro_sum_new_ik/hygro_sum_old_ik)
                end if
!               dqdt_tmp(i,k) = (water_new - water_old)/dt
                qaerwat(i,k,mm) = water_new
             end do
             end do

!            ptend%q(1:ncol,:,mm) = ptend%q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:)

!            call outfld( trim(cnst_name(mm))//'WET', dqdt_tmp(:,:), pcols, lchnk)

!            sflx(:)=0._r8
!            do k=1,pver
!               do i=1,ncol
!                  sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
!               enddo
!            enddo
!            call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)

          else   ! lphase == 2
             dqdt_tmp(:,:) = 0.0_r8
             qqcw_tmp(:,:) = 0.0_r8    ! rce 2010/05/01
             fldcw => qqcw_get_field(pbuf, mm,lchnk)
             call wetdepa_v2( state%t, state%pmid, state%q, state%pdel,  &
                  cldn, cldc, cmfdqr, evapc, conicw, prain, cme,                     &
                  evapr, totcond, fldcw, dt,            &
                  dqdt_tmp(:,:), iscavt, cldv, cldvcu, cldvst, dlf, fracis_cw(:,:), sol_factb, ncol, &
                  scavcoefnv(:,:,jnv), &
                  .true., rate1ord_cw2pr_st, qqcw_tmp, f_act_conv, &     ! rce 2010/05/01
                  icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
                  sol_facti_in=sol_facti, sol_factbi_in=sol_factbi, sol_factii_in=sol_factii, &   ! rce 2010/05/03
                  sol_factic_in=sol_factic, sol_factiic_in=sol_factiic )                          ! rce 2010/05/03

             fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt

             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name_cw(mm))//'SFWET', sflx, pcols, lchnk)
             aerdepwetcw(:ncol,mm) = sflx(:ncol)

             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name_cw(mm))//'SFSIC', sflx, pcols, lchnk)
             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name_cw(mm))//'SFSIS', sflx, pcols, lchnk)
             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name_cw(mm))//'SFSBC', sflx, pcols, lchnk)
             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name_cw(mm))//'SFSBS', sflx, pcols, lchnk)

          endif

          enddo   ! lspec = 0, nspec_amode(m)+1
       enddo   ! lphase = 1, 2
    enddo   ! m = 1, ntot_amode

    ! if the user has specified prescribed aerosol dep fluxes then 
    ! do not set cam_out dep fluxes according to the prognostic aerosols
    if (.not.aerodep_flx_prescribed()) then
       call set_srf_wetdep(aerdepwetis, aerdepwetcw, cam_out)
    endif

#else
       scavcoef(:ncol,:)=0.1_r8
       do m = 1, num_mz_aerosols

          mm = mz_aerosol_ids(m)

          if (mm > 0 .and. (.not. any(dust_names==cnst_name(mm)))  .and. (.not. any(progseasalts_names==cnst_name(mm))) ) then
             ptend%lq(mm) = .TRUE.

             sol_factb = 0.3_r8
             sol_facti = 0.3_r8
             call wetdepa_v1( state%t, state%pmid, state%q, state%pdel,  &
                  cldn, cldc, cmfdqr, conicw, prain, cme,                     &
                  evapr, totcond, state%q(:,:,mm), dt,            &
                  ptend%q(:,:,mm), iscavt, cldv, fracis(:,:,mm), sol_factb, ncol, &
                  scavcoef, icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
                  sol_facti_in=sol_facti)

             call outfld( trim(cnst_name(mm))//'WET', ptend%q(:,:,mm), pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SIC', icscavt , pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)

             sflx(:)=0._r8
             do k=1,pver
                do i=1,ncol
                   sflx(i)=sflx(i)+ptend%q(i,k,mm)*state%pdel(i,k)/gravit
                enddo
             enddo
             call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)

             ! if the user has specified prescribed aerosol dep fluxes then 
             ! do not set cam_out dep fluxes according to the prognostic aerosols
             if (.not.aerodep_flx_prescribed()) then
                ! export deposition fluxes to coupler
                select case (trim(cnst_name(mm)))
                case('CB2')
                   do i = 1, ncol
                      cam_out%bcphiwet(i) = max(-sflx(i), 0._r8)
                   end do
                case('OC2')
                   do i = 1, ncol
                      cam_out%ocphiwet(i) = max(-sflx(i), 0._r8)
                   end do
                end select
             endif

          endif

       end do
#endif
    endif

    if ( do_cam_sulfchem ) then

       so4(:ncol,:)  = state%q(:ncol,:,id_so4) + ptend%q(:ncol,:,id_so4)*dt
       so2(:ncol,:)  = state%q(:ncol,:,id_so2) + ptend%q(:ncol,:,id_so2)*dt
       dms(:ncol,:)  = state%q(:ncol,:,id_dms)
       h2o2(:ncol,:) = state%q(:ncol,:,id_h2o2)

       ! need oh,ho2,no3 in molec/cm3, o3 in kg/kg

       if ( inv_o3 ) then
         call get_cnst_data( 'O3', o3(:ncol,:),  ncol, lchnk, pbuf ) ! vmr
       else
         o3(:ncol,:) = state%q(:ncol,:,id_o3)*mwdry/fix_mass(id_o3) ! vmr
       endif

       ! vmr -> molec/cm3
       xhnm(:ncol,:) = 1.e-6_r8 * state%pmiddry(:ncol,:) / (boltz*state%t(:ncol,:))

       if ( inv_oh ) then
         call get_cnst_data( 'OH', oh(:ncol,:),  ncol, lchnk, pbuf ) ! vmr
         oh(:ncol,:) = oh(:ncol,:) * xhnm(:ncol,:)             ! molec/cm3
       else
         oh(:ncol,:) = state%q(:ncol,:,id_oh)* (mwdry/cnst_mw(id_oh )) * xhnm(:ncol,:)
       endif
       if ( inv_no3 ) then
         call get_cnst_data( 'NO3', no3(:ncol,:),  ncol, lchnk, pbuf ) ! vmr
         no3(:ncol,:) = no3(:ncol,:) * xhnm(:ncol,:)             ! molec/cm3
       else
         no3(:ncol,:) = state%q(:ncol,:,id_no3)* (mwdry/cnst_mw(id_no3)) * xhnm(:ncol,:)
       endif
       if ( inv_ho2 ) then
         call get_cnst_data( 'HO2', ho2(:ncol,:),  ncol, lchnk, pbuf ) ! vmr
         ho2(:ncol,:) = ho2(:ncol,:) * xhnm(:ncol,:)             ! molec/cm3
       else
         ho2(:ncol,:) = state%q(:ncol,:,id_ho2)* (mwdry/cnst_mw(id_ho2)) * xhnm(:ncol,:) 
       endif

       ! alter HO2 from diurnal average to instantaneous value
       call dicor( calday, state%lat, ncol, dicorfac)
       do k = 1, pver
          do i = 1, ncol
             ho2(i,k) = ho2(i,k)*dicorfac(i)
          end do
       end do

       icwmr2(:,:) = 0._r8
       call cldychmini( state%t,  state%pmid, totcond, rainmr, cldn, &
            cldv, conicw, icwmr2, so2, so4, &
            ph, hion, ekso2, ekh2o2, indcp, &
            ncldypts, ncol )
       call outfld( 'PH', ph, pcols, lchnk )

       call sulf_chemdr( state%t, state%pmid, totcond, rainmr, cldn, cldv, &
                         conicw, icwmr2, indcp, ncldypts, hion, &
                         so2, so4, dms,  h2o2, ekso2, ekh2o2, &
                         o3, h2o23d, oh, no3, ho2, state%q(:,:,1),  &
                         state%zm, dt, calday, state%lat, state%lon, ncol, lchnk  )

       ptend%q(:ncol,:,id_so2) = (so2(:ncol,:)-state%q(:ncol,:,id_so2))/dt
       ptend%lq(id_so2) = .true.

       ptend%q(:ncol,:,id_so4) = (so4(:ncol,:)-state%q(:ncol,:,id_so4))/dt
       ptend%lq(id_so4) = .true.

       ptend%q(:ncol,:,id_dms) = (dms(:ncol,:)-state%q(:ncol,:,id_dms))/dt
       ptend%lq(id_dms) = .true.

       ptend%q(:ncol,:,id_h2o2) = (h2o2(:ncol,:)-state%q(:ncol,:,id_h2o2))/dt
       ptend%lq(id_h2o2) = .true.

    endif
    return

  end subroutine mz_aero_wet_intr


#ifdef MODAL_AERO
  !===============================================================================
  subroutine mz_aero_dry_intr ( state, ptend, wvflx, dt, lat, clat, &
       fsds, obklen, ts, ustar, prect, snowh, pblh, hflx, month, landfrac, &
       icefrac, ocnfrac, fvin,ram1in, dgncur_awet, wetdens,  qaerwat, &
       cam_out, pbuf)
    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! Interface to dry deposition and sedimentation of aerosols
    !
    ! Author: Xiaohong Liu and Dick Easter
    !
    !-----------------------------------------------------------------------
    use cam_history,       only: outfld
    use ppgrid,            only: pverp
    use physics_types,     only: physics_state, physics_ptend
    use camsrfexch,        only: cam_out_t
    use physconst,         only: gravit, rair, rhoh2o
    use drydep_mod,        only: setdvel,  d3ddflux, calcram
    use dust_sediment_mod, only: dust_sediment_tend, dust_sediment_vel
    use modal_aero_data
    use modal_aero_deposition, only: set_srf_drydep
    use physics_buffer, only : physics_buffer_desc

    !-----------------------------------------------------------------------
    implicit none
    !-----------------------------------------------------------------------
    !
    ! Arguments:
    !
    real(r8),            intent(in)  :: dt             ! time step
    type(physics_state), intent(in ) :: state          ! Physics state variables
    integer, intent(in) :: lat(pcols)                  ! latitude index for S->N storage
    real(r8), intent(in) :: clat(pcols)                 ! latitude
    real(r8), intent(in) :: fsds(pcols)                 ! longwave down at sfc
    real(r8), intent(in) :: obklen(pcols)                 ! obklen
    real(r8), intent(in) :: ustar(pcols)                  ! sfc fric vel--used over oceans and sea ice.
    real(r8), intent(in) :: ts(pcols)                     ! sfc temp
    real(r8), intent(in) :: landfrac(pcols)               ! land fraction
    real(r8), intent(in) :: icefrac(pcols)                ! ice fraction
    real(r8), intent(in) :: ocnfrac(pcols)                ! ocean fraction
    real(r8), intent(in) :: hflx(pcols)                  ! sensible heat flux
    real(r8), intent(in) :: prect(pcols)                     ! prect
    real(r8), intent(in) :: snowh(pcols)                     ! snow depth
    real(r8), intent(in) :: pblh(pcols)                     ! pbl height
    integer, intent(in)  :: month
    real(r8), intent(in) :: wvflx(pcols)       ! water vapor flux
    real(r8), intent(in) :: fvin(pcols)        ! for dry dep velocities from land model
    real(r8), intent(in) :: ram1in(pcols)       ! for dry dep velocities from land model
    real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
    real(r8), intent(in) :: wetdens(pcols,pver,ntot_amode)
    real(r8), intent(inout) :: qaerwat(pcols,pver,ntot_amode)

    type(physics_ptend), intent(inout) :: ptend          ! indivdual parameterization tendencies
    type(cam_out_t),     intent(inout) :: cam_out
    type(physics_buffer_desc), pointer :: pbuf(:)

    !
    ! Local variables
    !
    integer :: i,k
    integer :: jvlc                    ! index for last dimension of vlc_xxx arrays
    integer :: lphase                  ! index for interstitial / cloudborne aerosol
    integer :: lspec                   ! index for aerosol number / chem-mass / water-mass
    integer :: m                       ! aerosol mode index
    integer :: mm                      ! tracer index
    integer :: lchnk                   ! chunk identifier
    integer :: ncol                    ! number of atmospheric columns
    real(r8) :: tvs(pcols,pver)
    real(r8) :: sflx(pcols)            ! deposition flux
    real(r8)::  dep_trb(pcols)       !kg/m2/s
    real(r8)::  dep_dry(pcols)       !kg/m2/s (total of grav and trb)
    real(r8)::  dep_grv(pcols)       !kg/m2/s (total of grav and trb)
    real (r8) :: rho(pcols,pver)                    ! air density in kg/m3
    real(r8) :: pvmzaer(pcols,pverp)    ! sedimentation velocity in Pa
    real(r8) :: fv(pcols)         ! for dry dep velocities, from land modified over ocean & ice
    real(r8) :: ram1(pcols)       ! for dry dep velocities, from land modified over ocean & ice
    real(r8) :: dqdt_tmp(pcols,pver)   ! temporary array to hold tendency for 1 species

    real(r8) :: rad_drop(pcols,pver)
    real(r8) :: dens_drop(pcols,pver)
    real(r8) :: sg_drop(pcols,pver)
    real(r8) :: rad_aer(pcols,pver)
    real(r8) :: dens_aer(pcols,pver)
    real(r8) :: sg_aer(pcols,pver)

    real(r8) :: vlc_dry(pcols,pver,4)       ! dep velocity
    real(r8) :: vlc_grv(pcols,pver,4)       ! dep velocity
    real(r8)::  vlc_trb(pcols,4)            ! dep velocity
    real(r8) :: aerdepdryis(pcols,pcnst)  ! aerosol dry deposition (interstitial)
    real(r8) :: aerdepdrycw(pcols,pcnst)  ! aerosol dry deposition (cloud water)
    real(r8), pointer :: fldcw(:,:)
    !-----------------------------------------------------------------------

    lchnk = state%lchnk
    ncol  = state%ncol
    tvs(:ncol,:) = state%t(:ncol,:)!*(1+state%q(:ncol,k)
    rho(:ncol,:)=  state%pmid(:ncol,:)/(rair*state%t(:ncol,:))

    ! Dry deposition of mozart aerosol species.
    ptend%name  = ptend%name//'+mz_aero_drydep'

    !   #################################
    !    call setdvel( ncol, landfrac, icefrac, ocnfrac, .001_r8, .001_r8, .001_r8, dvel )
    !  we get the ram1,fv from the land model as ram1in,fvin, but need to calculate it over oceans and ice.
    !  better if we got thse from the ocean and ice model
    !  for friction velocity, we use ustar (from taux and tauy), except over land, where we use fv from the land model.

    call calcram(ncol,landfrac,icefrac,ocnfrac,obklen,&
         ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
         state%pdel(:,pver),fvin,fv)
!   call outfld( 'airFV', fv(:), pcols, lchnk )
!   call outfld( 'RAM1', ram1(:), pcols, lchnk )
    !       call outfld( 'icefrc', icefrac(:), pcols, lchnk )


!
! calc settling/deposition velocities for cloud droplets (and cloud-borne aerosols)
!
! *** mean drop radius should eventually be computed from ndrop and qcldwtr
    rad_drop(:,:) = 5.0e-6_r8
    dens_drop(:,:) = rhoh2o
    sg_drop(:,:) = 1.46_r8
    jvlc = 3
    call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv,  &
                     vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
                     rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 0, lchnk)
    jvlc = 4
    call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv,  &
                     vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
                     rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 3, lchnk)



    do m = 1, ntot_amode   ! main loop over aerosol modes

       do lphase = 1, 2   ! loop over interstitial / cloud-borne forms

          if (lphase == 1) then   ! interstial aerosol - calc settling/dep velocities of mode

! rad_aer = volume mean wet radius (m)
! dgncur_awet = geometric mean wet diameter for number distribution (m)
             rad_aer(1:ncol,:top_lev-1) = 0._r8
             rad_aer(1:ncol,top_lev:) = 0.5_r8*dgncur_awet(1:ncol,top_lev:,m)   &
                                 *exp(1.5_r8*(alnsg_amode(m)**2))
! dens_aer(1:ncol,:) = wet density (kg/m3)
             dens_aer(1:ncol,:top_lev-1) = 0._r8
             dens_aer(1:ncol,top_lev:) = wetdens(1:ncol,top_lev:,m)
             sg_aer(1:ncol,:) = sigmag_amode(m)

             jvlc = 1
             call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv,  & 
                        vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
                        rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 0, lchnk)
             jvlc = 2
             call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv,  & 
                        vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc),  &
                        rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 3, lchnk)
          end if

          do lspec = 0, nspec_amode(m)+1   ! loop over number + constituents + water

             if (lspec == 0) then   ! number
                if (lphase == 1) then
                   mm = numptr_amode(m)
                   jvlc = 1
                else
                   mm = numptrcw_amode(m)
                   jvlc = 3
                endif
             else if (lspec <= nspec_amode(m)) then   ! non-water mass
                if (lphase == 1) then
                   mm = lmassptr_amode(lspec,m)
                   jvlc = 2
                else
                   mm = lmassptrcw_amode(lspec,m)
                   jvlc = 4
                endif
             else   ! water mass
!   bypass dry deposition of aerosol water
                cycle
                if (lphase == 1) then
                   mm = 0
!                  mm = lwaterptr_amode(m)
                   jvlc = 2
                else
                   mm = 0
                   jvlc = 4
                endif
             endif


          if (mm <= 0) cycle

!         if (lphase == 1) then
          if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then
             ptend%lq(mm) = .TRUE.

             ! use pvprogseasalts instead (means making the top level 0)
             pvmzaer(:ncol,1)=0._r8
             pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)

             call outfld( trim(cnst_name(mm))//'DDV', pvmzaer(:,2:pverp), pcols, lchnk )

             if(.true.) then ! use phil's method
             !      convert from meters/sec to pascals/sec
             !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
                pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit

             !      calculate the tendencies and sfc fluxes from the above velocities
                call dust_sediment_tend( &
                     ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
                     state%q(:,:,mm),  pvmzaer,  ptend%q(:,:,mm), sflx  )
             else   !use charlie's method
                call d3ddflux( ncol, vlc_dry(:,:,jvlc), state%q(:,:,mm), state%pmid, &
                               state%pdel, tvs, sflx, ptend%q(:,:,mm), dt )
             endif

             ! apportion dry deposition into turb and gravitational settling for tapes
             do i=1,ncol
                dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
                dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
             enddo

             call outfld( trim(cnst_name(mm))//'DDF', sflx, pcols, lchnk)
             call outfld( trim(cnst_name(mm))//'TBF', dep_trb, pcols, lchnk )
             call outfld( trim(cnst_name(mm))//'GVF', dep_grv, pcols, lchnk )
             call outfld( trim(cnst_name(mm))//'DTQ', ptend%q(:,:,mm), pcols, lchnk)
             aerdepdryis(:ncol,mm) = sflx(:ncol)

          else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then  ! aerosol water
             ! use pvprogseasalts instead (means making the top level 0)
             pvmzaer(:ncol,1)=0._r8
             pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)

             if(.true.) then ! use phil's method
             !      convert from meters/sec to pascals/sec
             !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
                pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit

             !      calculate the tendencies and sfc fluxes from the above velocities
                call dust_sediment_tend( &
                     ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
                     qaerwat(:,:,mm),  pvmzaer,  dqdt_tmp(:,:), sflx  )
             else   !use charlie's method
                call d3ddflux( ncol, vlc_dry(:,:,jvlc), qaerwat(:,:,mm), state%pmid, &
                               state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
             endif

             ! apportion dry deposition into turb and gravitational settling for tapes
             do i=1,ncol
                dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
                dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
             enddo

             qaerwat(1:ncol,:,mm) = qaerwat(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) * dt

          else  ! lphase == 2
             ! use pvprogseasalts instead (means making the top level 0)
             pvmzaer(:ncol,1)=0._r8
             pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc)
             fldcw => qqcw_get_field(pbuf, mm,lchnk)
             if(.true.) then ! use phil's method
             !      convert from meters/sec to pascals/sec
             !      pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion
                pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit

             !      calculate the tendencies and sfc fluxes from the above velocities
                call dust_sediment_tend( &
                     ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
                     fldcw(:,:),  pvmzaer,  dqdt_tmp(:,:), sflx  )
             else   !use charlie's method
                call d3ddflux( ncol, vlc_dry(:,:,jvlc), fldcw(:,:), state%pmid, &
                               state%pdel, tvs, sflx, dqdt_tmp(:,:), dt )
             endif

             ! apportion dry deposition into turb and gravitational settling for tapes
             do i=1,ncol
                dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc)
                dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc)
             enddo

             fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt

             call outfld( trim(cnst_name_cw(mm))//'DDF', sflx, pcols, lchnk)
             call outfld( trim(cnst_name_cw(mm))//'TBF', dep_trb, pcols, lchnk )
             call outfld( trim(cnst_name_cw(mm))//'GVF', dep_grv, pcols, lchnk )
             aerdepdrycw(:ncol,mm) = sflx(:ncol)

          endif

          enddo   ! lspec = 0, nspec_amode(m)+1
       enddo   ! lphase = 1, 2
    enddo   ! m = 1, ntot_amode

    ! if the user has specified prescribed aerosol dep fluxes then 
    ! do not set cam_out dep fluxes according to the prognostic aerosols
    if (.not.aerodep_flx_prescribed()) then
       call set_srf_drydep(aerdepdryis, aerdepdrycw, cam_out)
    endif

    return
  end subroutine mz_aero_dry_intr


  !===============================================================================
  subroutine modal_aero_depvel_part( ncol, t, pmid, ram1, fv, vlc_dry, vlc_trb, vlc_grv,  &
                                     radius_part, density_part, sig_part, moment, lchnk )

!    calculates surface deposition velocity of particles
!    L. Zhang, S. Gong, J. Padro, and L. Barrie
!    A size-seggregated particle dry deposition scheme for an atmospheric aerosol module
!    Atmospheric Environment, 35, 549-560, 2001.
!
!    Authors: X. Liu

    !
    ! !USES
    !
    use physconst,     only: pi,boltz, gravit, rair
    use mo_drydep,     only: n_land_type, fraction_landuse

    ! !ARGUMENTS:
    !
    implicit none
    !
    real(r8), intent(in) :: t(pcols,pver)       !atm temperature (K)
    real(r8), intent(in) :: pmid(pcols,pver)    !atm pressure (Pa)
    real(r8), intent(in) :: fv(pcols)           !friction velocity (m/s)
    real(r8), intent(in) :: ram1(pcols)         !aerodynamical resistance (s/m)
    real(r8), intent(in) :: radius_part(pcols,pver)    ! mean (volume/number) particle radius (m)
    real(r8), intent(in) :: density_part(pcols,pver)   ! density of particle material (kg/m3)
    real(r8), intent(in) :: sig_part(pcols,pver)       ! geometric standard deviation of particles
    integer,  intent(in) :: moment ! moment of size distribution (0 for number, 2 for surface area, 3 for volume)
    integer,  intent(in) :: ncol
    integer,  intent(in) :: lchnk

    real(r8), intent(out) :: vlc_trb(pcols)       !Turbulent deposn velocity (m/s)
    real(r8), intent(out) :: vlc_grv(pcols,pver)       !grav deposn velocity (m/s)
    real(r8), intent(out) :: vlc_dry(pcols,pver)       !dry deposn velocity (m/s)
    !------------------------------------------------------------------------

    !------------------------------------------------------------------------
    ! Local Variables
    integer  :: m,i,k,ix                !indices
    real(r8) :: rho     !atm density (kg/m**3)
    real(r8) :: vsc_dyn_atm(pcols,pver)   ![kg m-1 s-1] Dynamic viscosity of air
    real(r8) :: vsc_knm_atm(pcols,pver)   ![m2 s-1] Kinematic viscosity of atmosphere
    real(r8) :: shm_nbr       ![frc] Schmidt number
    real(r8) :: stk_nbr       ![frc] Stokes number
    real(r8) :: mfp_atm(pcols,pver)       ![m] Mean free path of air
    real(r8) :: dff_aer       ![m2 s-1] Brownian diffusivity of particle
    real(r8) :: slp_crc(pcols,pver) ![frc] Slip correction factor
    real(r8) :: rss_trb       ![s m-1] Resistance to turbulent deposition
    real(r8) :: rss_lmn       ![s m-1] Quasi-laminar layer resistance
    real(r8) :: brownian      ! collection efficiency for Browning diffusion
    real(r8) :: impaction     ! collection efficiency for impaction
    real(r8) :: interception  ! collection efficiency for interception
    real(r8) :: stickfrac     ! fraction of particles sticking to surface
    real(r8) :: radius_moment(pcols,pver) ! median radius (m) for moment
    real(r8) :: lnsig         ! ln(sig_part)
    real(r8) :: dispersion    ! accounts for influence of size dist dispersion on bulk settling velocity
                              ! assuming radius_part is number mode radius * exp(1.5 ln(sigma))

    integer  :: lt
    real(r8) :: lnd_frc
    real(r8) :: wrk1, wrk2, wrk3

    ! constants
    real(r8) gamma(11)      ! exponent of schmidt number
!   data gamma/0.54d+00,  0.56d+00,  0.57d+00,  0.54d+00,  0.54d+00, &
!              0.56d+00,  0.54d+00,  0.54d+00,  0.54d+00,  0.56d+00, &
!              0.50d+00/
    data gamma/0.56e+00_r8,  0.54e+00_r8,  0.54e+00_r8,  0.56e+00_r8,  0.56e+00_r8, &        
               0.56e+00_r8,  0.50e+00_r8,  0.54e+00_r8,  0.54e+00_r8,  0.54e+00_r8, &
               0.54e+00_r8/
    save gamma

    real(r8) alpha(11)      ! parameter for impaction
!   data alpha/50.00d+00,  0.95d+00,  0.80d+00,  1.20d+00,  1.30d+00, &
!               0.80d+00, 50.00d+00, 50.00d+00,  2.00d+00,  1.50d+00, &
!             100.00d+00/
    data alpha/1.50e+00_r8,   1.20e+00_r8,  1.20e+00_r8,  0.80e+00_r8,  1.00e+00_r8, &
               0.80e+00_r8, 100.00e+00_r8, 50.00e+00_r8,  2.00e+00_r8,  1.20e+00_r8, &
              50.00e+00_r8/
    save alpha

    real(r8) radius_collector(11) ! radius (m) of surface collectors
!   data radius_collector/-1.00d+00,  5.10d-03,  3.50d-03,  3.20d-03, 10.00d-03, &
!                          5.00d-03, -1.00d+00, -1.00d+00, 10.00d-03, 10.00d-03, &
!                         -1.00d+00/
    data radius_collector/10.00e-03_r8,  3.50e-03_r8,  3.50e-03_r8,  5.10e-03_r8,  2.00e-03_r8, &
                           5.00e-03_r8, -1.00e+00_r8, -1.00e+00_r8, 10.00e-03_r8,  3.50e-03_r8, &
                          -1.00e+00_r8/
    save radius_collector

    integer            :: iwet(11) ! flag for wet surface = 1, otherwise = -1
!   data iwet/1,   -1,   -1,   -1,   -1,  &
!            -1,   -1,   -1,    1,   -1,  &
!             1/
    data iwet/-1,  -1,   -1,   -1,   -1,  &
              -1,   1,   -1,    1,   -1,  &
              -1/
    save iwet


    !------------------------------------------------------------------------

    vlc_grv(:ncol,:top_lev-1) = 0._r8
    vlc_dry(:ncol,:top_lev-1) = 0._r8

    do k=top_lev,pver
       do i=1,ncol

          lnsig = log(sig_part(i,k))
! use a maximum radius of 50 microns when calculating deposition velocity
          radius_moment(i,k) = min(50.0e-6_r8,radius_part(i,k))*   &
                          exp((float(moment)-1.5_r8)*lnsig*lnsig)
          dispersion = exp(2._r8*lnsig*lnsig)

          rho=pmid(i,k)/rair/t(i,k)

          ! Quasi-laminar layer resistance: call rss_lmn_get
          ! Size-independent thermokinetic properties
          vsc_dyn_atm(i,k) = 1.72e-5_r8 * ((t(i,k)/273.0_r8)**1.5_r8) * 393.0_r8 / &
               (t(i,k)+120.0_r8)      ![kg m-1 s-1] RoY94 p. 102
          mfp_atm(i,k) = 2.0_r8 * vsc_dyn_atm(i,k) / &   ![m] SeP97 p. 455
               (pmid(i,k)*sqrt(8.0_r8/(pi*rair*t(i,k))))
          vsc_knm_atm(i,k) = vsc_dyn_atm(i,k) / rho ![m2 s-1] Kinematic viscosity of air

          slp_crc(i,k) = 1.0_r8 + mfp_atm(i,k) * &
                  (1.257_r8+0.4_r8*exp(-1.1_r8*radius_moment(i,k)/(mfp_atm(i,k)))) / &
                  radius_moment(i,k)   ![frc] Slip correction factor SeP97 p. 464
          vlc_grv(i,k) = (4.0_r8/18.0_r8) * radius_moment(i,k)*radius_moment(i,k)*density_part(i,k)* &
                  gravit*slp_crc(i,k) / vsc_dyn_atm(i,k) ![m s-1] Stokes' settling velocity SeP97 p. 466
          vlc_grv(i,k) = vlc_grv(i,k) * dispersion

          vlc_dry(i,k)=vlc_grv(i,k)
       enddo
    enddo
    k=pver  ! only look at bottom level for next part
    do i=1,ncol
       dff_aer = boltz * t(i,k) * slp_crc(i,k) / &    ![m2 s-1]
                 (6.0_r8*pi*vsc_dyn_atm(i,k)*radius_moment(i,k)) !SeP97 p.474
       shm_nbr = vsc_knm_atm(i,k) / dff_aer                        ![frc] SeP97 p.972

       wrk2 = 0._r8
       wrk3 = 0._r8
       do lt = 1,n_land_type
          lnd_frc = fraction_landuse(i,lt,lchnk)
          if ( lnd_frc /= 0._r8 ) then
             brownian = shm_nbr**(-gamma(lt))
             if (radius_collector(lt) > 0.0_r8) then
!       vegetated surface
                stk_nbr = vlc_grv(i,k) * fv(i) / (gravit*radius_collector(lt))
                interception = 2.0_r8*(radius_moment(i,k)/radius_collector(lt))**2.0_r8
             else
!       non-vegetated surface
                stk_nbr = vlc_grv(i,k) * fv(i) * fv(i) / (gravit*vsc_knm_atm(i,k))  ![frc] SeP97 p.965
                interception = 0.0_r8
             endif
             impaction = (stk_nbr/(alpha(lt)+stk_nbr))**2.0_r8   

             if (iwet(lt) > 0) then
                stickfrac = 1.0_r8
             else
                stickfrac = exp(-sqrt(stk_nbr))
                if (stickfrac < 1.0e-10_r8) stickfrac = 1.0e-10_r8
             endif
             rss_lmn = 1.0_r8 / (3.0_r8 * fv(i) * stickfrac * (brownian+interception+impaction))
             rss_trb = ram1(i) + rss_lmn + ram1(i)*rss_lmn*vlc_grv(i,k)

             wrk1 = 1.0_r8 / rss_trb
             wrk2 = wrk2 + lnd_frc*( wrk1 )
             wrk3 = wrk3 + lnd_frc*( wrk1 + vlc_grv(i,k) )
          endif
       enddo  ! n_land_type
       vlc_trb(i) = wrk2
       vlc_dry(i,k) = wrk3
    enddo !ncol

    return
  end subroutine modal_aero_depvel_part


  !===============================================================================
  subroutine modal_aero_bcscavcoef_get( m, ncol, isprx, dgn_awet, scavcoefnum, scavcoefvol )

    use modal_aero_data
    !-----------------------------------------------------------------------
    implicit none

	  integer,intent(in) :: m, ncol
	  logical,intent(in):: isprx(pcols,pver)
	  real(r8), intent(in) :: dgn_awet(pcols,pver,ntot_amode)
	  real(r8), intent(out) :: scavcoefnum(pcols,pver), scavcoefvol(pcols,pver)

	  integer i, k, jgrow
	  real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum


      scavcoefnum(:ncol,:top_lev-1) = 0._r8
      scavcoefvol(:ncol,:top_lev-1) = 0._r8

      do k = top_lev, pver
      do i = 1, ncol

!        do only if no precip
         if ( isprx(i,k)  ) then
!
!           interpolate table values using log of (actual-wet-size)/(base-dry-size)

	    dumdgratio = dgn_awet(i,k,m)/dgnum_amode(m)

	    if ((dumdgratio .ge. 0.99_r8) .and. (dumdgratio .le. 1.01_r8)) then
	        scavimpvol = scavimptblvol(0,m)
	        scavimpnum = scavimptblnum(0,m)
	    else
	        xgrow = log( dumdgratio ) / dlndg_nimptblgrow
	        jgrow = int( xgrow )
	        if (xgrow .lt. 0._r8) jgrow = jgrow - 1
	        if (jgrow .lt. nimptblgrow_mind) then
		    jgrow = nimptblgrow_mind
		    xgrow = jgrow
	        else
	            jgrow = min( jgrow, nimptblgrow_maxd-1 )
	        end if

	        dumfhi = xgrow - jgrow
	        dumflo = 1._r8 - dumfhi

	        scavimpvol = dumflo*scavimptblvol(jgrow,m) +   &
			  dumfhi*scavimptblvol(jgrow+1,m)
	        scavimpnum = dumflo*scavimptblnum(jgrow,m) +   &
			  dumfhi*scavimptblnum(jgrow+1,m)

	    end if

!           impaction scavenging removal amount for volume
	    scavcoefvol(i,k) = exp( scavimpvol  )
!           impaction scavenging removal amount to number
	    scavcoefnum(i,k) = exp( scavimpnum  )

!           scavcoef = impaction scav rate (1/h) for precip = 1 mm/h
!           scavcoef = impaction scav rate (1/s) for precip = pfx_inrain
!           (scavcoef/3600) = impaction scav rate (1/s) for precip = 1 mm/h
!           (pfx_inrain*3600) = in-rain-area precip rate (mm/h)
!           impactrate = (scavcoef/3600) * (pfx_inrain*3600)
          else
            scavcoefvol(i,k) = 0._r8
	    scavcoefnum(i,k) = 0._r8
          end if

      end do
      end do

        return
	end subroutine modal_aero_bcscavcoef_get


  !===============================================================================
  subroutine modal_aero_bcscavcoef_init
!-----------------------------------------------------------------------
!
! Purpose:
! Computes lookup table for aerosol impaction/interception scavenging rates
!
! Authors: R. Easter
!
!-----------------------------------------------------------------------

  use shr_kind_mod,only: r8 => shr_kind_r8
  use modal_aero_data
  use abortutils,  only: endrun

  implicit none


!   local variables
	integer nnfit_maxd
	parameter (nnfit_maxd=27)

	integer i, jgrow, jdens, jpress, jtemp, ll, mode, nnfit
        integer lunerr

	real(r8) dg0, dg0_cgs, press, &
	rhodryaero, rhowetaero, rhowetaero_cgs, rmserr, &
	scavratenum, scavratevol, sigmag,                &
	temp, wetdiaratio, wetvolratio
	real(r8) aafitnum(1), xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd)
	real(r8) aafitvol(1), xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd)


	lunerr = 6
	dlndg_nimptblgrow = log( 1.25_r8 )

	do 7900 mode = 1, ntot_amode

	sigmag = sigmag_amode(mode)

!BSINGH - "spectype" is no longer used
!	ll = lspectype_amode(1,mode)
!	rhodryaero = specdens_amode(ll)
    rhodryaero = specdens_amode(1,mode)

	do 7800 jgrow = nimptblgrow_mind, nimptblgrow_maxd

	wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
	dg0 = dgnum_amode(mode)*wetdiaratio

	wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 )
	rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio
	rhowetaero = min( rhowetaero, rhodryaero )

!
!   compute impaction scavenging rates at 1 temp-press pair and save
!
	nnfit = 0

	temp = 273.16_r8
	press = 0.75e6_r8   ! dynes/cm2
	rhowetaero = rhodryaero

	dg0_cgs = dg0*1.0e2_r8   ! m to cm
	rhowetaero_cgs = rhowetaero*1.0e-3_r8   ! kg/m3 to g/cm3
	call calc_1_impact_rate( &
     		dg0_cgs, sigmag, rhowetaero_cgs, temp, press, &
     		scavratenum, scavratevol, lunerr )

	nnfit = nnfit + 1
	if (nnfit .gt. nnfit_maxd) then
	    write(lunerr,9110)
	    call endrun()
	end if
9110	format( '*** subr. modal_aero_bcscavcoef_init -- nnfit too big' )

	xxfitnum(1,nnfit) = 1._r8
	yyfitnum(nnfit) = log( scavratenum )

	xxfitvol(1,nnfit) = 1._r8
	yyfitvol(nnfit) = log( scavratevol )

5900	continue

!
! skip mlinfit stuff because scav table no longer has dependencies on
!    air temp, air press, and particle wet density
! just load the log( scavrate--- ) values
!
!!
!!   do linear regression
!!	log(scavrate) = a1 + a2*log(wetdens)
!!
!	call mlinft( xxfitnum, yyfitnum, aafitnum, nnfit, 1, 1, rmserr )
!	call mlinft( xxfitvol, yyfitvol, aafitvol, nnfit, 1, 1, rmserr )
!
!	scavimptblnum(jgrow,mode) = aafitnum(1)
!	scavimptblvol(jgrow,mode) = aafitvol(1)

	scavimptblnum(jgrow,mode) = yyfitnum(1)
	scavimptblvol(jgrow,mode) = yyfitvol(1)



7800	continue
7900	continue

        return
        end subroutine modal_aero_bcscavcoef_init


  !===============================================================================
	subroutine calc_1_impact_rate(             &
     		dg0, sigmag, rhoaero, temp, press, &
     		scavratenum, scavratevol, lunerr )
!
!   routine computes a single impaction scavenging rate
!	for precipitation rate of 1 mm/h
!
!   dg0 = geometric mean diameter of aerosol number size distrib. (cm)
!   sigmag = geometric standard deviation of size distrib.
!   rhoaero = density of aerosol particles (g/cm^3)
!   temp = temperature (K)
!   press = pressure (dyne/cm^2)
!   scavratenum = number scavenging rate (1/h)
!   scavratevol = volume or mass scavenging rate (1/h)
!   lunerr = logical unit for error message
!
     use shr_kind_mod, only: r8 => shr_kind_r8

	implicit none

!   subr. parameters
	integer lunerr
	real(r8) dg0, sigmag, rhoaero, temp, press, scavratenum, scavratevol

!   local variables
	integer nrainsvmax
	parameter (nrainsvmax=50)
	real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
     		vfallrainsv(nrainsvmax)

	integer naerosvmax
	parameter (naerosvmax=51)
	real(r8) aaerosv(naerosvmax), &
     	ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)

	integer i, ja, jr, na, nr
	real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
     	real(r8) anumsum, avolsum, boltzmann, cair, chi
     	real(r8) d, dr, dum, dumfuchs, dx
     	real(r8) ebrown, eimpact, eintercept, etotal, freepath, gravity 
     	real(r8) pi, precip, precipmmhr, precipsum
     	real(r8) r, rainsweepout, reynolds, rhi, rhoair, rhowater, rlo, rnumsum
     	real(r8) scavsumnum, scavsumnumbb
     	real(r8) scavsumvol, scavsumvolbb
     	real(r8) schmidt, sqrtreynolds, sstar, stokes, sx              
     	real(r8) taurelax, vfall, vfallstp
     	real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair                     


	rlo = .005_r8
	rhi = .250_r8
	dr = 0.005_r8
	nr = 1 + nint( (rhi-rlo)/dr )
	if (nr .gt. nrainsvmax) then
	    write(lunerr,9110)
	    call endrun()
	end if
9110	format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' )

	precipmmhr = 1.0_r8
	precip = precipmmhr/36000._r8

	ag0 = dg0/2._r8
	sx = log( sigmag )
	xg0 = log( ag0 )
	xg3 = xg0 + 3._r8*sx*sx

	xlo = xg3 - 4._r8*sx
	xhi = xg3 + 4._r8*sx
	dx = 0.2_r8*sx

	dx = max( 0.2_r8*sx, 0.01_r8 )
	xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
	xhi = xg3 + max( 4._r8*sx, 2._r8*dx )

	na = 1 + nint( (xhi-xlo)/dx )
	if (na .gt. naerosvmax) then
	    write(lunerr,9120)
	    call endrun()
	end if
9120	format( '*** subr. calc_1_impact_rate -- na > naerosvmax' )

!   air molar density
	cair = press/(8.31436e7_r8*temp)
!   air mass density
	rhoair = 28.966_r8*cair
!   molecular freepath
	freepath = 2.8052e-10_r8/cair
!   boltzmann constant
	boltzmann = 1.3807e-16_r8
!   water density
	rhowater = 1.0_r8
!   gravity
	gravity = 980.616_r8
!   air dynamic viscosity
	airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) *    &
                                        ((temp/296.16_r8)**1.5_r8)
!   air kinemaic viscosity
	airkinvisc = airdynvisc/rhoair
!   ratio of water viscosity to air viscosity (from Slinn)
	xmuwaterair = 60.0_r8

	pi = 3.1415926536_r8

!
!   compute rain drop number concentrations
!	rrainsv = raindrop radius (cm)
!	xnumrainsv = raindrop number concentration (#/cm^3)
!		(number in the bin, not number density)
!	vfallrainsv = fall velocity (cm/s)
!
	precipsum = 0._r8
	do i = 1, nr
	    r = rlo + (i-1)*dr
	    rrainsv(i) = r
	    xnumrainsv(i) = exp( -r/2.7e-2_r8 )

	    d = 2._r8*r
	    if (d .le. 0.007_r8) then
		vfallstp = 2.88e5_r8 * d**2._r8
	    else if (d .le. 0.025_r8) then
		vfallstp = 2.8008e4_r8 * d**1.528_r8
	    else if (d .le. 0.1_r8) then
		vfallstp = 4104.9_r8 * d**1.008_r8
	    else if (d .le. 0.25_r8) then
		vfallstp = 1812.1_r8 * d**0.638_r8
	    else
		vfallstp = 1069.8_r8 * d**0.235_r8
	    end if

	    vfall = vfallstp * sqrt(1.204e-3_r8/rhoair)
	    vfallrainsv(i) = vfall
	    precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
	end do
	precipsum = precipsum*pi*1.333333_r8

	rnumsum = 0._r8
	do i = 1, nr
	    xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
	    rnumsum = rnumsum + xnumrainsv(i)
	end do

!
!   compute aerosol concentrations
!	aaerosv = particle radius (cm)
!	fnumaerosv = fraction of total number in the bin (--)
!	fvolaerosv = fraction of total volume in the bin (--)
!
	anumsum = 0._r8
	avolsum = 0._r8
	do i = 1, na
	    x = xlo + (i-1)*dx
	    a = exp( x )
	    aaerosv(i) = a
	    dum = (x - xg0)/sx
	    ynumaerosv(i) = exp( -0.5_r8*dum*dum )
	    yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a
	    anumsum = anumsum + ynumaerosv(i)
	    avolsum = avolsum + yvolaerosv(i)
	end do

	do i = 1, na
	    ynumaerosv(i) = ynumaerosv(i)/anumsum
	    yvolaerosv(i) = yvolaerosv(i)/avolsum
	end do


!
!   compute scavenging
!
	scavsumnum = 0._r8
	scavsumvol = 0._r8
!
!   outer loop for rain drop radius
!
	do 5900 jr = 1, nr

	r = rrainsv(jr)
	vfall = vfallrainsv(jr)

	reynolds = r * vfall / airkinvisc
	sqrtreynolds = sqrt( reynolds )

!
!   inner loop for aerosol particle radius
!
	scavsumnumbb = 0._r8
	scavsumvolbb = 0._r8

	do 5500 ja = 1, na

	a = aaerosv(ja)

	chi = a/r

	dum = freepath/a
	dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum)
	taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc)

	aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8
	aerodiffus = boltzmann*temp*taurelax/aeromass

	schmidt = airkinvisc/aerodiffus
	stokes = vfall*taurelax/r

	ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) /  &
     			(reynolds*schmidt)

	dum = (1._r8 + 2._r8*xmuwaterair*chi) /         &
     			(1._r8 + xmuwaterair/sqrtreynolds)
	eintercept = 4._r8*chi*(chi + dum)

	dum = log( 1._r8 + reynolds )
	sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum)
	eimpact = 0._r8
	if (stokes .gt. sstar) then
	    dum = stokes - sstar
	    eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8
	end if

	etotal = ebrown + eintercept + eimpact
	etotal = min( etotal, 1.0_r8 )

	rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall

	scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
	scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)

5500	continue

	scavsumnum = scavsumnum + scavsumnumbb
	scavsumvol = scavsumvol + scavsumvolbb
5900	continue

	scavratenum = scavsumnum*3600._r8
	scavratevol = scavsumvol*3600._r8


   
  return
  end subroutine calc_1_impact_rate

#endif 
! (defined MODAL_AERO)


!-------------------------------------------------------------------
!-------------------------------------------------------------------
  !===============================================================================

  function mz_prescribed_dust()

     ! Return true if MOZART is supplying prescribed dust.

     use dust_intr,    only: dust_names
     use tracer_cnst,  only: num_tracer_cnst, tracer_cnst_flds

     logical :: mz_prescribed_dust

     integer :: i
     !-----------------------------------------------------------------------

     mz_prescribed_dust = .false.
     do i = 1, num_tracer_cnst
        if ( trim(tracer_cnst_flds(i)) == dust_names(1) ) then
           mz_prescribed_dust = .true.
           exit
        end if
     end do

  end function mz_prescribed_dust

  !===============================================================================


end module mz_aerosols_intr