LCOV - code coverage report
Current view: | top level - chemistry/aerosol - aero_wetdep_cam.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 473 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 7 | 0.0 %

              Line data    Source code
           1             : module aero_wetdep_cam
           2             :   use shr_kind_mod,  only: r8 => shr_kind_r8
           3             :   use physics_types, only: physics_state, physics_ptend, physics_ptend_init
           4             :   use camsrfexch,    only: cam_out_t
           5             :   use physics_buffer,only: physics_buffer_desc, pbuf_get_index, pbuf_set_field, pbuf_get_field
           6             :   use constituents,  only: pcnst, cnst_name, cnst_get_ind
           7             :   use phys_control,  only: phys_getopts
           8             :   use ppgrid,        only: pcols, pver
           9             :   use physconst,     only: gravit
          10             :
          11             :   use cam_abortutils,only: endrun
          12             :   use cam_logfile,   only: iulog
          13             :   use spmd_utils,    only: masterproc
          14             :   use infnan,        only: nan, assignment(=)
          15             :
          16             :   use cam_history,   only: addfld, add_default, horiz_only, outfld
          17             :   use wetdep,        only: wetdep_init
          18             :
          19             :   use rad_constituents, only: rad_cnst_get_info
          20             :
          21             :   use aerosol_properties_mod, only: aero_name_len
          22             :   use aerosol_properties_mod, only: aerosol_properties
          23             :   use modal_aerosol_properties_mod, only: modal_aerosol_properties
          24             :
          25             :   use aerosol_state_mod, only: aerosol_state, ptr2d_t
          26             :   use modal_aerosol_state_mod, only: modal_aerosol_state
          27             :
          28             :   use aero_convproc, only: aero_convproc_readnl, aero_convproc_init, aero_convproc_intr
          29             :   use aero_convproc, only: convproc_do_evaprain_atonce
          30             :   use aero_convproc, only: deepconv_wetdep_history
          31             :
          32             :   use infnan,         only: nan, assignment(=)
          33             :   use perf_mod,       only: t_startf, t_stopf
          34             :
          35             :   implicit none
          36             :   private
          37             :
          38             :   public :: aero_wetdep_readnl
          39             :   public :: aero_wetdep_init
          40             :   public :: aero_wetdep_tend
          41             :
          42             :   real(r8), parameter :: NOTSET = -huge(1._r8)
          43             :   real(r8) :: sol_facti_cloud_borne   = NOTSET
          44             :   real(r8) :: sol_factb_interstitial  = NOTSET
          45             :   real(r8) :: sol_factic_interstitial = NOTSET
          46             :
          47             :   integer :: fracis_idx = -1
          48             :   integer :: rprddp_idx          = -1
          49             :   integer :: rprdsh_idx          = -1
          50             :   integer :: nevapr_shcu_idx     = -1
          51             :   integer :: nevapr_dpcu_idx     = -1
          52             :
          53             :   logical :: wetdep_active = .false.
          54             :   integer :: nwetdep = 0
          55             :   logical :: convproc_do_aer = .false.
          56             :   logical,allocatable :: aero_cnst_lq(:,:)
          57             :   integer,allocatable :: aero_cnst_id(:,:)
          58             :   logical, public, protected :: wetdep_lq(pcnst) ! set flags true for constituents with non-zero tendencies
          59             :
          60             :   ! variables for table lookup of aerosol impaction/interception scavenging rates
          61             :   integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
          62             :   real(r8) :: dlndg_nimptblgrow
          63             :   real(r8),allocatable :: scavimptblnum(:,:)
          64             :   real(r8),allocatable :: scavimptblvol(:,:)
          65             :
          66             :   integer :: nmodes=0
          67             :   integer :: nspec_max=0
          68             :   integer :: nele_tot            ! total number of aerosol elements
          69             :   class(aerosol_properties), pointer :: aero_props=>null()
          70             :
          71             : contains
          72             :
          73             :   !------------------------------------------------------------------------------
          74             :   !------------------------------------------------------------------------------
          75           0 :   subroutine aero_wetdep_readnl(nlfile)
          76             :     use namelist_utils,  only: find_group_name
          77             :     use spmd_utils,      only: mpicom, masterprocid, mpi_character, mpi_real8, mpi_integer, mpi_success
          78             :     use spmd_utils,      only: mpi_logical
          79             :
          80             :     character(len=*), intent(in)  :: nlfile  ! filepath for file containing namelist input
          81             :
          82             :     integer                       :: unitn, ierr
          83             :     character(len=*), parameter   :: subname = 'aero_wetdep_readnl'
          84             :
          85             :     ! ===================
          86             :     ! Namelist definition
          87             :     ! ===================
          88             :     namelist /aero_wetdep_nl/ sol_facti_cloud_borne, sol_factb_interstitial, sol_factic_interstitial
          89             :
          90             :     ! =============
          91             :     ! Read namelist
          92             :     ! =============
          93           0 :     if (masterproc) then
          94           0 :        open( newunit=unitn, file=trim(nlfile), status='old' )
          95           0 :        call find_group_name(unitn, 'aero_wetdep_nl', status=ierr)
          96           0 :        if (ierr == 0) then
          97           0 :           read(unitn, aero_wetdep_nl, iostat=ierr)
          98           0 :           if (ierr /= 0) then
          99           0 :              call endrun(subname // ':: ERROR reading namelist')
         100             :           end if
         101             :        end if
         102           0 :        close(unitn)
         103             :
         104             :        ! ============================
         105             :        ! Log namelist options
         106             :        ! ============================
         107           0 :        write(iulog,*) subname,' namelist settings: '
         108           0 :        write(iulog,*) '   sol_facti_cloud_borne  : ',sol_facti_cloud_borne
         109           0 :        write(iulog,*) '   sol_factb_interstitial : ',sol_factb_interstitial
         110           0 :        write(iulog,*) '   sol_factic_interstitial: ',sol_factic_interstitial
         111             :     end if
         112             :
         113             :     ! ============================
         114             :     ! Broadcast namelist variables
         115             :     ! ============================
         116           0 :     call mpi_bcast(sol_facti_cloud_borne, 1, mpi_real8, masterprocid, mpicom, ierr)
         117           0 :     if (ierr/=mpi_success) then
         118           0 :        call endrun(subname//': MPI_BCAST ERROR: sol_facti_cloud_borne')
         119             :     end if
         120           0 :     call mpi_bcast(sol_factb_interstitial, 1, mpi_real8, masterprocid, mpicom, ierr)
         121           0 :     if (ierr/=mpi_success) then
         122           0 :        call endrun(subname//': MPI_BCAST ERROR: sol_factb_interstitial')
         123             :     end if
         124           0 :     call mpi_bcast(sol_factic_interstitial, 1, mpi_real8, masterprocid, mpicom, ierr)
         125           0 :     if (ierr/=mpi_success) then
         126           0 :        call endrun(subname//': MPI_BCAST ERROR: sol_factic_interstitial')
         127             :     end if
         128             :
         129           0 :     call mpi_bcast(nwetdep, 1, mpi_integer, masterprocid, mpicom, ierr)
         130           0 :     if (ierr/=mpi_success) then
         131           0 :        call endrun(subname//': MPI_BCAST ERROR: nwetdep')
         132             :     end if
         133             :
         134           0 :     wetdep_active = .true. !nwetdep>0
         135             :
         136           0 :     if (masterproc) then
         137           0 :        write(iulog,*) subname,' wetdep_active = ',wetdep_active,' nwetdep = ',nwetdep
         138             :     endif
         139             :
         140           0 :     call aero_convproc_readnl(nlfile)
         141             :
         142           0 :   end subroutine aero_wetdep_readnl
         143             :
         144             :   !------------------------------------------------------------------------------
         145             :   !------------------------------------------------------------------------------
         146           0 :   subroutine aero_wetdep_init( )
         147             :
         148             :     character(len=*), parameter :: subrname = 'aero_wetdep_init'
         149             :
         150             :     character(len=2)  :: unit_basename  ! Units 'kg' or '1'
         151             :     character(len=aero_name_len) :: tmpname
         152             :     character(len=aero_name_len) :: tmpname_cw
         153             :
         154             :     logical  :: history_aerosol ! Output MAM or SECT aerosol tendencies
         155             :     logical  :: history_chemistry
         156             :
         157             :     integer :: l,m, id, astat
         158             :     character(len=2) :: binstr
         159             :
         160           0 :     fracis_idx = pbuf_get_index('FRACIS')
         161           0 :     rprddp_idx      = pbuf_get_index('RPRDDP')
         162           0 :     rprdsh_idx      = pbuf_get_index('RPRDSH')
         163           0 :     nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
         164           0 :     nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
         165             :
         166           0 :     if (.not.wetdep_active) return
         167             :
         168             :     call phys_getopts(history_aerosol_out = history_aerosol, &
         169             :                       history_chemistry_out=history_chemistry, &
         170             :                       convproc_do_aer_out = convproc_do_aer)
         171             :
         172           0 :     call rad_cnst_get_info(0, nmodes=nmodes)
         173             :
         174           0 :     if (nmodes>0) then
         175           0 :        aero_props => modal_aerosol_properties()
         176           0 :        if (.not.associated(aero_props)) then
         177           0 :           call endrun(subrname//' : construction of aero_props modal_aerosol_properties object failed')
         178             :        end if
         179             :     else
         180           0 :        call endrun(subrname//' : cannot determine aerosol model')
         181             :     endif
         182             :
         183           0 :     nele_tot = aero_props%ncnst_tot()
         184             :
         185           0 :     allocate(aero_cnst_lq(aero_props%nbins(),0:maxval(aero_props%nmasses())), stat=astat)
         186           0 :     if (astat/=0) then
         187           0 :        call endrun(subrname//' : not able to allocate aero_cnst_lq array')
         188             :     end if
         189           0 :     aero_cnst_lq(:,:) = .false.
         190             :
         191           0 :     allocate(aero_cnst_id(aero_props%nbins(),0:maxval(aero_props%nmasses())), stat=astat)
         192           0 :     if (astat/=0) then
         193           0 :        call endrun(subrname//' : not able to allocate aero_cnst_id array')
         194             :     end if
         195           0 :     aero_cnst_id(:,:) = -1
         196             :
         197           0 :     wetdep_lq = .false.
         198             :
         199           0 :     do m = 1, aero_props%nbins()
         200           0 :        write(binstr,'(i2.2)') m
         201           0 :        call addfld('SOLFACTB'//binstr,  (/ 'lev' /), 'A', '1', 'below cld sol fact')
         202             :
         203           0 :        do l = 0, aero_props%nmasses(m)
         204             :
         205           0 :           if (l == 0) then   ! number
         206           0 :              call aero_props%num_names( m, tmpname, tmpname_cw)
         207             :           else
         208           0 :              call aero_props%mmr_names( m,l, tmpname, tmpname_cw)
         209             :           end if
         210             :
         211           0 :           call cnst_get_ind(tmpname, id, abort=.false.)
         212           0 :           aero_cnst_id(m,l) = id
         213           0 :           aero_cnst_lq(m,l) = id > 0
         214           0 :           if (id > 0) then
         215           0 :              wetdep_lq(id) = .true.
         216             :           end if
         217             :
         218             :           ! units --
         219           0 :           if (l==0) then
         220           0 :              unit_basename = ' 1' ! for num
         221             :           else
         222           0 :              unit_basename = 'kg'
         223             :           endif
         224             :
         225           0 :           call add_hist_fields(tmpname, unit_basename)
         226           0 :           call add_hist_fields(tmpname_cw, unit_basename)
         227             :
         228             :           call addfld( trim(tmpname_cw)//'RSPTD', (/ 'lev' /), 'A', unit_basename//'/kg/s',   &
         229           0 :                        trim(tmpname_cw)//' resuspension tendency')
         230             :
         231             :        end do
         232             :     end do
         233             :
         234           0 :     allocate(scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, aero_props%nbins()), stat=astat)
         235           0 :     if (astat/=0) then
         236           0 :        call endrun(subrname//' : not able to allocate scavimptblnum array')
         237             :     end if
         238           0 :     allocate(scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, aero_props%nbins()), stat=astat)
         239           0 :     if (astat/=0) then
         240           0 :        call endrun(subrname//' : not able to allocate scavimptblvol array')
         241             :     end if
         242           0 :     scavimptblnum = nan
         243           0 :     scavimptblvol = nan
         244             :
         245           0 :     call wetdep_init()
         246             :
         247           0 :     nspec_max = maxval(aero_props%nspecies()) + 2
         248             :
         249           0 :     call init_bcscavcoef()
         250             :
         251           0 :     if (convproc_do_aer) then
         252           0 :        call aero_convproc_init(aero_props)
         253             :     end if
         254             :
         255             :   contains
         256             :
         257           0 :     subroutine add_hist_fields(name,baseunits)
         258             :       character(len=*), intent(in) :: name
         259             :       character(len=*), intent(in) :: baseunits
         260             :
         261             :       call addfld (trim(name)//'SFWET', &
         262           0 :            horiz_only,  'A',baseunits//'/m2/s ','Wet deposition flux at surface')
         263             :       call addfld (trim(name)//'SFSIC', &
         264           0 :            horiz_only,  'A',baseunits//'/m2/s ','Wet deposition flux (incloud, convective) at surface')
         265             :       call addfld (trim(name)//'SFSIS', &
         266           0 :            horiz_only,  'A',baseunits//'/m2/s ','Wet deposition flux (incloud, stratiform) at surface')
         267             :       call addfld (trim(name)//'SFSBC', &
         268           0 :            horiz_only,  'A',baseunits//'/m2/s ','Wet deposition flux (belowcloud, convective) at surface')
         269             :       call addfld (trim(name)//'SFSBS', &
         270           0 :            horiz_only,  'A',baseunits//'/m2/s ','Wet deposition flux (belowcloud, stratiform) at surface')
         271             :
         272           0 :       if (convproc_do_aer) then
         273             :          call addfld (trim(name)//'SFSEC', &
         274           0 :               horiz_only,  'A',unit_basename//'/m2/s','Wet deposition flux (precip evap, convective) at surface')
         275             :          call addfld (trim(name)//'SFSES', &
         276           0 :               horiz_only,  'A',unit_basename//'/m2/s','Wet deposition flux (precip evap, stratiform) at surface')
         277             :          call addfld (trim(name)//'SFSBD', &
         278           0 :               horiz_only,  'A',unit_basename//'/m2/s','Wet deposition flux (belowcloud, deep convective) at surface')
         279             :          call addfld (trim(name)//'WETC',  &
         280           0 :               (/ 'lev' /), 'A',unit_basename//'/kg/s ','wet deposition tendency')
         281             :          call addfld (trim(name)//'CONU',  &
         282           0 :               (/ 'lev' /), 'A',unit_basename//'/kg ','updraft mixing ratio')
         283             :       end if
         284             :
         285           0 :       call addfld (trim(name)//'WET',(/ 'lev' /), 'A',baseunits//'/kg/s ','wet deposition tendency')
         286           0 :       call addfld (trim(name)//'INS',(/ 'lev' /), 'A',baseunits//'/kg/s ','insol frac')
         287             :
         288             :       call addfld (trim(name)//'SIC',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
         289           0 :            trim(name)//' ic wet deposition')
         290             :       call addfld (trim(name)//'SIS',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
         291           0 :            trim(name)//' is wet deposition')
         292             :       call addfld (trim(name)//'SBC',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
         293           0 :            trim(name)//' bc wet deposition')
         294             :       call addfld (trim(name)//'SBS',(/ 'lev' /), 'A',baseunits//'/kg/s ', &
         295           0 :            trim(name)//' bs wet deposition')
         296             :
         297           0 :       if ( history_aerosol .or. history_chemistry ) then
         298           0 :          call add_default (trim(name)//'SFWET', 1, ' ')
         299             :       endif
         300           0 :       if ( history_aerosol ) then
         301           0 :          call add_default (trim(name)//'SFSEC', 1, ' ')
         302           0 :          call add_default (trim(name)//'SFSIC', 1, ' ')
         303           0 :          call add_default (trim(name)//'SFSIS', 1, ' ')
         304           0 :          call add_default (trim(name)//'SFSBC', 1, ' ')
         305           0 :          call add_default (trim(name)//'SFSBS', 1, ' ')
         306           0 :          if (convproc_do_aer) then
         307           0 :             call add_default (trim(name)//'SFSES', 1, ' ')
         308           0 :             call add_default (trim(name)//'SFSBD', 1, ' ')
         309             :          end if
         310             :       endif
         311             :
         312           0 :     end subroutine add_hist_fields
         313             :
         314             :   end subroutine aero_wetdep_init
         315             :
         316             :   !------------------------------------------------------------------------------
         317             :   !------------------------------------------------------------------------------
         318           0 :   subroutine aero_wetdep_tend( state, dt, dlf, cam_out, ptend, pbuf)
         319             :     use wetdep, only: wetdepa_v2, wetdep_inputs_set, wetdep_inputs_t
         320             :     use aerodep_flx, only: aerodep_flx_prescribed
         321             :     use aero_deposition_cam, only: aero_deposition_cam_setwet
         322             :
         323             :     type(physics_state), target, intent(in) :: state  ! Physics state variables
         324             :     real(r8),            intent(in)    :: dt          ! time step
         325             :     real(r8),            intent(in)    :: dlf(:,:)    ! shallow+deep convective detrainment [kg/kg/s]
         326             :     type(cam_out_t),     intent(inout) :: cam_out     ! export state
         327             :     type(physics_ptend), intent(out)   :: ptend       ! indivdual parameterization tendencies
         328             :     type(physics_buffer_desc), pointer :: pbuf(:)
         329             :
         330             :     character(len=*), parameter :: subrname = 'aero_wetdep_tend'
         331             :     type(wetdep_inputs_t) :: dep_inputs
         332           0 :     real(r8), pointer :: fracis(:,:,:)   ! fraction of transported species that are insoluble (pcols, pver, pcnst)
         333             :     real(r8), target :: fracis_nadv(pcols,pver)  ! fraction of not-transported aerosols
         334             :
         335             :     real(r8) :: scavcoefnv(pcols,pver,0:2) ! Dana and Hales coefficient (/mm) for
         336             :                                            ! cloud-borne num & vol (0),
         337             :                                            ! interstitial num (1), interstitial vol (2)
         338             :     integer :: jnv ! index for scavcoefnv 3rd dimension
         339             :     integer :: lphase ! index for interstitial / cloudborne aerosol
         340             :     integer :: strt_loop, end_loop, stride_loop !loop indices for the lphase loop
         341             :
         342             :     real(r8) :: sol_factb(pcols, pver)
         343             :     real(r8) :: sol_facti(pcols, pver)
         344             :     real(r8) :: sol_factic(pcols,pver)
         345             :
         346             :     real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species
         347             :     real(r8) :: rcscavt(pcols, pver)
         348             :     real(r8) :: rsscavt(pcols, pver)
         349             :     real(r8) :: iscavt(pcols, pver)
         350             :     real(r8) :: icscavt(pcols, pver)
         351             :     real(r8) :: isscavt(pcols, pver)
         352             :     real(r8) :: bcscavt(pcols, pver)
         353             :     real(r8) :: bsscavt(pcols, pver)
         354             :
         355           0 :     real(r8) :: diam_wet(state%ncol, pver)
         356             :     logical  :: isprx(pcols,pver) ! true if precipation
         357             :     real(r8) :: prec(pcols) ! precipitation rate
         358             :
         359           0 :     real(r8) :: rtscavt(pcols, pver, 0:nspec_max)
         360             :
         361             :     integer :: ncol, lchnk, m, ndx,mm, l
         362             :     integer :: i,k
         363             :
         364           0 :     real(r8), pointer :: insolfr_ptr(:,:)
         365             :     real(r8) :: q_tmp(pcols,pver) ! temporary array to hold "most current" mixing ratio for 1 species
         366             :     logical :: cldbrn
         367             :
         368           0 :     type(ptr2d_t) :: raer(nele_tot)
         369           0 :     type(ptr2d_t) :: qqcw(nele_tot)
         370             :
         371             :     real(r8) :: sflx(pcols)
         372             :     character(len=aero_name_len) :: aname, cname, name
         373             :
         374           0 :     real(r8) :: qqcw_in(pcols,pver), qqcw_sav(pcols,pver,0:nspec_max)
         375             :     real(r8) :: f_act_conv(pcols,pver) ! prescribed aerosol activation fraction for convective cloud ! rce 2010/05/01
         376             :
         377             :     character(len=2) :: binstr
         378             :     real(r8) :: aerdepwetcw(pcols,pcnst) ! aerosol wet deposition (cloud water)
         379             :     real(r8) :: aerdepwetis(pcols,pcnst)  ! aerosol wet deposition (interstitial)
         380           0 :     real(r8) :: dcondt_resusp3d(nele_tot,pcols,pver)
         381             :
         382             :     integer, parameter :: nsrflx_mzaer2cnvpr = 2
         383           0 :     real(r8) :: qsrflx_mzaer2cnvpr(pcols,nele_tot,nsrflx_mzaer2cnvpr)
         384             :
         385           0 :     real(r8), pointer :: rprddp(:,:)     ! rain production, deep convection
         386           0 :     real(r8), pointer :: rprdsh(:,:)     ! rain production, shallow convection
         387           0 :     real(r8), pointer :: evapcdp(:,:)    ! Evaporation rate of deep    convective precipitation >=0.
         388           0 :     real(r8), pointer :: evapcsh(:,:)    ! Evaporation rate of shallow convective precipitation >=0.
         389             :
         390             :     real(r8) :: rprddpsum(pcols)
         391             :     real(r8) :: rprdshsum(pcols)
         392             :     real(r8) :: evapcdpsum(pcols)
         393             :     real(r8) :: evapcshsum(pcols)
         394             :
         395             :     real(r8) :: tmp_resudp, tmp_resush
         396             :     real(r8) :: tmpa, tmpb
         397             :     real(r8) :: sflxec(pcols), sflxecdp(pcols)  ! deposition flux
         398             :     real(r8) :: sflxic(pcols), sflxicdp(pcols)  ! deposition flux
         399             :     real(r8) :: sflxbc(pcols), sflxbcdp(pcols)  ! deposition flux
         400             :
         401             :     class(aerosol_state), pointer :: aero_state
         402             :
         403           0 :     nullify(aero_state)
         404             :
         405           0 :     if (.not.wetdep_active) return
         406             :
         407           0 :     dcondt_resusp3d(:,:,:) = 0._r8
         408             :
         409           0 :     if (nmodes>0) then
         410           0 :        aero_state => modal_aerosol_state(state,pbuf)
         411           0 :        if (.not.associated(aero_state)) then
         412           0 :           call endrun(subrname//' : construction of aero_state modal_aerosol_state object failed')
         413             :        end if
         414             :     else
         415           0 :        call endrun(subrname//' : cannot determine aerosol model')
         416             :     endif
         417             :
         418           0 :     lchnk = state%lchnk
         419           0 :     ncol = state%ncol
         420             :
         421           0 :     call physics_ptend_init(ptend, state%psetcols, subrname, lq=wetdep_lq)
         422             :
         423           0 :     call wetdep_inputs_set( state, pbuf, dep_inputs )
         424             :
         425           0 :     call pbuf_get_field(pbuf, fracis_idx, fracis)
         426             :
         427           0 :     call aero_state%get_states( aero_props, raer, qqcw )
         428             :
         429           0 :     qsrflx_mzaer2cnvpr(:,:,:) = 0.0_r8
         430           0 :     aerdepwetis(:,:) = 0.0_r8
         431           0 :     aerdepwetcw(:,:) = 0.0_r8
         432             :
         433           0 :     if (convproc_do_aer) then
         434             :        !Do cloudborne first for unified convection scheme so that the resuspension of cloudborne
         435             :        !can be saved then applied to interstitial
         436             :        strt_loop   =  2
         437             :        end_loop    =  1
         438             :        stride_loop = -1
         439             :     else
         440             :        ! Counters for "without" unified convective treatment (i.e. default case)
         441           0 :        strt_loop   = 1
         442           0 :        end_loop    = 2
         443           0 :        stride_loop = 1
         444             :     endif
         445             :
         446           0 :     prec(:ncol)=0._r8
         447           0 :     do k=1,pver
         448           0 :        where (prec(:ncol) >= 1.e-7_r8)
         449             :           isprx(:ncol,k) = .true.
         450             :        elsewhere
         451           0 :           isprx(:ncol,k) = .false.
         452             :        endwhere
         453           0 :        prec(:ncol) = prec(:ncol) + (dep_inputs%prain(:ncol,k) + dep_inputs%cmfdqr(:ncol,k) - dep_inputs%evapr(:ncol,k)) &
         454           0 :                     *state%pdel(:ncol,k)/gravit
         455             :     end do
         456             :
         457           0 :     f_act_conv = 0._r8
         458           0 :     scavcoefnv = nan
         459           0 :     qqcw_sav = nan
         460             :
         461           0 :     if (convproc_do_aer) then
         462             :
         463           0 :        call t_startf('aero_convproc')
         464             :        call aero_convproc_intr( aero_props, aero_state, state, ptend, pbuf, dt, &
         465           0 :             nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr, aerdepwetis, dcondt_resusp3d )
         466             :
         467           0 :        if (convproc_do_evaprain_atonce) then
         468             :
         469           0 :           do m = 1,aero_props%nbins()
         470           0 :              do l = 0,aero_props%nmasses(m)
         471           0 :                 mm = aero_props%indexer(m,l)
         472             :
         473           0 :                 if (l == 0) then   ! number
         474           0 :                    call aero_props%num_names(m, aname, cname)
         475             :                 else
         476           0 :                    call aero_props%mmr_names(m,l, aname, cname)
         477             :                 end if
         478             :
         479           0 :                 call outfld( trim(cname)//'RSPTD', dcondt_resusp3d(mm,:ncol,:), ncol, lchnk )
         480             :
         481           0 :                 do k = 1,pver
         482           0 :                    do i = 1,ncol
         483           0 :                       qqcw(mm)%fld(i,k) = max(0._r8, qqcw(mm)%fld(i,k) + dcondt_resusp3d(mm,i,k)*dt)
         484             :                    end do
         485             :                 end do
         486             :
         487             :              end do
         488             :           end do
         489             :        end if
         490           0 :        call t_stopf('aero_convproc')
         491             :
         492             :     end if
         493             :
         494           0 :     bins_loop: do m = 1,aero_props%nbins()
         495             :
         496           0 :        phase_loop: do lphase = strt_loop,end_loop, stride_loop ! loop over interstitial (1) and cloud-borne (2) forms
         497             :
         498           0 :           cldbrn = lphase==2
         499             :
         500           0 :           sol_factb = nan
         501           0 :           sol_facti = nan
         502           0 :           sol_factic = nan
         503             :
         504           0 :           if (lphase == 1) then ! interstial aerosol
         505             :
         506           0 :              sol_facti = 0.0_r8 ! strat in-cloud scav totally OFF for institial
         507             :
         508           0 :              sol_factic = sol_factic_interstitial
         509             :
         510             :           else ! cloud-borne aerosol (borne by stratiform cloud drops)
         511             :
         512           0 :              sol_factb  = 0.0_r8   ! all below-cloud scav OFF (anything cloud-borne is located "in-cloud")
         513           0 :              sol_facti  = sol_facti_cloud_borne   ! strat  in-cloud scav cloud-borne tuning factor
         514           0 :              sol_factic = 0.0_r8   ! conv   in-cloud scav OFF (having this on would mean
         515             :              !        that conv precip collects strat droplets)
         516           0 :              f_act_conv = 0.0_r8   ! conv   in-cloud scav OFF (having this on would mean
         517             :
         518             :           end if
         519           0 :           if (convproc_do_aer .and. lphase == 1) then
         520             :              ! if modal aero convproc is turned on for aerosols, then
         521             :              !    turn off the convective in-cloud removal for interstitial aerosols
         522             :              !    (but leave the below-cloud on, as convproc only does in-cloud)
         523             :              !    and turn off the outfld SFWET, SFSIC, SFSID, SFSEC, and SFSED calls
         524             :              ! for (stratiform)-cloudborne aerosols, convective wet removal
         525             :              !    (all forms) is zero, so no action is needed
         526           0 :              sol_factic = 0.0_r8
         527             :           endif
         528             :
         529           0 :           diam_wet = aero_state%wet_diameter(m,ncol,pver)
         530             :
         531           0 :           scavcoefnv = 0.0_r8
         532             :
         533           0 :           if (lphase == 1) then ! interstial aerosol
         534           0 :              call get_bcscavcoefs( m, ncol, isprx, diam_wet, scavcoefnv(:,:,1), scavcoefnv(:,:,2) )
         535             :
         536           0 :              if ( sol_factb_interstitial /= NOTSET ) then
         537           0 :                 sol_factb(:ncol,:) = sol_factb_interstitial ! all below-cloud scav
         538             :              else
         539           0 :                 sol_factb(:ncol,:) = aero_state%sol_factb_interstitial( m, ncol, pver, aero_props )
         540             :              end if
         541             :
         542           0 :              write(binstr,'(i2.2)') m
         543           0 :              call outfld('SOLFACTB'//binstr,sol_factb, pcols, lchnk)
         544             :
         545             :           end if
         546             :
         547           0 :           masses_loop: do l = 0,aero_props%nmasses(m)
         548             :
         549           0 :              ndx = aero_cnst_id(m,l)
         550             :
         551           0 :              if (.not. cldbrn .and. ndx>0) then
         552           0 :                 insolfr_ptr => fracis(:,:,ndx)
         553             :              else
         554           0 :                 insolfr_ptr => fracis_nadv
         555             :              endif
         556             :
         557           0 :              mm = aero_props%indexer(m,l)
         558             :
         559           0 :              if (l == 0) then   ! number
         560           0 :                 call aero_props%num_names( m, aname, cname)
         561             :              else
         562           0 :                 call aero_props%mmr_names( m,l, aname, cname)
         563             :              end if
         564             :
         565           0 :              if (cldbrn) then
         566           0 :                 q_tmp(1:ncol,:) = qqcw(mm)%fld(1:ncol,:)
         567           0 :                 jnv = 0
         568           0 :                 if (convproc_do_aer) then
         569           0 :                    qqcw_sav(:ncol,:,l) = q_tmp(1:ncol,:)
         570             :                 endif
         571           0 :                 name = cname
         572           0 :                 qqcw_in = nan
         573           0 :                 f_act_conv = nan
         574             :              else ! interstial aerosol
         575           0 :                 q_tmp(1:ncol,:) = raer(mm)%fld(1:ncol,:) + ptend%q(1:ncol,:,ndx)*dt
         576           0 :                 if (l==0) then
         577             :                    jnv = 1
         578             :                 else
         579           0 :                    jnv = 2
         580             :                 end if
         581           0 :                 if(convproc_do_aer) then
         582             :                    !Feed in the saved cloudborne mixing ratios from phase 2
         583           0 :                    qqcw_in(:ncol,:) = qqcw_sav(:ncol,:,l)
         584             :                 else
         585           0 :                    qqcw_in(:ncol,:) = qqcw(mm)%fld(:ncol,:)
         586             :                 end if
         587             :
         588           0 :                 f_act_conv(:ncol,:) = aero_state%convcld_actfrac( m, l, ncol, pver)
         589           0 :                 name = aname
         590             :              end if
         591             :
         592           0 :              dqdt_tmp(1:ncol,:) = 0.0_r8
         593             :
         594             :              call wetdepa_v2(state%pmid, state%q(:,:,1), state%pdel, &
         595             :                   dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, &
         596             :                   dep_inputs%evapc, dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, &
         597             :                   dep_inputs%evapr, dep_inputs%totcond, q_tmp, dt, &
         598             :                   dqdt_tmp, iscavt, dep_inputs%cldvcu, dep_inputs%cldvst, &
         599             :                   dlf, insolfr_ptr, sol_factb, ncol, &
         600             :                   scavcoefnv(:,:,jnv), &
         601             :                   is_strat_cloudborne=cldbrn, &
         602             :                   qqcw=qqcw_in(:,:), f_act_conv=f_act_conv, &
         603             :                   icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, &
         604             :                   convproc_do_aer=convproc_do_aer, rcscavt=rcscavt, rsscavt=rsscavt,  &
         605             :                   sol_facti_in=sol_facti, sol_factic_in=sol_factic, &
         606             :                   convproc_do_evaprain_atonce_in=convproc_do_evaprain_atonce, &
         607           0 :                   bergso_in=dep_inputs%bergso )
         608             :
         609           0 :              if(convproc_do_aer) then
         610           0 :                 if(cldbrn) then
         611             :                    ! save resuspension of cloudborne species
         612           0 :                    rtscavt(1:ncol,:,l) = rcscavt(1:ncol,:) + rsscavt(1:ncol,:)
         613             :                    ! wetdepa_v2 adds the resuspension of cloudborne to the dqdt of cloudborne (as a source)
         614             :                    ! undo this, so the resuspension of cloudborne can be added to the dqdt of interstitial (above)
         615           0 :                    dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - rtscavt(1:ncol,:,l)
         616             :                 else
         617             :                    ! add resuspension of cloudborne species to dqdt of interstitial species
         618           0 :                    dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt(1:ncol,:,l)
         619             :                 end if
         620             :              endif
         621             :
         622           0 :              if (cldbrn .or. ndx<0) then
         623           0 :                 do k = 1,pver
         624           0 :                    do i = 1,ncol
         625           0 :                       if ( (qqcw(mm)%fld(i,k) + dqdt_tmp(i,k) * dt) .lt. 0.0_r8 )   then
         626           0 :                          dqdt_tmp(i,k) = - qqcw(mm)%fld(i,k) / dt
         627             :                       end if
         628             :                    end do
         629             :                 end do
         630             :
         631           0 :                 qqcw(mm)%fld(1:ncol,:) = qqcw(mm)%fld(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt
         632             :
         633             :              else
         634           0 :                 ptend%q(1:ncol,:,ndx) = ptend%q(1:ncol,:,ndx) + dqdt_tmp(1:ncol,:)
         635             :              end if
         636             :
         637           0 :              call outfld( trim(name)//'WET', dqdt_tmp(:,:), pcols, lchnk)
         638           0 :              call outfld( trim(name)//'SIC', icscavt, pcols, lchnk)
         639           0 :              call outfld( trim(name)//'SIS', isscavt, pcols, lchnk)
         640           0 :              call outfld( trim(name)//'SBC', bcscavt, pcols, lchnk)
         641           0 :              call outfld( trim(name)//'SBS', bsscavt, pcols, lchnk)
         642             :
         643           0 :              call outfld( trim(name)//'INS', insolfr_ptr, pcols, lchnk)
         644             :
         645           0 :              sflx(:)=0._r8
         646           0 :              do k=1,pver
         647           0 :                 do i=1,ncol
         648           0 :                    sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit
         649             :                 enddo
         650             :              enddo
         651           0 :              if (cldbrn) then
         652           0 :                 call outfld( trim(name)//'SFWET', sflx, pcols, lchnk)
         653           0 :                 if (ndx>0) aerdepwetcw(:ncol,ndx) = sflx(:ncol)
         654             :              else
         655           0 :                 if (.not.convproc_do_aer) call outfld( trim(name)//'SFWET', sflx, pcols, lchnk)
         656           0 :                 if (ndx>0) aerdepwetis(:ncol,ndx) = aerdepwetis(:ncol,ndx) + sflx(:ncol)
         657             :              end if
         658             :
         659           0 :              sflx(:)=0._r8
         660           0 :              do k=1,pver
         661           0 :                 do i=1,ncol
         662           0 :                    sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit
         663             :                 enddo
         664             :              enddo
         665           0 :              if (cldbrn) then
         666           0 :                 call outfld( trim(name)//'SFSIC', sflx, pcols, lchnk)
         667             :              else
         668           0 :                 if (.not.convproc_do_aer) call outfld( trim(name)//'SFSIC', sflx, pcols, lchnk)
         669           0 :                 if (convproc_do_aer) sflxic = sflx
         670             :              end if
         671             :
         672           0 :              sflx(:)=0._r8
         673           0 :              do k=1,pver
         674           0 :                 do i=1,ncol
         675           0 :                    sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit
         676             :                 enddo
         677             :              enddo
         678           0 :              call outfld( trim(name)//'SFSIS', sflx, pcols, lchnk)
         679             :
         680           0 :              sflx(:)=0._r8
         681           0 :              do k=1,pver
         682           0 :                 do i=1,ncol
         683           0 :                    sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit
         684             :                 enddo
         685             :              enddo
         686           0 :              call outfld( trim(name)//'SFSBC', sflx, pcols, lchnk)
         687           0 :              if (convproc_do_aer) sflxbc = sflx
         688             :
         689           0 :              sflx(:)=0._r8
         690           0 :              do k=1,pver
         691           0 :                 do i=1,ncol
         692           0 :                    sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit
         693             :                 enddo
         694             :              enddo
         695           0 :              call outfld( trim(name)//'SFSBS', sflx, pcols, lchnk)
         696             :
         697           0 :              if(convproc_do_aer) then
         698             :
         699           0 :                 sflx(:)=0.0_r8
         700           0 :                 do k=1,pver
         701           0 :                    do i=1,ncol
         702           0 :                       sflx(i)=sflx(i)+rsscavt(i,k)*state%pdel(i,k)/gravit
         703             :                    enddo
         704             :                 enddo
         705           0 :                 call outfld( trim(name)//'SFSES', sflx, pcols, lchnk)
         706             :
         707           0 :                 sflx(:)=0._r8
         708           0 :                 do k=1,pver
         709           0 :                    do i=1,ncol
         710           0 :                       sflx(i)=sflx(i)+rcscavt(i,k)*state%pdel(i,k)/gravit
         711             :                    enddo
         712             :                 enddo
         713           0 :                 if (.not.convproc_do_aer) call outfld( trim(name)//'SFSEC', sflx, pcols, lchnk)
         714           0 :                 sflxec = sflx
         715             :
         716           0 :                 if(.not.cldbrn) then
         717             :
         718             :                    ! apportion convective surface fluxes to deep and shallow conv
         719             :                    ! this could be done more accurately in subr wetdepa
         720             :                    ! since deep and shallow rarely occur simultaneously, and these
         721             :                    !    fields are just diagnostics, this approximate method is adequate
         722             :                    ! only do this for interstitial aerosol, because conv clouds to not
         723             :                    !    affect the stratiform-cloudborne aerosol
         724           0 :                    if ( deepconv_wetdep_history) then
         725             :
         726           0 :                       call pbuf_get_field(pbuf, rprddp_idx,      rprddp  )
         727           0 :                       call pbuf_get_field(pbuf, rprdsh_idx,      rprdsh  )
         728           0 :                       call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp )
         729           0 :                       call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh )
         730             :
         731           0 :                       rprddpsum(:)  = 0.0_r8
         732           0 :                       rprdshsum(:)  = 0.0_r8
         733           0 :                       evapcdpsum(:) = 0.0_r8
         734           0 :                       evapcshsum(:) = 0.0_r8
         735             :
         736           0 :                       do k = 1, pver
         737           0 :                          rprddpsum(:ncol)  = rprddpsum(:ncol)  +  rprddp(:ncol,k)*state%pdel(:ncol,k)/gravit
         738           0 :                          rprdshsum(:ncol)  = rprdshsum(:ncol)  +  rprdsh(:ncol,k)*state%pdel(:ncol,k)/gravit
         739           0 :                          evapcdpsum(:ncol) = evapcdpsum(:ncol) + evapcdp(:ncol,k)*state%pdel(:ncol,k)/gravit
         740           0 :                          evapcshsum(:ncol) = evapcshsum(:ncol) + evapcsh(:ncol,k)*state%pdel(:ncol,k)/gravit
         741             :                       end do
         742             :
         743           0 :                       do i = 1, ncol
         744           0 :                          rprddpsum(i)  = max( rprddpsum(i),  1.0e-35_r8 )
         745           0 :                          rprdshsum(i)  = max( rprdshsum(i),  1.0e-35_r8 )
         746           0 :                          evapcdpsum(i) = max( evapcdpsum(i), 0.1e-35_r8 )
         747           0 :                          evapcshsum(i) = max( evapcshsum(i), 0.1e-35_r8 )
         748             :
         749             :                          ! assume that in- and below-cloud removal are proportional to column precip production
         750           0 :                          tmpa = rprddpsum(i) / (rprddpsum(i) + rprdshsum(i))
         751           0 :                          tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
         752           0 :                          sflxicdp(i) = sflxic(i)*tmpa
         753           0 :                          sflxbcdp(i) = sflxbc(i)*tmpa
         754             :
         755             :                          ! assume that resuspension is proportional to (wet removal)*[(precip evap)/(precip production)]
         756           0 :                          tmp_resudp =           tmpa  * min( (evapcdpsum(i)/rprddpsum(i)), 1.0_r8 )
         757           0 :                          tmp_resush = (1.0_r8 - tmpa) * min( (evapcshsum(i)/rprdshsum(i)), 1.0_r8 )
         758           0 :                          tmpb = max( tmp_resudp, 1.0e-35_r8 ) / max( (tmp_resudp+tmp_resush), 1.0e-35_r8 )
         759           0 :                          tmpb = max( 0.0_r8, min( 1.0_r8, tmpb ) )
         760           0 :                          sflxecdp(i) = sflxec(i)*tmpb
         761             :                       end do
         762           0 :                       call outfld( trim(name)//'SFSBD', sflxbcdp, pcols, lchnk)
         763             :                    else
         764           0 :                       sflxec(1:ncol)   = 0.0_r8
         765           0 :                       sflxecdp(1:ncol) = 0.0_r8
         766             :                    end if
         767             :
         768             :                    ! when ma_convproc_intr is used, convective in-cloud wet removal is done there
         769             :                    ! the convective (total and deep) precip-evap-resuspension includes in- and below-cloud
         770             :                    ! contributions
         771             :                    ! so pass the below-cloud contribution to ma_convproc_intr
         772           0 :                    qsrflx_mzaer2cnvpr(1:ncol,mm,1) = sflxec(  1:ncol)
         773           0 :                    qsrflx_mzaer2cnvpr(1:ncol,mm,2) = sflxecdp(1:ncol)
         774             :
         775             :                 end if
         776             :              end if
         777             :
         778             :           end do masses_loop
         779             :        end do phase_loop
         780             :
         781             :     end do bins_loop
         782             :
         783           0 :     if (associated(aero_state)) then
         784           0 :        deallocate(aero_state)
         785             :        nullify(aero_state)
         786             :     end if
         787             :
         788             :     ! if the user has specified prescribed aerosol dep fluxes then
         789             :     ! do not set cam_out dep fluxes according to the prognostic aerosols
         790           0 :     if (.not. aerodep_flx_prescribed()) then
         791           0 :        call aero_deposition_cam_setwet(aerdepwetis, aerdepwetcw, cam_out)
         792             :     endif
         793             :
         794             :   contains
         795             :
         796             :     ! below cloud impaction scavenging coefs
         797           0 :     subroutine get_bcscavcoefs( m, ncol, isprx, diam_wet, scavcoefnum, scavcoefvol )
         798             :
         799             :       integer,intent(in) :: m, ncol
         800             :       logical,intent(in):: isprx(:,:)
         801             :       real(r8), intent(in) :: diam_wet(:,:)
         802             :       real(r8), intent(out) :: scavcoefnum(:,:), scavcoefvol(:,:)
         803             :
         804             :       integer i, k, jgrow
         805             :       real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum
         806             :
         807           0 :       do k = 1, pver
         808           0 :          do i = 1, ncol
         809             :
         810             :             ! do only if no precip
         811           0 :             if ( isprx(i,k) .and. diam_wet(i,k)>0.0_r8) then
         812             :                !
         813             :                ! interpolate table values using log of (actual-wet-size)/(base-dry-size)
         814             :
         815           0 :                dumdgratio = diam_wet(i,k)/aero_props%scav_diam(m)
         816           0 :                if ((dumdgratio >= 0.99_r8) .and. (dumdgratio <= 1.01_r8)) then
         817           0 :                   scavimpvol = scavimptblvol(0,m)
         818           0 :                   scavimpnum = scavimptblnum(0,m)
         819             :                else
         820           0 :                   xgrow = log( dumdgratio ) / dlndg_nimptblgrow
         821           0 :                   jgrow = int( xgrow )
         822           0 :                   if (xgrow < 0._r8) jgrow = jgrow - 1
         823           0 :                   if (jgrow < nimptblgrow_mind) then
         824             :                      jgrow = nimptblgrow_mind
         825             :                      xgrow = jgrow
         826             :                   else
         827           0 :                      jgrow = min( jgrow, nimptblgrow_maxd-1 )
         828             :                   end if
         829             :
         830           0 :                   dumfhi = xgrow - jgrow
         831           0 :                   dumflo = 1._r8 - dumfhi
         832             :
         833           0 :                   scavimpvol = dumflo*scavimptblvol(jgrow,m) + &
         834           0 :                                dumfhi*scavimptblvol(jgrow+1,m)
         835           0 :                   scavimpnum = dumflo*scavimptblnum(jgrow,m) + &
         836           0 :                                dumfhi*scavimptblnum(jgrow+1,m)
         837             :
         838             :                end if
         839             :
         840             :                ! impaction scavenging removal amount for volume
         841           0 :                scavcoefvol(i,k) = exp( scavimpvol )
         842             :                ! impaction scavenging removal amount to number
         843           0 :                scavcoefnum(i,k) = exp( scavimpnum )
         844             :
         845             :             else
         846           0 :                scavcoefvol(i,k) = 0._r8
         847           0 :                scavcoefnum(i,k) = 0._r8
         848             :             end if
         849             :
         850             :          end do
         851             :       end do
         852             :
         853           0 :     end subroutine get_bcscavcoefs
         854             :
         855             :   end subroutine aero_wetdep_tend
         856             :
         857             :   !------------------------------------------------------------------------------
         858             :   !------------------------------------------------------------------------------
         859           0 :   subroutine init_bcscavcoef( )
         860             :     !-----------------------------------------------------------------------
         861             :     !
         862             :     ! Purpose:
         863             :     ! Computes lookup table for aerosol impaction/interception scavenging rates
         864             :     !
         865             :     ! Authors: R. Easter
         866             :     ! Simone Tilmes Nov 2021
         867             :     ! added modifications for bin model, assuming sigma = 1.
         868             :     !
         869             :     !-----------------------------------------------------------------------
         870             :
         871             :     use mo_constants,   only: pi
         872             :
         873             :     !   local variables
         874             :     integer nnfit_maxd
         875             :     parameter (nnfit_maxd=27)
         876             :
         877             :     integer m, jgrow, nnfit
         878             :     integer lunerr
         879             :
         880             :     real(r8) dg0, dg0_cgs, press, dg0_base, &
         881             :          rhodryaero, rhowetaero, rhowetaero_cgs, &
         882             :          scavratenum, scavratevol, logsig,                &
         883             :          temp, wetdiaratio, wetvolratio
         884             :
         885             :     real(r8) :: xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd)
         886             :     real(r8) :: xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd)
         887             :
         888             :     character(len=*), parameter :: subname = 'aero_wetdep_cam::init_bcscavcoef'
         889             :
         890           0 :     lunerr = iulog
         891           0 :     dlndg_nimptblgrow = log( 1.25_r8 )
         892             :
         893             :     ! bin model: main loop over aerosol bins
         894             :
         895           0 :     modeloop: do m = 1, aero_props%nbins()
         896             :
         897             :        ! for setting up the lookup table, use the dry density of the first species
         898             :        ! -- assume the first species of the mode/bin is the dominate species
         899           0 :        call aero_props%get(m,1,density=rhodryaero)
         900             :
         901           0 :        dg0_base = aero_props%scav_diam(m)
         902             :
         903           0 :        logsig = aero_props%alogsig(m)
         904             :
         905           0 :        growloop: do jgrow = nimptblgrow_mind, nimptblgrow_maxd
         906             :
         907           0 :           wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
         908           0 :           dg0 = dg0_base*wetdiaratio
         909             :
         910             :           wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 )
         911           0 :           rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio
         912           0 :           rhowetaero = min( rhowetaero, rhodryaero )
         913             :
         914             :           !
         915             :           !   compute impaction scavenging rates at 1 temp-press pair and save
         916             :           !
         917           0 :           nnfit = 0
         918             :
         919           0 :           temp = 273.16_r8
         920           0 :           press = 0.75e6_r8   ! dynes/cm2
         921           0 :           rhowetaero = rhodryaero
         922             :
         923           0 :           dg0_cgs = dg0*1.0e2_r8   ! m to cm
         924             :
         925           0 :           rhowetaero_cgs = rhowetaero*1.0e-3_r8   ! kg/m3 to g/cm3
         926             :
         927             :           call calc_1_impact_rate( &
         928             :                dg0_cgs, logsig, rhowetaero_cgs, temp, press, &
         929           0 :                scavratenum, scavratevol, lunerr )
         930             :
         931           0 :           nnfit = nnfit + 1
         932             :           if (nnfit > nnfit_maxd) then
         933             :              write(lunerr,9110)
         934             :              call endrun(subname//' : nnfit > nnfit_maxd')
         935             :           end if
         936             : 9110      format( '*** subr. init_bcscavcoef -- nnfit too big' )
         937             :
         938           0 :           xxfitnum(1,nnfit) = 1._r8
         939           0 :           yyfitnum(nnfit) = log( scavratenum )
         940             :
         941           0 :           xxfitvol(1,nnfit) = 1._r8
         942           0 :           yyfitvol(nnfit) = log( scavratevol )
         943             :
         944             :          !depends on both bins and different species
         945           0 :           scavimptblnum(jgrow,m) = yyfitnum(1)
         946           0 :           scavimptblvol(jgrow,m) = yyfitvol(1)
         947             :
         948             :        enddo growloop
         949             :     enddo modeloop
         950             :
         951             :   contains
         952             :
         953             :     !===============================================================================
         954           0 :     subroutine calc_1_impact_rate(          &
         955             :          dg0, logsig, rhoaero, temp, press, &
         956             :          scavratenum, scavratevol, lunerr )
         957             :       !
         958             :       !   routine computes a single impaction scavenging rate
         959             :       ! for precipitation rate of 1 mm/h
         960             :       !
         961             :       !   dg0 = geometric mean diameter of aerosol number size distrib. (cm)
         962             :       !   sigmag = geometric standard deviation of size distrib.
         963             :       !   rhoaero = density of aerosol particles (g/cm^3)
         964             :       !   temp = temperature (K)
         965             :       !   press = pressure (dyne/cm^2)
         966             :       !   scavratenum = number scavenging rate (1/h)
         967             :       !   scavratevol = volume or mass scavenging rate (1/h)
         968             :       !   lunerr = logical unit for error message
         969             :       !
         970             :       use mo_constants, only: boltz_cgs, pi, rhowater => rhoh2o_cgs, rgas => rgas_cgs
         971             :
         972             :       implicit none
         973             :
         974             :       !   subr. parameters
         975             :       integer, intent(in) :: lunerr
         976             :       real(r8), intent(in) :: dg0, logsig, rhoaero, temp, press
         977             :       real(r8), intent(out) :: scavratenum, scavratevol
         978             :
         979             :       !   local variables
         980             :       integer nrainsvmax
         981             :       parameter (nrainsvmax=50)
         982             :       real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
         983             :            vfallrainsv(nrainsvmax)
         984             :
         985             :       integer naerosvmax
         986             :       parameter (naerosvmax=51)
         987             :       real(r8) aaerosv(naerosvmax), &
         988             :            ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)
         989             :
         990             :       integer i, ja, jr, na, nr
         991             :       real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
         992             :       real(r8) anumsum, avolsum, cair, chi
         993             :       real(r8) d, dr, dum, dumfuchs, dx
         994             :       real(r8) ebrown, eimpact, eintercept, etotal, freepath
         995             :       real(r8) precip, precipmmhr, precipsum
         996             :       real(r8) r, rainsweepout, reynolds, rhi, rhoair, rlo, rnumsum
         997             :       real(r8) scavsumnum, scavsumnumbb
         998             :       real(r8) scavsumvol, scavsumvolbb
         999             :       real(r8) schmidt, sqrtreynolds, sstar, stokes, sx
        1000             :       real(r8) taurelax, vfall, vfallstp
        1001             :       real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair
        1002             :
        1003           0 :       rlo = .005_r8
        1004           0 :       rhi = .250_r8
        1005           0 :       dr = 0.005_r8
        1006           0 :       nr = 1 + nint( (rhi-rlo)/dr )
        1007             :       if (nr > nrainsvmax) then
        1008             :          write(lunerr,9110)
        1009             :          call endrun(subname//' : nr > nrainsvmax')
        1010             :       end if
        1011             :
        1012             : 9110  format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' )
        1013             :
        1014           0 :       precipmmhr = 1.0_r8
        1015           0 :       precip = precipmmhr/36000._r8
        1016             :
        1017           0 :       ag0 = dg0/2._r8
        1018           0 :       sx = logsig
        1019           0 :       xg0 = log( ag0 )
        1020           0 :       xg3 = xg0 + 3._r8*sx*sx
        1021             :
        1022           0 :       xlo = xg3 - 4._r8*sx
        1023           0 :       xhi = xg3 + 4._r8*sx
        1024           0 :       dx = 0.2_r8*sx
        1025             :
        1026           0 :       dx = max( 0.2_r8*sx, 0.01_r8 )
        1027           0 :       xlo = xg3 - max( 4._r8*sx, 2._r8*dx )
        1028           0 :       xhi = xg3 + max( 4._r8*sx, 2._r8*dx )
        1029             :
        1030           0 :       na = 1 + nint( (xhi-xlo)/dx )
        1031           0 :       if (na > naerosvmax) then
        1032           0 :          write(lunerr,9120)
        1033           0 :          call endrun(subname//' : na > naerosvmax')
        1034             :       end if
        1035             :
        1036             : 9120  format( '*** subr. calc_1_impact_rate -- na > naerosvmax' )
        1037             :
        1038             :       !   air molar density
        1039           0 :       cair = press/(rgas*temp)
        1040             :       !   air mass density
        1041           0 :       rhoair = 28.966_r8*cair
        1042             :       !   molecular freepath
        1043           0 :       freepath = 2.8052e-10_r8/cair
        1044             :       !   air dynamic viscosity
        1045             :       airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) *    &
        1046           0 :            ((temp/296.16_r8)**1.5_r8)
        1047             :       !   air kinemaic viscosity
        1048           0 :       airkinvisc = airdynvisc/rhoair
        1049             :       !   ratio of water viscosity to air viscosity (from Slinn)
        1050           0 :       xmuwaterair = 60.0_r8
        1051             :
        1052             :       !
        1053             :       !   compute rain drop number concentrations
        1054             :       ! rrainsv = raindrop radius (cm)
        1055             :       ! xnumrainsv = raindrop number concentration (#/cm^3)
        1056             :       !         (number in the bin, not number density)
        1057             :       ! vfallrainsv = fall velocity (cm/s)
        1058             :       !
        1059           0 :       precipsum = 0._r8
        1060           0 :       do i = 1, nr
        1061           0 :          r = rlo + (i-1)*dr
        1062           0 :          rrainsv(i) = r
        1063           0 :          xnumrainsv(i) = exp( -r/2.7e-2_r8 )
        1064             :
        1065           0 :          d = 2._r8*r
        1066           0 :          if (d <= 0.007_r8) then
        1067           0 :             vfallstp = 2.88e5_r8 * d**2._r8
        1068           0 :          else if (d <= 0.025_r8) then
        1069           0 :             vfallstp = 2.8008e4_r8 * d**1.528_r8
        1070           0 :          else if (d <= 0.1_r8) then
        1071           0 :             vfallstp = 4104.9_r8 * d**1.008_r8
        1072           0 :          else if (d <= 0.25_r8) then
        1073           0 :             vfallstp = 1812.1_r8 * d**0.638_r8
        1074             :          else
        1075           0 :             vfallstp = 1069.8_r8 * d**0.235_r8
        1076             :          end if
        1077             :
        1078           0 :          vfall = vfallstp * sqrt(1.204e-3_r8/rhoair)
        1079           0 :          vfallrainsv(i) = vfall
        1080           0 :          precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
        1081             :       end do
        1082           0 :       precipsum = precipsum*pi*1.333333_r8
        1083             :
        1084           0 :       rnumsum = 0._r8
        1085           0 :       do i = 1, nr
        1086           0 :          xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
        1087           0 :          rnumsum = rnumsum + xnumrainsv(i)
        1088             :       end do
        1089             :
        1090             :       !
        1091             :       !   compute aerosol concentrations
        1092             :       ! aaerosv = particle radius (cm)
        1093             :       ! fnumaerosv = fraction of total number in the bin (--)
        1094             :       ! fvolaerosv = fraction of total volume in the bin (--)
        1095             :       !
        1096             :       anumsum = 0._r8
        1097             :       avolsum = 0._r8
        1098           0 :       do i = 1, na
        1099           0 :          x = xlo + (i-1)*dx
        1100           0 :          a = exp( x )
        1101           0 :          aaerosv(i) = a
        1102           0 :          dum = (x - xg0)/sx
        1103           0 :          ynumaerosv(i) = exp( -0.5_r8*dum*dum )
        1104           0 :          yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a
        1105           0 :          anumsum = anumsum + ynumaerosv(i)
        1106           0 :          avolsum = avolsum + yvolaerosv(i)
        1107             :       end do
        1108             :
        1109           0 :       do i = 1, na
        1110           0 :          ynumaerosv(i) = ynumaerosv(i)/anumsum
        1111           0 :          yvolaerosv(i) = yvolaerosv(i)/avolsum
        1112             :       end do
        1113             :
        1114             :       !
        1115             :       !   compute scavenging
        1116             :       !
        1117             :       scavsumnum = 0._r8
        1118             :       scavsumvol = 0._r8
        1119             :       !
        1120             :       !   outer loop for rain drop radius
        1121             :       !
        1122           0 :       jr_loop: do jr = 1, nr
        1123             :
        1124           0 :          r = rrainsv(jr)
        1125           0 :          vfall = vfallrainsv(jr)
        1126             :
        1127           0 :          reynolds = r * vfall / airkinvisc
        1128           0 :          sqrtreynolds = sqrt( reynolds )
        1129             :
        1130             :          !
        1131             :          !   inner loop for aerosol particle radius
        1132             :          !
        1133           0 :          scavsumnumbb = 0._r8
        1134           0 :          scavsumvolbb = 0._r8
        1135             :
        1136           0 :          ja_loop: do ja = 1, na
        1137             :
        1138           0 :             a = aaerosv(ja)
        1139             :
        1140           0 :             chi = a/r
        1141             :
        1142           0 :             dum = freepath/a
        1143           0 :             dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum)
        1144           0 :             taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc)
        1145             :
        1146           0 :             aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8
        1147           0 :             aerodiffus = boltz_cgs*temp*taurelax/aeromass
        1148             :
        1149           0 :             schmidt = airkinvisc/aerodiffus
        1150           0 :             stokes = vfall*taurelax/r
        1151             :
        1152             :             ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) /  &
        1153           0 :                  (reynolds*schmidt)
        1154             :
        1155             :             dum = (1._r8 + 2._r8*xmuwaterair*chi) /         &
        1156           0 :                  (1._r8 + xmuwaterair/sqrtreynolds)
        1157           0 :             eintercept = 4._r8*chi*(chi + dum)
        1158             :
        1159           0 :             dum = log( 1._r8 + reynolds )
        1160           0 :             sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum)
        1161           0 :             eimpact = 0._r8
        1162           0 :             if (stokes > sstar) then
        1163           0 :                dum = stokes - sstar
        1164           0 :                eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8
        1165             :             end if
        1166             :
        1167           0 :             etotal = ebrown + eintercept + eimpact
        1168           0 :             etotal = min( etotal, 1.0_r8 )
        1169             :
        1170           0 :             rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall
        1171             :
        1172           0 :             scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
        1173           0 :             scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
        1174             :
        1175             :          enddo ja_loop
        1176             :
        1177           0 :          scavsumnum = scavsumnum + scavsumnumbb
        1178           0 :          scavsumvol = scavsumvol + scavsumvolbb
        1179             :
        1180             :       enddo jr_loop
        1181             :
        1182           0 :       scavratenum = scavsumnum*3600._r8
        1183           0 :       scavratevol = scavsumvol*3600._r8
        1184             :
        1185           0 :     end subroutine calc_1_impact_rate
        1186             :
        1187             :   end subroutine init_bcscavcoef
        1188             :
        1189             : end module aero_wetdep_cam

Generated by: LCOV version 1.14