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 22656 : 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 768 : call phys_getopts( history_scwaccm_forcing_out = history_scwaccm_forcing )
76 :
77 768 : call phys_getopts( convproc_do_aer_out = convproc_do_aer, history_cesm_forcing_out=history_cesm_forcing )
78 :
79 768 : ndx_h2so4 = get_spc_ndx('H2SO4')
80 : !
81 : ! CCMI
82 : !
83 768 : st80_25_ndx = get_spc_ndx ('ST80_25')
84 768 : st80_25_tau_ndx = get_rxt_ndx ('ST80_25_tau')
85 768 : aoa_nh_ndx = get_spc_ndx ('AOA_NH')
86 768 : aoa_nh_ext_ndx = get_extfrc_ndx('AOA_NH')
87 768 : nh_5_ndx = get_spc_ndx('NH_5')
88 768 : nh_50_ndx = get_spc_ndx('NH_50')
89 768 : nh_50w_ndx = get_spc_ndx('NH_50W')
90 : !
91 768 : cb1_ndx = get_spc_ndx('CB1')
92 768 : cb2_ndx = get_spc_ndx('CB2')
93 768 : oc1_ndx = get_spc_ndx('OC1')
94 768 : oc2_ndx = get_spc_ndx('OC2')
95 768 : dst1_ndx = get_spc_ndx('DST01')
96 768 : dst2_ndx = get_spc_ndx('DST02')
97 768 : sslt1_ndx = get_spc_ndx('SSLT01')
98 768 : sslt2_ndx = get_spc_ndx('SSLT02')
99 768 : soa_ndx = get_spc_ndx('SOA')
100 768 : soam_ndx = get_spc_ndx('SOAM')
101 768 : soai_ndx = get_spc_ndx('SOAI')
102 768 : soat_ndx = get_spc_ndx('SOAT')
103 768 : soab_ndx = get_spc_ndx('SOAB')
104 768 : 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 768 : .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 768 : .and. soam_ndx>0 .and. soai_ndx>0 .and. soat_ndx>0 .and. soab_ndx>0 .and. soax_ndx>0
113 :
114 768 : 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 768 : call addfld('U_SRF',horiz_only,'I','m/s','bottom layer wind velocity' )
118 768 : call addfld('V_SRF',horiz_only,'I','m/s','bottom layer wind velocity' )
119 768 : call addfld('Q_SRF',horiz_only,'I','kg/kg','bottom layer specific humidity' )
120 : !
121 768 : het1_ndx= get_rxt_ndx('het1')
122 768 : o3_ndx = get_spc_ndx('O3')
123 768 : o3_inv_ndx = get_inv_ndx( 'O3' )
124 768 : o3s_ndx = get_spc_ndx('O3S')
125 768 : o_ndx = get_spc_ndx('O')
126 768 : o2_ndx = get_spc_ndx('O2')
127 768 : so4_ndx = get_spc_ndx('SO4')
128 768 : h2o_ndx = get_spc_ndx('H2O')
129 768 : hno3_ndx = get_spc_ndx('HNO3')
130 768 : hcl_ndx = get_spc_ndx('HCL')
131 768 : call cnst_get_ind( 'CLDICE', cldice_ndx )
132 768 : call cnst_get_ind( 'SNOWQM', snow_ndx, abort=.false. )
133 :
134 768 : if (o3_ndx>0 .or. o3_inv_ndx>0) then
135 768 : srf_ozone_pbf_ndx = pbuf_get_index('SRFOZONE')
136 : endif
137 :
138 19968 : do m = 1,extcnt
139 19200 : WRITE(UNIT=string, FMT='(I2.2)') m
140 19200 : extfrc_name(m) = 'extfrc_'// trim(string)
141 39168 : call addfld( extfrc_name(m), (/ 'lev' /), 'I', ' ', 'ext frcing' )
142 : end do
143 :
144 466176 : do n = 1,rxt_tag_cnt
145 465408 : tag_names(n) = trim(rxt_tag_lst(n))
146 465408 : if (n<=phtcnt) then
147 159744 : call addfld( tag_names(n), (/ 'lev' /), 'I', '/s', 'photolysis rate constant' )
148 : else
149 385536 : ii = n-phtcnt
150 385536 : select case(num_rnts(ii))
151 : case(1)
152 13824 : unitstr='/s'
153 : case(2)
154 329472 : unitstr='cm3/molecules/s'
155 : case(3)
156 42240 : unitstr='cm6/molecules2/s'
157 : case default
158 385536 : call endrun('gas_phase_chemdr_inti: invalid value in num_rnts used to set units in reaction rate constant')
159 : end select
160 771072 : call addfld( tag_names(n), (/ 'lev' /), 'I', unitstr, 'reaction rate constant' )
161 : endif
162 466176 : 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 1536 : call addfld( 'QDSAD', (/ 'lev' /), 'I', '/s', 'water vapor sad delta' )
171 1536 : call addfld( 'SAD_STRAT', (/ 'lev' /), 'I', 'cm2/cm3', 'stratospheric aerosol SAD' )
172 1536 : call addfld( 'SAD_SULFC', (/ 'lev' /), 'I', 'cm2/cm3', 'chemical sulfate aerosol SAD' )
173 1536 : call addfld( 'SAD_SAGE', (/ 'lev' /), 'I', 'cm2/cm3', 'SAGE sulfate aerosol SAD' )
174 1536 : call addfld( 'SAD_LNAT', (/ 'lev' /), 'I', 'cm2/cm3', 'large-mode NAT aerosol SAD' )
175 1536 : call addfld( 'SAD_ICE', (/ 'lev' /), 'I', 'cm2/cm3', 'water-ice aerosol SAD' )
176 1536 : call addfld( 'RAD_SULFC', (/ 'lev' /), 'I', 'cm', 'chemical sad sulfate' )
177 1536 : call addfld( 'RAD_LNAT', (/ 'lev' /), 'I', 'cm', 'large nat radius' )
178 1536 : call addfld( 'RAD_ICE', (/ 'lev' /), 'I', 'cm', 'sad ice' )
179 1536 : call addfld( 'SAD_TROP', (/ 'lev' /), 'I', 'cm2/cm3', 'tropospheric aerosol SAD' )
180 1536 : call addfld( 'SAD_AERO', (/ 'lev' /), 'I', 'cm2/cm3', 'aerosol surface area density' )
181 768 : if (history_cesm_forcing) then
182 0 : call add_default ('SAD_AERO',8,' ')
183 : endif
184 1536 : call addfld( 'REFF_AERO', (/ 'lev' /), 'I', 'cm', 'aerosol effective radius' )
185 1536 : call addfld( 'REFF_TROP', (/ 'lev' /), 'I', 'cm', 'tropospheric aerosol effective radius' )
186 1536 : call addfld( 'REFF_STRAT', (/ 'lev' /), 'I', 'cm', 'stratospheric aerosol effective radius' )
187 1536 : call addfld( 'SULF_TROP', (/ 'lev' /), 'I', 'mol/mol', 'tropospheric aerosol SAD' )
188 1536 : call addfld( 'QDSETT', (/ 'lev' /), 'I', '/s', 'water vapor settling delta' )
189 1536 : call addfld( 'QDCHEM', (/ 'lev' /), 'I', '/s', 'water vapor chemistry delta')
190 1536 : call addfld( 'HNO3_TOTAL', (/ 'lev' /), 'I', 'mol/mol', 'total HNO3' )
191 1536 : call addfld( 'HNO3_STS', (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HNO3' )
192 1536 : call addfld( 'HNO3_NAT', (/ 'lev' /), 'I', 'mol/mol', 'NAT condensed HNO3' )
193 1536 : call addfld( 'HNO3_GAS', (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hno3' )
194 1536 : call addfld( 'H2O_GAS', (/ 'lev' /), 'I', 'mol/mol', 'gas-phase h2o' )
195 1536 : call addfld( 'HCL_TOTAL', (/ 'lev' /), 'I', 'mol/mol', 'total hcl' )
196 1536 : call addfld( 'HCL_GAS', (/ 'lev' /), 'I', 'mol/mol', 'gas-phase hcl' )
197 1536 : call addfld( 'HCL_STS', (/ 'lev' /), 'I', 'mol/mol', 'STS condensed HCL' )
198 :
199 768 : if (het1_ndx>0) then
200 1536 : call addfld( 'het1_total', (/ 'lev' /), 'I', '/s', 'total N2O5 + H2O het rate constant' )
201 : endif
202 768 : call addfld( 'SZA', horiz_only, 'I', 'degrees', 'solar zenith angle' )
203 :
204 768 : call chm_diags_inti()
205 768 : call rate_diags_init()
206 :
207 : !-----------------------------------------------------------------------
208 : ! get pbuf indicies
209 : !-----------------------------------------------------------------------
210 768 : ndx_cldfr = pbuf_get_index('CLD')
211 768 : ndx_cmfdqr = pbuf_get_index('RPRDTOT')
212 768 : ndx_nevapr = pbuf_get_index('NEVAPR')
213 768 : ndx_prain = pbuf_get_index('PRAIN')
214 768 : ndx_cldtop = pbuf_get_index('CLDTOP')
215 :
216 768 : sad_pbf_ndx= pbuf_get_index('VOLC_SAD',errcode=err) ! prescribed strat aerosols (volcanic)
217 :
218 768 : ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index
219 768 : ion_temp_ndx = pbuf_get_index('TIon',errcode=err) ! ion temperature index
220 :
221 768 : jno2_pbuf_ndx = pbuf_get_index('JNO2',errcode=err)
222 768 : 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 1536 : call addfld( 'GAMMA_HET1', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
228 1536 : call addfld( 'GAMMA_HET2', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
229 1536 : call addfld( 'GAMMA_HET3', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
230 1536 : call addfld( 'GAMMA_HET4', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
231 1536 : call addfld( 'GAMMA_HET5', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
232 1536 : call addfld( 'GAMMA_HET6', (/ 'lev' /), 'I', '1', 'Reaction Probability' )
233 1536 : call addfld( 'WTPER', (/ 'lev' /), 'I', '%', 'H2SO4 Weight Percent' )
234 :
235 1536 : call addfld( 'O3S_LOSS', (/ 'lev' /), 'I', '1/sec', 'O3S loss rate const' )
236 :
237 768 : call chem_prod_loss_diags_init
238 :
239 : ! diagnostics for HEMCO ParaNOx extension
240 768 : hco_jno2_idx = pbuf_get_index('HCO_IN_JNO2',errcode=err)
241 768 : 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 768 : rxt_jno2_idx = get_rxt_ndx( 'jno2' )
254 768 : rxt_joh_idx = get_rxt_ndx( 'jo3_b' )
255 :
256 768 : end subroutine gas_phase_chemdr_inti
257 :
258 :
259 : !-----------------------------------------------------------------------
260 : !-----------------------------------------------------------------------
261 21888 : 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 768 : 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 21888 : real(r8), pointer :: prain(:,:)
383 21888 : real(r8), pointer :: nevapr(:,:)
384 21888 : real(r8), pointer :: cmfdqr(:,:)
385 21888 : real(r8), pointer :: cldfr(:,:)
386 21888 : 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 43776 : real(r8) :: invariants(ncol,pver,nfs)
393 43776 : real(r8) :: col_dens(ncol,pver,max(1,nabscol)) ! column densities (molecules/cm^2)
394 43776 : real(r8) :: col_delta(ncol,0:pver,max(1,nabscol)) ! layer column densities (molecules/cm^2)
395 43776 : real(r8) :: extfrc(ncol,pver,max(1,extcnt))
396 43776 : real(r8) :: vmr(ncol,pver,gas_pcnst) ! xported species (vmr)
397 43776 : real(r8) :: reaction_rates(ncol,pver,max(1,rxntot)) ! reaction rates
398 43776 : real(r8) :: depvel(ncol,gas_pcnst) ! dry deposition velocity (cm/s)
399 43776 : real(r8) :: het_rates(ncol,pver,max(1,gas_pcnst)) ! washout rate (1/s)
400 : real(r8), dimension(ncol,pver) :: &
401 43776 : h2ovmr, & ! water vapor volume mixing ratio
402 43776 : mbar, & ! mean wet atmospheric mass ( amu )
403 43776 : zmid, & ! midpoint geopotential in km
404 43776 : sulfate, & ! trop sulfate aerosols
405 43776 : pmb ! pressure at midpoints ( hPa )
406 : real(r8), dimension(ncol,pver) :: &
407 43776 : cwat, & ! cloud water mass mixing ratio (kg/kg)
408 43776 : wrk
409 : real(r8), dimension(ncol,pver+1) :: &
410 43776 : zintr ! interface geopotential in km realitive to surf
411 : real(r8), dimension(ncol,pver+1) :: &
412 43776 : zint ! interface geopotential in km
413 : real(r8), dimension(ncol) :: &
414 43776 : zen_angle, & ! solar zenith angles
415 43776 : zsurf, & ! surface height (m)
416 43776 : rlats, rlons ! chunk latitudes and longitudes (radians)
417 43776 : real(r8) :: sza(ncol) ! solar zenith angles (degrees)
418 : real(r8), parameter :: rad2deg = 180._r8/pi ! radians to degrees conversion factor
419 43776 : real(r8) :: relhum(ncol,pver) ! relative humidity
420 43776 : real(r8) :: satv(ncol,pver) ! wrk array for relative humidity
421 43776 : 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 21888 : 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 43776 : real(r8) :: o2mmr(ncol,pver) ! o2 concentration (kg/kg)
437 43776 : 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 43776 : real(r8) :: hno3_gas(ncol,pver) ! hno3 gas phase concentration (mol/mol)
441 43776 : real(r8) :: hno3_cond(ncol,pver,2) ! hno3 condensed phase concentration (mol/mol)
442 43776 : real(r8) :: hcl_gas(ncol,pver) ! hcl gas phase concentration (mol/mol)
443 43776 : real(r8) :: hcl_cond(ncol,pver) ! hcl condensed phase concentration (mol/mol)
444 43776 : real(r8) :: h2o_gas(ncol,pver) ! h2o gas phase concentration (mol/mol)
445 43776 : real(r8) :: h2o_cond(ncol,pver) ! h2o condensed phase concentration (mol/mol)
446 : real(r8) :: cldice(pcols,pver) ! cloud water "ice" (kg/kg)
447 43776 : real(r8) :: radius_strat(ncol,pver,3) ! radius of sulfate, nat, & ice ( cm )
448 43776 : 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 43776 : real(r8) :: del_h2so4_gasprod(ncol,pver)
455 43776 : real(r8) :: vmr0(ncol,pver,gas_pcnst)
456 :
457 : ! for HEMCO-CESM ... passing J-values to ParaNOx ship plume extension
458 21888 : real(r8), pointer :: hco_j_tmp_fld(:) ! J-value pointer (sfc only) [1/s]
459 :
460 : !
461 : ! CCMI
462 : !
463 : real(r8) :: xlat
464 43776 : real(r8) :: pm25(ncol)
465 :
466 43776 : real(r8) :: dlats(ncol)
467 :
468 : real(r8), dimension(ncol,pver) :: & ! aerosol reaction diagnostics
469 43776 : gprob_n2o5, &
470 43776 : gprob_cnt_hcl, &
471 43776 : gprob_cnt_h2o, &
472 43776 : gprob_bnt_h2o, &
473 43776 : gprob_hocl_hcl, &
474 43776 : gprob_hobr_hcl, &
475 43776 : wtper
476 :
477 21888 : real(r8), pointer :: ele_temp_fld(:,:) ! electron temperature pointer
478 21888 : real(r8), pointer :: ion_temp_fld(:,:) ! ion temperature pointer
479 43776 : real(r8) :: prod_out(ncol,pver,max(1,clscnt4))
480 43776 : real(r8) :: loss_out(ncol,pver,max(1,clscnt4))
481 :
482 43776 : real(r8) :: o3s_loss(ncol,pver)
483 21888 : real(r8), pointer :: srf_ozone_fld(:)
484 21888 : real(r8), pointer :: jno2_fld_ptr(:,:)
485 :
486 21888 : if ( ele_temp_ndx>0 .and. ion_temp_ndx>0 ) then
487 21888 : call pbuf_get_field(pbuf, ele_temp_ndx, ele_temp_fld)
488 21888 : call pbuf_get_field(pbuf, ion_temp_ndx, ion_temp_fld)
489 : else
490 0 : ele_temp_fld => tfld
491 0 : ion_temp_fld => tfld
492 : endif
493 :
494 : ! initialize to NaN to hopefully catch user defined rxts that go unset
495 21888 : reaction_rates(:,:,:) = nan
496 :
497 21888 : delt_inverse = 1._r8 / delt
498 : !-----------------------------------------------------------------------
499 : ! ... Get chunck latitudes and longitudes
500 : !-----------------------------------------------------------------------
501 21888 : call get_rlat_all_p( lchnk, ncol, rlats )
502 21888 : call get_rlon_all_p( lchnk, ncol, rlons )
503 21888 : tim_ndx = pbuf_old_tim_idx()
504 65664 : call pbuf_get_field(pbuf, ndx_prain, prain, start=(/1,1/), kount=(/ncol,pver/))
505 153216 : call pbuf_get_field(pbuf, ndx_cldfr, cldfr, start=(/1,1,tim_ndx/), kount=(/ncol,pver,1/) )
506 65664 : call pbuf_get_field(pbuf, ndx_cmfdqr, cmfdqr, start=(/1,1/), kount=(/ncol,pver/))
507 65664 : call pbuf_get_field(pbuf, ndx_nevapr, nevapr, start=(/1,1/), kount=(/ncol,pver/))
508 21888 : call pbuf_get_field(pbuf, ndx_cldtop, cldtop )
509 :
510 21888 : 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 21888 : reff_strat(:,:) = 0._r8
515 :
516 284544 : dlats(:) = rlats(:)*rad2deg ! convert to degrees
517 :
518 : !-----------------------------------------------------------------------
519 : ! ... Calculate cosine of zenith angle
520 : ! then cast back to angle (radians)
521 : !-----------------------------------------------------------------------
522 21888 : call zenith( calday, rlats, rlons, zen_angle, ncol )
523 284544 : zen_angle(:) = acos( zen_angle(:) )
524 :
525 284544 : sza(:) = zen_angle(:) * rad2deg
526 21888 : 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 284544 : zsurf(:ncol) = rga * phis(:ncol)
533 2867328 : do k = 1,pver
534 36990720 : zintr(:ncol,k) = m2km * zi(:ncol,k)
535 36990720 : zmid(:ncol,k) = m2km * (zm(:ncol,k) + zsurf(:ncol))
536 36990720 : zint(:ncol,k) = m2km * (zi(:ncol,k) + zsurf(:ncol))
537 37012608 : pmb(:ncol,k) = Pa2mb * pmid(:ncol,k)
538 : end do
539 284544 : zint(:ncol,pver+1) = m2km * (zi(:ncol,pver+1) + zsurf(:ncol))
540 284544 : zintr(:ncol,pver+1)= m2km * zi(:ncol,pver+1)
541 :
542 : !-----------------------------------------------------------------------
543 : ! ... map incoming concentrations to working array
544 : !-----------------------------------------------------------------------
545 2232576 : do m = 1,pcnst
546 2210688 : n = map2chm(m)
547 2232576 : if( n > 0 ) then
548 3331134720 : mmr(:ncol,:,n) = q(:ncol,:,m)
549 : end if
550 : end do
551 :
552 21888 : call get_short_lived_species( mmr, lchnk, ncol, pbuf )
553 :
554 : !-----------------------------------------------------------------------
555 : ! ... Set atmosphere mean mass
556 : !-----------------------------------------------------------------------
557 21888 : call set_mean_mass( ncol, mmr, mbar )
558 :
559 : !-----------------------------------------------------------------------
560 : ! ... Xform from mmr to vmr
561 : !-----------------------------------------------------------------------
562 5144774400 : 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 21888 : 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 21888 : 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 21888 : 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 21888 : 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 21888 : 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 21888 : if (h2o_ndx>0) then
611 : !-----------------------------------------------------------------------
612 : ! ... store water vapor in wrk variable
613 : !-----------------------------------------------------------------------
614 37012608 : qh2o(:ncol,:) = mmr(:ncol,:,h2o_ndx)
615 37012608 : 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 21888 : call charge_balance( ncol, vmr )
631 :
632 : !-----------------------------------------------------------------------
633 : ! ... Set the "invariants"
634 : !-----------------------------------------------------------------------
635 21888 : call setinv( invariants, tfld, h2ovmr, vmr, pmid, ncol, lchnk, pbuf )
636 :
637 : !-----------------------------------------------------------------------
638 : ! ... stratosphere aerosol surface area
639 : !-----------------------------------------------------------------------
640 21888 : if (sad_pbf_ndx>0) then
641 0 : call pbuf_get_field(pbuf, sad_pbf_ndx, strato_sad)
642 : else
643 21888 : allocate(strato_sad(pcols,pver))
644 48394368 : strato_sad(:,:) = 0._r8
645 :
646 : ! Prognostic modal stratospheric sulfate: compute dry strato_sad
647 21888 : call aero_model_strat_surfarea( state, ncol, mmr, pmid, tfld, troplevchem, pbuf, strato_sad, reff_strat )
648 :
649 : endif
650 :
651 21888 : stratochem: if ( has_strato_chem ) then
652 : !-----------------------------------------------------------------------
653 : ! ... initialize condensed and gas phases; all hno3 to gas
654 : !-----------------------------------------------------------------------
655 37012608 : hcl_cond(:,:) = 0.0_r8
656 37012608 : hcl_gas (:,:) = 0.0_r8
657 2867328 : do k = 1,pver
658 36990720 : hno3_gas(:,k) = vmr(:,k,hno3_ndx)
659 36990720 : h2o_gas(:,k) = h2ovmr(:,k)
660 36990720 : hcl_gas(:,k) = vmr(:,k,hcl_ndx)
661 36990720 : wrk(:,k) = h2ovmr(:,k)
662 2867328 : if (snow_ndx>0) then
663 36990720 : 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 65664 : do m = 1,2
669 5756544 : do k = 1,pver
670 74025216 : hno3_cond(:,k,m) = 0._r8
671 : end do
672 : end do
673 :
674 37012608 : 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 21888 : sad_strat, ncol, pbuf )
682 :
683 : ! NOTE: output of total HNO3 is before vmr is set to gas-phase.
684 21888 : call outfld( 'HNO3_TOTAL', vmr(:ncol,:,hno3_ndx), ncol ,lchnk )
685 :
686 :
687 2867328 : do k = 1,pver
688 36990720 : vmr(:,k,hno3_ndx) = hno3_gas(:,k)
689 36990720 : h2ovmr(:,k) = h2o_gas(:,k)
690 36990720 : vmr(:,k,h2o_ndx) = h2o_gas(:,k)
691 37012608 : wrk(:,k) = (h2ovmr(:,k) - wrk(:,k))*delt_inverse
692 : end do
693 :
694 21888 : call outfld( 'QDSAD', wrk(:,:), ncol, lchnk )
695 : !
696 21888 : call outfld( 'SAD_STRAT', strato_sad(:ncol,:), ncol, lchnk )
697 21888 : call outfld( 'SAD_SULFC', sad_strat(:,:,1), ncol, lchnk )
698 21888 : call outfld( 'SAD_LNAT', sad_strat(:,:,2), ncol, lchnk )
699 21888 : call outfld( 'SAD_ICE', sad_strat(:,:,3), ncol, lchnk )
700 : !
701 21888 : call outfld( 'RAD_SULFC', radius_strat(:,:,1), ncol, lchnk )
702 21888 : call outfld( 'RAD_LNAT', radius_strat(:,:,2), ncol, lchnk )
703 21888 : call outfld( 'RAD_ICE', radius_strat(:,:,3), ncol, lchnk )
704 : !
705 21888 : call outfld( 'HNO3_GAS', vmr(:ncol,:,hno3_ndx), ncol, lchnk )
706 21888 : call outfld( 'HNO3_STS', hno3_cond(:,:,1), ncol, lchnk )
707 21888 : call outfld( 'HNO3_NAT', hno3_cond(:,:,2), ncol, lchnk )
708 : !
709 21888 : call outfld( 'HCL_TOTAL', vmr(:ncol,:,hcl_ndx), ncol, lchnk )
710 21888 : call outfld( 'HCL_GAS', hcl_gas (:,:), ncol ,lchnk )
711 21888 : 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 21888 : gprob_hocl_hcl, gprob_hobr_hcl, wtper )
721 :
722 21888 : call outfld( 'GAMMA_HET1', gprob_n2o5 (:ncol,:), ncol, lchnk )
723 21888 : call outfld( 'GAMMA_HET2', gprob_cnt_h2o (:ncol,:), ncol, lchnk )
724 21888 : call outfld( 'GAMMA_HET3', gprob_bnt_h2o (:ncol,:), ncol, lchnk )
725 21888 : call outfld( 'GAMMA_HET4', gprob_cnt_hcl (:ncol,:), ncol, lchnk )
726 21888 : call outfld( 'GAMMA_HET5', gprob_hocl_hcl(:ncol,:), ncol, lchnk )
727 21888 : call outfld( 'GAMMA_HET6', gprob_hobr_hcl(:ncol,:), ncol, lchnk )
728 21888 : 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 21888 : if (hcl_ndx>0) then
735 37012608 : vmr(:,:,hcl_ndx) = hcl_gas(:,:)
736 : endif
737 :
738 : !-----------------------------------------------------------------------
739 : ! ... Set the column densities at the upper boundary
740 : !-----------------------------------------------------------------------
741 21888 : 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 21888 : call setrxt( reaction_rates, tfld, invariants(1,1,indexm), ncol )
747 :
748 37012608 : sulfate(:,:) = 0._r8
749 21888 : if ( .not. carma_hetchem_feedback ) then
750 21888 : if( so4_ndx < 1 ) then ! get offline so4 field if not prognostic
751 21888 : 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 2867328 : do k = 1, pver
761 37012608 : do i = 1, ncol
762 36990720 : if (k < troplevchem(i)) then
763 29345395 : sulfate(i,k) = 0.0_r8
764 : end if
765 : end do
766 : end do
767 :
768 21888 : call outfld( 'SULF_TROP', sulfate(:ncol,:), ncol, lchnk )
769 :
770 : !-----------------------------------------------------------------
771 : ! ... compute the relative humidity
772 : !-----------------------------------------------------------------
773 2867328 : do k = 1, pver
774 2867328 : call qsat(tfld(1:ncol,k), pmid(1:ncol,k), satv(1:ncol,k), satq(1:ncol,k), ncol)
775 : end do
776 :
777 2867328 : do k = 1,pver
778 36990720 : relhum(:,k) = .622_r8 * h2ovmr(:,k) / satq(:,k)
779 37012608 : relhum(:,k) = max( 0._r8,min( 1._r8,relhum(:,k) ) )
780 : end do
781 :
782 37012608 : 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 21888 : troplevchem, dlats, ncol, sad_trop, reff, cwat, mbar, pbuf )
787 :
788 37012608 : call outfld( 'SAD_TROP', sad_trop(:ncol,:), ncol, lchnk )
789 :
790 : ! Add trop/strat components of SAD for output
791 37012608 : sad_trop(:ncol,:)=sad_trop(:ncol,:)+strato_sad(:ncol,:)
792 37012608 : call outfld( 'SAD_AERO', sad_trop(:ncol,:), ncol, lchnk )
793 :
794 : ! Add trop/strat components of effective radius for output
795 37012608 : call outfld( 'REFF_TROP', reff(:ncol,:), ncol, lchnk )
796 37012608 : call outfld( 'REFF_STRAT', reff_strat(:ncol,:), ncol, lchnk )
797 37012608 : reff(:ncol,:)=reff(:ncol,:)+reff_strat(:ncol,:)
798 37012608 : call outfld( 'REFF_AERO', reff(:ncol,:), ncol, lchnk )
799 :
800 21888 : if (het1_ndx>0) then
801 21888 : call outfld( 'het1_total', reaction_rates(:,:,het1_ndx), ncol, lchnk )
802 : endif
803 :
804 21888 : if (ghg_chem) then
805 0 : call ghg_chem_set_rates( reaction_rates, latmapback, zen_angle, ncol, lchnk )
806 : endif
807 :
808 11009664 : do i = phtcnt+1,rxt_tag_cnt
809 11009664 : call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk )
810 : enddo
811 :
812 21888 : 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 21888 : call setcol( col_delta, col_dens )
821 :
822 : !-----------------------------------------------------------------------
823 : ! ... Calculate the photodissociation rates
824 : !-----------------------------------------------------------------------
825 :
826 21888 : esfact = 1._r8
827 : call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr , &
828 21888 : 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 21888 : esfact, vmr, invariants, ncol, lchnk, pbuf )
836 :
837 2298240 : do i = 1,phtcnt
838 2298240 : call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk )
839 : enddo
840 :
841 21888 : 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 21888 : call phtadj( reaction_rates, invariants, invariants(:,:,indexm), ncol,pver )
849 :
850 21888 : 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 21888 : if ( o2_ndx > 0 .and. o_ndx > 0 ) then
874 2867328 : do k = 1,pver
875 36990720 : o2mmr(:ncol,k) = mmr(:ncol,k,o2_ndx)
876 37012608 : 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 21888 : pmid, mbar, rlats, calday, ncol, rlons, pbuf )
882 : ! include forcings from fire emissions ...
883 21888 : call fire_emissions_vrt( ncol, lchnk, zint, fire_sflx, fire_ztop, extfrc )
884 :
885 569088 : do m = 1,extcnt
886 547200 : if( m /= aoa_nh_ext_ndx ) then
887 71683200 : do k = 1,pver
888 925315200 : extfrc(:ncol,k,m) = extfrc(:ncol,k,m) / invariants(:ncol,k,indexm)
889 : end do
890 : endif
891 569088 : call outfld( extfrc_name(m), extfrc(:ncol,:,m), ncol, lchnk )
892 : end do
893 :
894 : !-----------------------------------------------------------------------
895 : ! ... Form the washout rates
896 : !-----------------------------------------------------------------------
897 21888 : 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 5144774400 : het_rates = 0._r8
906 : end if
907 : !
908 : ! CCMI
909 : !
910 : ! set loss to below the tropopause only
911 : !
912 21888 : 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 284544 : ltrop_sol(:ncol) = 0 ! apply solver to all levels
919 :
920 : ! save h2so4 before gas phase chem (for later new particle nucleation)
921 21888 : if (ndx_h2so4 > 0) then
922 37012608 : del_h2so4_gasprod(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4)
923 : else
924 0 : del_h2so4_gasprod(:,:) = 0.0_r8
925 : endif
926 :
927 5144774400 : 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 37012608 : if ( has_strato_chem ) wrk(:,:) = vmr(:,:,h2o_ndx)
941 21888 : call t_startf('imp_sol')
942 : !
943 : call imp_sol( vmr, reaction_rates, het_rates, extfrc, delt, &
944 21888 : ncol,pver, lchnk, prod_out, loss_out )
945 :
946 21888 : call t_stopf('imp_sol')
947 :
948 21888 : call chem_prod_loss_diags_out( ncol, lchnk, vmr, reaction_rates, prod_out, loss_out, invariants(:ncol,:,indexm) )
949 21888 : 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 21888 : 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 21888 : if (convproc_do_aer) then
962 10289526912 : call vmr2mmr( vmr(:ncol,:,:), mmr_new(:ncol,:,:), mbar(:ncol,:), ncol )
963 : ! mmr_new = average of mmr values before and after imp_sol
964 5144774400 : mmr_new(:ncol,:,:) = 0.5_r8*( mmr(:ncol,:,:) + mmr_new(:ncol,:,:) )
965 5181765120 : 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 21888 : if (ndx_h2so4 > 0) then
970 37012608 : 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 21888 : vmr0, vmr, pbuf )
982 :
983 21888 : if ( has_strato_chem ) then
984 :
985 37012608 : wrk(:ncol,:) = (vmr(:ncol,:,h2o_ndx) - wrk(:ncol,:))*delt_inverse
986 21888 : 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 37012608 : 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 21888 : 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 2867328 : do k = 1,pver
1022 5690880 : vmr(:,k,hno3_ndx) = vmr(:,k,hno3_ndx) + hno3_cond(:,k,1) &
1023 42681600 : + hno3_cond(:,k,2)
1024 37012608 : vmr(:,k,hcl_ndx) = vmr(:,k,hcl_ndx) + hcl_cond(:,k)
1025 :
1026 : end do
1027 :
1028 37012608 : wrk(:,:) = (vmr(:,:,h2o_ndx) - wrk(:,:))*delt_inverse
1029 21888 : call outfld( 'QDSETT', wrk(:,:), ncol, lchnk )
1030 :
1031 : endif
1032 :
1033 : !-----------------------------------------------------------------------
1034 : ! ... Check for negative values and reset to zero
1035 : !-----------------------------------------------------------------------
1036 21888 : call negtrc( 'After chemistry ', vmr, ncol )
1037 :
1038 : !-----------------------------------------------------------------------
1039 : ! ... Set upper boundary mmr values
1040 : !-----------------------------------------------------------------------
1041 21888 : call set_fstrat_vals( vmr, pmid, pint, troplev, calday, ncol,lchnk )
1042 :
1043 : !-----------------------------------------------------------------------
1044 : ! ... Set fixed lower boundary mmr values
1045 : !-----------------------------------------------------------------------
1046 21888 : call flbc_set( vmr, ncol, lchnk, map2chm )
1047 :
1048 21888 : if ( ghg_chem ) then
1049 0 : 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 21888 : call charge_balance( ncol, vmr )
1056 :
1057 : !-----------------------------------------------------------------------
1058 : ! ... Xform from vmr to mmr
1059 : !-----------------------------------------------------------------------
1060 10289526912 : call vmr2mmr( vmr(:ncol,:,:), mmr_tend(:ncol,:,:), mbar(:ncol,:), ncol )
1061 :
1062 21888 : call set_short_lived_species( mmr_tend, lchnk, ncol, pbuf )
1063 :
1064 : !-----------------------------------------------------------------------
1065 : ! ... Form the tendencies
1066 : !-----------------------------------------------------------------------
1067 3064320 : do m = 1,gas_pcnst
1068 5144752512 : mmr_new(:ncol,:,m) = mmr_tend(:ncol,:,m)
1069 5144774400 : mmr_tend(:ncol,:,m) = (mmr_tend(:ncol,:,m) - mmr(:ncol,:,m))*delt_inverse
1070 : enddo
1071 :
1072 2232576 : do m = 1,pcnst
1073 2210688 : n = map2chm(m)
1074 2232576 : if( n > 0 ) then
1075 3331134720 : qtend(:ncol,:,m) = qtend(:ncol,:,m) + mmr_tend(:ncol,:,n)
1076 : end if
1077 : end do
1078 :
1079 284544 : tvs(:ncol) = tfld(:ncol,pver) * (1._r8 + qh2o(:ncol,pver))
1080 :
1081 21888 : sflx(:,:) = 0._r8
1082 284544 : wind_speed(:ncol) = sqrt( ufld(:ncol,pver)*ufld(:ncol,pver) + vfld(:ncol,pver)*vfld(:ncol,pver) )
1083 284544 : 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 21888 : tvs, ncol, lchnk )
1089 :
1090 21888 : drydepflx(:,:) = 0._r8
1091 21888 : wetdepflx_diag(:,:) = 0._r8
1092 2232576 : do m = 1,pcnst
1093 2210688 : n = map2chm( m )
1094 2232576 : if ( n > 0 ) then
1095 1969920 : if (cam_physpkg_is("cam7")) then
1096 : ! apply to qtend array
1097 0 : if (cnst_type(m).eq.'dry') then
1098 0 : qtend(:ncol,pver,m) = qtend(:ncol,pver,m) - sflx(:ncol,n)*rpdeldry(:ncol,pver)*gravit
1099 : else
1100 0 : 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 25608960 : cflx(:ncol,m) = cflx(:ncol,m) - sflx(:ncol,n)
1105 : end if
1106 25608960 : drydepflx(:ncol,m) = sflx(:ncol,n)
1107 25608960 : wetdepflx_diag(:ncol,n) = wetdepflx(:ncol,m)
1108 : endif
1109 : end do
1110 :
1111 21888 : if (srf_ozone_pbf_ndx>0) then
1112 21888 : call pbuf_get_field(pbuf, srf_ozone_pbf_ndx, srf_ozone_fld)
1113 21888 : if (o3_ndx>0) then
1114 284544 : srf_ozone_fld(:ncol) = vmr(:ncol,pver,o3_ndx)
1115 : else
1116 0 : 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 10442611584 : nhx_nitrogen_flx(:ncol), noy_nitrogen_flx(:ncol) )
1124 :
1125 21888 : call rate_diags_calc( reaction_rates(:,:,:), vmr(:,:,:), invariants(:,:,indexm), ncol, lchnk )
1126 : !
1127 : ! jfl
1128 : !
1129 : ! surface vmr
1130 : !
1131 21888 : 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 21888 : 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 21888 : call outfld('Q_SRF',qh2o(:ncol,pver) , ncol, lchnk )
1164 21888 : call outfld('U_SRF',ufld(:ncol,pver) , ncol, lchnk )
1165 21888 : call outfld('V_SRF',vfld(:ncol,pver) , ncol, lchnk )
1166 :
1167 : !
1168 21888 : if (.not.sad_pbf_ndx>0) then
1169 21888 : deallocate(strato_sad)
1170 : endif
1171 :
1172 43776 : end subroutine gas_phase_chemdr
1173 :
1174 : end module mo_gas_phase_chemdr
|
