/[MITgcm]/MITgcm_contrib/jscott/igsm/src_chem/chemtrop.F
ViewVC logotype

Contents of /MITgcm_contrib/jscott/igsm/src_chem/chemtrop.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Thu Sep 17 17:40:33 2009 UTC (15 years, 10 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
Error occurred while calculating annotation data.
chem module archive

1
2 #include "ctrparam.h"
3
4 ! ============================================================
5 !
6 ! CHEMTROP.F: Subroutine for calculating chemical
7 ! reactions in troposphere
8 ! of MIT Global Chemistry Model
9 !
10 ! ------------------------------------------------------------
11 !
12 ! Author: Chien Wang
13 ! MIT Joint Program on Science and Policy
14 ! of Global Change
15 !
16 ! ----------------------------------------------------------
17 !
18 ! Revision History:
19 !
20 ! When Who What
21 ! ---- ---------- -------
22 ! 081898 Chien Wang rev.
23 ! 080200 Chien Wang repack based on CliChem3 & add cpp
24 !
25 ! ==========================================================
26
27 !
28 subroutine chemtrop (dtr, ifss, ktrop, Temp, qv, den)
29 ! =====================================================
30
31 ! ==========================================================
32 ! Brief Description of Dummy Variables:
33 !
34 ! dtr: time step of integration in second
35 ! ifss: 1 for using steady-state assumption
36 ! 0 for not using
37 ! ktrop: layer label of tropopause
38 ! Temp: temperature in Kelvin
39 ! p: air pressure in hPa
40 ! qv: water vapor density ratio in kg/m^3
41 ! den: air density in kg/m^3
42 !
43 ! ==========================================================
44
45 #include "chem_para"
46 #include "chem_const1"
47 #include "chem_com"
48
49 dimension Temp (nlon,nlat,nlev)
50 dimension qv (nlon,nlat,nlev)
51 dimension den (nlon,nlat,nlev)
52
53 dimension ychem (nvaria)
54 dimension awchem (nvaria)
55 dimension revawchem(nvaria)
56 dimension ychem0 (2)
57
58 real ychem_hfc134a
59
60 ! ------------------------------------------------------------
61
62 #if ( defined CPL_CHEM )
63
64 awchem(1) = awO3
65 awchem(2) = awCO
66 awchem(3) = awCO2
67 awchem(4) = awHO
68 awchem(5) = awHO2
69 awchem(6) = awH2O2
70 awchem(7) = awNO
71 awchem(8) = awNO2
72 awchem(9) = awCH4
73 awchem(10)= awCH2O
74 awchem(11)= awCH3O2H
75 awchem(12)= awSO2
76 awchem(13)= awN2O5
77 awchem(14)= awHNO3
78 awchem(15)= awH2SO4
79
80 ! write(6,*)"Here are solarflux & coszangle"
81 ! write(6,*)solarflux
82 ! write(6,*)coszangle
83 ! write(6,*)
84
85 do ind = 1,15
86 revawchem(ind) = 1./awchem(ind)
87 enddo
88
89 yyy = Avogadro*1.e-3
90
91 yyy1 = yyy/28.964 !dry air
92 yyy2 = yyy/18.0152 !water vapor
93
94 do 1 k=1,ktrop
95 patm = airpress(k)/1013.27
96 do 2 j=1,n2dh
97
98 Rflux = solarflux(1,j,k)
99 & * 0.5678 !surface ratio of
100 ! Act(<= 800nm)/Act(total)
101 ! WMO 1981
102
103 if(Rflux.le.50.0) goto 2 !no photochemistry
104 !for nighttime
105
106 ! === Convert air density from kg/m^3 to mols/cm^3:
107
108 ychem0(1) = den(1,j,k)*yyy1 !dry air
109 ychem0(2) = qv (1,j,k)*yyy2 !water vapor
110
111 ! === Convert mixing ratio from ppb(m) to mols/cm^3:
112
113 ychem( 1) = o3 (1,j,k)
114 ychem( 2) = co (1,j,k)
115 ychem( 3) = zco2 (1,j,k)
116 ychem( 4) = ho (1,j,k)
117 ychem( 5) = ho2 (1,j,k)
118 ychem( 6) = h2o2 (1,j,k)
119 ychem( 7) = xno (1,j,k)
120 ychem( 8) = xno2 (1,j,k)
121 ychem( 9) = ch4 (1,j,k)
122 ychem(10) = ch2o (1,j,k)
123 ychem(11) = ch3o2h(1,j,k)
124 ychem(12) = so2 (1,j,k)
125
126 ychem(13) = xn2o5 (1,j,k)
127 ychem(14) = hno3 (1,j,k)
128 ychem(15) = h2so4 (1,j,k)
129
130 do iii = 1,nvaria
131 ychem(iii) = max(0.0, ychem(iii))
132 enddo
133
134 call ppbm2mcm (nvaria, ychem, den(1,j,k), revawchem)
135
136 #ifdef INC_3GASES
137 ! === if hfc, pfc, and sf6 are included:
138 ! ===
139 ! === 033098 add hfc134a:
140 ychem_hfc134a = hfc134a (1,j,k)
141
142 revawchemhfc134a = 1./awHFC134a
143 call ppbm2mcm(1,ychem_hfc134a,den(1,j,k),
144 & revawchemhfc134a)
145 #endif
146
147 ! === Calculate reaction rates:
148
149 ntemp = nint(max(1.0,
150 & min(300.0,(Temp(1,j,k) - 200.)*2.0)))
151
152 ! ===
153 ! === 042596 add zenith angle function:
154
155 call rateph (Temp(1,j,k), ntemp, Rflux, coszangle(1,j))
156
157 call ratebm (Temp(1,j,k), ntemp, patm)
158
159 call ratetm (Temp(1,j,k), ntemp, ychem0)
160
161 ! if(ifss.eq.1) call chemsteady(nvaria, ychem0,ychem)
162
163 ! =====
164 ! ===== New scheme for calculating NOy and S(VI)
165 ! ===== Chien Wang, 092195
166 ! =====
167
168 !
169 ! === weighting d[N2O5] and d[HNO3]:
170 !
171 xxxn2o5 = rk(14)*ychem(1) ! k14[O3]
172 xxxhno3 = rk(13)*ychem(4) ! k13[HO]
173 xxxtotal= 2*xxxn2o5 + xxxhno3 ! [NO2] = 1./
174 ! (2xxxn2o5 + xxxhno3)
175 !
176 ! === save initial NOx and SO2:
177 !
178 ynox0 = ychem(7) + ychem(8)
179 yso20 = ychem(12)
180
181 ! -------------------------------------------
182 ! call tropreact1(dtr, nvaria, ychem0, ychem)
183 call tropreact1(dtr, 12, ychem0, ychem)
184 ! -------------------------------------------
185
186 !
187 ! === Derive [N2O5], [HNO3], and [H2SO4]; 092195:
188 !
189 dnox = min(0.0, ychem(7) + ychem(8) - ynox0)
190 if(xxxtotal.gt.0.0)then
191 ychem(13) = ychem(13) - xxxn2o5/xxxtotal * dnox
192 ychem(14) = ychem(14) - xxxhno3/xxxtotal * dnox
193 endif
194
195 ychem(15) = ychem(15) - min(0.0, ychem(12) - yso20)
196
197 #ifdef INC_3GASES
198 ! === if hfc, pfc, and sf6 are included:
199 ! ===
200 ! === 033098 add hfc134a oxidation: CF3CH2F + OH -> CF3CHF + H2O
201 ! === here derive exp{-k[OH]t} first:
202
203 xxxhfc134a = exp(- 1.7e-12
204 & * exp(-1750.0/Temp(1,j,k))
205 & * dtr*ychem(4))
206 ychem_hfc134a = max(0.0, ychem_hfc134a*xxxhfc134a)
207
208 call mcm2ppbm (1,ychem_hfc134a,den(1,j,k),awHFC134a)
209
210 hfc134a(1,j,k)= ychem_hfc134a
211 #endif
212
213 ! === Convert concentration to ppb(m):
214
215 call mcm2ppbm (nvaria, ychem, den(1,j,k), awchem)
216
217 do 10 iii = 1,nvaria
218 ychem(iii) = max(0.0, ychem(iii))
219 10 continue
220
221 ! === Give value back to the common block:
222
223 o3 (1,j,k) = ychem(01)
224 co (1,j,k) = ychem(02)
225 zco2 (1,j,k) = ychem(03)
226 ho (1,j,k) = ychem(04)
227 ho2 (1,j,k) = ychem(05)
228 h2o2 (1,j,k) = ychem(06)
229 xno (1,j,k) = ychem(07)
230 xno2 (1,j,k) = ychem(08)
231 ch4 (1,j,k) = ychem(09)
232 ch2o (1,j,k) = ychem(10)
233 ch3o2h(1,j,k) = ychem(11)
234 so2 (1,j,k) = ychem(12)
235 xn2o5 (1,j,k) = ychem(13)
236 hno3 (1,j,k) = ychem(14)
237 h2so4 (1,j,k) = ychem(15)
238
239 2 continue
240 1 continue
241
242 #endif
243
244 return
245 end
246
247 !
248 ! --- won't need this portion for run excludes chemistry
249 !
250 #if ( defined CPL_CHEM )
251
252 c
253 subroutine tropreact1 (tout, neq, c00, y)
254 c ==============================================
255
256 c==================================================================c
257 c c
258 c TROPREACT1.F: Subroutine for calculating chemical c
259 c reactions in troposphere c
260 c of MIT Global Chemistry Model c
261 c partly based on LSODES c
262 c ------------------------------------------------- c
263 c Author: Chien Wang c
264 c MIT Joint Program on Science and Policy c
265 c of Global Change c
266 c Last Revised on: January 31, 1995 c
267 c c
268 c==================================================================c
269
270 external fex, jex
271
272 dimension y(neq), rwork(872), iwork(45)
273 dimension c00(2)
274
275 common /c00tmp/c00tmp(2)
276
277 c data lrw/500/, liw/30/
278
279 c00tmp(1) = c00(1)
280 c00tmp(2) = c00(2)
281
282 t = 0.
283 itol = 1
284 rtol = 1.e-3
285 atol = 1.e+3 !molecules/cm3
286
287 c atol(1) = 1.e-6
288 c atol(2) = 1.e-10
289 c atol(3) = 1.e-6
290
291 itask = 1
292 istate = 1
293 iopt = 0
294 lrw = 872
295 liw = 45
296 mf = 21
297
298 call lsodenew(fex,neq,y,t,tout,itol,rtol,atol,itask,istate,
299 & iopt,rwork,lrw,iwork,liw,jex,mf)
300 c call lsode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate,
301 c & iopt,rwork,lrw,iwork,liw,jex,mf)
302
303 return
304 end
305
306 c
307 subroutine fex (neq, t, y, ydot)
308 c ================================
309 c
310 c Subroutine for defining of ydot
311 c
312
313 #include "chem_para"
314 #include "chem_com"
315
316 common /c00tmp/c00(2)
317
318 dimension y(neq), ydot(neq)
319
320 bbb = 0.7809*c00(1)
321 ccc = 0.2095*c00(1)
322
323 c === redesigned scheme, 092195
324 c
325
326 xbeta1 = rk(3)*bbb + rk(4)*ccc !k3[N2] + k4[O2]
327 xbeta2 = rk(2)*c00(2) !k2[H2O]
328 xbeta = xbeta2/(xbeta1 + xbeta2)
329
330 rr1 = rk(1) *y(1) !j1[O3]
331 rr5 = rk(5) *y(2) *y(4) !k5[CO][HO]
332 rr7 = rk(7) *y(5) *y(7) !k7[HO2][NO]
333 rr8 = rk(8) *y(8) !k8[NO2]
334 rr10= rk(10)*y(1) *y(5) !k10[O3][HO2]
335 rr11= rk(11)*y(1) *y(4) !k11[O3][HO]
336 rr12= rk(12)*y(1) *y(7) !k12[O3][NO]
337 rr13= rk(13)*y(4) *y(8) !k13[HO][NO2]
338 rr14= rk(14)*y(1) *y(8) !k14[O3][NO2]
339 rr16= rk(16)*y(5) *y(5) !k16[HO2][HO2]
340 rr17= rk(17)*y(6) !j17[H2O2]
341 rr18= rk(18)*y(4) *y(6) !k18[HO][H2O2]
342 rr19= rk(19)*y(4) *y(9) !k19[HO][CH4]
343 rr21= rk(21)*y(7) !k21[NO] - for ss only
344 rr23= rk(23)*y(5) !k23[HO2]- for ss only
345 rr24= rk(24)*y(11) !j24[CH3O2H]
346 rr25= rk(25)*y(4) *y(11) !k25[HO][CH3O2H]
347 rr26= rk(26)*y(10) !j26[CH2O]
348 rr27= rk(27)*y(4) *y(10) !k27[HO][CH2O]
349 rr29= rk(29)*y(4) *y(12) !k29[HO][SO2]
350
351 c New reactions - 062995:
352
353 rr32= rk(32)*y(4)*y(5) !k32[HO][HO2]
354 rr33= rk(33)*y(4)*y(4) !k33[HO][HO]
355 rr34= rk(34)*y(4)*y(4) !k34[HO][HO]
356
357 xxx = (rr21 + rr23)
358
359 gamax = 0.0
360 if(rr21.gt.1.e-12)gamax = rr21/xxx
361
362 gamay = 0.0
363 if(xxx.gt.0.0)gamay = rr23/xxx
364
365 ydot(1) = rr8 + rr33 !O3
366 & -(xbeta*rr1 + rr10 + rr11 + rr12 + rr14)
367
368 ydot(2) = rr26 + rr27 !CO
369 & - rr5
370
371 ydot(3) = rr5 !CO2
372
373 ydot(4) = 2.0*xbeta*rr1 + rr7 + rr10 !HO
374 & + 2.0*rr17+ rr24
375 & -(rr5 + rr11 + rr13 + rr18 + rr19
376 & +rr25+ rr27 + rr29
377 & +rr32+ rr33 + rr34 )
378
379 ydot(5) = rr5 + rr11 + rr18 + rr27 + rr29 !HO2
380 & + rr24+ 2.0*rr26
381 & +(2.0*gamax - 1.0)*(rr19 + rr25)
382 & -(rr7 + rr10 + rr16 + rr32)
383
384 ydot(6) = rr16 + rr34 !H2O2
385 & -(rr17+ rr18)
386
387 ydot(7) = rr8 !NO
388 & -(rr7 + rr12 + gamax*(rr19 + rr25))
389
390 ydot(8) = rr7 + rr12 + gamax*(rr19 + rr25) !NO2
391 & -(rr8 + rr13+ rr14+ rr14)
392
393 ydot(9) =-rr19 !CH4
394
395 ydot(10) = rr24 + gamax*(rr19 + rr25) !CH2O
396 & -(rr26 + rr27)
397
398 ydot(11) = gamay*(rr19 + rr25) !CH3O2H
399 & -(rr24 + rr25)
400
401 ydot(12) =-rr29 !SO2
402
403 return
404 end
405
406 c
407 subroutine jex (neq, t, y, ml, mu, pd, nrpd)
408 c ==========================================
409 c
410 c Subroutine for calculating the Jacobian for f(i)
411 c
412 c P(m,n) here represents df(m)/dn
413 c or d ydot(m)/dn
414 c
415
416 #include "chem_para"
417 #include "chem_com"
418
419 common /c00tmp/c00(2)
420
421 dimension y(neq), pd(nrpd,neq)
422
423 bbb = 0.7809*c00(1)
424 ccc = 0.2095*c00(1)
425
426 xbeta1 = rk(3)*bbb + rk(4)*ccc !k3[N2] + k4[O2]
427 xbeta2 = rk(2)*c00(2) !k2[H2O]
428 xbeta = xbeta2/(xbeta1 + xbeta2)
429
430 rr21= rk(21)*y(7) !k21[NO] - for ss only
431 rr23= rk(23)*y(5) !k23[HO2]- for ss only
432 xxx = (rr21 + rr23)
433
434 gamax = 0.0
435 if(rr21.gt.1.e-12)gamax = rr21/xxx
436
437 gamay = 0.0
438 if(xxx.gt.0.0)gamay = rr23/xxx
439
440 do 1 j=1,neq
441 do 1 i=1,neq
442 pd(i,j) = 0.0
443 1 continue
444
445 c =====
446 c ===== pd(i,j) = d ydot(i)/dy(j):
447 c =====
448
449 rk5y2 = rk(5)*y(2)
450 rk5y4 = rk(5)*y(4)
451 rk7y5 = rk(7)*y(5)
452 rk7y7 = rk(7)*y(7)
453 rk10y1 = rk(10)*y(1)
454 rk10y5 = rk(10)*y(5)
455 rk11y1 = rk(11)*y(1)
456 rk11y4 = rk(11)*y(4)
457 rk12y1 = rk(12)*y(1)
458 rk12y7 = rk(12)*y(7)
459 rk13y4 = rk(13)*y(4)
460 rk13y8 = rk(13)*y(8)
461 rk14y1 = rk(14)*y(1)
462 rk14y8 = rk(14)*y(8)
463 rk18y4 = rk(18)*y(4)
464 rk18y6 = rk(18)*y(6)
465 rk19y4 = rk(19)*y(4)
466 rk19y9 = rk(19)*y(9)
467 rk25y4 = rk(25)*y(4)
468 rk25y11 = rk(25)*y(11)
469 rk27y4 = rk(27)*y(4)
470 rk27y10 = rk(27)*y(10)
471 rk29y4 = rk(29)*y(4)
472 rk29y12 = rk(29)*y(12)
473 rk32y4 = rk(32)*y(4)
474 rk32y5 = rk(32)*y(5)
475
476 pd(1,1) = -(xbeta*rk(1) + rk10y5 + rk11y4
477 & +rk12y7 + rk14y8)
478 pd(4,1) = 2.0*xbeta*rk(1) + rk10y5 - rk11y4
479 pd(5,1) = rk11y4 - rk10y5
480 pd(8,1) = rk12y7 - 2.*rk14y8
481
482 pd(2,2) =-rk5y4
483 pd(3,2) = rk5y4
484 pd(4,2) =-rk5y4
485 pd(5,2) = rk5y4
486
487 pd(1,4) = 2.0*rk(33)*y(4) - rk11y1
488 pd(2,4) = rk27y10 - rk5y2
489 pd(3,4) = rk5y2
490 pd(4,4) =-(
491 & rk5y2 + rk11y1 + rk(13)*y(8)
492 & + rk18y6 + rk19y9 + rk25y11
493 & + rk27y10+ rk29y12+ rk32y5
494 & +(rk(33) + rk(34))*y(4)
495 & )
496 pd(5,4) =(rk5y2 + rk11y1 + rk18y6
497 & + rk27y10+ rk29y12)
498 & +(2.0*gamax - 1.0)*(rk19y9 + rk25y11)
499 & - rk32y5
500 pd(6,4) =-rk18y6
501 pd(7,4) =-gamax*(rk19y9 + rk25y11)
502 pd(8,4) = gamax*(rk19y9 + rk25y11)
503 & - rk13y8
504 pd(9,4) =-rk19y9
505 pd(10,4) = gamax*(rk19y9 + rk25y11)
506 & - rk27y10
507 pd(11,4) = gamay*(rk19y9 + rk25y11)
508 & - rk25y11
509 pd(12,4) =-rk29y12
510
511 pd(1,5) = rk11y1
512 pd(4,5) = rk7y7 + rk10y1
513 & - rk32y4
514 pd(5,5) =
515 & -(rk7y7 + rk10y1 + rk(16)*y(5) + rk32y4)
516 pd(6,5) = 2.0* (rk(16) + rk(34))*y(5)
517 pd(7,5) =-rk7y7
518 pd(8,5) = rk7y7
519 c pd(11,5)
520
521 pd(4,6) = 2.0*rk(17)
522 & - rk18y4
523 pd(5,6) = rk18y4
524 pd(6,6) =-(rk(17) + rk18y4)
525
526 pd(1,7) =-rk12y1
527 pd(4,7) = rk7y5
528 pd(5,7) =-rk7y5
529 pd(7,7) =-(rk7y5 + rk12y1)
530 pd(8,7) = (rk7y5 + rk12y1)
531
532 pd(1,8) = rk(8)
533 & - rk14y1
534 pd(4,8) =-rk13y4
535 pd(7,8) = rk(8)
536 pd(8,8) =-(rk(8) + rk13y4 + 2.0*rk14y1)
537
538 pd(4,9) =-rk19y4
539 pd(5,9) =(2.0*gamax - 1.0)*rk19y4
540 pd(7,9) =-gamax*rk19y4
541 pd(8,9) = gamax*rk19y4
542 pd(9,9) =-rk19y4
543 pd(10,9) = gamax*rk19y4
544 pd(11,9) = gamay*rk19y4
545
546 pd(2,10) = rk(26) + rk27y4
547 pd(4,10) =-rk27y4
548 pd(5,10) = rk27y4 + 2.0*rk(26)
549 pd(10,10)=-(rk(26)+ rk27y4)
550
551 pd(4,11) = rk(24) - rk25y4
552 pd(5,11) = rk(24) + (2.0*gamax - 1.0)*rk25y4
553 pd(7,11) =-gamax*rk25y4
554 pd(8,11) = gamax*rk25y4
555 pd(10,11)= rk(24) + gamax*rk25y4
556 pd(11,11)=-(rk(24) + rk25y4*max(0.0,1.0 - gamay))
557
558 pd(4,12) =-rk(29)*y(4)
559 pd(5,12) = rk29y4
560 pd(12,12)=-rk29y4
561
562 return
563 end
564
565
566 c==================================================================c
567 c c
568 c CHEMFUNCPACK.F: Subroutine and function pack c
569 c of MIT Global Chemistry Model c
570 c ------------------------------------------------- c
571 c Author: Chien Wang c
572 c MIT Joint Program on Science and Policy c
573 c of Global Change c
574 c Last Revised on: February 22, 1995 c
575 c c
576 c==================================================================c
577
578 c
579 subroutine ppbm2mcm (np, xx, den, revaw)
580 c ==========================================
581
582 c------------------------------------------------
583 c A program for convert mixing ratio in
584 c ppb(m) to concentration in molecules/cm^3
585 c------------------------------------------------
586
587 dimension xx (np)
588 dimension revaw (np)
589
590 ddd = 6.02217e+11*den
591 do 1 i=1,np
592 xx (i) = xx (i)
593 & *ddd*revaw(i)
594 1 continue
595
596 return
597 end
598
599 c
600 subroutine mcm2ppbm (np, xx, den, aw)
601 c ==========================================
602
603 c------------------------------------------------
604 c A program for convert concentration in
605 c molecules/cm^3 to mixing ratio in ppb(m)
606 c------------------------------------------------
607
608 dimension xx (np)
609 dimension aw (np)
610
611 ddd = 1./(6.02217e+11*den)
612 do 1 i=1,np
613 xx (i) = xx (i)
614 & *ddd*aw(i)
615 1 continue
616
617 return
618 end
619
620 c
621 subroutine rateph (T, ntemp, Rflux, coszagcm)
622 c ===========================
623
624 c------------------------------------------------
625 c A program for calculating rate constants
626 c of photochemical reactions
627 c ----------------------------
628 c The time unit for all reaction rates is second
629 c
630 c New version revised at June 3, 1995
631 c
632 c------------------------------------------------
633
634 #include "chem_para"
635 #include "chem_com"
636
637 coszagcm = max(cosza4rk(10),
638 & min(cosza4rk(1),coszagcm))
639
640 do 1 nband=1,9
641 if(coszagcm.le.cosza4rk(nband).and.
642 & coszagcm.gt.cosza4rk(nband+1))then
643
644 wgtxx = 1.0/(cosza4rk(nband) - cosza4rk(nband+1))
645 wgt1 = (cosza4rk(nband) - coszagcm) *wgtxx
646 wgt2 = (coszagcm - cosza4rk(nband+1))*wgtxx
647 c
648 c +--- wgt2 ---+--- wgt1 ---+
649 c nband --- coszagcm --- nband+1
650
651 c-- R8: NO2 + hv -> NO + O
652 c Atkinson et al., 1992:
653 c
654 c rk(8) = (18096.7314929
655 c & -3.9704*(T-273.))
656 c & *5.0e-10
657 c & *Rflux
658
659 c rk(8) = rktable1(08,ntemp)*Rflux
660
661 rk(8) = ((rk08gama(nband)*wgt2 + rk08gama(nband+1)*wgt1)
662 & + (rk08aaa (nband)*wgt2 + rk08aaa (nband+1)*wgt1)
663 & *(T-273.))
664 & *Rflux
665
666 c-- R1: O3 + hv -> O(1D) + O2
667 c Atkinson et al., 1992:
668 c
669 c rk(1) = 0.0028*rk(8)
670 cc rk(1) = 3.467625419939578E-08*Rflux
671
672 c rk(1) = (rk01table(nband)*wgt2
673 c & + rk01table(nband+1)*wgt1)*Rflux
674 rk(1) = 0.0028*rk(8)
675
676 c-- R17: H2O2 + hv -> 2OH
677 c Atkinson et al., 1992:
678 c
679 c rk(17)= 6.263464249748237E-09*Rflux
680
681 rk(17) = (rk17table(nband)*wgt2
682 & + rk17table(nband+1)*wgt1)*Rflux
683
684 c-- R24: CH3O2H + hv -> CH3O + HO
685 c Atkinson et al., 1992:
686 c
687 c rk(24)= 4.474208962739173E-09*Rflux
688
689 rk(24) = (rk24table(nband)*wgt2
690 & + rk24table(nband+1)*wgt1)*Rflux
691
692 c-- R26: CH2O + hv -> CHO + H
693 c Atkinson et al., 1992:
694 c
695 c rk(26)= 9.410107812688823E-08*Rflux
696
697 rk(26) = (rk26table(nband)*wgt2
698 & + rk26table(nband+1)*wgt1)*Rflux
699
700 goto 2
701 endif
702
703 1 continue
704
705 c = For znith angle.ge.86:
706
707 rk(8) = (rk08gama(10)
708 & + rk08aaa (10)
709 & *(T-273.))
710 & *Rflux
711
712 c rk(1) = rk01table(10)
713 c & * Rflux
714 rk(1) = 0.0028*rk(8)
715
716 rk(17) = rk17table(10)
717 & * Rflux
718
719 rk(24) = rk24table(10)
720 & * Rflux
721
722 rk(26) = rk26table(10)
723 & *Rflux
724
725 2 continue
726
727 return
728 end
729
730 c
731 subroutine ratebm (T, ntemp, patm)
732 c ===========================
733
734 c------------------------------------------------
735 c A program for calculating rate constants
736 c of bimolecular reactions
737 c ----------------------------
738 c The time unit for all reaction rates is second
739 c
740 c------------------------------------------------
741
742 #include "chem_para"
743 #include "chem_com"
744
745 Trev = 1./T
746
747 c-- R2: O(1D) + H2O -> 2OH
748
749 rk(2) = 2.2e-10
750
751 c-- R3: O(1D) + N2 -> O + N2
752
753 c rk(3) = 1.8e-11 * exp( 107.0*Trev)
754 rk(3) = rktable1(03,ntemp)
755
756 c-- R4: O(1D) + O2 -> O + O2
757
758 c rk(4) = 3.2e-11 * exp( 67.0*Trev)
759 rk(4) = rktable1(04,ntemp)
760
761 c-- R5: CO + HO -> H + CO2
762
763 rk(5) = 1.5e-13 * (1.0 + 0.6*patm)
764
765 c-- R7: HO2 + NO -> HO + NO2
766
767 c rk(7) = 3.7e-12 * exp( 240.0*Trev)
768 rk(7) = rktable1(07,ntemp)
769
770 c-- R10: HO2 + O3 -> HO + 2O2
771
772 c rk(10) = 1.1e-14 * exp(-500.0*Trev)
773 rk(10) = rktable1(10,ntemp)
774
775 c-- R11: HO + O3 -> HO2 + O2
776
777 c rk(11) = 1.6e-12 * exp(-940.0*Trev)
778 rk(11) = rktable1(11,ntemp)
779 c rk(11) = 0.0
780
781 c-- R12: NO + O3 -> NO2 + O2
782
783 c rk(12) = 2.0e-12 * exp(-1400.0*Trev)
784 rk(12) = rktable1(12,ntemp)
785
786 c-- R14: NO2 + O3 -> NO3 + O2
787
788 c rk(14) = 1.2e-13 * exp(-2450.0*Trev)
789 rk(14) = rktable1(14,ntemp)
790
791 c-- R16: HO2 + HO2 -> H2O2 + O2
792
793 c rk(16) = 2.3e-13 * exp( 600.0*Trev)
794
795 c-- R18: H2O2 + HO -> HO2 + H2O
796
797 c rk(18) = 2.9e-12 * exp(-160.0*Trev)
798 rk(18) = rktable1(18,ntemp)
799
800 c-- R19: CH4 + HO -> CH3 + H2O
801
802 c rk(19) = 2.65e-12 * exp(-1800.0*Trev)
803 rk(19) = rktable1(19,ntemp)
804
805 c-- R21: CH3O2 + NO -> CH3O + NO2
806
807 c rk(21) = 4.2e-12 * exp( 180.0*Trev)
808 rk(21) = rktable1(21,ntemp)
809
810 c-- R22: CH3O + O2 -> CH2O + HO2
811
812 c rk(22) = 3.9e-14 * exp(-900.0*Trev)
813 rk(22) = rktable1(22,ntemp)
814
815 c-- R23: CH3O2 + HO2 -> CH3O2H + O2
816
817 c rk(23) = 3.8e-13 * exp( 780.0*Trev)
818 rk(23) = rktable1(23,ntemp)
819
820 c-- R25: CH3O2H + HO -> CH3O2 + H2O
821
822 c rk(25) = 1.9e-12 * exp( 190.0*Trev)
823 rk(25) = rktable1(25,ntemp)
824
825 c-- R27: CH2O + HO -> CHO + H2O
826
827 rk(27) = 1.0e-11
828
829 c-- R28: CHO + O2 -> CO + HO2
830
831 c rk(28) = 3.5e-12 * exp( 140.0*Trev)
832 rk(28) = rktable1(28,ntemp)
833
834 c-- R30: HOSO2 + O2 -> HO2 + SO3
835
836 c rk(30) = 1.3e-12 * exp(-330.0*Trev)
837 rk(30) = rktable1(30,ntemp)
838
839 c-- R31: SO3 + H2O -> H2SO4
840
841 rk(31) = 2.4e-15
842
843 c New reactions - 062995:
844
845 c-- R32: HO + HO2 -> H2O + O2
846
847 c rk(32) = 4.8e-11 * exp( 250.0*Trev)
848 rk(32) = rktable1(32,ntemp)
849
850 c-- R33: HO + HO -> H2O + O
851
852 c rk(33) = 4.2e-12 * exp(-240.0*Trev)
853 rk(33) = rktable1(33,ntemp)
854
855 return
856 end
857
858 c
859 subroutine ratetm (T, ntemp, ym)
860 c =========================
861
862 c------------------------------------------------
863 c A program for calculating rate constants
864 c of termolecular reactions
865 c ----------------------------
866 c The time unit for all reaction rates is second
867 c
868 c------------------------------------------------
869
870 #include "chem_para"
871 #include "chem_com"
872
873 dimension ym(2)
874
875 T300 = T/300.
876 xm = ym(1)
877 xn2 = 0.7809*xm
878 h2o = ym(2)
879
880 c-- R6: H + O2 + M -> HO2 + M
881
882 c rkt0 = 6.2e-32 * xm * T300**(-1.6)
883 rkt0 = rktable1(06,ntemp) * xm
884 rkt00 = 7.5e-11
885 rk(6) = calkmt(rkt0, rkt00)
886
887 c-- R9: O + O2 + M -> O3 + M (no high pres. limit)
888
889 c rkt0 = 5.6e-34 * xm * T300**(-2.8)
890 rkt0 = rktable1(09,ntemp) * xm
891 rk(9) = rkt0
892
893 c-- R13: NO2 + HO + M -> HNO3 + M
894
895 c rkt0 = 2.6e-30 * xm * T300**(-3.2)
896 c rkt00 = 2.4e-11 * T300**(-1.3)
897 rkt0 = rktable1(13,ntemp) * xm
898 rkt00 = rktable2(1, ntemp)
899 rk(13) = calkmt(rkt0, rkt00)
900
901 c-- R15: NO3 + NO2 + M -> N2O5 + M
902
903 c rkt0 = 2.2e-30 * xm * T300**(-3.9)
904 c rkt00 = 1.5e-12 * T300**(-0.7)
905 rkt0 = rktable1(15,ntemp) * xm
906 rkt00 = rktable2(2, ntemp)
907 rk(15) = calkmt(rkt0, rkt00)
908
909 c-- R16: HO2 + HO2 + M -> H2O2 + O2
910 c 062895:
911
912 c rk(16) = 1.7e-33 * xn2 * exp(1000.0/T)
913 rk(16) = rktable1(16,ntemp) * xn2
914 & *(1.+ 1.4e-21*h2o*exp(2200.0/T))
915
916 c-- R20: CH3 + O2 + M -> CH3O2 + M
917
918 c rkt0 = 4.5e-31 * xm * T300**(-3.0)
919 c rkt00 = 1.8e-12 * T300**(-1.7)
920 rkt0 = rktable1(20,ntemp) * xm
921 rkt00 = rktable2(3, ntemp)
922 rk(20) = calkmt(rkt0, rkt00)
923
924 c-- R29: SO2 + HO + M -> HOSO2 + M
925
926 c rkt0 = 3.0e-31 * xm * T300**(-3.3)
927 rkt0 = rktable1(29,ntemp) * xm
928 rkt00 = 1.5e-12
929 rk(29) = calkmt(rkt0, rkt00)
930
931 c New reactions, 063095:
932
933 c-- R34: HO + HO + M -> H2O2 + M
934
935 c rkt0 = 6.9e-31 * xm * T300**(-0.8)
936 rkt0 = rktable1(34,ntemp) * xm
937 rkt00 = 1.5e-11
938 rk(34) = calkmt(rkt0, rkt00)
939
940
941 return
942 end
943
944 c
945 function calkmt (rkt0, rkt00)
946 c =============================
947
948 c
949 c A function for calculate k(M,T) of termolecular
950 c reaction rate from low and high limit
951 c
952
953 aaa = rkt0/rkt00
954 bbb = log10(aaa)**2
955 ccc = 1./(1.+bbb)
956
957 calkmt = rkt0/(1.+aaa) * 0.6**ccc
958
959 return
960 end
961
962
963 c
964 Block data rktable3dat
965 c ======================
966
967 c----------------------------------------c
968 c Block data for holding rktable3 data c
969 c----------------------------------------c
970
971 #include "chem_para"
972 #include "chem_com"
973
974 data cosza4rk/
975 & 0.100000E+01, 0.984808E+00, 0.939693E+00, 0.866025E+00,
976 & 0.766044E+00, 0.642788E+00, 0.500000E+00, 0.342020E+00,
977 & 0.207912E+00, 0.697565E-01/
978
979 data rk08gama/
980 & 0.958260E-05, 0.956540E-05, 0.947786E-05, 0.933214E-05,
981 & 0.904837E-05, 0.853986E-05, 0.790570E-05, 0.644646E-05,
982 & 0.540693E-05, 0.640000E-05/
983
984 data rk08aaa/
985 & -0.213357E-08, -0.212833E-08, -0.210322E-08, -0.206239E-08,
986 & -0.198520E-08, -0.185147E-08, -0.169861E-08, -0.133796E-08,
987 & -0.109912E-08, -0.131294E-08/
988
989 data rk01table/
990 & 0.553973E-07, 0.541612E-07, 0.501300E-07, 0.436050E-07,
991 & 0.346763E-07, 0.241176E-07, 0.141663E-07, 0.593145E-08,
992 & 0.238321E-08, 0.148548E-08/
993
994 data rk17table/
995 & 0.761430E-08, 0.755268E-08, 0.730677E-08, 0.691454E-08,
996 & 0.626346E-08, 0.536943E-08, 0.430112E-08, 0.288139E-08,
997 & 0.199044E-08, 0.209279E-08/
998
999 data rk24table/
1000 & 0.535338E-08, 0.531281E-08, 0.515418E-08, 0.489785E-08,
1001 & 0.447421E-08, 0.387878E-08, 0.318768E-08, 0.217960E-08,
1002 & 0.155503E-08, 0.170283E-08/
1003
1004 data rk26table/
1005 & 0.109106E-06, 0.108473E-06, 0.105771E-06, 0.101587E-06,
1006 & 0.941011E-07, 0.831775E-07, 0.701839E-07, 0.494384E-07,
1007 & 0.364853E-07, 0.410051E-07/
1008
1009 end
1010
1011 #endif
1012

  ViewVC Help
Powered by ViewVC 1.1.22