LCOV - code coverage report
Current view: | top level - chemistry/bulk_aero - aero_model.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 63 | 427 | 14.8 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 4 | 9 | 44.4 %

              Line data    Source code
           1             : !===============================================================================
           2             : ! Bulk Aerosol Model
           3             : !===============================================================================
           4             : module aero_model
           5             :   use shr_kind_mod,      only: r8 => shr_kind_r8
           6             :   use constituents,      only: pcnst, cnst_name, cnst_get_ind
           7             :   use ppgrid,            only: pcols, pver, pverp
           8             :   use cam_abortutils,    only: endrun
           9             :   use cam_logfile,       only: iulog
          10             :   use perf_mod,          only: t_startf, t_stopf
          11             :   use camsrfexch,        only: cam_in_t, cam_out_t
          12             :   use aerodep_flx,       only: aerodep_flx_prescribed
          13             :   use physics_types,     only: physics_state, physics_ptend, physics_ptend_init
          14             :   use physics_buffer,    only: physics_buffer_desc
          15             :   use physconst,         only: gravit, rair
          16             :   use dust_model,        only: dust_active, dust_names, dust_nbin
          17             :   use seasalt_model,     only: sslt_active=>seasalt_active, seasalt_names, seasalt_nbin
          18             :   use spmd_utils,        only: masterproc
          19             :   use physics_buffer,    only: pbuf_get_field, pbuf_get_index
          20             :   use cam_history,       only: outfld
          21             :   use infnan,            only: nan, assignment(=)
          22             :
          23             :   implicit none
          24             :   private
          25             :
          26             :   public :: aero_model_readnl
          27             :   public :: aero_model_register
          28             :   public :: aero_model_init
          29             :   public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols.
          30             :   public :: aero_model_drydep     ! aerosol dry deposition and sediment
          31             :   public :: aero_model_wetdep     ! aerosol wet removal
          32             :   public :: aero_model_emissions  ! aerosol emissions
          33             :   public :: aero_model_surfarea    ! tropospheric aerosol wet surface area for chemistry
          34             :   public :: aero_model_strat_surfarea   ! stub
          35             :
          36             :   public :: wetdep_lq
          37             :
          38             :  ! Misc private data
          39             :
          40             :   integer :: so4_ndx, cb2_ndx, oc2_ndx, nit_ndx
          41             :   integer :: soa_ndx, soai_ndx, soam_ndx, soab_ndx, soat_ndx, soax_ndx
          42             :
          43             :   ! Namelist variables
          44             :   character(len=16) :: wetdep_list(pcnst) = ' '
          45             :   character(len=16) :: drydep_list(pcnst) = ' '
          46             :
          47             :   integer :: ndrydep = 0
          48             :   integer,allocatable :: drydep_indices(:)
          49             :   integer :: nwetdep = 0
          50             :   integer,allocatable :: wetdep_indices(:)
          51             :   logical :: drydep_lq(pcnst)
          52             :   logical, protected :: wetdep_lq(pcnst)
          53             :
          54             :   integer :: fracis_idx = 0
          55             :
          56             :   real(r8) :: aer_sol_facti(pcnst) ! in-cloud solubility factor
          57             :   real(r8) :: aer_sol_factb(pcnst) ! below-cloud solubility factor
          58             :   real(r8) :: aer_scav_coef(pcnst)
          59             :
          60             : contains
          61             :
          62             :   !=============================================================================
          63             :   ! reads aerosol namelist options
          64             :   !=============================================================================
          65           0 :   subroutine aero_model_readnl(nlfile)
          66             :
          67             :     use namelist_utils,  only: find_group_name
          68             :     use units,           only: getunit, freeunit
          69             :     use mpishorthand
          70             :
          71             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          72             :
          73             :     ! Local variables
          74             :     integer :: unitn, ierr
          75             :     character(len=*), parameter :: subname = 'aero_model_readnl'
          76             :
          77             :     ! Namelist variables
          78             :     character(len=16) :: aer_wetdep_list(pcnst) = ' '
          79             :     character(len=16) :: aer_drydep_list(pcnst) = ' '
          80             :
          81             :     namelist /aerosol_nl/ aer_wetdep_list, aer_drydep_list
          82             :     namelist /aerosol_nl/ aer_sol_facti, aer_sol_factb, aer_scav_coef
          83             :     !-----------------------------------------------------------------------------
          84           0 :     aer_sol_facti = nan
          85           0 :     aer_sol_factb = nan
          86           0 :     aer_scav_coef = nan
          87             :
          88             :     ! Read namelist
          89           0 :     if (masterproc) then
          90           0 :        unitn = getunit()
          91           0 :        open( unitn, file=trim(nlfile), status='old' )
          92           0 :        call find_group_name(unitn, 'aerosol_nl', status=ierr)
          93           0 :        if (ierr == 0) then
          94           0 :           read(unitn, aerosol_nl, iostat=ierr)
          95           0 :           if (ierr /= 0) then
          96           0 :              call endrun(subname // ':: ERROR reading namelist')
          97             :           end if
          98             :        end if
          99           0 :        close(unitn)
         100           0 :        call freeunit(unitn)
         101             :     end if
         102             :
         103             : #ifdef SPMD
         104             :     ! Broadcast namelist variables
         105           0 :     call mpibcast(aer_wetdep_list, len(aer_wetdep_list(1))*pcnst, mpichar, 0, mpicom)
         106           0 :     call mpibcast(aer_drydep_list, len(aer_drydep_list(1))*pcnst, mpichar, 0, mpicom)
         107           0 :     call mpibcast(aer_sol_facti, pcnst, mpir8, 0, mpicom)
         108           0 :     call mpibcast(aer_sol_factb, pcnst, mpir8, 0, mpicom)
         109           0 :     call mpibcast(aer_scav_coef, pcnst, mpir8, 0, mpicom)
         110             : #endif
         111             :
         112           0 :     wetdep_list = aer_wetdep_list
         113           0 :     drydep_list = aer_drydep_list
         114             :
         115           0 :   end subroutine aero_model_readnl
         116             :
         117             :   !=============================================================================
         118             :   !=============================================================================
         119        1536 :   subroutine aero_model_register()
         120             :     use mo_setsoa, only : soa_register
         121             :
         122        1536 :     call soa_register()
         123        1536 :   end subroutine aero_model_register
         124             :
         125             :   !=============================================================================
         126             :   !=============================================================================
         127        1536 :   subroutine aero_model_init( pbuf2d )
         128             :
         129             :     use mo_chem_utls,   only: get_inv_ndx, get_spc_ndx
         130             :     use cam_history,    only: addfld, add_default, horiz_only
         131             :     use phys_control,   only: phys_getopts
         132             :     use mo_aerosols,    only: aerosols_inti
         133             :     use mo_setsoa,      only: soa_inti
         134             :     use dust_model,     only: dust_init
         135             :     use seasalt_model,  only: seasalt_init
         136             :     use aer_drydep_mod, only: inidrydep
         137             :     use wetdep,         only: wetdep_init
         138             :     use mo_setsox,      only: has_sox
         139             :
         140             :     ! args
         141             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         142             :
         143             :     ! local vars
         144             :     character(len=12), parameter :: subrname = 'aero_model_init'
         145             :     integer :: m, id
         146             :     character(len=20) :: dummy
         147             :     logical  :: history_aerosol ! Output MAM or SECT aerosol tendencies
         148             :     logical  :: history_dust    ! Output dust
         149             :
         150             :     call phys_getopts( history_aerosol_out = history_aerosol,&
         151        1536 :                        history_dust_out    = history_dust   )
         152             :
         153        1536 :     call aerosols_inti()
         154        1536 :     call soa_inti(pbuf2d)
         155        1536 :     call dust_init()
         156        1536 :     call seasalt_init()
         157        1536 :     call wetdep_init()
         158             :
         159        1536 :     fracis_idx = pbuf_get_index('FRACIS')
         160             :
         161        1536 :     nwetdep = 0
         162        1536 :     ndrydep = 0
         163             :
         164       18432 :     count_species: do m = 1,pcnst
         165       16896 :        if ( len_trim(wetdep_list(m)) /= 0 ) then
         166           0 :           nwetdep = nwetdep+1
         167             :        endif
         168       18432 :        if ( len_trim(drydep_list(m)) /= 0 ) then
         169           0 :           ndrydep = ndrydep+1
         170             :        endif
         171             :     enddo count_species
         172             :
         173        1536 :     if (nwetdep>0) &
         174           0 :          allocate(wetdep_indices(nwetdep))
         175        1536 :     if (ndrydep>0) &
         176           0 :          allocate(drydep_indices(ndrydep))
         177             :
         178        1536 :     do m = 1,ndrydep
         179           0 :        call cnst_get_ind ( drydep_list(m), id, abort=.false. )
         180           0 :        if (id>0) then
         181           0 :           drydep_indices(m) = id
         182             :        else
         183           0 :           call endrun(subrname//': invalid drydep species: '//trim(drydep_list(m)) )
         184             :        endif
         185             :
         186        1536 :        if (masterproc) then
         187           0 :           write(iulog,*) subrname//': '//drydep_list(m)//' will have drydep applied'
         188             :        endif
         189             :     enddo
         190        1536 :     do m = 1,nwetdep
         191           0 :        call cnst_get_ind ( wetdep_list(m), id, abort=.false. )
         192           0 :        if (id>0) then
         193           0 :           wetdep_indices(m) = id
         194             :        else
         195           0 :           call endrun(subrname//': invalid wetdep species: '//trim(wetdep_list(m)) )
         196             :        endif
         197             :
         198        1536 :        if (masterproc) then
         199           0 :           write(iulog,*) subrname//': '//wetdep_list(m)//' will have wet removal'
         200             :        endif
         201             :     enddo
         202             :
         203             :     ! set flags for drydep tendencies
         204        1536 :     drydep_lq(:) = .false.
         205        1536 :     do m=1,ndrydep
         206           0 :        id = drydep_indices(m)
         207        1536 :        drydep_lq(id) =  .true.
         208             :     enddo
         209             :
         210             :     ! set flags for wetdep tendencies
         211        1536 :     wetdep_lq(:) = .false.
         212        1536 :     do m=1,nwetdep
         213           0 :        id = wetdep_indices(m)
         214        1536 :        wetdep_lq(id) = .true.
         215             :     enddo
         216             :
         217        1536 :     do m = 1,ndrydep
         218             :
         219           0 :        dummy = trim(drydep_list(m)) // 'TB'
         220           0 :        call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m))//' turbulent dry deposition flux')
         221           0 :        if ( history_aerosol ) then
         222           0 :           call add_default (dummy, 1, ' ')
         223             :        endif
         224           0 :        dummy = trim(drydep_list(m))  // 'GV'
         225           0 :        call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m)) //' gravitational dry deposition flux')
         226           0 :        if ( history_aerosol ) then
         227           0 :           call add_default (dummy, 1, ' ')
         228             :        endif
         229           0 :        dummy = trim(drydep_list(m))  // 'DD'
         230           0 :        call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m)) //' dry deposition flux at bottom (grav + turb)')
         231           0 :        if ( history_aerosol ) then
         232           0 :           call add_default (dummy, 1, ' ')
         233             :        endif
         234           0 :        dummy = trim(drydep_list(m)) // 'DT'
         235           0 :        call addfld (dummy,(/ 'lev' /), 'A','kg/kg/s',trim(drydep_list(m))//' dry deposition')
         236           0 :        if ( history_aerosol ) then
         237           0 :           call add_default (dummy, 1, ' ')
         238             :        endif
         239           0 :        dummy = trim(drydep_list(m)) // 'DV'
         240           0 :        call addfld (dummy,(/ 'lev' /), 'A','m/s',trim(drydep_list(m))//' deposition velocity')
         241        1536 :        if ( history_aerosol ) then
         242           0 :           call add_default (dummy, 1, ' ')
         243             :        endif
         244             :
         245             :     enddo
         246             :
         247        1536 :     if (ndrydep>0) then
         248             :
         249           0 :        call inidrydep(rair, gravit)
         250             :
         251           0 :        dummy = 'RAM1'
         252           0 :        call addfld (dummy,horiz_only, 'A','frac','RAM1')
         253           0 :        if ( history_aerosol ) then
         254           0 :           call add_default (dummy, 1, ' ')
         255             :        endif
         256           0 :        dummy = 'airFV'
         257           0 :        call addfld (dummy,horiz_only, 'A','frac','FV')
         258           0 :        if ( history_aerosol ) then
         259           0 :           call add_default (dummy, 1, ' ')
         260             :        endif
         261             :
         262           0 :        if (sslt_active) then
         263           0 :           dummy = 'SSTSFDRY'
         264           0 :           call addfld (dummy,horiz_only, 'A','kg/m2/s','Sea salt deposition flux at surface')
         265           0 :           if ( history_aerosol ) then
         266           0 :              call add_default (dummy, 1, ' ')
         267             :           endif
         268             :        endif
         269           0 :        if (dust_active) then
         270           0 :           dummy = 'DSTSFDRY'
         271           0 :           call addfld (dummy,horiz_only, 'A','kg/m2/s','Dust deposition flux at surface')
         272           0 :           if ( history_aerosol ) then
         273           0 :              call add_default (dummy, 1, ' ')
         274             :           endif
         275             :        endif
         276             :
         277             :     endif
         278             :
         279        1536 :     do m = 1,nwetdep
         280             :
         281           0 :        call addfld (trim(wetdep_list(m))//'SFWET', horiz_only,  'A','kg/m2/s', &
         282           0 :             'Wet deposition flux at surface')
         283             :        call addfld (trim(wetdep_list(m))//'SFSIC', horiz_only,  'A','kg/m2/s', &
         284           0 :             'Wet deposition flux (incloud, convective) at surface')
         285             :        call addfld (trim(wetdep_list(m))//'SFSIS', horiz_only,  'A','kg/m2/s', &
         286           0 :             'Wet deposition flux (incloud, stratiform) at surface')
         287             :        call addfld (trim(wetdep_list(m))//'SFSBC', horiz_only,  'A','kg/m2/s', &
         288           0 :             'Wet deposition flux (belowcloud, convective) at surface')
         289             :        call addfld (trim(wetdep_list(m))//'SFSBS', horiz_only,  'A','kg/m2/s', &
         290           0 :             'Wet deposition flux (belowcloud, stratiform) at surface')
         291             :        call addfld (trim(wetdep_list(m))//'WET',   (/ 'lev' /), 'A','kg/kg/s', &
         292           0 :             'wet deposition tendency')
         293             :        call addfld (trim(wetdep_list(m))//'SIC',   (/ 'lev' /), 'A','kg/kg/s', &
         294           0 :             trim(wetdep_list(m))//' ic wet deposition')
         295             :        call addfld (trim(wetdep_list(m))//'SIS',   (/ 'lev' /), 'A','kg/kg/s', &
         296           0 :             trim(wetdep_list(m))//' is wet deposition')
         297             :        call addfld (trim(wetdep_list(m))//'SBC',   (/ 'lev' /), 'A','kg/kg/s', &
         298           0 :             trim(wetdep_list(m))//' bc wet deposition')
         299             :        call addfld (trim(wetdep_list(m))//'SBS',   (/ 'lev' /), 'A','kg/kg/s', &
         300        1536 :             trim(wetdep_list(m))//' bs wet deposition')
         301             :     enddo
         302             :
         303        1536 :     if (nwetdep>0) then
         304           0 :        if (sslt_active) then
         305           0 :           dummy = 'SSTSFWET'
         306           0 :           call addfld (dummy,horiz_only, 'A','kg/m2/s','Sea salt wet deposition flux at surface')
         307           0 :           if ( history_aerosol ) then
         308           0 :              call add_default (dummy, 1, ' ')
         309             :           endif
         310             :        endif
         311           0 :        if (dust_active) then
         312           0 :           dummy = 'DSTSFWET'
         313           0 :           call addfld (dummy,horiz_only, 'A','kg/m2/s','Dust wet deposition flux at surface')
         314           0 :           if ( history_aerosol ) then
         315           0 :              call add_default (dummy, 1, ' ')
         316             :           endif
         317             :        endif
         318             :     endif
         319             :
         320        1536 :     if (dust_active) then
         321             :        ! emissions diagnostics ....
         322             :
         323           0 :        do m = 1, dust_nbin
         324           0 :           dummy = trim(dust_names(m)) // 'SF'
         325           0 :           call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(dust_names(m))//' dust surface emission')
         326           0 :           if (history_aerosol) then
         327           0 :              call add_default (dummy, 1, ' ')
         328             :           endif
         329             :        enddo
         330             :
         331           0 :        dummy = 'DSTSFMBL'
         332           0 :        call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
         333           0 :        if (history_aerosol .or. history_dust) then
         334           0 :           call add_default (dummy, 1, ' ')
         335             :        endif
         336             :
         337           0 :        dummy = 'LND_MBL'
         338           0 :        call addfld (dummy,horiz_only, 'A','frac','Soil erodibility factor')
         339           0 :        if (history_aerosol) then
         340           0 :           call add_default (dummy, 1, ' ')
         341             :        endif
         342             :
         343             :     endif
         344             :
         345        1536 :     if (sslt_active) then
         346             :
         347           0 :        dummy = 'SSTSFMBL'
         348           0 :        call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
         349           0 :        if (history_aerosol) then
         350           0 :           call add_default (dummy, 1, ' ')
         351             :        endif
         352             :
         353           0 :        do m = 1, seasalt_nbin
         354           0 :           dummy = trim(seasalt_names(m)) // 'SF'
         355           0 :           call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(seasalt_names(m))//' seasalt surface emission')
         356           0 :           if (history_aerosol) then
         357           0 :              call add_default (dummy, 1, ' ')
         358             :           endif
         359             :        enddo
         360             :
         361             :     endif
         362             :
         363        1536 :     if( has_sox ) then
         364        3072 :        call addfld( 'XPH_LWC',(/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc')
         365             :
         366        1536 :        if ( history_aerosol ) then
         367           0 :           call add_default ('XPH_LWC', 1, ' ')
         368             :        endif
         369             :     endif
         370             :
         371        1536 :     so4_ndx    = get_spc_ndx( 'SO4' )
         372        1536 :     soa_ndx    = get_spc_ndx( 'SOA' )
         373        1536 :     soai_ndx   = get_spc_ndx( 'SOAI' )
         374        1536 :     soam_ndx   = get_spc_ndx( 'SOAM' )
         375        1536 :     soab_ndx   = get_spc_ndx( 'SOAB' )
         376        1536 :     soat_ndx   = get_spc_ndx( 'SOAT' )
         377        1536 :     soax_ndx   = get_spc_ndx( 'SOAX' )
         378        1536 :     cb2_ndx    = get_spc_ndx( 'CB2' )
         379        1536 :     oc2_ndx    = get_spc_ndx( 'OC2' )
         380        1536 :     nit_ndx    = get_spc_ndx( 'NH4NO3' )
         381             :
         382        1536 :   end subroutine aero_model_init
         383             :
         384             :   !=============================================================================
         385             :   !=============================================================================
         386    19359288 :   subroutine aero_model_drydep  ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend )
         387             :
         388        1536 :     use dust_sediment_mod, only: dust_sediment_tend
         389             :     use aer_drydep_mod,    only: d3ddflux, calcram
         390             :     use dust_model,        only: dust_depvel, dust_nbin, dust_names
         391             :     use seasalt_model,     only: sslt_depvel=>seasalt_depvel, sslt_nbin=>seasalt_nbin, sslt_names=>seasalt_names
         392             :
         393             :     ! args
         394             :     type(physics_state),    intent(in)    :: state     ! Physics state variables
         395             :     real(r8),               intent(in)    :: obklen(:)
         396             :     real(r8),               intent(in)    :: ustar(:)  ! sfc fric vel
         397             :     type(cam_in_t), target, intent(in)    :: cam_in    ! import state
         398             :     real(r8),               intent(in)    :: dt        ! time step
         399             :     type(cam_out_t),        intent(inout) :: cam_out   ! export state
         400             :     type(physics_ptend),    intent(out)   :: ptend     ! indivdual parameterization tendencies
         401             :     type(physics_buffer_desc),    pointer :: pbuf(:)
         402             :
         403             :   ! local vars
         404     1489176 :     real(r8), pointer :: landfrac(:) ! land fraction
         405     1489176 :     real(r8), pointer :: icefrac(:)  ! ice fraction
         406     1489176 :     real(r8), pointer :: ocnfrac(:)  ! ocean fraction
         407     1489176 :     real(r8), pointer :: fvin(:)     !
         408     1489176 :     real(r8), pointer :: ram1in(:)   ! for dry dep velocities from land model for progseasalts
         409             :
         410             :     real(r8) :: fv(pcols)            ! for dry dep velocities, from land modified over ocean & ice
         411             :     real(r8) :: ram1(pcols)          ! for dry dep velocities, from land modified over ocean & ice
         412             :
         413             :      ! local decarations
         414             :
         415             :     integer, parameter :: naero = sslt_nbin+dust_nbin
         416             :     integer, parameter :: begslt = 1
         417             :     integer, parameter :: endslt = sslt_nbin
         418             :     integer, parameter :: begdst = sslt_nbin+1
         419             :     integer, parameter :: enddst = sslt_nbin+dust_nbin
         420             :
         421             :     integer :: ncol, lchnk
         422             :
         423             :     character(len=6) :: aeronames(naero) ! = (/ sslt_names, dust_names /)
         424             :
         425             :     real(r8) :: vlc_trb(pcols,naero)    !Turbulent deposn velocity (m/s)
         426             :     real(r8) :: vlc_grv(pcols,pver,naero)  !grav deposn velocity (m/s)
         427             :     real(r8) :: vlc_dry(pcols,pver,naero)  !dry deposn velocity (m/s)
         428             :
         429             :     real(r8) :: dep_trb(pcols)       !kg/m2/s
         430             :     real(r8) :: dep_grv(pcols)       !kg/m2/s (total of grav and trb)
         431             :
         432             :     real(r8) :: tsflx_dst(pcols)
         433             :     real(r8) :: tsflx_slt(pcols)
         434             :     real(r8) :: pvaeros(pcols,pverp)    ! sedimentation velocity in Pa
         435             :     real(r8) :: sflx(pcols)
         436             :
         437             :     real(r8) :: tvs(pcols,pver)
         438             :     real(r8) :: rho(pcols,pver)      ! air density in kg/m3
         439             :
         440             :     integer :: m,mm, i, im
         441             :
         442     1489176 :     if (ndrydep<1) return
         443             :
         444           0 :     landfrac => cam_in%landfrac(:)
         445           0 :     icefrac  => cam_in%icefrac(:)
         446           0 :     ocnfrac  => cam_in%ocnfrac(:)
         447           0 :     fvin     => cam_in%fv(:)
         448           0 :     ram1in   => cam_in%ram1(:)
         449             :
         450           0 :     lchnk = state%lchnk
         451           0 :     ncol  = state%ncol
         452             :
         453             :     ! calc ram and fv over ocean and sea ice ...
         454             :     call calcram( ncol,landfrac,icefrac,ocnfrac,obklen,&
         455             :                   ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
         456           0 :                   state%pdel(:,pver),fvin,fv)
         457             :
         458           0 :     call outfld( 'airFV', fv(:), pcols, lchnk )
         459           0 :     call outfld( 'RAM1', ram1(:), pcols, lchnk )
         460             :
         461             :     ! note that tendencies are not only in sfc layer (because of sedimentation)
         462             :     ! and that ptend is updated within each subroutine for different species
         463             :
         464           0 :     call physics_ptend_init(ptend, state%psetcols, 'aero_model_drydep', lq=drydep_lq)
         465             :
         466           0 :     aeronames(:sslt_nbin)   = sslt_names(:)
         467           0 :     aeronames(sslt_nbin+1:) = dust_names(:)
         468             :
         469           0 :     lchnk = state%lchnk
         470           0 :     ncol  = state%ncol
         471             :
         472           0 :     tvs(:ncol,:) = state%t(:ncol,:)
         473           0 :     rho(:ncol,:) = state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
         474             :
         475             :     ! compute dep velocities for sea salt and dust...
         476           0 :     if (sslt_active) then
         477             :        call sslt_depvel( state%t(:,:), state%pmid(:,:), state%q(:,:,1), ram1, fv, ncol, lchnk, &
         478           0 :                          vlc_dry(:,:,begslt:endslt), vlc_trb(:,begslt:endslt), vlc_grv(:,:,begslt:endslt))
         479             :     endif
         480           0 :     if (dust_active) then
         481             :        call dust_depvel( state%t(:,:), state%pmid(:,:),                 ram1, fv, ncol, &
         482           0 :                          vlc_dry(:,:,begdst:enddst), vlc_trb(:,begdst:enddst), vlc_grv(:,:,begdst:enddst) )
         483             :     endif
         484             :
         485           0 :     tsflx_dst(:)=0._r8
         486           0 :     tsflx_slt(:)=0._r8
         487             :
         488             :     ! do drydep for each of the bins of dust and seasalt
         489           0 :     do m=1,ndrydep
         490             :
         491           0 :        mm = drydep_indices(m)
         492           0 :        findindex: do im = 1,naero
         493           0 :          if (trim(cnst_name(mm))==trim(aeronames(im))) exit findindex
         494             :        enddo findindex
         495             :
         496           0 :        pvaeros(:ncol,1)=0._r8
         497           0 :        pvaeros(:ncol,2:pverp) = vlc_dry(:ncol,:,im)
         498             :
         499           0 :        call outfld( trim(cnst_name(mm))//'DV', pvaeros(:,2:pverp), pcols, lchnk )
         500             :
         501             :        if(.true.) then ! use phil's method
         502             :           !      convert from meters/sec to pascals/sec
         503             :           !      pvaeros(:,1) is assumed zero, use density from layer above in conversion
         504           0 :           pvaeros(:ncol,2:pverp) = pvaeros(:ncol,2:pverp) * rho(:ncol,:)*gravit
         505             :
         506             :           !      calculate the tendencies and sfc fluxes from the above velocities
         507             :           call dust_sediment_tend( &
         508             :                ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
         509           0 :                state%q(:,:,mm) , pvaeros  , ptend%q(:,:,mm), sflx  )
         510             :        else   !use charlie's method
         511             :           call d3ddflux(ncol, vlc_dry(:,:,im), state%q(:,:,mm),state%pmid,state%pdel, tvs,sflx,ptend%q(:,:,mm),dt)
         512             :        endif
         513             :        ! apportion dry deposition into turb and gravitational settling for tapes
         514           0 :        do i=1,ncol
         515           0 :           dep_trb(i)=sflx(i)*vlc_trb(i,im)/vlc_dry(i,pver,im)
         516           0 :           dep_grv(i)=sflx(i)*vlc_grv(i,pver,im)/vlc_dry(i,pver,im)
         517             :        enddo
         518             :
         519           0 :        if ( any( sslt_names(:)==trim(cnst_name(mm)) ) ) &
         520           0 :             tsflx_slt(:ncol)=tsflx_slt(:ncol)+sflx(:ncol)
         521           0 :        if ( any( dust_names(:)==trim(cnst_name(mm)) ) ) &
         522           0 :             tsflx_dst(:ncol)=tsflx_dst(:ncol)+sflx(:ncol)
         523             :
         524             :        ! if the user has specified prescribed aerosol dep fluxes then
         525             :        ! do not set cam_out dep fluxes according to the prognostic aerosols
         526           0 :        if (.not. aerodep_flx_prescribed()) then
         527             :           ! set deposition in export state
         528           0 :           if (im==begdst) then
         529           0 :              cam_out%dstdry1(:ncol) = max(sflx(:ncol), 0._r8)
         530           0 :           elseif(im==begdst+1) then
         531           0 :              cam_out%dstdry2(:ncol) = max(sflx(:ncol), 0._r8)
         532           0 :           elseif(im==begdst+2) then
         533           0 :              cam_out%dstdry3(:ncol) = max(sflx(:ncol), 0._r8)
         534           0 :           elseif(im==begdst+3) then
         535           0 :              cam_out%dstdry4(:ncol) = max(sflx(:ncol), 0._r8)
         536             :           endif
         537             :        endif
         538             :
         539           0 :        call outfld( trim(cnst_name(mm))//'DD', sflx, pcols, lchnk)
         540           0 :        call outfld( trim(cnst_name(mm))//'TB', dep_trb, pcols, lchnk )
         541           0 :        call outfld( trim(cnst_name(mm))//'GV', dep_grv, pcols, lchnk )
         542           0 :        call outfld( trim(cnst_name(mm))//'DT', ptend%q(:,:,mm), pcols, lchnk)
         543             :
         544             :     end do
         545             :
         546             :     ! output the total dry deposition
         547           0 :     if (sslt_active) then
         548           0 :        call outfld( 'SSTSFDRY', tsflx_slt, pcols, lchnk)
         549             :     endif
         550           0 :     if (dust_active) then
         551           0 :        call outfld( 'DSTSFDRY', tsflx_dst, pcols, lchnk)
         552             :     endif
         553             :
         554     1489176 :   endsubroutine aero_model_drydep
         555             :
         556             :   !=============================================================================
         557             :   !=============================================================================
         558    17870112 :   subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf)
         559             :
         560             :     use wetdep,        only : wetdepa_v1, wetdep_inputs_set, wetdep_inputs_t
         561             :     use dust_model,    only : dust_names
         562             :     use seasalt_model, only : sslt_names=>seasalt_names
         563             :
         564             :     ! args
         565             :
         566             :     type(physics_state), intent(in)    :: state       ! Physics state variables
         567             :     real(r8),            intent(in)    :: dt          ! time step
         568             :     real(r8),            intent(in)    :: dlf(:,:)    ! shallow+deep convective detrainment [kg/kg/s]
         569             :     type(cam_out_t),     intent(inout) :: cam_out     ! export state
         570             :     type(physics_ptend), intent(out)   :: ptend       ! indivdual parameterization tendencies
         571             :     type(physics_buffer_desc), pointer :: pbuf(:)
         572             :
         573             :     ! local vars
         574             :
         575             :     integer  :: ncol                     ! number of atmospheric columns
         576             :     integer  :: lchnk                    ! chunk identifier
         577             :     integer  :: m,mm, i,k
         578             :
         579             :     real(r8) :: sflx_tot_dst(pcols)
         580             :     real(r8) :: sflx_tot_slt(pcols)
         581             :
         582             :     real(r8) :: iscavt(pcols, pver)
         583             :     real(r8) :: scavt(pcols, pver)
         584             :     real(r8) :: scavcoef(pcols,pver)     ! Dana and Hales coefficient (/mm) (0.1)
         585             :     real(r8) :: sflx(pcols)              ! deposition flux
         586             :
         587             :     real(r8) :: icscavt(pcols, pver)
         588             :     real(r8) :: isscavt(pcols, pver)
         589             :     real(r8) :: bcscavt(pcols, pver)
         590             :     real(r8) :: bsscavt(pcols, pver)
         591             :
         592             :     real(r8) :: sol_factb, sol_facti
         593             :
         594             :     real(r8) :: rainmr(pcols,pver)       ! mixing ratio of rain within cloud volume
         595             :     real(r8) :: cldv(pcols,pver)         ! cloudy volume undergoing scavenging
         596             :     real(r8) :: cldvcu(pcols,pver)       ! Convective precipitation area at the top interface of current layer
         597             :     real(r8) :: cldvst(pcols,pver)       ! Stratiform precipitation area at the top interface of current layer
         598             :
         599     1489176 :     real(r8), pointer :: fracis(:,:,:)   ! fraction of transported species that are insoluble
         600             :
         601             :     type(wetdep_inputs_t) :: dep_inputs  ! obj that contains inputs to wetdepa routine
         602             :
         603     1489176 :     if (nwetdep<1) return
         604             :
         605           0 :     call pbuf_get_field(pbuf, fracis_idx, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/) )
         606             :
         607           0 :     call physics_ptend_init(ptend, state%psetcols, 'aero_model_wetdep', lq=wetdep_lq)
         608             :
         609           0 :     call wetdep_inputs_set( state, pbuf, dep_inputs )
         610             :
         611           0 :     lchnk = state%lchnk
         612           0 :     ncol  = state%ncol
         613             :
         614           0 :     sflx_tot_dst(:) = 0._r8
         615           0 :     sflx_tot_slt(:) = 0._r8
         616             :
         617           0 :     do m = 1, nwetdep
         618             :
         619           0 :        mm = wetdep_indices(m)
         620             :
         621           0 :        sol_factb = aer_sol_factb(m)
         622           0 :        sol_facti = aer_sol_facti(m)
         623             :
         624           0 :        scavcoef(:ncol,:) = aer_scav_coef(m)
         625             :
         626             :        call wetdepa_v1( state%t, state%pmid, state%q(:,:,1), state%pdel, &
         627             :             dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, &
         628             :             dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, &
         629             :             dep_inputs%evapr, dep_inputs%totcond, state%q(:,:,mm), dt, &
         630             :             scavt, iscavt, dep_inputs%cldv, &
         631             :             fracis(:,:,mm), sol_factb, ncol, &
         632             :             scavcoef, &
         633             :             sol_facti_in=sol_facti, &
         634           0 :             icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt )
         635             :
         636           0 :        ptend%q(:ncol,:,mm)=scavt(:ncol,:)
         637             :
         638           0 :        call outfld( trim(cnst_name(mm))//'WET', ptend%q(:,:,mm), pcols, lchnk)
         639           0 :        call outfld( trim(cnst_name(mm))//'SIC', icscavt , pcols, lchnk)
         640           0 :        call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
         641           0 :        call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
         642           0 :        call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)
         643             :
         644           0 :        sflx(:)=0._r8
         645             :
         646           0 :        do k=1,pver
         647           0 :           do i=1,ncol
         648           0 :              sflx(i)=sflx(i)+ptend%q(i,k,mm)*state%pdel(i,k)/gravit
         649             :           enddo
         650             :        enddo
         651           0 :        call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)
         652             :
         653           0 :        if ( any( sslt_names(:)==trim(cnst_name(mm)) ) ) &
         654           0 :             sflx_tot_slt(:ncol) = sflx_tot_slt(:ncol) + sflx(:ncol)
         655           0 :        if ( any( dust_names(:)==trim(cnst_name(mm)) ) ) &
         656           0 :             sflx_tot_dst(:ncol) = sflx_tot_dst(:ncol) + sflx(:ncol)
         657             :
         658             :        ! if the user has specified prescribed aerosol dep fluxes then
         659             :        ! do not set cam_out dep fluxes according to the prognostic aerosols
         660           0 :        if (.not.aerodep_flx_prescribed()) then
         661             :           ! export deposition fluxes to coupler ??? why "-" sign ???
         662           0 :           if (trim(cnst_name(mm))=='CB2') then
         663           0 :              cam_out%bcphiwet(:) = max(-sflx(:), 0._r8)
         664           0 :           elseif (trim(cnst_name(mm))=='OC2') then
         665           0 :              cam_out%ocphiwet(:) = max(-sflx(:), 0._r8)
         666           0 :           elseif (trim(cnst_name(mm))==trim(dust_names(1))) then
         667           0 :              cam_out%dstwet1(:) = max(-sflx(:), 0._r8)
         668           0 :           elseif (trim(cnst_name(mm))==trim(dust_names(2))) then
         669           0 :              cam_out%dstwet2(:) = max(-sflx(:), 0._r8)
         670           0 :           elseif (trim(cnst_name(mm))==trim(dust_names(3))) then
         671           0 :              cam_out%dstwet3(:) = max(-sflx(:), 0._r8)
         672           0 :           elseif (trim(cnst_name(mm))==trim(dust_names(4))) then
         673           0 :              cam_out%dstwet4(:) = max(-sflx(:), 0._r8)
         674             :           endif
         675             :        endif
         676             :
         677             :     enddo
         678             :
         679           0 :     if (sslt_active) then
         680           0 :        call outfld( 'SSTSFWET', sflx_tot_slt, pcols, lchnk)
         681             :     endif
         682           0 :     if (dust_active) then
         683           0 :        call outfld( 'DSTSFWET', sflx_tot_dst, pcols, lchnk)
         684             :     endif
         685             :
         686     1489176 :   endsubroutine aero_model_wetdep
         687             :
         688             :   !-------------------------------------------------------------------------
         689             :   ! provides aerosol surface area info for sectional aerosols
         690             :   ! called from mo_usrrxt
         691             :   !-------------------------------------------------------------------------
         692           0 :   subroutine aero_model_surfarea( &
         693           0 :                   mmr, radmean, relhum, pmid, temp, strato_sad, sulfate,  m, ltrop, &
         694           0 :                   dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_total, reff_trop )
         695             :
         696             :     use mo_constants, only : pi, avo => avogadro
         697             :
         698             :     ! dummy args
         699             :     real(r8), intent(in)    :: pmid(:,:)
         700             :     real(r8), intent(in)    :: temp(:,:)
         701             :     real(r8), intent(in)    :: mmr(:,:,:)
         702             :     real(r8), intent(in)    :: radmean      ! mean radii in cm
         703             :     real(r8), intent(in)    :: strato_sad(:,:)
         704             :     integer,  intent(in)    :: ncol
         705             :     integer,  intent(in)    :: ltrop(:)
         706             :     real(r8), intent(in)    :: dlat(:)                    ! degrees latitude
         707             :     integer,  intent(in)    :: het1_ndx
         708             :     real(r8), intent(in)    :: relhum(:,:)
         709             :     real(r8), intent(in)    :: m(:,:) ! total atm density (/cm^3)
         710             :     real(r8), intent(in)    :: sulfate(:,:)
         711             :     type(physics_buffer_desc), pointer :: pbuf(:)
         712             :
         713             :     real(r8), intent(inout) :: sfc(:,:,:)
         714             :     real(r8), intent(inout) :: dm_aer(:,:,:)
         715             :     real(r8), intent(inout) :: sad_total(:,:)
         716             :     real(r8), intent(out)   :: reff_trop(:,:)
         717             :
         718             :     ! local vars
         719             :
         720             :     integer  :: i,k
         721             :     real(r8) :: rho_air
         722             :     real(r8) :: v, n, n_exp, r_rd, r_sd
         723             :     real(r8) :: dm_sulf, dm_sulf_wet, log_sd_sulf, sfc_sulf, sfc_nit
         724             :     real(r8) :: dm_orgc, dm_orgc_wet, log_sd_orgc, sfc_oc, sfc_soa
         725             :     real(r8) :: sfc_soai, sfc_soam, sfc_soab, sfc_soat, sfc_soax
         726             :     real(r8) :: dm_bc, dm_bc_wet, log_sd_bc, sfc_bc
         727             :     real(r8) :: rxt_sulf, rxt_nit, rxt_oc, rxt_soa
         728             :     real(r8) :: c_n2o5, c_ho2, c_no2, c_no3
         729             :     real(r8) :: s_exp
         730             :
         731             :     !-----------------------------------------------------------------
         732             :     !   ... parameters for log-normal distribution by number
         733             :     ! references:
         734             :     !   Chin et al., JAS, 59, 461, 2003
         735             :     !   Liao et al., JGR, 108(D1), 4001, 2003
         736             :     !   Martin et al., JGR, 108(D3), 4097, 2003
         737             :     !-----------------------------------------------------------------
         738             :     real(r8), parameter :: rm_sulf  = 6.95e-6_r8        ! mean radius of sulfate particles (cm) (Chin)
         739             :     real(r8), parameter :: sd_sulf  = 2.03_r8           ! standard deviation of radius for sulfate (Chin)
         740             :     real(r8), parameter :: rho_sulf = 1.7e3_r8          ! density of sulfate aerosols (kg/m3) (Chin)
         741             :
         742             :     real(r8), parameter :: rm_orgc  = 2.12e-6_r8        ! mean radius of organic carbon particles (cm) (Chin)
         743             :     real(r8), parameter :: sd_orgc  = 2.20_r8           ! standard deviation of radius for OC (Chin)
         744             :     real(r8), parameter :: rho_orgc = 1.8e3_r8          ! density of OC aerosols (kg/m3) (Chin)
         745             :
         746             :     real(r8), parameter :: rm_bc    = 1.18e-6_r8        ! mean radius of soot/BC particles (cm) (Chin)
         747             :     real(r8), parameter :: sd_bc    = 2.00_r8           ! standard deviation of radius for BC (Chin)
         748             :     real(r8), parameter :: rho_bc   = 1.0e3_r8          ! density of BC aerosols (kg/m3) (Chin)
         749             :
         750             :     real(r8), parameter :: mw_so4 = 98.e-3_r8     ! so4 molecular wt (kg/mole)
         751             :
         752             :     integer  ::  irh, rh_l, rh_u
         753             :     real(r8) ::  factor, rfac_sulf, rfac_oc, rfac_bc, rfac_ss
         754             :     logical :: zero_aerosols
         755             :
         756             :     !-----------------------------------------------------------------
         757             :     !   ... table for hygroscopic growth effect on radius (Chin et al)
         758             :     !           (no growth effect for mineral dust)
         759             :     !-----------------------------------------------------------------
         760             :     real(r8), dimension(7) :: table_rh, table_rfac_sulf, table_rfac_bc, table_rfac_oc, table_rfac_ss
         761             :
         762             :     data table_rh(1:7)        / 0.0_r8, 0.5_r8, 0.7_r8, 0.8_r8, 0.9_r8, 0.95_r8, 0.99_r8/
         763             :     data table_rfac_sulf(1:7) / 1.0_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 1.9_r8,  2.2_r8/
         764             :     data table_rfac_oc(1:7)   / 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8,  2.2_r8/
         765             :     data table_rfac_bc(1:7)   / 1.0_r8, 1.0_r8, 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8,  1.9_r8/
         766             :     data table_rfac_ss(1:7)   / 1.0_r8, 1.6_r8, 1.8_r8, 2.0_r8, 2.4_r8, 2.9_r8,  4.8_r8/
         767             :
         768             :     !-----------------------------------------------------------------
         769             :     !   ... exponent for calculating number density
         770             :     !-----------------------------------------------------------------
         771           0 :     n_exp = exp( -4.5_r8*log(sd_sulf)*log(sd_sulf) )
         772             :
         773           0 :     dm_sulf = 2._r8 * rm_sulf
         774           0 :     dm_orgc = 2._r8 * rm_orgc
         775           0 :     dm_bc   = 2._r8 * rm_bc
         776             :
         777           0 :     log_sd_sulf = log(sd_sulf)
         778           0 :     log_sd_orgc = log(sd_orgc)
         779           0 :     log_sd_bc   = log(sd_bc)
         780             :
         781           0 :     reff_trop(:,:) = 0._r8
         782             :
         783           0 :     ver_loop: do k = 1,pver
         784           0 :        col_loop: do i = 1,ncol
         785             :           !-------------------------------------------------------------------------
         786             :           !     ... air density (kg/m3)
         787             :           !-------------------------------------------------------------------------
         788           0 :           rho_air = pmid(i,k)/(temp(i,k)*287.04_r8)
         789             :           !-------------------------------------------------------------------------
         790             :           !       ... aerosol growth interpolated from M.Chin's table
         791             :           !-------------------------------------------------------------------------
         792           0 :           if (relhum(i,k) >= table_rh(7)) then
         793           0 :              rfac_sulf = table_rfac_sulf(7)
         794           0 :              rfac_oc = table_rfac_oc(7)
         795           0 :              rfac_bc = table_rfac_bc(7)
         796             :           else
         797           0 :              do irh = 2,7
         798           0 :                 if (relhum(i,k) <= table_rh(irh)) then
         799             :                    exit
         800             :                 end if
         801             :              end do
         802           0 :              rh_l = irh-1
         803           0 :              rh_u = irh
         804             :
         805           0 :              factor = (relhum(i,k) - table_rh(rh_l))/(table_rh(rh_u) - table_rh(rh_l))
         806             :
         807           0 :              rfac_sulf = table_rfac_sulf(rh_l) + factor*(table_rfac_sulf(rh_u) - table_rfac_sulf(rh_l))
         808           0 :              rfac_oc = table_rfac_oc(rh_u) + factor*(table_rfac_oc(rh_u) - table_rfac_oc(rh_l))
         809           0 :              rfac_bc = table_rfac_bc(rh_u) + factor*(table_rfac_bc(rh_u) - table_rfac_bc(rh_l))
         810             :           end if
         811             :
         812           0 :           dm_sulf_wet = dm_sulf * rfac_sulf
         813           0 :           dm_orgc_wet = dm_orgc * rfac_oc
         814           0 :           dm_bc_wet = dm_bc * rfac_bc
         815             :
         816           0 :           dm_bc_wet   = min(dm_bc_wet  ,50.e-6_r8) ! maximum size is 0.5 micron (Chin)
         817           0 :           dm_orgc_wet = min(dm_orgc_wet,50.e-6_r8) ! maximum size is 0.5 micron (Chin)
         818             :
         819             :
         820             :           !-------------------------------------------------------------------------
         821             :           !     ... sulfate aerosols
         822             :           !-------------------------------------------------------------------------
         823           0 :           zero_aerosols = k < ltrop(i)
         824           0 :           if ( abs( dlat(i) ) > 50._r8 ) then
         825           0 :              zero_aerosols = pmid(i,k) < 30000._r8
         826             :           endif
         827             :           !-------------------------------------------------------------------------
         828             :           !       ... use ubvals climatology for stratospheric sulfate surface area density
         829             :           !-------------------------------------------------------------------------
         830           0 :           if( zero_aerosols ) then
         831           0 :              sfc_sulf = strato_sad(i,k)
         832           0 :              if ( het1_ndx > 0 ) then
         833           0 :                 sfc_sulf = 0._r8        ! reaction already taken into account in mo_strato_rates.F90
         834             :              end if
         835             :              sfc_nit = 0._r8
         836             :              sfc_soa = 0._r8
         837             :              sfc_oc  = 0._r8
         838             :              sfc_bc  = 0._r8
         839             :           else
         840             :
         841           0 :              if( so4_ndx > 0 ) then
         842             :                 !-------------------------------------------------------------------------
         843             :                 ! convert mass mixing ratio of aerosol to cm3/cm3 (cm^3_aerosol/cm^3_air)
         844             :                 ! v=volume density (m^3/m^3)
         845             :                 ! rho_aer=density of aerosol (kg/m^3)
         846             :                 ! v=m*rho_air/rho_aer   [kg/kg * (kg/m3)_air/(kg/m3)_aer]
         847             :                 !-------------------------------------------------------------------------
         848           0 :                 v = mmr(i,k,so4_ndx) * rho_air/rho_sulf
         849             :                 !-------------------------------------------------------------------------
         850             :                 ! calculate the number density of aerosol (aerosols/cm3)
         851             :                 ! assuming a lognormal distribution
         852             :                 ! n  = (aerosols/cm3)
         853             :                 ! dm = geometric mean diameter
         854             :                 !
         855             :                 ! because only the dry mass of the aerosols is known, we
         856             :                 ! use the mean dry radius
         857             :                 !-------------------------------------------------------------------------
         858           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
         859             :                 !-------------------------------------------------------------------------
         860             :                 ! find surface area of aerosols using dm_wet, log_sd
         861             :                 !  (increase of sd due to RH is negligible)
         862             :                 ! and number density calculated above as distribution
         863             :                 ! parameters
         864             :                 ! sfc = surface area of wet aerosols (cm^2/cm^3)
         865             :                 !-------------------------------------------------------------------------
         866           0 :                 s_exp    = exp(2._r8*log_sd_sulf*log_sd_sulf)
         867           0 :                 sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp
         868             :
         869             :              else
         870             :                 !-------------------------------------------------------------------------
         871             :                 !  if so4 not simulated, use off-line sulfate and calculate as above
         872             :                 !  convert sulfate vmr to volume density of aerosol (cm^3_aerosol/cm^3_air)
         873             :                 !-------------------------------------------------------------------------
         874           0 :                 v = sulfate(i,k) * m(i,k) * mw_so4 / (avo * rho_sulf) *1.e6_r8
         875           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
         876           0 :                 s_exp    = exp(2._r8*log_sd_sulf*log_sd_sulf)
         877           0 :                 sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp
         878             :
         879             :              end if
         880             :
         881             :              !-------------------------------------------------------------------------
         882             :              ! ammonium nitrate (follow same procedure as sulfate, using size and density of sulfate)
         883             :              !-------------------------------------------------------------------------
         884           0 :              if( nit_ndx > 0 ) then
         885           0 :                 v = mmr(i,k,nit_ndx) * rho_air/rho_sulf
         886           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
         887           0 :                 s_exp   = exp(2._r8*log_sd_sulf*log_sd_sulf)
         888           0 :                 sfc_nit = n * pi * (dm_sulf_wet**2._r8) * s_exp
         889             :              else
         890             :                 sfc_nit = 0._r8
         891             :              end if
         892             :
         893             :              !-------------------------------------------------------------------------
         894             :              ! hydrophylic organic carbon (follow same procedure as sulfate)
         895             :              !-------------------------------------------------------------------------
         896           0 :              if( oc2_ndx > 0 ) then
         897           0 :                 v = mmr(i,k,oc2_ndx) * rho_air/rho_orgc
         898           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3))*n_exp
         899           0 :                 s_exp    = exp(2._r8*log_sd_orgc*log_sd_orgc)
         900           0 :                 sfc_oc   = n * pi * (dm_orgc_wet**2._r8) * s_exp
         901             :              else
         902             :                 sfc_oc = 0._r8
         903             :              end if
         904             :
         905             :              !-------------------------------------------------------------------------
         906             :              ! secondary organic carbon (follow same procedure as sulfate)
         907             :              !-------------------------------------------------------------------------
         908           0 :              if( soa_ndx > 0 ) then
         909           0 :                 v = mmr(i,k,soa_ndx) * rho_air/rho_orgc
         910           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
         911           0 :                 s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
         912           0 :                 sfc_soa   = n * pi * (dm_orgc_wet**2._r8) * s_exp
         913             :              else
         914             :                 sfc_soa = 0._r8
         915             :              end if
         916             :
         917             :              !-------------------------------------------------------------------------
         918             :              ! black carbon (follow same procedure as sulfate)
         919             :              !-------------------------------------------------------------------------
         920           0 :              if( cb2_ndx > 0 ) then
         921           0 :                 v = mmr(i,k,cb2_ndx) * rho_air/rho_bc
         922           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_bc**3._r8))*n_exp
         923           0 :                 s_exp     = exp(2._r8*log_sd_bc*log_sd_bc)
         924           0 :                 sfc_bc   = n * pi * (dm_bc_wet**2._r8) * s_exp
         925             :              else
         926             :                 sfc_bc = 0._r8
         927             :              end if
         928           0 :              if( soai_ndx > 0 ) then
         929           0 :                 v = mmr(i,k,soai_ndx) * rho_air/rho_orgc
         930           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
         931           0 :                 s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
         932           0 :                 sfc_soai   = n * pi * (dm_orgc_wet**2._r8) * s_exp
         933             :              else
         934             :                 sfc_soai = 0._r8
         935             :              end if
         936           0 :              if( soam_ndx > 0 ) then
         937           0 :                 v = mmr(i,k,soam_ndx) * rho_air/rho_orgc
         938           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
         939           0 :                 s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
         940           0 :                 sfc_soam   = n * pi * (dm_orgc_wet**2._r8) * s_exp
         941             :              else
         942             :                 sfc_soam = 0._r8
         943             :              end if
         944           0 :              if( soab_ndx > 0 ) then
         945           0 :                 v = mmr(i,k,soab_ndx) * rho_air/rho_orgc
         946           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
         947           0 :                 s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
         948           0 :                 sfc_soab   = n * pi * (dm_orgc_wet**2._r8) * s_exp
         949             :              else
         950             :                 sfc_soab = 0._r8
         951             :              end if
         952           0 :              if( soat_ndx > 0 ) then
         953           0 :                 v = mmr(i,k,soat_ndx) * rho_air/rho_orgc
         954           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
         955           0 :                 s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
         956           0 :                 sfc_soat   = n * pi * (dm_orgc_wet**2._r8) * s_exp
         957             :              else
         958             :                 sfc_soat = 0._r8
         959             :              end if
         960           0 :              if( soax_ndx > 0 ) then
         961           0 :                 v = mmr(i,k,soax_ndx) * rho_air/rho_orgc
         962           0 :                 n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
         963           0 :                 s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
         964           0 :                 sfc_soax   = n * pi * (dm_orgc_wet**2._r8) * s_exp
         965             :              else
         966             :                 sfc_soax = 0._r8
         967             :              end if
         968           0 :              sfc_soa = sfc_soa + sfc_soai + sfc_soam + sfc_soab + sfc_soat + sfc_soax
         969             :
         970             :           end if
         971             :
         972           0 :           sfc(i,k,:) = (/ sfc_sulf, sfc_nit, sfc_oc, sfc_soa, sfc_bc /)
         973           0 :           dm_aer(i,k,:) = (/ dm_sulf_wet,dm_sulf_wet,dm_orgc_wet,dm_orgc_wet,dm_bc_wet /)
         974             :
         975             :           !-------------------------------------------------------------------------
         976             :           !     ... add up total surface area density for output
         977             :           !-------------------------------------------------------------------------
         978           0 :           sad_total(i,k) = sfc_sulf + sfc_nit + sfc_oc + sfc_soa + sfc_bc
         979             :
         980             :        enddo col_loop
         981             :     enddo ver_loop
         982             :
         983           0 :   end subroutine aero_model_surfarea
         984             :
         985             :   !-------------------------------------------------------------------------
         986             :   ! stub
         987             :   !-------------------------------------------------------------------------
         988           0 :   subroutine aero_model_strat_surfarea( ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat )
         989             :
         990             :     ! dummy args
         991             :     integer,  intent(in)    :: ncol
         992             :     real(r8), intent(in)    :: mmr(:,:,:)
         993             :     real(r8), intent(in)    :: pmid(:,:)
         994             :     real(r8), intent(in)    :: temp(:,:)
         995             :     integer,  intent(in)    :: ltrop(:) ! tropopause level indices
         996             :     type(physics_buffer_desc), pointer :: pbuf(:)
         997             :     real(r8), intent(out)   :: strato_sad(:,:)
         998             :     real(r8), intent(out)   :: reff_strat(:,:)
         999             :
        1000           0 :     strato_sad(:,:) = 0._r8
        1001           0 :     reff_strat(:,:) = 0._r8
        1002             :
        1003           0 :   end subroutine aero_model_strat_surfarea
        1004             :
        1005             :   !=============================================================================
        1006             :   !=============================================================================
        1007           0 :   subroutine aero_model_gasaerexch( loffset, ncol, lchnk, troplev, delt, reaction_rates, &
        1008           0 :                                     tfld, pmid, pdel, mbar, relhum, &
        1009           0 :                                     zm,  qh2o, cwat, cldfr, cldnum, &
        1010           0 :                                     airdens, invariants, del_h2so4_gasprod,  &
        1011           0 :                                     vmr0, vmr, pbuf )
        1012             :
        1013             :     use chem_mods,   only : gas_pcnst
        1014             :     use mo_aerosols, only : aerosols_formation, has_aerosols
        1015             :     use mo_setsox,   only : setsox, has_sox
        1016             :     use mo_setsoa,   only : setsoa, has_soa
        1017             :
        1018             :     !-----------------------------------------------------------------------
        1019             :     !      ... dummy arguments
        1020             :     !-----------------------------------------------------------------------
        1021             :     integer,  intent(in) :: loffset                ! offset applied to modal aero "pointers"
        1022             :     integer,  intent(in) :: ncol                   ! number columns in chunk
        1023             :     integer,  intent(in) :: lchnk                  ! chunk index
        1024             :     integer,  intent(in) :: troplev(:)
        1025             :     real(r8), intent(in) :: delt                   ! time step size (sec)
        1026             :     real(r8), intent(in) :: reaction_rates(:,:,:)  ! reaction rates
        1027             :     real(r8), intent(in) :: tfld(:,:)              ! temperature (K)
        1028             :     real(r8), intent(in) :: pmid(:,:)              ! pressure at model levels (Pa)
        1029             :     real(r8), intent(in) :: pdel(:,:)              ! pressure thickness of levels (Pa)
        1030             :     real(r8), intent(in) :: mbar(:,:)              ! mean wet atmospheric mass ( amu )
        1031             :     real(r8), intent(in) :: relhum(:,:)            ! relative humidity
        1032             :     real(r8), intent(in) :: airdens(:,:)           ! total atms density (molec/cm**3)
        1033             :     real(r8), intent(in) :: invariants(:,:,:)
        1034             :     real(r8), intent(in) :: del_h2so4_gasprod(:,:)
        1035             :     real(r8), intent(in) :: zm(:,:)
        1036             :     real(r8), intent(in) :: qh2o(:,:)
        1037             :     real(r8), intent(in) :: cwat(:,:)          ! cloud liquid water content (kg/kg)
        1038             :     real(r8), intent(in) :: cldfr(:,:)
        1039             :     real(r8), intent(in) :: cldnum(:,:)       ! droplet number concentration (#/kg)
        1040             :     real(r8), intent(in) :: vmr0(:,:,:)       ! initial mixing ratios (before gas-phase chem changes)
        1041             :     real(r8), intent(inout) :: vmr(:,:,:)         ! mixing ratios ( vmr )
        1042             :
        1043             :     type(physics_buffer_desc), pointer :: pbuf(:)
        1044             :
        1045             :     ! local vars
        1046             :
        1047           0 :     real(r8) :: vmrcw(ncol,pver,gas_pcnst)            ! cloud-borne aerosol (vmr)
        1048             :
        1049           0 :     real(r8) ::  aqso4(ncol,1)               ! aqueous phase chemistry
        1050           0 :     real(r8) ::  aqh2so4(ncol,1)             ! aqueous phase chemistry
        1051           0 :     real(r8) ::  aqso4_h2o2(ncol)            ! SO4 aqueous phase chemistry due to H2O2
        1052           0 :     real(r8) ::  aqso4_o3(ncol)              ! SO4 aqueous phase chemistry due to O3
        1053           0 :     real(r8) ::  xphlwc(ncol,pver)           ! pH value multiplied by lwc
        1054             :
        1055             :
        1056             :   ! aqueous chemistry ...
        1057             :
        1058           0 :     if( has_sox ) then
        1059             :        call setsox(   &
        1060             :             ncol,     &
        1061             :             lchnk,    &
        1062             :             loffset,  &
        1063             :             delt,     &
        1064             :             pmid,     &
        1065             :             pdel,     &
        1066             :             tfld,     &
        1067             :             mbar,     &
        1068             :             cwat,     &
        1069             :             cldfr,    &
        1070             :             cldnum,   &
        1071             :             airdens,  &
        1072             :             invariants, &
        1073             :             vmrcw,    &
        1074             :             vmr,      &
        1075             :             xphlwc,   &
        1076             :             aqso4,    &
        1077             :             aqh2so4,  &
        1078             :             aqso4_h2o2,&
        1079             :             aqso4_o3  &
        1080           0 :             )
        1081           0 :         call outfld( 'XPH_LWC',xphlwc(:ncol,:), ncol , lchnk )
        1082             :     endif
        1083             :
        1084           0 :     if( has_soa ) then
        1085           0 :        call setsoa( ncol, lchnk, delt, reaction_rates, tfld, airdens, vmr, pbuf)
        1086             :     endif
        1087             :
        1088           0 :     if( has_aerosols ) then
        1089           0 :        call aerosols_formation( ncol, lchnk, tfld, relhum, vmr )
        1090             :     endif
        1091             :
        1092             :
        1093           0 :   end subroutine aero_model_gasaerexch
        1094             :
        1095             :   !=============================================================================
        1096             :   !=============================================================================
        1097           0 :   subroutine aero_model_emissions( state, cam_in )
        1098             :     use seasalt_model, only: seasalt_emis, seasalt_indices
        1099             :     use dust_model,    only: dust_emis, dust_indices
        1100             :     use physics_types, only: physics_state
        1101             :
        1102             :     ! Arguments:
        1103             :
        1104             :     type(physics_state),    intent(in)    :: state   ! Physics state variables
        1105             :     type(cam_in_t),         intent(inout) :: cam_in  ! import state
        1106             :
        1107             :     ! local vars
        1108             :
        1109             :     integer :: lchnk, ncol
        1110             :     integer :: m, mm
        1111             :     real(r8) :: soil_erod_tmp(pcols)
        1112             :     real(r8) :: sflx(pcols)   ! accumulate over all bins for output
        1113             :     real(r8) :: u10cubed(pcols)
        1114             :     real (r8), parameter :: z0=0.0001_r8  ! m roughness length over oceans--from ocean model
        1115             :
        1116           0 :     lchnk = state%lchnk
        1117           0 :     ncol = state%ncol
        1118             :
        1119           0 :     if (dust_active) then
        1120             :
        1121           0 :        call dust_emis( ncol, lchnk, cam_in%dstflx, cam_in%cflx, soil_erod_tmp )
        1122             :
        1123             :        ! some dust emis diagnostics ...
        1124           0 :        sflx(:)=0._r8
        1125           0 :        do m=1,dust_nbin
        1126           0 :           mm = dust_indices(m)
        1127           0 :           sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
        1128           0 :           call outfld(trim(dust_names(m))//'SF',cam_in%cflx(:,mm),pcols, lchnk)
        1129             :        enddo
        1130           0 :        call outfld('DSTSFMBL',sflx(:),pcols,lchnk)
        1131           0 :        call outfld('LND_MBL',soil_erod_tmp(:),pcols, lchnk )
        1132             :     endif
        1133             :
        1134           0 :     if (sslt_active) then
        1135           0 :        u10cubed(:ncol)=sqrt(state%u(:ncol,pver)**2+state%v(:ncol,pver)**2)
        1136             :        ! move the winds to 10m high from the midpoint of the gridbox:
        1137             :        ! follows Tie and Seinfeld and Pandis, p.859 with math.
        1138             :
        1139           0 :        u10cubed(:ncol)=u10cubed(:ncol)*log(10._r8/z0)/log(state%zm(:ncol,pver)/z0)
        1140             :
        1141             :        ! we need them to the 3.41 power, according to Gong et al., 1997:
        1142           0 :        u10cubed(:ncol)=u10cubed(:ncol)**3.41_r8
        1143             :
        1144           0 :        sflx(:)=0._r8
        1145             :
        1146           0 :        call seasalt_emis( u10cubed, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx )
        1147             :
        1148           0 :        do m=1,seasalt_nbin
        1149           0 :           mm = seasalt_indices(m)
        1150           0 :           sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
        1151           0 :           call outfld(trim(seasalt_names(m))//'SF',cam_in%cflx(:,mm),pcols,lchnk)
        1152             :        enddo
        1153           0 :        call outfld('SSTSFMBL',sflx(:),pcols,lchnk)
        1154             :     endif
        1155             :
        1156           0 :   end subroutine aero_model_emissions
        1157             :
        1158             : end module aero_model

Generated by: LCOV version 1.14