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