Line data Source code
1 :
2 : module mo_strato_rates
3 : !=======================================================================
4 : ! ROUTINE
5 : ! ratecon_sfstrat.f
6 : !
7 : ! Date...
8 : ! 15 August 2002
9 : ! 11 April 2008
10 : ! 15 December 2014
11 : !
12 : ! Programmed by...
13 : ! Douglas E. Kinnison
14 : !
15 : ! DESCRIPTION
16 : !
17 : ! Derivation of the rate constant for reactions on
18 : ! sulfate, NAT, and ICE aerosols.
19 : !
20 : !
21 : ! Sulfate Aerosol Reactions Rxn# Gamma
22 : ! N2O5 + H2O(l) => 2HNO3 (1) f(wt%)
23 : ! ClONO2 + H2O(l) => HOCl + HNO3 (2) f(T,P,HCl,H2O,r)
24 : ! BrONO2 + H2O(l) => HOBr + HNO3 (3) f(T,P,H2O,r)
25 : ! ClONO2 + HCl(l) => Cl2 + HNO3 (4) f(T,P,HCl,H2O,r)
26 : ! HOCl + HCl(l) => Cl2 + H2O (5) f(T,P,HCl,HOCl,H2O,r)
27 : ! HOBr + HCl(l) => BrCl + H2O (6) f(T,P,HCl,HOBr,H2O,r)
28 : !
29 : ! Nitric Acid Tri-hydrate Reactions Rxn# Gamma Reference
30 : ! N2O5 + H2O(s) => 2HNO3 (7) 4e-4 JPL10-6
31 : ! ClONO2 + H2O(s) => HOCl + HNO3 (8) 4e-3 JPL10-6
32 : ! ClONO2 + HCl(s) => Cl2 + HNO3 (9) 0.2 JPL10-6
33 : ! HOCl + HCl(s) => Cl2 + H2O (10) 0.1 JPL10-6
34 : ! BrONO2 + H2O(s) => HOBr + HNO3 (11) 0.006 Davies et JGR, 2003
35 : !
36 : ! WATER-ICE Aersol Reactions Rxn# Gamma
37 : ! N2O5 + H2O(s) => 2HNO3 (12) 0.02 JPL10-6
38 : ! ClONO2 + H2O(s) => HOCl + HNO3 (13) 0.3 JPL10-6
39 : ! BrONO2 + H2O(s) => HOBr + HNO3 (14) 0.3 JPL10-6
40 : ! ClONO2 + HCl(s) => Cl2 + HNO3 (15) 0.3 JPL10-6
41 : ! HOCl + HCl(s) => Cl2 + H2O (16) 0.2 JPL10-6
42 : ! HOBr + HCl(s) => BrCl + H2O (17) 0.3 JPL10-6
43 : !
44 : ! NOTE: The rate constants derived from species reacting with H2O are
45 : ! first order (i.e., sec-1 units) - an example is N2O5 + H2O = 2HNO3.
46 : ! Other reactions, e.g., ClONO2 + HCl have rate constants that
47 : ! are second order (i.e., cm+3 molecules-1 sec-1 units). In all
48 : ! of these types of reactions the derived first order rate constant
49 : ! {0.25*(mean Velocity)*SAD*gamma} is divided by the HCl abundance
50 : ! to derive the correct second order units.
51 : !
52 : ! NOTE: Liquid Sulfate Aerosols...
53 : ! See coding for references on how the Sulfate Aerosols were handled.
54 : ! Approach follows Shi et al., JGR, 106, D20, 24259, 2001.
55 : !
56 : !
57 : ! INPUT:
58 : ! ad . .... air density, molec. cm-3
59 : ! pmid ..... pressures, hPa
60 : ! temp ..... temperatures, K
61 : ! rad_sulfate ..... Surface area density, cm2 cm-3
62 : ! sad_sulfate ..... Surface area density, cm2 cm-3
63 : ! sad_nat ..... Surface area density, cm2 cm-3
64 : ! sad_ice ..... Surface area density, cm2 cm-3
65 : ! brono2mv ..... BrONO2 Volume Mixing Ratio
66 : ! clono2mvr ..... ClONO2 Volume Mixing Ratio
67 : ! h2omvr ..... H2O Volume Mixing Ratio
68 : ! hclmvr ..... HCl Volume Mixing Ratio
69 : ! hobrmvr ..... HOBr Volume Mixing Ratio
70 : ! hoclmvr ..... HOCl Volume Mixing Ratio
71 : ! n2o5mvr ..... N2O5 Volume Mixing Ratio
72 : !
73 : ! OUTPUT:
74 : !
75 : ! rxt ..... Rate constant (s-1 and cm3 sec-1 molec-1)
76 : !=======================================================================
77 :
78 : private
79 : public :: ratecon_sfstrat, init_strato_rates, has_strato_chem
80 :
81 : integer :: id_brono2, id_clono2, id_hcl, id_hocl, &
82 : id_hobr, id_n2o5
83 : integer :: rid_het1, rid_het2, rid_het3, rid_het4, rid_het5, &
84 : rid_het6, rid_het7, rid_het8, rid_het9, rid_het10, &
85 : rid_het11, rid_het12, rid_het13, rid_het14, rid_het15, &
86 : rid_het16, rid_het17
87 :
88 : logical :: has_strato_chem
89 :
90 : contains
91 :
92 1536 : subroutine init_strato_rates
93 :
94 : use mo_chem_utls, only : get_rxt_ndx, get_spc_ndx
95 : use mo_aero_settling, only : strat_aer_settl_init
96 : use ppgrid, only : pcols, pver
97 :
98 : implicit none
99 :
100 : integer :: ids(23)
101 :
102 1536 : rid_het1 = get_rxt_ndx( 'het1' )
103 1536 : rid_het2 = get_rxt_ndx( 'het2' )
104 1536 : rid_het3 = get_rxt_ndx( 'het3' )
105 1536 : rid_het4 = get_rxt_ndx( 'het4' )
106 1536 : rid_het5 = get_rxt_ndx( 'het5' )
107 1536 : rid_het6 = get_rxt_ndx( 'het6' )
108 1536 : rid_het7 = get_rxt_ndx( 'het7' )
109 1536 : rid_het8 = get_rxt_ndx( 'het8' )
110 1536 : rid_het9 = get_rxt_ndx( 'het9' )
111 1536 : rid_het10 = get_rxt_ndx( 'het10' )
112 1536 : rid_het11 = get_rxt_ndx( 'het11' )
113 1536 : rid_het12 = get_rxt_ndx( 'het12' )
114 1536 : rid_het13 = get_rxt_ndx( 'het13' )
115 1536 : rid_het14 = get_rxt_ndx( 'het14' )
116 1536 : rid_het15 = get_rxt_ndx( 'het15' )
117 1536 : rid_het16 = get_rxt_ndx( 'het16' )
118 1536 : rid_het17 = get_rxt_ndx( 'het17' )
119 :
120 1536 : id_brono2 = get_spc_ndx( 'BRONO2' )
121 1536 : id_clono2 = get_spc_ndx( 'CLONO2' )
122 1536 : id_hcl = get_spc_ndx( 'HCL' )
123 1536 : id_hocl = get_spc_ndx( 'HOCL' )
124 1536 : id_hobr = get_spc_ndx( 'HOBR' )
125 1536 : id_n2o5 = get_spc_ndx( 'N2O5' )
126 :
127 : ids(:) = (/ rid_het1, rid_het2, rid_het3, rid_het4, rid_het5, rid_het6, rid_het7, rid_het8, &
128 : rid_het9, rid_het10, rid_het11, rid_het12, rid_het13, rid_het14, rid_het15, &
129 36864 : rid_het16, rid_het17, id_brono2, id_clono2, id_hcl, id_hocl, id_hobr, id_n2o5 /)
130 :
131 1536 : has_strato_chem = all( ids(:) > 0 )
132 :
133 1536 : if (.not. has_strato_chem) return
134 :
135 0 : call strat_aer_settl_init
136 :
137 : endsubroutine init_strato_rates
138 :
139 0 : subroutine ratecon_sfstrat( ncol, ad, pmid, temp, rad_sulfate, sad_sulfate, &
140 0 : sad_nat, sad_ice, h2ovmr, vmr, rxt, &
141 0 : gprob_n2o5, gprob_cnt_hcl, gprob_cnt_h2o, gprob_bnt_h2o, &
142 0 : gprob_hocl_hcl, gprob_hobr_hcl, wtper )
143 :
144 : use shr_kind_mod, only : r8 => shr_kind_r8
145 : use chem_mods, only : adv_mass, rxntot, gas_pcnst
146 : use ppgrid, only : pcols, pver
147 : use mo_sad, only : sad_top
148 : use cam_logfile, only : iulog
149 :
150 : implicit none
151 :
152 : !-----------------------------------------------------------------------
153 : ! ... dummy arguments
154 : !-----------------------------------------------------------------------
155 : integer, intent(in) :: ncol ! columns in chunk
156 : real(r8), dimension(ncol,pver,gas_pcnst), intent(in) :: & ! species concentrations (mol/mol)
157 : vmr
158 : real(r8), dimension(ncol,pver), intent(in) :: &
159 : ad, & ! Air Density (molec. cm-3)
160 : rad_sulfate, & ! Radius of Sulfate Aerosol (cm)
161 : sad_ice, & ! ICE Surface Area Density (cm-1)
162 : sad_nat, & ! NAT Surface Area Density (cm-1)
163 : sad_sulfate, & ! Sulfate Surface Area Density (cm-1)
164 : h2ovmr ! water vapor volume mixing ratio( gas phase )
165 : real(r8), dimension(pcols,pver), intent(in) :: &
166 : pmid, & ! pressure (Pa)
167 : temp ! temperature (K)
168 :
169 : real(r8), intent(out) :: &
170 : rxt(ncol,pver,rxntot) ! rate constants
171 :
172 : real(r8), dimension(ncol,pver), intent(out) :: & ! diagnostics
173 : gprob_n2o5, &
174 : gprob_cnt_hcl, &
175 : gprob_cnt_h2o, &
176 : gprob_bnt_h2o, &
177 : gprob_hocl_hcl, &
178 : gprob_hobr_hcl, &
179 : wtper
180 :
181 : !-----------------------------------------------------------------------
182 : ! ... local variables
183 : !-----------------------------------------------------------------------
184 : real(r8), parameter :: small_div = 1.e-16_r8 ! for divid by excess species
185 : real(r8), parameter :: av_const = 2.117265e4_r8 ! (8*8.31448*1000 / PI)
186 : real(r8), parameter :: pa2mb = 1.e-2_r8 ! Pa to mb
187 : real(r8), parameter :: m2cm = 100._r8 ! meters to cms
188 :
189 : integer :: &
190 : i, & ! altitude loop index
191 : k, & ! level loop index
192 : m ! species index
193 :
194 : !-----------------------------------------------------------------------
195 : ! ... variables for gamma calculations
196 : !-----------------------------------------------------------------------
197 : real(r8) :: &
198 : brono2vmr, & ! BrONO2 Volume Mixing Ratio
199 : clono2vmr, & ! ClONO2 Volume Mixing Ratio
200 : hclvmr, & ! HCl Volume Mixing Ratio
201 : hcldeni, & ! inverse of HCl density
202 : cntdeni, & ! inverse of ClONO2 density
203 : hocldeni, & ! inverse of HOCl density
204 : hobrdeni, & ! inverse of HOBr density
205 : hoclvmr, & ! HOCl Volume Mixing Ratio
206 : hobrvmr, & ! HOBr Volume Mixing Ratio
207 : n2o5vmr ! N2O5 Volume Mixing Ratio
208 :
209 : real(r8) :: &
210 : av_n2o5, & ! N2O5 Mean Velocity (cm s-1)
211 : av_clono2, & ! ClONO2 Mean Velocity (cm s-1)
212 : av_brono2, & ! BrONO2Mean Velocity (cm s-1)
213 : av_hocl, & ! HOCl Mean Velocity (cm s-1)
214 : av_hobr ! HOBr Mean Velocity (cm s-1)
215 :
216 : real(r8) :: &
217 : pzero_h2o, & ! H2O sat vapor press (mbar)
218 : e0, e1, e2, e3, & ! coefficients for H2O sat vapor press.
219 : aw, & ! Water activity
220 : m_h2so4, & ! H2SO4 molality (mol/kg)
221 : wt, & ! wt % H2SO4
222 : y1, y2, & ! used in H2SO4 molality
223 : & a1, b1, c1, d1, & ! used in H2SO4 molality
224 : a2, b2, c2, d2 ! used in H2SO4 molality
225 :
226 : real(r8) :: &
227 : z1, z2, z3, & ! used in H2SO4 soln density
228 : den_h2so4, & ! H2SO4 soln density, g/cm3
229 : mol_h2so4, & ! Molality of H2SO4, mol / kg
230 : molar_h2so4, & ! Molarity of H2SO4, mol / l
231 : x_h2so4, & ! H2SO4 mole fraction
232 : aconst, tzero, & ! used in viscosity of H2SO4
233 : vis_h2so4, & ! H2SO4 viscosity
234 : ah, & ! Acid activity, molarity units
235 : term1,term2,term3,term4, & ! used in ah
236 : term5,term6,term7,term0, &
237 : T_limit, & ! temporary variable for temp (185-260K range)
238 : T_limiti, & ! 1./T_limit
239 : T_limitsq, & ! sqrt( T_limit )
240 : rad_sulf, & ! temporary variable for sulfate radius (cm)
241 : sadsulf, & ! temporary variable for sulfate radius (cm)
242 : sadice, & ! temporary variable for sulfate radius (cm)
243 : sadnat ! temporary variable for sulfate radius (cm)
244 :
245 : real(r8) :: &
246 : C_cnt, S_cnt, & ! used in H_cnt
247 : H_cnt, & ! Henry's law coeff. for ClONO2
248 : H_hcl, & ! Henry's law coeff. for HCl
249 : D_cnt, &
250 : k_hydr, &
251 : k_h2o, &
252 : k_h, &
253 : k_hcl, &
254 : rdl_cnt, &
255 : f_cnt, &
256 : M_hcl, &
257 : atmos
258 :
259 : real(r8) :: &
260 : Gamma_b_h2o, &
261 : Gamma_cnt_rxn, &
262 : Gamma_b_hcl, &
263 : Gamma_s, &
264 : Fhcl, &
265 : Gamma_s_prime, &
266 : Gamma_b_hcl_prime, &
267 : Gamma_b, &
268 : gprob_rxn, &
269 : gprob_tot, &
270 : gprob_cnt
271 :
272 : real(r8) :: &
273 : D_hocl, &
274 : k_hocl_hcl, &
275 : C_hocl, &
276 : S_hocl, &
277 : H_hocl, &
278 : Gamma_hocl_rxn, &
279 : rdl_hocl, &
280 : f_hocl
281 :
282 : real(r8) :: &
283 : h1, h2, h3, &
284 : alpha
285 :
286 : real(r8) :: &
287 : C_hobr, &
288 : D_hobr, &
289 : aa, bb, cc, dd, &
290 : k_hobr_hcl, &
291 : k_dl, &
292 : k_wasch, &
293 : H_hobr, &
294 : rdl_hobr, &
295 : Gamma_hobr_rxn, &
296 : f_hobr
297 :
298 : real(r8) :: &
299 : pmb,& ! Pressure, mbar (hPa)
300 : pH2O_atm,& ! Partial press. H2O (atm)
301 : pH2O_hPa,& ! Partial press. H2O (hPa)
302 : pHCl_atm,& ! Partial press. HCl (atm)
303 : pCNT_atm ! Partial press. ClONO2 (atm)
304 :
305 : !-----------------------------------------------------------------------
306 : ! ... Used in pzero h2o calculation
307 : !-----------------------------------------------------------------------
308 : real(r8), parameter :: wt_e0 = 18.452406985_r8
309 : real(r8), parameter :: wt_e1 = 3505.1578807_r8
310 : real(r8), parameter :: wt_e2 = 330918.55082_r8
311 : real(r8), parameter :: wt_e3 = 12725068.262_r8
312 :
313 : real(r8) :: &
314 : wrk, tmp
315 :
316 :
317 :
318 0 : if (.not. has_strato_chem) return
319 :
320 : !-----------------------------------------------------------------------
321 : ! ... intialize rate constants
322 : !-----------------------------------------------------------------------
323 0 : do k = 1,pver
324 0 : rxt(:,k,rid_het1) = 0._r8
325 0 : rxt(:,k,rid_het2) = 0._r8
326 0 : rxt(:,k,rid_het3) = 0._r8
327 0 : rxt(:,k,rid_het4) = 0._r8
328 0 : rxt(:,k,rid_het5) = 0._r8
329 0 : rxt(:,k,rid_het6) = 0._r8
330 0 : rxt(:,k,rid_het7) = 0._r8
331 0 : rxt(:,k,rid_het8) = 0._r8
332 0 : rxt(:,k,rid_het9) = 0._r8
333 0 : rxt(:,k,rid_het10) = 0._r8
334 0 : rxt(:,k,rid_het11) = 0._r8
335 0 : rxt(:,k,rid_het12) = 0._r8
336 0 : rxt(:,k,rid_het13) = 0._r8
337 0 : rxt(:,k,rid_het14) = 0._r8
338 0 : rxt(:,k,rid_het15) = 0._r8
339 0 : rxt(:,k,rid_het16) = 0._r8
340 0 : rxt(:,k,rid_het17) = 0._r8
341 :
342 0 : gprob_n2o5(:,k) = 0._r8
343 0 : gprob_cnt_h2o(:,k) = 0._r8
344 0 : gprob_cnt_hcl(:,k) = 0._r8
345 0 : gprob_bnt_h2o(:,k) = 0._r8
346 0 : gprob_hocl_hcl(:,k)= 0._r8
347 0 : gprob_hobr_hcl(:,k)= 0._r8
348 0 : wtper(:,k) = 0._r8
349 : end do
350 :
351 : !-----------------------------------------------------------------------
352 : ! ... set rate constants
353 : !-----------------------------------------------------------------------
354 : Level_loop : &
355 0 : do k = sad_top+1,pver
356 : column_loop : &
357 0 : do i = 1,ncol
358 : !-----------------------------------------------------------------------
359 : ! ... set species, pmb, and atmos
360 : !-----------------------------------------------------------------------
361 0 : brono2vmr = vmr(i,k,id_brono2)
362 0 : clono2vmr = vmr(i,k,id_clono2)
363 0 : hclvmr = vmr(i,k,id_hcl)
364 0 : hoclvmr = vmr(i,k,id_hocl)
365 0 : hobrvmr = vmr(i,k,id_hobr)
366 0 : if( hclvmr > 0._r8 ) then
367 0 : hcldeni = 1._r8/(hclvmr*ad(i,k))
368 : end if
369 0 : if( clono2vmr > 0._r8 ) then
370 0 : cntdeni = 1._r8/(clono2vmr*ad(i,k))
371 : end if
372 0 : if( hoclvmr > 0._r8 ) then
373 0 : hocldeni = 1._r8/(hoclvmr*ad(i,k))
374 : end if
375 0 : if( hobrvmr > 0._r8 ) then
376 0 : hobrdeni = 1._r8/(hobrvmr*ad(i,k))
377 : end if
378 0 : n2o5vmr = vmr(i,k,id_n2o5)
379 0 : sadsulf = sad_sulfate(i,k)
380 0 : sadnat = sad_nat(i,k)
381 0 : sadice = sad_ice(i,k)
382 0 : pmb = pa2mb*pmid(i,k)
383 0 : atmos = pmb/1013.25_r8
384 :
385 : !-----------------------------------------------------------------------
386 : ! ... setup for stratospheric aerosols
387 : ! data range set: 185K - 240K; Tabazedeh GRL, 24, 1931, 1997
388 : !-----------------------------------------------------------------------
389 0 : T_limit = max( temp(i,k),185._r8 )
390 0 : T_limit = min( T_limit,240._r8 )
391 0 : T_limiti = 1._r8/T_limit
392 0 : T_limitsq = sqrt( T_limit )
393 :
394 : !-----------------------------------------------------------------------
395 : ! .... Average velocity (8RT*1000/(PI*MW))**1/2 * 100.(units cm s-1)
396 : ! .... or (av_const*T/M2)**1/2
397 : !-----------------------------------------------------------------------
398 0 : wrk = av_const*T_limit
399 0 : av_n2o5 = sqrt( wrk/adv_mass(id_n2o5) )*m2cm
400 0 : av_clono2 = sqrt( wrk/adv_mass(id_clono2) )*m2cm
401 0 : av_brono2 = sqrt( wrk/adv_mass(id_brono2) )*m2cm
402 0 : av_hocl = sqrt( wrk/adv_mass(id_hocl) )*m2cm
403 0 : av_hobr = sqrt( wrk/adv_mass(id_hobr) )*m2cm
404 : has_sadsulf : &
405 0 : if( sadsulf > 0._r8 ) then
406 : !-----------------------------------------------------------------------
407 : ! .... Partial Pressure of H2O, ClONO2, and HCl in atmospheres
408 : !-----------------------------------------------------------------------
409 0 : if( hclvmr > 0._r8 ) then
410 0 : pHCl_atm = hclvmr*atmos
411 : else
412 : pHCl_atm = 0._r8
413 : end if
414 :
415 0 : if( clono2vmr > 0._r8 ) then
416 0 : pCNT_atm = clono2vmr*atmos
417 : else
418 : pCNT_atm = 0._r8
419 : end if
420 :
421 0 : if( h2ovmr(i,k) > 0._r8 ) then
422 : pH2O_atm = h2ovmr(i,k)*atmos
423 : else
424 : pH2O_atm = 0._r8
425 : end if
426 : !-----------------------------------------------------------------------
427 : ! .... Partial Pressure of H2O in hPa
428 : !-----------------------------------------------------------------------
429 0 : pH2O_hpa = h2ovmr(i,k)*pmb
430 : !-----------------------------------------------------------------------
431 : ! .... Calculate the h2so4 Wt% and Activity of H2O -
432 : !-----------------------------------------------------------------------
433 : !-----------------------------------------------------------------------
434 : ! ... Saturation Water Vapor Pressure (mbar)
435 : !-----------------------------------------------------------------------
436 0 : pzero_h2o = exp( wt_e0 - T_limiti*(wt_e1 + T_limiti*(wt_e2 - T_limiti*wt_e3)) )
437 : !-----------------------------------------------------------------------
438 : ! ... H2O activity
439 : ! ... if the activity of H2O goes above 1.0, wt% can go negative
440 : !-----------------------------------------------------------------------
441 0 : aw = ph2o_hpa / pzero_h2o
442 0 : aw = min( aw,1._r8 )
443 0 : aw = max( aw,.001_r8 )
444 : !-----------------------------------------------------------------------
445 : ! ... h2so4 Molality (mol/kg)
446 : !-----------------------------------------------------------------------
447 0 : if( aw <= .05_r8 ) then
448 : a1 = 12.37208932_r8
449 : b1 = -0.16125516114_r8
450 : c1 = -30.490657554_r8
451 : d1 = -2.1133114241_r8
452 : a2 = 13.455394705_r8
453 : b2 = -0.1921312255_r8
454 : c2 = -34.285174607_r8
455 : d2 = -1.7620073078_r8
456 0 : else if( aw > .05_r8 .and. aw < .85_r8 ) then
457 : a1 = 11.820654354_r8
458 : b1 = -0.20786404244_r8
459 : c1 = -4.807306373_r8
460 : d1 = -5.1727540348_r8
461 : a2 = 12.891938068_r8
462 : b2 = -0.23233847708_r8
463 : c2 = -6.4261237757_r8
464 : d2 = -4.9005471319_r8
465 : else
466 0 : a1 = -180.06541028_r8
467 0 : b1 = -0.38601102592_r8
468 0 : c1 = -93.317846778_r8
469 0 : d1 = 273.88132245_r8
470 0 : a2 = -176.95814097_r8
471 0 : b2 = -0.36257048154_r8
472 0 : c2 = -90.469744201_r8
473 0 : d2 = 267.45509988_r8
474 : end if
475 : !-----------------------------------------------------------------------
476 : ! ... h2so4 mole fraction
477 : !-----------------------------------------------------------------------
478 0 : y1 = a1*(aw**b1) + c1*aw + d1
479 0 : y2 = a2*(aw**b2) + c2*aw + d2
480 0 : m_h2so4 = y1 + ((T_limit - 190._r8)*(y2 - y1)) / 70._r8
481 : !-----------------------------------------------------------------------
482 : ! ... h2so4 Weight Percent
483 : !-----------------------------------------------------------------------
484 0 : wt = 9800._r8*m_h2so4 / (98._r8*m_h2so4 + 1000._r8)
485 0 : wtper(i,k) = wt
486 : !-----------------------------------------------------------------------
487 : ! .... Parameters for h2so4 Solution
488 : !-----------------------------------------------------------------------
489 : ! ... h2so4 Solution Density (g/cm3)
490 : !-----------------------------------------------------------------------
491 0 : wrk = T_limit*T_limit
492 0 : z1 = .12364_r8 - 5.6e-7_r8*wrk
493 0 : z2 = -.02954_r8 + 1.814e-7_r8*wrk
494 0 : z3 = 2.343e-3_r8 - T_limit*1.487e-6_r8 - 1.324e-8_r8*wrk
495 : !-----------------------------------------------------------------------
496 : ! ... where mol_h2so4 is molality in mol/kg
497 : !-----------------------------------------------------------------------
498 0 : den_h2so4 = 1._r8 + m_h2so4*(z1 + z2*sqrt(m_h2so4) + z3*m_h2so4)
499 : !-----------------------------------------------------------------------
500 : ! ... h2so4 Molarity, mol / l
501 : !-----------------------------------------------------------------------
502 0 : molar_h2so4 = den_h2so4*wt/9.8_r8
503 : !-----------------------------------------------------------------------
504 : ! ... h2so4 Mole fraction
505 : !-----------------------------------------------------------------------
506 0 : x_h2so4 = wt / (wt + ((100._r8 - wt)*98._r8/18._r8))
507 0 : term1 = .094_r8 - x_h2so4*(.61_r8 - 1.2_r8*x_h2so4)
508 0 : term2 = (8515._r8 - 10718._r8*(x_h2so4**.7_r8))*T_limiti
509 0 : H_hcl = term1 * exp( -8.68_r8 + term2 )
510 0 : M_hcl = H_hcl*pHCl_atm
511 : !-----------------------------------------------------------------------
512 : ! ... h2so4 solution viscosity
513 : !-----------------------------------------------------------------------
514 0 : aconst = 169.5_r8 + wt*(5.18_r8 - wt*(.0825_r8 - 3.27e-3_r8*wt))
515 0 : tzero = 144.11_r8 + wt*(.166_r8 - wt*(.015_r8 - 2.18e-4_r8*wt))
516 0 : vis_h2so4 = aconst/(T_limit**1.43_r8) * exp( 448._r8/(T_limit - tzero) )
517 : !-----------------------------------------------------------------------
518 : ! ... Acid activity in molarity
519 : !-----------------------------------------------------------------------
520 0 : term1 = 60.51_r8
521 0 : term2 = .095_r8*wt
522 0 : wrk = wt*wt
523 0 : term3 = .0077_r8*wrk
524 0 : term4 = 1.61e-5_r8*wt*wrk
525 0 : term5 = (1.76_r8 + 2.52e-4_r8*wrk) * T_limitsq
526 0 : term6 = -805.89_r8 + (253.05_r8*(wt**.076_r8))
527 0 : term7 = T_limitsq
528 0 : ah = exp( term1 - term2 + term3 - term4 - term5 + term6/term7 )
529 0 : if( ah <= 0._r8 ) then
530 0 : write(iulog,*) 'ratecon: ah <= 0 at i,k, = ',i,k
531 0 : write(iulog,*) 'ratecon: term1,term2,term3,term4,term5,term6,term7,wt,T_limit,ah = ', &
532 0 : term1,term2,term3,term4,term5,term6,term7,wt,T_limit,ah
533 : end if
534 :
535 0 : wrk = .25_r8*sadsulf
536 0 : rad_sulf = max( rad_sulfate(i,k),1.e-6_r8 )
537 : !-----------------------------------------------------------------------
538 : ! N2O5 + H2O(liq) => 2.00*HNO3 Sulfate Aerosol Reaction
539 : !-----------------------------------------------------------------------
540 0 : term0 = -25.5265_r8 - wt*(.133188_r8 - wt*(.00930846_r8 - 9.0194e-5_r8*wt))
541 0 : term1 = 9283.76_r8 + wt*(115.345_r8 - wt*(5.19258_r8 - .0483464_r8*wt))
542 0 : term2 = -851801._r8 - wt*(22191.2_r8 - wt*(766.916_r8 - 6.85427_r8*wt))
543 0 : gprob_n2o5(i,k) = exp( term0 + T_limiti*(term1 + term2*T_limiti) )
544 0 : rxt(i,k,rid_het1) = max( 0._r8,wrk*av_n2o5*gprob_n2o5(i,k) )
545 :
546 : !-----------------------------------------------------------------------
547 : ! ClONO2 + H2O(liq) = HOCl + HNO3 Sulfate Aerosol Reaction
548 : !-----------------------------------------------------------------------
549 : ! ... NOTE: Aerosol radius in units of cm.
550 : !-----------------------------------------------------------------------
551 : !-----------------------------------------------------------------------
552 : ! ... Radius sulfate set (from sad module)
553 : ! Set min radius to 0.01 microns (1e-6 cm)
554 : ! Typical radius is 0.1 microns (1e-5 cm)
555 : ! f_cnt may go negative under if not set.
556 : !-----------------------------------------------------------------------
557 0 : C_cnt = 1474._r8*T_limitsq
558 0 : S_cnt = .306_r8 + 24._r8*T_limiti
559 0 : term1 = exp( -S_cnt*molar_h2so4 )
560 0 : H_cnt = 1.6e-6_r8 * exp( 4710._r8*T_limiti )*term1
561 0 : D_cnt = 5.e-8_r8*T_limit / vis_h2so4
562 0 : k_h = 1.22e12_r8*exp( -6200._r8*T_limiti )
563 0 : k_h2o = 1.95e10_r8*exp( -2800._r8*T_limiti )
564 0 : k_hydr = (k_h2o + k_h*ah)*aw
565 0 : k_hcl = 7.9e11_r8*ah*D_cnt*M_hcl
566 0 : rdl_cnt = sqrt( D_cnt/(k_hydr + k_hcl) )
567 0 : term1 = 1._r8/tanh( rad_sulf/rdl_cnt )
568 0 : term2 = rdl_cnt/rad_sulf
569 0 : f_cnt = term1 - term2
570 0 : if( f_cnt > 0._r8 ) then
571 0 : term1 = 4._r8*H_cnt*.082_r8*T_limit
572 0 : term2 = sqrt( D_cnt*k_hydr )
573 0 : Gamma_b_h2o = term1*term2/C_cnt
574 0 : term1 = sqrt( 1._r8 + k_hcl/k_hydr )
575 0 : Gamma_cnt_rxn = f_cnt*Gamma_b_h2o*term1
576 0 : Gamma_b_hcl = Gamma_cnt_rxn*k_hcl/(k_hcl + k_hydr)
577 0 : term1 = exp( -1374._r8*T_limiti )
578 0 : Gamma_s = 66.12_r8*H_cnt*M_hcl*term1
579 0 : if( pHCl_atm > 0._r8 ) then
580 0 : term1 = .612_r8*(Gamma_s+Gamma_b_hcl)* pCNT_atm/pHCl_atm
581 0 : Fhcl = 1._r8/(1._r8 + term1)
582 : else
583 : Fhcl = 1._r8
584 : end if
585 0 : Gamma_s_prime = Fhcl*Gamma_s
586 0 : Gamma_b_hcl_prime = Fhcl*Gamma_b_hcl
587 0 : term1 = Gamma_cnt_rxn*k_hydr
588 0 : term2 = k_hcl + k_hydr
589 0 : Gamma_b = Gamma_b_hcl_prime + (term1/term2)
590 0 : term1 = 1._r8 / (Gamma_s_prime + Gamma_b)
591 0 : gprob_cnt = 1._r8 / (1._r8 + term1)
592 0 : term1 = Gamma_s_prime + Gamma_b_hcl_prime
593 0 : term2 = Gamma_s_prime + Gamma_b
594 0 : gprob_cnt_hcl(i,k) = gprob_cnt * term1/term2
595 0 : gprob_cnt_h2o(i,k) = gprob_cnt - gprob_cnt_hcl(i,k)
596 : else
597 0 : gprob_cnt_h2o(i,k) = 0._r8
598 0 : gprob_cnt_hcl(i,k) = 0._r8
599 0 : Fhcl = 1._r8
600 : end if
601 :
602 0 : rxt(i,k,rid_het2) = max( 0._r8,wrk*av_clono2*gprob_cnt_h2o(i,k) )
603 :
604 : !-----------------------------------------------------------------------
605 : ! ... BrONO2 + H2O(liq) = HOBr + HNO3 Sulfate Aerosol Reaction
606 : !-----------------------------------------------------------------------
607 0 : h1 = 29.24_r8
608 0 : h2 = -.396_r8
609 0 : h3 = .114_r8
610 0 : alpha = .805_r8
611 0 : gprob_rxn = exp( h1 + h2*wt ) + h3
612 0 : term1 = 1._r8/alpha
613 0 : term2 = 1._r8/gprob_rxn
614 0 : gprob_bnt_h2o(i,k) = 1._r8 / (term1 + term2)
615 0 : rxt(i,k,rid_het3) = max( 0._r8,wrk*av_brono2*gprob_bnt_h2o(i,k) )
616 :
617 : !-----------------------------------------------------------------------
618 : ! ... ClONO2 + HCl(liq) = Cl2 + HNO3 Sulfate Aerosol Reaction
619 : !-----------------------------------------------------------------------
620 0 : if( hclvmr > small_div .and. clono2vmr > small_div ) then
621 0 : if ( hclvmr > clono2vmr ) then
622 0 : rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*hcldeni
623 : else
624 0 : rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*cntdeni
625 : end if
626 : end if
627 :
628 : !-----------------------------------------------------------------------
629 : ! ... HOCl + HCl(liq) = Cl2 + H2O Sulfate Aerosol Reaction
630 : !-----------------------------------------------------------------------
631 : !-----------------------------------------------------------------------
632 : ! ... Radius sulfate set (from sad module)
633 : ! Set min radius to 0.01 microns (1e-6 cm)
634 : ! Typical radius is 0.1 microns (1e-5 cm)
635 : ! f_hocl may go negative under if not set.
636 : !-----------------------------------------------------------------------
637 0 : if( pCNT_atm > 0._r8 ) then
638 0 : D_hocl = 6.4e-8_r8*T_limit/vis_h2so4
639 0 : k_hocl_hcl = 1.25e9_r8*ah*D_hocl*M_hcl
640 0 : C_hocl = 2009._r8*T_limitsq
641 0 : S_hocl = .0776_r8 + 59.18_r8*T_limiti
642 0 : term1 = exp( -S_hocl*molar_h2so4 )
643 0 : H_hocl = 1.91e-6_r8 * exp( 5862.4_r8*T_limiti )*term1
644 0 : term1 = 4._r8*H_hocl*.082_r8*T_limit
645 0 : term2 = sqrt( D_hocl*k_hocl_hcl )
646 0 : Gamma_hocl_rxn = term1*term2/C_hocl
647 0 : rdl_hocl = sqrt( D_hocl/k_hocl_hcl )
648 0 : term1 = 1._r8/tanh( rad_sulf/rdl_hocl )
649 0 : term2 = rdl_hocl/rad_sulf
650 0 : f_hocl = term1 - term2
651 0 : if( f_hocl > 0._r8 ) then
652 0 : term1 = 1._r8 / (f_hocl*Gamma_hocl_rxn*Fhcl)
653 0 : gprob_hocl_hcl(i,k) = 1._r8 / (1._r8 + term1)
654 : else
655 0 : gprob_hocl_hcl(i,k) = 0._r8
656 : end if
657 :
658 0 : if( hclvmr > small_div .and. hoclvmr > small_div ) then
659 0 : if ( hclvmr > hoclvmr ) then
660 0 : rxt(i,k,rid_het5) = max( 0._r8,wrk*av_hocl*gprob_hocl_hcl(i,k) )*hcldeni
661 : else
662 0 : rxt(i,k,rid_het5) = max( 0._r8,wrk*av_hocl*gprob_hocl_hcl(i,k) )*hocldeni
663 : end if
664 : end if
665 : end if
666 :
667 : !-----------------------------------------------------------------------
668 : ! ... HOBr + HCl(liq) = BrCl + H2O Sulfate Aerosol Reaction
669 : !-----------------------------------------------------------------------
670 : !-----------------------------------------------------------------------
671 : ! ... Radius sulfate set (from sad module)
672 : ! Set min radius to 0.01 microns (1e-6 cm)
673 : ! Typical radius is 0.1 microns (1e-5 cm)
674 : ! f_hobr may go negative under if not set.
675 : !-----------------------------------------------------------------------
676 0 : C_hobr = 1477._r8*T_limitsq
677 0 : D_hobr = 9.e-9_r8
678 : !-----------------------------------------------------------------------
679 : ! ... Taken from Waschewsky and Abbat
680 : ! Dave Hanson (PC) suggested we divide this rc by eight to agee
681 : ! with his data (Hanson, 108, D8, 4239, JGR, 2003).
682 : ! k1=k2*Mhcl for gamma(HOBr)
683 : !-----------------------------------------------------------------------
684 0 : k_wasch = .125_r8 * exp( .542_r8*wt - 6440._r8*T_limiti + 10.3_r8)
685 : !-----------------------------------------------------------------------
686 : ! ... Taken from Hanson 2002.
687 : !-----------------------------------------------------------------------
688 0 : H_hobr = exp( -9.86_r8 + 5427._r8*T_limiti )
689 0 : k_dl = 7.5e14_r8*D_hobr*2._r8 ! or 7.5e14*D *(2nm)
690 : !-----------------------------------------------------------------------
691 : ! ... If k_wasch is GE than the diffusion limit...
692 : !-----------------------------------------------------------------------
693 0 : if( M_hcl > 0._r8 ) then
694 0 : if( k_wasch >= k_dl ) then
695 0 : k_hobr_hcl = k_dl * M_hcl
696 : else
697 0 : k_hobr_hcl = k_wasch * M_hcl
698 : end if
699 0 : term1 = 4._r8*H_hobr*.082_r8*T_limit
700 0 : term2 = sqrt( D_hobr*k_hobr_hcl )
701 0 : tmp = rad_sulf/term2
702 0 : Gamma_hobr_rxn = term1*term2/C_hobr
703 0 : rdl_hobr = sqrt( D_hobr/k_hobr_hcl )
704 0 : if( tmp < 1.e2_r8 ) then
705 0 : term1 = 1._r8/tanh( rad_sulf/rdl_hobr )
706 : else
707 0 : term1 = 1._r8
708 : end if
709 0 : term2 = rdl_hobr/rad_sulf
710 0 : f_hobr = term1 - term2
711 0 : if( f_hobr > 0._r8 ) then
712 0 : term1 = 1._r8 / (f_hobr*Gamma_hobr_rxn)
713 0 : gprob_hobr_hcl(i,k)= 1._r8 / (1._r8 + term1)
714 : else
715 0 : gprob_hobr_hcl(i,k)= 0._r8
716 : end if
717 0 : if( hclvmr > small_div .and. hobrvmr > small_div ) then
718 0 : if ( hclvmr > hobrvmr ) then
719 0 : rxt(i,k,rid_het6) = max( 0._r8,wrk*av_hobr*gprob_hobr_hcl(i,k) )*hcldeni
720 : else
721 0 : rxt(i,k,rid_het6) = max( 0._r8,wrk*av_hobr*gprob_hobr_hcl(i,k) )*hobrdeni
722 : end if
723 : end if
724 : end if
725 :
726 : end if has_sadsulf
727 :
728 : has_sadnat : &
729 0 : if( sadnat > 0._r8 ) then
730 0 : wrk = .25_r8*sadnat
731 : !-----------------------------------------------------------------------
732 : ! ... N2O5 + H2O(s) => 2HNO3 NAT Aerosol Reaction
733 : !-----------------------------------------------------------------------
734 : !-----------------------------------------------------------------------
735 : ! ... gprob based on JPL10-6 for NAT.
736 : ! also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
737 : ! gprob_tot = 4.e-4
738 : !-----------------------------------------------------------------------
739 0 : rxt(i,k,rid_het7) = wrk*av_n2o5*4.e-4_r8
740 :
741 : !-----------------------------------------------------------------------
742 : ! ClONO2 + H2O(s) => HNO3 + HOCl NAT Aerosol Reaction
743 : !-----------------------------------------------------------------------
744 : !-----------------------------------------------------------------------
745 : ! ... gprob based on JPL10-6 for NAT.
746 : ! also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
747 : ! gprob_tot = 0.004
748 : !-----------------------------------------------------------------------
749 0 : rxt(i,k,rid_het8) = wrk*av_clono2*4.0e-3_r8
750 :
751 : !-----------------------------------------------------------------------
752 : ! ... ClONO2 + HCl(s) => HNO3 + Cl2, NAT Aerosol Reaction
753 : !-----------------------------------------------------------------------
754 : !-----------------------------------------------------------------------
755 : ! ... gprob based on JPL10-6 for NAT.
756 : ! also see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
757 : ! gprob_tot = 0.2
758 : !-----------------------------------------------------------------------
759 0 : if( hclvmr > small_div .and. clono2vmr > small_div ) then
760 0 : if ( hclvmr > clono2vmr ) then
761 0 : rxt(i,k,rid_het9) = wrk*av_clono2*0.2_r8*hcldeni
762 : else
763 0 : rxt(i,k,rid_het9) = wrk*av_clono2*0.2_r8*cntdeni
764 : end if
765 : end if
766 :
767 : !-----------------------------------------------------------------------
768 : ! ... HOCl + HCl(s) => H2O + Cl2 NAT Aerosol Reaction
769 : !-----------------------------------------------------------------------
770 : !-----------------------------------------------------------------------
771 : ! ... gprob based on JPL10-6 for NAT.
772 : ! see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
773 : ! and Abbatt and Molina, GRL, 19, 461-464, 1992.
774 : ! gprob_tot = 0.1
775 : !-----------------------------------------------------------------------
776 0 : if( hclvmr > small_div .and. hoclvmr > small_div ) then
777 0 : if ( hclvmr > hoclvmr ) then
778 0 : rxt(i,k,rid_het10) = wrk*av_hocl*0.1_r8*hcldeni
779 : else
780 0 : rxt(i,k,rid_het10) = wrk*av_hocl*0.1_r8*hocldeni
781 : end if
782 : end if
783 :
784 : !-----------------------------------------------------------------------
785 : ! ... BrONO2 + H2O(s) => HOBr + HNO3 NAT Aerosol Reaction
786 : !-----------------------------------------------------------------------
787 : !-----------------------------------------------------------------------
788 : ! ... Davies et al.,
789 : ! JGR, 108, NO. D5, 8322, doi:10.1029/2001JD000445, 2003
790 : ! gprob_tot = 0.006
791 : !-----------------------------------------------------------------------
792 0 : rxt(i,k,rid_het11) = wrk*av_brono2*0.006_r8
793 :
794 : end if has_sadnat
795 :
796 : has_sadice : &
797 0 : if( sadice > 0._r8 ) then
798 0 : wrk = .25_r8*sadice
799 : !-----------------------------------------------------------------------
800 : ! N2O5 + H2O(s) => 2HNO3 ICE Aerosol Reaction
801 : !-----------------------------------------------------------------------
802 : !-----------------------------------------------------------------------
803 : ! ... gprob based on JPL10-6 for ICE.
804 : ! also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
805 : ! gprob_tot = 0.02
806 : !-----------------------------------------------------------------------
807 0 : rxt(i,k,rid_het12) = wrk*av_n2o5*0.02_r8
808 :
809 : !-----------------------------------------------------------------------
810 : ! ... ClONO2 + H2O(s) => HNO3 + HOCl ICE Aerosol Reaction
811 : !-----------------------------------------------------------------------
812 : !-----------------------------------------------------------------------
813 : ! ... gprob based on JPL10-6 for ICE.
814 : ! also see Hanson and Ravi, JGR, 96, 17307-17314, 1991.
815 : ! gprob_tot = 0.3
816 : !-----------------------------------------------------------------------
817 0 : rxt(i,k,rid_het13) = wrk*av_clono2*0.3_r8
818 :
819 : !-----------------------------------------------------------------------
820 : ! ... BrONO2 + H2O(s) => HNO3 + HOBr ICE Aerosol Reaction
821 : !-----------------------------------------------------------------------
822 : !-----------------------------------------------------------------------
823 : ! ... gprob based on JPL10-6 for ICE.
824 : ! also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
825 : ! could be as high as 1.0
826 : ! gprob_tot = 0.3
827 : !-----------------------------------------------------------------------
828 0 : rxt(i,k,rid_het14) = wrk*av_brono2*0.3_r8
829 :
830 : !-----------------------------------------------------------------------
831 : ! ClONO2 + HCl(s) => HNO3 + Cl2, ICE Aerosol Reaction
832 : !-----------------------------------------------------------------------
833 : !-----------------------------------------------------------------------
834 : ! ... gprob based on JPL10-6 for ICE.
835 : ! also see Hanson and Ravi, GRL, 15, 17-20, 1988.
836 : ! also see Lue et al.,
837 : ! gprob_tot = 0.3
838 : !-----------------------------------------------------------------------
839 0 : if( hclvmr > small_div .and. clono2vmr > small_div ) then
840 0 : if ( hclvmr > clono2vmr ) then
841 0 : rxt(i,k,rid_het15) = wrk*av_clono2*0.3_r8*hcldeni
842 : else
843 0 : rxt(i,k,rid_het15) = wrk*av_clono2*0.3_r8*cntdeni
844 : end if
845 : end if
846 : !
847 : !-----------------------------------------------------------------------
848 : ! ... HOCl + HCl(s) => H2O + Cl2, ICE Aerosol Reaction
849 : !-----------------------------------------------------------------------
850 : !-----------------------------------------------------------------------
851 : ! ... gprob based on JPL10-6 for ICE.
852 : ! also see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
853 : ! also see Abbatt and Molina, GRL, 19, 461-464, 1992.
854 : ! gprob_tot = 0.2
855 : !-----------------------------------------------------------------------
856 0 : if( hoclvmr > small_div .and. hclvmr > small_div ) then
857 0 : if ( hclvmr > hoclvmr ) then
858 0 : rxt(i,k,rid_het16) = wrk*av_hocl*0.2_r8*hcldeni
859 : else
860 0 : rxt(i,k,rid_het16) = wrk*av_hocl*0.2_r8*hocldeni
861 : end if
862 : end if
863 :
864 : !-----------------------------------------------------------------------
865 : ! HOBr + HCl(s) => H2O + BrCl, ICE Aerosol Reaction
866 : !-----------------------------------------------------------------------
867 : !-----------------------------------------------------------------------
868 : ! ... gprob based on JPL10-6 for ICE.
869 : ! Abbatt GRL, 21, 665-668, 1994.
870 : ! gprob_tot = 0.3
871 : !-----------------------------------------------------------------------
872 0 : if( hobrvmr > small_div .and. hclvmr > small_div ) then
873 0 : if ( hclvmr > hobrvmr ) then
874 0 : rxt(i,k,rid_het17) = wrk*av_hobr*0.3_r8*hcldeni
875 : else
876 0 : rxt(i,k,rid_het17) = wrk*av_hobr*0.3_r8*hobrdeni
877 : end if
878 : end if
879 :
880 : end if has_sadice
881 : end do column_loop
882 : end do Level_loop
883 :
884 : end subroutine ratecon_sfstrat
885 :
886 : end module mo_strato_rates
|
