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