LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_gas_phase_chemdr.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 331 | 468 | 70.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 2 | 100.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     1490712 :   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        1536 :     call phys_getopts( history_scwaccm_forcing_out = history_scwaccm_forcing )
          75             :
          76        1536 :     call phys_getopts( convproc_do_aer_out = convproc_do_aer, history_cesm_forcing_out=history_cesm_forcing )
          77             :
          78        1536 :     ndx_h2so4 = get_spc_ndx('H2SO4')
          79             : !
          80             : ! CCMI
          81             : !
          82        1536 :     st80_25_ndx     = get_spc_ndx   ('ST80_25')
          83        1536 :     st80_25_tau_ndx = get_rxt_ndx   ('ST80_25_tau')
          84        1536 :     aoa_nh_ndx      = get_spc_ndx   ('AOA_NH')
          85        1536 :     aoa_nh_ext_ndx  = get_extfrc_ndx('AOA_NH')
          86        1536 :     nh_5_ndx        = get_spc_ndx('NH_5')
          87        1536 :     nh_50_ndx       = get_spc_ndx('NH_50')
          88        1536 :     nh_50w_ndx      = get_spc_ndx('NH_50W')
          89             : !
          90        1536 :     cb1_ndx         = get_spc_ndx('CB1')
          91        1536 :     cb2_ndx         = get_spc_ndx('CB2')
          92        1536 :     oc1_ndx         = get_spc_ndx('OC1')
          93        1536 :     oc2_ndx         = get_spc_ndx('OC2')
          94        1536 :     dst1_ndx        = get_spc_ndx('DST01')
          95        1536 :     dst2_ndx        = get_spc_ndx('DST02')
          96        1536 :     sslt1_ndx       = get_spc_ndx('SSLT01')
          97        1536 :     sslt2_ndx       = get_spc_ndx('SSLT02')
          98        1536 :     soa_ndx         = get_spc_ndx('SOA')
          99        1536 :     soam_ndx         = get_spc_ndx('SOAM')
         100        1536 :     soai_ndx         = get_spc_ndx('SOAI')
         101        1536 :     soat_ndx         = get_spc_ndx('SOAT')
         102        1536 :     soab_ndx         = get_spc_ndx('SOAB')
         103        1536 :     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        1536 :               .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        1536 :               .and. soam_ndx>0 .and. soai_ndx>0 .and. soat_ndx>0 .and. soab_ndx>0 .and. soax_ndx>0
         112             :
         113        1536 :     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        1536 :     call addfld('U_SRF',horiz_only,'I','m/s','bottom layer wind velocity' )
         117        1536 :     call addfld('V_SRF',horiz_only,'I','m/s','bottom layer wind velocity' )
         118        1536 :     call addfld('Q_SRF',horiz_only,'I','kg/kg','bottom layer specific humidity' )
         119             : !
         120        1536 :     het1_ndx= get_rxt_ndx('het1')
         121        1536 :     o3_ndx  = get_spc_ndx('O3')
         122        1536 :     o3_inv_ndx  = get_inv_ndx( 'O3' )
         123        1536 :     o3s_ndx = get_spc_ndx('O3S')
         124        1536 :     o_ndx   = get_spc_ndx('O')
         125        1536 :     o2_ndx  = get_spc_ndx('O2')
         126        1536 :     so4_ndx = get_spc_ndx('SO4')
         127        1536 :     h2o_ndx = get_spc_ndx('H2O')
         128        1536 :     hno3_ndx = get_spc_ndx('HNO3')
         129        1536 :     hcl_ndx  = get_spc_ndx('HCL')
         130        1536 :     call cnst_get_ind( 'CLDICE', cldice_ndx )
         131        1536 :     call cnst_get_ind( 'SNOWQM', snow_ndx, abort=.false. )
         132             :
         133        1536 :     if (o3_ndx>0 .or. o3_inv_ndx>0) then
         134        1536 :        srf_ozone_pbf_ndx = pbuf_get_index('SRFOZONE')
         135             :     endif
         136             :
         137       18432 :     do m = 1,extcnt
         138       16896 :        WRITE(UNIT=string, FMT='(I2.2)') m
         139       16896 :        extfrc_name(m) = 'extfrc_'// trim(string)
         140       35328 :        call addfld( extfrc_name(m), (/ 'lev' /), 'I', ' ', 'ext frcing' )
         141             :     end do
         142             :
         143       24576 :     do n = 1,rxt_tag_cnt
         144       23040 :        tag_names(n) = trim(rxt_tag_lst(n))
         145       23040 :        if (n<=phtcnt) then
         146        9216 :           call addfld( tag_names(n), (/ 'lev' /), 'I', '/s', 'photolysis rate constant' )
         147             :        else
         148       18432 :           ii = n-phtcnt
         149       18432 :           select case(num_rnts(ii))
         150             :           case(1)
         151        9216 :              unitstr='/s'
         152             :           case(2)
         153        7680 :              unitstr='cm3/molecules/s'
         154             :           case(3)
         155        1536 :              unitstr='cm6/molecules2/s'
         156             :           case default
         157       18432 :              call endrun('gas_phase_chemdr_inti: invalid value in num_rnts used to set units in reaction rate constant')
         158             :           end select
         159       36864 :           call addfld( tag_names(n), (/ 'lev' /), 'I', unitstr, 'reaction rate constant' )
         160             :        endif
         161       24576 :        if (history_scwaccm_forcing) then
         162           0 :           select case (trim(tag_names(n)))
         163             :           case ('jh2o_a', 'jh2o_b', 'jh2o_c' )
         164           0 :              call add_default( tag_names(n), 1, ' ')
         165             :           end select
         166             :        endif
         167             :     enddo
         168             :
         169        3072 :     call addfld( 'QDSAD',      (/ 'lev' /), 'I', '/s',      'water vapor sad delta' )
         170        3072 :     call addfld( 'SAD_STRAT',  (/ 'lev' /), 'I', 'cm2/cm3', 'stratospheric aerosol SAD' )
         171        3072 :     call addfld( 'SAD_SULFC',  (/ 'lev' /), 'I', 'cm2/cm3', 'chemical sulfate aerosol SAD' )
         172        3072 :     call addfld( 'SAD_SAGE',   (/ 'lev' /), 'I', 'cm2/cm3', 'SAGE sulfate aerosol SAD' )
         173        3072 :     call addfld( 'SAD_LNAT',   (/ 'lev' /), 'I', 'cm2/cm3', 'large-mode NAT aerosol SAD' )
         174        3072 :     call addfld( 'SAD_ICE',    (/ 'lev' /), 'I', 'cm2/cm3', 'water-ice aerosol SAD' )
         175        3072 :     call addfld( 'RAD_SULFC',  (/ 'lev' /), 'I', 'cm',      'chemical sad sulfate' )
         176        3072 :     call addfld( 'RAD_LNAT',   (/ 'lev' /), 'I', 'cm',      'large nat radius' )
         177        3072 :     call addfld( 'RAD_ICE',    (/ 'lev' /), 'I', 'cm',      'sad ice' )
         178        3072 :     call addfld( 'SAD_TROP',   (/ 'lev' /), 'I', 'cm2/cm3', 'tropospheric aerosol SAD' )
         179        3072 :     call addfld( 'SAD_AERO',   (/ 'lev' /), 'I', 'cm2/cm3', 'aerosol surface area density' )
         180        1536 :     if (history_cesm_forcing) then
         181           0 :        call add_default ('SAD_AERO',8,' ')
         182             :     endif
         183        3072 :     call addfld( 'REFF_AERO',  (/ 'lev' /), 'I', 'cm',      'aerosol effective radius' )
         184        3072 :     call addfld( 'SULF_TROP',  (/ 'lev' /), 'I', 'mol/mol', 'tropospheric aerosol SAD' )
         185        3072 :     call addfld( 'QDSETT',     (/ 'lev' /), 'I', '/s',      'water vapor settling delta' )
         186        3072 :     call addfld( 'QDCHEM',     (/ 'lev' /), 'I', '/s',      'water vapor chemistry delta')
         187        3072 :     call addfld( 'HNO3_TOTAL', (/ 'lev' /), 'I', 'mol/mol', 'total HNO3' )
         188        3072 :     call addfld( 'HNO3_STS',   (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HNO3' )
         189        3072 :     call addfld( 'HNO3_NAT',   (/ 'lev' /), 'I', 'mol/mol', 'NAT condensed HNO3' )
         190        3072 :     call addfld( 'HNO3_GAS',   (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hno3' )
         191        3072 :     call addfld( 'H2O_GAS',    (/ 'lev' /), 'I', 'mol/mol', 'gas-phase h2o' )
         192        3072 :     call addfld( 'HCL_TOTAL',  (/ 'lev' /), 'I', 'mol/mol', 'total hcl' )
         193        3072 :     call addfld( 'HCL_GAS',    (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hcl' )
         194        3072 :     call addfld( 'HCL_STS',    (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HCL' )
         195             :
         196        1536 :     if (het1_ndx>0) then
         197           0 :        call addfld( 'het1_total', (/ 'lev' /), 'I', '/s', 'total N2O5 + H2O het rate constant' )
         198             :     endif
         199        1536 :     call addfld( 'SZA', horiz_only, 'I', 'degrees', 'solar zenith angle' )
         200             :
         201        1536 :     call chm_diags_inti()
         202        1536 :     call rate_diags_init()
         203             :
         204             : !-----------------------------------------------------------------------
         205             : ! get pbuf indicies
         206             : !-----------------------------------------------------------------------
         207        1536 :     ndx_cldfr  = pbuf_get_index('CLD')
         208        1536 :     ndx_cmfdqr = pbuf_get_index('RPRDTOT')
         209        1536 :     ndx_nevapr = pbuf_get_index('NEVAPR')
         210        1536 :     ndx_prain  = pbuf_get_index('PRAIN')
         211        1536 :     ndx_cldtop = pbuf_get_index('CLDTOP')
         212             :
         213        1536 :     sad_pbf_ndx= pbuf_get_index('VOLC_SAD',errcode=err) ! prescribed  strat aerosols (volcanic)
         214        1536 :     if (.not.sad_pbf_ndx>0) sad_pbf_ndx = pbuf_get_index('SADSULF',errcode=err) ! CARMA's version of strat aerosols
         215             :
         216        1536 :     ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index
         217        1536 :     ion_temp_ndx = pbuf_get_index('TIon',errcode=err) ! ion temperature index
         218             :
         219             :     ! diagnostics for stratospheric heterogeneous reactions
         220        3072 :     call addfld( 'GAMMA_HET1', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         221        3072 :     call addfld( 'GAMMA_HET2', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         222        3072 :     call addfld( 'GAMMA_HET3', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         223        3072 :     call addfld( 'GAMMA_HET4', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         224        3072 :     call addfld( 'GAMMA_HET5', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         225        3072 :     call addfld( 'GAMMA_HET6', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
         226        3072 :     call addfld( 'WTPER',      (/ 'lev' /), 'I', '%', 'H2SO4 Weight Percent' )
         227             :
         228        3072 :     call addfld( 'O3S_LOSS',      (/ 'lev' /), 'I', '1/sec', 'O3S loss rate const' )
         229             :
         230        1536 :     call chem_prod_loss_diags_init
         231             :
         232             :     ! diagnostics for HEMCO ParaNOx extension
         233        1536 :     hco_jno2_idx = pbuf_get_index('HCO_IN_JNO2',errcode=err)
         234        1536 :     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        1536 :     rxt_jno2_idx  = get_rxt_ndx( 'jno2' )
         247        1536 :     rxt_joh_idx   = get_rxt_ndx( 'jo3_b' )
         248             :
         249        1536 :   end subroutine gas_phase_chemdr_inti
         250             :
         251             :
         252             : !-----------------------------------------------------------------------
         253             : !-----------------------------------------------------------------------
         254     1489176 :   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        1536 :     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     1489176 :     real(r8),       pointer    :: prain(:,:)
         374     1489176 :     real(r8),       pointer    :: nevapr(:,:)
         375     1489176 :     real(r8),       pointer    :: cmfdqr(:,:)
         376     1489176 :     real(r8),       pointer    :: cldfr(:,:)
         377     1489176 :     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     2978352 :     real(r8)     ::  invariants(ncol,pver,nfs)
         384     2978352 :     real(r8)     ::  col_dens(ncol,pver,max(1,nabscol))    ! column densities (molecules/cm^2)
         385     2978352 :     real(r8)     ::  col_delta(ncol,0:pver,max(1,nabscol)) ! layer column densities (molecules/cm^2)
         386     2978352 :     real(r8)     ::  extfrc(ncol,pver,max(1,extcnt))
         387     2978352 :     real(r8)     ::  vmr(ncol,pver,gas_pcnst)         ! xported species (vmr)
         388     2978352 :     real(r8)     ::  reaction_rates(ncol,pver,max(1,rxntot))      ! reaction rates
         389     2978352 :     real(r8)     ::  depvel(ncol,gas_pcnst)                ! dry deposition velocity (cm/s)
         390     2978352 :     real(r8)     ::  het_rates(ncol,pver,max(1,gas_pcnst)) ! washout rate (1/s)
         391             :     real(r8), dimension(ncol,pver) :: &
         392     2978352 :          h2ovmr, &                                         ! water vapor volume mixing ratio
         393     2978352 :          mbar, &                                           ! mean wet atmospheric mass ( amu )
         394     2978352 :          zmid, &                                           ! midpoint geopotential in km
         395     2978352 :          sulfate, &                                        ! trop sulfate aerosols
         396     2978352 :          pmb                                               ! pressure at midpoints ( hPa )
         397             :     real(r8), dimension(ncol,pver) :: &
         398     2978352 :          cwat, &                                           ! cloud water mass mixing ratio (kg/kg)
         399     2978352 :          wrk
         400             :     real(r8), dimension(ncol,pver+1) :: &
         401     2978352 :          zintr                                              ! interface geopotential in km realitive to surf
         402             :     real(r8), dimension(ncol,pver+1) :: &
         403     2978352 :          zint                                              ! interface geopotential in km
         404             :     real(r8), dimension(ncol)  :: &
         405     2978352 :          zen_angle, &                                      ! solar zenith angles
         406     2978352 :          zsurf, &                                          ! surface height (m)
         407     2978352 :          rlats, rlons                                      ! chunk latitudes and longitudes (radians)
         408     2978352 :     real(r8) :: sza(ncol)                                  ! solar zenith angles (degrees)
         409             :     real(r8), parameter :: rad2deg = 180._r8/pi                ! radians to degrees conversion factor
         410     2978352 :     real(r8) :: relhum(ncol,pver)                          ! relative humidity
         411     2978352 :     real(r8) :: satv(ncol,pver)                            ! wrk array for relative humidity
         412     2978352 :     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     1489176 :     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     2978352 :     real(r8) :: o2mmr(ncol,pver)               ! o2 concentration (kg/kg)
         428     2978352 :     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     2978352 :     real(r8) :: hno3_gas(ncol,pver)            ! hno3 gas phase concentration (mol/mol)
         432     2978352 :     real(r8) :: hno3_cond(ncol,pver,2)         ! hno3 condensed phase concentration (mol/mol)
         433     2978352 :     real(r8) :: hcl_gas(ncol,pver)             ! hcl gas phase concentration (mol/mol)
         434     2978352 :     real(r8) :: hcl_cond(ncol,pver)            ! hcl condensed phase concentration (mol/mol)
         435     2978352 :     real(r8) :: h2o_gas(ncol,pver)             ! h2o gas phase concentration (mol/mol)
         436     2978352 :     real(r8) :: h2o_cond(ncol,pver)            ! h2o condensed phase concentration (mol/mol)
         437             :     real(r8) :: cldice(pcols,pver)             ! cloud water "ice" (kg/kg)
         438     2978352 :     real(r8) :: radius_strat(ncol,pver,3)      ! radius of sulfate, nat, & ice ( cm )
         439     2978352 :     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     2978352 :     real(r8) :: del_h2so4_gasprod(ncol,pver)
         446     2978352 :     real(r8) :: vmr0(ncol,pver,gas_pcnst)
         447             :
         448             :   ! for HEMCO-CESM ... passing J-values to ParaNOx ship plume extension
         449     1489176 :     real(r8), pointer :: hco_j_tmp_fld(:)    ! J-value pointer (sfc only) [1/s]
         450             :
         451             : !
         452             : ! CCMI
         453             : !
         454             :     real(r8) :: xlat
         455     2978352 :     real(r8) :: pm25(ncol)
         456             :
         457     2978352 :     real(r8) :: dlats(ncol)
         458             :
         459             :     real(r8), dimension(ncol,pver) :: &      ! aerosol reaction diagnostics
         460     2978352 :         gprob_n2o5, &
         461     2978352 :         gprob_cnt_hcl, &
         462     2978352 :         gprob_cnt_h2o, &
         463     2978352 :         gprob_bnt_h2o, &
         464     2978352 :         gprob_hocl_hcl, &
         465     2978352 :         gprob_hobr_hcl, &
         466     2978352 :         wtper
         467             :
         468     1489176 :     real(r8), pointer :: ele_temp_fld(:,:) ! electron temperature pointer
         469     1489176 :     real(r8), pointer :: ion_temp_fld(:,:) ! ion temperature pointer
         470     2978352 :     real(r8) :: prod_out(ncol,pver,max(1,clscnt4))
         471     2978352 :     real(r8) :: loss_out(ncol,pver,max(1,clscnt4))
         472             :
         473     2978352 :     real(r8) :: o3s_loss(ncol,pver)
         474     1489176 :     real(r8), pointer :: srf_ozone_fld(:)
         475             :
         476     1489176 :     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     1489176 :        ele_temp_fld => tfld
         481     1489176 :        ion_temp_fld => tfld
         482             :     endif
         483             :
         484             :     ! initialize to NaN to hopefully catch user defined rxts that go unset
         485     1489176 :     reaction_rates(:,:,:) = nan
         486             :
         487     1489176 :     delt_inverse = 1._r8 / delt
         488             :     !-----------------------------------------------------------------------
         489             :     !        ... Get chunck latitudes and longitudes
         490             :     !-----------------------------------------------------------------------
         491     1489176 :     call get_rlat_all_p( lchnk, ncol, rlats )
         492     1489176 :     call get_rlon_all_p( lchnk, ncol, rlons )
         493     1489176 :     tim_ndx = pbuf_old_tim_idx()
         494     4467528 :     call pbuf_get_field(pbuf, ndx_prain,      prain,  start=(/1,1/), kount=(/ncol,pver/))
         495    10424232 :     call pbuf_get_field(pbuf, ndx_cldfr,        cldfr, start=(/1,1,tim_ndx/), kount=(/ncol,pver,1/) )
         496     4467528 :     call pbuf_get_field(pbuf, ndx_cmfdqr,     cmfdqr, start=(/1,1/), kount=(/ncol,pver/))
         497     4467528 :     call pbuf_get_field(pbuf, ndx_nevapr,     nevapr, start=(/1,1/), kount=(/ncol,pver/))
         498     1489176 :     call pbuf_get_field(pbuf, ndx_cldtop,     cldtop )
         499             :
         500     1489176 :     reff_strat(:,:) = 0._r8
         501             :
         502    24865776 :     dlats(:) = rlats(:)*rad2deg ! convert to degrees
         503             :
         504             :     !-----------------------------------------------------------------------
         505             :     !        ... Calculate cosine of zenith angle
         506             :     !            then cast back to angle (radians)
         507             :     !-----------------------------------------------------------------------
         508     1489176 :     call zenith( calday, rlats, rlons, zen_angle, ncol )
         509    24865776 :     zen_angle(:) = acos( zen_angle(:) )
         510             :
         511    24865776 :     sza(:) = zen_angle(:) * rad2deg
         512     1489176 :     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    24865776 :     zsurf(:ncol) = rga * phis(:ncol)
         519   139982544 :     do k = 1,pver
         520  2312517168 :        zintr(:ncol,k) = m2km * zi(:ncol,k)
         521  2312517168 :        zmid(:ncol,k) = m2km * (zm(:ncol,k) + zsurf(:ncol))
         522  2312517168 :        zint(:ncol,k) = m2km * (zi(:ncol,k) + zsurf(:ncol))
         523  2314006344 :        pmb(:ncol,k)  = Pa2mb * pmid(:ncol,k)
         524             :     end do
         525    24865776 :     zint(:ncol,pver+1) = m2km * (zi(:ncol,pver+1) + zsurf(:ncol))
         526    24865776 :     zintr(:ncol,pver+1)= m2km *  zi(:ncol,pver+1)
         527             :
         528             :     !-----------------------------------------------------------------------
         529             :     !        ... map incoming concentrations to working array
         530             :     !-----------------------------------------------------------------------
         531    62545392 :     do m = 1,pcnst
         532    61056216 :        n = map2chm(m)
         533    62545392 :        if( n > 0 ) then
         534 71734196664 :           mmr(:ncol,:,n) = q(:ncol,:,m)
         535             :        end if
         536             :     end do
         537             :
         538     1489176 :     call get_short_lived_species( mmr, lchnk, ncol, pbuf )
         539             :
         540             :     !-----------------------------------------------------------------------
         541             :     !        ... Set atmosphere mean mass
         542             :     !-----------------------------------------------------------------------
         543     1489176 :     call set_mean_mass( ncol, mmr, mbar )
         544             :
         545             :     !-----------------------------------------------------------------------
         546             :     !        ... Xform from mmr to vmr
         547             :     !-----------------------------------------------------------------------
         548 16820023584 :     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     1489176 :     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     1489176 :     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     1489176 :     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     1489176 :     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     1489176 :     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     1489176 :     if (h2o_ndx>0) then
         597             :        !-----------------------------------------------------------------------
         598             :        !        ... store water vapor in wrk variable
         599             :        !-----------------------------------------------------------------------
         600  2314006344 :        qh2o(:ncol,:) = mmr(:ncol,:,h2o_ndx)
         601  2314006344 :        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     1489176 :     call charge_balance( ncol, vmr )
         617             :
         618             :     !-----------------------------------------------------------------------
         619             :     !        ... Set the "invariants"
         620             :     !-----------------------------------------------------------------------
         621     1489176 :     call setinv( invariants, tfld, h2ovmr, vmr, pmid, ncol, lchnk, pbuf )
         622             :
         623             :     !-----------------------------------------------------------------------
         624             :     !        ... stratosphere aerosol surface area
         625             :     !-----------------------------------------------------------------------
         626     1489176 :     if (sad_pbf_ndx>0) then
         627     1489176 :        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     1489176 :     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           0 :        do k = 1,pver
         644           0 :           hno3_gas(:,k)   = vmr(:,k,hno3_ndx)
         645           0 :           h2o_gas(:,k)    = h2ovmr(:,k)
         646           0 :           hcl_gas(:,k)    = vmr(:,k,hcl_ndx)
         647           0 :           wrk(:,k)        = h2ovmr(:,k)
         648           0 :           if (snow_ndx>0) then
         649           0 :              cldice(:ncol,k) = q(:ncol,k,cldice_ndx) + q(:ncol,k,snow_ndx)
         650             :           else
         651           0 :              cldice(:ncol,k) = q(:ncol,k,cldice_ndx)
         652             :           endif
         653             :        end do
         654           0 :        do m = 1,2
         655           0 :           do k = 1,pver
         656           0 :              hno3_cond(:,k,m) = 0._r8
         657             :           end do
         658             :        end do
         659             :
         660           0 :        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           0 :             hno3_cond, h2o_gas, h2o_cond, hcl_gas, hcl_cond, strato_sad(:ncol,:), radius_strat, &
         667           0 :             sad_strat, ncol, pbuf )
         668             :
         669             : !      NOTE: output of total HNO3 is before vmr is set to gas-phase.
         670           0 :        call outfld( 'HNO3_TOTAL', vmr(:ncol,:,hno3_ndx), ncol ,lchnk )
         671             :
         672             :
         673           0 :        do k = 1,pver
         674           0 :           vmr(:,k,hno3_ndx) = hno3_gas(:,k)
         675           0 :           h2ovmr(:,k)       = h2o_gas(:,k)
         676           0 :           vmr(:,k,h2o_ndx)  = h2o_gas(:,k)
         677           0 :           wrk(:,k)          = (h2ovmr(:,k) - wrk(:,k))*delt_inverse
         678             :        end do
         679             :
         680           0 :        call outfld( 'QDSAD', wrk(:,:), ncol, lchnk )
         681             : !
         682           0 :        call outfld( 'SAD_STRAT',  strato_sad(:ncol,:), ncol, lchnk )
         683           0 :        call outfld( 'SAD_SULFC',  sad_strat(:,:,1), ncol, lchnk )
         684           0 :        call outfld( 'SAD_LNAT',   sad_strat(:,:,2), ncol, lchnk )
         685           0 :        call outfld( 'SAD_ICE',    sad_strat(:,:,3), ncol, lchnk )
         686             : !
         687           0 :        call outfld( 'RAD_SULFC',  radius_strat(:,:,1), ncol, lchnk )
         688           0 :        call outfld( 'RAD_LNAT',   radius_strat(:,:,2), ncol, lchnk )
         689           0 :        call outfld( 'RAD_ICE',    radius_strat(:,:,3), ncol, lchnk )
         690             : !
         691           0 :        call outfld( 'HNO3_GAS',   vmr(:ncol,:,hno3_ndx), ncol, lchnk )
         692           0 :        call outfld( 'HNO3_STS',   hno3_cond(:,:,1), ncol, lchnk )
         693           0 :        call outfld( 'HNO3_NAT',   hno3_cond(:,:,2), ncol, lchnk )
         694             : !
         695           0 :        call outfld( 'HCL_TOTAL',  vmr(:ncol,:,hcl_ndx), ncol, lchnk )
         696           0 :        call outfld( 'HCL_GAS',    hcl_gas (:,:), ncol ,lchnk )
         697           0 :        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           0 :             gprob_hocl_hcl, gprob_hobr_hcl, wtper )
         707             :
         708           0 :        call outfld( 'GAMMA_HET1', gprob_n2o5    (:ncol,:), ncol, lchnk )
         709           0 :        call outfld( 'GAMMA_HET2', gprob_cnt_h2o (:ncol,:), ncol, lchnk )
         710           0 :        call outfld( 'GAMMA_HET3', gprob_bnt_h2o (:ncol,:), ncol, lchnk )
         711           0 :        call outfld( 'GAMMA_HET4', gprob_cnt_hcl (:ncol,:), ncol, lchnk )
         712           0 :        call outfld( 'GAMMA_HET5', gprob_hocl_hcl(:ncol,:), ncol, lchnk )
         713           0 :        call outfld( 'GAMMA_HET6', gprob_hobr_hcl(:ncol,:), ncol, lchnk )
         714           0 :        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     1489176 :     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     1489176 :     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     1489176 :     call setrxt( reaction_rates, tfld, invariants(1,1,indexm), ncol )
         733             :
         734  2314006344 :     sulfate(:,:) = 0._r8
         735     1489176 :     if ( .not. carma_hetchem_feedback ) then
         736     1489176 :        if( so4_ndx < 1 ) then ! get offline so4 field if not prognostic
         737     1489176 :           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   139982544 :     do k = 1, pver
         747  2314006344 :        do i = 1, ncol
         748  2312517168 :           if (k < troplevchem(i)) then
         749  1198140589 :              sulfate(i,k) = 0.0_r8
         750             :           end if
         751             :        end do
         752             :     end do
         753             :
         754     1489176 :     call outfld( 'SULF_TROP', sulfate(:ncol,:), ncol, lchnk )
         755             :
         756             :     !-----------------------------------------------------------------
         757             :     !   ... compute the relative humidity
         758             :     !-----------------------------------------------------------------
         759   139982544 :     do k = 1, pver
         760   139982544 :        call qsat(tfld(1:ncol,k), pmid(1:ncol,k), satv(1:ncol,k), satq(1:ncol,k), ncol)
         761             :     end do
         762             :
         763   139982544 :     do k = 1,pver
         764  2312517168 :        relhum(:,k) = .622_r8 * h2ovmr(:,k) / satq(:,k)
         765  2314006344 :        relhum(:,k) = max( 0._r8,min( 1._r8,relhum(:,k) ) )
         766             :     end do
         767             :
         768  2314006344 :     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     1489176 :                  troplevchem, dlats, ncol, sad_trop, reff, cwat, mbar, pbuf )
         773             :
         774   543653136 :     call outfld( 'SAD_TROP', sad_trop(:ncol,:), ncol, lchnk )
         775             :
         776             :     ! Add trop/strat components of SAD for output
         777  2314006344 :     sad_trop(:ncol,:)=sad_trop(:ncol,:)+strato_sad(:ncol,:)
         778   543653136 :     call outfld( 'SAD_AERO', sad_trop(:ncol,:), ncol, lchnk )
         779             :
         780             :     ! Add trop/strat components of effective radius for output
         781  2314006344 :     reff(:ncol,:)=reff(:ncol,:)+reff_strat(:ncol,:)
         782   543653136 :     call outfld( 'REFF_AERO', reff(:ncol,:), ncol, lchnk )
         783             :
         784     1489176 :     if (het1_ndx>0) then
         785           0 :        call outfld( 'het1_total', reaction_rates(:,:,het1_ndx), ncol, lchnk )
         786             :     endif
         787             :
         788     1489176 :     if (ghg_chem) then
         789     1489176 :        call ghg_chem_set_rates( reaction_rates, latmapback, zen_angle, ncol, lchnk )
         790             :     endif
         791             :
         792    19359288 :     do i = phtcnt+1,rxt_tag_cnt
         793    19359288 :        call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk )
         794             :     enddo
         795             :
         796     1489176 :     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     1489176 :     call setcol( col_delta, col_dens )
         805             :
         806             :     !-----------------------------------------------------------------------
         807             :     !           ... Calculate the photodissociation rates
         808             :     !-----------------------------------------------------------------------
         809             :
         810     1489176 :     esfact = 1._r8
         811             :     call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr  , &
         812     1489176 :          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     1489176 :                       esfact, vmr, invariants, ncol, lchnk, pbuf )
         820             :
         821     5956704 :     do i = 1,phtcnt
         822     5956704 :        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     1489176 :     call phtadj( reaction_rates, invariants, invariants(:,:,indexm), ncol,pver )
         829             :
         830     1489176 :     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     1489176 :     if ( o2_ndx > 0 .and. o_ndx > 0 ) then
         854           0 :        do k = 1,pver
         855           0 :           o2mmr(:ncol,k) = mmr(:ncol,k,o2_ndx)
         856           0 :           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     1489176 :                  pmid, mbar, rlats, calday, ncol, rlons, pbuf )
         862             :     ! include forcings from fire emissions ...
         863     1489176 :     call fire_emissions_vrt( ncol, lchnk, zint, fire_sflx, fire_ztop, extfrc )
         864             :
         865    17870112 :     do m = 1,extcnt
         866    16380936 :        if( m /= aoa_nh_ext_ndx ) then
         867  1539807984 :           do k = 1,pver
         868 25454069784 :              extfrc(:ncol,k,m) = extfrc(:ncol,k,m) / invariants(:ncol,k,indexm)
         869             :           end do
         870             :        endif
         871    17870112 :        call outfld( extfrc_name(m), extfrc(:ncol,:,m), ncol, lchnk )
         872             :     end do
         873             :
         874             :     !-----------------------------------------------------------------------
         875             :     !        ... Form the washout rates
         876             :     !-----------------------------------------------------------------------
         877     1489176 :     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 71735685840 :        het_rates = 0._r8
         886             :     end if
         887             : !
         888             : ! CCMI
         889             : !
         890             : ! set loss to below the tropopause only
         891             : !
         892     1489176 :     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    24865776 :     ltrop_sol(:ncol) = 0 ! apply solver to all levels
         899             :
         900             :     ! save h2so4 before gas phase chem (for later new particle nucleation)
         901     1489176 :     if (ndx_h2so4 > 0) then
         902  2314006344 :        del_h2so4_gasprod(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4)
         903             :     else
         904           0 :        del_h2so4_gasprod(:,:) = 0.0_r8
         905             :     endif
         906             :
         907 71735685840 :     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     1489176 :     if ( has_strato_chem ) wrk(:,:) = vmr(:,:,h2o_ndx)
         921     1489176 :     call t_startf('imp_sol')
         922             :     !
         923             :     call imp_sol( vmr, reaction_rates, het_rates, extfrc, delt, &
         924     1489176 :                   ncol,pver, lchnk,  prod_out, loss_out )
         925             :
         926     1489176 :     call t_stopf('imp_sol')
         927             :
         928     1489176 :     call chem_prod_loss_diags_out( ncol, lchnk, vmr, reaction_rates, prod_out, loss_out, invariants(:ncol,:,indexm) )
         929     1489176 :     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     1489176 :     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     1489176 :     if (convproc_do_aer) then
         942 33638557992 :        call vmr2mmr( vmr(:ncol,:,:), mmr_new(:ncol,:,:), mbar(:ncol,:), ncol )
         943             :        ! mmr_new = average of mmr values before and after imp_sol
         944 71735685840 :        mmr_new(:ncol,:,:) = 0.5_r8*( mmr(:ncol,:,:) + mmr_new(:ncol,:,:) )
         945 17362187544 :        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     1489176 :     if (ndx_h2so4 > 0) then
         950  2314006344 :        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     1489176 :                                 vmr0, vmr, pbuf )
         962             :
         963     1489176 :     if ( has_strato_chem ) then
         964             :
         965           0 :        wrk(:ncol,:) = (vmr(:ncol,:,h2o_ndx) - wrk(:ncol,:))*delt_inverse
         966           0 :        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           0 :        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           0 :             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           0 :        do k = 1,pver
        1002           0 :           vmr(:,k,hno3_ndx) = vmr(:,k,hno3_ndx) + hno3_cond(:,k,1) &
        1003           0 :                + hno3_cond(:,k,2)
        1004           0 :           vmr(:,k,hcl_ndx)  = vmr(:,k,hcl_ndx)  + hcl_cond(:,k)
        1005             :
        1006             :        end do
        1007             :
        1008           0 :        wrk(:,:) = (vmr(:,:,h2o_ndx) - wrk(:,:))*delt_inverse
        1009           0 :        call outfld( 'QDSETT', wrk(:,:), ncol, lchnk )
        1010             :
        1011             :     endif
        1012             :
        1013             :     !-----------------------------------------------------------------------
        1014             :     !         ... Check for negative values and reset to zero
        1015             :     !-----------------------------------------------------------------------
        1016     1489176 :     call negtrc( 'After chemistry ', vmr, ncol )
        1017             :
        1018             :     !-----------------------------------------------------------------------
        1019             :     !         ... Set upper boundary mmr values
        1020             :     !-----------------------------------------------------------------------
        1021     1489176 :     call set_fstrat_vals( vmr, pmid, pint, troplev, calday, ncol,lchnk )
        1022             :
        1023             :     !-----------------------------------------------------------------------
        1024             :     !         ... Set fixed lower boundary mmr values
        1025             :     !-----------------------------------------------------------------------
        1026     1489176 :     call flbc_set( vmr, ncol, lchnk, map2chm )
        1027             :
        1028     1489176 :     if ( ghg_chem ) then
        1029     1489176 :        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     1489176 :     call charge_balance( ncol, vmr )
        1036             :
        1037             :     !-----------------------------------------------------------------------
        1038             :     !         ... Xform from vmr to mmr
        1039             :     !-----------------------------------------------------------------------
        1040 33638557992 :     call vmr2mmr( vmr(:ncol,:,:), mmr_tend(:ncol,:,:), mbar(:ncol,:), ncol )
        1041             :
        1042     1489176 :     call set_short_lived_species( mmr_tend, lchnk, ncol, pbuf )
        1043             :
        1044             :     !-----------------------------------------------------------------------
        1045             :     !         ... Form the tendencies
        1046             :     !-----------------------------------------------------------------------
        1047    47653632 :     do m = 1,gas_pcnst
        1048 71734196664 :        mmr_new(:ncol,:,m) = mmr_tend(:ncol,:,m)
        1049 71735685840 :        mmr_tend(:ncol,:,m) = (mmr_tend(:ncol,:,m) - mmr(:ncol,:,m))*delt_inverse
        1050             :     enddo
        1051             :
        1052    62545392 :     do m = 1,pcnst
        1053    61056216 :        n = map2chm(m)
        1054    62545392 :        if( n > 0 ) then
        1055 71734196664 :           qtend(:ncol,:,m) = qtend(:ncol,:,m) + mmr_tend(:ncol,:,n)
        1056             :        end if
        1057             :     end do
        1058             :
        1059    24865776 :     tvs(:ncol) = tfld(:ncol,pver) * (1._r8 + qh2o(:ncol,pver))
        1060             :
        1061     1489176 :     sflx(:,:) = 0._r8
        1062    24865776 :     wind_speed(:ncol) = sqrt( ufld(:ncol,pver)*ufld(:ncol,pver) + vfld(:ncol,pver)*vfld(:ncol,pver) )
        1063    24865776 :     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     1489176 :                  tvs, ncol, lchnk )
        1069             :
        1070     1489176 :     drydepflx(:,:) = 0._r8
        1071     1489176 :     wetdepflx_diag(:,:) = 0._r8
        1072    62545392 :     do m = 1,pcnst
        1073    61056216 :        n = map2chm( m )
        1074    62545392 :        if ( n > 0 ) then
        1075    46164456 :          if (cam_physpkg_is("cam7")) then
        1076             :            ! apply to qtend array
        1077    46164456 :            if (cnst_type(m).eq.'dry') then
        1078   745973280 :              qtend(:ncol,pver,m) = qtend(:ncol,pver,m) - sflx(:ncol,n)*rpdeldry(:ncol,pver)*gravit
        1079             :            else
        1080    24865776 :              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   770839056 :          drydepflx(:ncol,m) = sflx(:ncol,n)
        1087   770839056 :          wetdepflx_diag(:ncol,n) = wetdepflx(:ncol,m)
        1088             :        endif
        1089             :     end do
        1090             :
        1091     1489176 :     if (srf_ozone_pbf_ndx>0) then
        1092     1489176 :        call pbuf_get_field(pbuf, srf_ozone_pbf_ndx, srf_ozone_fld)
        1093     1489176 :        if (o3_ndx>0) then
        1094           0 :           srf_ozone_fld(:ncol) = vmr(:ncol,pver,o3_ndx)
        1095             :        else
        1096    24865776 :           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             :                     mmr_tend(:ncol,:,:), pdel(:ncol,:), pmid(:ncol,:), troplev(:ncol), wetdepflx_diag(:ncol,:), &
        1103 35084328552 :                     nhx_nitrogen_flx(:ncol), noy_nitrogen_flx(:ncol) )
        1104             :
        1105     1489176 :     call rate_diags_calc( reaction_rates(:,:,:), vmr(:,:,:), invariants(:,:,indexm), ncol, lchnk )
        1106             : !
        1107             : ! jfl
        1108             : !
        1109             : ! surface vmr
        1110             : !
        1111     1489176 :     if ( pm25_srf_diag ) then
        1112           0 :        pm25(:ncol) = mmr_new(:ncol,pver,cb1_ndx)   &
        1113           0 :             + mmr_new(:ncol,pver,cb2_ndx)   &
        1114           0 :             + mmr_new(:ncol,pver,oc1_ndx)   &
        1115           0 :             + mmr_new(:ncol,pver,oc2_ndx)   &
        1116           0 :             + mmr_new(:ncol,pver,dst1_ndx)  &
        1117           0 :             + mmr_new(:ncol,pver,dst2_ndx)  &
        1118           0 :             + mmr_new(:ncol,pver,sslt1_ndx) &
        1119           0 :             + mmr_new(:ncol,pver,sslt2_ndx) &
        1120           0 :             + mmr_new(:ncol,pver,soa_ndx)   &
        1121           0 :             + mmr_new(:ncol,pver,so4_ndx)
        1122           0 :        call outfld('PM25_SRF',pm25(:ncol) , ncol, lchnk )
        1123             :     endif
        1124     1489176 :     if ( pm25_srf_diag_soa ) then
        1125           0 :        pm25(:ncol) = mmr_new(:ncol,pver,cb1_ndx)   &
        1126           0 :             + mmr_new(:ncol,pver,cb2_ndx)   &
        1127           0 :             + mmr_new(:ncol,pver,oc1_ndx)   &
        1128           0 :             + mmr_new(:ncol,pver,oc2_ndx)   &
        1129           0 :             + mmr_new(:ncol,pver,dst1_ndx)  &
        1130           0 :             + mmr_new(:ncol,pver,dst2_ndx)  &
        1131           0 :             + mmr_new(:ncol,pver,sslt1_ndx) &
        1132           0 :             + mmr_new(:ncol,pver,sslt2_ndx) &
        1133           0 :             + mmr_new(:ncol,pver,soam_ndx)   &
        1134           0 :             + mmr_new(:ncol,pver,soai_ndx)   &
        1135           0 :             + mmr_new(:ncol,pver,soat_ndx)   &
        1136           0 :             + mmr_new(:ncol,pver,soab_ndx)   &
        1137           0 :             + mmr_new(:ncol,pver,soax_ndx)   &
        1138           0 :             + mmr_new(:ncol,pver,so4_ndx)
        1139           0 :        call outfld('PM25_SRF',pm25(:ncol) , ncol, lchnk )
        1140             :     endif
        1141             : !
        1142             : !
        1143     1489176 :     call outfld('Q_SRF',qh2o(:ncol,pver) , ncol, lchnk )
        1144     1489176 :     call outfld('U_SRF',ufld(:ncol,pver) , ncol, lchnk )
        1145     1489176 :     call outfld('V_SRF',vfld(:ncol,pver) , ncol, lchnk )
        1146             :
        1147             : !
        1148     1489176 :     if (.not.sad_pbf_ndx>0) then
        1149           0 :        deallocate(strato_sad)
        1150             :     endif
        1151             :
        1152     2978352 :   end subroutine gas_phase_chemdr
        1153             :
        1154             : end module mo_gas_phase_chemdr

Generated by: LCOV version 1.14