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