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.
7 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************
10 module module_cmu_bulkaqchem
13 use module_peg_util, only: peg_error_fatal, peg_message
23 !**********************************************************************
25 ! revised aqueous-phase chemistry module for the 3-d acid deposition model
27 !**********************************************************************
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
34 ! (2) the water content is assumed to be constant (an input to the code)
36 ! ****** for bulk ******
38 !-----------------------------------------------------------------------
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.
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
49 ! this version is completely double precision, including arguments
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)
61 ! deleted the call to hybrd, so iradical_in>0 results in iradical=100
62 !-----------------------------------------------------------------------
67 !----------------------------------------------------------------------
68 subroutine aqoperator1( &
69 istat_fatal, istat_warn, &
71 gas, aerosol, gas_aqfrac, &
73 co2_mixrat_in, photo_in, iradical_in, idecomp_hmsa_hso5, iaq, &
74 yaq_beg, yaq_end, ph_cmuaq_cur )
76 use module_data_cmu_bulkaqchem, only: &
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, &
87 wh2o2, wh2so4, whcho, whcl, whcooh, whno3, &
93 ! tbeg_sec : start of integration (in sec)
94 ! tend_sec : end of integration time (in sec)
96 ! gas(ngas) : gas-phase concentration vector (in ppm)
97 ! aerosol(naers) : aqueous (aerosol) species concentration (in ug/m3)
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
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
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
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
121 ! gcon(ngas) : gas-phase concentrations (in ppm) (through rpar to rates)
124 !.....differential equations are solved for the following species
127 !-----------------------------------------------------------------------
128 ! 1. formaldehyde total 1. so2 1. s(iv)
129 ! 2. formic acid total 2. h2o2 2. h2o2(aq)
138 !.....y matrix structure everything is (ug/m3)
139 ! only 11 odes are solved for the whole distribution
141 ! yaq(1) = total formaldehyde
142 ! yaq(2) = total formic acid
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)
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
168 integer icount, iprint, negconc_count
171 ammonold, chlorold, crustal, &
172 dammon, dchlor, dhmsa, dhso5, dnit, dsulf, &
173 fammon, fh2o2, fh2so4, fhcho, fhcl, fhcooh, fhno3, fso2, &
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
181 duma, dumb, dumhion, vlwc
182 double precision gascon(ngas)
183 double precision yaq(meqn1max)
186 double precision rpar(8+ngas)
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))
205 ! calc temperature dependent parameters
207 pres = p ! pressure in atm
211 ! zero the yaq matrix
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
230 co2_mixrat = co2_mixrat_in
233 ! if (iradical .eq. 100) meqn1 = 20
234 if (iradical .gt. 0) meqn1 = 20
240 ! transfer of all gas-phase concentrations to gascon and then
241 ! through rpar to rates
248 ! unit conversion factors
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
257 ! *** beginning of aqueous-phase calculation
258 ! loading of the concentrations in the yaq vector
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)
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
272 ! aqueous-phase species
274 ! "aerosol" array holds the bulk aqueous chemical concentrations
275 ! (in ug/m3-of-air), so there is no "activation" or summing over sections
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
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
292 ! save the old aqueous-phase concentrations before the integration
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)
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
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)
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
323 !.....the yaq matrix has been initialized
327 ! write(27,*)' initial values of the yaqs'
328 ! write(27,*)lwc,rh,temp
330 ! write(27,*)i, yaq(i)
334 ! variables for integration
336 stime = tbeg_sec ! in seconds
337 stout = tend_sec ! in seconds
339 ! calculation of sodium (needed for the integration routine)
341 sodium = aerosol(nas)
343 ! calculate crustal species concentration inside droplets
344 ! (it is used in the aqueous-phase chemistry calculation to
345 ! estimate fe and mn concentrations
347 crustal = aerosol(nar)
351 ! save yaq to yaq_beg
376 rpar(8) = ph_cmuaq_cur
378 rpar(8+i) = gascon(i)
384 call aqintegr1( meqn1, yaq, stime, stout, rpar, ipar )
389 ph_cmuaq_cur = rpar(8)
392 ! calculate the differences
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
402 ! transfer new aqueous concentrations back to aerosol array
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
412 ! check if the algorithm resulted in negative concentrations
413 ! report the corrections at the warning file
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)
431 ! return gas yaq results to their original matrices
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")
438 gas(nghcho) = yaq(1)/fhcho ! total hcho (ppm)
439 gas(nghcooh) = yaq(2)/fhcooh ! total hcooh (ppm)
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)
446 gas(ngh2o2) = (yaq(4)+yaq(7))/fh2o2 ! h2o2(g) (ppm)
447 gas_aqfrac(ngh2o2) = yaq(7) / (yaq(4)+yaq(7))
449 gas(nga) = yaq(5)/fammon ! nh3(g) (ppm)
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)
456 vlwc = lwc*1.0e-6 ! (liter-h2o)/(liter-air)
458 duma = rideal*temp*vlwc*akhen(7)
459 gas_aqfrac(nghcho) = duma/(duma + 1.0)
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)
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=
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
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))
493 ! save yaq to yaq_end
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
508 if (ipar(3) .ne. 0) istat_fatal = istat_fatal - 10
509 if (ipar(7) .ne. 0) istat_fatal = istat_fatal - 200
512 if (ipar(5) .ne. 0) istat_warn = istat_warn + 10
513 if (negconc_count .ne. 0) istat_warn = istat_warn + 200
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
526 end subroutine aqoperator1
530 !----------------------------------------------------------------------
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
537 ! last modification: june 13, 1998 by s.p.
539 !----------------------------------------------------------------------
540 subroutine aqintegr1( neqa, y, stime, stout, rpar, ipar )
542 use module_data_cmu_bulkaqchem, only: &
543 iopt, itask, itol, liw1, lrw1, &
544 mdiag_svode, meqn1, meqn1max, mf, &
545 tola, tolr, worki, workr
547 use module_cmu_dvode_solver, only: dvode
550 integer neqa, ipar(*)
551 double precision y(meqn1max), stime, stout, rpar(*)
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)
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
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
571 ! ready for the call to svode
575 call dvode( fex1, neqa, y, stime, stout, itol, rtol, atol, itask, &
576 istate, iopt, rwork, lrw, iwork, liw, jex, mf, &
579 if (istate .ne. 2) then
580 ! write(*,*) '*** aqintegr1 -- svode failed'
581 ! write(*,*) ' stime, istate =', stime, istate
584 if (mdiag_svode .gt. 0) then
586 '*** module_cmuaq_bulk, aubr aqintegr1 -- ' // &
587 'svode failed, istate =', istate
589 ipar(3) = ipar(3) + 1
594 if (y(i) .lt. 0.0) y(i) = 1.e-20
598 end subroutine aqintegr1
602 !----------------------------------------------------------------------
603 subroutine fex1( neq, t, y, f, rpar, ipar )
605 use module_data_cmu_bulkaqchem, only: meqn1, meqn1max
609 double precision t, y(meqn1max), f(meqn1max), rpar(*)
612 double precision a(meqn1max),b(meqn1max)
615 call aqratesa( meqn1, t, y, a, b, f, rpar, ipar )
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 ***' )
632 !----------------------------------------------------------------------
634 ! last modification : feb. 15, 1995 by s. pandis
635 !*************************************************************************
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
643 subroutine aqratesa( neq, t, yaq, aqprod, aqdest, yaqprime, &
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, &
659 wh2o2, whcho, whcooh, wnh3, wso2, wmol, wvol, &
664 double precision t, yaq(meqn1max), yaqprime(meqn1max)
665 double precision aqprod(meqn1max), aqdest(meqn1max)
666 double precision rpar(*)
669 integer nfev, info, i, j, n
670 integer icount, iprint
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)
681 double precision, parameter :: one=1.
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)
695 double precision gascon(ngas)
699 ! unload ipar, rpar values
712 gascon(i) = rpar(8+i)
715 temp_cmuaq_cur = temp
716 pres_cmuaq_cur = pres
719 ! calculation of the temperature for this time (temp)
720 ! (assuming linear temperature change for each operator stepp)
725 call qsaturation (temp, qsat) ! qsat is in ug/m3
727 ! zero the rate of change vectors
735 ! set dummy ph to zero for printing only
741 tlwc = watcont*1.e6 ! in ug/m3
742 vlwc=tlwc*1.e-12 ! vol/vol
744 ! check for negative values
747 !sp if (yaq(i) .le. 0.) yaq(i)=1.e-20
750 ! reconstruct the matrices using the available yaq values
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
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)
783 ! calculation of gcon vector (gas-phase concentrations in mole/l)
786 ! gcon(i) = spres(i)*1.e-6/(0.08206*temp)
787 gcon(i) = spres(i)*1.e-6/(rideal*temp)
790 ! ** radius and lwc for the section
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
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.
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
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)
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
813 ! loading of the main aqueous concentrations (in m)
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)
846 ! set a minimum concentration (to avoid divisions by zero)
849 if (con(i) .lt. 1.0e-20) con(i)=1.0e-20
853 ! 27-oct-2005 rce - previously there was a bug here and the code
854 ! would continue on even if fullequil failed
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
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
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)
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 )
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
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
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
908 ! we calculate the ionic species concentrations
910 call values(hyd, con, cmet, akeq, c)
912 ! bypass the radical calculation by hybrd if necessary
914 if (iradical .eq. 0) go to 270
915 if (iabs(iradical-101) .le. 1) go to 280
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
923 ! set the radical concentrations to zero
926 ! if (info .eq. 2 .or. info .eq. 3 &
927 ! .or. info .eq. 4 .or. info .eq. 5 .or. iradical .eq. 0) then
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
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
952 !sp con(22)=pro/destr
957 call values(hyd, con, cmet, akeq,c)
959 ! final calculation of reaction rates
961 call react(c,cmet,con,akre,rr,arytm)
963 ! calculation of net production and consumption rates
965 call addit(rr, arytm, rp, rl)
967 ! calculation of mass transfer rates
969 call mass(wvol,rad,temp,gcon,con,c,akeq,akhen,fgl,flg,gfgl,gflg)
971 ! calculation of net rates of change
973 call differ(rp,rl,fgl,flg,gfgl,gflg,dp,dl)
975 ! calculation of right hand sides of the derivatives
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)
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)
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)
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)
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)
999 ! ** aqueous-phase species **
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)
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)
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)
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)
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)
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)
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)
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)
1033 if (iabs(iradical-101) .le. 1) then
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)
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)
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)
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)
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)
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)
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)
1064 ! calculation of appropriate destruction rate
1067 if (yaq(i) .le. 1.e-20) then
1071 aqdest(i) = aqdest(i)/yaq(i)
1074 ! change to avoid divisions by zero in integration
1077 if (aqdest(i) .le. 1.e-18) aqdest(i)=1.e-18
1081 ! calculation of characteristic times (used for debugging)
1084 !db do 110 i=1, meqn1
1086 !db if (aqdest(i) .le. 1.e-10) go to 110
1087 !db tchar=1./aqdest(i)
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)
1094 !db if (tchar .lt. tsm) then
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
1106 ! sulfrate=yaqprime(3)/gmol(1)+yaqprime(6)/wmol(1)+ &
1107 ! yaqprime(11)/wmol(2)
1108 ! rsrate=sulfrate/(abs(yaqprime(3))+abs(yaqprime(11)) + &
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, *)'************************'
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
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))
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)*******'
1150 !sp write(3,*)i,rr(i)*1.e6*3600.
1154 !120 format(1x,2(1x,f8.4))
1155 120 format( 'aqratesa - tmin, yaq(11)=so4', 2(1x,f8.4) )
1158 ! load ipar, rpar values
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)
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)
1173 end subroutine aqratesa
1178 !************************************************************************
1179 ! this routine calculates the the steady state species concentrations
1180 !************************************************************************
1181 subroutine steady(radius,temp,c,con,gcon,akeq,akhen,akre)
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
1190 ! x(8) the values of the steady state species concentrations
1194 double precision radius, temp
1195 double precision c(46),gcon(22),akeq(17),akhen(21),akre(120)
1196 double precision con(28)
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)
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
1210 ! kn is the knudsen number
1213 ! acc is the accomodation coefficient assumed the same here for
1216 ! n is the coefficient entering the flux expression
1217 n=1.0/(1.+((1.33+0.71*ikn)/(1.+ikn)+4.*(1.-acc) &
1219 ! dg is the gas phase diffusivity assumed here the same for all
1220 ! the gases. dg=1.x10-5 m**2/sec
1222 ! rideal is the gas constant in [atm/K/(mol/liter)]
1224 kmt=(3.0*n*dg)/(radius*radius)
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))
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))
1250 a1=c(46)/(akeq(15)+c(46))
1251 a2=akeq(15)/(akeq(15)+c(46))
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))
1277 ho2=(x(2)*c(46))/(akeq(15)+c(46))
1278 o2=(x(2)*akeq(15))/(akeq(15)+c(46))
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))
1288 x(1)=(2.*akre(1)*c(14)+akre(15)*c(15)*c(22)+akre(30)*c(45)* &
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)
1305 x(4)=(akre(21)*c(27)*x(1)+akre(24)*c(37))/(akre(22)+ &
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)+ &
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))
1332 end subroutine steady
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 ! ***********************************************************************
1348 subroutine fullequil(acon,aspres,acmet,aakeq,aakhen,awv, &
1351 use module_data_cmu_bulkaqchem, only: mdiag_fullequil
1355 double precision acon(28), aspres(21), acmet(4), aakeq(17), aakhen(21)
1356 double precision awv, atemp, axsol
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)
1365 ! change variables to double precision to avoid low ph errors
1391 ! we find the initial interval [aa,bb] for the bisection method
1392 ! new version (31/10/87)
1394 call electro(x,fa,con,spres,cmet,akeq,akhen,wv,temp)
1398 do 1035 i = -(14+idum_01), (1+idum_01)
1400 call electro(x,f,con,spres,cmet,akeq,akhen,wv,temp)
1401 if (f*fa .ge. 0.0d0) then
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.
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
1420 '*** module_cmuaq_bulk - mistake in fullequil with ipass_01 = 1'
1424 else if (ipass_01 .eq. 2) then
1426 '*** module_cmuaq_bulk - mistake in fullequil with ipass_01 = 2'
1432 ! otherwise, report the error and exit
1433 ipar(7) = ipar(7) + 1
1434 if (mdiag_fullequil .gt. 0) then
1436 '*** module_cmuaq_bulk - mistake in fullequil - con, cmet ='
1437 write(6,*) con, cmet
1443 ! bisection method for the solution of the equation
1444 ! rtol : relative tolerance
1446 ! 02-nov-2005 rce - smaller rtol for h+ makes dvode run faster
1449 1050 error= dabs(bb-aa)/aa
1451 if (error .le. rtol) then
1453 axsol=xsol ! single precision
1457 call electro(xm,fm,con,spres,cmet,akeq,akhen,wv,temp)
1458 if (fa*fm .gt. 0.0d0) then
1466 end subroutine fullequil
1471 ! ***********************************************************************
1472 ! routine that gives values of the electroneutrality equation
1473 ! called by fullequil
1474 ! ***********************************************************************
1476 subroutine electro(x,f,con,spres,cmet,zkeq,zkhen,wv,temp)
1478 use module_data_cmu_bulkaqchem, only: &
1479 mprescribe_ph, rideal, xprescribe_ph
1483 ! original subr arguments were akeq & akhen
1484 ! renamed them to zkeq & zkhen
1485 ! to avoid conflict with module_data_cmu_bulkaqchem
1487 double precision x, f, wv, temp
1488 double precision con(28),spres(21),cmet(4),zkeq(17),zkhen(21)
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)
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--
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
1505 cc(10)=(zkeq(6)*con(4))/(x+zkeq(6)) ! no3-
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
1511 cc(15)=(zkeq(5)*con(6))/(x+zkeq(5)) ! ho2-
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
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-
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-
1537 cc(33)=(zkeq(10)*x*con(19))/(zkeq(11)+zkeq(10)*x) ! nh4+
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)
1548 if (mprescribe_ph .gt. 0) then
1549 f = 10.0**(-xprescribe_ph) - x
1553 end subroutine electro
1557 !----------------------------------------------------------------------
1559 ! routines used by the aqeous-phase module
1561 ! 1. dropinit : initialization
1565 use module_data_cmu_bulkaqchem, only: &
1567 wmol, wh2o2, wh2so4, whcho, whcl, whcooh, whno3, wnh3, wso2
1574 ! loading of molecular weights
1593 wmol(7)= 48.0e0 ! was previously 60.0e0
1627 ! 09-nov-2005 rce - set gmol(6) == wh2o2 to conserve h2o2
1647 end subroutine dropinit
1651 !----------------------------------------------------------------------
1652 subroutine qsaturation(temp,qsat)
1654 ! this routine calculates the saturation mass concentration (in ug/m3)
1655 ! over liquid water for a temperature temp (k)
1659 double precision temp,qsat
1662 double precision psat, t ! these should be double precision ?
1663 double precision rideal,a0,a1,a2,a3,a4,a5,a6,esat,csat
1665 t = temp-273.15 ! in c
1666 rideal = 0.082058d0 ! gas constant in [atm/K/(mol/liter)]
1673 a6 = 6.136820929d-11
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
1681 end subroutine qsaturation
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 !************************************************************************
1695 subroutine constants(temp)
1697 use module_data_cmu_bulkaqchem, only: akeq, akhen, akre, &
1698 maqurxn_all, maqurxn_sulf1, mopt_eqrt_cons
1701 double precision temp
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)
1709 data dheq/1960.,1500.,0.,2720.,-3730.,8700.,-1260., &
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, &
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, &
1720 data bkhen/1.23e0,0.0e0,49.e0,2.1e5,3.4e-2,7.45e4,6.3e3, &
1722 0.01e0,1.13e-2,2.9e0,473.e0,227.e0,727.e0,25.e0,2000.e0,2.1e5, &
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/
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/
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
1787 akeq(i)=bkeq(i)*exp(dheq(i)*(1.0/temp-1.0/298.0))
1790 akhen(i)=bkhen(i)*exp(dhhen(i)*(1.0/temp-1.0/298.0))
1793 akre(i)=bkre(i)*exp(dhre(i)*(1.0/temp-1.0/298.0))
1796 ! turn reactions on/off selectively
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
1803 if (iusei .le. 0) akre(i) = 0.0
1804 if (iusei .le. 0) bkre(i) = 0.0
1808 end subroutine constants
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 !*************************************************************************
1820 subroutine values(x,con,cmet,akeq,cc)
1824 double precision con(28),cmet(4),akeq(17),cc(46)
1827 double precision bparam,cparam,diak,cl,hcl
1829 ! species in the aqueous phase mechanism
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(-)
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(+)
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(-)
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))
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))
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))
1887 cc(7)=(x*con(3))/(x+akeq(7))
1888 cc(8)=(akeq(7)*con(3))/(x+akeq(7))
1890 cc(9)=(x*con(4))/(x+akeq(6))
1891 cc(10)=(akeq(6)*con(4))/(x+akeq(6))
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))
1897 cc(14)=(x*con(6))/(x+akeq(5))
1898 cc(15)=(akeq(5)*con(6))/(x+akeq(5))
1900 cc(16)=con(7)/(1.0d0+akeq(12))
1901 cc(17)=(akeq(12)*con(7))/(1.0d0+akeq(12))
1903 cc(18)=(x*con(8))/(x+akeq(13))
1904 cc(19)=(akeq(13)*con(8))/(x+akeq(13))
1919 cc(27)=(akeq(14)*hcl)/x
1923 cc(29)=(x*con(17))/(x+akeq(15))
1924 cc(30)=(akeq(15)*con(17))/(x+akeq(15))
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)
1935 cc(36)=(akeq(14)*cl*hcl)/(akeq(16)*x)
1942 cc(42)=(x*con(27))/(x+akeq(17))
1943 cc(43)=(akeq(17)*con(27))/(x+akeq(17))
1949 end subroutine values
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 !************************************************************************
1962 subroutine react(c,cmet,con,zkre,rr,arytm)
1964 use module_data_cmu_bulkaqchem, only: chlorine, kiron, photo, &
1969 ! original argument was akre
1970 ! renamed it to zkre
1971 ! to avoid conflict with module_data_cmu_bulkaqchem
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
1978 double precision ph, r1, r2, r3, r4, r5, sn
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)
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
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)
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))
2064 ! rate expressions for the metal catalysed oxidation of s(iv)
2066 ! ** phenomenological expression by martin et al. (1991) **
2068 if (kiron .eq. 1) then
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)
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)))
2088 if (cmet(1)*cmet(2) .lt. 1.0d-15) then
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))
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)
2106 rr(74)=zkre(78)*cmet(2)*c(2)
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
2115 ! ** end of martin's expression **
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)
2148 if (c(46) .ge. 1.0d-3) then
2149 rr(105)=zkre(111)*con(3)*con(1)*c(46)**0.5d0
2152 rr(105)=zkre(112)*c(8)*c(2)*c(46)
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
2163 end subroutine react
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 !************************************************************************
2175 subroutine mass(wl,radius,temp,gcon,con,c,akeq,akhen,fgl,flg, &
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)
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
2192 ekhen(1)=akhen(1)*(1.d0+akeq(1)/c(46)+akeq(1)*akeq(2)/c(46)**2)
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))
2206 ekhen(15)=akhen(15)*(1.d0+akeq(14)/c(46)+(akeq(14)*c(37))/ &
2209 ekhen(17)=akhen(17)*(1.d0+akeq(15)/c(46))
2211 ekhen(19)=akhen(19)*(1.d0+akeq(10)/c(45))
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
2219 ! kn is the knudsen number
2222 ! acc is the accomodation coefficient assumed the same here for
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) &
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
2232 ! rideal is the gas constant in [atm/K/(mol/liter)]
2234 kmt=(3.0d0*n*dg)/(radius*radius)
2238 flg(i)=(kmt*con(i))/(ekhen(i)*rideal*temp)
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 !***********************************************************************
2259 subroutine differ(rp,rl,fgl,flg,gfgl,gflg,dp,dl)
2262 double precision rp(28),rl(28),fgl(21),flg(21),gfgl(21),gflg(21)
2263 double precision dp(49),dl(49)
2286 end subroutine differ
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 ! ***********************************************************************
2299 subroutine addit(rr,arytm,rp,rl)
2302 double precision arytm
2303 double precision rr(120),rp(28),rl(28)
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) &
2317 rp(2)= rr(72)+rr(73)+rr(74)+rr(98)+rr(101)+rr(105)*arytm &
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)
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)
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)
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)+ &
2341 rl(5)=rr(17)+rr(18)+rr(94)
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)
2351 rp(7)= rr(65)+rr(67)+rr(68)+rr(69)+rr(70)+rr(71)+rr(102)+rr(107) &
2353 rl(7)= rr(106)+rr(50)+rr(51)
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)
2363 rp(9)=rr(35)+rr(36)+rr(45)
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)
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)
2380 rl(12)=rr(62)+rr(98)
2389 rp(14)=rr(63)+rr(64)
2390 rl(14)=rr(65)+rr(66)+rr(70)+rr(99)
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)
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)+ &
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)
2421 rl(18)= rr(43)+rr(45)+rr(46)+rr(47)+rr(48)+rr(49)+rr(54)+ &
2422 rr(59)+rr(71)+rr(103)
2432 rl(20)=rr(63)+rr(64)
2437 rl(21)=rr(67)+rr(68)+rr(69)+rr(71)+rr(102)
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)
2447 rp(23)=rr(21)+rr(24)
2448 rl(23)=rr(22)+rr(23)
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)
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)
2463 rp(26)=rr(77)+rr(78)+rr(79)+rr(80)
2464 rl(26)=+rr(82)+rr(83)+rr(84)+rr(85)+rr(86)
2469 rl(27)=rr(107)+rr(108)
2473 rp(28)=rr(17)+rr(18)+rr(94)
2474 rl(28)=rr(19)+rr(20)+rr(41)+rr(60)+rr(68)
2477 end subroutine addit
2478 !----------------------------------------------------------------------
2482 end module module_cmu_bulkaqchem