module_cmu_bulkaqchem.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 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************
9
10 module module_cmu_bulkaqchem
11
12
13 use module_peg_util, only: peg_error_fatal, peg_message
14
15
16 implicit none
17
18
19 contains
20
21
22
23 !**********************************************************************
24 !
25 ! revised aqueous-phase chemistry module for the 3-d acid deposition model
26 !
27 !**********************************************************************
28 !
29 ! last modification : june 2, 1998 by s. pandis
30 ! the algorithm was modified to speed up execution. two major changes
31 ! (1) a bulk aqueous-phase chemistry approach is used and the
32 ! results are distributed over the particle spectrum according
33 ! to water content
34 ! (2) the water content is assumed to be constant (an input to the code)
35 !
36 ! ****** for bulk ******
37 !
38 !-----------------------------------------------------------------------
39 !
40 ! This code was obtained from S. Pandis in July 2003.
41 ! It was converted to Fortran-90 and adapted for use in WRF-chem with
42 ! the MOSAIC aerosol modules by R. Easter (PNNL) in July 2005.
43 !
44 ! 27-oct-2005 rce
45 ! switched from svode (single prec) to dvode (double prec) solver
46 ! changes in fullequil & aqratesa to make the ph calc more robust,
47 ! and do a graceful exit (?) when it fails
48 ! 01-nov-2005 rce
49 ! this version is completely double precision, including arguments
50 ! 09-dec-2005 rce
51 ! when iradical_in=100/101/102, the aqueous radical species are
52 ! calculated directly by dvode, and hybrd is not used
53 ! when idecomp_hmsa_hso5 > 0 & iradical_in = 100/1/2
54 ! aqueous so5- is transferred to gas so2
55 ! aqueous so4- is transferred to aerosol so4=
56 ! (this is in addition to the previous actions
57 ! when idecomp_hmsa_hso5 > 0:
58 ! aqueous hso5 is transferred to gas so2
59 ! aqueous hmsa is transferred to gas so2 & hcho)
60 ! 15-may-2006 rce
61 ! deleted the call to hybrd, so iradical_in>0 results in iradical=100
62 !-----------------------------------------------------------------------
63
64
65
66
67 !----------------------------------------------------------------------
68 subroutine aqoperator1( &
69 istat_fatal, istat_warn, &
70 tbeg_sec, tend_sec, &
71 gas, aerosol, gas_aqfrac, &
72 temp, p, lwc, rh, &
73 co2_mixrat_in, photo_in, iradical_in, idecomp_hmsa_hso5, iaq, &
74 yaq_beg, yaq_end, ph_cmuaq_cur )
75
76 use module_data_cmu_bulkaqchem, only: &
77 akeq, akhen, &
78 co2_mixrat, iradical, &
79 kaqx_hmsa, kaqx_hso5m, kaqx_siv, &
80 mdiag_fullequil, mdiag_hybrd, &
81 mdiag_negconc, mdiag_rsrate, mdiag_svode, &
82 meqn1, meqn1max, naers, ngas, &
83 na4, naa, nac, nahmsa, nahso5, nan, nar, nas, &
84 ng4, nga, ngc, ngh2o2, nghcho, nghcooh, ngn, ngso2, &
85 photo, &
86 rideal, &
87 wh2o2, wh2so4, whcho, whcl, whcooh, whno3, &
88 wnh3, wso2, wmol
89
90 !
91 !.....inputs
92 !
93 ! tbeg_sec : start of integration (in sec)
94 ! tend_sec : end of integration time (in sec)
95 !
96 ! gas(ngas) : gas-phase concentration vector (in ppm)
97 ! aerosol(naers) : aqueous (aerosol) species concentration (in ug/m3)
98 !
99 ! temp : temperature (in k)
100 ! p : pressure (in atm)
101 ! lwc : liquid water content (in g/m3)
102 ! rh : relative humidity (in 0-1 scale)
103 ! iradical_in : flag: 1/0 means aqueous free radical chemistry is on/off
104 ! photo_in : factor applied to aqueous photochemical rates
105 !
106 ! idecomp_hmsa_hso5 : flag -- if > 0, then hso5 is decomposed to so2,
107 ! and hmsa is decomposed to so2 & hcho, at end of integration.
108 ! If 0, hso5 and hmsa are left as is.
109 ! iaq : flag: 1 at first call, 0 there after
110 !
111 !.....outputs
112 ! gas(ngas) : updated gas-phase concentrations (in ppm)
113 ! aerosol(naers) : updated aqueous species concentration (in ug/m3)
114 ! iaq : flag set to zero
115 !
116 ! yaq_beg(meqn1max) : initial yaq values computed from initial gas & aerosol
117 ! yaq_end(meqn1max) : final yaq values
118 ! ph_cmuaq_cur : final ph values
119 !
120 !.....variables
121 ! gcon(ngas) : gas-phase concentrations (in ppm) (through rpar to rates)
122 !
123 !
124 !.....differential equations are solved for the following species
125 !
126 ! total gas aqueous
127 !-----------------------------------------------------------------------
128 ! 1. formaldehyde total 1. so2 1. s(iv)
129 ! 2. formic acid total 2. h2o2 2. h2o2(aq)
130 ! 3. nh3 3. nitrate
131 ! 4. chloride
132 ! 5. ammonium
133 ! 6. sulfate
134 ! 7. hso5-
135 ! 8. hmsa
136 !
137 !
138 !.....y matrix structure everything is (ug/m3)
139 ! only 11 odes are solved for the whole distribution
140 !
141 ! yaq(1) = total formaldehyde
142 ! yaq(2) = total formic acid
143 ! yaq(3) = so2 (g)
144 ! yaq(4) = h2o2 (g)
145 ! yaq(5) = nh3 (g)
146 ! yaq(6) = so2 (aq)
147 ! yaq(7) = h2o2 (aq)
148 ! yaq(8) = nitrate (aq)
149 ! yaq(9) = chloride (aq)
150 ! yaq(10) = ammonium (aq)
151 ! yaq(11) = sulfate (aq)
152 ! yaq(12) = hso5- (aq)
153 ! yaq(13) = hmsa (aq)
154 !
155
156 ! arguments
157 integer istat_fatal, istat_warn
158 integer iradical_in, idecomp_hmsa_hso5, iaq
159 double precision tbeg_sec, tend_sec
160 double precision temp, p, lwc, rh
161 double precision co2_mixrat_in, photo_in
162 double precision gas(ngas), aerosol(naers)
163 double precision gas_aqfrac(ngas)
164 double precision yaq_beg(meqn1max), yaq_end(meqn1max), ph_cmuaq_cur
165
166 ! local variables
167 integer i, isp
168 integer icount, iprint, negconc_count
169
170 double precision &
171 ammonold, chlorold, crustal, &
172 dammon, dchlor, dhmsa, dhso5, dnit, dsulf, &
173 fammon, fh2o2, fh2so4, fhcho, fhcl, fhcooh, fhno3, fso2, &
174 hmsaold, hso5old, &
175 pres, &
176 salt, sodium, stime, stout, sulfateold, &
177 tnitold, watcont, watvap, &
178 whchowhmsa, wso2whmsa, wso2whso5, wso2wsiv, &
179 yaq_h2so4_g, yaq_hcl_g, yaq_hno3_g
180 double precision &
181 duma, dumb, dumhion, vlwc
182 double precision gascon(ngas)
183 double precision yaq(meqn1max)
184
185 integer ipar(8)
186 double precision rpar(8+ngas)
187
188 !
189 ! initialization (only on first call)
190 ! (calculation of sectional diameters (which isn't really needed),
191 ! and loading of molec. weights (which is needed))
192 !
193 if (iaq .eq. 1) then
194 call dropinit
195 iaq = 0
196 endif
197
198 mdiag_fullequil = 1
199 mdiag_hybrd = 1
200 mdiag_negconc = 1
201 mdiag_rsrate = 1
202 mdiag_svode = 1
203
204 !
205 ! calc temperature dependent parameters
206 !
207 pres = p ! pressure in atm
208 call constants(temp)
209
210 !
211 ! zero the yaq matrix
212 !
213 do i=1, meqn1max
214 yaq(i)=0.0
215 enddo
216
217 watcont = lwc
218
219 ! for iradical_in<0, set iradical=0
220 ! for iradical_in>0, set iradical=100 except when iradical_in=101,102
221 iradical = iradical_in
222 if (iradical .gt. 0) then
223 if ((iradical .ne. 101) .and. &
224 (iradical .ne. 102)) iradical = 100
225 else
226 iradical = 0
227 end if
228
229 photo = photo_in
230 co2_mixrat = co2_mixrat_in
231
232 meqn1 = 13
233 ! if (iradical .eq. 100) meqn1 = 20
234 if (iradical .gt. 0) meqn1 = 20
235
236 icount=0
237 ! iprint=20
238
239 !
240 ! transfer of all gas-phase concentrations to gascon and then
241 ! through rpar to rates
242 !
243 do i=1, ngas
244 gascon(i) = gas(i)
245 gas_aqfrac(i) = 0.0
246 enddo
247 !
248 ! unit conversion factors
249 !
250 fso2=1000.*pres*wso2/(rideal*temp) !so2 conver. factor from ppm to ug/m3
251 fh2o2=1000.*pres*wh2o2/(rideal*temp) !h2o2 conver. factor from ppm to ug/m3
252 fhcho=1000.*pres*whcho/(rideal*temp) !hcho conver. factor from ppm to ug/m3
253 fhcooh=1000.*pres*whcooh/(rideal*temp) !hcooh conver. factor from ppm to ug/m3
254 fammon=1000.*pres*wnh3/(rideal*temp) !nh3 conver. factor from ppm to ug/m3
255 !
256 !
257 ! *** beginning of aqueous-phase calculation
258 ! loading of the concentrations in the yaq vector
259 !
260 ! the total formaldehyde/formic acid are transferred as gas-phase species
261 ! (they are not included in the aqueous or aerosol matrices)
262 yaq(1) = gas(nghcho)*fhcho ! total hcho (ug/m3)
263 yaq(2) = gas(nghcooh)*fhcooh ! total hcooh (ug/m3)
264 !
265 ! gas-phase species
266 !
267 yaq(3) = gas(ngso2)*fso2 ! total so2 (ug/m3)
268 yaq(4) = gas(ngh2o2)*fh2o2 ! total h2o2 (ug/m3)
269 yaq(5) = gas(nga) *fammon ! nh3(g) in ug/m3
270
271 !
272 ! aqueous-phase species
273 !
274 ! "aerosol" array holds the bulk aqueous chemical concentrations
275 ! (in ug/m3-of-air), so there is no "activation" or summing over sections
276 !
277 ! ** the droplets are assumed to by without so2 or
278 ! h2o2 in the beginning of a timestep **
279 ! (this is done to avoid carrying the s(iv)/h2o2 concentrations
280 ! in the 3d simulation **
281 yaq(6) = 1.e-10 ! so2 (aq) in ug/m3
282 yaq(7) = 1.e-10 ! h2o2 (aq) in ug/m3
283
284 yaq(8) = aerosol(nan) ! nitrate (aq) in ug/m3
285 yaq(9) = aerosol(nac) ! chloride (aq) in ug/m3
286 yaq(10) = aerosol(naa) ! ammonium (aq) in ug/m3
287 yaq(11) = aerosol(na4) ! sulfate (aq) in ug/m3
288 yaq(12) = aerosol(nahso5) ! hso5- in ug/m3
289 yaq(13) = aerosol(nahmsa) ! hmsa in ug/m3
290
291 !
292 ! save the old aqueous-phase concentrations before the integration
293 !
294 tnitold = yaq(8)
295 chlorold = yaq(9)
296 ammonold = yaq(10)
297 sulfateold = yaq(11)
298 hso5old = yaq(12)
299 hmsaold = yaq(13)
300 !
301 ! calculation of the dissolution of the available h2so4, hno3 and hcl
302 ! to the droplets/particles. we assume that the timescale of dissolution
303 ! is small and that all hno3 and hcl are dissolved (zero vapor pressure)
304 !
305 fh2so4=1000.*wh2so4*pres/(rideal*temp) !h2so4 conver. factor from ppm to ug/m3
306 fhno3 =1000.*whno3*pres/(rideal*temp) !hno3 conver. factor from ppm to ug/m3
307 fhcl =1000.*whcl*pres/(rideal*temp) !hcl conver. factor from ppm to ug/m3
308 !
309 yaq_h2so4_g = gas(ng4)*fh2so4 ! h2so4(g) (in ug/m3)
310 yaq_hno3_g = gas(ngn)*fhno3 ! hno3(g) (in ug/m3)
311 yaq_hcl_g = gas(ngc)*fhcl ! hcl(g) (in ug/m3)
312 !
313 ! ** transfer the gas-phase h2so4, hno3 and hcl to aqueous phase **
314 ! (and account for molec-weight differences)
315 gas(ng4) = 0.0 ! all h2so4 is dissolved
316 gas(ngn) = 0.0 ! all hno3 is dissolved
317 gas(ngc) = 0.0 ! all hcl is dissolved
318 yaq(11) = yaq(11) + (yaq_h2so4_g/wh2so4)*wmol(2) ! so4-- increase
319 yaq(8) = yaq(8) + (yaq_hno3_g/whno3)*wmol(4) ! no3- increase
320 yaq(9) = yaq(9) + (yaq_hcl_g/whcl)*wmol(15) ! cl- increase
321
322 !
323 !.....the yaq matrix has been initialized
324 !
325
326 ! test to be removed
327 ! write(27,*)' initial values of the yaqs'
328 ! write(27,*)lwc,rh,temp
329 ! do i = 1, 11
330 ! write(27,*)i, yaq(i)
331 ! enddo
332
333 !
334 ! variables for integration
335 !
336 stime = tbeg_sec ! in seconds
337 stout = tend_sec ! in seconds
338 !
339 ! calculation of sodium (needed for the integration routine)
340 !
341 sodium = aerosol(nas)
342 !
343 ! calculate crustal species concentration inside droplets
344 ! (it is used in the aqueous-phase chemistry calculation to
345 ! estimate fe and mn concentrations
346 !
347 crustal = aerosol(nar)
348 salt = aerosol(nas)
349
350 !
351 ! save yaq to yaq_beg
352 !
353 do i = 1, 13
354 yaq_beg(i) = yaq(i)
355 end do
356
357 !
358 ! load ipar, rpar
359 !
360 ipar(1) = icount
361 ipar(2) = iprint
362 ipar(3) = 0
363 ipar(4) = 0
364 ipar(5) = 0
365 ipar(6) = 0
366 ipar(7) = 0
367 ipar(8) = 0
368
369 rpar(1) = temp
370 rpar(2) = pres
371 rpar(3) = watcont
372 rpar(4) = watvap
373 rpar(5) = crustal
374 rpar(6) = salt
375 rpar(7) = sodium
376 rpar(8) = ph_cmuaq_cur
377 do i = 1, ngas
378 rpar(8+i) = gascon(i)
379 end do
380
381 !
382 ! integrate
383 !
384 call aqintegr1( meqn1, yaq, stime, stout, rpar, ipar )
385
386 !
387 ! unload ipar, rpar
388 !
389 ph_cmuaq_cur = rpar(8)
390
391 !
392 ! calculate the differences
393 !
394 dnit = yaq(8) - tnitold
395 dchlor = yaq(9) - chlorold
396 dammon = yaq(10) - ammonold
397 dsulf = yaq(11) - sulfateold
398 dhso5 = yaq(12) - hso5old
399 dhmsa = yaq(13) - hmsaold
400
401 !
402 ! transfer new aqueous concentrations back to aerosol array
403 !
404 aerosol(nan) = yaq(8) ! nitrate (aq) in ug/m3
405 aerosol(nac) = yaq(9) ! chloride (aq) in ug/m3
406 aerosol(naa) = yaq(10) ! ammonium (aq) in ug/m3
407 aerosol(na4) = yaq(11) ! sulfate (aq) in ug/m3
408 aerosol(nahso5) = yaq(12) ! hso5 (aq) in ug/m3
409 aerosol(nahmsa) = yaq(13) ! hmsa (aq) in ug/m3
410
411 !
412 ! check if the algorithm resulted in negative concentrations
413 ! report the corrections at the warning file
414 !
415 negconc_count = 0
416 do isp=1, naers
417 if (aerosol(isp) .lt. 0.0) then
418 negconc_count = negconc_count + 1
419 if (mdiag_negconc .gt. 0) then
420 ! write(2,*)'negative concentration at', stime
421 ! write(2,*)'species=',isp, aerosol(isp)
422 if (negconc_count .eq. 1) write(6,*) &
423 '*** module_cmuaq_bulk aqoperator1 neg conc at t', stime
424 write(6,*) ' isp, conc', isp, aerosol(isp)
425 end if
426 aerosol(isp)=1.e-12
427 endif
428 enddo
429
430 !
431 ! return gas yaq results to their original matrices
432 !
433 ! ** gas-phase species **
434 ! important : the dissolved s(iv) and h2o2 are transferred
435 ! back to the gas phase
436 ! (this change is applied to "gas" but not to "yaq")
437 !
438 gas(nghcho) = yaq(1)/fhcho ! total hcho (ppm)
439 gas(nghcooh) = yaq(2)/fhcooh ! total hcooh (ppm)
440
441 wso2wsiv = wso2/wmol(kaqx_siv)
442 ! gas(ngso2) = (yaq(3)+yaq(6)*0.7901)/fso2 ! so2(g) (ppm)
443 gas(ngso2) = (yaq(3)+yaq(6)*wso2wsiv)/fso2 ! so2(g) (ppm)
444 gas_aqfrac(ngso2) = (yaq(6)*wso2wsiv) / (yaq(3)+yaq(6)*wso2wsiv)
445
446 gas(ngh2o2) = (yaq(4)+yaq(7))/fh2o2 ! h2o2(g) (ppm)
447 gas_aqfrac(ngh2o2) = yaq(7) / (yaq(4)+yaq(7))
448
449 gas(nga) = yaq(5)/fammon ! nh3(g) (ppm)
450
451 !
452 ! gas fractions [gas/(gas+aqueous)] for hcho and hcooh
453 ! hcho expression is from subr aqratesa (duma here = hc1 there)
454 ! hcooh expression is from subr electro (duma here = dform there)
455 !
456 vlwc = lwc*1.0e-6 ! (liter-h2o)/(liter-air)
457
458 duma = rideal*temp*vlwc*akhen(7)
459 gas_aqfrac(nghcho) = duma/(duma + 1.0)
460
461 dumb = max( 0.0d0, min( 14.0d0, ph_cmuaq_cur ) )
462 dumhion = 10.0**(-dumb)
463 duma = rideal*temp*vlwc*akhen(8)*(1.+akeq(13)/dumhion)
464 gas_aqfrac(nghcooh) = duma/(duma + 1.0)
465
466
467 !
468 ! if idecomp_hmsa_hso5 > 0, then
469 ! aqueous hso5 is transferred to gas so2
470 ! aqueous hmsa is transferred to gas so2 & hcho
471 ! (this change is applied to 'gas' & 'aerosol' but not to 'yaq')
472 ! also, if idecomp_hmsa_hso5 > 0 & iradical = 100/1/2
473 ! aqueous so5- is transferred to gas so2
474 ! aqueous so4- is transferred to aerosol so4=
475 !
476 if (idecomp_hmsa_hso5 .gt. 0) then
477 wso2whso5 = wso2/wmol(kaqx_hso5m)
478 wso2whmsa = wso2/wmol(kaqx_hmsa)
479 whchowhmsa = whcho/wmol(kaqx_hmsa)
480 gas(ngso2) = gas(ngso2) + &
481 (yaq(12)*wso2whso5 + yaq(13)*wso2whmsa)/fso2
482 gas(nghcho) = gas(nghcho) + yaq(13)*whchowhmsa/fhcho
483 aerosol(nahso5) = 0.0
484 aerosol(nahmsa) = 0.0
485 end if
486 if ( (idecomp_hmsa_hso5 .gt. 0) .and. &
487 (iabs(iradical-101) .le. 1) ) then
488 gas(ngso2) = gas(ngso2) + (yaq(19)*wso2/wmol(25))/fso2
489 aerosol(na4) = aerosol(na4) + (yaq(18)*wmol(2)/wmol(24))
490 end if
491
492 !
493 ! save yaq to yaq_end
494 !
495 do i = 1, 13
496 yaq_end(i) = yaq(i)
497 end do
498
499 !
500 ! istat_fatal = fatal error status
501 ! (ipar(3) .ne. 0) --> svode had problems --> 10s digit = -1
502 ! (ipar(7) .ne. 0) --> fullequil had problems --> 100s digit = -2
503 ! istat_warn = warning status
504 ! (ipar(5) .ne. 0) --> hybrd had problems --> 10s digit = +1
505 ! (negconc_count .ne. 0) --> 100s digit = +2
506 !
507 istat_fatal = 0
508 if (ipar(3) .ne. 0) istat_fatal = istat_fatal - 10
509 if (ipar(7) .ne. 0) istat_fatal = istat_fatal - 200
510
511 istat_warn = 0
512 if (ipar(5) .ne. 0) istat_warn = istat_warn + 10
513 if (negconc_count .ne. 0) istat_warn = istat_warn + 200
514
515 !
516 ! the code neglects the change in the size distribution shape
517 ! of the nonvolatile aerosol components because of the sulfate
518 ! production. the change in the sulfate distribution is calculated
519 ! using the distribution functions while the change in the
520 ! distribution of the volatile inorganic aerosol components will
521 ! calculated by the aerosol module after cloud evaporation
522 !
523 !
524 !.....end of code
525 return
526 end subroutine aqoperator1
527
528
529
530 !----------------------------------------------------------------------
531 !
532 ! this routine prepares the necessary variables for the call
533 ! to the stiff ode integrator (svode)
534 ! its interface to the main code is the same as with the
535 ! assymptotic solver so one can interchange solvers easily
536 !
537 ! last modification: june 13, 1998 by s.p.
538 !
539 !----------------------------------------------------------------------
540 subroutine aqintegr1( neqa, y, stime, stout, rpar, ipar )
541
542 use module_data_cmu_bulkaqchem, only: &
543 iopt, itask, itol, liw1, lrw1, &
544 mdiag_svode, meqn1, meqn1max, mf, &
545 tola, tolr, worki, workr
546
547 use module_cmu_dvode_solver, only: dvode
548
549 ! arguments
550 integer neqa, ipar(*)
551 double precision y(meqn1max), stime, stout, rpar(*)
552
553 ! local variables
554 integer i, istate, liw, lrw
555 integer iwork(30+meqn1max)
556 double precision rwork(22+9*meqn1max+2*meqn1max**2)
557 double precision atol(meqn1max),rtol(meqn1max)
558
559 !
560 do i=1, meqn1
561 atol(i) = tola ! absolute tolerance in ug/m3
562 rtol(i) = tolr ! relative tolerance
563 if (i .gt. 13) atol(i) = 1.0e-20
564 enddo
565 lrw = lrw1 ! dimension of double precision work vector
566 liw = liw1 ! dimension of integer work vector
567 iwork(6) = worki ! steps allowed
568 rwork(6) = workr ! maximum step in seconds
569
570 !
571 ! ready for the call to svode
572 !
573 istate = 1
574
575 call dvode( fex1, neqa, y, stime, stout, itol, rtol, atol, itask, &
576 istate, iopt, rwork, lrw, iwork, liw, jex, mf, &
577 rpar, ipar )
578
579 if (istate .ne. 2) then
580 ! write(*,*) '*** aqintegr1 -- svode failed'
581 ! write(*,*) ' stime, istate =', stime, istate
582 ! stop
583
584 if (mdiag_svode .gt. 0) then
585 write(6,*) &
586 '*** module_cmuaq_bulk, aubr aqintegr1 -- ' // &
587 'svode failed, istate =', istate
588 end if
589 ipar(3) = ipar(3) + 1
590 ipar(4) = istate
591 endif
592 !
593 do i=1, meqn1
594 if (y(i) .lt. 0.0) y(i) = 1.e-20
595 enddo
596 !
597 return
598 end subroutine aqintegr1
599
600
601
602 !----------------------------------------------------------------------
603 subroutine fex1( neq, t, y, f, rpar, ipar )
604
605 use module_data_cmu_bulkaqchem, only: meqn1, meqn1max
606
607 ! arguments
608 integer neq, ipar(*)
609 double precision t, y(meqn1max), f(meqn1max), rpar(*)
610
611 ! local variables
612 double precision a(meqn1max),b(meqn1max)
613
614
615 call aqratesa( meqn1, t, y, a, b, f, rpar, ipar )
616
617 return
618 end subroutine fex1
619
620
621
622 !----------------------------------------------------------------------
623 subroutine jex( neq, t, y, ml, mu, pd, nrowpd, rpar, ipar )
624 integer neq, ml, mu, nrowpd, ipar(*)
625 double precision t, y(neq), pd(nrowpd,neq), rpar(*)
626 call peg_error_fatal( -1, &
627 '*** module_cmuaq_bulk, subr jex -- should not be here ***' )
628 end subroutine jex
629
630
631
632 !----------------------------------------------------------------------
633 !
634 ! last modification : feb. 15, 1995 by s. pandis
635 !*************************************************************************
636 ! aqrates.for
637 !*************************************************************************
638 ! this routine calculates the rates of change of the y's at time tmin
639 ! for the aqueous-phase species
640 ! it does bulk-phase chemistry
641 !
642 !
643 subroutine aqratesa( neq, t, yaq, aqprod, aqdest, yaqprime, &
644 rpar, ipar )
645
646 use module_data_cmu_bulkaqchem, only: &
647 akeq, akhen, akre, avdiam, &
648 caratio, chyd, cmet, co2_mixrat, con, &
649 epsfcn, factor, firon, fman, gcon, gmol, &
650 iradical, ldfjac, lr, &
651 maxfev, meqn1, meqn1max, ml, mode, mu, &
652 mdiag_hybrd, mdiag_rsrate, &
653 ngas, ngch3o2, nghno2, ngho2, &
654 ngno, ngno2, ngno3, ngo3, ngoh, ngpan, &
655 nprint, numfunc, &
656 pres_cmuaq_cur, &
657 rad, rideal, &
658 temp_cmuaq_cur, &
659 wh2o2, whcho, whcooh, wnh3, wso2, wmol, wvol, &
660 xtol
661
662 ! arguments
663 integer neq, ipar(*)
664 double precision t, yaq(meqn1max), yaqprime(meqn1max)
665 double precision aqprod(meqn1max), aqdest(meqn1max)
666 double precision rpar(*)
667
668 ! local variables
669 integer nfev, info, i, j, n
670 integer icount, iprint
671
672 double precision spres(22)
673 double precision c(46)
674 double precision fgl(21), flg(21), gfgl(21), gflg(21), rp(28), rl(28)
675 double precision dp(49), dl(49), rr(120)
676 double precision x(numfunc)
677 double precision fvec(numfunc),diag(numfunc),fjac(ldfjac,numfunc),r(lr)
678 double precision qtf(numfunc),wa1(numfunc),wa2(numfunc),wa3(numfunc)
679 double precision wa4(numfunc)
680
681 double precision, parameter :: one=1.
682 !
683 double precision hc1, hc2, wl, wlm, tlwc, vlwc, hyd, arytm, hno2, af
684 double precision ah, chen
685 double precision rsrate, sulfrate, sulfrateb
686 double precision form
687 double precision ph, pres
688 double precision qsat
689 double precision tmin, temp
690 double precision watvap, watcont, crustal, salt, sodium
691 ! lwc,wat. vapor in g/m3
692 ! crustal species concentration inside droplets (ug/m3)
693 ! na concentration inside droplets (ug/m3)
694 double precision dum
695 double precision gascon(ngas)
696
697
698 !
699 ! unload ipar, rpar values
700 !
701 icount = ipar(1)
702 iprint = ipar(2)
703 temp = rpar(1)
704 pres = rpar(2)
705 watcont = rpar(3)
706 watvap = rpar(4)
707 crustal = rpar(5)
708 salt = rpar(6)
709 sodium = rpar(7)
710 ph = rpar(8)
711 do i = 1, ngas
712 gascon(i) = rpar(8+i)
713 end do
714
715 temp_cmuaq_cur = temp
716 pres_cmuaq_cur = pres
717
718 !
719 ! calculation of the temperature for this time (temp)
720 ! (assuming linear temperature change for each operator stepp)
721 !
722 tmin = t !in seconds
723 n = numfunc
724
725 call qsaturation (temp, qsat) ! qsat is in ug/m3
726 !
727 ! zero the rate of change vectors
728 !
729 do i=1, meqn1
730 yaqprime(i)=0.0
731 aqprod(i)=0.0
732 aqdest(i)=0.0
733 enddo
734 !
735 ! set dummy ph to zero for printing only
736 !
737 ph=0.0
738 !
739 ! find total lwc
740 !
741 tlwc = watcont*1.e6 ! in ug/m3
742 vlwc=tlwc*1.e-12 ! vol/vol
743 !
744 ! check for negative values
745 !
746 !sp do i=1, meqn1
747 !sp if (yaq(i) .le. 0.) yaq(i)=1.e-20
748 !sp enddo
749 !
750 ! reconstruct the matrices using the available yaq values
751 !
752 ! ** gas phase concentrations (in ppm) **
753 ! (some are assumed to remain constant during the aqueous-phase
754 ! chemical processes, the rest are updated)
755 ! note : all hcho and hcooh are still in the gas-phase
756 !
757 spres(1) = yaq(3)*rideal*temp/(1000.*wso2*pres) ! so2 (g)
758 spres(2) = 1.e-20 ! h2so4 (g)
759 spres(3) = gascon(nghno2) ! hno2 (g)
760 spres(4) = 1.e-20 ! hno3 (g) [it has already dissolved]
761 ! spres(5) = 350. ! co2 (g)
762 spres(5) = co2_mixrat ! 2004-nov-15 rce
763 spres(6) = yaq(4)*rideal*temp/(1000.*wh2o2*pres) ! h2o2 (g)
764 hc1=rideal*temp*vlwc*akhen(7)
765 hc2=rideal*temp/(1000.*whcho*pres)
766 spres(7) = yaq(1)*hc2/(hc1+1.) ! hcho (g)
767 spres(8) = yaq(2)*rideal*temp/(1000.*whcooh*pres) ! hcooh (g)
768 spres(9) = gascon(ngno) ! no (g)
769 spres(10) = gascon(ngno2) ! no2 (g)
770 spres(11) = gascon(ngo3) ! o3 (g)
771 spres(12) = gascon(ngpan) ! pan (g)
772 spres(13) = 1.e-20 ! ch3coooh (g)
773 spres(14) = 1.e-20 ! ch3ooh (g)
774 spres(15) = 1.e-20 ! hcl (g) [it has already dissolved]
775 spres(16) = gascon(ngoh) ! oh (g)
776 spres(17) = gascon(ngho2) ! ho2 (g)
777 spres(18) = gascon(ngno3) ! no3 (g)
778 spres(19) = yaq(5)*rideal*temp/(1000.*wnh3*pres) ! nh3 (g)
779 spres(20) = gascon(ngch3o2) ! ch3o2 (g)
780 spres(21) = 1.e-20 ! ch3oh (g)
781 spres(22) = watvap*rideal*temp/(1000.*18.*pres) ! h2o (g)
782 !
783 ! calculation of gcon vector (gas-phase concentrations in mole/l)
784 !
785 do 5 i=1,22
786 ! gcon(i) = spres(i)*1.e-6/(0.08206*temp)
787 gcon(i) = spres(i)*1.e-6/(rideal*temp)
788 5 continue
789 !
790 ! ** radius and lwc for the section
791 !
792 rad = 0.5*avdiam ! in m
793 wl = watcont ! in g/m3
794 wvol= wl*1.e-6 ! in vol/vol
795 wlm = wl*1.e6 ! in ug/m3
796 !
797 ! loading of all metal species
798 ! note : na+ is an input. the rest of the metal species are
799 ! calculated as a fraction of the crustal aerosol mass.
800 !
801 cmet(1) = firon*crustal*1000./(55.85*wlm) ! fe(3+) mol/l
802 cmet(2) = fman*crustal*1000./(54.9*wlm) ! mn(2+) mol/l
803 cmet(3) = salt*1000./(23.*wlm) ! na(+) mol/l
804 cmet(4) = caratio*crustal*1000./(40.08*wlm) ! ca(2+) mol/l
805 !
806 ! do not let the fe(3+) and mn(2+) concentrations exceed a
807 ! certain limit because they cause extreme stiffness due to
808 ! the reaction s(iv)->s(vi)
809 !
810 ! if (cmet(1) .gt. 1.0e-5) cmet(1)=1.0e-5
811 ! if (cmet(2) .gt. 1.0e-5) cmet(2)=1.0e-5
812 !
813 ! loading of the main aqueous concentrations (in m)
814 !
815 con(1) = yaq(6)*1000./(wmol(1)*wlm) ! s(iv)
816 if (con(1) .lt. 1.e-20) con(1)=1.e-20
817 con(2) = yaq(11)*1000./(wmol(2)*wlm) ! s(vi)
818 con(3) = 0. ! n(iii) (determined later)
819 con(4) = yaq(8)*1000./(wmol(4)*wlm) ! n(v)
820 con(5) = 0. ! co2 (determined later)
821 con(6) = yaq(7)*1000./(wmol(6)*wlm) ! h2o2
822 if (con(6) .lt. 1.e-20) con(6)=1.e-20
823 con(7) = akhen(7)*spres(7)*1.e-6 ! hcho
824 con(8) = 0. ! hcooh (determined later)
825 con(9) = 1.0e-6*akhen(9)*spres(9) ! no
826 con(10) = 1.0e-6*akhen(10)*spres(10) ! no2
827 con(11) = 1.0e-6*akhen(11)*spres(11) ! o3
828 con(12) = 1.0e-6*akhen(12)*spres(12) ! pan
829 con(13) = 1.0e-6*akhen(13)*spres(13) ! ch3coooh
830 con(14) = 1.0e-6*akhen(14)*spres(14) ! ch3ooh
831 con(15) = yaq(9)*1000./(wmol(15)*wlm) ! hcl
832 con(16) = 0. ! oh (determined later)
833 con(17) = 0. ! ho2 (determined later)
834 con(18) = 0. ! no3 (determined later)
835 con(19) = yaq(10)*1000./(wmol(19)*wlm) ! nh3
836 con(20) = 1.0e-6*akhen(20)*spres(20) ! ch3o2
837 con(21) = 1.0e-6*akhen(21)*spres(21) ! ch3oh
838 con(22) = 0. ! cl (determined later)
839 con(23) = 0. ! cloh- (determined later)
840 con(24) = 0. ! so4- (determined later)
841 con(25) = 0. ! so5- (determined later)
842 con(26) = yaq(12)*1000./(wmol(26)*wlm) ! hso5-
843 con(27) = yaq(13)*1000./(wmol(27)*wlm) ! hoch2so3-
844 con(28) = 0. ! co3- (determined later)
845 !
846 ! set a minimum concentration (to avoid divisions by zero)
847 !
848 do i=1, 28
849 if (con(i) .lt. 1.0e-20) con(i)=1.0e-20
850 enddo
851
852 !
853 ! 27-oct-2005 rce - previously there was a bug here and the code
854 ! would continue on even if fullequil failed
855 !
856 ! calculation of ph and volatile concentrations (co2, n(iii), hcooh)
857 ! (solve the system of equations)
858 ! if ipar(7)>0, this means that fullequil has failed
859 ! and it's time to shut down the integration
860 ! do this by setting the yaqprime=0, which hopefully will allow
861 ! the integrator to complete
862 !
863 if (ipar(7) .le. 0) &
864 call fullequil(con,spres,cmet,akeq,akhen,vlwc,temp,hyd,ipar)
865 if (ipar(7) .gt. 0) then
866 do i = 1, meqn1
867 yaqprime(i) = 0.0
868 end do
869 ph=30.0
870 return
871 end if
872 ph=-log10(hyd)
873
874 ! when iradical = 100/101/102, the radical species are computed
875 ! by directly by dvode rather than by hybrd
876 ! the con's for radicals are loaded here (after call to fullequil)
877 ! as that more closely follows the approach with hybrd
878 if (iabs(iradical-101) .le. 1) then
879 con(16) = yaq(14)*1000./(wmol(16)*wlm) ! oh (mw = 17)
880 con(17) = yaq(15)*1000./(wmol(17)*wlm) ! ho2 (mw = 33)
881 con(18) = yaq(16)*1000./(wmol(18)*wlm) ! no3 (mw = 62)
882 con(23) = yaq(17)*1000./(wmol(23)*wlm) ! cloh- (mw = 52.5)
883 con(24) = yaq(18)*1000./(wmol(24)*wlm) ! so4- (mw = 96)
884 con(25) = yaq(19)*1000./(wmol(25)*wlm) ! so5- (mw = 122)
885 con(28) = yaq(20)*1000./(wmol(28)*wlm) ! co3- (mw = 60)
886 dum = 1.0d-35
887 con(16) = max( con(16), dum )
888 con(17) = max( con(17), dum )
889 con(18) = max( con(18), dum )
890 con(23) = max( con(23), dum )
891 con(24) = max( con(24), dum )
892 con(25) = max( con(25), dum )
893 con(28) = max( con(28), dum )
894 end if
895
896 !
897 ah = rideal*temp*vlwc*akhen(3)*(1.+akeq(7)/hyd)/pres
898 hno2=spres(3)/(1.+ah)
899 con(3)=akhen(3)*(1.+akeq(7)/hyd)*1.e-6*hno2
900 !
901 chen=akhen(5)*(1.+akeq(8)/hyd+akeq(8)*akeq(9)/hyd**2)
902 con(5)=chen*spres(5)*1.e-6 ! [co2 t]aq m
903 !
904 af=rideal*temp*vlwc*akhen(8)*(1.+akeq(13)/hyd)/pres
905 form=spres(8)/(1.+af) ! new hcooh(g) ppm
906 con(8)=akhen(8)*(1.+akeq(13)/hyd)*1.e-6*form
907 !
908 ! we calculate the ionic species concentrations
909 !
910 call values(hyd, con, cmet, akeq, c)
911 !
912 ! bypass the radical calculation by hybrd if necessary
913 !
914 if (iradical .eq. 0) go to 270
915 if (iabs(iradical-101) .le. 1) go to 280
916 !
917 ! this is where the call to hybrd was made when
918 ! the aqueous radical species were treated as steady state
919 ! now we treate iradical.gt.0 same as iradical=100
920 !
921 go to 280
922 !
923 ! set the radical concentrations to zero
924 !
925 270 continue
926 ! if (info .eq. 2 .or. info .eq. 3 &
927 ! .or. info .eq. 4 .or. info .eq. 5 .or. iradical .eq. 0) then
928 con(16)=1.e-25
929 con(17)=1.e-25
930 con(18)=1.e-25
931 con(23)=1.e-25
932 con(24)=1.e-25
933 con(25)=1.e-25
934 con(28)=1.e-25
935 ! endif
936 !
937 ! pseudo-steady-state approximation for cl radical
938 ! not used because it is of secondary importance
939 !sp call values(hyd, con, cmet, akeq,c)
940 !sp call react(c,cmet,con,akre,rr,arytm)
941 !sp pro=rr(23)+rr(49)+rr(96)
942 !sp if (con(22) .le. 0.0)then
943 !sp con(22) = 1.e-20
944 !sp go to 20
945 !sp endif
946 !sp destr=(rr(24)+rr(25)+rr(26)+rr(27)+rr(28)+rr(29)+rr(30)+rr(42)+
947 !sp & rr(56)+rr(61)+rr(69)+rr(109))/con(22)
948 !sp if (destr .eq. 0.0)then
949 !sp con(22)= 1.e-10
950 !sp go to 20
951 !sp endif
952 !sp con(22)=pro/destr
953 !sp 20 continue
954 !
955
956 280 continue
957 call values(hyd, con, cmet, akeq,c)
958 !
959 ! final calculation of reaction rates
960 !
961 call react(c,cmet,con,akre,rr,arytm)
962 !
963 ! calculation of net production and consumption rates
964 !
965 call addit(rr, arytm, rp, rl)
966 !
967 ! calculation of mass transfer rates
968 !
969 call mass(wvol,rad,temp,gcon,con,c,akeq,akhen,fgl,flg,gfgl,gflg)
970 !
971 ! calculation of net rates of change
972 !
973 call differ(rp,rl,fgl,flg,gfgl,gflg,dp,dl)
974 !
975 ! calculation of right hand sides of the derivatives
976 !
977 ! ** gas-phase species **
978 ! (rates in ug/m3 air s)
979 yaqprime(1) = 1.e9*wvol*wmol(7)*(rp(7)-rl(7)) ! hcho t
980 aqprod(1) = 1.e9*wvol*wmol(7)*rp(7)
981 aqdest(1) = 1.e9*wvol*wmol(7)*rl(7)
982 !
983 yaqprime(2) = 1.e9*wvol*wmol(8)*(rp(8)-rl(8)) ! hcooh t
984 aqprod(2) = 1.e9*wvol*wmol(8)*rp(8)
985 aqdest(2) = 1.e9*wvol*wmol(8)*rl(8)
986 !
987 yaqprime(3) = 1.e9*gmol(1)*(dp(29)-dl(29)) ! so2(g)
988 aqprod(3) = 1.e9*gmol(1)*dp(29)
989 aqdest(3) = 1.e9*gmol(1)*dl(29)
990 !
991 yaqprime(4) = 1.e9*gmol(6)*(dp(34)-dl(34)) ! h2o2(g)
992 aqprod(4) = 1.e9*gmol(6)*dp(34)
993 aqdest(4) = 1.e9*gmol(6)*dl(34)
994 !
995 yaqprime(5) = 1.e9*gmol(19)*(dp(47)-dl(47)) ! nh3(g)
996 aqprod(5) = 1.e9*gmol(19)*dp(47)
997 aqdest(5) = 1.e9*gmol(19)*dl(47)
998 !
999 ! ** aqueous-phase species **
1000 !
1001 yaqprime(6)= 1.e9*wvol*wmol(1)*(dp(1)-dl(1)) ! s(iv)
1002 aqprod(6) = 1.e9*wvol*wmol(1)*dp(1)
1003 aqdest(6) = 1.e9*wvol*wmol(1)*dl(1)
1004 !
1005 yaqprime(7)= 1.e9*wvol*wmol(6)*(dp(6)-dl(6)) ! h2o2
1006 aqprod(7) = 1.e9*wvol*wmol(6)*dp(6)
1007 aqdest(7) = 1.e9*wvol*wmol(6)*dl(6)
1008 !
1009 yaqprime(8)= 1.e9*wvol*wmol(4)*(dp(4)-dl(4)) ! n(v)
1010 aqprod(8) = 1.e9*wvol*wmol(4)*dp(4)
1011 aqdest(8) = 1.e9*wvol*wmol(4)*dl(4)
1012 !
1013 yaqprime(9)= 1.e9*wvol*wmol(15)*(dp(15)-dl(15)) ! cl-
1014 aqprod(9) = 1.e9*wvol*wmol(15)*dp(15)
1015 aqdest(9) = 1.e9*wvol*wmol(15)*dl(15)
1016 !
1017 yaqprime(10)= 1.e9*wvol*wmol(19)*(dp(19)-dl(19)) ! nh4+
1018 aqprod(10) = 1.e9*wvol*wmol(19)*dp(19)
1019 aqdest(10) = 1.e9*wvol*wmol(19)*dl(19)
1020 !
1021 yaqprime(11)= 1.e9*wvol*wmol(2)*(dp(2)-dl(2)) ! s(vi)
1022 aqprod(11) = 1.e9*wvol*wmol(2)*dp(2)
1023 aqdest(11) = 1.e9*wvol*wmol(2)*dl(2)
1024 !
1025 yaqprime(12)= 1.e9*wvol*wmol(26)*(dp(26)-dl(26)) ! hso5-
1026 aqprod(12) = 1.e9*wvol*wmol(26)*dp(26)
1027 aqdest(12) = 1.e9*wvol*wmol(26)*dl(26)
1028 !
1029 yaqprime(13)= 1.e9*wvol*wmol(27)*(dp(27)-dl(27)) ! hmsa
1030 aqprod(13) = 1.e9*wvol*wmol(27)*dp(27)
1031 aqdest(13) = 1.e9*wvol*wmol(27)*dl(27)
1032 !
1033 if (iabs(iradical-101) .le. 1) then
1034 !
1035 yaqprime(14)= 1.e9*wvol*wmol(16)*(dp(16)-dl(16)) ! oh(aq)
1036 aqprod(14) = 1.e9*wvol*wmol(16)*dp(16)
1037 aqdest(14) = 1.e9*wvol*wmol(16)*dl(16)
1038 !
1039 yaqprime(15)= 1.e9*wvol*wmol(17)*(dp(17)-dl(17)) ! ho2(aq)
1040 aqprod(15) = 1.e9*wvol*wmol(17)*dp(17)
1041 aqdest(15) = 1.e9*wvol*wmol(17)*dl(17)
1042 !
1043 yaqprime(16)= 1.e9*wvol*wmol(18)*(dp(18)-dl(18)) ! no3(aq)
1044 aqprod(16) = 1.e9*wvol*wmol(18)*dp(18)
1045 aqdest(16) = 1.e9*wvol*wmol(18)*dl(18)
1046 !
1047 yaqprime(17)= 1.e9*wvol*wmol(23)*(dp(23)-dl(23)) ! cloh-(aq)
1048 aqprod(17) = 1.e9*wvol*wmol(23)*dp(23)
1049 aqdest(17) = 1.e9*wvol*wmol(23)*dl(23)
1050 !
1051 yaqprime(18)= 1.e9*wvol*wmol(24)*(dp(24)-dl(24)) ! so4-(aq)
1052 aqprod(18) = 1.e9*wvol*wmol(24)*dp(24)
1053 aqdest(18) = 1.e9*wvol*wmol(24)*dl(24)
1054 !
1055 yaqprime(19)= 1.e9*wvol*wmol(25)*(dp(25)-dl(25)) ! so5-(aq)
1056 aqprod(19) = 1.e9*wvol*wmol(25)*dp(25)
1057 aqdest(19) = 1.e9*wvol*wmol(25)*dl(25)
1058 !
1059 yaqprime(20)= 1.e9*wvol*wmol(28)*(dp(28)-dl(28)) ! co3-(aq)
1060 aqprod(20) = 1.e9*wvol*wmol(28)*dp(28)
1061 aqdest(20) = 1.e9*wvol*wmol(28)*dl(28)
1062 end if
1063 !
1064 ! calculation of appropriate destruction rate
1065 !
1066 do 50 i=1, meqn1
1067 if (yaq(i) .le. 1.e-20) then
1068 aqdest(i) = 0.0
1069 go to 50
1070 endif
1071 aqdest(i) = aqdest(i)/yaq(i)
1072 50 continue
1073 !
1074 ! change to avoid divisions by zero in integration
1075 !
1076 do i=1,meqn1
1077 if (aqdest(i) .le. 1.e-18) aqdest(i)=1.e-18
1078 enddo
1079
1080 !
1081 ! calculation of characteristic times (used for debugging)
1082 !
1083 !db tsm=100.
1084 !db do 110 i=1, meqn1
1085 !
1086 !db if (aqdest(i) .le. 1.e-10) go to 110
1087 !db tchar=1./aqdest(i)
1088 !
1089 !db if (tchar .lt. 0.01)then
1090 !db write(80,*)tmin,i,yaq(i),yaqprime(i),tchar
1091 !db write(6,*) i,yaq(i),yprime(i)
1092 !db endif
1093 !
1094 !db if (tchar .lt. tsm) then
1095 !db tsm=tchar
1096 !db endif
1097 !
1098 !db110 continue
1099
1100 !
1101 ! mass balance for sulfur
1102 ! original sulfrate calc includes so2(g), so2(aq), so4=, hso5-, hmsa
1103 ! and does not always balance
1104 ! sulfrateb calc also includes so4-, so5- and gives a closer balance
1105 !
1106 ! sulfrate=yaqprime(3)/gmol(1)+yaqprime(6)/wmol(1)+ &
1107 ! yaqprime(11)/wmol(2)
1108 ! rsrate=sulfrate/(abs(yaqprime(3))+abs(yaqprime(11)) + &
1109 ! abs(yaqprime(6)))
1110 ! if (abs(rsrate) .ge. 0.01) then
1111 ! write(80, *)'problem at ',tmin/60.
1112 ! write(80, *) yaqprime(3),yaqprime(6),yaqprime(11)
1113 ! write(80, *) rsrate
1114 ! write(80, *)'************************'
1115 ! endif
1116 sulfrate = (yaqprime( 3)/gmol( 1)) + (yaqprime( 6)/wmol( 1)) + &
1117 (yaqprime(11)/wmol( 2)) + (yaqprime(12)/wmol(26)) + &
1118 (yaqprime(13)/wmol(27))
1119 sulfrateb = sulfrate + &
1120 1.0e9*wvol*(rp(24) - rl(24) + rp(25) - rl(25))
1121 rsrate = abs(yaqprime( 3)/gmol( 1)) + abs(yaqprime( 6)/wmol( 1)) + &
1122 abs(yaqprime(11)/wmol( 2)) + abs(yaqprime(12)/wmol(26)) + &
1123 abs(yaqprime(13)/wmol(27))
1124 rsrate = max(rsrate, 1.0d-30)
1125 if (mdiag_rsrate .gt. 0) then
1126 if (abs(sulfrateb/rsrate) .ge. 1.0e-5) then
1127 write(6,*)
1128 write(6,'(a,1p,3e11.2)') &
1129 'aqratesa sulfbal prob - rerr, rerrb, t =', &
1130 (sulfrate/rsrate), (sulfrateb/rsrate), tmin
1131 write(6,'(a,1p,e15.6/4e15.6)') &
1132 'yaqprime/wmol so2,siv,svi,hso5-,hmsa =', &
1133 (yaqprime(3)/gmol(1)), (yaqprime(6)/wmol(1)), &
1134 (yaqprime(11)/wmol(2)), (yaqprime(12)/wmol(26)), &
1135 (yaqprime(13)/wmol(27))
1136 write(6,*)
1137 end if
1138 end if
1139
1140 !
1141 ! diagnostic output
1142 !
1143 icount=icount+1
1144 if (icount .ge. iprint)then
1145 !rce write(6,120)tmin/60.,yaq(11) !,ph,rsrate,x(1)*1.e12,yaq(13)
1146 !rce write(79,*)tmin/60.,ph
1147 ! printing of all reaction rates for debugging
1148 !sp write(3,*)tmin/60.,'****(um/hr)*******'
1149 !sp do i=1,109
1150 !sp write(3,*)i,rr(i)*1.e6*3600.
1151 !sp enddo
1152 icount=0
1153 endif
1154 !120 format(1x,2(1x,f8.4))
1155 120 format( 'aqratesa - tmin, yaq(11)=so4', 2(1x,f8.4) )
1156
1157 !
1158 ! load ipar, rpar values
1159 !
1160 ipar(1) = icount
1161 rpar(8) = ph
1162
1163 ! write(*,'(a,1p,8e10.2 )') 'xxx t,yaq1-7 ', t, (yaq(i), i=1,7)
1164 ! write(*,'(a,1p,8e10.2 )') 'xxx t,yaq8-13', t, (yaq(i), i=8,13)
1165 ! write(*,'(a,1p,8e10.2 )') 'xxx t,rad-con', t, &
1166 ! (con(i), i=16,18), (con(i), i=23,25), (con(i), i=28,28)
1167 ! dum = 1.e9*wvol
1168 ! write(*,'(a,1p,8e10.2/)') 'xxx t,rad-yaq', t, &
1169 ! (con(i)*wmol(i)*dum, i=16,18), (con(i)*wmol(i)*dum, i=23,25), &
1170 ! (con(i)*wmol(i)*dum, i=28,28)
1171
1172 return
1173 end subroutine aqratesa
1174
1175
1176
1177
1178 !************************************************************************
1179 ! this routine calculates the the steady state species concentrations
1180 !************************************************************************
1181 subroutine steady(radius,temp,c,con,gcon,akeq,akhen,akre)
1182 !
1183 !..inputs:
1184 ! radius : droplet radius in m
1185 ! temp : temperature (in k)
1186 ! c(46) : the concentrations of the rest of the aqueous-phase species
1187 ! gcon(22) : gas-phase concentrations
1188 ! akeq,akhen,akre : reaction constants
1189 !..outputs:
1190 ! x(8) the values of the steady state species concentrations
1191 !
1192
1193 ! arguments
1194 double precision radius, temp
1195 double precision c(46),gcon(22),akeq(17),akhen(21),akre(120)
1196 double precision con(28)
1197
1198 ! local variables
1199 integer icount
1200 double precision a1, a2, a3, a4, acc, airl
1201 double precision dg, ho2, o2
1202 double precision kn,n,ikn,kmt
1203 double precision rideal
1204 double precision x(8)
1205 !
1206 ! airl is the mean free path of air. later we have to substitute
1207 ! the numerical value given here by a function of temperature
1208 ! airl=65x10-9 m
1209 airl=65.e-9
1210 ! kn is the knudsen number
1211 kn=airl/radius
1212 ikn=1.0/kn
1213 ! acc is the accomodation coefficient assumed the same here for
1214 ! all the species
1215 acc=0.01
1216 ! n is the coefficient entering the flux expression
1217 n=1.0/(1.+((1.33+0.71*ikn)/(1.+ikn)+4.*(1.-acc) &
1218 /(3.*acc))*kn)
1219 ! dg is the gas phase diffusivity assumed here the same for all
1220 ! the gases. dg=1.x10-5 m**2/sec
1221 dg=1.0e-5
1222 ! rideal is the gas constant in [atm/K/(mol/liter)]
1223 rideal=0.082058
1224 kmt=(3.0*n*dg)/(radius*radius)
1225 !
1226 ! iteration loop
1227 !
1228 do icount=1,2
1229
1230 !
1231 ! no3 calculation
1232 !
1233 x(3)=(kmt*gcon(18))/(akre(43)*c(8)+akre(45)+akre(46)*c(29)+ &
1234 akre(47)*c(30) +akre(48)*c(14)+ &
1235 akre(49)*c(27)+akre(54)*c(18)+akre(59)*c(19)+akre(71)*c(35)+ &
1236 akre(109)*c(2)+kmt/(akhen(18)*rideal*temp))
1237 con(18)=x(3)
1238 !
1239 ! so4- calculation
1240 !
1241 x(5)=(akre(109)*c(2)*x(3)+2.*akre(86)*c(40)*c(40)) &
1242 /(akre(89)*c(41)+akre(92)*c(2)+ &
1243 akre(93)*c(3)+akre(94)*c(29)+akre(95)*c(30)+ &
1244 akre(96)*c(45)+akre(97)*c(14)+akre(98)*c(8)+ &
1245 akre(99)*c(12)+akre(100)*c(19)+akre(101)*c(27)+ &
1246 akre(102)*c(18)+akre(108)*c(35))
1247 c(39)=x(5)
1248 con(24)=c(39)
1249 !
1250 a1=c(46)/(akeq(15)+c(46))
1251 a2=akeq(15)/(akeq(15)+c(46))
1252 !
1253 ! ho2 calculation
1254 !
1255 x(2)=((akre(48)*c(14)+akre(54)*c(18)+akre(59)*c(19)+ &
1256 akre(71)*c(35)) * x(3)+ &
1257 (akre(97)*c(14)+akre(100)*c(19)+akre(102)*c(18)+ &
1258 akre(108)*c(35))*x(5)+ &
1259 2.0*akre(14)*c(45)*c(22)+ &
1260 akre(28)*c(14)*c(36)+akre(29)*c(14)*c(37)+akre(55)*c(18)*c(22)+ &
1261 akre(56)*c(18)*c(36)+akre(61)*c(19)*c(36)+akre(69)*c(35)*c(36)+ &
1262 akre(65)*c(25)+akre(15)*c(22)*c(15)+akre(58)*c(19)*c(22)+ &
1263 kmt*gcon(17) +akre(5)*c(14)*c(28)+akre(11)*c(22)*c(28)+ &
1264 akre(20)*c(14)*c(44)+akre(50)*c(17)*c(28)+akre(52)*c(18)*c(28)+ &
1265 akre(57)*c(19)*c(28)+akre(60)*c(19)*c(44)+akre(67)*c(35)*c(28)+ &
1266 akre(68)*c(35)*c(44)+akre(84)*c(18)*c(40)+ &
1267 akre(85)*c(19)*c(40))/ &
1268 (a1*(akre(3)*c(28)+2.*akre(6)*c(29)+2.*akre(7)*c(30)+ &
1269 akre(9)*c(14)+akre(12)*c(22)+akre(25)*c(36)+ &
1270 akre(27)*c(37)+akre(46)*x(3)+akre(63)*c(34)+ &
1271 akre(94)*c(39)+akre(107)*c(3))+ &
1272 a2*(akre(4)*c(28)+2.*akre(8)*c(30)+akre(10)*c(14)+ &
1273 akre(13)*c(22)+akre(18)*c(12)+akre(19)*c(44)+akre(26)*c(36)+ &
1274 akre(47)*x(3)+akre(64)*c(34)+akre(83)*c(40)+akre(95)*c(39)) &
1275 +(kmt*a1)/(akhen(17)*rideal*temp))
1276 !
1277 ho2=(x(2)*c(46))/(akeq(15)+c(46))
1278 o2=(x(2)*akeq(15))/(akeq(15)+c(46))
1279 c(29)=ho2
1280 c(30)=o2
1281 con(17)=c(29)+c(30)
1282 !
1283 a3=(akre(21)*akre(22)*c(27))/(akre(22)+akre(23)*c(46))
1284 a4=(akre(22)*akre(24)*c(37))/(akre(22)+akre(23)*c(46))
1285 !
1286 ! oh calculation
1287 !
1288 x(1)=(2.*akre(1)*c(14)+akre(15)*c(15)*c(22)+akre(30)*c(45)* &
1289 c(36)+ &
1290 akre(35)*c(7)+akre(36)*c(8)+akre(44)*c(10)+akre(55)*c(18)*c(22)+ &
1291 akre(58)*c(19)*c(22)+akre(65)*c(25)+kmt*gcon(16)+a4+ &
1292 (akre(9)*c(14)+akre(12)*c(22)+akre(107)*c(3))*ho2+ &
1293 (akre(10)*c(14)+akre(13)*c(22))*o2+akre(96)*c(45)*x(5))/ &
1294 (akre(3)*ho2+akre(4)*o2+akre(5)*c(14)+akre(11)*c(22)+ &
1295 akre(17)*c(12)+akre(21)*c(27)+akre(33)*c(20)+akre(34)*c(21)+ &
1296 akre(37)*c(7)+akre(38)*c(8)+akre(50)*c(17)+akre(52)*c(18)+ &
1297 akre(57)*c(19)+akre(66)*c(25)+akre(67)*c(35)+akre(80)*c(3)+ &
1298 akre(81)*c(2)+akre(88)*c(41)+akre(115)*c(42)+ &
1299 kmt/(akhen(16)*rideal*temp)-a3)
1300 c(28)=x(1)
1301 con(16)=c(28)
1302 !
1303 ! cloh- calculation
1304 !
1305 x(4)=(akre(21)*c(27)*x(1)+akre(24)*c(37))/(akre(22)+ &
1306 akre(23)*c(46))
1307 c(38)=x(4)
1308 con(23)=c(38)
1309 !
1310 ! co3- calculation
1311 !
1312 x(7)=(akre(17)*c(12)*x(1)+akre(99)*c(12)*x(5)+akre(18)*c(12)*o2)/ &
1313 (akre(19)*o2+akre(20)*c(14)+akre(41)*c(8)+akre(60)*c(19)+ &
1314 akre(68)*c(35))
1315 c(44)=x(7)
1316 con(28)=c(44)
1317 !
1318 ! so5- calculation
1319 !
1320 x(6)=(akre(116)*c(2)*c(36)+(akre(80)*c(3)+akre(81)*c(2)+ &
1321 akre(88)*c(41)+akre(115)*c(42))*x(1)+(akre(89)*c(41)+ &
1322 akre(92)*c(2)+akre(93)*c(3))*x(5))/ &
1323 (akre(83)*o2+akre(84)*c(18)+akre(85)*c(19)+2.*akre(86)*c(40))
1324 c(40)=x(6)
1325 con(25)=c(40)
1326
1327 enddo
1328 !
1329
1330
1331 return
1332 end subroutine steady
1333
1334
1335
1336
1337 ! ***********************************************************************
1338 ! the routine fullequil solves the electroneutrality equation
1339 ! using bisection method when the concentrations of [c], n(iii)
1340 ! and hcooh are unknown.
1341 ! inputs in the sub are con(28),spres(21),cmet(4),akeq(17),akhen(21).
1342 ! output is the value of [h+]
1343 ! the routine electro gives values of the electroneutrality equation.
1344 ! inputs in the sub are x=[h+],con(28),cmet(4),akeq(17).
1345 ! output is the value of f ( f(x)=0.0 at the solution)
1346 ! ***********************************************************************
1347
1348 subroutine fullequil(acon,aspres,acmet,aakeq,aakhen,awv, &
1349 atemp,axsol,ipar)
1350
1351 use module_data_cmu_bulkaqchem, only: mdiag_fullequil
1352
1353 ! arguments
1354 integer ipar(*)
1355 double precision acon(28), aspres(21), acmet(4), aakeq(17), aakhen(21)
1356 double precision awv, atemp, axsol
1357
1358 ! local variables
1359 integer i, k, ntry
1360 integer ipass_01, idum_01
1361 double precision aa, bb, error, f, fa, fm, rtol, x, xm
1362 double precision wv, temp, xsol
1363 double precision con(28), spres(21), cmet(4), akeq(17), akhen(21)
1364 !
1365 ! change variables to double precision to avoid low ph errors
1366 !
1367 ipass_01 = 1
1368 idum_01 = 0
1369 300 continue
1370
1371 do k=1,28
1372 con(k)=acon(k)
1373 enddo
1374 !
1375 do k=1,21
1376 spres(k)=aspres(k)
1377 akhen(k)=aakhen(k)
1378 enddo
1379 !
1380 do k=1,4
1381 cmet(k)=acmet(k)
1382 enddo
1383 !
1384 do k=1,17
1385 akeq(k)=aakeq(k)
1386 enddo
1387 !
1388 wv=awv
1389 temp=atemp
1390 !
1391 ! we find the initial interval [aa,bb] for the bisection method
1392 ! new version (31/10/87)
1393 x=10.0d0**(-14)
1394 call electro(x,fa,con,spres,cmet,akeq,akhen,wv,temp)
1395 aa=x
1396
1397 ! do 1035 i=-14,1
1398 do 1035 i = -(14+idum_01), (1+idum_01)
1399 x=10.0d0**i
1400 call electro(x,f,con,spres,cmet,akeq,akhen,wv,temp)
1401 if (f*fa .ge. 0.0d0) then
1402 aa=x
1403 fa=f
1404 else
1405 bb=x
1406 ! fb=f
1407 go to 1040
1408 end if
1409 1035 continue
1410
1411 ! 27-oct-2005 rce - previously there was a bug here and the code
1412 ! continued on to label 1040 after reporting the "mistake in fullequil"
1413 ! Now the code tries a greater range of initial ph values,
1414 ! then gives up if it fails.
1415 !
1416 ! unable to find 2 initial hion values that bracket the "solution"
1417 ! if ipass_01 = 1, try again
1418 if (ipass_01 .eq. 1) then
1419 write(6,*) &
1420 '*** module_cmuaq_bulk - mistake in fullequil with ipass_01 = 1'
1421 ipass_01 = 2
1422 idum_01 = 5
1423 goto 300
1424 else if (ipass_01 .eq. 2) then
1425 write(6,*) &
1426 '*** module_cmuaq_bulk - mistake in fullequil with ipass_01 = 2'
1427 ipass_01 = 3
1428 idum_01 = 10
1429 goto 300
1430 end if
1431
1432 ! otherwise, report the error and exit
1433 ipar(7) = ipar(7) + 1
1434 if (mdiag_fullequil .gt. 0) then
1435 write(6,*) &
1436 '*** module_cmuaq_bulk - mistake in fullequil - con, cmet ='
1437 write(6,*) con, cmet
1438 return
1439 end if
1440
1441 1040 continue
1442 !
1443 ! bisection method for the solution of the equation
1444 ! rtol : relative tolerance
1445 ntry=0
1446 ! 02-nov-2005 rce - smaller rtol for h+ makes dvode run faster
1447 ! rtol=0.00001d0
1448 rtol=1.0d-8
1449 1050 error= dabs(bb-aa)/aa
1450 ntry=ntry+1
1451 if (error .le. rtol) then
1452 xsol=(aa+bb)/2.0d0
1453 axsol=xsol ! single precision
1454 return
1455 end if
1456 xm=(aa+bb)/2.0d0
1457 call electro(xm,fm,con,spres,cmet,akeq,akhen,wv,temp)
1458 if (fa*fm .gt. 0.0d0) then
1459 aa=xm
1460 fa=fm
1461 else
1462 bb=xm
1463 ! fb=fm
1464 end if
1465 go to 1050
1466 end subroutine fullequil
1467
1468
1469
1470
1471 ! ***********************************************************************
1472 ! routine that gives values of the electroneutrality equation
1473 ! called by fullequil
1474 ! ***********************************************************************
1475
1476 subroutine electro(x,f,con,spres,cmet,zkeq,zkhen,wv,temp)
1477
1478 use module_data_cmu_bulkaqchem, only: &
1479 mprescribe_ph, rideal, xprescribe_ph
1480
1481 ! arguments
1482 !
1483 ! original subr arguments were akeq & akhen
1484 ! renamed them to zkeq & zkhen
1485 ! to avoid conflict with module_data_cmu_bulkaqchem
1486 !
1487 double precision x, f, wv, temp
1488 double precision con(28),spres(21),cmet(4),zkeq(17),zkhen(21)
1489
1490 ! local variables
1491 double precision bparam, cparam, cl, dfac, dform, diak
1492 double precision f1, f2, f3, f4, f5, form, hcl, hno2
1493 double precision cc(46)
1494 !
1495 cc(2)=(zkeq(1)*con(1)*x)/(x*x+zkeq(1)*x+zkeq(1)*zkeq(2)) ! hso3-
1496 cc(3)=(zkeq(1)*zkeq(2)*con(1))/(x*x+zkeq(1)*x+zkeq(1)*zkeq(2)) ! so3--
1497 cc(5)=(zkeq(3)*con(2)*x)/(x*x+zkeq(3)*x+zkeq(3)*zkeq(4)) ! hso4-
1498 cc(6)=(zkeq(3)*zkeq(4)*con(2))/(x*x+zkeq(3)*x+zkeq(3)*zkeq(4)) ! so4--
1499 !
1500 ! ** no2- calculation from equilibrium **
1501 dfac=rideal*temp*wv*zkhen(3)*(1.+zkeq(7)/x)
1502 hno2 = spres(3)/(1.+dfac) ! new hno2(g) in ppm
1503 cc(8)=zkhen(3)*1.e-6*(zkeq(7)/x)*hno2
1504 !
1505 cc(10)=(zkeq(6)*con(4))/(x+zkeq(6)) ! no3-
1506 !
1507 ! ** co2 equilibrium (constant gas co2 concentration) **
1508 cc(12)= zkeq(8)*zkhen(5)*spres(5)*1.e-6/x
1509 cc(13)= zkeq(9)*cc(12)/x
1510 !
1511 cc(15)=(zkeq(5)*con(6))/(x+zkeq(5)) ! ho2-
1512 !
1513 ! ** hcoo- equilibrium (partitioning with gas-phase) **
1514 dform=rideal*temp*wv*zkhen(8)*(1.+zkeq(13)/x)
1515 form=spres(8)/(1.+dform) ! new hcooh
1516 cc(19)=zkhen(8)*1.e-6*(zkeq(13)/x)*form
1517 !
1518 cc(30)=(zkeq(15)*con(17))/(x+zkeq(15)) ! o2-
1519 cc(38)=con(23) ! cloh-
1520 cc(39)=con(24) ! so4-
1521 cc(40)=con(25) ! so5-
1522 cc(41)=con(26) ! hso5-
1523 cc(42)=(x*con(27))/(x+zkeq(17)) ! hoch2so3-
1524 cc(43)=(zkeq(17)*con(27))/(x+zkeq(17)) ! -och2so3-
1525 cc(44)=con(28) ! co3-
1526 cc(45)=zkeq(11)/x ! oh-
1527
1528 bparam=zkeq(16)+con(15)-con(22)
1529 cparam=zkeq(16)*con(22)
1530 diak=bparam*bparam+4.0*cparam
1531 if (diak .le. 0.) diak=1.0e-20
1532 cl=(-bparam+(diak)**0.5)/2.0
1533 hcl=(x*(con(15)-con(22)+cl))/(x+zkeq(14))
1534 cc(27)=(zkeq(14)*hcl)/x ! cl-
1535 cc(36)=(zkeq(14)*cl*hcl)/(zkeq(16)*x) ! cl2-
1536
1537 cc(33)=(zkeq(10)*x*con(19))/(zkeq(11)+zkeq(10)*x) ! nh4+
1538 cc(46)=x ! h+
1539
1540 f1=cc(2)+2.0*cc(3)+cc(5)+2.0*cc(6)+cc(8)+cc(10)
1541 f2=cc(12)+2.0*cc(13)+cc(15)+cc(19)+cc(27)+cc(30)
1542 f3=cc(36)+cc(38)+cc(39)+cc(40)+cc(41)+cc(42)
1543 f4=2.0*cc(43)+cc(44)+cc(45)-cc(33)-cc(46)
1544 f5=-3.0*cmet(1)-2.0*cmet(2)-cmet(3)-2.0*cmet(4)
1545
1546 f=f1+f2+f3+f4+f5
1547
1548 if (mprescribe_ph .gt. 0) then
1549 f = 10.0**(-xprescribe_ph) - x
1550 end if
1551
1552 return
1553 end subroutine electro
1554
1555
1556
1557 !----------------------------------------------------------------------
1558 !
1559 ! routines used by the aqeous-phase module
1560 !
1561 ! 1. dropinit : initialization
1562
1563 subroutine dropinit
1564
1565 use module_data_cmu_bulkaqchem, only: &
1566 amol, gmol, &
1567 wmol, wh2o2, wh2so4, whcho, whcl, whcooh, whno3, wnh3, wso2
1568
1569
1570 ! local variables
1571
1572
1573 !
1574 ! loading of molecular weights
1575 !
1576 wso2 = 64.
1577 wh2o2 = 34.
1578 whcho = 30.
1579 whcooh = 46.
1580 wnh3 = 17.
1581 whno3 = 63.
1582 whcl = 36.5
1583 wh2so4 = 98.
1584 !
1585 ! molecular weights
1586 !
1587 wmol(1)= 81.0e0
1588 wmol(2)= 96.0e0
1589 wmol(3)= 47.0e0
1590 wmol(4)= 62.0e0
1591 wmol(5)= 62.0e0
1592 wmol(6)= 34.0e0
1593 wmol(7)= 48.0e0 ! was previously 60.0e0
1594 wmol(8)= 46.0e0
1595 wmol(9)= 30.0e0
1596 wmol(10)=46.0e0
1597 wmol(11)=48.0e0
1598 wmol(12)=121.0e0
1599 wmol(13)=76.0e0
1600 wmol(14)=48.0e0
1601 wmol(15)=35.5e0
1602 wmol(16)=17.0e0
1603 wmol(17)=33.0e0
1604 wmol(18)=62.0e0
1605 wmol(19)=18.0e0
1606 wmol(20)=47.0e0
1607 wmol(21)=32.0e0
1608 wmol(22)=35.5e0
1609 wmol(23)=52.50e0
1610 wmol(24)=96.0e0
1611 wmol(25)=112.0e0
1612 wmol(26)=113.0e0
1613 wmol(27)=111.0e0
1614 wmol(28)=60.00e0
1615 wmol(29)=18.0e0
1616 !
1617 amol(1)= 55.85e0
1618 amol(2)= 55.0e0
1619 amol(3)= 23.0e0
1620 !
1621 gmol(1)=64.0
1622 gmol(2)=98.08
1623 gmol(3)=47.02
1624 gmol(4)=63.02
1625 gmol(5)=44.01
1626 gmol(6)=34.02
1627 ! 09-nov-2005 rce - set gmol(6) == wh2o2 to conserve h2o2
1628 gmol(6)=34.0
1629 gmol(7)=30.03
1630 gmol(8)=46.00
1631 gmol(9)=30.01
1632 gmol(10)=46.01
1633 gmol(11)=48.00
1634 gmol(12)=121.05
1635 gmol(13)=76.00
1636 gmol(14)=48.00
1637 gmol(15)=36.50
1638 gmol(16)=17.00
1639 gmol(17)=33.01
1640 gmol(18)=62.01
1641 gmol(19)=17.00
1642 gmol(20)=47.00
1643 gmol(21)=32.00
1644 gmol(22)=18.00
1645
1646 return
1647 end subroutine dropinit
1648
1649
1650
1651 !----------------------------------------------------------------------
1652 subroutine qsaturation(temp,qsat)
1653 !
1654 ! this routine calculates the saturation mass concentration (in ug/m3)
1655 ! over liquid water for a temperature temp (k)
1656 !
1657
1658 ! arguments
1659 double precision temp,qsat
1660
1661 ! local variables
1662 double precision psat, t ! these should be double precision ?
1663 double precision rideal,a0,a1,a2,a3,a4,a5,a6,esat,csat
1664 !
1665 t = temp-273.15 ! in c
1666 rideal = 0.082058d0 ! gas constant in [atm/K/(mol/liter)]
1667 a0 = 6.107799961d-0
1668 a1 = 4.436518521d-1
1669 a2 = 1.428945805d-2
1670 a3 = 2.650648471d-4
1671 a4 = 3.031240396d-6
1672 a5 = 2.034080948d-8
1673 a6 = 6.136820929d-11
1674 !
1675 esat=a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+a6*t))))) ! in mb
1676 psat = esat/(1000.0*1.01325) ! in atm
1677 csat = psat/(rideal*temp) ! in mole/l
1678 qsat = 18000.0d0*csat*1.e6 ! in ug/m3
1679 ! write(6,*)t,esat/1000.,qsat
1680 return
1681 end subroutine qsaturation
1682
1683
1684
1685
1686 ! s u b r o u t i n e s f o r d r o p l e t p r o g r a m
1687 !************************************************************************
1688 ! this routine calculates the 20 equilibrium costants {akeq}, the 20
1689 ! henry's law costants {akhen} and the 120 reaction rate costants {akre}
1690 ! with only input the temperature (temp).
1691 ! included in the routine are the corresponding enthalpies
1692 ! {dheq,dhhen,dhre} and the costants at 298 k {bkeq,bkhen,bkre}.
1693 !************************************************************************
1694
1695 subroutine constants(temp)
1696
1697 use module_data_cmu_bulkaqchem, only: akeq, akhen, akre, &
1698 maqurxn_all, maqurxn_sulf1, mopt_eqrt_cons
1699
1700 ! arguments
1701 double precision temp
1702
1703 ! local variable
1704 integer i, iusei
1705 ! dimension akeq(17),akhen(21),akre(120)
1706 double precision, save :: dheq(17),dhhen(21),dhre(120)
1707 double precision, save :: bkeq(17),bkhen(21),bkre(120)
1708
1709 data dheq/1960.,1500.,0.,2720.,-3730.,8700.,-1260., &
1710 -1000.,-1760., &
1711 -450.,-6710.,4020.,-20.,6900.,0.e0,0.,0./
1712 data bkeq/1.23e-2,6.61e-8,1.e3,1.02e-2,2.2e-12,15.4,5.1e-4, &
1713 4.46e-7,4.68e-11,1.75e-5,1.0e-14,1.82e3,1.78e-4,1.74e6,3.5e-5, &
1714 5.26e-6,2.0e-12/
1715 data dhhen/3120.,0.0,4780.,8700.,2420.,6620., &
1716 6460.e0,5740.e0,1480.e0, &
1717 2500.e0,2300.e0,5910.e0,6170.e0,5610.e0,2020.e0,5280.e0, &
1718 6640.e0,8700.e0,3400.e0, &
1719 5600.e0,4900.e0/
1720 data bkhen/1.23e0,0.0e0,49.e0,2.1e5,3.4e-2,7.45e4,6.3e3, &
1721 3.5e3,1.9e-3, &
1722 0.01e0,1.13e-2,2.9e0,473.e0,227.e0,727.e0,25.e0,2000.e0,2.1e5, &
1723 75.e0,6.e0,220.e0/
1724
1725
1726 data dhre/0.0e0,0.0e0,-1500.e0,-1500.e0,-1700.e0,-2365.e0, &
1727 -1500.e0,0.0e0,0.0e0, &
1728 0.0e0,0.0e0,0.0e0,-1500.e0,0.0e0,-2520.e0,0.0e0,-1910.e0, &
1729 0.0e0,-1500.e0,-2820.e0, &
1730 -1500.e0,0.0e0,0.0e0,0.0e0,-1500.e0,-1500.e0,-1500.e0, &
1731 -3370.e0,0.0e0,-2160.e0, &
1732 -1500.e0,-1500.e0,-1500.e0,-1500.e0,0.0e0,0.0e0,-1500.e0, &
1733 -1500.e0,-6693.e0,-6950.e0, &
1734 0.0e0,-1500.e0,-1500.e0,0.0e0,0.0e0,-1500.e0,-1500.e0, &
1735 -2800.e0,-1500.e0,-1500.e0, &
1736 0.0e0,-1500.e0,-5180.e0,-3200.e0,0.0e0,-4300.e0,-1500.e0, &
1737 0.0e0,-1500.e0,-3400.e0, &
1738 -2600.e0,0.0e0,-3000.e0,-1600.e0,0.0e0,-1700.e0,-1500.e0, &
1739 -4500.e0,-4400.e0,-1800.e0, &
1740 -2800.e0,0.0e0,-5530.e0,-5280.e0,-4430.e0,-13700.e0, &
1741 -11000.e0,-13700e0,-11000.e0, &
1742 -1500.e0,-1500.e0,-3100.e0,-1500.e0,-5300.e0,-4000.e0, &
1743 -1500.e0,-4755.e0,-1900.e0, &
1744 0.0e0,-6650.e0,-7050.e0,-1500.e0,-1500.e0,-1500.e0, &
1745 -1500.e0,-1500.e0,-2000.e0, &
1746 -1500.e0,-2100.e0,-1500.e0,-1500.e0,-2700.e0,0.0e0, &
1747 -3800.e0,-4000.e0,0.0e0,0.0e0, &
1748 -1800.e0,0.0e0,0.0e0,0.0e0,-6100.e0,-4900.e0, &
1749 -4500.e0,-1500.e0,-1500.e0, &
1750 -2000.e0,0.e0,-1800.e0,120.e0/
1751
1752
1753 data bkre/2.5e-6,2.0e-5,7.0e9,1.0e10,2.7e7, &
1754 8.6e5,1.0e8,0.3e0,0.5e0, &
1755 0.13e0,2.0e9,1.0e4,1.5e9,70.e0,2.8e6,7.8e-3,1.5e7, &
1756 1.5e6,4.0e8,8.0e5, &
1757 4.3e9,6.1e9,2.1e10,1.3e3,4.5e9,1.0e9,3.1e9, &
1758 1.4e5,4.5e7,7.3e6, &
1759 2.0e8,1.0e8,2.0e10,1.3e9,3.7e-5,6.3e-6,1.0e9, &
1760 1.0e10,6.3e3,5.0e5, &
1761 4.0e5,2.5e8,1.2e9,1.0e-7,1.0e-5,4.5e9,1.0e9, &
1762 1.0e6,1.0e8,2.0e9, &
1763 0.1e0,1.6e8,4.6e-6,2.1e5,5.0e0,6.7e3,2.5e9,100.0e0, &
1764 6.0e7,1.1e5,1.9e6, &
1765 4.0e-4,7.6e5,5.0e7,5.4e-7,2.7e7,4.5e8,2.6e3, &
1766 3.5e3,1.9e7,1.0e6, &
1767 2.4e4,3.7e5,1.5e9,1.3e6,4.7e0,0.82e0,5.0e3,1.0e7, &
1768 4.6e9,4.2e9,3.0e5, &
1769 1.0e8,200.e0,1.4e4,2.e8,7.5e7,1.7e7,1.e5,0.31e0, &
1770 1.8e-3,1.3e9,5.3e8, &
1771 5.0e9,5.0e9,8.0e7,1.2e7,8.8e8,9.1e6,1.7e8, &
1772 2.0e8,1.4e6,6.7e-3, &
1773 1.9e7,5.0e7,6.0e2,1.0e6,2.5e7,1.0e8,2.0e6, &
1774 1.42e2,4.77e3,2.94e2, &
1775 3.6e3,1.4e9,3.4e8, &
1776 2.5e4,1.0e5,2.5e7,120.e0/
1777
1778
1779 ! when mopt_eqrt_cons=20, set s(iv)+h2o2 rxn rate constant
1780 ! to that used in mtcrm and testaqu22
1781 if (mopt_eqrt_cons .eq. 20) then
1782 bkre(75) = 4.19e7
1783 dhre(75) = -1950.0
1784 end if
1785
1786 do 1020 i=1,17
1787 akeq(i)=bkeq(i)*exp(dheq(i)*(1.0/temp-1.0/298.0))
1788 1020 continue
1789 do 1025 i=1,21
1790 akhen(i)=bkhen(i)*exp(dhhen(i)*(1.0/temp-1.0/298.0))
1791 1025 continue
1792 do 1030 i=1,120
1793 akre(i)=bkre(i)*exp(dhre(i)*(1.0/temp-1.0/298.0))
1794 1030 continue
1795
1796 ! turn reactions on/off selectively
1797 do i = 1, 120
1798 iusei = 0
1799 if (maqurxn_all .gt. 0) iusei = 1
1800 if (maqurxn_sulf1 .gt. 0) then
1801 if ((i .ge. 72) .and. (i .le. 75)) iusei = 1
1802 end if
1803 if (iusei .le. 0) akre(i) = 0.0
1804 if (iusei .le. 0) bkre(i) = 0.0
1805 end do
1806
1807 return
1808 end subroutine constants
1809
1810
1811
1812
1813 !*************************************************************************
1814 ! this routine calculates the values of the concentrations
1815 ! of all the 46 species that appear in our aqueous phase mechanism.
1816 ! it has as inputs the values of [h+],con(28),cmet(3),akeq(17) and as
1817 ! outputs the values of cc(46)
1818 !*************************************************************************
1819
1820 subroutine values(x,con,cmet,akeq,cc)
1821
1822 ! arguments
1823 double precision x
1824 double precision con(28),cmet(4),akeq(17),cc(46)
1825
1826 ! local variables
1827 double precision bparam,cparam,diak,cl,hcl
1828 !
1829 ! species in the aqueous phase mechanism
1830 ! cc (1 - 46)
1831 ! 1.) so2*h2o 24.) ch3c(o)ooh
1832 ! 2.) hso3(-) 25.) ch3ooh
1833 ! 3.) so3(2-) 26.) hcl
1834 ! 4.) h2so4 27.) cl(-)
1835 ! 5.) hso4(-) 28.) oh
1836 ! 6.) so4(2-) 29.) ho2
1837 ! 7.) hno2 30.) o2(-)
1838 ! 8.) no2(-) 31.) no3
1839 ! 9.) hno3 32.) nh4oh
1840 ! 10.) no3(-) 33.) nh4(+)
1841 ! 11.) co2*h2o 34.) ch3o2
1842 ! 12.) hco3(-) 35.) ch3oh
1843 ! 13.) co3(2-) 36.) cl2(-)
1844 ! 14.) h2o2 37.) cl
1845 ! 15.) ho2(-) 38.) cloh(-)
1846 ! 16.) hcho 39.) so4(-)
1847 ! 17.) h2c(oh)2 40.) so5(-)
1848 ! 18.) hcooh 41.) hso5(-)
1849 ! 19.) hcoo(-) 42.) hoch2so3(-)
1850 ! 20.) no 43.) och2so3(2-)
1851 ! 21.) no2 44.) co3(-)
1852 ! 22.) o3 45.) oh(-)
1853 ! 23.) pan 46.) h(+)
1854 !
1855 ! con(1-28)
1856 !
1857 ! 1.) so2(g) 15.) hcl(g)
1858 ! 2.) h2so4(g) 16.) oh(g)
1859 ! 3.) hno2(g) 17.) ho2(g)
1860 ! 4.) hno3(g) 18.) no3(g)
1861 ! 5.) co2(g) 19.) nh3(g)
1862 ! 6.) h2o2(g) 20.) ch3o2(g)
1863 ! 7.) hcho(g) 21.) ch3oh(g)
1864 ! 8.) hcooh(g) 22.) cl2(-), cl
1865 ! 9.) no(g) 23.) cloh(-)
1866 ! 10.) no2(g) 24.) so4(-)
1867 ! 11.) o3(g) 25.) so5(-)
1868 ! 12.) pan(g) 26.) hso5(-)
1869 ! 13.) ch3c(o)ooh(g) 27.) hoch2so3(-),och2so3(2-)
1870 ! 14.) ch3ooh(g) 28.) co3(-)
1871 !
1872 bparam=akeq(16)+con(15)-con(22)
1873 cparam=akeq(16)*con(22)
1874 diak=bparam*bparam+4.0d0*cparam
1875 if (diak .le. 0.0d0) diak=1.0d-30
1876 cl=(-bparam+(diak)**0.5d0)/2.0d0
1877 hcl=(x*(con(15)-con(22)+cl))/(x+akeq(14))
1878
1879 cc(1)=(con(1)*x*x)/(x*x+akeq(1)*x+akeq(1)*akeq(2))
1880 cc(2)=(akeq(1)*con(1)*x)/(x*x+akeq(1)*x+akeq(1)*akeq(2))
1881 cc(3)=(akeq(1)*akeq(2)*con(1))/(x*x+akeq(1)*x+akeq(1)*akeq(2))
1882
1883 cc(4)=(con(2)*x*x)/(x*x+akeq(3)*x+akeq(3)*akeq(4))
1884 cc(5)=(akeq(3)*con(2)*x)/(x*x+akeq(3)*x+akeq(3)*akeq(4))
1885 cc(6)=(akeq(3)*akeq(4)*con(2))/(x*x+akeq(3)*x+akeq(3)*akeq(4))
1886
1887 cc(7)=(x*con(3))/(x+akeq(7))
1888 cc(8)=(akeq(7)*con(3))/(x+akeq(7))
1889
1890 cc(9)=(x*con(4))/(x+akeq(6))
1891 cc(10)=(akeq(6)*con(4))/(x+akeq(6))
1892
1893 cc(11)=(x*x*con(5))/(x*x+akeq(8)*x+akeq(8)*akeq(9))
1894 cc(12)=(akeq(8)*con(5)*x)/(x*x+akeq(8)*x+akeq(8)*akeq(9))
1895 cc(13)=(akeq(8)*akeq(9)*con(5))/(x*x+akeq(8)*x+akeq(8)*akeq(9))
1896
1897 cc(14)=(x*con(6))/(x+akeq(5))
1898 cc(15)=(akeq(5)*con(6))/(x+akeq(5))
1899
1900 cc(16)=con(7)/(1.0d0+akeq(12))
1901 cc(17)=(akeq(12)*con(7))/(1.0d0+akeq(12))
1902
1903 cc(18)=(x*con(8))/(x+akeq(13))
1904 cc(19)=(akeq(13)*con(8))/(x+akeq(13))
1905
1906 cc(20)=con(9)
1907
1908 cc(21)=con(10)
1909
1910 cc(22)=con(11)
1911
1912 cc(23)=con(12)
1913
1914 cc(24)=con(13)
1915
1916 cc(25)=con(14)
1917
1918 cc(26)=hcl
1919 cc(27)=(akeq(14)*hcl)/x
1920
1921 cc(28)=con(16)
1922
1923 cc(29)=(x*con(17))/(x+akeq(15))
1924 cc(30)=(akeq(15)*con(17))/(x+akeq(15))
1925
1926 cc(31)=con(18)
1927
1928 cc(32)=(akeq(11)*con(19))/(akeq(11)+akeq(10)*x)
1929 cc(33)=(akeq(10)*x*con(19))/(akeq(11)+akeq(10)*x)
1930
1931 cc(34)=con(20)
1932
1933 cc(35)=con(21)
1934
1935 cc(36)=(akeq(14)*cl*hcl)/(akeq(16)*x)
1936 cc(37)=cl
1937
1938 cc(38)=con(23)
1939 cc(39)=con(24)
1940 cc(40)=con(25)
1941 cc(41)=con(26)
1942 cc(42)=(x*con(27))/(x+akeq(17))
1943 cc(43)=(akeq(17)*con(27))/(x+akeq(17))
1944 cc(44)=con(28)
1945 cc(45)=akeq(11)/x
1946 cc(46)=x
1947
1948 return
1949 end subroutine values
1950
1951
1952
1953
1954 !************************************************************************
1955 ! this program contains the routine react which calculates
1956 ! the rates of all the reactions taking place in the aqueous phase.
1957 ! inputs in the sub are the 46 concentrations ,the 3 metal concentrations
1958 ! the 28 main species concentrations and the 120 reaction constants.
1959 ! output is the matrix of the 120 reaction rates.
1960 !************************************************************************
1961
1962 subroutine react(c,cmet,con,zkre,rr,arytm)
1963
1964 use module_data_cmu_bulkaqchem, only: chlorine, kiron, photo, &
1965 mopt_eqrt_cons
1966
1967 ! arguments
1968 !
1969 ! original argument was akre
1970 ! renamed it to zkre
1971 ! to avoid conflict with module_data_cmu_bulkaqchem
1972 !
1973 ! dimension c(46),cmet(4),con(28),zkre(120),rr(120)
1974 double precision c(46),cmet(4),con(28),zkre(120),rr(120)
1975 double precision arytm
1976
1977 ! local variables
1978 double precision ph, r1, r2, r3, r4, r5, sn
1979
1980
1981 rr(1)=zkre(1)*c(14)*photo
1982 rr(2)=zkre(2)*c(22)*photo
1983 rr(3)=zkre(3)*c(28)*c(29)
1984 rr(4)=zkre(4)*c(28)*c(30)
1985 rr(5)=zkre(5)*c(28)*c(14)
1986 rr(6)=zkre(6)*c(29)*c(29)
1987 rr(7)=zkre(7)*c(29)*c(30)
1988 rr(8)=zkre(8)*c(30)*c(30)
1989 rr(9)=zkre(9)*c(29)*c(14)
1990 rr(10)=zkre(10)*c(30)*c(14)
1991 rr(11)=zkre(11)*c(28)*c(22)
1992 rr(12)=zkre(12)*c(29)*c(22)
1993 rr(13)=zkre(13)*c(30)*c(22)
1994 rr(14)=zkre(14)*c(45)*c(22)
1995 rr(15)=zkre(15)*c(15)*c(22)
1996 if (c(22) .le. 0.0d0) c(22)=1.0d-30
1997 rr(16)=zkre(16)*c(14)*(c(22)**0.5)
1998
1999 rr(17)=zkre(17)*c(12)*c(28)
2000 rr(18)=zkre(18)*c(12)*c(30)
2001 rr(19)=zkre(19)*c(44)*c(30)
2002 rr(20)=zkre(20)*c(44)*c(14)
2003 rr(21)=zkre(21)*c(27)*c(28)*chlorine
2004 rr(22)=zkre(22)*c(38)*chlorine
2005 rr(23)=zkre(23)*c(46)*c(38)*chlorine
2006 rr(24)=zkre(24)*c(37)*chlorine
2007 rr(25)=zkre(25)*c(29)*c(36)*chlorine
2008 rr(26)=zkre(26)*c(30)*c(36)*chlorine
2009 rr(27)=zkre(27)*c(29)*c(37)*chlorine
2010 rr(28)=zkre(28)*c(14)*c(36)*chlorine
2011 rr(29)=zkre(29)*c(37)*c(14)*chlorine
2012 rr(30)=zkre(30)*c(45)*c(36)*chlorine
2013
2014 rr(31)=zkre(31)*c(20)*c(21)
2015 rr(32)=zkre(32)*c(21)*c(21)
2016 rr(33)=zkre(33)*c(20)*c(28)
2017 rr(34)=zkre(34)*c(21)*c(28)
2018 rr(35)=zkre(35)*c(7)*photo
2019 rr(36)=zkre(36)*c(8)*photo
2020 rr(37)=zkre(37)*c(7)*c(28)
2021 rr(38)=zkre(38)*c(8)*c(28)
2022 rr(39)=zkre(39)*c(46)*c(14)*c(7)
2023 rr(40)=zkre(40)*c(8)*c(22)
2024 rr(41)=zkre(41)*c(8)*c(44)
2025 rr(42)=zkre(42)*c(8)*c(36)*chlorine
2026 rr(43)=zkre(43)*c(8)*c(31)
2027 rr(44)=zkre(44)*c(10)*photo
2028 rr(45)=zkre(45)*c(31)*photo
2029 rr(46)=zkre(46)*c(31)*c(29)
2030 rr(47)=zkre(47)*c(31)*c(30)
2031 rr(48)=zkre(48)*c(31)*c(14)
2032 rr(49)=zkre(49)*c(31)*c(27)*chlorine
2033 rr(50)=zkre(50)*c(17)*c(28)
2034 rr(51)=zkre(51)*c(17)*c(22)
2035 rr(52)=zkre(52)*c(18)*c(28)
2036 rr(53)=zkre(53)*c(18)*c(14)
2037 rr(54)=zkre(54)*c(18)*c(31)
2038 rr(55)=zkre(55)*c(18)*c(22)
2039 rr(56)=zkre(56)*c(18)*c(36)*chlorine
2040 rr(57)=zkre(57)*c(19)*c(28)
2041 rr(58)=zkre(58)*c(19)*c(22)
2042 rr(59)=zkre(59)*c(19)*c(31)
2043 rr(60)=zkre(60)*c(19)*c(44)
2044 rr(61)=zkre(61)*c(19)*c(36)*chlorine
2045 rr(62)=zkre(62)*c(23)
2046 rr(63)=zkre(63)*c(34)*c(29)
2047 rr(64)=zkre(64)*c(34)*c(30)
2048 rr(65)=zkre(65)*c(25)*photo
2049 rr(66)=zkre(66)*c(25)*c(28)
2050 rr(67)=zkre(67)*c(35)*c(28)
2051 rr(68)=zkre(68)*c(35)*c(44)
2052 rr(69)=zkre(69)*c(35)*c(36)*chlorine
2053 rr(70)=zkre(70)*c(25)*c(28)
2054 rr(71)=zkre(71)*c(35)*c(31)
2055
2056 rr(72)=(zkre(72)*c(1)+zkre(73)*c(2)+zkre(74)*c(3))*c(22)
2057 rr(73)=(zkre(75)*c(14)*c(1))/(1.0d0+16.0d0*c(46))
2058 ! when mopt_eqrt_cons=20, calc s(iv)+h2o2 rxn rate
2059 ! same as in mtcrm and testaqu22
2060 if (mopt_eqrt_cons .eq. 20) then
2061 rr(73)=(zkre(75)*c(14)*c(2)*c(46))/(1.0d0+16.0d0*c(46))
2062 end if
2063 !
2064 ! rate expressions for the metal catalysed oxidation of s(iv)
2065 !
2066 ! ** phenomenological expression by martin et al. (1991) **
2067 ph=-log10(c(46))
2068 if (kiron .eq. 1) then
2069 !
2070 if (ph .le. 3.0) rr(74)=6.*cmet(1)*con(1)/c(46)
2071 if (ph .gt. 3.0 .and. ph .le. 4.5) &
2072 rr(74) = 1.e9*con(1)*cmet(1)*cmet(1)
2073 if (ph .gt. 4.5 .and. ph .le. 6.5) rr(74) = 1.0e-3*con(1)
2074 if (ph .gt. 6.5) rr(74)=1.0e-4*con(1)
2075 endif
2076
2077 ! ** expression by martin (1984) **
2078 if (kiron .eq. 2) then
2079 if ((c(46) .ge. 1.0d-4).and.(con(1) .ge. 1.0d-5)) then
2080 r1=(zkre(76)*cmet(2)*cmet(2))/c(46)
2081 r2=(zkre(77)*cmet(1)*con(1)/c(46))
2082 if (cmet(2) .le. 0.0d0) cmet(2)=1.0d-30
2083 r3=r2*(1.0d0+(1.7d3*cmet(2)**1.5)/(6.3d-6+cmet(1)))
2084 rr(74)=r1+r3
2085 go to 1300
2086 end if
2087
2088 if (cmet(1)*cmet(2) .lt. 1.0d-15) then
2089 sn=1.0d0
2090 else
2091 sn=3.0d0
2092 end if
2093
2094 if ((c(46) .ge. 1.0d-4).and.(con(1) .lt. 1.0d-5)) then
2095 rr(74)=sn*(zkre(78)*cmet(2)*c(2)+zkre(77)*cmet(1)*con(1)/c(46))
2096 go to 1300
2097 end if
2098
2099 if ((c(46) .lt. 1.0d-4).and.(con(1) .ge. 1.0d-5)) then
2100 r4=zkre(76)*cmet(2)*cmet(2)/c(46)
2101 r5=zkre(79)*cmet(1)*con(1)*con(1)
2102 rr(74)=r4+r5
2103 go to 1300
2104 end if
2105
2106 rr(74)=zkre(78)*cmet(2)*c(2)
2107
2108 1300 continue
2109 endif
2110
2111 ! 09-nov-2005 rce - if rate constants 76,77,78,79 are all zero, set rr(74)=0.
2112 ! This allows ANY/ALL rxns to be turned on/off in subr constants.
2113 if (abs(zkre(76)+zkre(77)+zkre(78)+zkre(79)) .le. 1.0e-37) rr(74)=0.0
2114
2115 ! ** end of martin's expression **
2116
2117 rr(75)=zkre(80)*c(3)*c(28)
2118 rr(76)=zkre(81)*c(2)*c(28)
2119 rr(77)=zkre(82)*c(40)*c(2)+zkre(117)*c(40)*c(3)
2120 rr(78)=zkre(83)*c(40)*c(30)
2121 rr(79)=zkre(84)*c(40)*c(18)
2122 rr(80)=zkre(85)*c(40)*c(19)
2123 rr(81)=zkre(86)*c(40)*c(40)
2124 rr(82)=zkre(87)*c(41)*c(2)*c(46)
2125 rr(83)=zkre(88)*c(41)*c(28)
2126 rr(84)=zkre(89)*c(41)*c(39)
2127 rr(85)=zkre(90)*c(41)*c(8)
2128 rr(86)=zkre(91)*c(41)*c(27)
2129 rr(87)=zkre(92)*c(39)*c(2)
2130 rr(88)=zkre(93)*c(39)*c(3)
2131 rr(89)=zkre(94)*c(39)*c(29)
2132 rr(90)=zkre(95)*c(39)*c(30)
2133 rr(91)=zkre(96)*c(39)*c(45)
2134 rr(92)=zkre(97)*c(39)*c(14)
2135 rr(93)=zkre(98)*c(39)*c(8)
2136 rr(94)=zkre(99)*c(39)*c(12)
2137 rr(95)=zkre(100)*c(39)*c(19)
2138 rr(96)=zkre(101)*c(39)*c(27)
2139 rr(97)=zkre(102)*c(39)*c(18)
2140 rr(98)=zkre(103)*c(23)*c(2)/c(46)
2141 rr(99)=zkre(104)*c(2)*c(25)*c(46)
2142 rr(100)=(zkre(105)*c(46)+zkre(106))*c(2)*c(24)
2143 rr(101)=zkre(107)*c(29)*c(3)+zkre(118)*c(3)*c(30)
2144 rr(102)=zkre(108)*c(39)*c(35)
2145 rr(103)=zkre(109)*c(2)*c(31)
2146 rr(104)=zkre(110)*con(1)*c(21)
2147
2148 if (c(46) .ge. 1.0d-3) then
2149 rr(105)=zkre(111)*con(3)*con(1)*c(46)**0.5d0
2150 arytm=1.0d0
2151 else
2152 rr(105)=zkre(112)*c(8)*c(2)*c(46)
2153 arytm=0.0d0
2154 end if
2155
2156 rr(106)=zkre(113)*c(16)*c(2)+zkre(119)*c(16)*c(3)
2157 rr(107)=zkre(114)*c(42)*c(45)
2158 rr(108)=zkre(115)*c(42)*c(28)
2159 rr(109)=zkre(116)*c(2)*c(36)*chlorine+ &
2160 zkre(116)*c(3)*c(36)*chlorine
2161
2162 return
2163 end subroutine react
2164
2165
2166
2167
2168 !************************************************************************
2169 ! this program contains the routine mass which calculates the mass
2170 ! fluxes for the mass balances.
2171 ! inputs : wl, radius, temp, gcon(21), con(28), akeq(17),akhen(21)
2172 ! outputs : fgl(21),flg(21)
2173 !************************************************************************
2174
2175 subroutine mass(wl,radius,temp,gcon,con,c,akeq,akhen,fgl,flg, &
2176 gfgl,gflg)
2177
2178 ! arguments
2179 double precision wl, radius, temp
2180 double precision gcon(22), con(28), c(46), akeq(17), akhen(21)
2181 double precision fgl(21), flg(21), gfgl(21), gflg(21)
2182
2183 ! local variables
2184 integer i
2185 double precision acc, airl, dg, rideal
2186 ! ekhen(i) is the effective henry's law constant
2187 ! dimension ekhen(21)
2188 double precision ekhen(21)
2189 double precision kn,n,ikn,kmt
2190
2191
2192 ekhen(1)=akhen(1)*(1.d0+akeq(1)/c(46)+akeq(1)*akeq(2)/c(46)**2)
2193 ekhen(2)=1.0d30
2194 ekhen(3)=akhen(3)*(1.d0+akeq(7)/c(46))
2195 ekhen(4)=akhen(4)*(1.d0+akeq(6)/c(46))
2196 ekhen(5)=akhen(5)*(1.d0+akeq(8)/c(46)+akeq(8)*akeq(9)/c(46)**2)
2197 ekhen(6)=akhen(6)*(1.d0+akeq(5)/c(46))
2198 ekhen(7)=akhen(7)*((1.d0+akeq(12))/akeq(12))
2199 ekhen(8)=akhen(8)*(1.d0+akeq(13)/c(46))
2200 ekhen(9)=akhen(9)
2201 ekhen(10)=akhen(10)
2202 ekhen(11)=akhen(11)
2203 ekhen(12)=akhen(12)
2204 ekhen(13)=akhen(13)
2205 ekhen(14)=akhen(14)
2206 ekhen(15)=akhen(15)*(1.d0+akeq(14)/c(46)+(akeq(14)*c(37))/ &
2207 (akeq(16)*c(46)))
2208 ekhen(16)=akhen(16)
2209 ekhen(17)=akhen(17)*(1.d0+akeq(15)/c(46))
2210 ekhen(18)=akhen(18)
2211 ekhen(19)=akhen(19)*(1.d0+akeq(10)/c(45))
2212 ekhen(20)=akhen(20)
2213 ekhen(21)=akhen(21)
2214
2215 ! airl is the mean free path of air. later we have to substitute
2216 ! the numerical value given here by a function of temperature
2217 ! airl=65x10-9 m
2218 airl=65.d-9
2219 ! kn is the knudsen number
2220 kn=airl/radius
2221 ikn=1.d0/kn
2222 ! acc is the accomodation coefficient assumed the same here for
2223 ! all the species
2224 acc=0.1d0
2225 ! n is the coefficient entering the flux expression
2226 n=1.d0/(1.d0+((1.33d0+0.71d0*ikn)/(1.d0+ikn)+4.d0*(1.d0-acc) &
2227 /(3.d0*acc))*kn)
2228 ! dg is the gas phase diffusivity assumed here the same for all
2229 ! the gases.we shall probably have to change it later.
2230 ! dg=1.x10-5 m**2/sec
2231 dg=1.0d-5
2232 ! rideal is the gas constant in [atm/K/(mol/liter)]
2233 rideal=0.082058d0
2234 kmt=(3.0d0*n*dg)/(radius*radius)
2235
2236 do 1500 i=1,21
2237 fgl(i)=kmt*gcon(i)
2238 flg(i)=(kmt*con(i))/(ekhen(i)*rideal*temp)
2239 gfgl(i)=fgl(i)*wl
2240 gflg(i)=flg(i)*wl
2241 1500 continue
2242
2243
2244
2245 return
2246 end subroutine mass
2247
2248
2249
2250
2251 !***********************************************************************
2252 ! this routine calculates the differentials dc(i)/dt for the 28
2253 ! main species in the aqueous phase and the 21 species in the gas phase.
2254 ! units for all rates are (mol/lt.s)
2255 ! note that there are no reaction terms for the gas phase.
2256 ! revised 23 nov 1988
2257 !***********************************************************************
2258 !
2259 subroutine differ(rp,rl,fgl,flg,gfgl,gflg,dp,dl)
2260
2261 ! arguments
2262 double precision rp(28),rl(28),fgl(21),flg(21),gfgl(21),gflg(21)
2263 double precision dp(49),dl(49)
2264
2265 ! local variables
2266 integer i
2267
2268
2269 do 1510 i=1,21
2270 dp(i)=rp(i)+fgl(i)
2271 dl(i)=rl(i)+flg(i)
2272 1510 continue
2273
2274 do 1520 i=22,28
2275 dp(i)=rp(i)
2276 dl(i)=rl(i)
2277 1520 continue
2278
2279 do 1530 i=29,49
2280 dp(i)=gflg(i-28)
2281 dl(i)=gfgl(i-28)
2282 1530 continue
2283
2284
2285 return
2286 end subroutine differ
2287
2288
2289
2290
2291 ! ***********************************************************************
2292 ! the routine addit sums up the rates of
2293 ! the 120 reactions to give the rates for the 28 main species.
2294 ! input : the 120 reaction rates from the sub react
2295 ! output : the 28 formation and destruction rates
2296 ! revised 23 nov 1988
2297 ! ***********************************************************************
2298
2299 subroutine addit(rr,arytm,rp,rl)
2300
2301 ! arguments
2302 double precision arytm
2303 double precision rr(120),rp(28),rl(28)
2304
2305
2306 !
2307 ! ** s(iv) **
2308 !
2309 rp(1)=rr(107)
2310 rl(1)=+rr(72)+rr(73)+rr(74)+rr(98)+rr(101)+rr(105)*arytm &
2311 +rr(76)+rr(77)+rr(82)+rr(87)+rr(99)+rr(100)+2.0d0*rr(103) &
2312 +rr(104)+2.0d0*rr(105)*(1.0d0-arytm)+rr(106)+rr(109) &
2313 +rr(75)+rr(88)
2314 !
2315 ! ** s(vi) **
2316 !
2317 rp(2)= rr(72)+rr(73)+rr(74)+rr(98)+rr(101)+rr(105)*arytm &
2318 +rr(85) &
2319 +2.d0*rr(82)+rr(84)+rr(86)+rr(87)+rr(88)+rr(89)+rr(90) &
2320 +rr(91)+rr(92)+rr(93)+rr(94)+rr(95)+rr(96)+rr(97)+rr(99) &
2321 +rr(100)+rr(102)+rr(103)+rr(104)
2322 rl(2)=0.0d0
2323 !
2324 ! ** n(iii) **
2325 !
2326 rp(3)=2.0d0*rr(31)+rr(32)+rr(33)+2.0d0*rr(104)
2327 rl(3)=rr(35)+rr(37)+rr(39)+rr(105)*arytm+rr(36)+rr(38)+rr(40) &
2328 +rr(41)+rr(42)+rr(43)+rr(85)+rr(93)+rr(105)*(1.0d0-arytm)
2329 !
2330 ! ** n(v) **
2331 !
2332 rp(4)=rr(32)+rr(34)+rr(39)+rr(40)+rr(43)+rr(46)+rr(47)+rr(48)+ &
2333 rr(49)+rr(54)+rr(59)+rr(62)+rr(71)+rr(85)+rr(103)
2334 rl(4)=rr(44)
2335 !
2336 ! ** co2 **
2337 !
2338 rp(5)=rr(52)+rr(54)+rr(55)+rr(56)+rr(57)+rr(58)+rr(59)+ &
2339 rr(60)+rr(61)+rr(79)+rr(80)+rr(95)+rr(97)+rr(19)+rr(20)+rr(60)+ &
2340 rr(68)+rr(41)
2341 rl(5)=rr(17)+rr(18)+rr(94)
2342 !
2343 ! ** h2o2 **
2344 !
2345 rp(6)=rr(2)+rr(6)+rr(7)+rr(8)+rr(18)
2346 rl(6)=rr(1)+rr(5)+rr(9)+rr(10)+rr(16)+rr(20)+rr(29)+rr(39)+ &
2347 rr(48)+rr(53)+rr(73)+rr(92)+rr(15)+rr(28)
2348 !
2349 ! ** hcho **
2350 !
2351 rp(7)= rr(65)+rr(67)+rr(68)+rr(69)+rr(70)+rr(71)+rr(102)+rr(107) &
2352 +rr(108)
2353 rl(7)= rr(106)+rr(50)+rr(51)
2354 !
2355 ! ** hcooh **
2356 !
2357 rp(8)=rr(50)
2358 rl(8)=rr(52)+rr(53)+rr(54)+rr(55)+rr(56)+rr(79)+rr(97)+rr(57)+ &
2359 rr(58)+rr(59)+rr(60)+rr(61)+rr(80)+rr(95)
2360 !
2361 ! ** no **
2362 !
2363 rp(9)=rr(35)+rr(36)+rr(45)
2364 rl(9)=rr(31)+rr(33)
2365 !
2366 ! ** no2 **
2367 !
2368 rp(10)=rr(37)+rr(38)+rr(41)+rr(42)+rr(43)+rr(44)+rr(93)
2369 rl(10)=rr(31)+2.0d0*rr(32)+rr(34)+2.0d0*rr(104)
2370 !
2371 ! ** o3 **
2372 !
2373 rp(11)=0.0d0
2374 rl(11)=rr(2)+rr(11)+rr(12)+rr(13)+rr(14)+rr(15)+rr(16)+rr(40)+ &
2375 rr(51)+rr(55)+rr(58)+rr(72)
2376 !
2377 ! ** pan **
2378 !
2379 rp(12)=0.0d0
2380 rl(12)=rr(62)+rr(98)
2381 !
2382 ! ** ch3coooh **
2383 !
2384 rp(13)=0.0d0
2385 rl(13)=rr(100)
2386 !
2387 ! ** ch3ooh **
2388 !
2389 rp(14)=rr(63)+rr(64)
2390 rl(14)=rr(65)+rr(66)+rr(70)+rr(99)
2391 !
2392 ! ** hcl **
2393 !
2394 rp(15)=rr(22)+2.d0*rr(25)+2.d0*rr(26)+rr(27)+rr(29)+2.d0*rr(30) &
2395 +2.d0*rr(42)+2.d0*rr(56)+2.d0*rr(61)+ &
2396 2.d0*rr(69)+2.d0*rr(109)+2.d0*rr(28)
2397 rl(15)=rr(21)+rr(49)+rr(86)+rr(96)
2398 !
2399 ! ** oh **
2400 !
2401 rp(16)=2.0d0*rr(1)+rr(9)+rr(10)+rr(12)+ &
2402 rr(13)+rr(15)+rr(22)+rr(30)+rr(35)+rr(36)+rr(44)+rr(55)+rr(58)+ &
2403 rr(65)+rr(91)+rr(101)
2404 rl(16)=rr(3)+rr(4)+rr(5)+rr(11)+rr(17)+rr(21)+rr(33)+rr(34)+ &
2405 rr(37)+rr(38)+rr(50)+rr(52)+rr(57)+rr(66)+rr(67)+rr(75)+rr(76)+ &
2406 rr(83)+rr(108)
2407 !
2408 ! ** ho2 **
2409 !
2410 rp(17)=rr(5)+rr(11)+rr(20)+rr(28)+rr(29)+rr(48)+rr(50)+rr(52)+ &
2411 rr(54)+rr(55)+rr(56)+rr(57)+rr(59)+rr(60)+rr(61)+rr(65)+ &
2412 rr(67)+rr(68)+rr(69)+rr(71)+rr(79)+rr(92)+rr(95)+rr(97)+ &
2413 rr(102)+rr(14)+rr(14)+rr(15)+rr(58)+rr(80)
2414 rl(17)=rr(3)+2.0d0*rr(6)+rr(7)+rr(9)+rr(12)+rr(25)+rr(27)+rr(46) &
2415 +rr(63)+rr(89)+rr(101)+rr(4)+rr(7)+2.0d0*rr(8)+rr(10)+rr(13) &
2416 +rr(18)+rr(19)+rr(26)+rr(47)+rr(64)+rr(78)+rr(90)
2417 !
2418 ! ** no3 **
2419 !
2420 rp(18)=0.0d0
2421 rl(18)= rr(43)+rr(45)+rr(46)+rr(47)+rr(48)+rr(49)+rr(54)+ &
2422 rr(59)+rr(71)+rr(103)
2423 !
2424 ! ** nh3 **
2425 !
2426 rp(19)=0.0d0
2427 rl(19)=0.0d0
2428 !
2429 ! ** ch3o2 **
2430 !
2431 rp(20)=rr(66)
2432 rl(20)=rr(63)+rr(64)
2433 !
2434 ! ** ch3oh **
2435 !
2436 rp(21)=0.0d0
2437 rl(21)=rr(67)+rr(68)+rr(69)+rr(71)+rr(102)
2438 !
2439 ! ** cl2-, cl **
2440 !
2441 rp(22)=rr(49)+rr(96)+rr(23)
2442 rl(22)=rr(25)+rr(26)+rr(28)+rr(30)+rr(42)+rr(56)+rr(61)+ &
2443 rr(69)+rr(109)+rr(24)+rr(27)+rr(29)
2444 !
2445 ! ** cloh- **
2446 !
2447 rp(23)=rr(21)+rr(24)
2448 rl(23)=rr(22)+rr(23)
2449 !
2450 ! ** so4- **
2451 !
2452 rp(24)=2.d0*rr(81)+rr(103)
2453 rl(24)= rr(84)+rr(87)+rr(88)+rr(89)+rr(90)+rr(91)+ &
2454 rr(92)+rr(93)+rr(94)+rr(95)+rr(96)+rr(97)+rr(102)
2455 !
2456 ! ** so5- **
2457 !
2458 rp(25)=rr(75)+rr(76)+rr(83)+rr(84)+rr(87)+rr(88)+rr(108)+rr(109)
2459 rl(25)=rr(78)+rr(79)+rr(80)+2.0d0*rr(81)
2460 !
2461 ! ** hso5- **
2462 !
2463 rp(26)=rr(77)+rr(78)+rr(79)+rr(80)
2464 rl(26)=+rr(82)+rr(83)+rr(84)+rr(85)+rr(86)
2465 !
2466 ! ** hoch2so3- **
2467 !
2468 rp(27)=rr(106)
2469 rl(27)=rr(107)+rr(108)
2470 !
2471 ! ** co3- **
2472 !
2473 rp(28)=rr(17)+rr(18)+rr(94)
2474 rl(28)=rr(19)+rr(20)+rr(41)+rr(60)+rr(68)
2475
2476 return
2477 end subroutine addit
2478 !----------------------------------------------------------------------
2479
2480
2481
2482 end module module_cmu_bulkaqchem