module_cbmz.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! Chemistry Option: CBMZ (Carbon Bond Mechanism IV - Zaveri)
8 ! * Primary investigator: Rahul A. Zaveri
9 ! * Co-investigator: Richard C. Easter, William I. Gustafson Jr.
10 ! Last update: September 2005
11 !
12 ! Contacts:
13 ! Rahul A. Zaveri, PhD Jerome D. Fast, PhD
14 ! Senior Research Scientist Staff Scientist
15 ! Pacific Northwest National Laboratory Pacific Northwest National Laboratory
16 ! P.O. Box 999, MSIN K9-30 P.O. Box 999, MSIN K9-30
17 ! Richland, WA 99352 Richland, WA, 99352
18 ! Phone: (509) 372-6159 Phone: (509) 372-6116
19 ! Email: Rahul.Zaveri@pnl.gov Email: Jerome.Fast@pnl.gov
20 !
21 ! Please report any bugs or problems to Rahul Zaveri, the primary author of the
22 ! code, or Jerome Fast, the WRF-chem implementation team leader
23 !
24 !Terms of Use:
25 ! 1) CBMZ and its sub-modules may not be included in any commerical package,
26 ! or used for any commercial applications without the primary author's
27 ! prior consent.
28 ! 2) The CBMZ source code is provided to the WRF modeling community; however,
29 ! no portion of CBMZ can be used separately or in another code without the
30 ! primary author's prior consent.
31 ! 3) The CBMZ source code may be used for research, educational, and non-profit
32 ! purposes only. Any other usage must be first approved by the primary author.
33 ! 4) Publications resulting from the usage of CBMZ must use one or more of the
34 ! references below (depending on the application) for proper acknowledgment.
35 !
36 ! References:
37 ! 1) Zaveri R.A., and L.K. Peters (1999), A new lumped structure photochemical
38 ! mechanism for large-scale applications, J. Geophys. Res., 104, 30387-30415.
39 ! 2) Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C. Barnard, E.G.
40 ! Chapman, G.A. Grell, and S.E. Peckham (2005), Evolution of ozone, particulates,
41 ! and aerosol direct radiative forcing in the vicinity of Houston using a fully-
42 ! coupled meteorology-chemistry-aerosol model, Submitted to J. Geophys. Res.
43 !
44 ! Contact Jerome Fast for updates on the status of manuscripts under review.
45 !
46 ! Additional information:
47 ! 1) www.pnl.gov/atmos_sciences/raz
48 ! 2) www.pnl.gov/atmos_sciences/Jdf/wrfchem.html
49 !
50 ! Support:
51 ! Funding for developing and evaluating CBMZ was provided by the U.S. Department
52 ! of Energy under the auspices of Atmospheric Science Program of the Office of
53 ! Biological and Environmental Research the the PNNL Laboratory Research and
54 ! Directed Research and Development program.
55 !**********************************************************************************
56 module module_cbmz
57
58
59
60 use module_peg_util
61
62 contains
63
64
65 !***********************************************************************
66 ! < 1.> subr cbmz_driver
67 !
68 ! purpose: serves as an interface between subr. gas_chemistry and
69 ! the actual solver subr such as lsodes, rodas, etc.
70 !
71 ! grid : fixed i,j,k (box-model)
72 !
73 ! author : Rahul A. Zaveri
74 ! date : November 1998
75 !
76 !-----------------------------------------------------------------------
77
78 subroutine cbmz_driver( &
79 id, ktau, dtstep, ktauc, dtstepc, config_flags, &
80 gmt, julday, t_phy, moist, p8w, t8w, &
81 p_phy, chem, rho_phy, dz8w, z, z_at_w, vdrog3, &
82 ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
83 ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
84 ph_ch3o2h, ph_n2o5, &
85 ids,ide, jds,jde, kds,kde, &
86 ims,ime, jms,jme, kms,kme, &
87 its,ite, jts,jte, kts,kte )
88
89 USE module_configure, only: grid_config_rec_type, num_moist, num_chem, &
90 p_qv, p_so2, p_ho2, p_so4aj, p_corn, p_hcl, p_mtf
91 USE module_data_sorgam, only: ldrog
92 USE module_data_cbmz
93 IMPLICIT NONE
94
95
96 !-----------------------------------------------------------------------
97 ! subr arguments
98
99 INTEGER, INTENT(IN ) :: id, julday, &
100 ids,ide, jds,jde, kds,kde, &
101 ims,ime, jms,jme, kms,kme, &
102 its,ite, jts,jte, kts,kte
103
104 INTEGER, INTENT(IN ) :: ktau, ktauc
105
106 REAL, INTENT(IN ) :: dtstep, dtstepc, gmt
107 !
108 ! advected moisture variables
109 !
110 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
111 INTENT(IN ) :: moist
112 !
113 ! advected chemical tracers
114 !
115 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
116 INTENT(INOUT ) :: chem
117 !
118 ! arrays that hold photolysis rates
119 !
120 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
121 INTENT(INOUT ) :: &
122 ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
123 ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
124 ph_ch3o2h, ph_n2o5
125 !
126 ! on input from met model
127 !
128 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
129 INTENT(IN ) :: &
130 t_phy, & ! temperature
131 rho_phy, & ! air density (kg/m3)
132 p_phy, & ! NOT USED
133 z, z_at_w, & ! NOT USED
134 dz8w, & ! NOT USED
135 t8w, p8w ! NOT USED
136 !
137 ! for interaction with aerosols (really is output)
138 !
139 REAL, DIMENSION( ims:ime , kms:kme-0 , jms:jme , ldrog ) , &
140 INTENT(INOUT ) :: &
141 vdrog3 ! NOT USED
142
143 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
144
145
146 !-----------------------------------------------------------------------
147
148
149 ! local variables
150 integer :: idum, iok
151 integer :: iregime
152 integer :: igaschem_allowed_m1, igaschem_allowed_m2, &
153 igaschem_allowed_m3, igaschem_allowed_m4
154 integer :: igas_solver, iregime_forced
155 integer :: i_boxtest_units_convert
156 integer :: i_print_gasode_stats
157 integer :: i_force_dump, mode_force_dump
158 integer :: it, jt, kt
159 integer :: jsolver
160 integer :: lunerr, lunout, levdbg_err, levdbg_info
161 integer :: mgaschem
162
163
164 real :: abs_error, rel_error, trun
165 real :: tchem, dtchem
166 real :: tstart, tstop
167 real :: airdenbox, pressbox, tempbox
168 real :: cair_mlc
169 real :: h2o, ch4, oxygen, nitrogen, hydrogen
170 real :: cboxnew(ngas_z), cboxold(ngas_z)
171 real :: Aperox(nperox,nperox), Bperox(nperox,nperox)
172 real :: rk_param(nperox), rk_photo(nphoto)
173 real :: rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
174
175 integer, dimension(2,6), save :: inforodas=0
176 integer, dimension(6), save :: iodestatus_count=0, ioderegime_count=0
177
178 #ifdef CHEM_DBG_I
179 !rcetestb diagnostics --------------------------------------------------
180 print 93010, ' '
181 print 93010, 'rcetestb diagnostics from cbmz_driver'
182 print 93010, 'id, chem_opt, ktau, ktauc, julday ', &
183 id, config_flags%chem_opt, ktau, ktauc, julday
184 print 93020, 'dtstep, dtstepc, gmt ', &
185 dtstep, dtstepc, gmt
186 print 93010, 'ids/e, j, k', ids, ide, jds, jde, kds, kde
187 print 93010, 'ims/e, j, k', ims, ime, jms, jme, kms, kme
188 print 93010, 'its/e, j, k', its, ite, jts, jte, kts, kte
189 print 93010, 'num_moist, p_qv ', num_moist, p_qv
190 print 93010, 'num_chem, p_so2, p_ho2 ', num_chem, p_so2, p_ho2
191 print 93010, 'p_so4aj, p_corn, p_hcl, p_mtf', p_so4aj, p_corn, p_hcl, p_mtf
192 93010 format( a, 8(1x,i6) )
193 93020 format( a, 8(1p,e14.6) )
194 !rcetestb diagnostics --------------------------------------------------
195 #endif
196
197
198 ! set some control variables to their "standard for wrf-chem" values
199 igas_solver = 1
200 iregime_forced = -1
201 mgaschem = +1
202 i_boxtest_units_convert = 0
203
204 i_print_gasode_stats = 1
205 mode_force_dump = 0
206 lunerr = -1
207 lunout = -1
208 levdbg_err = 0
209 levdbg_info = 15
210
211 abs_error = 1.0e1 ! solver absolute tolerance (molecules/cm3)
212 rel_error = 1.0e-3 ! solver relative tolerance
213
214 ! set some control variables to non-standard values for testing
215 ! force dumps for center column, every 3rd level
216 ! mode_force_dump = +77
217 ! force dumps for center column, 1st level
218 ! mode_force_dump = +7
219 ! do levdbg_info output always
220 levdbg_info = 0
221
222
223 ! following call is for boxwrf testing only
224 ! it must be commented out for actual wrf applications
225 ! call boxtest_get_extra_args( &
226 ! igas_solver, iregime_forced, &
227 ! i_boxtest_units_convert, lunerr, lunout, &
228 ! abs_error, rel_error, trun )
229
230
231 ! currently nothing is done with vdrog3
232 ! vdrog3(its:ite,kts:kte,jts:jte,:) = 0.0 !This is already set to zero in chem_driver.
233
234 ! initialize rate constant arrays to 0.0 since not all elements are set
235 ! to valid values before some were used in loops, wig 16-Nov-2007
236 rk_m1(:) = 0.
237 rk_m2(:) = 0.
238 rk_m3(:) = 0.
239 rk_m4(:) = 0.
240
241 ! determine which regimes are allowed
242 ! based on which gas species are "active"
243 call set_gaschem_allowed_regimes( lunerr, &
244 igaschem_allowed_m1, igaschem_allowed_m2, &
245 igaschem_allowed_m3, igaschem_allowed_m4 )
246
247 !
248 ! main loop -- do gas chemistry at each i,j,k
249 !
250 do 2900 jt = jts, jte
251 do 2900 kt = kts, kte
252 do 2900 it = its, ite
253
254 trun = dtstep*(ktau-1) ! run time in s
255 tchem = gmt*3600.0 + dtstep*(ktau-1)
256 tchem = mod( tchem, 86400.0 ) ! time from 00 UTC in s
257 dtchem = dtstepc
258 tstart = tchem ! s
259 tstop = tstart + dtchem ! s
260
261 ! skip integration for very small dtchem
262 if ((tstop-tstart) .le. 1.0e-5) goto 2900
263
264
265 ! initial species mapping from host array
266 call mapgas_tofrom_host( 0, &
267 i_boxtest_units_convert, &
268 it,jt,kt, ims,ime, jms,jme, kms,kme, &
269 num_moist, num_chem, moist, chem, &
270 t_phy, p_phy, rho_phy, &
271 cboxold, tempbox, pressbox, airdenbox, &
272 cair_mlc, &
273 h2o, ch4, oxygen, nitrogen, hydrogen )
274 cboxnew(:) = cboxold(:)
275
276 ! determine regime
277 call selectgasregime( iregime, iregime_forced, cboxold, &
278 igaschem_allowed_m1, igaschem_allowed_m2, &
279 igaschem_allowed_m3, igaschem_allowed_m4 )
280 idum = iregime
281 if ((idum .lt. 1) .or. (idum .ge. 6)) idum = 6
282 ioderegime_count(idum) = ioderegime_count(idum) + 1
283 iodestatus_count(6) = iodestatus_count(6) + 1
284
285 ! compute rate constants
286 ! transfer/map incoming photolysis rate contants to local array
287 call gasphotoconstants( rk_photo, &
288 i_boxtest_units_convert, &
289 it,jt,kt, ims,ime, jms,jme, kms,kme, &
290 ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
291 ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
292 ph_ch3o2h, ph_n2o5 )
293 ! loads Aperox and Bperox
294 call loadperoxyparameters( Aperox, Bperox )
295 ! calculate parameterized rate constants
296 call peroxyrateconstants( tempbox, cboxold, &
297 Aperox, Bperox, rk_param )
298 ! calculate thermal rate constants
299 call gasrateconstants( iregime, tempbox, cair_mlc, &
300 rk_photo, rk_param, rk_m1, rk_m2, rk_m3, rk_m4 )
301
302 ! mode_force_dump selects a detailed dump of gaschem at either
303 ! first ijk grid, first ij column, all ijk, or no ijk
304 i_force_dump = 0
305 if (mode_force_dump .eq. 1) then
306 if ((it.eq.its) .and. (jt.eq.jts) &
307 .and. (kt.eq.kts)) i_force_dump = 1
308 else if (mode_force_dump .eq. 10) then
309 if ((it.eq.its) .and. (jt.eq.jts)) i_force_dump = 1
310 else if (mode_force_dump .eq. 100) then
311 i_force_dump = 1
312 else if (mode_force_dump .eq. 7) then
313 if ( (it .eq. (its+ite)/2) .and. &
314 (jt .eq. (jts+jte)/2) .and. &
315 (kt .eq. kts) ) i_force_dump = 1
316 else if (mode_force_dump .eq. 77) then
317 if ( (it .eq. (its+ite)/2) .and. &
318 (jt .eq. (jts+jte)/2) .and. &
319 (mod(kt-kts,3) .eq. 0) ) i_force_dump = 1
320 end if
321
322
323 ! rodas
324 iok = 0
325 jsolver = 0
326 if (igas_solver .eq. 1) then
327 jsolver = 1
328 call gasodesolver_rodas( tstart, tstop, iok, &
329 it, jt, kt, iregime, &
330 mgaschem, lunerr, lunout, levdbg_err, levdbg_info, &
331 i_force_dump, inforodas, iodestatus_count, &
332 abs_error, rel_error, trun, &
333 tempbox, pressbox, airdenbox, cboxnew, cboxold, &
334 rk_m1, rk_m2, rk_m3, rk_m4 )
335 endif
336
337 ! lsodes
338 if (igas_solver.eq.2 .or. iok.le.0) then
339 jsolver = 2
340 call gasodesolver_lsodes( tstart, tstop, iok, &
341 it, jt, kt, iregime, &
342 mgaschem, lunerr, lunout, levdbg_err, levdbg_info, &
343 i_force_dump, iodestatus_count, &
344 abs_error, rel_error, trun, &
345 tempbox, pressbox, airdenbox, cboxnew, cboxold, &
346 rk_m1, rk_m2, rk_m3, rk_m4 )
347 endif
348
349 ! final species mapping back to host array -- only when iok > 0
350 if (iok .gt. 0) then
351 call mapgas_tofrom_host( 1, &
352 i_boxtest_units_convert, &
353 it,jt,kt, ims,ime, jms,jme, kms,kme, &
354 num_moist, num_chem, moist, chem, &
355 t_phy, p_phy, rho_phy, &
356 cboxnew, tempbox, pressbox, airdenbox, &
357 cair_mlc, &
358 h2o, ch4, oxygen, nitrogen, hydrogen )
359 end if
360
361 ! following call is for boxwrf testing only
362 ! it must be commented out for actual wrf applications
363 ! call boxtest_set_extra_args( iregime, it, jt, kt )
364
365 2900 continue
366
367 if (i_print_gasode_stats .gt. 0) &
368 call print_gasode_stats( lunout, levdbg_info, &
369 inforodas, iodestatus_count, ioderegime_count )
370 return
371 end subroutine cbmz_driver
372
373
374
375 !***********************************************************************
376 ! < xx.> subr print_gasode_stats
377 !
378 ! purpose: writes some statistics on ode solver performance to unit lunout
379 !
380 !-----------------------------------------------------------------------
381
382 subroutine print_gasode_stats( lunout, levdbg, &
383 inforodas, iodestatus_count, ioderegime_count )
384
385 implicit none
386
387 ! subr arguments
388 integer lunout, levdbg
389 integer inforodas(2,6), iodestatus_count(6), ioderegime_count(6)
390
391 ! local variables
392 integer i, j
393 character*80 msg
394
395
396 msg = ' '
397 call peg_debugmsg( lunout, levdbg, msg )
398 msg = 'output from dump_cbmz_gasodeinfo'
399 call peg_debugmsg( lunout, levdbg, msg )
400 write(msg,9100) 'oderegime(1-6)', (ioderegime_count(i), i=1,6)
401 call peg_debugmsg( lunout, levdbg, msg )
402 write(msg,9100) 'odestatus(1-6)', (iodestatus_count(i), i=1,6)
403 call peg_debugmsg( lunout, levdbg, msg )
404
405 write(msg,9200) &
406 'inforodas(1-3)', ((inforodas(j,i), j=1,2), i=1,3)
407 call peg_debugmsg( lunout, levdbg, msg )
408 write(msg,9200) &
409 'inforodas(4-6)', ((inforodas(j,i), j=1,2), i=4,6)
410 call peg_debugmsg( lunout, levdbg, msg )
411
412 9100 format( a, 6i11 )
413 9200 format( a, 3( i11, '--', i9.9 ) )
414
415 return
416 end subroutine print_gasode_stats
417
418
419
420 !***********************************************************************
421 ! < xx.> subr gasodesolver_rodas
422 !
423 ! purpose: interfaces to rodas ode solver
424 !
425 !-----------------------------------------------------------------------
426
427 subroutine gasodesolver_rodas( tstart, tstop, iok, &
428 isvode, jsvode, ksvode, iregime, &
429 mgaschem, lunerr, lunout, levdbg_err, levdbg_info, &
430 i_force_dump, inforodas, iodestatus_count, &
431 abs_error, rel_error, trun, &
432 tempbox, pressbox, airdenbox, cboxnew, cboxold, &
433 rk_m1, rk_m2, rk_m3, rk_m4 )
434
435 use module_data_cbmz
436 use module_cbmz_rodas_prep, only: &
437 cbmz_v02r01_mapconcs, cbmz_v02r01_maprates, cbmz_v02r01_torodas, &
438 cbmz_v02r02_mapconcs, cbmz_v02r02_maprates, cbmz_v02r02_torodas, &
439 cbmz_v02r03_mapconcs, cbmz_v02r03_maprates, cbmz_v02r03_torodas, &
440 cbmz_v02r04_mapconcs, cbmz_v02r04_maprates, cbmz_v02r04_torodas, &
441 cbmz_v02r05_mapconcs, cbmz_v02r05_maprates, cbmz_v02r05_torodas, &
442 cbmz_v02r06_mapconcs, cbmz_v02r06_maprates, cbmz_v02r06_torodas
443
444 implicit none
445
446 ! subr arguments
447 integer iok, isvode, jsvode, ksvode, i_force_dump, iregime, &
448 mgaschem, lunerr, lunout, levdbg_err, levdbg_info
449 integer inforodas(2,6), iodestatus_count(6)
450 real tstart, tstop, abs_error, rel_error, trun
451 real tempbox, pressbox, airdenbox
452 real cboxnew(ngas_z), cboxold(ngas_z)
453 real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
454
455 ! local variables
456 integer :: ia, idum, idydt_sngldble, ig, l, ntot
457 integer, save :: nrodas_failures = 0
458 integer, dimension(6) :: inforodas_cur
459
460 real hmin, hstart, taa, tzz
461 real atolvec(ngas_z), rtolvec(ngas_z), &
462 stot(ngas_z), &
463 yposlimit(ngas_z), yneglimit(ngas_z)
464 real sfixedkpp(nfixed_kppmax), rconstkpp(nreact_kppmax)
465
466 character*80 msg
467
468 ! map reaction rate constants (pegasus --> kpp)
469 ! map concentrations (cboxold --> stot)
470 ! dump rates (for debugging)
471 if (iregime .eq. 1) then
472 call cbmz_v02r01_maprates( rk_m1, rk_m2, rk_m3, rk_m4, &
473 rconstkpp )
474 call cbmz_v02r01_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
475
476 else if (iregime .eq. 2) then
477 call cbmz_v02r02_maprates( rk_m1, rk_m2, rk_m3, rk_m4, &
478 rconstkpp )
479 call cbmz_v02r02_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
480
481 else if (iregime .eq. 3) then
482 call cbmz_v02r03_maprates( rk_m1, rk_m2, rk_m3, rk_m4, &
483 rconstkpp )
484 call cbmz_v02r03_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
485
486 else if (iregime .eq. 4) then
487 call cbmz_v02r04_maprates( rk_m1, rk_m2, rk_m3, rk_m4, &
488 rconstkpp )
489 call cbmz_v02r04_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
490
491 else if (iregime .eq. 5) then
492 call cbmz_v02r05_maprates( rk_m1, rk_m2, rk_m3, rk_m4, &
493 rconstkpp )
494 call cbmz_v02r05_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
495
496 else
497 call cbmz_v02r06_maprates( rk_m1, rk_m2, rk_m3, rk_m4, &
498 rconstkpp )
499 call cbmz_v02r06_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
500 end if
501
502 ! set parameters for rodas call
503 do l = 1, ntot
504 atolvec(l) = abs_error
505 rtolvec(l) = rel_error
506 yposlimit(l) = 1.0e20
507 yneglimit(l) = -1.0e8
508 end do
509
510 taa = tstart
511 tzz = tstop
512 hmin = 1.0e-5
513 hstart = 60.0
514 idydt_sngldble = 1
515
516 ! call rodas integrator
517 ! subr cbmz_v02r06_torodas(
518 ! + ngas, taa, tzz,
519 ! + stot, atol, rtol, yposlimit, yneglimit,
520 ! + hmin, hstart,
521 ! + inforodas_cur, iok, lunerr, idydt_sngldble )
522
523 if (iregime .eq. 1) then
524 call cbmz_v02r01_torodas( &
525 ntot, taa, tzz, &
526 stot, atolvec, rtolvec, yposlimit, yneglimit, &
527 sfixedkpp, rconstkpp, &
528 hmin, hstart, &
529 inforodas_cur, iok, lunerr, idydt_sngldble )
530
531 else if (iregime .eq. 2) then
532 call cbmz_v02r02_torodas( &
533 ntot, taa, tzz, &
534 stot, atolvec, rtolvec, yposlimit, yneglimit, &
535 sfixedkpp, rconstkpp, &
536 hmin, hstart, &
537 inforodas_cur, iok, lunerr, idydt_sngldble )
538
539 else if (iregime .eq. 3) then
540 call cbmz_v02r03_torodas( &
541 ntot, taa, tzz, &
542 stot, atolvec, rtolvec, yposlimit, yneglimit, &
543 sfixedkpp, rconstkpp, &
544 hmin, hstart, &
545 inforodas_cur, iok, lunerr, idydt_sngldble )
546
547 else if (iregime .eq. 4) then
548 call cbmz_v02r04_torodas( &
549 ntot, taa, tzz, &
550 stot, atolvec, rtolvec, yposlimit, yneglimit, &
551 sfixedkpp, rconstkpp, &
552 hmin, hstart, &
553 inforodas_cur, iok, lunerr, idydt_sngldble )
554
555 else if (iregime .eq. 5) then
556 call cbmz_v02r05_torodas( &
557 ntot, taa, tzz, &
558 stot, atolvec, rtolvec, yposlimit, yneglimit, &
559 sfixedkpp, rconstkpp, &
560 hmin, hstart, &
561 inforodas_cur, iok, lunerr, idydt_sngldble )
562
563 else
564 call cbmz_v02r06_torodas( &
565 ntot, taa, tzz, &
566 stot, atolvec, rtolvec, yposlimit, yneglimit, &
567 sfixedkpp, rconstkpp, &
568 hmin, hstart, &
569 inforodas_cur, iok, lunerr, idydt_sngldble )
570 end if
571
572
573 ! increment odeinfo counters
574 if (iok .gt. 0) then
575 if (inforodas_cur(6) .le. 0) then
576 ia = 1
577 else
578 ia = 2
579 end if
580 else
581 ia = 3
582 end if
583 iodestatus_count(ia) = iodestatus_count(ia) + 1
584 ! do following to avoid overflow of the "inforodas" numbers
585 ! inforodas(2,i) contains rightmost 9 digits of each inforodas number
586 ! inforodas(1,i) contains any higher digits of each inforodas number
587 do ia = 1, 6
588 idum = inforodas(2,ia) + inforodas_cur(ia)
589 inforodas(1,ia) = inforodas(1,ia) + (idum/1000000000)
590 inforodas(2,ia) = mod(idum, 1000000000)
591 end do
592
593
594 ! map concentrations (stot --> cboxnew)
595 if (iregime .eq. 1) then
596 call cbmz_v02r01_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
597 else if (iregime .eq. 2) then
598 call cbmz_v02r02_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
599 else if (iregime .eq. 3) then
600 call cbmz_v02r03_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
601 else if (iregime .eq. 4) then
602 call cbmz_v02r04_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
603 else if (iregime .eq. 5) then
604 call cbmz_v02r05_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
605 else
606 call cbmz_v02r06_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
607 end if
608
609
610 ! diagnostic output if integration fails OR if i_force_dump > 0
611 if (iok .gt. 0) then
612 if (i_force_dump .le. 0) goto 20000
613 else
614 nrodas_failures = nrodas_failures + 1
615 if (nrodas_failures .gt. 100) goto 20000
616 end if
617
618 msg = ' '
619 call peg_debugmsg( lunout, levdbg_err, msg )
620 if (iok .gt. 0) then
621 msg = '*** gasodesolver_rodas forced dump'
622 else
623 write(msg,*) '*** gasodesolver_rodas failure no.', &
624 nrodas_failures
625 end if
626 call peg_debugmsg( lunout, levdbg_err, msg )
627 msg = 'iregime, iok, i, j, k / t'
628 call peg_debugmsg( lunout, levdbg_err, msg )
629 write(msg,97010) iregime, iok, isvode, jsvode, ksvode
630 call peg_debugmsg( lunout, levdbg_err, msg )
631 write(msg,97020) trun
632 call peg_debugmsg( lunout, levdbg_err, msg )
633 msg = 'inforodas_cur(1-6) ='
634 call peg_debugmsg( lunout, levdbg_err, msg )
635 write(msg,97010) inforodas_cur
636 call peg_debugmsg( lunout, levdbg_err, msg )
637 msg = &
638 'tstart, tstop, abs_error, rel_error / temp, press, cair, cos_sza ='
639 call peg_debugmsg( lunout, levdbg_err, msg )
640 write(msg,97020) tstart, tstop, abs_error, rel_error
641 call peg_debugmsg( lunout, levdbg_err, msg )
642 write(msg,97020) tempbox, pressbox, airdenbox, -99.0
643 call peg_debugmsg( lunout, levdbg_err, msg )
644
645 idum = 0
646 do ig = nreact_kppmax, 1, -1
647 if ((idum .eq. 0) .and. (rconstkpp(ig) .ne. 0.0)) idum = ig
648 end do
649 msg = 'ngas_z, nrconst_nonzero ='
650 call peg_debugmsg( lunout, levdbg_err, msg )
651 write(msg,97010) ngas_z, idum
652 call peg_debugmsg( lunout, levdbg_err, msg )
653 msg = 'l, name, cboxold, cboxnew for l=1,ngas_z'
654 call peg_debugmsg( lunout, levdbg_err, msg )
655 do l = 1, ngas_z
656 write(msg,97030) l, name_z(l), cboxold(l), cboxnew(l)
657 call peg_debugmsg( lunout, levdbg_err, msg )
658 end do
659 msg = 'rconst for i=1,nrconst_nonzero'
660 call peg_debugmsg( lunout, levdbg_err, msg )
661 do ia = 1, idum, 4
662 write(msg,97020) ( rconstkpp(ig), ig = ia, min(ia+3,idum) )
663 call peg_debugmsg( lunout, levdbg_err, msg )
664 end do
665
666 97010 format( 6i12 )
667 97020 format( 4(1pe18.10) )
668 97030 format(( i3, 1x, a, 2(1pe18.10) ))
669
670
671 ! force non-negative values
672 20000 do l = 1, ngas_z
673 cboxnew(l) = max( cboxnew(l), 0.0 )
674 end do
675
676 return
677 end subroutine gasodesolver_rodas
678
679
680
681 !***********************************************************************
682 ! < xx.> subr gasodesolver_lsodes
683 !
684 ! purpose: interface to lsodes ode solver
685 !
686 ! author : Rahul A. Zaveri
687 ! date : May, 2000
688 !
689 !-----------------------------------------------------------------------
690
691 subroutine gasodesolver_lsodes( tstart, tstop, iok, &
692 isvode, jsvode, ksvode, iregime, &
693 mgaschem, lunerr, lunout, levdbg_err, levdbg_info, &
694 i_force_dump, iodestatus_count, &
695 abs_error, rel_error, trun, &
696 tempbox, pressbox, airdenbox, cboxnew, cboxold, &
697 rk_m1, rk_m2, rk_m3, rk_m4 )
698
699 use module_data_cbmz
700 use module_cbmz_lsodes_solver, only: lsodes_solver, xsetf, &
701 set_lsodes_common_vars
702 implicit none
703
704 ! subr arguments
705 integer i, iok, isvode, jsvode, ksvode, i_force_dump, iregime, &
706 mgaschem, lunerr, lunout, levdbg_err, levdbg_info
707 integer iodestatus_count(6)
708 real tstart, tstop, abs_error, rel_error, trun
709 real tempbox, pressbox, airdenbox
710 real cboxnew(ngas_z), cboxold(ngas_z)
711 real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
712
713 ! lsodes parameters and local variables
714 integer itoler, itask, iopt, mf, lwm, nrdim, nidim
715 integer nruserpar, niuserpar
716 parameter( itoler = 1, itask = 1, iopt = 1, mf= 222 )
717 parameter( lwm = 3*ngas_tot*ngas_tot + 12*ngas_tot )
718 parameter( nrdim = 20 + 9*ngas_tot + lwm )
719 parameter( nidim = 31 + ngas_tot + ngas_tot*ngas_tot )
720 parameter( nruserpar = 5 + nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4)
721 parameter( niuserpar = ngas_z + 1 )
722
723 integer ia, idum, ig, ioffset, istate, iwork(nidim), l
724 integer ntotvec(1), iuserpar(niuserpar)
725 integer indx(ngas_z)
726 integer, save :: iflagout = 0
727 integer, save :: nlsodes_failures = 0
728
729 real dtchem, rwork(nrdim), stot(ngas_tot)
730 real atolvec(1), rtolvec(1), ruserpar(nruserpar)
731
732 character*80 msg
733
734
735
736 iok = 1 ! reset
737
738 call set_lsodes_common_vars()
739
740 ! sets gas species indices for iregime
741 call setgasindices( iregime, indx )
742
743 ! map cboxold --> stot
744 call mapgasspecies( cboxold, stot, 0, iregime, indx )
745
746 !----------------------------------------------------------------------
747 ! set number of species (ntot) for the selected regime for LSODES
748 if (iregime .eq. 1) then
749 ntotvec(1) = ngas_r1
750 else if (iregime .eq. 2) then
751 ntotvec(1) = ngas_r2
752 else if (iregime .eq. 3) then
753 ntotvec(1) = ngas_r3
754 else if (iregime .eq. 4) then
755 ntotvec(1) = ngas_r4
756 else if (iregime .eq. 5) then
757 ntotvec(1) = ngas_r5
758 else
759 ntotvec(1) = ngas_r6
760 end if
761
762 100 continue
763
764 ! set other LSODES parameters...
765 iwork(6) = 1000 ! max iterations for a time step
766 iwork(7) = 1
767 istate = 1
768 rwork(6) = dtchem
769 if(iflagout.eq.0)then
770 call xsetf(iflagout)
771 endif
772
773 atolvec(1) = abs_error
774 rtolvec(1) = rel_error
775
776 do ig = 1, 5
777 ruserpar(ig) = ig*7.0
778 end do
779 ruserpar(1) = cboxold(ih2o_z)
780 ruserpar(2) = cboxold(ich4_z)
781 ruserpar(3) = cboxold(io2_z)
782 ruserpar(4) = cboxold(in2_z)
783 ruserpar(5) = cboxold(ih2_z)
784 ioffset = 5
785 do ig = 1, nrxn_m1
786 ruserpar(ioffset+ig) = rk_m1(ig)
787 end do
788 ioffset = ioffset + nrxn_m1
789 do ig = 1, nrxn_m2
790 ruserpar(ioffset+ig) = rk_m2(ig)
791 end do
792 ioffset = ioffset + nrxn_m2
793 do ig = 1, nrxn_m3
794 ruserpar(ioffset+ig) = rk_m3(ig)
795 end do
796 ioffset = ioffset + nrxn_m3
797 do ig = 1, nrxn_m4
798 ruserpar(ioffset+ig) = rk_m4(ig)
799 end do
800
801 iuserpar(1) = iregime
802 do ig = 1, ngas_z
803 iuserpar(1+ig) = indx(ig)
804 end do
805
806 call lsodes_solver( &
807 gasode_cbmz, ntotvec, stot, tstart, tstop, &
808 itoler, rtolvec, atolvec, itask, istate, iopt, &
809 rwork, nrdim, iwork, nidim, jcs, mf, &
810 ruserpar, nruserpar, iuserpar, niuserpar )
811
812 if (istate .le. 0) iok = -1
813
814
815 ! increment odeinfo counters
816 if (iok .gt. 0) then
817 ia = 4
818 else
819 ia = 5
820 end if
821 iodestatus_count(ia) = iodestatus_count(ia) + 1
822
823
824 ! map stot --> cboxnew
825 call mapgasspecies( cboxnew, stot, 1, iregime, indx )
826
827
828 ! do diagnostic output if integration fails OR if i_force_dump > 0
829 if (iok .gt. 0) then
830 if (i_force_dump .le. 0) goto 20000
831 else
832 nlsodes_failures = nlsodes_failures + 1
833 end if
834
835 msg = ' '
836 call peg_debugmsg( lunout, levdbg_err, msg )
837 if (iok .gt. 0) then
838 msg = '*** gasodesolver_lsodes forced dump'
839 else
840 write(msg,*) '*** gasodesolver_lsodes failure no.', &
841 nlsodes_failures
842 end if
843 call peg_debugmsg( lunout, levdbg_err, msg )
844 msg = 'iregime, iok, i, j, k / t'
845 call peg_debugmsg( lunout, levdbg_err, msg )
846 write(msg,97010) iregime, iok, isvode, jsvode, ksvode
847 call peg_debugmsg( lunout, levdbg_err, msg )
848 write(msg,97020) trun
849 call peg_debugmsg( lunout, levdbg_err, msg )
850 if (nlsodes_failures .gt. 1000) then
851 write(msg,*) '*** exceeded lsodes failure limit =', 1000
852 call peg_debugmsg( lunout, levdbg_err, msg )
853 call peg_error_fatal( lunerr, msg )
854 end if
855 if (nlsodes_failures .gt. 100) goto 20000
856
857 write(msg,*) 'istate -', istate
858 call peg_debugmsg( lunout, levdbg_err, msg )
859 msg = &
860 'tstart, tstop, abs_error, rel_error / temp, press, cair, cos_sza ='
861 call peg_debugmsg( lunout, levdbg_err, msg )
862 write(msg,97020) tstart, tstop, abs_error, rel_error
863 call peg_debugmsg( lunout, levdbg_err, msg )
864 write(msg,97020) tempbox, pressbox, airdenbox, -99.0
865 call peg_debugmsg( lunout, levdbg_err, msg )
866
867 idum = nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4
868 msg = 'ngas_z, nrconst_m1+m2+m3+m4 ='
869 call peg_debugmsg( lunout, levdbg_err, msg )
870 write(msg,97010) ngas_z, idum
871 call peg_debugmsg( lunout, levdbg_err, msg )
872 msg = 'l, name, cboxold, cboxnew for l=1,ngas_z'
873 call peg_debugmsg( lunout, levdbg_err, msg )
874 do l = 1, ngas_z
875 write(msg,97030) l, name_z(l), cboxold(l), cboxnew(l)
876 call peg_debugmsg( lunout, levdbg_err, msg )
877 end do
878 msg = 'rconst for i=1,nrconst_nonzero'
879 call peg_debugmsg( lunout, levdbg_err, msg )
880 do ia = 1, nrxn_m1, 4
881 write(msg,97020) ( rk_m1(ig), ig = ia, min(ia+3,nrxn_m1) )
882 call peg_debugmsg( lunout, levdbg_err, msg )
883 end do
884 do ia = 1, nrxn_m2, 4
885 write(msg,97020) ( rk_m2(ig), ig = ia, min(ia+3,nrxn_m2) )
886 call peg_debugmsg( lunout, levdbg_err, msg )
887 end do
888 do ia = 1, nrxn_m3, 4
889 write(msg,97020) ( rk_m3(ig), ig = ia, min(ia+3,nrxn_m3) )
890 call peg_debugmsg( lunout, levdbg_err, msg )
891 end do
892 do ia = 1, nrxn_m4, 4
893 write(msg,97020) ( rk_m4(ig), ig = ia, min(ia+3,nrxn_m4) )
894 call peg_debugmsg( lunout, levdbg_err, msg )
895 end do
896
897 97010 format( 6i12 )
898 97020 format( 4(1pe18.10) )
899 97030 format(( i3, 1x, a, 2(1pe18.10) ))
900
901
902 ! force non-negative values
903 20000 do l = 1, ngas_z
904 cboxnew(l) = max( cboxnew(l), 0.0 )
905 end do
906
907 return
908 end subroutine gasodesolver_lsodes
909
910
911
912 !***********************************************************************
913 ! < 2.> subr selectgasregime
914 !
915 ! purpose: selects an optimum combination of gas-phase
916 ! mechanisms based on sensitivity of [OH]
917 ! concentrations to some lumped structure
918 ! hydrocarbon groups concentrations and [DMS]
919 ! concentration.
920 !
921 ! input : cbox = full species concentrations array (mol/cc)
922 !
923 ! output: iregime = 1 : com
924 ! = 2 : com + urb
925 ! = 3 : com + urb + bio
926 ! = 4 : com + mar
927 ! = 5 : com + urb + mar
928 ! = 6 : com + urb + bio + mar
929 ! ntot = number of gas-phase species in the selected mechanism
930 !
931 ! author: Rahul A. Zaveri
932 ! date : April 2000
933 !
934 !---------------------------------------------------------------------
935
936 subroutine selectgasregime( iregime, iregime_forced, cbox, &
937 igaschem_allowed_m1, igaschem_allowed_m2, &
938 igaschem_allowed_m3, igaschem_allowed_m4 )
939
940 use module_data_cbmz
941 implicit none
942
943 ! subr arguments
944 integer iregime, iregime_forced
945 integer igaschem_allowed_m1, igaschem_allowed_m2, &
946 igaschem_allowed_m3, igaschem_allowed_m4
947 real cbox(ngas_z)
948
949 ! local variables
950 integer iwork(6)
951 integer m_m1, m_m2, m_m3, m_m4
952 real cut_molecpcc
953
954
955 cut_molecpcc = 5.e+6 ! molecules/cc
956
957 ! initialize mechanism flags
958 m_m1 = 1 ! 1 (always)
959 m_m2 = 0 ! 0 or 1
960 m_m3 = 0 ! 0 or 2
961 m_m4 = 0 ! 0 or 3
962
963 if (igaschem_allowed_m2 .gt. 0) then
964 if (cbox(ipar_z ) .gt. cut_molecpcc) m_m2 = 1
965 if (cbox(iaone_z ) .gt. cut_molecpcc) m_m2 = 1
966 if (cbox(imgly_z ) .gt. cut_molecpcc) m_m2 = 1
967 if (cbox(ieth_z ) .gt. cut_molecpcc) m_m2 = 1
968 if (cbox(iolet_z ) .gt. cut_molecpcc) m_m2 = 1
969 if (cbox(iolei_z ) .gt. cut_molecpcc) m_m2 = 1
970 if (cbox(ixyl_z ) .gt. cut_molecpcc) m_m2 = 1
971 if (cbox(icres_z ) .gt. cut_molecpcc) m_m2 = 1
972 if (cbox(ito2_z ) .gt. cut_molecpcc) m_m2 = 1
973 if (cbox(icro_z ) .gt. cut_molecpcc) m_m2 = 1
974 if (cbox(iopen_z ) .gt. cut_molecpcc) m_m2 = 1
975 if (cbox(ionit_z ) .gt. cut_molecpcc) m_m2 = 1
976 if (cbox(irooh_z ) .gt. cut_molecpcc) m_m2 = 1
977 if (cbox(iro2_z ) .gt. cut_molecpcc) m_m2 = 1
978 if (cbox(iano2_z ) .gt. cut_molecpcc) m_m2 = 1
979 if (cbox(inap_z ) .gt. cut_molecpcc) m_m2 = 1
980 if (cbox(ixo2_z ) .gt. cut_molecpcc) m_m2 = 1
981 if (cbox(ixpar_z ) .gt. cut_molecpcc) m_m2 = 1
982 end if
983
984 if (igaschem_allowed_m3 .gt. 0) then
985 if (cbox(iisop_z ) .gt. cut_molecpcc) m_m3 = 2
986 if (cbox(iisoprd_z ) .gt. cut_molecpcc) m_m3 = 2
987 if (cbox(iisopp_z ) .gt. cut_molecpcc) m_m3 = 2
988 if (cbox(iisopn_z ) .gt. cut_molecpcc) m_m3 = 2
989 if (cbox(iisopo2_z ) .gt. cut_molecpcc) m_m3 = 2
990 end if
991
992 if (igaschem_allowed_m4 .gt. 0) then
993 if (cbox(idms_z ) .gt. cut_molecpcc) m_m4 = 3
994 if (cbox(imsa_z ) .gt. cut_molecpcc) m_m4 = 3
995 if (cbox(idmso_z ) .gt. cut_molecpcc) m_m4 = 3
996 if (cbox(idmso2_z ) .gt. cut_molecpcc) m_m4 = 3
997 if (cbox(ich3so2h_z ) .gt. cut_molecpcc) m_m4 = 3
998 if (cbox(ich3sch2oo_z ) .gt. cut_molecpcc) m_m4 = 3
999 if (cbox(ich3so2_z ) .gt. cut_molecpcc) m_m4 = 3
1000 if (cbox(ich3so3_z ) .gt. cut_molecpcc) m_m4 = 3
1001 if (cbox(ich3so2oo_z ) .gt. cut_molecpcc) m_m4 = 3
1002 if (cbox(ich3so2ch2oo_z) .gt. cut_molecpcc) m_m4 = 3
1003 if (cbox(imtf_z ) .gt. cut_molecpcc) m_m4 = 3
1004 end if
1005
1006 iregime = m_m1 + m_m2*((2-m_m3)/2) + m_m3 + m_m4
1007
1008 ! force iregime = iregime_forced
1009 if ((iregime_forced .ge. 1) .and. (iregime_forced .le. 6)) &
1010 iregime = iregime_forced
1011
1012 return
1013 end subroutine selectgasregime
1014
1015
1016
1017 !***********************************************************************
1018 ! < 3.> subr setgasindices
1019 !
1020 ! purpose: sets gas species indices
1021 !
1022 ! author : Rahul A. Zaveri
1023 ! date : May, 2000
1024 !
1025 !-----------------------------------------------------------------------
1026
1027 subroutine setgasindices( iregime, indx )
1028
1029 use module_data_cbmz
1030 implicit none
1031
1032 ! subr arguments
1033 integer iregime, indx(ngas_z)
1034
1035 ! local variables
1036 integer ilast
1037
1038 ilast = 0
1039 indx(:) = -999888777
1040
1041 goto (1,2,3,4,5,6), iregime
1042
1043 1 call setgasindex_m1( ilast, indx ) ! regime 1
1044 return
1045
1046 2 call setgasindex_m1( ilast, indx ) ! regime 2
1047 call setgasindex_m2( ilast, indx )
1048 return
1049
1050 3 call setgasindex_m1( ilast, indx ) ! regime 3
1051 call setgasindex_m2( ilast, indx )
1052 call setgasindex_m3( ilast, indx )
1053 return
1054
1055 4 call setgasindex_m1( ilast, indx ) ! regime 4
1056 call setgasindex_m4( ilast, indx )
1057 return
1058
1059 5 call setgasindex_m1( ilast, indx ) ! regime 5
1060 call setgasindex_m2( ilast, indx )
1061 call setgasindex_m4( ilast, indx )
1062 return
1063
1064 6 call setgasindex_m1( ilast, indx ) ! regime 6
1065 call setgasindex_m2( ilast, indx )
1066 call setgasindex_m3( ilast, indx )
1067 call setgasindex_m4( ilast, indx )
1068 return
1069
1070 end subroutine setgasindices
1071
1072
1073
1074
1075 !***********************************************************************
1076 ! < 4.> subr gasrateconstants
1077 !
1078 ! purpose: calls regime-dependent subrs for calculating
1079 ! gas-phase thermal reaction rate constants
1080 !
1081 ! author : Rahul A. Zaveri
1082 ! date : May, 2000
1083 !
1084 !-----------------------------------------------------------------------
1085
1086 subroutine gasrateconstants( iregime, tempbox, cair_mlc, &
1087 rk_photo, rk_param, rk_m1, rk_m2, rk_m3, rk_m4 )
1088
1089 use module_data_cbmz
1090 implicit none
1091
1092 ! subr arguments
1093 integer iregime
1094 real tempbox, cair_mlc
1095 real rk_photo(nphoto), rk_param(nperox)
1096 real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
1097
1098
1099 ! iregime=1 --> do m1
1100 ! iregime=2 --> do m1, m2
1101 ! iregime=3 --> do m1, m2, m3
1102 ! iregime=4 --> do m1, --, --, m4
1103 ! iregime=5 --> do m1, m2, --, m4
1104 ! iregime=6 --> do m1, m2, m3, m4
1105
1106 call gasthermrk_m1( tempbox, cair_mlc, &
1107 rk_photo, rk_param, rk_m1, rk_m2 )
1108
1109 if ((iregime .eq. 2) .or. &
1110 (iregime .eq. 3) .or. &
1111 (iregime .eq. 5) .or. &
1112 (iregime .eq. 6)) &
1113 call gasthermrk_m2( tempbox, cair_mlc, &
1114 rk_photo, rk_param, rk_m2 )
1115
1116 if ((iregime .eq. 3) .or. &
1117 (iregime .eq. 6)) &
1118 call gasthermrk_m3( tempbox, cair_mlc, &
1119 rk_photo, rk_param, rk_m3 )
1120
1121 if ((iregime .eq. 4) .or. &
1122 (iregime .eq. 5) .or. &
1123 (iregime .eq. 6)) &
1124 call gasthermrk_m4( tempbox, cair_mlc, &
1125 rk_photo, rk_param, rk_m4 )
1126
1127 return
1128 end subroutine gasrateconstants
1129
1130
1131
1132 !***********************************************************************
1133 ! < 3.> subr setgasindex_m1
1134 !
1135 ! purpose: defines gas species indices for regime 1
1136 !
1137 ! author : Rahul A. Zaveri
1138 ! date : December, 1998
1139 !
1140 !-----------------------------------------------------------------------
1141
1142 subroutine setgasindex_m1( ilast, indx )
1143
1144 use module_data_cbmz
1145 implicit none
1146
1147 ! subr arguments
1148 integer ilast, indx(ngas_z)
1149
1150
1151 indx(ino_z) = 1
1152 indx(ino2_z) = 2
1153 indx(ino3_z) = 3
1154 indx(in2o5_z) = 4
1155 indx(ihono_z) = 5
1156 indx(ihno3_z) = 6
1157 indx(ihno4_z) = 7
1158 indx(io3_z) = 8
1159 indx(io1d_z) = 9
1160 indx(io3p_z) = 10
1161 indx(ioh_z) = 11
1162 indx(iho2_z) = 12
1163 indx(ih2o2_z) = 13
1164 indx(ico_z) = 14
1165 indx(iso2_z) = 15
1166 indx(ih2so4_z) = 16
1167 indx(inh3_z) = 17
1168 indx(ihcl_z) = 18
1169 indx(ic2h6_z) = 19
1170 indx(ich3o2_z) = 20
1171 indx(iethp_z) = 21
1172 indx(ihcho_z) = 22
1173 indx(ich3oh_z) = 23
1174 indx(ic2h5oh_z) = 24
1175 indx(ich3ooh_z) = 25
1176 indx(iethooh_z) = 26
1177 indx(iald2_z) = 27
1178 indx(ihcooh_z) = 28
1179 indx(ircooh_z) = 29
1180 indx(ic2o3_z) = 30
1181 indx(ipan_z) = 31
1182
1183 ilast = indx(ipan_z)
1184
1185 return
1186 end subroutine setgasindex_m1
1187
1188
1189
1190 !***********************************************************************
1191 ! < 4.> subr setgasindex_m2
1192 !
1193 ! purpose: defines gas species indices for regime 2
1194 !
1195 ! author : Rahul A. Zaveri
1196 ! date : December, 1998
1197 !
1198 !-----------------------------------------------------------------------
1199 subroutine setgasindex_m2( ilast, indx )
1200
1201 use module_data_cbmz
1202 implicit none
1203
1204 ! subr arguments
1205 integer ilast, indx(ngas_z)
1206
1207
1208 indx(ipar_z) = ilast + 1
1209 indx(iaone_z) = ilast + 2
1210 indx(imgly_z) = ilast + 3
1211 indx(ieth_z) = ilast + 4
1212 indx(iolet_z) = ilast + 5
1213 indx(iolei_z) = ilast + 6
1214 indx(itol_z) = ilast + 7
1215 indx(ixyl_z) = ilast + 8
1216 indx(icres_z) = ilast + 9
1217 indx(ito2_z) = ilast + 10
1218 indx(icro_z) = ilast + 11
1219 indx(iopen_z) = ilast + 12
1220 indx(ionit_z) = ilast + 13
1221 ! indx(ipan_z) = ilast + 14
1222 ! indx(ircooh_z) = ilast + 15
1223 indx(irooh_z) = ilast + 14
1224 ! indx(ic2o3_z) = ilast + 17
1225 indx(iro2_z) = ilast + 15
1226 indx(iano2_z) = ilast + 16
1227 indx(inap_z) = ilast + 17
1228 indx(ixo2_z) = ilast + 18
1229 indx(ixpar_z) = ilast + 19
1230
1231 ilast = indx(ixpar_z)
1232
1233 return
1234 end subroutine setgasindex_m2
1235
1236
1237
1238 !***********************************************************************
1239 ! < 5.> subr setgasindex_m3
1240 !
1241 ! purpose: defines gas species indices for regime 3
1242 !
1243 ! author : Rahul A. Zaveri
1244 ! date : December, 1998
1245 !
1246 !-----------------------------------------------------------------------
1247 subroutine setgasindex_m3( ilast, indx )
1248
1249 use module_data_cbmz
1250 implicit none
1251
1252 ! subr arguments
1253 integer ilast, indx(ngas_z)
1254
1255
1256 indx(iisop_z) = ilast + 1
1257 indx(iisoprd_z) = ilast + 2
1258 indx(iisopp_z) = ilast + 3
1259 indx(iisopn_z) = ilast + 4
1260 indx(iisopo2_z) = ilast + 5
1261
1262 ilast = indx(iisopo2_z)
1263
1264 return
1265 end subroutine setgasindex_m3
1266
1267
1268
1269 !***********************************************************************
1270 ! < 6.> subr setgasindex_m4
1271 !
1272 ! purpose: defines gas species indices for regime 4
1273 !
1274 ! author : Rahul A. Zaveri
1275 ! date : December, 1998
1276 !
1277 !-----------------------------------------------------------------------
1278
1279 subroutine setgasindex_m4( ilast, indx )
1280
1281 use module_data_cbmz
1282 implicit none
1283
1284 ! subr arguments
1285 integer ilast, indx(ngas_z)
1286
1287
1288 !
1289 indx(idms_z) = ilast + 1
1290 indx(imsa_z) = ilast + 2
1291 indx(idmso_z) = ilast + 3
1292 indx(idmso2_z) = ilast + 4
1293 indx(ich3so2h_z) = ilast + 5
1294 indx(ich3sch2oo_z) = ilast + 6
1295 indx(ich3so2_z) = ilast + 7
1296 indx(ich3so3_z) = ilast + 8
1297 indx(ich3so2oo_z) = ilast + 9
1298 indx(ich3so2ch2oo_z) = ilast + 10
1299 indx(imtf_z) = ilast + 11
1300
1301 ilast = indx(imtf_z)
1302
1303 return
1304 end subroutine setgasindex_m4
1305
1306
1307
1308 !***********************************************************************
1309 ! < 9.> subr mapgasspecies
1310 !
1311 ! purpose: map gas species between stot and cbox arrays
1312 !
1313 ! author : Rahul A. Zaveri
1314 ! date : December, 1998
1315 !
1316 ! ----------------------------------------------------------------------
1317
1318 subroutine mapgasspecies( cbox, stot, imap, iregime, indx )
1319
1320 use module_data_cbmz
1321 implicit none
1322
1323 ! subr arguments
1324 integer imap, iregime, indx(ngas_z)
1325 real cbox(ngas_z)
1326 real stot(ngas_tot)
1327 !
1328 !
1329 goto (1,2,3,4,5,6), iregime
1330 !
1331 !
1332 1 call mapgas_m1( cbox, stot, imap, indx )
1333 return
1334 !
1335 !
1336 2 call mapgas_m1( cbox, stot, imap, indx )
1337 call mapgas_m2( cbox, stot, imap, indx )
1338 return
1339 !
1340 !
1341 3 call mapgas_m1( cbox, stot, imap, indx )
1342 call mapgas_m2( cbox, stot, imap, indx )
1343 call mapgas_m3( cbox, stot, imap, indx )
1344 return
1345 !
1346 !
1347 4 call mapgas_m1( cbox, stot, imap, indx )
1348 call mapgas_m4( cbox, stot, imap, indx )
1349 return
1350 !
1351 !
1352 5 call mapgas_m1( cbox, stot, imap, indx )
1353 call mapgas_m2( cbox, stot, imap, indx )
1354 call mapgas_m4( cbox, stot, imap, indx )
1355 return
1356 !
1357 !
1358 6 call mapgas_m1( cbox, stot, imap, indx )
1359 call mapgas_m2( cbox, stot, imap, indx )
1360 call mapgas_m3( cbox, stot, imap, indx )
1361 call mapgas_m4( cbox, stot, imap, indx )
1362 return
1363
1364 end subroutine mapgasspecies
1365
1366
1367
1368 !***********************************************************************
1369 ! <10.> subr mapgas_m1
1370 !
1371 ! purpose: maps gas species between stot and cbox arrays for mechanism 1
1372 !
1373 ! author : Rahul A. Zaveri
1374 ! date : December, 1998
1375 !
1376 ! ----------------------------------------------------------------------
1377
1378 subroutine mapgas_m1( cbox, stot, imap, indx )
1379
1380 use module_data_cbmz
1381 implicit none
1382
1383 ! subr arguments
1384 integer imap, indx(ngas_z)
1385 real cbox(ngas_z)
1386 real stot(ngas_tot)
1387
1388 if(imap.eq.0)then ! map cbox --> stot (both molec/cc)
1389 stot(indx(ino_z)) = cbox(ino_z)
1390 stot(indx(ino2_z)) = cbox(ino2_z)
1391 stot(indx(ino3_z)) = cbox(ino3_z)
1392 stot(indx(in2o5_z)) = cbox(in2o5_z)
1393 stot(indx(ihono_z)) = cbox(ihono_z)
1394 stot(indx(ihno3_z)) = cbox(ihno3_z)
1395 stot(indx(ihno4_z)) = cbox(ihno4_z)
1396 stot(indx(io3_z)) = cbox(io3_z)
1397 stot(indx(io1d_z)) = cbox(io1d_z)
1398 stot(indx(io3p_z)) = cbox(io3p_z)
1399 stot(indx(ioh_z)) = cbox(ioh_z)
1400 stot(indx(iho2_z)) = cbox(iho2_z)
1401 stot(indx(ih2o2_z)) = cbox(ih2o2_z)
1402 stot(indx(ico_z)) = cbox(ico_z)
1403 stot(indx(iso2_z)) = cbox(iso2_z)
1404 stot(indx(ih2so4_z)) = cbox(ih2so4_z)
1405 stot(indx(inh3_z)) = cbox(inh3_z)
1406 stot(indx(ihcl_z)) = cbox(ihcl_z)
1407 stot(indx(ic2h6_z)) = cbox(ic2h6_z)
1408 stot(indx(ich3o2_z)) = cbox(ich3o2_z)
1409 stot(indx(iethp_z)) = cbox(iethp_z)
1410 stot(indx(ihcho_z)) = cbox(ihcho_z)
1411 stot(indx(ich3oh_z)) = cbox(ich3oh_z)
1412 stot(indx(ic2h5oh_z)) = cbox(ic2h5oh_z)
1413 stot(indx(ich3ooh_z)) = cbox(ich3ooh_z)
1414 stot(indx(iethooh_z)) = cbox(iethooh_z)
1415 stot(indx(iald2_z)) = cbox(iald2_z)
1416 stot(indx(ihcooh_z)) = cbox(ihcooh_z)
1417 stot(indx(ircooh_z)) = cbox(ircooh_z)
1418 stot(indx(ic2o3_z)) = cbox(ic2o3_z)
1419 stot(indx(ipan_z)) = cbox(ipan_z)
1420 !
1421 else ! map stot --> cbox (both molec/cc)
1422 cbox(ino_z) = stot(indx(ino_z))
1423 cbox(ino2_z) = stot(indx(ino2_z))
1424 cbox(ino3_z) = stot(indx(ino3_z))
1425 cbox(in2o5_z) = stot(indx(in2o5_z))
1426 cbox(ihono_z) = stot(indx(ihono_z))
1427 cbox(ihno3_z) = stot(indx(ihno3_z))
1428 cbox(ihno4_z) = stot(indx(ihno4_z))
1429 cbox(io3_z) = stot(indx(io3_z))
1430 cbox(io1d_z) = stot(indx(io1d_z))
1431 cbox(io3p_z) = stot(indx(io3p_z))
1432 cbox(ioh_z) = stot(indx(ioh_z))
1433 cbox(iho2_z) = stot(indx(iho2_z))
1434 cbox(ih2o2_z) = stot(indx(ih2o2_z))
1435 cbox(ico_z) = stot(indx(ico_z))
1436 cbox(iso2_z) = stot(indx(iso2_z))
1437 cbox(ih2so4_z) = stot(indx(ih2so4_z))
1438 cbox(inh3_z) = stot(indx(inh3_z))
1439 cbox(ihcl_z) = stot(indx(ihcl_z))
1440 cbox(ic2h6_z) = stot(indx(ic2h6_z))
1441 cbox(ich3o2_z) = stot(indx(ich3o2_z))
1442 cbox(iethp_z) = stot(indx(iethp_z))
1443 cbox(ihcho_z) = stot(indx(ihcho_z))
1444 cbox(ich3oh_z) = stot(indx(ich3oh_z))
1445 cbox(ic2h5oh_z) = stot(indx(ic2h5oh_z))
1446 cbox(ich3ooh_z) = stot(indx(ich3ooh_z))
1447 cbox(iethooh_z) = stot(indx(iethooh_z))
1448 cbox(iald2_z) = stot(indx(iald2_z))
1449 cbox(ihcooh_z) = stot(indx(ihcooh_z))
1450 cbox(ircooh_z) = stot(indx(ircooh_z))
1451 cbox(ic2o3_z) = stot(indx(ic2o3_z))
1452 cbox(ipan_z) = stot(indx(ipan_z))
1453 endif
1454
1455 return
1456 end subroutine mapgas_m1
1457
1458
1459
1460 !***********************************************************************
1461 ! <11.> subr mapgas_m2
1462 !
1463 ! purpose: maps gas species between stot and cbox arrays for mechanism 2
1464 !
1465 ! author : Rahul A. Zaveri
1466 ! date : December, 1998
1467 !
1468 ! ----------------------------------------------------------------------
1469
1470 subroutine mapgas_m2( cbox, stot, imap, indx )
1471
1472 use module_data_cbmz
1473 implicit none
1474
1475 ! subr arguments
1476 integer imap, indx(ngas_z)
1477 real cbox(ngas_z)
1478 real stot(ngas_tot)
1479
1480 if(imap.eq.0)then ! map cbox --> stot (both molec/cc)
1481 stot(indx(ipar_z)) = cbox(ipar_z)
1482 stot(indx(iaone_z)) = cbox(iaone_z)
1483 stot(indx(imgly_z)) = cbox(imgly_z)
1484 stot(indx(ieth_z)) = cbox(ieth_z)
1485 stot(indx(iolet_z)) = cbox(iolet_z)
1486 stot(indx(iolei_z)) = cbox(iolei_z)
1487 stot(indx(itol_z)) = cbox(itol_z)
1488 stot(indx(ixyl_z)) = cbox(ixyl_z)
1489 stot(indx(icres_z)) = cbox(icres_z)
1490 stot(indx(ito2_z)) = cbox(ito2_z)
1491 stot(indx(icro_z)) = cbox(icro_z)
1492 stot(indx(iopen_z)) = cbox(iopen_z)
1493 stot(indx(ionit_z)) = cbox(ionit_z)
1494 ! stot(indx(ipan_z)) = cbox(ipan_z)
1495 ! stot(indx(ircooh_z)) = cbox(ircooh_z)
1496 stot(indx(irooh_z)) = cbox(irooh_z)
1497 ! stot(indx(ic2o3_z)) = cbox(ic2o3_z)
1498 stot(indx(iro2_z)) = cbox(iro2_z)
1499 stot(indx(iano2_z)) = cbox(iano2_z)
1500 stot(indx(inap_z)) = cbox(inap_z)
1501 stot(indx(ixo2_z)) = cbox(ixo2_z)
1502 stot(indx(ixpar_z)) = cbox(ixpar_z)
1503 !
1504 else ! map stot --> cbox (both molec/cc)
1505 cbox(ipar_z) = stot(indx(ipar_z))
1506 cbox(iaone_z) = stot(indx(iaone_z))
1507 cbox(imgly_z) = stot(indx(imgly_z))
1508 cbox(ieth_z) = stot(indx(ieth_z))
1509 cbox(iolet_z) = stot(indx(iolet_z))
1510 cbox(iolei_z) = stot(indx(iolei_z))
1511 cbox(itol_z) = stot(indx(itol_z))
1512 cbox(ixyl_z) = stot(indx(ixyl_z))
1513 cbox(icres_z) = stot(indx(icres_z))
1514 cbox(ito2_z) = stot(indx(ito2_z))
1515 cbox(icro_z) = stot(indx(icro_z))
1516 cbox(iopen_z) = stot(indx(iopen_z))
1517 cbox(ionit_z) = stot(indx(ionit_z))
1518 ! cbox(ipan_z) = stot(indx(ipan_z))
1519 ! cbox(ircooh_z) = stot(indx(ircooh_z))
1520 cbox(irooh_z) = stot(indx(irooh_z))
1521 ! cbox(ic2o3_z) = stot(indx(ic2o3_z))
1522 cbox(iro2_z) = stot(indx(iro2_z))
1523 cbox(iano2_z) = stot(indx(iano2_z))
1524 cbox(inap_z) = stot(indx(inap_z))
1525 cbox(ixo2_z) = stot(indx(ixo2_z))
1526 cbox(ixpar_z) = stot(indx(ixpar_z))
1527 endif
1528
1529 return
1530 end subroutine mapgas_m2
1531
1532
1533
1534 !***********************************************************************
1535 ! <12.> subr mapgas_m3
1536 !
1537 ! purpose: maps gas species between stot and cbox arrays for mechanism 3
1538 !
1539 ! author : Rahul A. Zaveri
1540 ! date : December, 1998
1541 !
1542 ! ----------------------------------------------------------------------
1543
1544 subroutine mapgas_m3( cbox, stot, imap, indx )
1545
1546 use module_data_cbmz
1547 implicit none
1548
1549 ! subr arguments
1550 integer imap, indx(ngas_z)
1551 real cbox(ngas_z)
1552 real stot(ngas_tot)
1553
1554 if(imap.eq.0)then ! map cbox --> stot (both molec/cc)
1555 stot(indx(iisop_z)) = cbox(iisop_z)
1556 stot(indx(iisoprd_z)) = cbox(iisoprd_z)
1557 stot(indx(iisopp_z)) = cbox(iisopp_z)
1558 stot(indx(iisopn_z)) = cbox(iisopn_z)
1559 stot(indx(iisopo2_z)) = cbox(iisopo2_z)
1560 !
1561 else ! map stot --> cbox (both molec/cc)
1562 cbox(iisop_z) = stot(indx(iisop_z))
1563 cbox(iisoprd_z) = stot(indx(iisoprd_z))
1564 cbox(iisopp_z) = stot(indx(iisopp_z))
1565 cbox(iisopn_z) = stot(indx(iisopn_z))
1566 cbox(iisopo2_z) = stot(indx(iisopo2_z))
1567 endif
1568
1569 return
1570 end subroutine mapgas_m3
1571
1572
1573
1574 !***********************************************************************
1575 ! <13.> subr mapgas_m4
1576 !
1577 ! purpose: maps gas species between stot and cbox arrays for mechanism 4
1578 !
1579 ! author : Rahul A. Zaveri
1580 ! date : December, 1998
1581 !
1582 ! ----------------------------------------------------------------------
1583
1584 subroutine mapgas_m4( cbox, stot, imap, indx )
1585
1586 use module_data_cbmz
1587 implicit none
1588
1589 ! subr arguments
1590 integer imap, indx(ngas_z)
1591 real cbox(ngas_z)
1592 real stot(ngas_tot)
1593
1594 if(imap.eq.0)then ! map cbox --> stot (both molec/cc)
1595 stot(indx(idms_z)) = cbox(idms_z)
1596 stot(indx(imsa_z)) = cbox(imsa_z)
1597 stot(indx(idmso_z)) = cbox(idmso_z)
1598 stot(indx(idmso2_z)) = cbox(idmso2_z)
1599 stot(indx(ich3so2h_z)) = cbox(ich3so2h_z)
1600 stot(indx(ich3sch2oo_z)) = cbox(ich3sch2oo_z)
1601 stot(indx(ich3so2_z)) = cbox(ich3so2_z)
1602 stot(indx(ich3so3_z)) = cbox(ich3so3_z)
1603 stot(indx(ich3so2oo_z)) = cbox(ich3so2oo_z)
1604 stot(indx(ich3so2ch2oo_z))= cbox(ich3so2ch2oo_z)
1605 stot(indx(imtf_z)) = cbox(imtf_z)
1606 !
1607 else ! map stot --> cbox (both molec/cc)
1608 cbox(idms_z) = stot(indx(idms_z))
1609 cbox(imsa_z) = stot(indx(imsa_z))
1610 cbox(idmso_z) = stot(indx(idmso_z))
1611 cbox(idmso2_z) = stot(indx(idmso2_z))
1612 cbox(ich3so2h_z) = stot(indx(ich3so2h_z))
1613 cbox(ich3sch2oo_z) = stot(indx(ich3sch2oo_z))
1614 cbox(ich3so2_z) = stot(indx(ich3so2_z))
1615 cbox(ich3so3_z) = stot(indx(ich3so3_z))
1616 cbox(ich3so2oo_z) = stot(indx(ich3so2oo_z))
1617 cbox(ich3so2ch2oo_z)= stot(indx(ich3so2ch2oo_z))
1618 cbox(imtf_z) = stot(indx(imtf_z))
1619 endif
1620
1621 return
1622 end subroutine mapgas_m4
1623
1624
1625
1626 !***********************************************************************
1627 ! <xx.> subr check_userpar
1628 !
1629 ! purpose: called by lsodes (external)
1630 ! computes time derivatives of species concentrations ds/dt.
1631 ! by calling subr. gasrate and ode
1632 !
1633 ! author : Rahul A. Zaveri
1634 ! date : December 1998
1635 !
1636 !---------------------------------------------------------------------------
1637
1638 subroutine check_userpar( ruserpar, nruserpar, iuserpar, niuserpar )
1639
1640 use module_data_cbmz
1641 implicit none
1642
1643 ! subr arguments
1644 integer nruserpar, niuserpar
1645 integer iuserpar(niuserpar)
1646 real ruserpar(nruserpar)
1647
1648 ! local variables
1649 integer i
1650 real dum
1651 character*80 msg
1652
1653 if (nruserpar .ne. (5 + nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4)) then
1654 write(msg,9010) 'nruserpar', -1, nruserpar
1655 call wrf_error_fatal( msg )
1656 end if
1657
1658 if (niuserpar .ne. (ngas_z + 1)) then
1659 write(msg,9010) 'niuserpar', -1, niuserpar
1660 call wrf_error_fatal( msg )
1661 end if
1662
1663 9010 format( '*** check_userpar error -- ', a, 1x, i8, 1x, i8 )
1664
1665
1666 return
1667 end subroutine check_userpar
1668
1669
1670
1671 !***********************************************************************
1672 ! <14.> subr gasode_cbmz
1673 !
1674 ! purpose: called by lsodes (external)
1675 ! computes time derivatives of species concentrations ds/dt.
1676 ! by calling subr. gasrate and ode
1677 !
1678 ! author : Rahul A. Zaveri
1679 ! date : December 1998
1680 !
1681 !---------------------------------------------------------------------------
1682
1683 subroutine gasode_cbmz( ntot, tt, s, sdot, &
1684 ruserpar, nruserpar, iuserpar, niuserpar )
1685
1686 use module_data_cbmz
1687 implicit none
1688
1689 ! subr arguments
1690 integer ntot, nruserpar, niuserpar
1691 integer iuserpar(niuserpar)
1692 real tt
1693 real s(ngas_tot), sdot(ngas_tot), ruserpar(nruserpar)
1694
1695 ! local variables
1696 integer ig, ioffset, iregime, irxn
1697 integer indx(ngas_z)
1698 real h2o, ch4, oxygen, nitrogen, hydrogen
1699 real rk_m1(nrxn_m1), r_m1(nrxn_m1)
1700 real rk_m2(nrxn_m2), r_m2(nrxn_m2)
1701 real rk_m3(nrxn_m3), r_m3(nrxn_m3)
1702 real rk_m4(nrxn_m4), r_m4(nrxn_m4)
1703 real p_m1(ngas_tot), d_m1(ngas_tot)
1704 real p_m2(ngas_tot), d_m2(ngas_tot)
1705 real p_m3(ngas_tot), d_m3(ngas_tot)
1706 real p_m4(ngas_tot), d_m4(ngas_tot)
1707
1708
1709 ! test on userpar
1710 call check_userpar( ruserpar, nruserpar, iuserpar, niuserpar )
1711
1712 iregime = iuserpar(1)
1713 do ig = 1, ngas_z
1714 indx(ig) = iuserpar(ig+1)
1715 end do
1716
1717 h2o = ruserpar(1)
1718 ch4 = ruserpar(2)
1719 oxygen = ruserpar(3)
1720 nitrogen = ruserpar(4)
1721 hydrogen = ruserpar(5)
1722 ioffset = 5
1723 do ig = 1, nrxn_m1
1724 rk_m1(ig) = ruserpar(ioffset+ig)
1725 end do
1726 ioffset = ioffset + nrxn_m1
1727 do ig = 1, nrxn_m2
1728 rk_m2(ig) = ruserpar(ioffset+ig)
1729 end do
1730 ioffset = ioffset + nrxn_m2
1731 do ig = 1, nrxn_m3
1732 rk_m3(ig) = ruserpar(ioffset+ig)
1733 end do
1734 ioffset = ioffset + nrxn_m3
1735 do ig = 1, nrxn_m4
1736 rk_m4(ig) = ruserpar(ioffset+ig)
1737 end do
1738
1739
1740 ! initialize to zero
1741 do irxn=1,nrxn_m1
1742 r_m1(irxn) = 0.
1743 enddo
1744
1745 do irxn=1,nrxn_m2
1746 r_m2(irxn) = 0.
1747 enddo
1748
1749 do irxn=1,nrxn_m3
1750 r_m3(irxn) = 0.
1751 enddo
1752
1753 do irxn=1,nrxn_m4
1754 r_m4(irxn) = 0.
1755 enddo
1756 !
1757 !
1758 ! initialize to zero
1759 do ig=1,ngas_tot
1760 p_m1(ig) = 0.
1761 p_m2(ig) = 0.
1762 p_m3(ig) = 0.
1763 p_m4(ig) = 0.
1764 !
1765 d_m1(ig) = 0.
1766 d_m2(ig) = 0.
1767 d_m3(ig) = 0.
1768 d_m4(ig) = 0.
1769 enddo
1770
1771
1772 ! iregime=1 --> do m1
1773 ! iregime=2 --> do m1, m2
1774 ! iregime=3 --> do m1, m2, m3
1775 ! iregime=4 --> do m1, --, --, m4
1776 ! iregime=5 --> do m1, m2, --, m4
1777 ! iregime=6 --> do m1, m2, m3, m4
1778
1779 call gasrate_m1( indx, s, r_m1, r_m2, rk_m1, rk_m2, &
1780 h2o, ch4, oxygen, nitrogen, hydrogen )
1781
1782 if ((iregime .eq. 2) .or. &
1783 (iregime .eq. 3) .or. &
1784 (iregime .eq. 5) .or. &
1785 (iregime .eq. 6)) &
1786 call gasrate_m2( indx, s, r_m2, rk_m2 )
1787
1788 if ((iregime .eq. 3) .or. &
1789 (iregime .eq. 6)) &
1790 call gasrate_m3( indx, s, r_m3, rk_m3 )
1791
1792 if ((iregime .eq. 4) .or. &
1793 (iregime .eq. 5) .or. &
1794 (iregime .eq. 6)) &
1795 call gasrate_m4( indx, s, r_m4, rk_m4, oxygen )
1796
1797 call gasode_m1( indx, r_m1, p_m1, d_m1, r_m2 )
1798
1799 if ((iregime .eq. 2) .or. &
1800 (iregime .eq. 3) .or. &
1801 (iregime .eq. 5) .or. &
1802 (iregime .eq. 6)) &
1803 call gasode_m2( indx, r_m2, p_m2, d_m2 )
1804
1805 if ((iregime .eq. 3) .or. &
1806 (iregime .eq. 6)) &
1807 call gasode_m3( indx, r_m3, p_m3, d_m3 )
1808
1809 if ((iregime .eq. 4) .or. &
1810 (iregime .eq. 5) .or. &
1811 (iregime .eq. 6)) &
1812 call gasode_m4( indx, r_m4, p_m4, d_m4 )
1813
1814
1815 if (iregime .eq. 1) then
1816 ! regime = 1
1817 do ig = 1, ngas_r1
1818 sdot(ig) = real( dble(p_m1(ig)) - &
1819 dble(d_m1(ig)) )
1820 end do
1821
1822 else if (iregime .eq. 2) then
1823 ! regime = 2
1824 do ig = 1, ngas_r2
1825 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)) - &
1826 dble(d_m1(ig)+d_m2(ig)) )
1827 end do
1828
1829 else if (iregime .eq. 3) then
1830 ! regime = 3
1831 do ig = 1, ngas_r3
1832 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m3(ig)) - &
1833 dble(d_m1(ig)+d_m2(ig)+d_m3(ig)) )
1834 end do
1835
1836 else if (iregime .eq. 4) then
1837 ! regime = 4
1838 do ig = 1, ngas_r4
1839 sdot(ig) = real( dble(p_m1(ig)+p_m4(ig)) - &
1840 dble(d_m1(ig)+d_m4(ig)) )
1841 end do
1842
1843 else if (iregime .eq. 5) then
1844 ! regime = 5
1845 do ig = 1, ngas_r5
1846 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m4(ig)) - &
1847 dble(d_m1(ig)+d_m2(ig)+d_m4(ig)) )
1848 end do
1849
1850 else if (iregime .eq. 6) then
1851 ! regime = 6
1852 do ig = 1, ngas_r6
1853 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m3(ig)+p_m4(ig)) - &
1854 dble(d_m1(ig)+d_m2(ig)+d_m3(ig)+d_m4(ig)) )
1855 end do
1856
1857 end if
1858
1859 return
1860 end subroutine gasode_cbmz
1861
1862
1863
1864 !***********************************************************************
1865 ! <15.> subr jcs
1866 !
1867 ! purpose: external dummy jacobian evaluation for LSODES (when mf=222)
1868 !
1869 !-----------------------------------------------------------------------
1870
1871 subroutine jcs( ngas, tt, s, j, ian, jan, pdj, &
1872 ruserpar, nruserpar, iuserpar, niuserpar )
1873
1874 implicit none
1875 integer ngas, j, ian(*), jan(*), nruserpar, niuserpar
1876 integer iuserpar(niuserpar)
1877 real tt, s(*), pdj(*)
1878 real ruserpar(nruserpar)
1879
1880 ! test on userpar
1881 call check_userpar( ruserpar, nruserpar, iuserpar, niuserpar )
1882
1883 return
1884 end subroutine jcs
1885
1886
1887
1888 !***********************************************************************
1889 ! <16.> subr gasode_m1
1890 !
1891 ! purpose: updates production and destruction rates for mechanism 1
1892 ! background troposphere
1893 !
1894 ! author : Rahul A. Zaveri
1895 ! date : December, 1998
1896 !
1897 ! ----------------------------------------------------------------------
1898
1899 subroutine gasode_m1( indx, r_m1, p_m1, d_m1, r_m2 )
1900
1901 use module_data_cbmz
1902 implicit none
1903
1904 ! subr arguments
1905 integer indx(ngas_z)
1906 real r_m1(nrxn_m1), p_m1(ngas_tot), d_m1(ngas_tot)
1907 real r_m2(nrxn_m2)
1908
1909
1910 p_m1(indx(ino_z))= r_m1(1)+0.11*r_m1(2) &
1911 +r_m1(3)+r_m1(15)+r_m1(38)
1912 d_m1(indx(ino_z))= r_m1(17)+r_m1(18)+r_m1(23) &
1913 +r_m1(33)+r_m1(37)+r_m1(57) &
1914 +r_m1(58) &
1915 +r_m2(34)
1916 !
1917 p_m1(indx(ino2_z))= 0.89*r_m1(2)+r_m1(4) &
1918 +r_m1(5)+r_m1(6)+r_m1(17) &
1919 +r_m1(18)+r_m1(25) &
1920 +r_m1(26)+r_m1(28) &
1921 +r_m1(33)+r_m1(36) &
1922 +r_m1(37)+r_m1(37) &
1923 +r_m1(38)+r_m1(40) &
1924 +r_m1(40)+.7*r_m1(41) &
1925 +r_m1(43)+r_m1(57) &
1926 +r_m1(58)+r_m1(59) &
1927 +r_m1(60) &
1928 +r_m2(32)+r_m2(34)+r_m2(39)
1929 d_m1(indx(ino2_z))= r_m1(1)+r_m1(15)+r_m1(16) &
1930 +r_m1(19)+r_m1(24) &
1931 +r_m1(34)+r_m1(35) &
1932 +r_m1(38)+r_m1(39) &
1933 +r_m2(31)
1934 !
1935 p_m1(indx(ino3_z))= r_m1(6)+r_m1(16)+r_m1(19) &
1936 +r_m1(27)+r_m1(43)
1937 d_m1(indx(ino3_z))= r_m1(2)+r_m1(25)+r_m1(37) &
1938 +r_m1(38)+r_m1(39) &
1939 +r_m1(40)+r_m1(40) &
1940 +r_m1(41)+r_m1(52) &
1941 +r_m1(59)+r_m1(60) &
1942 +r_m2(4)+r_m2(39)
1943 !
1944 p_m1(indx(in2o5_z))= r_m1(39)
1945 d_m1(indx(in2o5_z))= r_m1(6)+r_m1(42) &
1946 +r_m1(43)
1947 !
1948 p_m1(indx(ihono_z))= r_m1(23)+r_m1(35)
1949 d_m1(indx(ihono_z))= r_m1(3)+r_m1(26)
1950 !
1951 p_m1(indx(ihno3_z))= r_m1(24)+.3*r_m1(41) &
1952 +r_m1(42)+r_m1(42) &
1953 +r_m1(52) &
1954 +r_m2(4)
1955 d_m1(indx(ihno3_z))= r_m1(4)+r_m1(27)
1956 !
1957 p_m1(indx(ihno4_z))= r_m1(34)
1958 d_m1(indx(ihno4_z))= r_m1(5)+r_m1(28) &
1959 +r_m1(36)
1960 !
1961 p_m1(indx(io3_z))= r_m1(13) &
1962 +.4*r_m2(44)
1963 d_m1(indx(io3_z))= r_m1(7)+r_m1(8)+r_m1(14) &
1964 +r_m1(18)+r_m1(19)+r_m1(20) &
1965 +r_m1(21)
1966 !
1967 p_m1(indx(io1d_z))= r_m1(8)
1968 d_m1(indx(io1d_z))= r_m1(10)+r_m1(11) &
1969 +r_m1(12)
1970 !
1971 p_m1(indx(io3p_z))= r_m1(1)+0.89*r_m1(2) &
1972 +r_m1(7)+r_m1(10)+r_m1(11)
1973 d_m1(indx(io3p_z))= r_m1(13)+r_m1(14) &
1974 +r_m1(15)+r_m1(16) &
1975 +r_m1(17)
1976 !
1977 p_m1(indx(ioh_z))= r_m1(3)+r_m1(4)+2*r_m1(9) &
1978 +2*r_m1(12)+r_m1(21) &
1979 +r_m1(33)+.7*r_m1(41) &
1980 +r_m1(53)+r_m1(54)+.3*r_m1(55) &
1981 +.5*r_m1(56)
1982 d_m1(indx(ioh_z))= r_m1(20)+r_m1(22)+r_m1(23) &
1983 +r_m1(24)+r_m1(25)+r_m1(26) &
1984 +r_m1(27)+r_m1(28)+r_m1(29) &
1985 +r_m1(30)+r_m1(44)+r_m1(45) &
1986 +r_m1(46)+r_m1(47)+r_m1(48) &
1987 +r_m1(51)+r_m1(55)+r_m1(56) &
1988 +r_m1(65) &
1989 +r_m2(3)
1990 !
1991 p_m1(indx(iho2_z))= r_m1(5)+r_m1(20)+r_m1(22) &
1992 +r_m1(25)+r_m1(30) &
1993 +r_m1(36)+r_m1(44) &
1994 +r_m1(45)+r_m1(48) &
1995 +2*r_m1(49)+r_m1(51) &
1996 +r_m1(52)+r_m1(53) &
1997 +r_m1(54)+r_m1(57) &
1998 +r_m1(58)+r_m1(59) &
1999 +r_m1(60)+.32*r_m1(63) &
2000 +.6*r_m1(64)+r_m1(65) &
2001 +r_m2(2)
2002 d_m1(indx(iho2_z))= r_m1(21)+r_m1(29) &
2003 +r_m1(31)+r_m1(31) &
2004 +r_m1(32)+r_m1(32) &
2005 +r_m1(33)+r_m1(34) &
2006 +r_m1(35)+r_m1(41) &
2007 +r_m1(61)+r_m1(62) &
2008 +r_m2(44)
2009 !
2010 p_m1(indx(ih2o2_z))= r_m1(31)+r_m1(32)
2011 d_m1(indx(ih2o2_z))= r_m1(9)+r_m1(30)
2012 !
2013 p_m1(indx(ico_z))= r_m1(49)+r_m1(50)+r_m1(51) &
2014 +r_m1(52) &
2015 +r_m2(2)
2016 d_m1(indx(ico_z))= r_m1(44)
2017 !
2018 p_m1(indx(iso2_z))= 0.0
2019 d_m1(indx(iso2_z))= r_m1(45)
2020 !
2021 p_m1(indx(ih2so4_z))= r_m1(45)
2022 d_m1(indx(ih2so4_z))= 0.0
2023 !
2024 p_m1(indx(inh3_z))= 0.0
2025 d_m1(indx(inh3_z))= 0.0
2026 !
2027 p_m1(indx(ihcl_z))= 0.0
2028 d_m1(indx(ihcl_z))= 0.0
2029 !
2030 p_m1(indx(ic2h6_z))= .2*r_m1(64)
2031 d_m1(indx(ic2h6_z))= r_m1(47)
2032 !
2033 p_m1(indx(ich3o2_z))= r_m1(46)+.7*r_m1(55) &
2034 +r_m2(2)+r_m2(34)+r_m2(39)+r_m2(49)
2035 d_m1(indx(ich3o2_z))= r_m1(57)+r_m1(59) &
2036 +r_m1(61)+r_m1(63)
2037 !
2038 p_m1(indx(iethp_z))= r_m1(47)+.5*r_m1(56)
2039 d_m1(indx(iethp_z))= r_m1(58)+r_m1(60) &
2040 +r_m1(62)+r_m1(64)
2041 !
2042 p_m1(indx(ihcho_z))= r_m1(48)+r_m1(53) &
2043 +.3*r_m1(55)+r_m1(57) &
2044 +r_m1(59)+.66*r_m1(63)
2045 d_m1(indx(ihcho_z))= r_m1(49)+r_m1(50) &
2046 +r_m1(51)+r_m1(52)
2047 !
2048 p_m1(indx(ich3oh_z))= .34*r_m1(63)
2049 d_m1(indx(ich3oh_z))= r_m1(48)
2050 !
2051 p_m1(indx(ic2h5oh_z))= 0.0
2052 d_m1(indx(ic2h5oh_z))= r_m1(65)
2053 !
2054 p_m1(indx(ich3ooh_z))= r_m1(61)
2055 d_m1(indx(ich3ooh_z))= r_m1(53)+r_m1(55)
2056 !
2057 p_m1(indx(iethooh_z))= r_m1(62)
2058 d_m1(indx(iethooh_z))= r_m1(54)+r_m1(56)
2059 !
2060 p_m1(indx(iald2_z))= r_m1(54)+.5*r_m1(56) &
2061 +r_m1(58)+r_m1(60) &
2062 +.8*r_m1(64)+r_m1(65)
2063 d_m1(indx(iald2_z))= r_m2(2)+r_m2(3)+r_m2(4)
2064 !
2065 p_m1(indx(ihcooh_z))= 0.0
2066 d_m1(indx(ihcooh_z))= 0.0
2067 !
2068 p_m1(indx(ircooh_z))= .4*r_m2(44)
2069 d_m1(indx(ircooh_z))= 0.0
2070 !
2071 p_m1(indx(ic2o3_z))= r_m2(3)+r_m2(4)+r_m2(32)
2072 d_m1(indx(ic2o3_z))= r_m2(31)+r_m2(34)+r_m2(39)+r_m2(44)+r_m2(49)
2073 !
2074 p_m1(indx(ipan_z))= r_m2(31)
2075 d_m1(indx(ipan_z))= r_m2(32)
2076
2077 return
2078 end subroutine gasode_m1
2079
2080
2081
2082 !***********************************************************************
2083 ! <17.> subr gasode_m2
2084 !
2085 ! purpose: updates production and destruction rates for mechanism 2
2086 ! anthropogenic hydrocarbons
2087 !
2088 ! author : Rahul A. Zaveri
2089 ! date : December, 1998
2090 !
2091 ! ----------------------------------------------------------------------
2092
2093 subroutine gasode_m2( indx, r_m2, p_m2, d_m2 )
2094
2095 use module_data_cbmz
2096 implicit none
2097
2098 ! subr arguments
2099 integer indx(ngas_z)
2100 real r_m2(nrxn_m2), p_m2(ngas_tot), d_m2(ngas_tot)
2101
2102
2103 p_m2(indx(ino_z))= 0.0
2104 d_m2(indx(ino_z))= r_m2(20)+r_m2(33) +r_m2(35) &
2105 +r_m2(36)+r_m2(37)
2106 !
2107 p_m2(indx(ino2_z))= .95*r_m2(20)+r_m2(30) +.84*r_m2(33) &
2108 +r_m2(35)+1.5*r_m2(36)+r_m2(37) &
2109 +r_m2(38) +r_m2(40)+1.5*r_m2(41)+r_m2(42) &
2110 +.5*r_m2(51)
2111 d_m2(indx(ino2_z))= r_m2(23)
2112 !
2113 p_m2(indx(ino3_z))= 0.0
2114 d_m2(indx(ino3_z))= +r_m2(9)+r_m2(16)+r_m2(17)+r_m2(22) &
2115 +r_m2(38) +r_m2(40)+r_m2(41)+r_m2(42)
2116 !
2117 p_m2(indx(ihno3_z))= +r_m2(9)+r_m2(22)
2118 d_m2(indx(ihno3_z))= 0.0
2119 !
2120 p_m2(indx(io3_z))= 0.0
2121 d_m2(indx(io3_z))= r_m2(10)+r_m2(12)+r_m2(13)+r_m2(26)
2122 !
2123 p_m2(indx(ioh_z))= .12*r_m2(10)+.33*r_m2(12)+.60*r_m2(13) &
2124 +.08*r_m2(26)+r_m2(27)+.23*r_m2(28)
2125 d_m2(indx(ioh_z))= r_m2(1) +r_m2(6)+r_m2(8)+r_m2(11) &
2126 +r_m2(14)+r_m2(15)+r_m2(18)+r_m2(19)+r_m2(21) &
2127 +r_m2(24)+r_m2(28)+r_m2(29)
2128 !
2129 p_m2(indx(iho2_z))= +r_m2(7)+.22*r_m2(10)+r_m2(11) &
2130 +.26*r_m2(12)+.22*r_m2(13)+r_m2(14)+r_m2(15) &
2131 +.2*r_m2(18)+.55*r_m2(19)+.95*r_m2(20) &
2132 +.6*r_m2(21)+2*r_m2(24)+r_m2(25)+.76*r_m2(26) &
2133 +.9*r_m2(27)+.9*r_m2(30)+.76*r_m2(33)+.5*r_m2(36) &
2134 +.9*r_m2(38)+.5*r_m2(41)+.54*r_m2(48)
2135 d_m2(indx(iho2_z))= r_m2(43) +r_m2(45)+r_m2(46)+r_m2(47)
2136 !
2137 p_m2(indx(ico_z))= +r_m2(7)+r_m2(9)+.24*r_m2(10) &
2138 +.31*r_m2(12)+.30*r_m2(13)+2*r_m2(24)+r_m2(25) &
2139 +.69*r_m2(26)
2140 d_m2(indx(ico_z))= 0.0
2141 !
2142 p_m2(indx(ipar_z))= 1.1*r_m2(19)
2143 d_m2(indx(ipar_z))= r_m2(1) + r_m2(53)
2144 !
2145 p_m2(indx(ich3oh_z))= .03*r_m2(12)+.04*r_m2(13)
2146 d_m2(indx(ich3oh_z))= 0.0
2147 !
2148 p_m2(indx(ihcho_z))= r_m2(10)+1.56*r_m2(11)+.57*r_m2(12)+r_m2(14) &
2149 +r_m2(24)+.7*r_m2(26)+r_m2(35)+.5*r_m2(36) &
2150 +r_m2(40)+.5*r_m2(41)+.7*r_m2(50)+.5*r_m2(51)
2151 d_m2(indx(ihcho_z))= 0.0
2152 !
2153 p_m2(indx(iald2_z))= .22*r_m2(11)+.47*r_m2(12)+1.03*r_m2(13) &
2154 +r_m2(14)+1.77*r_m2(15)+.03*r_m2(26)+.3*r_m2(27) &
2155 +.04*r_m2(28)+.3*r_m2(30)+.25*r_m2(33)+.5*r_m2(36) &
2156 +.3*r_m2(38)+.5*r_m2(41)+.21*r_m2(48)+.5*r_m2(51)
2157 d_m2(indx(iald2_z))= 0.0
2158 !
2159 p_m2(indx(ihcooh_z))= .52*r_m2(10)+.22*r_m2(12)
2160 d_m2(indx(ihcooh_z))= 0.0
2161 !
2162 p_m2(indx(iaone_z))= .07*r_m2(13)+.23*r_m2(15)+.74*r_m2(27) &
2163 +.74*r_m2(30)+.62*r_m2(33)+.74*r_m2(38) &
2164 +.57*r_m2(48)+.15*r_m2(50)
2165 d_m2(indx(iaone_z))= r_m2(5)+r_m2(6)
2166 !
2167 p_m2(indx(imgly_z))= .04*r_m2(12)+.07*r_m2(13)+.8*r_m2(19) &
2168 +.2*r_m2(26)+.19*r_m2(28)+.15*r_m2(50)
2169 d_m2(indx(imgly_z))= r_m2(7)+r_m2(8)+r_m2(9)
2170 !
2171 p_m2(indx(ieth_z))= 0.0
2172 d_m2(indx(ieth_z))= r_m2(10)+r_m2(11)
2173 !
2174 p_m2(indx(iolet_z))= 0.0
2175 d_m2(indx(iolet_z))= r_m2(12)+r_m2(14)+r_m2(16)
2176 !
2177 p_m2(indx(iolei_z))= 0.0
2178 d_m2(indx(iolei_z))= r_m2(13)+r_m2(15)+r_m2(17)
2179 !
2180 p_m2(indx(itol_z))= 0.0
2181 d_m2(indx(itol_z))= r_m2(18)
2182 !
2183 p_m2(indx(ixyl_z))= 0.0
2184 d_m2(indx(ixyl_z))= r_m2(19)
2185 !
2186 p_m2(indx(icres_z))= .12*r_m2(18)+.05*r_m2(19)
2187 d_m2(indx(icres_z))= r_m2(21)+r_m2(22)
2188 !
2189 p_m2(indx(ito2_z))= .8*r_m2(18)+.45*r_m2(19)
2190 d_m2(indx(ito2_z))= r_m2(20)
2191 !
2192 p_m2(indx(icro_z))= .4*r_m2(21)+r_m2(22)
2193 d_m2(indx(icro_z))= r_m2(23)
2194 !
2195 p_m2(indx(iopen_z))= .95*r_m2(20)+.3*r_m2(21)
2196 d_m2(indx(iopen_z))= r_m2(24)+r_m2(25)+r_m2(26)
2197 !
2198 p_m2(indx(ionit_z))= .05*r_m2(20)+r_m2(23)+.16*r_m2(33) &
2199 +.5*r_m2(36)+.5*r_m2(41)+r_m2(46)+.5*r_m2(51)
2200 d_m2(indx(ionit_z))= r_m2(29)+r_m2(30)
2201 !
2202 p_m2(indx(ipan_z))= 0.0
2203 d_m2(indx(ipan_z))= 0.0
2204 !
2205 p_m2(indx(ircooh_z))= .09*r_m2(12)+.16*r_m2(13)
2206 d_m2(indx(ircooh_z))= 0.0
2207 !
2208 p_m2(indx(irooh_z))= r_m2(43)+r_m2(45)
2209 d_m2(indx(irooh_z))= r_m2(27)+r_m2(28)
2210 !
2211 p_m2(indx(ich3o2_z))= +r_m2(5)+.07*r_m2(12)+.10*r_m2(13)
2212 d_m2(indx(ich3o2_z))= 0.0
2213 !
2214 p_m2(indx(iethp_z))= .06*r_m2(12)+.05*r_m2(13)+.1*r_m2(27) &
2215 +.1*r_m2(30)+.08*r_m2(33)+.1*r_m2(38)+.06*r_m2(48)
2216 d_m2(indx(iethp_z))= 0.0
2217 !
2218 p_m2(indx(ic2o3_z))= +r_m2(5)+r_m2(7)+r_m2(8) &
2219 +r_m2(9)+.13*r_m2(12)+.19*r_m2(13)+r_m2(24) &
2220 +r_m2(25)+.62*r_m2(26) +r_m2(35) &
2221 +r_m2(40)+.7*r_m2(50)
2222 d_m2(indx(ic2o3_z))= 0.0
2223 !
2224 p_m2(indx(iro2_z))= r_m2(1)+.03*r_m2(12)+.09*r_m2(13)+.77*r_m2(28)
2225 d_m2(indx(iro2_z))= r_m2(33)+r_m2(38)+r_m2(43)+r_m2(48)
2226 !
2227 p_m2(indx(iano2_z))= r_m2(6)+.11*r_m2(13)
2228 d_m2(indx(iano2_z))= r_m2(35)+r_m2(40)+r_m2(45)+r_m2(50)
2229 !
2230 p_m2(indx(inap_z))= r_m2(16)+r_m2(17)+r_m2(29)
2231 d_m2(indx(inap_z))= r_m2(36)+r_m2(41)+r_m2(46)+r_m2(51)
2232 !
2233 p_m2(indx(ixo2_z))= r_m2(8)+r_m2(11)+r_m2(14)+r_m2(15) &
2234 +.08*r_m2(18)+.5*r_m2(19)+.6*r_m2(21) &
2235 +r_m2(24)+.03*r_m2(26)+.4*r_m2(27)+.4*r_m2(30) &
2236 +.34*r_m2(33)+.4*r_m2(38)+.24*r_m2(48)
2237 d_m2(indx(ixo2_z))= r_m2(37)+r_m2(42)+r_m2(47)+r_m2(52)
2238
2239 p_m2(indx(ixpar_z))= 1.06*r_m2(12)+2.26*r_m2(13)+r_m2(14) &
2240 +2.23*r_m2(15)+1.98*r_m2(27)+.42*r_m2(28) &
2241 +1.98*r_m2(30)+1.68*r_m2(33)+r_m2(36) &
2242 +1.98*r_m2(38)+r_m2(41)+1.25*r_m2(48)+r_m2(51)
2243 d_m2(indx(ixpar_z))= r_m2(53)
2244 !
2245 return
2246 end subroutine gasode_m2
2247
2248
2249
2250 !***********************************************************************
2251 ! <18.> subr gasode_m3
2252 !
2253 ! purpose: updates production and destruction rates for mechanism 3
2254 ! isoprene
2255 !
2256 ! author : Rahul A. Zaveri
2257 ! date : December, 1998
2258 !
2259 ! ----------------------------------------------------------------------
2260
2261 subroutine gasode_m3( indx, r_m3, p_m3, d_m3 )
2262
2263 use module_data_cbmz
2264 implicit none
2265
2266 ! subr arguments
2267 integer indx(ngas_z)
2268 real r_m3(nrxn_m3), p_m3(ngas_tot), d_m3(ngas_tot)
2269
2270
2271 p_m3(indx(ino_z))= 0.0
2272 d_m3(indx(ino_z))= r_m3(8)+r_m3(9)+r_m3(10)
2273 !
2274 p_m3(indx(ino2_z))= .91*r_m3(8)+1.2*r_m3(9)+r_m3(10)
2275 d_m3(indx(ino2_z))= 0.0
2276 !
2277 p_m3(indx(ino3_z))= 0.0
2278 d_m3(indx(ino3_z))= r_m3(3)+r_m3(7)
2279 !
2280 p_m3(indx(ihno3_z))= .07*r_m3(7)
2281 d_m3(indx(ihno3_z))= 0.0
2282 !
2283 p_m3(indx(io3_z))= 0.0
2284 d_m3(indx(io3_z))= r_m3(2)+r_m3(6)
2285 !
2286 p_m3(indx(ioh_z))= .27*r_m3(2)+.27*r_m3(6)
2287 d_m3(indx(ioh_z))= r_m3(1)+r_m3(5)
2288 !
2289 p_m3(indx(iho2_z))= .07*r_m3(2)+.33*r_m3(4)+.1*r_m3(6)+.93*r_m3(7) &
2290 +.91*r_m3(8)+.8*r_m3(9)+r_m3(10)
2291 d_m3(indx(iho2_z))= r_m3(11)+r_m3(12)+r_m3(13)
2292 !
2293 p_m3(indx(ico_z))= .07*r_m3(2)+.33*r_m3(4)+.16*r_m3(6)+.64*r_m3(7) &
2294 +.59*r_m3(10)
2295 d_m3(indx(ico_z))= 0.0
2296 !
2297 p_m3(indx(ipar_z))= 1.86*r_m3(7)+0.18*r_m3(8)+1.6*r_m3(9)+2*r_m3(12) &
2298 +2*r_m3(15)
2299 d_m3(indx(ipar_z))= 0.0
2300 !
2301 p_m3(indx(ihcho_z))= .6*r_m3(2)+.2*r_m3(4)+.15*r_m3(6)+.28*r_m3(7) &
2302 +.63*r_m3(8)+.25*r_m3(10)
2303 d_m3(indx(ihcho_z))= 0.0
2304 !
2305 p_m3(indx(iald2_z))= .15*r_m3(2)+.07*r_m3(4)+.02*r_m3(6)+.28*r_m3(7) &
2306 +.8*r_m3(9)+.55*r_m3(10)+r_m3(15)+.5*r_m3(16)
2307 d_m3(indx(iald2_z))= 0.0
2308 !
2309 p_m3(indx(iaone_z))= .03*r_m3(4)+.09*r_m3(6)+.63*r_m3(10)+.5*r_m3(16)
2310 d_m3(indx(iaone_z))= 0.0
2311 !
2312 p_m3(indx(imgly_z))= .85*r_m3(6)+.34*r_m3(10)
2313 d_m3(indx(imgly_z))= 0.0
2314 !
2315 p_m3(indx(ionit_z))= .93*r_m3(7)+.09*r_m3(8)+.8*r_m3(9)+r_m3(12) &
2316 +r_m3(15)
2317 d_m3(indx(ionit_z))= 0.0
2318 !
2319 p_m3(indx(ihcooh_z))= .39*r_m3(2)+.46*r_m3(6)
2320 d_m3(indx(ihcooh_z))= 0.0
2321 !
2322 p_m3(indx(irooh_z))= r_m3(11)+r_m3(13)
2323 d_m3(indx(irooh_z))= 0.0
2324 !
2325 p_m3(indx(ich3o2_z))= .7*r_m3(4)+.05*r_m3(6)
2326 d_m3(indx(ich3o2_z))= 0.0
2327 !
2328 p_m3(indx(ic2o3_z))= .2*r_m3(2)+.97*r_m3(4)+.5*r_m3(5)+.11*r_m3(6) &
2329 +.07*r_m3(7)
2330 d_m3(indx(ic2o3_z))= 0.0
2331 !
2332 p_m3(indx(ixo2_z))= .08*r_m3(1)+.2*r_m3(2)+.2*r_m3(5)+.07*r_m3(6) &
2333 +.93*r_m3(7)
2334 d_m3(indx(ixo2_z))= 0.0
2335 !
2336 p_m3(indx(iisop_z))= 0.0
2337 d_m3(indx(iisop_z))= r_m3(1)+r_m3(2)+r_m3(3)
2338 !
2339 p_m3(indx(iisoprd_z))= .65*r_m3(2)+.91*r_m3(8)+.2*r_m3(9)+r_m3(14)
2340 d_m3(indx(iisoprd_z))= r_m3(4)+r_m3(5)+r_m3(6)+r_m3(7)
2341 !
2342 p_m3(indx(iisopp_z))= r_m3(1)
2343 d_m3(indx(iisopp_z))= r_m3(8)+r_m3(11)+r_m3(14)
2344 !
2345 p_m3(indx(iisopn_z))= r_m3(3)
2346 d_m3(indx(iisopn_z))= r_m3(9)+r_m3(12)+r_m3(15)
2347 !
2348 p_m3(indx(iisopo2_z))= .5*r_m3(5)
2349 d_m3(indx(iisopo2_z))= r_m3(10)+r_m3(13)+r_m3(16)
2350 !
2351 return
2352 end subroutine gasode_m3
2353
2354
2355
2356 !***********************************************************************
2357 ! <19.> subr gasode_m4
2358 !
2359 ! purpose: updates production and destruction rates for mechanism 4
2360 ! dimethylsulfide
2361 !
2362 ! author : Rahul A. Zaveri
2363 ! date : December, 1998
2364 !
2365 ! ----------------------------------------------------------------------
2366
2367 subroutine gasode_m4( indx, r_m4, p_m4, d_m4 )
2368
2369 use module_data_cbmz
2370 implicit none
2371
2372 ! subr arguments
2373 integer indx(ngas_z)
2374 real r_m4(nrxn_m4), p_m4(ngas_tot), d_m4(ngas_tot)
2375
2376
2377 p_m4(indx(ino_z))= r_m4(19)
2378 d_m4(indx(ino_z))= r_m4(5)+r_m4(11)+r_m4(26)+r_m4(30)
2379 !
2380 p_m4(indx(ino2_z))= r_m4(5)+r_m4(11)+r_m4(26)
2381 d_m4(indx(ino2_z))= r_m4(19)+r_m4(29)
2382 !
2383 p_m4(indx(ino3_z))= 0.0
2384 d_m4(indx(ino3_z))= r_m4(2)+r_m4(14)
2385 !
2386 p_m4(indx(ihono_z))= r_m4(30)
2387 d_m4(indx(ihono_z))= 0.0
2388 !
2389 p_m4(indx(ihno3_z))= r_m4(2)+r_m4(14)+r_m4(29)
2390 d_m4(indx(ihno3_z))= 0.0
2391 !
2392 p_m4(indx(io3_z))= 0.0
2393 d_m4(indx(io3_z))= r_m4(20)
2394 !
2395 p_m4(indx(io3p_z))= 0.0
2396 d_m4(indx(io3p_z))= r_m4(3)
2397 !
2398 p_m4(indx(ioh_z))= r_m4(21)
2399 d_m4(indx(ioh_z))= r_m4(1)+r_m4(4)+r_m4(9)+r_m4(10)+r_m4(16)+r_m4(23)
2400 !
2401 p_m4(indx(iho2_z))= .965*r_m4(4)+r_m4(6)+.27*r_m4(9)+r_m4(12)+r_m4(22) &
2402 +r_m4(27)+r_m4(32)
2403 d_m4(indx(iho2_z))= r_m4(13)+r_m4(21)+r_m4(31)
2404 !
2405 p_m4(indx(ih2o2_z))= r_m4(13)
2406 d_m4(indx(ih2o2_z))= 0.0
2407 !
2408 p_m4(indx(ico_z))= r_m4(32)
2409 d_m4(indx(ico_z))= 0.0
2410 !
2411 p_m4(indx(iso2_z))= r_m4(18)
2412 d_m4(indx(iso2_z))= 0.0
2413 !
2414 p_m4(indx(ih2so4_z))= r_m4(28)
2415 d_m4(indx(ih2so4_z))= 0.0
2416 !
2417 p_m4(indx(ihcho_z))= r_m4(5)+2*r_m4(6)+r_m4(7)+r_m4(11)+2*r_m4(12) &
2418 +r_m4(22)+r_m4(27)
2419 d_m4(indx(ihcho_z))= r_m4(32)
2420 !
2421 p_m4(indx(ich3o2_z))= r_m4(3)+.035*r_m4(4)+.73*r_m4(9)+r_m4(18)+r_m4(28)
2422 d_m4(indx(ich3o2_z))= r_m4(6)+r_m4(12)+r_m4(15)+r_m4(22)+r_m4(27)
2423 !
2424 p_m4(indx(ich3ooh_z))= r_m4(15)
2425 d_m4(indx(ich3ooh_z))= 0.0
2426 !
2427 p_m4(indx(idms_z))= 0.0
2428 d_m4(indx(idms_z))= r_m4(1)+r_m4(2)+r_m4(3)+r_m4(4)
2429 !
2430 p_m4(indx(imsa_z))= r_m4(17)+r_m4(23)+r_m4(29)+r_m4(30)+r_m4(31)+r_m4(32)
2431 d_m4(indx(imsa_z))= 0.0
2432 !
2433 p_m4(indx(idmso_z))= .965*r_m4(4)
2434 d_m4(indx(idmso_z))= r_m4(9)
2435 !
2436 p_m4(indx(idmso2_z))= .27*r_m4(9)
2437 d_m4(indx(idmso2_z))= r_m4(10)
2438 !
2439 p_m4(indx(ich3so2h_z))= .73*r_m4(9)
2440 d_m4(indx(ich3so2h_z))= r_m4(13)+r_m4(14)+r_m4(15)+r_m4(16)+r_m4(17)
2441 !
2442 p_m4(indx(ich3sch2oo_z))= r_m4(1)+r_m4(2)
2443 d_m4(indx(ich3sch2oo_z))= r_m4(5)+r_m4(6)+r_m4(7)+r_m4(8)+r_m4(8)
2444 !
2445 p_m4(indx(ich3so2_z))= r_m4(3)+.035*r_m4(4)+r_m4(5)+r_m4(6)+r_m4(7) &
2446 +1.85*r_m4(8) &
2447 +r_m4(11)+r_m4(12)+r_m4(13)+r_m4(14)+r_m4(15) &
2448 +r_m4(16)+r_m4(17)+r_m4(25)
2449 d_m4(indx(ich3so2_z))= r_m4(7)+r_m4(18)+r_m4(19)+r_m4(20)+r_m4(21) &
2450 +r_m4(22)+r_m4(23)+r_m4(24)
2451 !
2452 p_m4(indx(ich3so3_z))= r_m4(7)+r_m4(19)+r_m4(20)+r_m4(21)+r_m4(22) &
2453 +r_m4(26)+r_m4(27)
2454 d_m4(indx(ich3so3_z))= r_m4(17)+r_m4(28)+r_m4(29)+r_m4(30)+r_m4(31) &
2455 +r_m4(32)
2456 !
2457 p_m4(indx(ich3so2oo_z))= r_m4(24)
2458 d_m4(indx(ich3so2oo_z))= r_m4(25)+r_m4(26)+r_m4(27)
2459 !
2460 p_m4(indx(ich3so2ch2oo_z))= r_m4(10)
2461 d_m4(indx(ich3so2ch2oo_z))= r_m4(11)+r_m4(12)
2462 !
2463 p_m4(indx(imtf_z))= .15*r_m4(8)
2464 d_m4(indx(imtf_z))= 0.0
2465 !
2466 return
2467 end subroutine gasode_m4
2468
2469
2470
2471 !***********************************************************************
2472 ! <20.> subr gasrate_m1
2473 !
2474 ! purpose: computes reaction rates for mechanism 1
2475 !
2476 ! author : Rahul A. Zaveri
2477 ! date : December 1998
2478 !
2479 !-------------------------------------------------------------------------
2480
2481 subroutine gasrate_m1( indx, s, r_m1, r_m2, rk_m1, rk_m2, &
2482 h2o, ch4, oxygen, nitrogen, hydrogen )
2483
2484 use module_data_cbmz
2485 implicit none
2486
2487 ! subr arguments
2488 integer indx(ngas_z)
2489 real s(ngas_tot), r_m1(nrxn_m1), r_m2(nrxn_m2)
2490 real rk_m1(nrxn_m1), rk_m2(nrxn_m2)
2491 real h2o, ch4, oxygen, nitrogen, hydrogen
2492
2493 r_m1(1) = rk_m1(1)*s(indx(ino2_z))
2494 r_m1(2) = rk_m1(2)*s(indx(ino3_z))
2495 r_m1(3) = rk_m1(3)*s(indx(ihono_z))
2496 r_m1(4) = rk_m1(4)*s(indx(ihno3_z))
2497 r_m1(5) = rk_m1(5)*s(indx(ihno4_z))
2498 r_m1(6) = rk_m1(6)*s(indx(in2o5_z))
2499 r_m1(7) = rk_m1(7)*s(indx(io3_z))
2500 r_m1(8) = rk_m1(8)*s(indx(io3_z))
2501 r_m1(9) = rk_m1(9)*s(indx(ih2o2_z))
2502 r_m1(10) = rk_m1(10)*s(indx(io1d_z))*oxygen
2503 r_m1(11) = rk_m1(11)*s(indx(io1d_z))*nitrogen
2504 r_m1(12) = rk_m1(12)*s(indx(io1d_z))*h2o
2505 r_m1(13) = rk_m1(13)*s(indx(io3p_z))*oxygen
2506 r_m1(14) = rk_m1(14)*s(indx(io3p_z))*s(indx(io3_z))
2507 r_m1(15) = rk_m1(15)*s(indx(io3p_z))*s(indx(ino2_z))
2508 r_m1(16) = rk_m1(16)*s(indx(io3p_z))*s(indx(ino2_z))
2509 r_m1(17) = rk_m1(17)*s(indx(io3p_z))*s(indx(ino_z))
2510 r_m1(18) = rk_m1(18)*s(indx(io3_z))*s(indx(ino_z))
2511 r_m1(19) = rk_m1(19)*s(indx(io3_z))*s(indx(ino2_z))
2512 r_m1(20) = rk_m1(20)*s(indx(io3_z))*s(indx(ioh_z))
2513 r_m1(21) = rk_m1(21)*s(indx(io3_z))*s(indx(iho2_z))
2514 r_m1(22) = rk_m1(22)*s(indx(ioh_z))*hydrogen
2515 r_m1(23) = rk_m1(23)*s(indx(ioh_z))*s(indx(ino_z))
2516 r_m1(24) = rk_m1(24)*s(indx(ioh_z))*s(indx(ino2_z))
2517 r_m1(25) = rk_m1(25)*s(indx(ioh_z))*s(indx(ino3_z))
2518 r_m1(26) = rk_m1(26)*s(indx(ioh_z))*s(indx(ihono_z))
2519 r_m1(27) = rk_m1(27)*s(indx(ioh_z))*s(indx(ihno3_z))
2520 r_m1(28) = rk_m1(28)*s(indx(ioh_z))*s(indx(ihno4_z))
2521 r_m1(29) = rk_m1(29)*s(indx(ioh_z))*s(indx(iho2_z))
2522 r_m1(30) = rk_m1(30)*s(indx(ioh_z))*s(indx(ih2o2_z))
2523 r_m1(31) = rk_m1(31)*s(indx(iho2_z))*s(indx(iho2_z))
2524 r_m1(32) = rk_m1(32)*s(indx(iho2_z))*s(indx(iho2_z))*h2o
2525 r_m1(33) = rk_m1(33)*s(indx(iho2_z))*s(indx(ino_z))
2526 r_m1(34) = rk_m1(34)*s(indx(iho2_z))*s(indx(ino2_z))
2527 r_m1(35) = rk_m1(35)*s(indx(iho2_z))*s(indx(ino2_z))
2528 r_m1(36) = rk_m1(36)*s(indx(ihno4_z))
2529 r_m1(37) = rk_m1(37)*s(indx(ino3_z))*s(indx(ino_z))
2530 r_m1(38) = rk_m1(38)*s(indx(ino3_z))*s(indx(ino2_z))
2531 r_m1(39) = rk_m1(39)*s(indx(ino3_z))*s(indx(ino2_z))
2532 r_m1(40) = rk_m1(40)*s(indx(ino3_z))*s(indx(ino3_z))
2533 r_m1(41) = rk_m1(41)*s(indx(ino3_z))*s(indx(iho2_z))
2534 r_m1(42) = rk_m1(42)*s(indx(in2o5_z))*h2o
2535 r_m1(43) = rk_m1(43)*s(indx(in2o5_z))
2536 r_m1(44) = rk_m1(44)*s(indx(ico_z))*s(indx(ioh_z))
2537 r_m1(45) = rk_m1(45)*s(indx(iso2_z))*s(indx(ioh_z))
2538 r_m1(46) = rk_m1(46)*s(indx(ioh_z))*ch4
2539 r_m1(47) = rk_m1(47)*s(indx(ic2h6_z))*s(indx(ioh_z))
2540 r_m1(48) = rk_m1(48)*s(indx(ich3oh_z))*s(indx(ioh_z))
2541 r_m1(49) = rk_m1(49)*s(indx(ihcho_z))
2542 r_m1(50) = rk_m1(50)*s(indx(ihcho_z))
2543 r_m1(51) = rk_m1(51)*s(indx(ihcho_z))*s(indx(ioh_z))
2544 r_m1(52) = rk_m1(52)*s(indx(ihcho_z))*s(indx(ino3_z))
2545 r_m1(53) = rk_m1(53)*s(indx(ich3ooh_z))
2546 r_m1(54) = rk_m1(54)*s(indx(iethooh_z))
2547 r_m1(55) = rk_m1(55)*s(indx(ich3ooh_z))*s(indx(ioh_z))
2548 r_m1(56) = rk_m1(56)*s(indx(iethooh_z))*s(indx(ioh_z))
2549 r_m1(57) = rk_m1(57)*s(indx(ich3o2_z))*s(indx(ino_z))
2550 r_m1(58) = rk_m1(58)*s(indx(iethp_z))*s(indx(ino_z))
2551 r_m1(59) = rk_m1(59)*s(indx(ich3o2_z))*s(indx(ino3_z))
2552 r_m1(60) = rk_m1(60)*s(indx(iethp_z))*s(indx(ino3_z))
2553 r_m1(61) = rk_m1(61)*s(indx(ich3o2_z))*s(indx(iho2_z))
2554 r_m1(62) = rk_m1(62)*s(indx(iethp_z))*s(indx(iho2_z))
2555 r_m1(63) = rk_m1(63)*s(indx(ich3o2_z))
2556 r_m1(64) = rk_m1(64)*s(indx(iethp_z))
2557 r_m1(65) = rk_m1(65)*s(indx(ic2h5oh_z))*s(indx(ioh_z))
2558
2559 r_m2(2) = rk_m2(2)*s(indx(iald2_z))
2560 r_m2(3) = rk_m2(3)*s(indx(iald2_z))*s(indx(ioh_z))
2561 r_m2(4) = rk_m2(4)*s(indx(iald2_z))*s(indx(ino3_z))
2562 r_m2(31) = rk_m2(31)*s(indx(ic2o3_z))*s(indx(ino2_z))
2563 r_m2(32) = rk_m2(32)*s(indx(ipan_z))
2564 r_m2(34) = rk_m2(34)*s(indx(ic2o3_z))*s(indx(ino_z))
2565 r_m2(39) = rk_m2(39)*s(indx(ic2o3_z))*s(indx(ino3_z))
2566 r_m2(44) = rk_m2(44)*s(indx(ic2o3_z))*s(indx(iho2_z))
2567 r_m2(49) = rk_m2(49)*s(indx(ic2o3_z))
2568
2569 return
2570 end subroutine gasrate_m1
2571
2572
2573
2574 !***********************************************************************
2575 ! <21.> subr gasrate_m2
2576 !
2577 ! purpose: computes reaction rates for mechanism 2
2578 !
2579 ! author : Rahul A. Zaveri
2580 ! date : December 1998
2581 !
2582 !-------------------------------------------------------------------------
2583 !
2584 subroutine gasrate_m2( indx, s, r_m2, rk_m2 )
2585
2586 use module_data_cbmz
2587 implicit none
2588
2589 ! subr arguments
2590 integer indx(ngas_z)
2591 real s(ngas_tot), r_m2(nrxn_m2), rk_m2(nrxn_m2)
2592
2593 r_m2(1) = rk_m2(1)*s(indx(ipar_z))*s(indx(ioh_z))
2594
2595 r_m2(5) = rk_m2(5)*s(indx(iaone_z))
2596 r_m2(6) = rk_m2(6)*s(indx(iaone_z))*s(indx(ioh_z))
2597 r_m2(7) = rk_m2(7)*s(indx(imgly_z))
2598 r_m2(8) = rk_m2(8)*s(indx(imgly_z))*s(indx(ioh_z))
2599 r_m2(9) = rk_m2(9)*s(indx(imgly_z))*s(indx(ino3_z))
2600 r_m2(10) = rk_m2(10)*s(indx(ieth_z))*s(indx(io3_z))
2601 r_m2(11) = rk_m2(11)*s(indx(ieth_z))*s(indx(ioh_z))
2602 r_m2(12) = rk_m2(12)*s(indx(iolet_z))*s(indx(io3_z))
2603 r_m2(13) = rk_m2(13)*s(indx(iolei_z))*s(indx(io3_z))
2604 r_m2(14) = rk_m2(14)*s(indx(iolet_z))*s(indx(ioh_z))
2605 r_m2(15) = rk_m2(15)*s(indx(iolei_z))*s(indx(ioh_z))
2606 r_m2(16) = rk_m2(16)*s(indx(iolet_z))*s(indx(ino3_z))
2607 r_m2(17) = rk_m2(17)*s(indx(iolei_z))*s(indx(ino3_z))
2608 r_m2(18) = rk_m2(18)*s(indx(itol_z))*s(indx(ioh_z))
2609 r_m2(19) = rk_m2(19)*s(indx(ixyl_z))*s(indx(ioh_z))
2610 r_m2(20) = rk_m2(20)*s(indx(ito2_z))*s(indx(ino_z))
2611 r_m2(21) = rk_m2(21)*s(indx(icres_z))*s(indx(ioh_z))
2612 r_m2(22) = rk_m2(22)*s(indx(icres_z))*s(indx(ino3_z))
2613 r_m2(23) = rk_m2(23)*s(indx(icro_z))*s(indx(ino2_z))
2614 r_m2(24) = rk_m2(24)*s(indx(iopen_z))*s(indx(ioh_z))
2615 r_m2(25) = rk_m2(25)*s(indx(iopen_z))
2616 r_m2(26) = rk_m2(26)*s(indx(iopen_z))*s(indx(io3_z))
2617 r_m2(27) = rk_m2(27)*s(indx(irooh_z))
2618 r_m2(28) = rk_m2(28)*s(indx(irooh_z))*s(indx(ioh_z))
2619 r_m2(29) = rk_m2(29)*s(indx(ionit_z))*s(indx(ioh_z))
2620 r_m2(30) = rk_m2(30)*s(indx(ionit_z))
2621
2622 r_m2(33) = rk_m2(33)*s(indx(iro2_z))*s(indx(ino_z))
2623
2624 r_m2(35) = rk_m2(35)*s(indx(iano2_z))*s(indx(ino_z))
2625 r_m2(36) = rk_m2(36)*s(indx(inap_z))*s(indx(ino_z))
2626 r_m2(37) = rk_m2(37)*s(indx(ixo2_z))*s(indx(ino_z))
2627 r_m2(38) = rk_m2(38)*s(indx(iro2_z))*s(indx(ino3_z))
2628
2629 r_m2(40) = rk_m2(40)*s(indx(iano2_z))*s(indx(ino3_z))
2630 r_m2(41) = rk_m2(41)*s(indx(inap_z))*s(indx(ino3_z))
2631 r_m2(42) = rk_m2(42)*s(indx(ixo2_z))*s(indx(ino3_z))
2632 r_m2(43) = rk_m2(43)*s(indx(iro2_z))*s(indx(iho2_z))
2633
2634 r_m2(45) = rk_m2(45)*s(indx(iano2_z))*s(indx(iho2_z))
2635 r_m2(46) = rk_m2(46)*s(indx(inap_z))*s(indx(iho2_z))
2636 r_m2(47) = rk_m2(47)*s(indx(ixo2_z))*s(indx(iho2_z))
2637 r_m2(48) = rk_m2(48)*s(indx(iro2_z))
2638
2639 r_m2(50) = rk_m2(50)*s(indx(iano2_z))
2640 r_m2(51) = rk_m2(51)*s(indx(inap_z))
2641 r_m2(52) = rk_m2(52)*s(indx(ixo2_z))
2642 r_m2(53) = rk_m2(53)*s(indx(ipar_z))*s(indx(ixpar_z))
2643 !
2644 return
2645 end subroutine gasrate_m2
2646
2647
2648
2649 !***********************************************************************
2650 ! <22.> subr gasrate_m3
2651 !
2652 ! purpose: computes reaction rates for mechanism 3
2653 !
2654 ! author : Rahul A. Zaveri
2655 ! date : December 1998
2656 !
2657 !-------------------------------------------------------------------------
2658 !
2659 subroutine gasrate_m3( indx, s, r_m3, rk_m3 )
2660
2661 use module_data_cbmz
2662 implicit none
2663
2664 ! subr arguments
2665 integer indx(ngas_z)
2666 real s(ngas_tot), r_m3(nrxn_m3), rk_m3(nrxn_m3)
2667
2668 r_m3(1) = rk_m3(1)*s(indx(iisop_z))*s(indx(ioh_z))
2669 r_m3(2) = rk_m3(2)*s(indx(iisop_z))*s(indx(io3_z))
2670 r_m3(3) = rk_m3(3)*s(indx(iisop_z))*s(indx(ino3_z))
2671 r_m3(4) = rk_m3(4)*s(indx(iisoprd_z))
2672 r_m3(5) = rk_m3(5)*s(indx(iisoprd_z))*s(indx(ioh_z))
2673 r_m3(6) = rk_m3(6)*s(indx(iisoprd_z))*s(indx(io3_z))
2674 r_m3(7) = rk_m3(7)*s(indx(iisoprd_z))*s(indx(ino3_z))
2675 r_m3(8) = rk_m3(8)*s(indx(iisopp_z))*s(indx(ino_z))
2676 r_m3(9) = rk_m3(9)*s(indx(iisopn_z))*s(indx(ino_z))
2677 r_m3(10) = rk_m3(10)*s(indx(iisopo2_z))*s(indx(ino_z))
2678 r_m3(11) = rk_m3(11)*s(indx(iisopp_z))*s(indx(iho2_z))
2679 r_m3(12) = rk_m3(12)*s(indx(iisopn_z))*s(indx(iho2_z))
2680 r_m3(13) = rk_m3(13)*s(indx(iisopo2_z))*s(indx(iho2_z))
2681 r_m3(14) = rk_m3(14)*s(indx(iisopp_z))
2682 r_m3(15) = rk_m3(15)*s(indx(iisopn_z))
2683 r_m3(16) = rk_m3(16)*s(indx(iisopo2_z))
2684
2685 return
2686 end subroutine gasrate_m3
2687
2688
2689
2690 !***********************************************************************
2691 ! <23.> subr gasrate_m4
2692 !
2693 ! purpose: computes reaction rates for mechanism 4
2694 !
2695 ! author : Rahul A. Zaveri
2696 ! date : December 1998
2697 !
2698 !-------------------------------------------------------------------------
2699 !
2700 subroutine gasrate_m4( indx, s, r_m4, rk_m4, oxygen )
2701
2702 use module_data_cbmz
2703 implicit none
2704
2705 ! subr arguments
2706 integer indx(ngas_z)
2707 real s(ngas_tot), r_m4(nrxn_m4), rk_m4(nrxn_m4)
2708 real oxygen
2709
2710 r_m4(1) = rk_m4(1)*s(indx(idms_z))*s(indx(ioh_z))
2711 r_m4(2) = rk_m4(2)*s(indx(idms_z))*s(indx(ino3_z))
2712 r_m4(3) = rk_m4(3)*s(indx(idms_z))*s(indx(io3p_z))
2713 r_m4(4) = rk_m4(4)*s(indx(idms_z))*s(indx(ioh_z))
2714 r_m4(5) = rk_m4(5)*s(indx(ich3sch2oo_z))*s(indx(ino_z))
2715 r_m4(6) = rk_m4(6)*s(indx(ich3sch2oo_z))*s(indx(ich3o2_z))
2716 r_m4(7) = rk_m4(7)*s(indx(ich3sch2oo_z))*s(indx(ich3so2_z))
2717 r_m4(8) = rk_m4(8)*s(indx(ich3sch2oo_z))*s(indx(ich3sch2oo_z))
2718 r_m4(9) = rk_m4(9)*s(indx(idmso_z))*s(indx(ioh_z))
2719 r_m4(10) = rk_m4(10)*s(indx(idmso2_z))*s(indx(ioh_z))
2720 r_m4(11) = rk_m4(11)*s(indx(ich3so2ch2oo_z))*s(indx(ino_z))
2721 r_m4(12) = rk_m4(12)*s(indx(ich3so2ch2oo_z))*s(indx(ich3o2_z))
2722 r_m4(13) = rk_m4(13)*s(indx(ich3so2h_z))*s(indx(iho2_z))
2723 r_m4(14) = rk_m4(14)*s(indx(ich3so2h_z))*s(indx(ino3_z))
2724 r_m4(15) = rk_m4(15)*s(indx(ich3so2h_z))*s(indx(ich3o2_z))
2725 r_m4(16) = rk_m4(16)*s(indx(ich3so2h_z))*s(indx(ioh_z))
2726 r_m4(17) = rk_m4(17)*s(indx(ich3so2h_z))*s(indx(ich3so3_z))
2727 r_m4(18) = rk_m4(18)*s(indx(ich3so2_z))
2728 r_m4(19) = rk_m4(19)*s(indx(ich3so2_z))*s(indx(ino2_z))
2729 r_m4(20) = rk_m4(20)*s(indx(ich3so2_z))*s(indx(io3_z))
2730 r_m4(21) = rk_m4(21)*s(indx(ich3so2_z))*s(indx(iho2_z))
2731 r_m4(22) = rk_m4(22)*s(indx(ich3so2_z))*s(indx(ich3o2_z))
2732 r_m4(23) = rk_m4(23)*s(indx(ich3so2_z))*s(indx(ioh_z))
2733 r_m4(24) = rk_m4(24)*s(indx(ich3so2_z))*oxygen
2734 r_m4(25) = rk_m4(25)*s(indx(ich3so2oo_z))
2735 r_m4(26) = rk_m4(26)*s(indx(ich3so2oo_z))*s(indx(ino_z))
2736 r_m4(27) = rk_m4(27)*s(indx(ich3so2oo_z))*s(indx(ich3o2_z))
2737 r_m4(28) = rk_m4(28)*s(indx(ich3so3_z))
2738 r_m4(29) = rk_m4(29)*s(indx(ich3so3_z))*s(indx(ino2_z))
2739 r_m4(30) = rk_m4(30)*s(indx(ich3so3_z))*s(indx(ino_z))
2740 r_m4(31) = rk_m4(31)*s(indx(ich3so3_z))*s(indx(iho2_z))
2741 r_m4(32) = rk_m4(32)*s(indx(ich3so3_z))*s(indx(ihcho_z))
2742
2743 return
2744 end subroutine gasrate_m4
2745
2746
2747
2748 !**************************************************************************
2749 ! <24.> subr loadperoxyparameters
2750 !
2751 ! purpose: loads thermal rate coefficients for peroxy-peroxy
2752 ! permutation reactions
2753 !
2754 ! author : Rahul A. Zaveri
2755 ! date : June 1998
2756 !
2757 ! nomenclature:
2758 ! Aperox = Pre-exponential factor (molec-cc-s)
2759 ! Bperox = activation energy (-E/R) (K)
2760 !
2761 !-------------------------------------------------------------------------
2762 subroutine loadperoxyparameters( Aperox, Bperox )
2763
2764 use module_data_cbmz
2765 implicit none
2766
2767 ! subr arguments
2768 real Aperox(nperox,nperox), Bperox(nperox,nperox)
2769
2770 ! local variables
2771 integer i, j
2772
2773 Aperox(jch3o2,jch3o2) = 2.5e-13
2774 Aperox(jethp,jethp) = 6.8e-14
2775 Aperox(jc2o3,jc2o3) = 2.9e-12
2776 Aperox(jano2,jano2) = 8.0e-12
2777 Aperox(jnap,jnap) = 1.0e-12
2778 Aperox(jro2,jro2) = 5.3e-16
2779 Aperox(jisopp,jisopp) = 3.1e-14
2780 Aperox(jisopn,jisopn) = 3.1e-14
2781 Aperox(jisopo2,jisopo2) = 3.1e-14
2782 Aperox(jxo2,jxo2) = 3.1e-14
2783
2784 Bperox(jch3o2,jch3o2) = 190.
2785 Bperox(jethp,jethp) = 0.0
2786 Bperox(jc2o3,jc2o3) = 500.
2787 Bperox(jano2,jano2) = 0.0
2788 Bperox(jnap,jnap) = 0.0
2789 Bperox(jro2,jro2) = 1980.
2790 Bperox(jisopp,jisopp) = 1000.
2791 Bperox(jisopn,jisopn) = 1000.
2792 Bperox(jisopo2,jisopo2) = 1000.
2793 Bperox(jxo2,jxo2) = 1000.
2794
2795 do i = 1, nperox
2796 do j = 1, nperox
2797 if(i.ne.j)then
2798 Aperox(i,j) = 2.0*sqrt(Aperox(i,i)*Aperox(j,j))
2799 Bperox(i,j) = 0.5*(Bperox(i,i) + Bperox(j,j))
2800 endif
2801 enddo
2802 enddo
2803
2804 ! except for
2805 Aperox(jc2o3,jch3o2) = 1.3e-12
2806 Aperox(jch3o2,jc2o3) = 1.3e-12
2807 Bperox(jc2o3,jch3o2) = 640.
2808 Bperox(jch3o2,jc2o3) = 640.
2809
2810 return
2811 end subroutine loadperoxyparameters
2812
2813
2814
2815
2816 !**************************************************************************
2817 ! <25.> subr peroxyrateconstants
2818 !
2819 ! purpose: computes parameterized thermal rate coefficients
2820 ! for the alkylperoxy radical permutation reactions
2821 ! for the entire mechanism.
2822 !
2823 ! author : Rahul A. Zaveri
2824 ! date : June 1998
2825 !
2826 ! nomenclature:
2827 ! rk_param = parameterized reaction rate constants (1/s)
2828 ! rk_perox = individual permutation reaction rate constants (molec-cc-s)
2829 ! te = ambient atmospheric temperature (K)
2830 !
2831 !-------------------------------------------------------------------------
2832 subroutine peroxyrateconstants( tempbox, cbox, &
2833 Aperox, Bperox, rk_param )
2834
2835 use module_data_cbmz
2836 implicit none
2837
2838 ! subr arguments
2839 real tempbox, cbox(ngas_z)
2840 real Aperox(nperox,nperox), Bperox(nperox,nperox), rk_param(nperox)
2841
2842 ! local variables
2843 integer i, j
2844 real te
2845 real sperox(nperox), rk_perox(nperox,nperox)
2846
2847
2848 te = tempbox
2849
2850 sperox(jch3o2) = cbox(ich3o2_z)
2851 sperox(jethp) = cbox(iethp_z)
2852 sperox(jro2) = cbox(iro2_z)
2853 sperox(jc2o3) = cbox(ic2o3_z)
2854 sperox(jano2) = cbox(iano2_z)
2855 sperox(jnap) = cbox(inap_z)
2856 sperox(jisopp) = cbox(iisopp_z)
2857 sperox(jisopn) = cbox(iisopn_z)
2858 sperox(jisopo2) = cbox(iisopo2_z)
2859 sperox(jxo2) = cbox(ixo2_z)
2860
2861 !
2862 ! initialize to zero
2863 do i = 1, nperox
2864 rk_param(i) = 0.0
2865 enddo
2866
2867 do i = 1, nperox
2868 do j = 1, nperox
2869 rk_perox(i,j) = arr( Aperox(i,j), Bperox(i,j), te )
2870 rk_param(i) = rk_param(i) + rk_perox(i,j)*sperox(j)
2871 enddo
2872 enddo
2873
2874 return
2875 end subroutine peroxyrateconstants
2876
2877
2878
2879 !***********************************************************************
2880 ! <26.> subr gasthermrk_m1
2881 !
2882 ! purpose: computes thermal reaction rate coefficients for
2883 ! mechanism 1
2884 !
2885 ! author : Rahul A. Zaveri
2886 ! date : December 1998
2887 !
2888 !-------------------------------------------------------------------------
2889
2890 subroutine gasthermrk_m1( tempbox, cair_mlc, &
2891 rk_photo, rk_param, rk_m1, rk_m2 )
2892
2893 use module_data_cbmz
2894 implicit none
2895
2896 ! subr arguments
2897 real tempbox, cair_mlc
2898 real rk_photo(nphoto), rk_param(nperox)
2899 real rk_m1(nrxn_m1), rk_m2(nrxn_m2)
2900 ! local variables
2901 integer i
2902 real rk0, rk2, rk3, rki, rko, rmm, rnn, te
2903 ! real arr, troe
2904
2905
2906 te = tempbox
2907
2908 rk_m1(1) = rk_photo(jphoto_no2)
2909 rk_m1(2) = rk_photo(jphoto_no3)
2910 rk_m1(3) = rk_photo(jphoto_hono)
2911 rk_m1(4) = rk_photo(jphoto_hno3)
2912 rk_m1(5) = rk_photo(jphoto_hno4)
2913 rk_m1(6) = rk_photo(jphoto_n2o5)
2914 rk_m1(7) = rk_photo(jphoto_o3a)
2915 rk_m1(8) = rk_photo(jphoto_o3b)
2916 rk_m1(9) = rk_photo(jphoto_h2o2)
2917 rk_m1(10) = arr(3.2e-11, 70., te)
2918 rk_m1(11) = arr(1.8e-11, 110., te)
2919 rk_m1(12) = 2.2e-10
2920 rk_m1(13) = cair_mlc*6.e-34*(te/300.)**(-2.3)
2921 rk_m1(14) = arr(8.0e-12, -2060., te)
2922 rk_m1(15) = arr(6.5e-12, -120., te)
2923 !
2924 rk0 = 9.0e-32
2925 rnn = 2.0
2926 rki = 2.2e-11
2927 rmm = 0.0
2928 rk_m1(16) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
2929 !
2930 rk0 = 9.0e-32
2931 rnn = 1.5
2932 rki = 3.0e-11
2933 rmm = 0.0
2934 rk_m1(17) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
2935 !
2936 rk_m1(18) = arr(2.0e-12, -1400., te)
2937 rk_m1(19) = arr(1.2e-13, -2450., te)
2938 rk_m1(20) = arr(1.6e-12, -940., te)
2939 rk_m1(21) = arr(1.1e-14, -500., te)
2940 rk_m1(22) = arr(5.5e-12, -2000., te)
2941 !
2942 rk0 = 7.0e-31
2943 rnn = 2.6
2944 rki = 3.6e-11
2945 rmm = 0.1
2946 rk_m1(23) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
2947 !
2948 rk0 = 2.5e-30
2949 rnn = 4.4
2950 rki = 1.6e-11
2951 rmm = 1.7
2952 rk_m1(24) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
2953 !
2954 rk_m1(25) = 2.2e-11
2955 rk_m1(26) = arr(1.8e-11, -390., te)
2956 rko = 7.2e-15 * exp(785./te)
2957 rk2 = 4.1e-16 * exp(1440./te)
2958 rk3 = 1.9e-33 * exp(725./te)*cair_mlc
2959 rk_m1(27) = rko + rk3/(1.+rk3/rk2)
2960 rk_m1(28) = arr(1.3e-12, 380., te)
2961 rk_m1(29) = arr(4.8e-11, 250., te)
2962 rk_m1(30) = arr(2.9e-12, -160., te)
2963 rk_m1(31) = 2.3e-13 * exp(600./te) + &
2964 1.7e-33 * exp(1000./te)*cair_mlc ! ho2 + ho2 --> h2o2
2965 rk_m1(32) = rk_m1(31)*1.4e-21*exp(2200./te) ! ho2 + ho2 + h2o --> h2o2
2966 rk_m1(33) = arr(3.5e-12, 250., te)
2967 !
2968 rk0 = 1.8e-31
2969 rnn = 3.2
2970 rki = 4.7e-12
2971 rmm = 1.4
2972 rk_m1(34) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
2973 !
2974 rk_m1(35) = 5.0e-16
2975 rk_m1(36) = rk_m1(34)*arr(4.8e26, -10900., te)
2976 rk_m1(37) = arr(1.5e-11, 170., te)
2977 rk_m1(38) = arr(4.5e-14, -1260., te)
2978 !
2979 rk0 = 2.2e-30
2980 rnn = 3.9
2981 rki = 1.5e-12
2982 rmm = 0.7
2983 rk_m1(39) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
2984 !
2985 rk_m1(40) = arr(8.5e-13, -2450., te)
2986 rk_m1(41) = 3.5e-12
2987 rk_m1(42) = 2.0e-21
2988 rk_m1(43) = rk_m1(39)*arr(3.7e26, -11000., te)
2989 rk_m1(44) = 1.5e-13 * (1.+8.18e-23*te*cair_mlc) ! co + oh --> ho2
2990
2991 rk0 = 3.0e-31
2992 rnn = 3.3
2993 rki = 1.5e-12
2994 rmm = 0.0
2995 rk_m1(45) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
2996
2997 rk_m1(46) = te**.667*arr(2.8e-14, -1575., te)
2998 rk_m1(47) = te**2*arr(1.5e-17, -492., te)
2999 rk_m1(48) = arr(6.7e-12, -600., te)
3000 rk_m1(49) = rk_photo(jphoto_hchoa) ! hcho + hv --> 2ho2 + co
3001 rk_m1(50) = rk_photo(jphoto_hchob) ! hcho + hv --> co
3002 rk_m1(51) = 1.0e-11
3003 rk_m1(52) = arr(3.4e-13, -1900., te)
3004 rk_m1(53) = rk_photo(jphoto_ch3ooh)
3005 rk_m1(54) = rk_photo(jphoto_ethooh)
3006 rk_m1(55) = arr(3.8e-12, 200., te)
3007 rk_m1(56) = arr(3.8e-12, 200., te)
3008 rk_m1(57) = arr(3.0e-12, 280., te)
3009 rk_m1(58) = arr(2.6e-12, 365., te)
3010 rk_m1(59) = 1.1e-12
3011 rk_m1(60) = 2.5e-12
3012 rk_m1(61) = arr(3.8e-13, 800., te)
3013 rk_m1(62) = arr(7.5e-13, 700., te)
3014 rk_m1(63) = rk_param(jch3o2)
3015 rk_m1(64) = rk_param(jethp)
3016 rk_m1(65) = arr(7.0e-12, -235.,te)
3017
3018 rk_m2(2) = rk_photo(jphoto_ald2)
3019 rk_m2(3) = arr(5.6e-12, 270., te)
3020 rk_m2(4) = arr(1.4e-12, -1900., te)
3021 !
3022 rk0 = 9.7e-29
3023 rnn = 5.6
3024 rki = 9.3e-12
3025 rmm = 1.5
3026 rk_m2(31) = troe(cair_mlc,te,rk0,rnn,rki,rmm)
3027 !
3028 rk_m2(32) = rk_m2(31)*arr(1.1e28, -14000., te)
3029 rk_m2(34) = arr(5.3e-12, 360., te)
3030 rk_m2(39) = 4.0e-12
3031 rk_m2(44) = arr(4.5e-13, 1000., te)
3032 rk_m2(49) = rk_param(jc2o3)
3033
3034 ! Heterogeneous reactions
3035 ! rk_m1(65) = rk_het(1) ! O3 -->
3036 ! rk_m1(66) = rk_het(2) ! HO2 --> 0.5H2O2
3037 ! rk_m1(67) = rk_het(3) ! NO2 --> 0.5HONO + 0.5HNO3
3038 ! rk_m1(68) = rk_het(4) ! N2O5 --> 2HNO3
3039 ! rk_m1(69) = rk_het(5) ! HNO3 --> NO2
3040 ! rk_m1(70) = rk_het(6) ! HNO3 --> NO
3041 ! rk_m1(71) = rk_het(7) ! NO3 --> NO + O2
3042
3043 !!$Turn off this check since initializing to 0's; wig 16-Nov-2007
3044 !!$! all rate constants but be >= 0
3045 !!$ do i = 1, nrxn_m1
3046 !!$ rk_m1(i) = max( rk_m1(i), 0.0 )
3047 !!$ end do
3048 !!$ do i = 1, nrxn_m2
3049 !!$ rk_m2(i) = max( rk_m2(i), 0.0 )
3050 !!$ end do
3051
3052 return
3053 end subroutine gasthermrk_m1
3054
3055
3056
3057 !***********************************************************************
3058 ! <27.> subr gasthermrk_m2
3059 !
3060 ! purpose: computes thermal reaction rate coefficients for
3061 ! mechanism 2
3062 !
3063 ! author : Rahul A. Zaveri
3064 ! date : December 1998
3065 !
3066 !-------------------------------------------------------------------------
3067
3068 subroutine gasthermrk_m2( tempbox, cair_mlc, &
3069 rk_photo, rk_param, rk_m2 )
3070
3071 use module_data_cbmz
3072 implicit none
3073
3074 ! subr arguments
3075 real tempbox, cair_mlc
3076 real rk_photo(nphoto), rk_param(nperox), rk_m2(nrxn_m2)
3077 ! local variables
3078 integer i
3079 real rk0, rki, rmm, rnn, te
3080 ! real arr, troe
3081
3082
3083 te = tempbox
3084
3085 rk_m2(1) = 8.1e-13
3086
3087 rk_m2(5) = rk_photo(jphoto_aone)
3088 rk_m2(6) = te**2*arr(5.3e-18, -230., te)
3089 rk_m2(7) = rk_photo(jphoto_mgly)
3090 rk_m2(8) = 1.7e-11
3091 rk_m2(9) = arr(1.4e-12, -1900., te)
3092 rk_m2(10) = arr(1.2e-14, -2630., te)
3093 !
3094 rk0 = 1.0e-28
3095 rnn = 0.8
3096 rki = 8.8e-12
3097 rmm = 0.0
3098 rk_m2(11) = troe(cair_mlc,te,rk0,rnn,rki,rmm)
3099 !
3100 rk_m2(12) = arr(4.2e-15, -1800., te)
3101 rk_m2(13) = arr(8.9e-16, -392., te)
3102 rk_m2(14) = arr(5.8e-12, 478., te)
3103 rk_m2(15) = arr(2.9e-11, 255., te)
3104 rk_m2(16) = arr(3.1e-13, -1010., te)
3105 rk_m2(17) = 2.5e-12
3106 rk_m2(18) = arr(2.1e-12, 322., te)
3107 rk_m2(19) = arr(1.7e-11, 116., te)
3108 rk_m2(20) = 8.1e-12
3109 rk_m2(21) = 4.1e-11
3110 rk_m2(22) = 2.2e-11
3111 rk_m2(23) = 1.4e-11
3112 rk_m2(24) = 3.0e-11
3113 rk_m2(25) = rk_photo(jphoto_open)
3114 rk_m2(26) = arr(5.4e-17, -500., te)
3115 rk_m2(27) = rk_photo(jphoto_rooh)
3116 rk_m2(28) = arr(3.8e-12, 200., te)
3117 rk_m2(29) = arr(1.6e-11, -540., te)
3118 rk_m2(30) = rk_photo(jphoto_onit)
3119
3120 rk_m2(33) = 4.0e-12
3121
3122 rk_m2(35) = 4.0e-12
3123 rk_m2(36) = 4.0e-12
3124 rk_m2(37) = 4.0e-12
3125 rk_m2(38) = 2.5e-12
3126
3127 rk_m2(40) = 1.2e-12
3128 rk_m2(41) = 4.0e-12
3129 rk_m2(42) = 2.5e-12
3130 rk_m2(43) = arr(1.7e-13, 1300., te)
3131
3132 rk_m2(45) = arr(1.2e-13, 1300., te)
3133 rk_m2(46) = arr(1.7e-13, 1300., te)
3134 rk_m2(47) = arr(1.7e-13, 1300., te)
3135 rk_m2(48) = rk_param(jro2)
3136
3137 rk_m2(50) = rk_param(jano2)
3138 rk_m2(51) = rk_param(jnap)
3139 rk_m2(52) = rk_param(jxo2)
3140 rk_m2(53) = 1.0e-11 ! XPAR + PAR -->
3141
3142 !!$Turn off this check since initializing to 0's; wig 16-Nov-2007
3143 !!$! all rate constants but be >= 0
3144 !!$ do i = 1, nrxn_m2
3145 !!$ rk_m2(i) = max( rk_m2(i), 0.0 )
3146 !!$ end do
3147
3148 return
3149 end subroutine gasthermrk_m2
3150
3151
3152
3153 !***********************************************************************
3154 ! <28.> subr gasthermrk_m3
3155 !
3156 ! purpose: computes thermal reaction rate coefficients for
3157 ! mechanism 3
3158 !
3159 ! author : Rahul A. Zaveri
3160 ! date : December 1998
3161 !
3162 !-------------------------------------------------------------------------
3163
3164 subroutine gasthermrk_m3( tempbox, cair_mlc, &
3165 rk_photo, rk_param, rk_m3 )
3166
3167 use module_data_cbmz
3168 implicit none
3169
3170 ! subr arguments
3171 real tempbox, cair_mlc
3172 real rk_photo(nphoto), rk_param(nperox), rk_m3(nrxn_m3)
3173 ! local variables
3174 integer i
3175 real te
3176 ! real arr
3177
3178
3179 te = tempbox
3180 !
3181 rk_m3(1) = arr(2.6e-11, 409., te)
3182 rk_m3(2) = arr(1.2e-14, -2013., te)
3183 rk_m3(3) = arr(3.0e-12, -446., te)
3184 rk_m3(4) = rk_photo(jphoto_isoprd)
3185 rk_m3(5) = 3.3e-11
3186 rk_m3(6) = 7.0e-18
3187 rk_m3(7) = 1.0e-15
3188 rk_m3(8) = 4.0e-12
3189 rk_m3(9) = 4.0e-12
3190 rk_m3(10) = 4.0e-12
3191 rk_m3(11) = arr(1.7e-13, 1300., te)
3192 rk_m3(12) = arr(1.7e-13, 1300., te)
3193 rk_m3(13) = arr(1.7e-13, 1300., te)
3194 rk_m3(14) = rk_param(jisopp)
3195 rk_m3(15) = rk_param(jisopn)
3196 rk_m3(16) = rk_param(jisopo2)
3197
3198 ! all rate constants but be >= 0
3199 do i = 1, nrxn_m3
3200 rk_m3(i) = max( rk_m3(i), 0.0 )
3201 end do
3202
3203 return
3204 end subroutine gasthermrk_m3
3205
3206
3207
3208 !***********************************************************************
3209 ! <29.> subr gasthermrk_m4
3210 !
3211 ! purpose: computes thermal reaction rate coefficients for
3212 ! mechanism 4
3213 !
3214 ! author : Rahul A. Zaveri
3215 ! date : December 1998
3216 !
3217 !-------------------------------------------------------------------------
3218
3219 subroutine gasthermrk_m4( tempbox, cair_mlc, &
3220 rk_photo, rk_param, rk_m4 )
3221
3222 use module_data_cbmz
3223 implicit none
3224
3225 ! subr arguments
3226 real tempbox, cair_mlc
3227 real rk_photo(nphoto), rk_param(nperox), rk_m4(nrxn_m4)
3228 ! local variables
3229 integer i
3230 real B_abs, B_add, rk_tot, rk_tot_den, rk_tot_num, te
3231 ! real arr
3232
3233
3234 te = tempbox
3235 !
3236 rk_m4(1) = arr(9.6e-12, -234., te) ! ch3sch3 + oh --> ch3sch2
3237 rk_m4(2) = arr(1.4e-13, 500., te)
3238 rk_m4(3) = arr(1.3e-11, 409., te)
3239
3240 ! Hynes et al. (1986)
3241 rk_tot_num = te * exp(-234./te) + &
3242 8.46e-10 * exp(7230./te) + &
3243 2.68e-10 * exp(7810./te)
3244 rk_tot_den = 1.04e+11 * te + 88.1 * exp(7460./te)
3245 rk_tot = rk_tot_num/rk_tot_den
3246 B_abs = rk_m4(1)/rk_tot
3247 B_add = 1. - B_abs
3248
3249 rk_m4(4) = B_add*rk_tot ! ch3sch3 + oh --> ch3s(oh)ch3
3250 rk_m4(5) = 8.0e-12
3251 rk_m4(6) = 1.8e-13
3252 rk_m4(7) = 2.5e-13
3253 rk_m4(8) = 8.6e-14
3254 rk_m4(9) = 5.8e-11
3255 rk_m4(10) = 1.0e-14
3256 rk_m4(11) = 5.0e-12
3257 rk_m4(12) = 1.8e-13
3258 rk_m4(13) = 1.0e-15
3259 rk_m4(14) = 1.0e-13
3260 rk_m4(15) = 1.0e-15
3261 rk_m4(16) = 1.6e-11
3262 rk_m4(17) = 1.0e-13
3263 rk_m4(18) = arr(2.5e-13, -8686., te)
3264 rk_m4(19) = 1.0e-14
3265 rk_m4(20) = 5.0e-15
3266 rk_m4(21) = 2.5e-13
3267 rk_m4(22) = 2.5e-13
3268 rk_m4(23) = 5.0e-11
3269 rk_m4(24) = 2.6e-18
3270 rk_m4(25) = 3.3
3271 rk_m4(26) = 1.0e-11
3272 rk_m4(27) = 5.5e-12
3273 rk_m4(28) = arr(2.0e17, -12626., te)
3274 rk_m4(29) = 3.0e-15
3275 rk_m4(30) = 3.0e-15
3276 rk_m4(31) = 5.0e-11
3277 rk_m4(32) = 1.6e-15
3278
3279 ! all rate constants but be >= 0
3280 do i = 1, nrxn_m4
3281 rk_m4(i) = max( rk_m4(i), 0.0 )
3282 end do
3283
3284 return
3285 end subroutine gasthermrk_m4
3286
3287
3288
3289
3290 !***********************************************************************
3291 ! <26.> subr hetrateconstants
3292 !
3293 ! purpose: computes heterogeneous reaction rate coefficients
3294 !
3295 ! author : Rahul A. Zaveri
3296 ! date : May 2000
3297 !
3298 !-------------------------------------------------------------------------
3299
3300 subroutine hetrateconstants
3301 implicit none
3302
3303 return
3304 end subroutine hetrateconstants
3305
3306
3307
3308 !***********************************************************************
3309 ! <31.> func troe
3310 !
3311 ! purpose: calculates Troe reaction rate coefficient
3312 !
3313 ! author : Rahul A. Zaveri
3314 ! date : December 1998
3315 !-----------------------------------------------------------------------
3316
3317 real function troe( cairmlc, te, rk0, rnn, rki, rmm )
3318 implicit none
3319 ! func parameters
3320 real cairmlc, te, rk0, rnn, rki, rmm
3321 ! local variables
3322 real expo
3323
3324 rk0 = rk0*cairmlc*(te/300.)**(-rnn)
3325 rki = rki*(te/300.)**(-rmm)
3326 expo= 1./(1. + (ALOG10(rk0/rki))**2)
3327 troe = (rk0*rki/(rk0+rki))*.6**expo
3328 return
3329 end function troe
3330
3331
3332
3333 !***********************************************************************
3334 ! <32.> func arr
3335 !
3336 ! purpose: calculates arrhenius rate coefficient
3337 !
3338 ! author : Rahul A. Zaveri
3339 ! date : December 1998
3340 !-----------------------------------------------------------------------
3341
3342 real function arr( aa, bb, te )
3343 implicit none
3344 ! func parameters
3345 real aa, bb, te
3346
3347 arr = aa*exp(bb/te)
3348 return
3349 end function arr
3350
3351
3352
3353 !***********************************************************************
3354 ! <xx.> subr mapgas_tofrom_host
3355 !
3356 ! purpose: maps gas species between cboxold/new and host arrays
3357 !
3358 ! author : R. C. Easter
3359 ! date : November, 2003
3360 !
3361 ! ----------------------------------------------------------------------
3362
3363 subroutine mapgas_tofrom_host( imap, &
3364 i_boxtest_units_convert, &
3365 it,jt,kt, ims,ime, jms,jme, kms,kme, &
3366 num_moist, num_chem, moist, chem, &
3367 t_phy, p_phy, rho_phy, &
3368 cbox, tempbox, pressbox, airdenbox, &
3369 cair_mlc, &
3370 h2o, ch4, oxygen, nitrogen, hydrogen )
3371
3372 use module_configure, only: &
3373 p_qv, &
3374 p_so2, p_sulf, p_no2, p_no, p_o3, &
3375 p_hno3, p_h2o2, p_ald, p_hcho, p_op1, &
3376 p_op2, p_paa, p_ora1, p_ora2, p_nh3, &
3377 p_n2o5, p_no3, p_pan, p_hc3, p_hc5, &
3378 p_hc8, p_eth, p_co, p_ol2, p_olt, &
3379 p_oli, p_tol, p_xyl, p_aco3, p_tpan, &
3380 p_hono, p_hno4, p_ket, p_gly, p_mgly, &
3381 p_dcb, p_onit, p_csl, p_iso, p_ho, &
3382 p_ho2, &
3383 p_hcl, p_ch3o2, p_ethp, p_ch3oh, p_c2h5oh, &
3384 p_par, p_to2, p_cro, p_open, p_op3, &
3385 p_c2o3, p_ro2, p_ano2, p_nap, p_xo2, &
3386 p_xpar, p_isoprd, p_isopp, p_isopn, p_isopo2, &
3387 p_dms, p_msa, p_dmso, p_dmso2, p_ch3so2h, &
3388 p_ch3sch2oo, p_ch3so2, p_ch3so3, p_ch3so2oo, p_ch3so2ch2oo, &
3389 p_mtf
3390 use module_data_cbmz
3391 implicit none
3392
3393 ! subr arguments
3394 INTEGER, INTENT(IN) :: imap, it,jt,kt, ims,ime, jms,jme, kms,kme, &
3395 num_moist, num_chem, i_boxtest_units_convert
3396 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
3397 INTENT(IN) :: moist
3398 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
3399 INTENT(INOUT) :: chem
3400 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3401 INTENT(IN) :: t_phy, & ! temperature
3402 p_phy, & ! air pressure (Pa)
3403 rho_phy ! air density (kg/m3)
3404 REAL, INTENT(INOUT) :: cbox(ngas_z)
3405 REAL, INTENT(INOUT) :: tempbox, pressbox, airdenbox
3406 REAL, INTENT(INOUT) :: cair_mlc
3407 REAL, INTENT(INOUT) :: h2o, ch4, oxygen, nitrogen, hydrogen
3408
3409 ! local variables
3410 integer l
3411 real factoraa
3412 real, parameter :: eps=0.622
3413
3414
3415 tempbox = t_phy(it,kt,jt)
3416 ! p_phy = (Pa); pressbox = (dynes/cm2)
3417 pressbox = p_phy(it,kt,jt)*10.0
3418 ! rho_phy = (kg_air/m3); airdenbox = (mole_air/cm3)
3419 airdenbox = rho_phy(it,kt,jt)/28.966e3
3420 if (i_boxtest_units_convert .eq. 10) then
3421 airdenbox = rho_phy(it,kt,jt)
3422 end if
3423
3424 if (imap .gt. 0) goto 2000
3425
3426 !
3427 ! imap==0 -- initial species mapping from host array to cboxold
3428 ! chem --> czz --> cbox
3429 !
3430 ! note: do not map nh3, hcl
3431 !
3432 cbox(:) = 0.0
3433
3434 ! cair_mlc = (molecules_air/cm3)
3435 cair_mlc = airdenbox*avognumkpp
3436
3437 ! moist = (kg_h2o/kg_air); czz = (mole_h2o/cm3); h2o = (molecules_h2o/cm3)
3438 ! czz(ih2o_z) = (moist(it,kt,jt,p_qv)/eps)*airdenbox
3439 ! if (i_boxtest_units_convert .eq. 10) then
3440 ! czz(ih2o_z) = moist(it,kt,jt,p_qv)*airdenbox
3441 ! end if
3442 ! h2o = czz(ih2o_z)*avognumkpp
3443 h2o = (moist(it,kt,jt,p_qv)/eps)*airdenbox
3444 if (i_boxtest_units_convert .eq. 10) then
3445 h2o = moist(it,kt,jt,p_qv)*airdenbox
3446 end if
3447 h2o = h2o*avognumkpp
3448
3449 ! czz(ich4_z) = 1.7e-6*airdenbox ! ch4 conc. in mol/cc
3450 ! ch4 = czz(ich4_z)*avognumkpp ! ch4 conc. in molec/cc
3451 ch4 = 1.7e-6*airdenbox*avognumkpp ! ch4 conc. in molec/cc
3452
3453 oxygen = 0.21*cair_mlc ! o2 conc. in molec/cc
3454 nitrogen = 0.79*cair_mlc ! n2 conc. in molec/cc
3455 hydrogen = 0.58e-6*cair_mlc ! h2 conc. in molec/cc
3456
3457 ! chem units = (ppm); czz units = (mole/cm3); cbox units = (molecules/cm3)
3458 factoraa = airdenbox*1.0e-6
3459 if (i_boxtest_units_convert .eq. 10) factoraa = airdenbox
3460 factoraa = factoraa*avognumkpp
3461
3462 cbox(iso2_z) = chem(it,kt,jt,p_so2)*factoraa
3463 cbox(ih2so4_z) = chem(it,kt,jt,p_sulf)*factoraa
3464 cbox(ino2_z) = chem(it,kt,jt,p_no2)*factoraa
3465 cbox(ino_z) = chem(it,kt,jt,p_no)*factoraa
3466 cbox(io3_z) = chem(it,kt,jt,p_o3)*factoraa
3467 cbox(ihno3_z) = chem(it,kt,jt,p_hno3)*factoraa
3468 cbox(ih2o2_z) = chem(it,kt,jt,p_h2o2)*factoraa
3469 cbox(iald2_z) = chem(it,kt,jt,p_ald)*factoraa
3470 cbox(ihcho_z) = chem(it,kt,jt,p_hcho)*factoraa
3471 cbox(ich3ooh_z) = chem(it,kt,jt,p_op1)*factoraa
3472 cbox(iethooh_z) = chem(it,kt,jt,p_op2)*factoraa
3473 cbox(ihcooh_z) = chem(it,kt,jt,p_ora1)*factoraa
3474 cbox(ircooh_z) = chem(it,kt,jt,p_ora2)*factoraa
3475 cbox(inh3_z) = chem(it,kt,jt,p_nh3)*factoraa
3476 cbox(in2o5_z) = chem(it,kt,jt,p_n2o5)*factoraa
3477 cbox(ino3_z) = chem(it,kt,jt,p_no3)*factoraa
3478 cbox(ipan_z) = chem(it,kt,jt,p_pan)*factoraa
3479 cbox(ic2h6_z) = chem(it,kt,jt,p_eth)*factoraa
3480 cbox(ico_z) = chem(it,kt,jt,p_co)*factoraa
3481 cbox(ieth_z) = chem(it,kt,jt,p_ol2)*factoraa
3482 cbox(iolet_z) = chem(it,kt,jt,p_olt)*factoraa
3483 cbox(iolei_z) = chem(it,kt,jt,p_oli)*factoraa
3484 cbox(itol_z) = chem(it,kt,jt,p_tol)*factoraa
3485 cbox(ixyl_z) = chem(it,kt,jt,p_xyl)*factoraa
3486 cbox(ihono_z) = chem(it,kt,jt,p_hono)*factoraa
3487 cbox(ihno4_z) = chem(it,kt,jt,p_hno4)*factoraa
3488 cbox(iaone_z) = chem(it,kt,jt,p_ket)*factoraa
3489 cbox(imgly_z) = chem(it,kt,jt,p_mgly)*factoraa
3490 cbox(ionit_z) = chem(it,kt,jt,p_onit)*factoraa
3491 cbox(icres_z) = chem(it,kt,jt,p_csl)*factoraa
3492 cbox(iisop_z) = chem(it,kt,jt,p_iso)*factoraa
3493 cbox(ioh_z) = chem(it,kt,jt,p_ho)*factoraa
3494 cbox(iho2_z) = chem(it,kt,jt,p_ho2)*factoraa
3495
3496 cbox(ihcl_z) = chem(it,kt,jt,p_hcl)*factoraa
3497 cbox(ich3o2_z) = chem(it,kt,jt,p_ch3o2)*factoraa
3498 cbox(iethp_z) = chem(it,kt,jt,p_ethp)*factoraa
3499 cbox(ich3oh_z) = chem(it,kt,jt,p_ch3oh)*factoraa
3500 cbox(ic2h5oh_z) = chem(it,kt,jt,p_c2h5oh)*factoraa
3501 cbox(ipar_z) = chem(it,kt,jt,p_par)*factoraa
3502 cbox(ito2_z) = chem(it,kt,jt,p_to2)*factoraa
3503 cbox(icro_z) = chem(it,kt,jt,p_cro)*factoraa
3504 cbox(iopen_z) = chem(it,kt,jt,p_open)*factoraa
3505 cbox(irooh_z) = chem(it,kt,jt,p_op3)*factoraa
3506 cbox(ic2o3_z) = chem(it,kt,jt,p_c2o3)*factoraa
3507 cbox(iro2_z) = chem(it,kt,jt,p_ro2)*factoraa
3508 cbox(iano2_z) = chem(it,kt,jt,p_ano2)*factoraa
3509 cbox(inap_z) = chem(it,kt,jt,p_nap)*factoraa
3510 cbox(ixo2_z) = chem(it,kt,jt,p_xo2)*factoraa
3511 cbox(ixpar_z) = chem(it,kt,jt,p_xpar)*factoraa
3512 cbox(iisoprd_z) = chem(it,kt,jt,p_isoprd)*factoraa
3513 cbox(iisopp_z) = chem(it,kt,jt,p_isopp)*factoraa
3514 cbox(iisopn_z) = chem(it,kt,jt,p_isopn)*factoraa
3515 cbox(iisopo2_z) = chem(it,kt,jt,p_isopo2)*factoraa
3516 cbox(idms_z) = chem(it,kt,jt,p_dms)*factoraa
3517 cbox(imsa_z) = chem(it,kt,jt,p_msa)*factoraa
3518 cbox(idmso_z) = chem(it,kt,jt,p_dmso)*factoraa
3519 cbox(idmso2_z) = chem(it,kt,jt,p_dmso2)*factoraa
3520 cbox(ich3so2h_z) = chem(it,kt,jt,p_ch3so2h)*factoraa
3521 cbox(ich3sch2oo_z) = chem(it,kt,jt,p_ch3sch2oo)*factoraa
3522 cbox(ich3so2_z) = chem(it,kt,jt,p_ch3so2)*factoraa
3523 cbox(ich3so3_z) = chem(it,kt,jt,p_ch3so3)*factoraa
3524 cbox(ich3so2oo_z) = chem(it,kt,jt,p_ch3so2oo)*factoraa
3525 cbox(ich3so2ch2oo_z) = chem(it,kt,jt,p_ch3so2ch2oo)*factoraa
3526 cbox(imtf_z) = chem(it,kt,jt,p_mtf)*factoraa
3527
3528 cbox(ih2o_z) = h2o
3529 cbox(ich4_z) = ch4
3530 cbox(io2_z) = oxygen
3531 cbox(in2_z) = nitrogen
3532 cbox(ih2_z) = hydrogen
3533
3534 return
3535
3536 !
3537 ! imap==1 -- final species mapping from cbox back to host array
3538 ! cbox --> czz --> chem
3539 !
3540 ! note1: do not map nh3, hcl, ch4
3541 !
3542 2000 continue
3543 ! chem = (ppm); czz = (mole/cm3); cbox = (molecules/cm3)
3544 factoraa = airdenbox*1.0e-6
3545 if (i_boxtest_units_convert .eq. 10) factoraa = airdenbox
3546 factoraa = factoraa*avognumkpp
3547
3548 chem(it,kt,jt,p_so2) = cbox(iso2_z)/factoraa
3549 chem(it,kt,jt,p_sulf) = cbox(ih2so4_z)/factoraa
3550 chem(it,kt,jt,p_no2) = cbox(ino2_z)/factoraa
3551 chem(it,kt,jt,p_no) = cbox(ino_z)/factoraa
3552 chem(it,kt,jt,p_o3) = cbox(io3_z)/factoraa
3553 chem(it,kt,jt,p_hno3) = cbox(ihno3_z)/factoraa
3554 chem(it,kt,jt,p_h2o2) = cbox(ih2o2_z)/factoraa
3555 chem(it,kt,jt,p_ald) = cbox(iald2_z)/factoraa
3556 chem(it,kt,jt,p_hcho) = cbox(ihcho_z)/factoraa
3557 chem(it,kt,jt,p_op1) = cbox(ich3ooh_z)/factoraa
3558 chem(it,kt,jt,p_op2) = cbox(iethooh_z)/factoraa
3559 chem(it,kt,jt,p_ora1) = cbox(ihcooh_z)/factoraa
3560 chem(it,kt,jt,p_ora2) = cbox(ircooh_z)/factoraa
3561 chem(it,kt,jt,p_nh3) = cbox(inh3_z)/factoraa
3562 chem(it,kt,jt,p_n2o5) = cbox(in2o5_z)/factoraa
3563 chem(it,kt,jt,p_no3) = cbox(ino3_z)/factoraa
3564 chem(it,kt,jt,p_pan) = cbox(ipan_z)/factoraa
3565 chem(it,kt,jt,p_eth) = cbox(ic2h6_z)/factoraa
3566 chem(it,kt,jt,p_co) = cbox(ico_z)/factoraa
3567 chem(it,kt,jt,p_ol2) = cbox(ieth_z)/factoraa
3568 chem(it,kt,jt,p_olt) = cbox(iolet_z)/factoraa
3569 chem(it,kt,jt,p_oli) = cbox(iolei_z)/factoraa
3570 chem(it,kt,jt,p_tol) = cbox(itol_z)/factoraa
3571 chem(it,kt,jt,p_xyl) = cbox(ixyl_z)/factoraa
3572 chem(it,kt,jt,p_hono) = cbox(ihono_z)/factoraa
3573 chem(it,kt,jt,p_hno4) = cbox(ihno4_z)/factoraa
3574 chem(it,kt,jt,p_ket) = cbox(iaone_z)/factoraa
3575 chem(it,kt,jt,p_mgly) = cbox(imgly_z)/factoraa
3576 chem(it,kt,jt,p_onit) = cbox(ionit_z)/factoraa
3577 chem(it,kt,jt,p_csl) = cbox(icres_z)/factoraa
3578 chem(it,kt,jt,p_iso) = cbox(iisop_z)/factoraa
3579 chem(it,kt,jt,p_ho) = cbox(ioh_z)/factoraa
3580 chem(it,kt,jt,p_ho2) = cbox(iho2_z)/factoraa
3581
3582 chem(it,kt,jt,p_hcl) = cbox(ihcl_z)/factoraa
3583 chem(it,kt,jt,p_ch3o2) = cbox(ich3o2_z)/factoraa
3584 chem(it,kt,jt,p_ethp) = cbox(iethp_z)/factoraa
3585 chem(it,kt,jt,p_ch3oh) = cbox(ich3oh_z)/factoraa
3586 chem(it,kt,jt,p_c2h5oh) = cbox(ic2h5oh_z)/factoraa
3587 chem(it,kt,jt,p_par) = cbox(ipar_z)/factoraa
3588 chem(it,kt,jt,p_to2) = cbox(ito2_z)/factoraa
3589 chem(it,kt,jt,p_cro) = cbox(icro_z)/factoraa
3590 chem(it,kt,jt,p_open) = cbox(iopen_z)/factoraa
3591 chem(it,kt,jt,p_op3) = cbox(irooh_z)/factoraa
3592 chem(it,kt,jt,p_c2o3) = cbox(ic2o3_z)/factoraa
3593 chem(it,kt,jt,p_ro2) = cbox(iro2_z)/factoraa
3594 chem(it,kt,jt,p_ano2) = cbox(iano2_z)/factoraa
3595 chem(it,kt,jt,p_nap) = cbox(inap_z)/factoraa
3596 chem(it,kt,jt,p_xo2) = cbox(ixo2_z)/factoraa
3597 chem(it,kt,jt,p_xpar) = cbox(ixpar_z)/factoraa
3598 chem(it,kt,jt,p_isoprd) = cbox(iisoprd_z)/factoraa
3599 chem(it,kt,jt,p_isopp) = cbox(iisopp_z)/factoraa
3600 chem(it,kt,jt,p_isopn) = cbox(iisopn_z)/factoraa
3601 chem(it,kt,jt,p_isopo2) = cbox(iisopo2_z)/factoraa
3602 chem(it,kt,jt,p_dms) = cbox(idms_z)/factoraa
3603 chem(it,kt,jt,p_msa) = cbox(imsa_z)/factoraa
3604 chem(it,kt,jt,p_dmso) = cbox(idmso_z)/factoraa
3605 chem(it,kt,jt,p_dmso2) = cbox(idmso2_z)/factoraa
3606 chem(it,kt,jt,p_ch3so2h) = cbox(ich3so2h_z)/factoraa
3607 chem(it,kt,jt,p_ch3sch2oo) = cbox(ich3sch2oo_z)/factoraa
3608 chem(it,kt,jt,p_ch3so2) = cbox(ich3so2_z)/factoraa
3609 chem(it,kt,jt,p_ch3so3) = cbox(ich3so3_z)/factoraa
3610 chem(it,kt,jt,p_ch3so2oo) = cbox(ich3so2oo_z)/factoraa
3611 chem(it,kt,jt,p_ch3so2ch2oo) = cbox(ich3so2ch2oo_z)/factoraa
3612 chem(it,kt,jt,p_mtf) = cbox(imtf_z)/factoraa
3613
3614 return
3615 end subroutine mapgas_tofrom_host
3616
3617
3618
3619 !***********************************************************************
3620 ! <xx.> subr set_gaschem_allowed_regimes
3621 !
3622 ! purpose: determines which gas-phase chemistry regimes are allowed based
3623 ! on which species are active in the simulation
3624 !
3625 ! author :
3626 ! date :
3627 !
3628 ! ----------------------------------------------------------------------
3629
3630 subroutine set_gaschem_allowed_regimes( lunerr, &
3631 igaschem_allowed_m1, igaschem_allowed_m2, &
3632 igaschem_allowed_m3, igaschem_allowed_m4 )
3633 !
3634 ! determines which gas-phase chemistry regimes are allowed based
3635 ! on which species are active in the simulation
3636 !
3637 use module_configure, only: &
3638 p_qv, &
3639 p_so2, p_sulf, p_no2, p_no, p_o3, &
3640 p_hno3, p_h2o2, p_ald, p_hcho, p_op1, &
3641 p_op2, p_paa, p_ora1, p_ora2, p_nh3, &
3642 p_n2o5, p_no3, p_pan, p_hc3, p_hc5, &
3643 p_hc8, p_eth, p_co, p_ol2, p_olt, &
3644 p_oli, p_tol, p_xyl, p_aco3, p_tpan, &
3645 p_hono, p_hno4, p_ket, p_gly, p_mgly, &
3646 p_dcb, p_onit, p_csl, p_iso, p_ho, &
3647 p_ho2, &
3648 p_hcl, p_ch3o2, p_ethp, p_ch3oh, p_c2h5oh, &
3649 p_par, p_to2, p_cro, p_open, p_op3, &
3650 p_c2o3, p_ro2, p_ano2, p_nap, p_xo2, &
3651 p_xpar, p_isoprd, p_isopp, p_isopn, p_isopo2, &
3652 p_dms, p_msa, p_dmso, p_dmso2, p_ch3so2h, &
3653 p_ch3sch2oo, p_ch3so2, p_ch3so3, p_ch3so2oo, p_ch3so2ch2oo, &
3654 p_mtf
3655 use module_state_description, only: param_first_scalar
3656 use module_data_cbmz
3657 implicit none
3658
3659 ! subr arguments
3660 integer lunerr
3661 integer igaschem_allowed_m1, igaschem_allowed_m2, &
3662 igaschem_allowed_m3, igaschem_allowed_m4
3663
3664 ! local variables
3665 integer nactive, ndum, p1st
3666 character*80 msg
3667
3668
3669 ! index for first "active" scalar (= 2)
3670 p1st = param_first_scalar
3671
3672 ! determine if regime 1 is allowed
3673 ! (note: p_xxx>1 if xxx is active, p_xxx=1 if inactive)
3674 if (p_qv .lt. p1st) then
3675 msg = '*** subr set_gaschem_allowed_regimes'
3676 call peg_message( lunerr, msg )
3677 msg = '*** water vapor IS NOT ACTIVE'
3678 call peg_message( lunerr, msg )
3679 call peg_error_fatal( lunerr, msg )
3680 end if
3681
3682 ! determine if regime 1 is allowed
3683 ! (note: p_xxx>1 if xxx is active, p_xxx=1 if inactive)
3684 nactive = 0
3685 ndum = 27
3686 if (p_no .ge. p1st) nactive = nactive + 1
3687 if (p_no2 .ge. p1st) nactive = nactive + 1
3688 if (p_no3 .ge. p1st) nactive = nactive + 1
3689 if (p_n2o5 .ge. p1st) nactive = nactive + 1
3690 if (p_hono .ge. p1st) nactive = nactive + 1
3691 if (p_hno3 .ge. p1st) nactive = nactive + 1
3692 if (p_hno4 .ge. p1st) nactive = nactive + 1
3693 if (p_o3 .ge. p1st) nactive = nactive + 1
3694 ! if (p_o1d .ge. p1st) nactive = nactive + 1
3695 ! if (p_o3p .ge. p1st) nactive = nactive + 1
3696 if (p_ho .ge. p1st) nactive = nactive + 1
3697 if (p_ho2 .ge. p1st) nactive = nactive + 1
3698 if (p_h2o2 .ge. p1st) nactive = nactive + 1
3699 if (p_co .ge. p1st) nactive = nactive + 1
3700 if (p_so2 .ge. p1st) nactive = nactive + 1
3701 if (p_sulf .ge. p1st) nactive = nactive + 1
3702 ! if (p_nh3 .ge. p1st) nactive = nactive + 1
3703 ! if (p_hcl .ge. p1st) nactive = nactive + 1
3704 if (p_eth .ge. p1st) nactive = nactive + 1
3705 if (p_ch3o2 .ge. p1st) nactive = nactive + 1
3706 if (p_ethp .ge. p1st) nactive = nactive + 1
3707 if (p_hcho .ge. p1st) nactive = nactive + 1
3708 if (p_ch3oh .ge. p1st) nactive = nactive + 1
3709 if (p_c2h5oh .ge. p1st) nactive = nactive + 1
3710 if (p_op1 .ge. p1st) nactive = nactive + 1
3711 if (p_op2 .ge. p1st) nactive = nactive + 1
3712 if (p_ald .ge. p1st) nactive = nactive + 1
3713 if (p_ora1 .ge. p1st) nactive = nactive + 1
3714 if (p_pan .ge. p1st) nactive = nactive + 1
3715 if (p_ora2 .ge. p1st) nactive = nactive + 1
3716 if (p_c2o3 .ge. p1st) nactive = nactive + 1
3717
3718 if (nactive .le. 0) then
3719 igaschem_allowed_m1 = 0
3720 else if (nactive .eq. ndum) then
3721 igaschem_allowed_m1 = 1
3722 else
3723 msg = '*** subr set_gaschem_allowed_regimes'
3724 call peg_message( lunerr, msg )
3725 write(msg,90200) 1, nactive, ndum
3726 call peg_message( lunerr, msg )
3727 call peg_error_fatal( lunerr, msg )
3728 end if
3729 90200 format( ' error for regime ', i1, ', nactive, nexpected = ', 2i5 )
3730
3731 ! determine if regime 2 is allowed
3732 nactive = 0
3733 ndum = 19
3734 if (p_par .ge. p1st) nactive = nactive + 1
3735 if (p_ket .ge. p1st) nactive = nactive + 1
3736 if (p_mgly .ge. p1st) nactive = nactive + 1
3737 if (p_ol2 .ge. p1st) nactive = nactive + 1
3738 if (p_olt .ge. p1st) nactive = nactive + 1
3739 if (p_oli .ge. p1st) nactive = nactive + 1
3740 if (p_tol .ge. p1st) nactive = nactive + 1
3741 if (p_xyl .ge. p1st) nactive = nactive + 1
3742 if (p_csl .ge. p1st) nactive = nactive + 1
3743 if (p_to2 .ge. p1st) nactive = nactive + 1
3744 if (p_cro .ge. p1st) nactive = nactive + 1
3745 if (p_open .ge. p1st) nactive = nactive + 1
3746 if (p_onit .ge. p1st) nactive = nactive + 1
3747 if (p_op3 .ge. p1st) nactive = nactive + 1
3748 if (p_ro2 .ge. p1st) nactive = nactive + 1
3749 if (p_ano2 .ge. p1st) nactive = nactive + 1
3750 if (p_nap .ge. p1st) nactive = nactive + 1
3751 if (p_xo2 .ge. p1st) nactive = nactive + 1
3752 if (p_xpar .ge. p1st) nactive = nactive + 1
3753 if (nactive .le. 0) then
3754 igaschem_allowed_m2 = 0
3755 else if (nactive .eq. ndum) then
3756 igaschem_allowed_m2 = 2
3757 else
3758 msg = '*** subr set_gaschem_allowed_regimes'
3759 call peg_message( lunerr, msg )
3760 write(msg,90200) 2, nactive, ndum
3761 call peg_message( lunerr, msg )
3762 call peg_error_fatal( lunerr, msg )
3763 end if
3764
3765 ! determine if regime 3 is allowed
3766 nactive = 0
3767 ndum = 5
3768 if (p_iso .ge. p1st) nactive = nactive + 1
3769 if (p_isoprd .ge. p1st) nactive = nactive + 1
3770 if (p_isopp .ge. p1st) nactive = nactive + 1
3771 if (p_isopn .ge. p1st) nactive = nactive + 1
3772 if (p_isopo2 .ge. p1st) nactive = nactive + 1
3773 if (nactive .le. 0) then
3774 igaschem_allowed_m3 = 0
3775 else if (nactive .eq. ndum) then
3776 igaschem_allowed_m3 = 3
3777 else
3778 msg = '*** subr set_gaschem_allowed_regimes'
3779 call peg_message( lunerr, msg )
3780 write(msg,90200) 3, nactive, ndum
3781 call peg_message( lunerr, msg )
3782 call peg_error_fatal( lunerr, msg )
3783 end if
3784
3785 ! determine if regime 4 is allowed
3786 nactive = 0
3787 ndum = 11
3788 if (p_dms .ge. p1st) nactive = nactive + 1
3789 if (p_msa .ge. p1st) nactive = nactive + 1
3790 if (p_dmso .ge. p1st) nactive = nactive + 1
3791 if (p_dmso2 .ge. p1st) nactive = nactive + 1
3792 if (p_ch3so2h .ge. p1st) nactive = nactive + 1
3793 if (p_ch3sch2oo .ge. p1st) nactive = nactive + 1
3794 if (p_ch3so2 .ge. p1st) nactive = nactive + 1
3795 if (p_ch3so3 .ge. p1st) nactive = nactive + 1
3796 if (p_ch3so2oo .ge. p1st) nactive = nactive + 1
3797 if (p_ch3so2ch2oo .ge. p1st) nactive = nactive + 1
3798 if (p_mtf .ge. p1st) nactive = nactive + 1
3799 if (nactive .le. 0) then
3800 igaschem_allowed_m4 = 0
3801 else if (nactive .eq. ndum) then
3802 igaschem_allowed_m4 = 4
3803 else
3804 msg = '*** subr set_gaschem_allowed_regimes'
3805 call peg_message( lunerr, msg )
3806 write(msg,90200) 4, nactive, ndum
3807 call peg_message( lunerr, msg )
3808 call peg_error_fatal( lunerr, msg )
3809 end if
3810
3811 ! regime 1 must always be allowed
3812 if (igaschem_allowed_m1 .le. 0) then
3813 msg = '*** subr set_gaschem_allowed_regimes'
3814 call peg_message( lunerr, msg )
3815 write(msg,90300) 1
3816 call peg_message( lunerr, msg )
3817 call peg_error_fatal( lunerr, msg )
3818 end if
3819 90300 format( ' regime ', i1, ' must always be allowed' )
3820
3821 ! if regime 2 is allowed, then regime 1 must be allowed
3822 if (igaschem_allowed_m2 .gt. 0) then
3823 if (igaschem_allowed_m1 .le. 0) then
3824 msg = '*** subr set_gaschem_allowed_regimes'
3825 call peg_message( lunerr, msg )
3826 write(msg,90400) 2, 1
3827 call peg_message( lunerr, msg )
3828 call peg_error_fatal( lunerr, msg )
3829 end if
3830 end if
3831 90400 format( ' regime ', i1, ' allowed BUT regime ', i1, ' unallowed' )
3832
3833 ! if regime 3 is allowed, then regimes 1&2 must be allowed
3834 if (igaschem_allowed_m3 .gt. 0) then
3835 if (igaschem_allowed_m1 .le. 0) then
3836 msg = '*** subr set_gaschem_allowed_regimes'
3837 call peg_message( lunerr, msg )
3838 write(msg,90400) 3, 1
3839 call peg_message( lunerr, msg )
3840 call peg_error_fatal( lunerr, msg )
3841 else if (igaschem_allowed_m2 .le. 0) then
3842 msg = '*** subr set_gaschem_allowed_regimes'
3843 call peg_message( lunerr, msg )
3844 write(msg,90400) 3, 2
3845 call peg_message( lunerr, msg )
3846 call peg_error_fatal( lunerr, msg )
3847 end if
3848 end if
3849
3850 ! if regime 4 is allowed, then regime 1 must be allowed
3851 if (igaschem_allowed_m4 .gt. 0) then
3852 if (igaschem_allowed_m1 .le. 0) then
3853 msg = '*** subr set_gaschem_allowed_regimes'
3854 call peg_message( lunerr, msg )
3855 write(msg,90400) 4, 1
3856 call peg_message( lunerr, msg )
3857 call peg_error_fatal( lunerr, msg )
3858 end if
3859 end if
3860
3861 return
3862 end subroutine set_gaschem_allowed_regimes
3863
3864
3865
3866 !***********************************************************************
3867 ! <xx.> subr gasphotoconstants
3868 !
3869 ! purpose: copy photolytic rate constants from host arrays to local array
3870 !
3871 !-----------------------------------------------------------------------
3872 subroutine gasphotoconstants( rk_photo, &
3873 i_boxtest_units_convert, &
3874 it,jt,kt, ims,ime, jms,jme, kms,kme, &
3875 ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
3876 ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
3877 ph_ch3o2h, ph_n2o5 )
3878 !
3879 ! copies photolytic rate constants from host arrays to local arrays
3880 ! note1: currently 8 rate constants are scaled to other rate constants
3881 ! as is done in zz06gasphotolysis.f
3882 ! note2: currently the n2o5 rate is set to zero
3883 !
3884 use module_data_cbmz
3885 implicit none
3886
3887 ! subr arguments
3888 integer it,jt,kt, ims,ime, jms,jme, kms,kme
3889 integer i_boxtest_units_convert
3890 real rk_photo(nphoto)
3891 real, dimension( ims:ime, kms:kme, jms:jme ) :: &
3892 ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
3893 ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
3894 ph_ch3o2h, ph_n2o5
3895
3896 ! local variables
3897 real ft
3898
3899
3900 rk_photo(:) = 0.0
3901
3902 ! these from wrf/madronnich rate constants
3903 rk_photo(jphoto_no2) = ph_no2(it,kt,jt)
3904 rk_photo(jphoto_no3) = ph_no3o(it,kt,jt) &
3905 + ph_no3o2(it,kt,jt)
3906 rk_photo(jphoto_o3a) = ph_o33p(it,kt,jt)
3907 rk_photo(jphoto_o3b) = ph_o31d(it,kt,jt)
3908 rk_photo(jphoto_hono) = ph_hno2(it,kt,jt)
3909 rk_photo(jphoto_hno3) = ph_hno3(it,kt,jt)
3910 rk_photo(jphoto_hno4) = ph_hno4(it,kt,jt)
3911 rk_photo(jphoto_h2o2) = ph_h2o2(it,kt,jt)
3912 rk_photo(jphoto_ch3ooh) = ph_ch3o2h(it,kt,jt)
3913 rk_photo(jphoto_hchoa) = ph_ch2or(it,kt,jt)
3914 rk_photo(jphoto_hchob) = ph_ch2om(it,kt,jt)
3915 rk_photo(jphoto_n2o5) = ph_n2o5(it,kt,jt)
3916
3917 ! these scaled to other rate constants
3918 rk_photo(jphoto_ethooh) = 0.7 *rk_photo(jphoto_h2o2)
3919 rk_photo(jphoto_ald2) = 4.6e-4*rk_photo(jphoto_no2)
3920 rk_photo(jphoto_aone) = 7.8e-5*rk_photo(jphoto_no2)
3921 rk_photo(jphoto_mgly) = 9.64 *rk_photo(jphoto_hchoa)
3922 rk_photo(jphoto_open) = 9.04 *rk_photo(jphoto_hchoa)
3923 rk_photo(jphoto_rooh) = 0.7 *rk_photo(jphoto_h2o2)
3924 rk_photo(jphoto_onit) = 1.0e-4*rk_photo(jphoto_no2)
3925 rk_photo(jphoto_isoprd) = .025 *rk_photo(jphoto_hchob)
3926
3927 ! convert from (1/min) to (1/s)
3928 ! (except when i_boxtest_units_convert = 10 or 20)
3929 ft = 60.0
3930 if (i_boxtest_units_convert .eq. 10) ft = 1.0
3931 if (i_boxtest_units_convert .eq. 20) ft = 1.0
3932 if (ft .ne. 1.0) then
3933 rk_photo(:) = rk_photo(:)/ft
3934 end if
3935
3936
3937 return
3938 end subroutine gasphotoconstants
3939
3940
3941
3942 !-----------------------------------------------------------------------
3943 end module module_cbmz