LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_gas_phase_chemdr.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 363 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 2 | 0.0 %

              Line data    Source code
           1             : module mo_gas_phase_chemdr
           2             :
           3             :   use shr_kind_mod,     only : r8 => shr_kind_r8
           4             :   use shr_const_mod,    only : pi => shr_const_pi
           5             :   use constituents,     only : pcnst
           6             :   use cam_history,      only : fieldname_len
           7             :   use chem_mods,        only : phtcnt, rxntot, gas_pcnst
           8             :   use chem_mods,        only : rxt_tag_cnt, rxt_tag_lst, rxt_tag_map, extcnt, num_rnts
           9             :   use ppgrid,           only : pcols, pver
          10             :   use phys_control,     only : phys_getopts
          11             :   use carma_flags_mod,  only : carma_hetchem_feedback
          12             :   use chem_prod_loss_diags, only: chem_prod_loss_diags_init, chem_prod_loss_diags_out
          13             :
          14             :   implicit none
          15             :   save
          16             :
          17             :   private
          18             :   public :: gas_phase_chemdr, gas_phase_chemdr_inti
          19             :   public :: map2chm
          20             :
          21             :   integer :: map2chm(pcnst) = 0           ! index map to/from chemistry/constituents list
          22             :
          23             :   integer :: so4_ndx, h2o_ndx, o2_ndx, o_ndx, hno3_ndx, hcl_ndx, cldice_ndx, snow_ndx
          24             :   integer :: o3_ndx, o3s_ndx, o3_inv_ndx, srf_ozone_pbf_ndx=-1
          25             :   integer :: het1_ndx
          26             :   integer :: ndx_cldfr, ndx_cmfdqr, ndx_nevapr, ndx_cldtop, ndx_prain
          27             :   integer :: ndx_h2so4
          28             : !
          29             : ! CCMI
          30             : !
          31             :   integer :: st80_25_ndx
          32             :   integer :: st80_25_tau_ndx
          33             :   integer :: aoa_nh_ndx
          34             :   integer :: aoa_nh_ext_ndx
          35             :   integer :: nh_5_ndx
          36             :   integer :: nh_50_ndx
          37             :   integer :: nh_50w_ndx
          38             :   integer :: sad_pbf_ndx
          39             :   integer :: cb1_ndx,cb2_ndx,oc1_ndx,oc2_ndx,dst1_ndx,dst2_ndx,sslt1_ndx,sslt2_ndx
          40             :   integer :: soa_ndx,soai_ndx,soam_ndx,soat_ndx,soab_ndx,soax_ndx
          41             :
          42             :   character(len=fieldname_len),dimension(rxt_tag_cnt)   :: tag_names
          43             :   character(len=fieldname_len),dimension(extcnt)        :: extfrc_name
          44             :
          45             :   logical :: pm25_srf_diag
          46             :   logical :: pm25_srf_diag_soa
          47             :
          48             :   logical :: convproc_do_aer
          49             :   integer :: ele_temp_ndx, ion_temp_ndx
          50             :
          51             :   ! for HEMCO-CESM ... passing J-values to ParaNOx ship plume extension
          52             :   integer :: hco_jno2_idx, hco_joh_idx
          53             :   integer :: rxt_jno2_idx, rxt_joh_idx
          54             :
          55             : contains
          56             :
          57           0 :   subroutine gas_phase_chemdr_inti()
          58             :
          59             :     use mo_chem_utls,      only : get_spc_ndx, get_extfrc_ndx, get_rxt_ndx, get_inv_ndx
          60             :     use cam_history,       only : addfld,add_default,horiz_only
          61             :     use mo_chm_diags,      only : chm_diags_inti
          62             :     use constituents,      only : cnst_get_ind
          63             :     use physics_buffer,    only : pbuf_get_index
          64             :     use rate_diags,        only : rate_diags_init
          65             :     use cam_abortutils,    only : endrun
          66             :
          67             :     character(len=3) :: string
          68             :     integer          :: n, m, err, ii
          69             :     logical :: history_cesm_forcing
          70             :     character(len=16) :: unitstr
          71             :     !-----------------------------------------------------------------------
          72             :     logical :: history_scwaccm_forcing
          73             :
          74           0 :     call phys_getopts( history_scwaccm_forcing_out = history_scwaccm_forcing )
          75             :
          76           0 :     call phys_getopts( convproc_do_aer_out = convproc_do_aer, history_cesm_forcing_out=history_cesm_forcing )
          77             :
          78           0 :     ndx_h2so4 = get_spc_ndx('H2SO4')
          79             : !
          80             : ! CCMI
          81             : !
          82           0 :     st80_25_ndx     = get_spc_ndx   ('ST80_25')
          83           0 :     st80_25_tau_ndx = get_rxt_ndx   ('ST80_25_tau')
          84           0 :     aoa_nh_ndx      = get_spc_ndx   ('AOA_NH')
          85           0 :     aoa_nh_ext_ndx  = get_extfrc_ndx('AOA_NH')
          86           0 :     nh_5_ndx        = get_spc_ndx('NH_5')
          87           0 :     nh_50_ndx       = get_spc_ndx('NH_50')
          88           0 :     nh_50w_ndx      = get_spc_ndx('NH_50W')
          89             : !
          90           0 :     cb1_ndx         = get_spc_ndx('CB1')
          91           0 :     cb2_ndx         = get_spc_ndx('CB2')
          92           0 :     oc1_ndx         = get_spc_ndx('OC1')
          93           0 :     oc2_ndx         = get_spc_ndx('OC2')
          94           0 :     dst1_ndx        = get_spc_ndx('DST01')
          95           0 :     dst2_ndx        = get_spc_ndx('DST02')
          96           0 :     sslt1_ndx       = get_spc_ndx('SSLT01')
          97           0 :     sslt2_ndx       = get_spc_ndx('SSLT02')
          98           0 :     soa_ndx         = get_spc_ndx('SOA')
          99           0 :     soam_ndx         = get_spc_ndx('SOAM')
         100           0 :     soai_ndx         = get_spc_ndx('SOAI')
         101           0 :     soat_ndx         = get_spc_ndx('SOAT')
         102           0 :     soab_ndx         = get_spc_ndx('SOAB')
         103           0 :     soax_ndx         = get_spc_ndx('SOAX')
         104             :
         105             :     pm25_srf_diag = cb1_ndx>0 .and. cb2_ndx>0 .and. oc1_ndx>0 .and. oc2_ndx>0 &
         106             :               .and. dst1_ndx>0 .and. dst2_ndx>0 .and. sslt1_ndx>0 .and. sslt2_ndx>0 &
         107           0 :               .and. soa_ndx>0
         108             :
         109             :     pm25_srf_diag_soa = cb1_ndx>0 .and. cb2_ndx>0 .and. oc1_ndx>0 .and. oc2_ndx>0 &
         110             :               .and. dst1_ndx>0 .and. dst2_ndx>0 .and. sslt1_ndx>0 .and. sslt2_ndx>0 &
         111           0 :               .and. soam_ndx>0 .and. soai_ndx>0 .and. soat_ndx>0 .and. soab_ndx>0 .and. soax_ndx>0
         112             :
         113           0 :     if ( pm25_srf_diag .or. pm25_srf_diag_soa) then
         114           0 :        call addfld('PM25_SRF',horiz_only,'I','kg/kg','bottom layer PM2.5 mixing ratio' )
         115             :     endif
         116           0 :     call addfld('U_SRF',horiz_only,'I','m/s','bottom layer wind velocity' )
         117           0 :     call addfld('V_SRF',horiz_only,'I','m/s','bottom layer wind velocity' )
         118           0 :     call addfld('Q_SRF',horiz_only,'I','kg/kg','bottom layer specific humidity' )
         119             : !
         120           0 :     het1_ndx= get_rxt_ndx('het1')
         121           0 :     o3_ndx  = get_spc_ndx('O3')
         122           0 :     o3_inv_ndx  = get_inv_ndx( 'O3' )
         123           0 :     o3s_ndx = get_spc_ndx('O3S')
         124           0 :     o_ndx   = get_spc_ndx('O')
         125           0 :     o2_ndx  = get_spc_ndx('O2')
         126           0 :     so4_ndx = get_spc_ndx('SO4')
         127           0 :     h2o_ndx = get_spc_ndx('H2O')
         128           0 :     hno3_ndx = get_spc_ndx('HNO3')
         129           0 :     hcl_ndx  = get_spc_ndx('HCL')
         130           0 :     call cnst_get_ind( 'CLDICE', cldice_ndx )
         131           0 :     call cnst_get_ind( 'SNOWQM', snow_ndx, abort=.false. )
         132             :
         133           0 :     if (o3_ndx>0 .or. o3_inv_ndx>0) then
         134           0 :        srf_ozone_pbf_ndx = pbuf_get_index('SRFOZONE')
         135             :     endif
         136             :
         137           0 :     do m = 1,extcnt
         138             :        WRITE(UNIT=string, FMT='(I2.2)') m
         139             :        extfrc_name(m) = 'extfrc_'// trim(string)
         140             :        call addfld( extfrc_name(m), (/ 'lev' /), 'I', ' ', 'ext frcing' )
         141             :     end do
         142             :
         143             :     do n = 1,rxt_tag_cnt
         144             :        tag_names(n) = trim(rxt_tag_lst(n))
         145             :        if (n<=phtcnt) then
         146             :           call addfld( tag_names(n), (/ 'lev' /), 'I', '/s', 'photolysis rate constant' )
         147             :        else
         148             :           ii = n-phtcnt
         149             :           select case(num_rnts(ii))
         150             :           case(1)
         151             :              unitstr='/s'
         152             :           case(2)
         153             :              unitstr='cm3/molecules/s'
         154             :           case(3)
         155             :              unitstr='cm6/molecules2/s'
         156             :           case default
         157             :              call endrun('gas_phase_chemdr_inti: invalid value in num_rnts used to set units in reaction rate constant')
         158             :           end select
         159             :           call addfld( tag_names(n), (/ 'lev' /), 'I', unitstr, 'reaction rate constant' )
         160             :        endif
         161             :        if (history_scwaccm_forcing) then
         162             :           select case (trim(tag_names(n)))
         163             :           case ('jh2o_a', 'jh2o_b', 'jh2o_c' )
         164             :              call add_default( tag_names(n), 1, ' ')
         165             :           end select
         166             :        endif
         167             :     enddo
         168             :
         169           0 :     call addfld( 'QDSAD',      (/ 'lev' /), 'I', '/s',      'water vapor sad delta' )
         170           0 :     call addfld( 'SAD_STRAT',  (/ 'lev' /), 'I', 'cm2/cm3', 'stratospheric aerosol SAD' )
         171           0 :     call addfld( 'SAD_SULFC',  (/ 'lev' /), 'I', 'cm2/cm3', 'chemical sulfate aerosol SAD' )
         172           0 :     call addfld( 'SAD_SAGE',   (/ 'lev' /), 'I', 'cm2/cm3', 'SAGE sulfate aerosol SAD' )
         173           0 :     call addfld( 'SAD_LNAT',   (/ 'lev' /), 'I', 'cm2/cm3', 'large-mode NAT aerosol SAD' )
         174           0 :     call addfld( 'SAD_ICE',    (/ 'lev' /), 'I', 'cm2/cm3', 'water-ice aerosol SAD' )
         175           0 :     call addfld( 'RAD_SULFC',  (/ 'lev' /), 'I', 'cm',      'chemical sad sulfate' )
         176           0 :     call addfld( 'RAD_LNAT',   (/ 'lev' /), 'I', 'cm',      'large nat radius' )
         177           0 :     call addfld( 'RAD_ICE',    (/ 'lev' /), 'I', 'cm',      'sad ice' )
         178           0 :     call addfld( 'SAD_TROP',   (/ 'lev' /), 'I', 'cm2/cm3', 'tropospheric aerosol SAD' )
         179           0 :     call addfld( 'SAD_AERO',   (/ 'lev' /), 'I', 'cm2/cm3', 'aerosol surface area density' )
         180           0 :     if (history_cesm_forcing) then
         181           0 :        call add_default ('SAD_AERO',8,' ')
         182             :     endif
         183           0 :     call addfld( 'REFF_AERO',  (/ 'lev' /), 'I', 'cm',      'aerosol effective radius' )
         184           0 :     call addfld( 'SULF_TROP',  (/ 'lev' /), 'I', 'mol/mol', 'tropospheric aerosol SAD' )
         185           0 :     call addfld( 'QDSETT',     (/ 'lev' /), 'I', '/s',      'water vapor settling delta' )
         186           0 :     call addfld( 'QDCHEM',     (/ 'lev' /), 'I', '/s',      'water vapor chemistry delta')
         187           0 :     call addfld( 'HNO3_TOTAL', (/ 'lev' /), 'I', 'mol/mol', 'total HNO3' )
         188           0 :     call addfld( 'HNO3_STS',   (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HNO3' )
         189           0 :     call addfld( 'HNO3_NAT',   (/ 'lev' /), 'I', 'mol/mol', 'NAT condensed HNO3' )
         190           0 :     call addfld( 'HNO3_GAS',   (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hno3' )
         191           0 :     call addfld( 'H2O_GAS',    (/ 'lev' /), 'I', 'mol/mol', 'gas-phase h2o' )
         192           0 :     call addfld( 'HCL_TOTAL',  (/ 'lev' /), 'I', 'mol/mol', 'total hcl' )
         193           0 :     call addfld( 'HCL_GAS',    (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hcl' )
         194           0 :     call addfld( 'HCL_STS',    (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HCL' )
         195             :
         196           0 :     if (het1_ndx>0) then
         197           0 :        call addfld( 'het1_total', (/ 'lev' /), 'I', '/s', 'total N2O5 + H2O het rate constant' )
         198             :     endif
         199           0 :     call addfld( 'SZA', horiz_only, 'I', 'degrees', 'solar zenith angle' )
         200             :
         201           0 :     call chm_diags_inti()
         202           0 :     call rate_diags_init()
         203             :
         204             : !-----------------------------------------------------------------------
         205             : ! get pbuf indicies
         206             : !-----------------------------------------------------------------------
         207           0 :     ndx_cldfr  = pbuf_get_index('CLD')
         208           0 :     ndx_cmfdqr = pbuf_get_index('RPRDTOT')
         209           0 :     ndx_nevapr = pbuf_get_index('NEVAPR')
         210           0 :     ndx_prain  = pbuf_get_index('PRAIN')
         211           0 :     ndx_cldtop = pbuf_get_index('CLDTOP')
         212             :
         213           0 :     sad_pbf_ndx= pbuf_get_index('VOLC_SAD',errcode=err) ! prescribed  strat aerosols (volcanic)
         214           0 :     if (.not.sad_pbf_ndx>0) sad_pbf_ndx = pbuf_get_index('SADSULF',errcode=err) ! CARMA's version of strat aerosols
         215             :
         216           0 :     ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index
         217           0 :     ion_temp_ndx = pbuf_get_index('TIon',errcode=err) ! ion temperature index
         218             :
         219             :     ! diagnostics for stratospheric heterogeneous reactions
         220           0 :     call addfld( 'GAMMA_HET1', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         221           0 :     call addfld( 'GAMMA_HET2', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         222           0 :     call addfld( 'GAMMA_HET3', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         223           0 :     call addfld( 'GAMMA_HET4', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         224           0 :     call addfld( 'GAMMA_HET5', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         225           0 :     call addfld( 'GAMMA_HET6', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         226           0 :     call addfld( 'WTPER',      (/ 'lev' /), 'I', '%', 'H2SO4 Weight Percent' )
         227             :
         228           0 :     call addfld( 'O3S_LOSS',      (/ 'lev' /), 'I', '1/sec', 'O3S loss rate const' )
         229             :
         230           0 :     call chem_prod_loss_diags_init
         231             :
         232             :     ! diagnostics for HEMCO ParaNOx extension
         233           0 :     hco_jno2_idx = pbuf_get_index('HCO_IN_JNO2',errcode=err)
         234           0 :     hco_joh_idx  = pbuf_get_index('HCO_IN_JOH',errcode=err)
         235             :
         236             :     !-------------------------- HEMCO_CESM ---------------------------------
         237             :     !  ... save photo rxn rates for HEMCO ParaNOx, chem_mech rxns:
         238             :     !    jo3_b            (  8)   O3 + hv ->  O + O2
         239             :     !    jno2             ( 16)   NO2 + hv ->  NO + O
         240             :     !
         241             :     ! The reactions jo2 and jo3_b exist in the mechanisms that would use
         242             :     ! the ParaNOx ship plume extension. If they do not exist, then the indices
         243             :     ! would not be found and the HCO_IN_JNO2 and HCO_IN_JOH fields would not
         244             :     ! be zero and the extension would have no effect.
         245             :     !-----------------------------------------------------------------------
         246           0 :     rxt_jno2_idx  = get_rxt_ndx( 'jno2' )
         247           0 :     rxt_joh_idx   = get_rxt_ndx( 'jo3_b' )
         248             :
         249           0 :   end subroutine gas_phase_chemdr_inti
         250             :
         251             :
         252             : !-----------------------------------------------------------------------
         253             : !-----------------------------------------------------------------------
         254           0 :   subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, &
         255             :                               phis, zm, zi, calday, &
         256             :                               tfld, pmid, pdel, pint, rpdel, rpdeldry, &
         257             :                               cldw, troplev, troplevchem, &
         258             :                               ncldwtr, ufld, vfld,  &
         259             :                               delt, ps, &
         260             :                               fsds, ts, asdir, ocnfrac, icefrac, &
         261             :                               precc, precl, snowhland, ghg_chem, latmapback, &
         262             :                               drydepflx, wetdepflx, cflx, fire_sflx, fire_ztop, nhx_nitrogen_flx, noy_nitrogen_flx, &
         263             :                               use_hemco, qtend, pbuf)
         264             :
         265             :     !-----------------------------------------------------------------------
         266             :     !     ... Chem_solver advances the volumetric mixing ratio
         267             :     !         forward one time step via a combination of explicit,
         268             :     !         ebi, hov, fully implicit, and/or rodas algorithms.
         269             :     !-----------------------------------------------------------------------
         270             :
         271           0 :     use phys_control,      only : cam_physpkg_is
         272             :     use chem_mods,         only : nabscol, nfs, indexm, clscnt4
         273             :     use physconst,         only : rga, gravit
         274             :     use mo_photo,          only : set_ub_col, setcol, table_photo
         275             :     use mo_exp_sol,        only : exp_sol
         276             :     use mo_imp_sol,        only : imp_sol
         277             :     use mo_setrxt,         only : setrxt
         278             :     use mo_adjrxt,         only : adjrxt
         279             :     use mo_phtadj,         only : phtadj
         280             :     use mo_usrrxt,         only : usrrxt
         281             :     use mo_setinv,         only : setinv
         282             :     use mo_negtrc,         only : negtrc
         283             :     use mo_sulf,           only : sulf_interp
         284             :     use mo_setext,         only : setext
         285             :     use fire_emissions,    only : fire_emissions_vrt
         286             :     use mo_sethet,         only : sethet
         287             :     use mo_drydep,         only : drydep
         288             :     use mo_fstrat,         only : set_fstrat_vals, set_fstrat_h2o
         289             :     use mo_flbc,           only : flbc_set
         290             :     use phys_grid,         only : get_rlat_all_p, get_rlon_all_p
         291             :     use mo_mean_mass,      only : set_mean_mass
         292             :     use cam_history,       only : outfld
         293             :     use wv_saturation,     only : qsat
         294             :     use constituents,      only : cnst_mw, cnst_type
         295             :     use mo_ghg_chem,       only : ghg_chem_set_rates, ghg_chem_set_flbc
         296             :     use mo_sad,            only : sad_strat_calc
         297             :     use charge_neutrality, only : charge_balance
         298             :     use mo_strato_rates,   only : ratecon_sfstrat
         299             :     use mo_aero_settling,  only : strat_aer_settling
         300             :     use shr_orb_mod,       only : shr_orb_decl
         301             :     use cam_control_mod,   only : lambm0, eccen, mvelpp, obliqr
         302             :     use mo_strato_rates,   only : has_strato_chem
         303             :     use short_lived_species,only: set_short_lived_species,get_short_lived_species
         304             :     use mo_chm_diags,      only : chm_diags, het_diags
         305             :     use perf_mod,          only : t_startf, t_stopf
         306             :     use gas_wetdep_opts,   only : gas_wetdep_method
         307             :     use physics_buffer,    only : physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx
         308             :     use infnan,            only : nan, assignment(=)
         309             :     use rate_diags,        only : rate_diags_calc, rate_diags_o3s_loss
         310             :     use mo_mass_xforms,    only : mmr2vmr, vmr2mmr, h2o_to_vmr, mmr2vmri
         311             :     use orbit,             only : zenith
         312             :
         313             : !
         314             : ! for aqueous chemistry and aerosol growth
         315             : !
         316             :     use aero_model,        only : aero_model_gasaerexch
         317             :
         318             :     use aero_model,        only : aero_model_strat_surfarea
         319             :
         320             :     !-----------------------------------------------------------------------
         321             :     !        ... Dummy arguments
         322             :     !-----------------------------------------------------------------------
         323             :     integer,        intent(in)    :: lchnk                          ! chunk index
         324             :     integer,        intent(in)    :: ncol                           ! number columns in chunk
         325             :     integer,        intent(in)    :: imozart                        ! gas phase start index in q
         326             :     real(r8),       intent(in)    :: delt                           ! timestep (s)
         327             :     real(r8),       intent(in)    :: calday                         ! day of year
         328             :     real(r8),       intent(in)    :: ps(pcols)                      ! surface pressure
         329             :     real(r8),       intent(in)    :: phis(pcols)                    ! surface geopotential
         330             :     real(r8),target,intent(in)    :: tfld(pcols,pver)               ! midpoint temperature (K)
         331             :     real(r8),       intent(in)    :: pmid(pcols,pver)               ! midpoint pressures (Pa)
         332             :     real(r8),       intent(in)    :: pdel(pcols,pver)               ! pressure delta about midpoints (Pa)
         333             :     real(r8),       intent(in)    :: rpdel(pcols,pver)              ! reciprocal pressure delta about midpoints (Pa)
         334             :     real(r8),       intent(in)    :: rpdeldry(pcols,pver)           ! reciprocal dry pressure delta about midpoints (Pa)
         335             :     real(r8),       intent(in)    :: ufld(pcols,pver)               ! zonal velocity (m/s)
         336             :     real(r8),       intent(in)    :: vfld(pcols,pver)               ! meridional velocity (m/s)
         337             :     real(r8),       intent(in)    :: cldw(pcols,pver)               ! cloud water (kg/kg)
         338             :     real(r8),       intent(in)    :: ncldwtr(pcols,pver)            ! droplet number concentration (#/kg)
         339             :     real(r8),       intent(in)    :: zm(pcols,pver)                 ! midpoint geopotential height above the surface (m)
         340             :     real(r8),       intent(in)    :: zi(pcols,pver+1)               ! interface geopotential height above the surface (m)
         341             :     real(r8),       intent(in)    :: pint(pcols,pver+1)             ! interface pressures (Pa)
         342             :     real(r8),       intent(in)    :: q(pcols,pver,pcnst)            ! species concentrations (kg/kg)
         343             :     real(r8),pointer, intent(in)  :: fire_sflx(:,:)                 ! fire emssions surface flux (kg/m^2/s)
         344             :     real(r8),pointer, intent(in)  :: fire_ztop(:)                   ! top of vertical distribution of fire emssions (m)
         345             :     real(r8),       intent(in)    :: fsds(pcols)                    ! longwave down at sfc
         346             :     real(r8),       intent(in)    :: icefrac(pcols)                 ! sea-ice areal fraction
         347             :     real(r8),       intent(in)    :: ocnfrac(pcols)                 ! ocean areal fraction
         348             :     real(r8),       intent(in)    :: asdir(pcols)                   ! albedo: shortwave, direct
         349             :     real(r8),       intent(in)    :: ts(pcols)                      ! sfc temp (merged w/ocean if coupled)
         350             :     real(r8),       intent(in)    :: precc(pcols)                   !
         351             :     real(r8),       intent(in)    :: precl(pcols)                   !
         352             :     real(r8),       intent(in)    :: snowhland(pcols)               !
         353             :     logical,        intent(in)    :: ghg_chem
         354             :     integer,        intent(in)    :: latmapback(pcols)
         355             :     integer,        intent(in)    :: troplev(pcols)                 ! trop/strat separation vertical index
         356             :     integer,        intent(in)    :: troplevchem(pcols)             ! trop/strat chemistry separation vertical index
         357             :     real(r8),       intent(inout) :: qtend(pcols,pver,pcnst)        ! species tendencies (kg/kg/s)
         358             :     real(r8),       intent(inout) :: cflx(pcols,pcnst)              ! constituent surface flux (kg/m^2/s)
         359             :     real(r8),       intent(out)   :: drydepflx(pcols,pcnst)         ! dry deposition flux (kg/m^2/s)
         360             :     real(r8),       intent(in)    :: wetdepflx(pcols,pcnst)         ! wet deposition flux (kg/m^2/s)
         361             :     real(r8), intent(out) :: nhx_nitrogen_flx(pcols)
         362             :     real(r8), intent(out) :: noy_nitrogen_flx(pcols)
         363             :     logical,        intent(in)    :: use_hemco                      ! use Harmonized Emissions Component (HEMCO)
         364             :
         365             :     type(physics_buffer_desc), pointer :: pbuf(:)
         366             :
         367             :     !-----------------------------------------------------------------------
         368             :     !           ... Local variables
         369             :     !-----------------------------------------------------------------------
         370             :     real(r8), parameter :: m2km  = 1.e-3_r8
         371             :     real(r8), parameter :: Pa2mb = 1.e-2_r8
         372             :
         373           0 :     real(r8),       pointer    :: prain(:,:)
         374           0 :     real(r8),       pointer    :: nevapr(:,:)
         375           0 :     real(r8),       pointer    :: cmfdqr(:,:)
         376           0 :     real(r8),       pointer    :: cldfr(:,:)
         377           0 :     real(r8),       pointer    :: cldtop(:)
         378             :
         379             :     integer      ::  i, k, m, n
         380             :     integer      ::  tim_ndx
         381             :     real(r8)     ::  delt_inverse
         382             :     real(r8)     ::  esfact
         383           0 :     real(r8)     ::  invariants(ncol,pver,nfs)
         384           0 :     real(r8)     ::  col_dens(ncol,pver,max(1,nabscol))    ! column densities (molecules/cm^2)
         385           0 :     real(r8)     ::  col_delta(ncol,0:pver,max(1,nabscol)) ! layer column densities (molecules/cm^2)
         386           0 :     real(r8)     ::  extfrc(ncol,pver,max(1,extcnt))
         387           0 :     real(r8)     ::  vmr(ncol,pver,gas_pcnst)         ! xported species (vmr)
         388           0 :     real(r8)     ::  reaction_rates(ncol,pver,max(1,rxntot))      ! reaction rates
         389           0 :     real(r8)     ::  depvel(ncol,gas_pcnst)                ! dry deposition velocity (cm/s)
         390           0 :     real(r8)     ::  het_rates(ncol,pver,max(1,gas_pcnst)) ! washout rate (1/s)
         391             :     real(r8), dimension(ncol,pver) :: &
         392           0 :          h2ovmr, &                                         ! water vapor volume mixing ratio
         393           0 :          mbar, &                                           ! mean wet atmospheric mass ( amu )
         394           0 :          zmid, &                                           ! midpoint geopotential in km
         395           0 :          sulfate, &                                        ! trop sulfate aerosols
         396           0 :          pmb                                               ! pressure at midpoints ( hPa )
         397             :     real(r8), dimension(ncol,pver) :: &
         398           0 :          cwat, &                                           ! cloud water mass mixing ratio (kg/kg)
         399           0 :          wrk
         400             :     real(r8), dimension(ncol,pver+1) :: &
         401           0 :          zintr                                              ! interface geopotential in km realitive to surf
         402             :     real(r8), dimension(ncol,pver+1) :: &
         403           0 :          zint                                              ! interface geopotential in km
         404             :     real(r8), dimension(ncol)  :: &
         405           0 :          zen_angle, &                                      ! solar zenith angles
         406           0 :          zsurf, &                                          ! surface height (m)
         407           0 :          rlats, rlons                                      ! chunk latitudes and longitudes (radians)
         408           0 :     real(r8) :: sza(ncol)                                  ! solar zenith angles (degrees)
         409             :     real(r8), parameter :: rad2deg = 180._r8/pi                ! radians to degrees conversion factor
         410           0 :     real(r8) :: relhum(ncol,pver)                          ! relative humidity
         411           0 :     real(r8) :: satv(ncol,pver)                            ! wrk array for relative humidity
         412           0 :     real(r8) :: satq(ncol,pver)                            ! wrk array for relative humidity
         413             :
         414             :     integer                   :: j
         415             :     integer                   ::  ltrop_sol(pcols)         ! tropopause vertical index used in chem solvers
         416           0 :     real(r8), pointer         ::  strato_sad(:,:)          ! stratospheric sad (1/cm)
         417             :
         418             :     real(r8)                  ::  sad_trop(pcols,pver)    ! total tropospheric sad (cm^2/cm^3)
         419             :     real(r8)                  ::  reff(pcols,pver)        ! aerosol effective radius (cm)
         420             :     real(r8)                  ::  reff_strat(pcols,pver)  ! stratospheric aerosol effective radius (cm)
         421             :
         422             :     real(r8) :: tvs(pcols)
         423             :     real(r8) :: wind_speed(pcols)        ! surface wind speed (m/s)
         424             :     real(r8) :: prect(pcols)
         425             :     real(r8) :: sflx(pcols,gas_pcnst)
         426             :     real(r8) :: wetdepflx_diag(pcols,gas_pcnst)
         427           0 :     real(r8) :: o2mmr(ncol,pver)               ! o2 concentration (kg/kg)
         428           0 :     real(r8) :: ommr(ncol,pver)                ! o concentration (kg/kg)
         429             :     real(r8) :: mmr(pcols,pver,gas_pcnst)      ! chem working concentrations (kg/kg)
         430             :     real(r8) :: mmr_new(pcols,pver,gas_pcnst)      ! chem working concentrations (kg/kg)
         431           0 :     real(r8) :: hno3_gas(ncol,pver)            ! hno3 gas phase concentration (mol/mol)
         432           0 :     real(r8) :: hno3_cond(ncol,pver,2)         ! hno3 condensed phase concentration (mol/mol)
         433           0 :     real(r8) :: hcl_gas(ncol,pver)             ! hcl gas phase concentration (mol/mol)
         434           0 :     real(r8) :: hcl_cond(ncol,pver)            ! hcl condensed phase concentration (mol/mol)
         435           0 :     real(r8) :: h2o_gas(ncol,pver)             ! h2o gas phase concentration (mol/mol)
         436           0 :     real(r8) :: h2o_cond(ncol,pver)            ! h2o condensed phase concentration (mol/mol)
         437             :     real(r8) :: cldice(pcols,pver)             ! cloud water "ice" (kg/kg)
         438           0 :     real(r8) :: radius_strat(ncol,pver,3)      ! radius of sulfate, nat, & ice ( cm )
         439           0 :     real(r8) :: sad_strat(ncol,pver,3)         ! surf area density of sulfate, nat, & ice ( cm^2/cm^3 )
         440             :     real(r8) :: mmr_tend(pcols,pver,gas_pcnst) ! chemistry species tendencies (kg/kg/s)
         441             :     real(r8) :: qh2o(pcols,pver)               ! specific humidity (kg/kg)
         442             :     real(r8) :: delta
         443             :
         444             :   ! for aerosol formation....
         445           0 :     real(r8) :: del_h2so4_gasprod(ncol,pver)
         446           0 :     real(r8) :: vmr0(ncol,pver,gas_pcnst)
         447             :
         448             :   ! for HEMCO-CESM ... passing J-values to ParaNOx ship plume extension
         449           0 :     real(r8), pointer :: hco_j_tmp_fld(:)    ! J-value pointer (sfc only) [1/s]
         450             :
         451             : !
         452             : ! CCMI
         453             : !
         454             :     real(r8) :: xlat
         455           0 :     real(r8) :: pm25(ncol)
         456             :
         457           0 :     real(r8) :: dlats(ncol)
         458             :
         459             :     real(r8), dimension(ncol,pver) :: &      ! aerosol reaction diagnostics
         460           0 :         gprob_n2o5, &
         461           0 :         gprob_cnt_hcl, &
         462           0 :         gprob_cnt_h2o, &
         463           0 :         gprob_bnt_h2o, &
         464           0 :         gprob_hocl_hcl, &
         465           0 :         gprob_hobr_hcl, &
         466           0 :         wtper
         467             :
         468           0 :     real(r8), pointer :: ele_temp_fld(:,:) ! electron temperature pointer
         469           0 :     real(r8), pointer :: ion_temp_fld(:,:) ! ion temperature pointer
         470           0 :     real(r8) :: prod_out(ncol,pver,max(1,clscnt4))
         471           0 :     real(r8) :: loss_out(ncol,pver,max(1,clscnt4))
         472             :
         473           0 :     real(r8) :: o3s_loss(ncol,pver)
         474           0 :     real(r8), pointer :: srf_ozone_fld(:)
         475             :
         476           0 :     if ( ele_temp_ndx>0 .and. ion_temp_ndx>0 ) then
         477           0 :        call pbuf_get_field(pbuf, ele_temp_ndx, ele_temp_fld)
         478           0 :        call pbuf_get_field(pbuf, ion_temp_ndx, ion_temp_fld)
         479             :     else
         480           0 :        ele_temp_fld => tfld
         481           0 :        ion_temp_fld => tfld
         482             :     endif
         483             :
         484             :     ! initialize to NaN to hopefully catch user defined rxts that go unset
         485           0 :     reaction_rates(:,:,:) = nan
         486             :
         487           0 :     delt_inverse = 1._r8 / delt
         488             :     !-----------------------------------------------------------------------
         489             :     !        ... Get chunck latitudes and longitudes
         490             :     !-----------------------------------------------------------------------
         491           0 :     call get_rlat_all_p( lchnk, ncol, rlats )
         492           0 :     call get_rlon_all_p( lchnk, ncol, rlons )
         493           0 :     tim_ndx = pbuf_old_tim_idx()
         494           0 :     call pbuf_get_field(pbuf, ndx_prain,      prain,  start=(/1,1/), kount=(/ncol,pver/))
         495           0 :     call pbuf_get_field(pbuf, ndx_cldfr,        cldfr, start=(/1,1,tim_ndx/), kount=(/ncol,pver,1/) )
         496           0 :     call pbuf_get_field(pbuf, ndx_cmfdqr,     cmfdqr, start=(/1,1/), kount=(/ncol,pver/))
         497           0 :     call pbuf_get_field(pbuf, ndx_nevapr,     nevapr, start=(/1,1/), kount=(/ncol,pver/))
         498           0 :     call pbuf_get_field(pbuf, ndx_cldtop,     cldtop )
         499             :
         500           0 :     reff_strat(:,:) = 0._r8
         501             :
         502           0 :     dlats(:) = rlats(:)*rad2deg ! convert to degrees
         503             :
         504             :     !-----------------------------------------------------------------------
         505             :     !        ... Calculate cosine of zenith angle
         506             :     !            then cast back to angle (radians)
         507             :     !-----------------------------------------------------------------------
         508           0 :     call zenith( calday, rlats, rlons, zen_angle, ncol )
         509           0 :     zen_angle(:) = acos( zen_angle(:) )
         510             :
         511           0 :     sza(:) = zen_angle(:) * rad2deg
         512           0 :     call outfld( 'SZA',   sza,    ncol, lchnk )
         513             :
         514             :     !-----------------------------------------------------------------------
         515             :     !        ... Xform geopotential height from m to km
         516             :     !            and pressure from Pa to mb
         517             :     !-----------------------------------------------------------------------
         518           0 :     zsurf(:ncol) = rga * phis(:ncol)
         519           0 :     do k = 1,pver
         520           0 :        zintr(:ncol,k) = m2km * zi(:ncol,k)
         521           0 :        zmid(:ncol,k) = m2km * (zm(:ncol,k) + zsurf(:ncol))
         522           0 :        zint(:ncol,k) = m2km * (zi(:ncol,k) + zsurf(:ncol))
         523           0 :        pmb(:ncol,k)  = Pa2mb * pmid(:ncol,k)
         524             :     end do
         525           0 :     zint(:ncol,pver+1) = m2km * (zi(:ncol,pver+1) + zsurf(:ncol))
         526           0 :     zintr(:ncol,pver+1)= m2km *  zi(:ncol,pver+1)
         527             :
         528             :     !-----------------------------------------------------------------------
         529             :     !        ... map incoming concentrations to working array
         530             :     !-----------------------------------------------------------------------
         531           0 :     do m = 1,pcnst
         532           0 :        n = map2chm(m)
         533           0 :        if( n > 0 ) then
         534           0 :           mmr(:ncol,:,n) = q(:ncol,:,m)
         535             :        end if
         536             :     end do
         537             :
         538           0 :     call get_short_lived_species( mmr, lchnk, ncol, pbuf )
         539             :
         540             :     !-----------------------------------------------------------------------
         541             :     !        ... Set atmosphere mean mass
         542             :     !-----------------------------------------------------------------------
         543           0 :     call set_mean_mass( ncol, mmr, mbar )
         544             :
         545             :     !-----------------------------------------------------------------------
         546             :     !        ... Xform from mmr to vmr
         547             :     !-----------------------------------------------------------------------
         548           0 :     call mmr2vmr( mmr(:ncol,:,:), vmr(:ncol,:,:), mbar(:ncol,:), ncol )
         549             :
         550             : !
         551             : ! CCMI
         552             : !
         553             : ! reset STE tracer to specific vmr of 200 ppbv
         554             : !
         555           0 :     if ( st80_25_ndx > 0 ) then
         556           0 :        where ( pmid(:ncol,:) < 80.e+2_r8 )
         557           0 :           vmr(:ncol,:,st80_25_ndx) = 200.e-9_r8
         558             :        end where
         559             :     end if
         560             : !
         561             : ! reset AOA_NH, NH_5, NH_50, NH_50W surface mixing ratios between 30N and 50N
         562             : !
         563           0 :     if ( aoa_nh_ndx>0 ) then
         564           0 :       do j=1,ncol
         565           0 :         xlat = dlats(j)
         566           0 :         if ( xlat >= 30._r8 .and. xlat <= 50._r8 ) then
         567           0 :            vmr(j,pver,aoa_nh_ndx) = 0._r8
         568             :         end if
         569             :       end do
         570             :     end if
         571           0 :     if ( nh_5_ndx>0 ) then
         572           0 :       do j=1,ncol
         573           0 :         xlat = dlats(j)
         574           0 :         if ( xlat >= 30._r8 .and. xlat <= 50._r8 ) then
         575           0 :            vmr(j,pver,nh_5_ndx)   = 100.e-9_r8
         576             :         end if
         577             :       end do
         578             :     end if
         579           0 :     if ( nh_50_ndx>0 ) then
         580           0 :       do j=1,ncol
         581           0 :         xlat = dlats(j)
         582           0 :         if ( xlat >= 30._r8 .and. xlat <= 50._r8 ) then
         583           0 :            vmr(j,pver,nh_50_ndx)  = 100.e-9_r8
         584             :         end if
         585             :       end do
         586             :     end if
         587           0 :     if ( nh_50w_ndx>0 ) then
         588           0 :       do j=1,ncol
         589           0 :         xlat = dlats(j)
         590           0 :         if ( xlat >= 30._r8 .and. xlat <= 50._r8 ) then
         591           0 :            vmr(j,pver,nh_50w_ndx) = 100.e-9_r8
         592             :         end if
         593             :       end do
         594             :     end if
         595             :
         596           0 :     if (h2o_ndx>0) then
         597             :        !-----------------------------------------------------------------------
         598             :        !        ... store water vapor in wrk variable
         599             :        !-----------------------------------------------------------------------
         600           0 :        qh2o(:ncol,:) = mmr(:ncol,:,h2o_ndx)
         601             :        h2ovmr(:ncol,:) = vmr(:ncol,:,h2o_ndx)
         602             :     else
         603           0 :        qh2o(:ncol,:) = q(:ncol,:,1)
         604             :        !-----------------------------------------------------------------------
         605             :        !        ... Xform water vapor from mmr to vmr and set upper bndy values
         606             :        !-----------------------------------------------------------------------
         607           0 :        call h2o_to_vmr( q(:ncol,:,1), h2ovmr(:ncol,:), mbar(:ncol,:), ncol )
         608             :
         609           0 :        call set_fstrat_h2o( h2ovmr, pmid, troplev, calday, ncol, lchnk )
         610             :
         611             :     endif
         612             :
         613             :     !-----------------------------------------------------------------------
         614             :     !        ... force ion/electron balance
         615             :     !-----------------------------------------------------------------------
         616           0 :     call charge_balance( ncol, vmr )
         617             :
         618             :     !-----------------------------------------------------------------------
         619             :     !        ... Set the "invariants"
         620             :     !-----------------------------------------------------------------------
         621           0 :     call setinv( invariants, tfld, h2ovmr, vmr, pmid, ncol, lchnk, pbuf )
         622             :
         623             :     !-----------------------------------------------------------------------
         624             :     !        ... stratosphere aerosol surface area
         625             :     !-----------------------------------------------------------------------
         626           0 :     if (sad_pbf_ndx>0) then
         627           0 :        call pbuf_get_field(pbuf, sad_pbf_ndx, strato_sad)
         628             :     else
         629           0 :        allocate(strato_sad(pcols,pver))
         630           0 :        strato_sad(:,:) = 0._r8
         631             :
         632             :        ! Prognostic modal stratospheric sulfate: compute dry strato_sad
         633           0 :        call aero_model_strat_surfarea( ncol, mmr, pmid, tfld, troplevchem, pbuf, strato_sad, reff_strat )
         634             :
         635             :     endif
         636             :
         637           0 :     stratochem: if ( has_strato_chem ) then
         638             :        !-----------------------------------------------------------------------
         639             :        !        ... initialize condensed and gas phases; all hno3 to gas
         640             :        !-----------------------------------------------------------------------
         641           0 :        hcl_cond(:,:)      = 0.0_r8
         642           0 :        hcl_gas (:,:)      = 0.0_r8
         643             :        do k = 1,pver
         644           0 :           hno3_gas(:,k)   = vmr(:,k,hno3_ndx)
         645             :           h2o_gas(:,k)    = h2ovmr(:,k)
         646             :           hcl_gas(:,k)    = vmr(:,k,hcl_ndx)
         647             :           wrk(:,k)        = h2ovmr(:,k)
         648             :           if (snow_ndx>0) then
         649             :              cldice(:ncol,k) = q(:ncol,k,cldice_ndx) + q(:ncol,k,snow_ndx)
         650             :           else
         651             :              cldice(:ncol,k) = q(:ncol,k,cldice_ndx)
         652             :           endif
         653             :        end do
         654             :        do m = 1,2
         655             :           do k = 1,pver
         656             :              hno3_cond(:,k,m) = 0._r8
         657             :           end do
         658             :        end do
         659             :
         660             :        call mmr2vmri( cldice(:ncol,:), h2o_cond(:ncol,:), mbar(:ncol,:), cnst_mw(cldice_ndx), ncol )
         661             :
         662             :        !-----------------------------------------------------------------------
         663             :        !        ... call SAD routine
         664             :        !-----------------------------------------------------------------------
         665             :        call sad_strat_calc( lchnk, invariants(:ncol,:,indexm), pmb, tfld, hno3_gas, &
         666             :             hno3_cond, h2o_gas, h2o_cond, hcl_gas, hcl_cond, strato_sad(:ncol,:), radius_strat, &
         667             :             sad_strat, ncol, pbuf )
         668             :
         669             : !      NOTE: output of total HNO3 is before vmr is set to gas-phase.
         670             :        call outfld( 'HNO3_TOTAL', vmr(:ncol,:,hno3_ndx), ncol ,lchnk )
         671             :
         672             :
         673             :        do k = 1,pver
         674             :           vmr(:,k,hno3_ndx) = hno3_gas(:,k)
         675             :           h2ovmr(:,k)       = h2o_gas(:,k)
         676             :           vmr(:,k,h2o_ndx)  = h2o_gas(:,k)
         677             :           wrk(:,k)          = (h2ovmr(:,k) - wrk(:,k))*delt_inverse
         678             :        end do
         679             :
         680             :        call outfld( 'QDSAD', wrk(:,:), ncol, lchnk )
         681             : !
         682             :        call outfld( 'SAD_STRAT',  strato_sad(:ncol,:), ncol, lchnk )
         683             :        call outfld( 'SAD_SULFC',  sad_strat(:,:,1), ncol, lchnk )
         684             :        call outfld( 'SAD_LNAT',   sad_strat(:,:,2), ncol, lchnk )
         685             :        call outfld( 'SAD_ICE',    sad_strat(:,:,3), ncol, lchnk )
         686             : !
         687             :        call outfld( 'RAD_SULFC',  radius_strat(:,:,1), ncol, lchnk )
         688             :        call outfld( 'RAD_LNAT',   radius_strat(:,:,2), ncol, lchnk )
         689             :        call outfld( 'RAD_ICE',    radius_strat(:,:,3), ncol, lchnk )
         690             : !
         691             :        call outfld( 'HNO3_GAS',   vmr(:ncol,:,hno3_ndx), ncol, lchnk )
         692             :        call outfld( 'HNO3_STS',   hno3_cond(:,:,1), ncol, lchnk )
         693             :        call outfld( 'HNO3_NAT',   hno3_cond(:,:,2), ncol, lchnk )
         694             : !
         695             :        call outfld( 'HCL_TOTAL',  vmr(:ncol,:,hcl_ndx), ncol, lchnk )
         696             :        call outfld( 'HCL_GAS',    hcl_gas (:,:), ncol ,lchnk )
         697             :        call outfld( 'HCL_STS',    hcl_cond(:,:), ncol ,lchnk )
         698             :
         699             :        !-----------------------------------------------------------------------
         700             :        !        ... call aerosol reaction rates
         701             :        !-----------------------------------------------------------------------
         702             :        call ratecon_sfstrat( ncol, invariants(:,:,indexm), pmid, tfld, &
         703             :             radius_strat(:,:,1), sad_strat(:,:,1), sad_strat(:,:,2), &
         704             :             sad_strat(:,:,3), h2ovmr, vmr, reaction_rates, &
         705             :             gprob_n2o5, gprob_cnt_hcl, gprob_cnt_h2o, gprob_bnt_h2o, &
         706             :             gprob_hocl_hcl, gprob_hobr_hcl, wtper )
         707             :
         708             :        call outfld( 'GAMMA_HET1', gprob_n2o5    (:ncol,:), ncol, lchnk )
         709             :        call outfld( 'GAMMA_HET2', gprob_cnt_h2o (:ncol,:), ncol, lchnk )
         710             :        call outfld( 'GAMMA_HET3', gprob_bnt_h2o (:ncol,:), ncol, lchnk )
         711             :        call outfld( 'GAMMA_HET4', gprob_cnt_hcl (:ncol,:), ncol, lchnk )
         712             :        call outfld( 'GAMMA_HET5', gprob_hocl_hcl(:ncol,:), ncol, lchnk )
         713             :        call outfld( 'GAMMA_HET6', gprob_hobr_hcl(:ncol,:), ncol, lchnk )
         714             :        call outfld( 'WTPER',      wtper         (:ncol,:), ncol, lchnk )
         715             :
         716             :     endif stratochem
         717             :
         718             : !      NOTE: For gas-phase solver only.
         719             : !            ratecon_sfstrat needs total hcl.
         720           0 :     if (hcl_ndx>0) then
         721           0 :        vmr(:,:,hcl_ndx)  = hcl_gas(:,:)
         722             :     endif
         723             :
         724             :     !-----------------------------------------------------------------------
         725             :     !        ... Set the column densities at the upper boundary
         726             :     !-----------------------------------------------------------------------
         727           0 :     call set_ub_col( col_delta, vmr, invariants, pint(:,1), pdel, ncol, lchnk)
         728             :
         729             :     !-----------------------------------------------------------------------
         730             :     !       ...  Set rates for "tabular" and user specified reactions
         731             :     !-----------------------------------------------------------------------
         732           0 :     call setrxt( reaction_rates, tfld, invariants(1,1,indexm), ncol )
         733             :
         734           0 :     sulfate(:,:) = 0._r8
         735           0 :     if ( .not. carma_hetchem_feedback ) then
         736           0 :        if( so4_ndx < 1 ) then ! get offline so4 field if not prognostic
         737           0 :           call sulf_interp( ncol, lchnk, sulfate )
         738             :        else
         739           0 :           sulfate(:,:) = vmr(:,:,so4_ndx)
         740             :        endif
         741             :     endif
         742             :
         743             :     !-----------------------------------------------------------------
         744             :     ! ... zero out sulfate above tropopause
         745             :     !-----------------------------------------------------------------
         746           0 :     do k = 1, pver
         747           0 :        do i = 1, ncol
         748           0 :           if (k < troplevchem(i)) then
         749           0 :              sulfate(i,k) = 0.0_r8
         750             :           end if
         751             :        end do
         752             :     end do
         753             :
         754           0 :     call outfld( 'SULF_TROP', sulfate(:ncol,:), ncol, lchnk )
         755             :
         756             :     !-----------------------------------------------------------------
         757             :     !   ... compute the relative humidity
         758             :     !-----------------------------------------------------------------
         759           0 :     do k = 1, pver
         760           0 :        call qsat(tfld(1:ncol,k), pmid(1:ncol,k), satv(1:ncol,k), satq(1:ncol,k), ncol)
         761             :     end do
         762             :
         763           0 :     do k = 1,pver
         764           0 :        relhum(:,k) = .622_r8 * h2ovmr(:,k) / satq(:,k)
         765           0 :        relhum(:,k) = max( 0._r8,min( 1._r8,relhum(:,k) ) )
         766             :     end do
         767             :
         768           0 :     cwat(:ncol,:pver) = cldw(:ncol,:pver)
         769             :
         770             :     call usrrxt( reaction_rates, tfld, ion_temp_fld, ele_temp_fld, invariants, h2ovmr, &
         771             :                  pmid, invariants(:,:,indexm), sulfate, mmr, relhum, strato_sad, &
         772           0 :                  troplevchem, dlats, ncol, sad_trop, reff, cwat, mbar, pbuf )
         773             :
         774           0 :     call outfld( 'SAD_TROP', sad_trop(:ncol,:), ncol, lchnk )
         775             :
         776             :     ! Add trop/strat components of SAD for output
         777           0 :     sad_trop(:ncol,:)=sad_trop(:ncol,:)+strato_sad(:ncol,:)
         778           0 :     call outfld( 'SAD_AERO', sad_trop(:ncol,:), ncol, lchnk )
         779             :
         780             :     ! Add trop/strat components of effective radius for output
         781           0 :     reff(:ncol,:)=reff(:ncol,:)+reff_strat(:ncol,:)
         782           0 :     call outfld( 'REFF_AERO', reff(:ncol,:), ncol, lchnk )
         783             :
         784           0 :     if (het1_ndx>0) then
         785           0 :        call outfld( 'het1_total', reaction_rates(:,:,het1_ndx), ncol, lchnk )
         786             :     endif
         787             :
         788           0 :     if (ghg_chem) then
         789           0 :        call ghg_chem_set_rates( reaction_rates, latmapback, zen_angle, ncol, lchnk )
         790             :     endif
         791             :
         792             :     do i = phtcnt+1,rxt_tag_cnt
         793             :        call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk )
         794             :     enddo
         795             :
         796           0 :     call adjrxt( reaction_rates, invariants, invariants(1,1,indexm), ncol,pver )
         797             :
         798             :     !-----------------------------------------------------------------------
         799             :     !        ... Compute the photolysis rates at time = t(n+1)
         800             :     !-----------------------------------------------------------------------
         801             :     !-----------------------------------------------------------------------
         802             :     !           ... Set the column densities
         803             :     !-----------------------------------------------------------------------
         804           0 :     call setcol( col_delta, col_dens )
         805             :
         806             :     !-----------------------------------------------------------------------
         807             :     !           ... Calculate the photodissociation rates
         808             :     !-----------------------------------------------------------------------
         809             :
         810           0 :     esfact = 1._r8
         811             :     call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr  , &
         812           0 :          delta, esfact )
         813             :
         814             :     !-----------------------------------------------------------------
         815             :     !   ... lookup the photolysis rates from table
         816             :     !-----------------------------------------------------------------
         817             :     call table_photo( reaction_rates, pmid, pdel, tfld, zmid, zint, &
         818             :                       col_dens, zen_angle, asdir, cwat, cldfr, &
         819           0 :                       esfact, vmr, invariants, ncol, lchnk, pbuf )
         820             :
         821             :     do i = 1,phtcnt
         822             :        call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk )
         823             :     enddo
         824             :
         825             :     !-----------------------------------------------------------------------
         826             :     !           ... Adjust the photodissociation rates
         827             :     !-----------------------------------------------------------------------
         828           0 :     call phtadj( reaction_rates, invariants, invariants(:,:,indexm), ncol,pver )
         829             :
         830           0 :     if ( use_hemco ) then
         831             :        !-------------------------- HEMCO_CESM ---------------------------------
         832             :        !  ... save photo rxn rates for HEMCO ParaNOx, chem_mech rxns:
         833             :        !    jo3_b            (  8)   O3 + hv ->  O + O2
         834             :        !    jno2             ( 16)   NO2 + hv ->  NO + O
         835             :        !-----------------------------------------------------------------------
         836             :        ! get the rxn rate [1/s] and write to pbuf
         837           0 :        if(rxt_jno2_idx > 0 .and. hco_jno2_idx > 0) then
         838           0 :          call pbuf_get_field(pbuf, hco_jno2_idx, hco_j_tmp_fld)
         839             :          ! this is already in chunk, write /pcols/ at surface
         840           0 :          hco_j_tmp_fld(:ncol) = reaction_rates(:ncol,pver,rxt_jno2_idx)
         841             :        endif
         842             :
         843           0 :        if(rxt_joh_idx > 0 .and. hco_joh_idx > 0) then
         844           0 :          call pbuf_get_field(pbuf, hco_joh_idx, hco_j_tmp_fld)
         845             :          ! this is already in chunk, write /pcols/ at surface
         846           0 :          hco_j_tmp_fld(:ncol) = reaction_rates(:ncol,pver,rxt_joh_idx)
         847             :        endif
         848             :     endif
         849             :
         850             :     !-----------------------------------------------------------------------
         851             :     !        ... Compute the extraneous frcing at time = t(n+1)
         852             :     !-----------------------------------------------------------------------
         853           0 :     if ( o2_ndx > 0 .and. o_ndx > 0 ) then
         854             :        do k = 1,pver
         855           0 :           o2mmr(:ncol,k) = mmr(:ncol,k,o2_ndx)
         856             :           ommr(:ncol,k)  = mmr(:ncol,k,o_ndx)
         857             :        end do
         858             :     endif
         859             :     call setext( extfrc, zint, zintr, cldtop, &
         860             :                  zmid, lchnk, tfld, o2mmr, ommr, &
         861           0 :                  pmid, mbar, rlats, calday, ncol, rlons, pbuf )
         862             :     ! include forcings from fire emissions ...
         863           0 :     call fire_emissions_vrt( ncol, lchnk, zint, fire_sflx, fire_ztop, extfrc )
         864             :
         865             :     do m = 1,extcnt
         866             :        if( m /= aoa_nh_ext_ndx ) then
         867             :           do k = 1,pver
         868             :              extfrc(:ncol,k,m) = extfrc(:ncol,k,m) / invariants(:ncol,k,indexm)
         869             :           end do
         870             :        endif
         871             :        call outfld( extfrc_name(m), extfrc(:ncol,:,m), ncol, lchnk )
         872             :     end do
         873             :
         874             :     !-----------------------------------------------------------------------
         875             :     !        ... Form the washout rates
         876             :     !-----------------------------------------------------------------------
         877           0 :     if ( gas_wetdep_method=='MOZ' ) then
         878             :        call sethet( het_rates, pmid, zmid, phis, tfld, &
         879             :                     cmfdqr, prain, nevapr, delt, invariants(:,:,indexm), &
         880           0 :                     vmr, ncol, lchnk )
         881           0 :        if (.not. convproc_do_aer) then
         882           0 :           call het_diags( het_rates(:ncol,:,:), mmr(:ncol,:,:), pdel(:ncol,:), lchnk, ncol )
         883             :        endif
         884             :     else
         885           0 :        het_rates = 0._r8
         886             :     end if
         887             : !
         888             : ! CCMI
         889             : !
         890             : ! set loss to below the tropopause only
         891             : !
         892           0 :     if ( st80_25_tau_ndx > 0 ) then
         893           0 :        do i = 1,ncol
         894           0 :           reaction_rates(i,1:troplev(i),st80_25_tau_ndx) = 0._r8
         895             :        enddo
         896             :     end if
         897             :
         898           0 :     ltrop_sol(:ncol) = 0 ! apply solver to all levels
         899             :
         900             :     ! save h2so4 before gas phase chem (for later new particle nucleation)
         901           0 :     if (ndx_h2so4 > 0) then
         902           0 :        del_h2so4_gasprod(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4)
         903             :     else
         904           0 :        del_h2so4_gasprod(:,:) = 0.0_r8
         905             :     endif
         906             :
         907           0 :     vmr0(:ncol,:,:) = vmr(:ncol,:,:) ! mixing ratios before chemistry changes
         908             :
         909             :     !=======================================================================
         910             :     !        ... Call the class solution algorithms
         911             :     !=======================================================================
         912             :     !-----------------------------------------------------------------------
         913             :     !   ... Solve for "Explicit" species
         914             :     !-----------------------------------------------------------------------
         915             :     call exp_sol( vmr, reaction_rates, het_rates, extfrc, delt, invariants(1,1,indexm), ncol, lchnk, ltrop_sol )
         916             :
         917             :     !-----------------------------------------------------------------------
         918             :     !   ... Solve for "Implicit" species
         919             :     !-----------------------------------------------------------------------
         920           0 :     if ( has_strato_chem ) wrk(:,:) = vmr(:,:,h2o_ndx)
         921           0 :     call t_startf('imp_sol')
         922             :     !
         923             :     call imp_sol( vmr, reaction_rates, het_rates, extfrc, delt, &
         924           0 :                   ncol,pver, lchnk,  prod_out, loss_out )
         925             :
         926           0 :     call t_stopf('imp_sol')
         927             :
         928           0 :     call chem_prod_loss_diags_out( ncol, lchnk, vmr, reaction_rates, prod_out, loss_out, invariants(:ncol,:,indexm) )
         929           0 :     if( h2o_ndx>0) call outfld( 'H2O_GAS',  vmr(1,1,h2o_ndx),  ncol ,lchnk )
         930             :
         931             :     ! reset O3S to O3 in the stratosphere ...
         932           0 :     if ( o3_ndx > 0 .and. o3s_ndx > 0 ) then
         933           0 :        o3s_loss = rate_diags_o3s_loss( reaction_rates, vmr, ncol )
         934           0 :        do i = 1,ncol
         935           0 :           vmr(i,1:troplev(i),o3s_ndx) = vmr(i,1:troplev(i),o3_ndx)
         936           0 :           vmr(i,troplev(i)+1:pver,o3s_ndx) = vmr(i,troplev(i)+1:pver,o3s_ndx) * exp(-delt*o3s_loss(i,troplev(i)+1:pver))
         937             :        end do
         938           0 :        call outfld( 'O3S_LOSS',  o3s_loss,  ncol ,lchnk )
         939             :     end if
         940             :
         941           0 :     if (convproc_do_aer) then
         942           0 :        call vmr2mmr( vmr(:ncol,:,:), mmr_new(:ncol,:,:), mbar(:ncol,:), ncol )
         943             :        ! mmr_new = average of mmr values before and after imp_sol
         944             :        mmr_new(:ncol,:,:) = 0.5_r8*( mmr(:ncol,:,:) + mmr_new(:ncol,:,:) )
         945           0 :        call het_diags( het_rates(:ncol,:,:), mmr_new(:ncol,:,:), pdel(:ncol,:), lchnk, ncol )
         946             :     endif
         947             :
         948             :     ! save h2so4 change by gas phase chem (for later new particle nucleation)
         949           0 :     if (ndx_h2so4 > 0) then
         950           0 :        del_h2so4_gasprod(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) - del_h2so4_gasprod(1:ncol,:)
         951             :     endif
         952             :
         953             : !
         954             : ! Aerosol processes ...
         955             : !
         956             :
         957             :     call aero_model_gasaerexch( imozart-1, ncol, lchnk, troplevchem, delt, reaction_rates, &
         958             :                                 tfld, pmid, pdel, mbar, relhum, &
         959             :                                 zm,  qh2o, cwat, cldfr, ncldwtr, &
         960             :                                 invariants(:,:,indexm), invariants, del_h2so4_gasprod,  &
         961           0 :                                 vmr0, vmr, pbuf )
         962             :
         963           0 :     if ( has_strato_chem ) then
         964             :
         965           0 :        wrk(:ncol,:) = (vmr(:ncol,:,h2o_ndx) - wrk(:ncol,:))*delt_inverse
         966             :        call outfld( 'QDCHEM',   wrk(:ncol,:),         ncol, lchnk )
         967             :
         968             :        !-----------------------------------------------------------------------
         969             :        !         ... aerosol settling
         970             :        !             first settle hno3(2) using radius ice
         971             :        !             secnd settle hno3(3) using radius large nat
         972             :        !-----------------------------------------------------------------------
         973             :        wrk(:,:) = vmr(:,:,h2o_ndx)
         974             : #ifdef ALT_SETTL
         975             :        where( h2o_cond(:,:) > 0._r8 )
         976             :           settl_rad(:,:) = radius_strat(:,:,3)
         977             :        elsewhere
         978             :           settl_rad(:,:) = 0._r8
         979             :        endwhere
         980             :        call strat_aer_settling( invariants(1,1,indexm), pmid, delt, zmid, tfld, &
         981             :             hno3_cond(1,1,2), settl_rad, ncol, lchnk, 1 )
         982             :
         983             :        where( h2o_cond(:,:) == 0._r8 )
         984             :           settl_rad(:,:) = radius_strat(:,:,2)
         985             :        elsewhere
         986             :           settl_rad(:,:) = 0._r8
         987             :        endwhere
         988             :        call strat_aer_settling( invariants(1,1,indexm), pmid, delt, zmid, tfld, &
         989             :             hno3_cond(1,1,2), settl_rad, ncol, lchnk, 2 )
         990             : #else
         991             :        call strat_aer_settling( invariants(1,1,indexm), pmid, delt, zmid, tfld, &
         992             :             hno3_cond(1,1,2), radius_strat(1,1,2), ncol, lchnk, 2 )
         993             : #endif
         994             :
         995             :        !-----------------------------------------------------------------------
         996             :        !        ... reform total hno3 and hcl = gas + all condensed
         997             :        !-----------------------------------------------------------------------
         998             : !      NOTE: vmr for hcl and hno3 is gas-phase at this point.
         999             : !            hno3_cond(:,k,1) = STS; hno3_cond(:,k,2) = NAT
        1000             :
        1001             :        do k = 1,pver
        1002             :           vmr(:,k,hno3_ndx) = vmr(:,k,hno3_ndx) + hno3_cond(:,k,1) &
        1003             :                + hno3_cond(:,k,2)
        1004             :           vmr(:,k,hcl_ndx)  = vmr(:,k,hcl_ndx)  + hcl_cond(:,k)
        1005             :
        1006             :        end do
        1007             :
        1008             :        wrk(:,:) = (vmr(:,:,h2o_ndx) - wrk(:,:))*delt_inverse
        1009             :        call outfld( 'QDSETT', wrk(:,:), ncol, lchnk )
        1010             :
        1011             :     endif
        1012             :
        1013             :     !-----------------------------------------------------------------------
        1014             :     !         ... Check for negative values and reset to zero
        1015             :     !-----------------------------------------------------------------------
        1016           0 :     call negtrc( 'After chemistry ', vmr, ncol )
        1017             :
        1018             :     !-----------------------------------------------------------------------
        1019             :     !         ... Set upper boundary mmr values
        1020             :     !-----------------------------------------------------------------------
        1021           0 :     call set_fstrat_vals( vmr, pmid, pint, troplev, calday, ncol,lchnk )
        1022             :
        1023             :     !-----------------------------------------------------------------------
        1024             :     !         ... Set fixed lower boundary mmr values
        1025             :     !-----------------------------------------------------------------------
        1026           0 :     call flbc_set( vmr, ncol, lchnk, map2chm )
        1027             :
        1028           0 :     if ( ghg_chem ) then
        1029           0 :        call ghg_chem_set_flbc( vmr, ncol )
        1030             :     endif
        1031             :
        1032             :     !-----------------------------------------------------------------------
        1033             :     ! force ion/electron balance -- ext forcings likely do not conserve charge
        1034             :     !-----------------------------------------------------------------------
        1035           0 :     call charge_balance( ncol, vmr )
        1036             :
        1037             :     !-----------------------------------------------------------------------
        1038             :     !         ... Xform from vmr to mmr
        1039             :     !-----------------------------------------------------------------------
        1040           0 :     call vmr2mmr( vmr(:ncol,:,:), mmr_tend(:ncol,:,:), mbar(:ncol,:), ncol )
        1041             :
        1042           0 :     call set_short_lived_species( mmr_tend, lchnk, ncol, pbuf )
        1043             :
        1044             :     !-----------------------------------------------------------------------
        1045             :     !         ... Form the tendencies
        1046             :     !-----------------------------------------------------------------------
        1047             :     do m = 1,gas_pcnst
        1048             :        mmr_new(:ncol,:,m) = mmr_tend(:ncol,:,m)
        1049           0 :        mmr_tend(:ncol,:,m) = (mmr_tend(:ncol,:,m) - mmr(:ncol,:,m))*delt_inverse
        1050             :     enddo
        1051             :
        1052           0 :     do m = 1,pcnst
        1053           0 :        n = map2chm(m)
        1054           0 :        if( n > 0 ) then
        1055           0 :           qtend(:ncol,:,m) = qtend(:ncol,:,m) + mmr_tend(:ncol,:,n)
        1056             :        end if
        1057             :     end do
        1058             :
        1059           0 :     tvs(:ncol) = tfld(:ncol,pver) * (1._r8 + qh2o(:ncol,pver))
        1060             :
        1061           0 :     sflx(:,:) = 0._r8
        1062           0 :     wind_speed(:ncol) = sqrt( ufld(:ncol,pver)*ufld(:ncol,pver) + vfld(:ncol,pver)*vfld(:ncol,pver) )
        1063           0 :     prect(:ncol) = precc(:ncol) + precl(:ncol)
        1064             :
        1065             :     call drydep( ocnfrac, icefrac, ts, ps,  &
        1066             :                  wind_speed, qh2o(:,pver), tfld(:,pver), pmid(:,pver), prect, &
        1067             :                  snowhland, fsds, depvel, sflx, mmr, &
        1068           0 :                  tvs, ncol, lchnk )
        1069             :
        1070           0 :     drydepflx(:,:) = 0._r8
        1071           0 :     wetdepflx_diag(:,:) = 0._r8
        1072           0 :     do m = 1,pcnst
        1073           0 :        n = map2chm( m )
        1074           0 :        if ( n > 0 ) then
        1075           0 :          if (cam_physpkg_is("cam7")) then
        1076             :            ! apply to qtend array
        1077           0 :            if (cnst_type(m).eq.'dry') then
        1078           0 :              qtend(:ncol,pver,m) = qtend(:ncol,pver,m) - sflx(:ncol,n)*rpdeldry(:ncol,pver)*gravit
        1079             :            else
        1080           0 :              qtend(:ncol,pver,m) = qtend(:ncol,pver,m) - sflx(:ncol,n)*rpdel(:ncol,pver)*gravit
        1081             :            end if
        1082             :          else
        1083             :            ! apply to emissions array
        1084           0 :            cflx(:ncol,m) = cflx(:ncol,m) - sflx(:ncol,n)
        1085             :          end if
        1086             :          drydepflx(:ncol,m) = sflx(:ncol,n)
        1087             :          wetdepflx_diag(:ncol,n) = wetdepflx(:ncol,m)
        1088             :        endif
        1089             :     end do
        1090             :
        1091           0 :     if (srf_ozone_pbf_ndx>0) then
        1092           0 :        call pbuf_get_field(pbuf, srf_ozone_pbf_ndx, srf_ozone_fld)
        1093           0 :        if (o3_ndx>0) then
        1094           0 :           srf_ozone_fld(:ncol) = vmr(:ncol,pver,o3_ndx)
        1095             :        else
        1096           0 :           srf_ozone_fld(:ncol) = invariants(:ncol,pver,o3_inv_ndx)/invariants(:ncol,pver,indexm)
        1097             :        endif
        1098             :     endif
        1099             :
        1100             :     call chm_diags( lchnk, ncol, vmr(:ncol,:,:), mmr_new(:ncol,:,:), &
        1101             :                     reaction_rates(:ncol,:,:), invariants(:ncol,:,:), depvel(:ncol,:),  sflx(:ncol,:), &
        1102           0 :                     mmr_tend(:ncol,:,:), pdel(:ncol,:), pmid(:ncol,:), troplev(:ncol), wetdepflx_diag(:ncol,:), &
        1103           0 :                     nhx_nitrogen_flx(:ncol), noy_nitrogen_flx(:ncol) )
        1104             :
        1105           0 :     call rate_diags_calc( reaction_rates(:,:,:), vmr(:,:,:), invariants(:,:,indexm), ncol, lchnk )
        1106             : !
        1107             : ! jfl
        1108             : !
        1109             : ! surface vmr
        1110             : !
        1111           0 :     if ( pm25_srf_diag ) then
        1112           0 :        pm25(:ncol) = mmr_new(:ncol,pver,cb1_ndx)   &
        1113             :             + mmr_new(:ncol,pver,cb2_ndx)   &
        1114             :             + mmr_new(:ncol,pver,oc1_ndx)   &
        1115             :             + mmr_new(:ncol,pver,oc2_ndx)   &
        1116             :             + mmr_new(:ncol,pver,dst1_ndx)  &
        1117             :             + mmr_new(:ncol,pver,dst2_ndx)  &
        1118             :             + mmr_new(:ncol,pver,sslt1_ndx) &
        1119             :             + mmr_new(:ncol,pver,sslt2_ndx) &
        1120             :             + mmr_new(:ncol,pver,soa_ndx)   &
        1121           0 :             + mmr_new(:ncol,pver,so4_ndx)
        1122             :        call outfld('PM25_SRF',pm25(:ncol) , ncol, lchnk )
        1123             :     endif
        1124           0 :     if ( pm25_srf_diag_soa ) then
        1125           0 :        pm25(:ncol) = mmr_new(:ncol,pver,cb1_ndx)   &
        1126             :             + mmr_new(:ncol,pver,cb2_ndx)   &
        1127             :             + mmr_new(:ncol,pver,oc1_ndx)   &
        1128             :             + mmr_new(:ncol,pver,oc2_ndx)   &
        1129             :             + mmr_new(:ncol,pver,dst1_ndx)  &
        1130             :             + mmr_new(:ncol,pver,dst2_ndx)  &
        1131             :             + mmr_new(:ncol,pver,sslt1_ndx) &
        1132             :             + mmr_new(:ncol,pver,sslt2_ndx) &
        1133             :             + mmr_new(:ncol,pver,soam_ndx)   &
        1134             :             + mmr_new(:ncol,pver,soai_ndx)   &
        1135             :             + mmr_new(:ncol,pver,soat_ndx)   &
        1136             :             + mmr_new(:ncol,pver,soab_ndx)   &
        1137             :             + mmr_new(:ncol,pver,soax_ndx)   &
        1138           0 :             + mmr_new(:ncol,pver,so4_ndx)
        1139             :        call outfld('PM25_SRF',pm25(:ncol) , ncol, lchnk )
        1140             :     endif
        1141             : !
        1142             : !
        1143           0 :     call outfld('Q_SRF',qh2o(:ncol,pver) , ncol, lchnk )
        1144           0 :     call outfld('U_SRF',ufld(:ncol,pver) , ncol, lchnk )
        1145           0 :     call outfld('V_SRF',vfld(:ncol,pver) , ncol, lchnk )
        1146             :
        1147             : !
        1148           0 :     if (.not.sad_pbf_ndx>0) then
        1149           0 :        deallocate(strato_sad)
        1150             :     endif
        1151             :
        1152           0 :   end subroutine gas_phase_chemdr
        1153             :
        1154             : end module mo_gas_phase_chemdr

Generated by: LCOV version 1.14