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 |
|