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