Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : !! Calculates particle production rates due to nucleation <rhompe>:
6 : !! binary homogeneous nucleation of sulfuric acid and water only
7 : !! Numerical method follows one of the following:
8 : !! 1. Zhao & Turco, JAS, V.26, No.5, 1995.
9 : !! 2. Vehkamaki, H., M. Kulmala, I. Napari, K.E.J. Lehtinen,
10 : !! C. Timmreck, M. Noppel and A. Laaksonen, 2002,
11 : !! An improved parameterization for sulfuric acid-water nucleation
12 : !! rates for tropospheric and stratospheric conditions,
13 : !! J. Geophys. Res., 107, 4622, doi:10.1029/2002jd002184
14 : !!
15 : !!
16 : !! @author Mike Mills, Chuck Bardeen
17 : !! @version May-2022
18 1418703437 : subroutine sulfnucrate(carma, cstate, iz, igroup, h2o, h2so4, beta1, beta2, radius_cluster, nucbin, nucrate, rc)
19 : use carma_precision_mod
20 : use carma_enums_mod
21 : use carma_constants_mod
22 : use carma_types_mod
23 : use carmastate_mod
24 : use carma_mod
25 : use sulfate_utils
26 :
27 : implicit none
28 :
29 : type(carma_type), intent(in) :: carma !! the carma object
30 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
31 : integer, intent(in) :: iz !! level index
32 : integer, intent(in) :: igroup !! group index
33 : real(kind=f), intent(out) :: h2o !! H2O concentrations in molec/cm3
34 : real(kind=f), intent(out) :: h2so4 !! H2SO4 concentrations in molec/cm3
35 : real(kind=f), intent(out) :: beta1
36 : real(kind=f), intent(out) :: beta2
37 : real(kind=f), intent(out) :: radius_cluster !! critical radius (cm)
38 : integer, intent(out) :: nucbin !! bin in which nucleation occurs
39 : real(kind=f), intent(out) :: nucrate !! nucleation rate #/x/y/z/s
40 : integer, intent(inout) :: rc !! return code, negative indicates failure
41 :
42 : ! Local declarations
43 : integer :: i, ibin, ie
44 : real(kind=f) :: nucrate_cgs ! binary nucleation rate, j (# cm-3 s-1)
45 : real(kind=f) :: cnum_h2so4 ! number of h2so4 molecules in the critical nucleus
46 : real(kind=f) :: cnum_tot ! total number of molecules in the critical nucleus
47 : real(kind=f) :: rb ! [erg/mol]
48 : real(kind=f) :: h2so4_cgs ! H2SO4 densities in g/cm3
49 : real(kind=f) :: h2o_cgs ! H2O densities in g/cm3
50 : real(kind=f) :: rh ! relative humidity (0-1)
51 : real(kind=f) :: mass_cluster_dry ! dry mass of the cluster ()
52 : real(kind=f) :: h2so4_bb ! bounded value of H2SO4 concentration in molec/cm3
53 : real(kind=f) :: temp_bb ! bounded value of temperature in Kelvins
54 : real(kind=f) :: rh_bb ! bounded value of relative humidity
55 : real(kind=f) :: ftry
56 :
57 : 5 format(/,'microfast::WARNING - nucleation rate exceeds 5.e1: ie=', i2,', iz=',i4,',lat=', &
58 : f7.2,',lon=',f7.2, ', rhompe=', e10.3)
59 :
60 : ! default values for outputs
61 1418703437 : nucbin = 1
62 1418703437 : nucrate = 0.0_f
63 1418703437 : nucrate_cgs = 0.0_f
64 1418703437 : radius_cluster = 0.0_f
65 1418703437 : mass_cluster_dry = 0.0_f
66 1418703437 : ftry = NOTSET
67 :
68 : !--------------------------------------------------------------
69 :
70 : ! beta1 and beta2 are calculated and used in sulfnucrate, and output for use in sulfhetnucrate
71 : ! kT/(2*Pi*M) = [erg/mol/K]*[K]/[g/mol] = [erg/g] = [cm2/s2]
72 : ! RB[erg/mol] = RGAS[erg/mol/K] * T[K] / (2Pi)
73 1418703437 : rb = RGAS * t(iz) / 2._f / PI
74 :
75 : ! Beta[cm/s] = sqrt(RB[erg/mol] / WTMOL[g/mol])
76 1418703437 : beta1 = sqrt(rb / gwtmol(igash2so4)) ! H2SO4
77 1418703437 : beta2 = sqrt(rb / gwtmol(igash2o)) ! H2O
78 :
79 : !--------------------------------------------------------------
80 :
81 : ! Compute H2SO4 densities in g/cm3
82 1418703437 : h2so4_cgs = gc(iz, igash2so4) / zmet(iz)
83 :
84 : ! Compute H2O densities in g/cm3
85 1418703437 : h2o_cgs = gc(iz, igash2o) / zmet(iz)
86 :
87 : ! Compute H2SO4 concentrations in molec/cm3
88 1418703437 : h2so4 = h2so4_cgs * AVG / gwtmol(igash2so4)
89 :
90 : ! Compute H2O concentrations in molec/cm3
91 1418703437 : h2o = h2o_cgs * AVG / gwtmol(igash2o)
92 :
93 : ! Compute relative humidity of water wrt liquid water
94 1418703437 : rh = (supsatl(iz, igash2o) + 1._f) !* 100._f
95 :
96 : ! Select nucleation method
97 1418703437 : select case (sulfnuclmethod)
98 : case('ZhaoTurco')
99 0 : call binary_nuc_zhao1995( carma, cstate, t(iz), wtpct(iz), rh, h2so4, h2so4_cgs, h2o, h2o_cgs, beta1, &
100 1418703437 : nucrate_cgs, mass_cluster_dry, radius_cluster, ftry, rc )
101 : case('Vehkamaki')
102 0 : if (h2so4 >= 1.0e4_f) then
103 0 : temp_bb = max( 230.0_f, min( 305.0_f, t(iz) ) )
104 0 : rh_bb = max( 1.0e-4_f, min( 1.0_f, rh ) )
105 0 : h2so4_bb = max( 1.0e4_f, min( 1.0e11_f, h2so4 ) )
106 :
107 : call binary_nuc_vehk2002( carma, temp_bb, rh_bb, h2so4_bb, nucrate_cgs, &
108 0 : mass_cluster_dry, radius_cluster )
109 : end if
110 : case default
111 0 : write(LUNOPRT,*)'sulfnucrate: '//trim(sulfnuclmethod)//' nucleation method no recognized'
112 0 : rc = RC_ERROR
113 2837406874 : return
114 : end select
115 :
116 : ! Calc bin # of crit nucleus
117 1418703437 : if (mass_cluster_dry.lt.rmassup(1,igroup)) then
118 565299250 : nucbin = 1
119 : else
120 853404187 : nucbin = 2 + int(log(mass_cluster_dry / rmassup(1,igroup)) / log(rmrat(igroup)))
121 : endif
122 :
123 : ! If none of the bins are large enough for the critical radius, then
124 : ! no nucleation will occur.
125 1418703437 : if (nucbin <= NBIN) then
126 : ! Scale to #z/s
127 1418462923 : nucrate = nucrate_cgs * zmet(iz)
128 : endif
129 :
130 : return
131 1418703437 : end subroutine sulfnucrate
132 :
133 :
134 : !----------------------------------------------------------------------
135 : !-----------------------------------------------------------------------
136 :
137 0 : subroutine binary_nuc_vehk2002( carma, temp, rh, h2so4, &
138 : nucrate_cgs, mass_cluster_dry, radius_cluster )
139 : !
140 : ! calculates binary nucleation rate and critical cluster size
141 : ! using the parameterization in
142 : ! Vehkamäki, H., M. Kulmala, I. Napari, K.E.J. Lehtinen,
143 : ! C. Timmreck, M. Noppel and A. Laaksonen, 2002,
144 : ! An improved parameterization for sulfuric acid-water nucleation
145 : ! rates for tropospheric and stratospheric conditions,
146 : ! J. Geophys. Res., 107, 4622, doi:10.1029/2002jd002184
147 : !
148 : use carma_precision_mod
149 : use carma_enums_mod
150 : use carma_constants_mod
151 : use carma_types_mod
152 : use carmastate_mod
153 : use carma_mod
154 : use sulfate_utils
155 :
156 : implicit none
157 :
158 : ! subr arguments (in)
159 : type(carma_type), intent(in) :: carma !! the carma object
160 : real(kind=f), intent(in) :: temp ! temperature (k)
161 : real(kind=f), intent(in) :: rh ! relative humidity (0-1)
162 : real(kind=f), intent(in) :: h2so4 ! concentration of h2so4 (molecules cm-3)
163 :
164 : ! subr arguments (out)
165 : real(kind=f), intent(out) :: nucrate_cgs ! binary nucleation rate, j (# cm-3 s-1)
166 : real(kind=f), intent(out) :: mass_cluster_dry ! the mass of cluster (g)
167 : real(kind=f), intent(out) :: radius_cluster ! the radius of cluster (cm)
168 :
169 : ! local variables
170 : real(kind=f) :: crit_x
171 : real(kind=f) :: acoe, bcoe, ccoe, dcoe, ecoe, fcoe, gcoe, hcoe, icoe, jcoe
172 : real(kind=f) :: tmpa, tmpb
173 : real(kind=f) :: cnum_h2so4 ! number of h2so4 molecules in the critical nucleus
174 : real(kind=f) :: cnum_tot ! total number of molecules in the critical nucleus
175 :
176 : ! executable
177 :
178 : ! calc sulfuric acid mole fraction in critical cluster
179 : crit_x = 0.740997_f - 0.00266379_f * temp &
180 : - 0.00349998_f * log (h2so4) &
181 : + 0.0000504022_f * temp * log (h2so4) &
182 : + 0.00201048_f * log (rh) &
183 : - 0.000183289_f * temp * log (rh) &
184 : + 0.00157407_f * (log (rh)) ** 2.0_f &
185 : - 0.0000179059_f * temp * (log (rh)) ** 2.0_f &
186 : + 0.000184403_f * (log (rh)) ** 3.0_f &
187 0 : - 1.50345e-6_f * temp * (log (rh)) ** 3.0_f
188 :
189 : ! calc nucleation rate
190 : acoe = 0.14309_f+2.21956_f*temp &
191 : - 0.0273911_f * temp**2.0_f &
192 0 : + 0.0000722811_f * temp**3.0_f + 5.91822_f/crit_x
193 :
194 : bcoe = 0.117489_f + 0.462532_f *temp &
195 : - 0.0118059_f * temp**2.0_f &
196 0 : + 0.0000404196_f * temp**3.0_f + 15.7963_f/crit_x
197 :
198 : ccoe = -0.215554_f-0.0810269_f * temp &
199 : + 0.00143581_f * temp**2.0_f &
200 : - 4.7758e-6_f * temp**3.0_f &
201 0 : - 2.91297_f/crit_x
202 :
203 : dcoe = -3.58856_f+0.049508_f * temp &
204 : - 0.00021382_f * temp**2.0_f &
205 : + 3.10801e-7_f * temp**3.0_f &
206 0 : - 0.0293333_f/crit_x
207 :
208 : ecoe = 1.14598_f - 0.600796_f * temp &
209 : + 0.00864245_f * temp**2.0_f &
210 : - 0.0000228947_f * temp**3.0_f &
211 0 : - 8.44985_f/crit_x
212 :
213 : fcoe = 2.15855_f + 0.0808121_f * temp &
214 : -0.000407382_f * temp**2.0_f &
215 : -4.01957e-7_f * temp**3.0_f &
216 0 : + 0.721326_f/crit_x
217 :
218 : gcoe = 1.6241_f - 0.0160106_f * temp &
219 : + 0.0000377124_f * temp**2.0_f &
220 : + 3.21794e-8_f * temp**3.0_f &
221 0 : - 0.0113255_f/crit_x
222 :
223 : hcoe = 9.71682_f - 0.115048_f * temp &
224 : + 0.000157098_f * temp**2.0_f &
225 : + 4.00914e-7_f * temp**3.0_f &
226 0 : + 0.71186_f/crit_x
227 :
228 : icoe = -1.05611_f + 0.00903378_f * temp &
229 : - 0.0000198417_f * temp**2.0_f &
230 : + 2.46048e-8_f * temp**3.0_f &
231 0 : - 0.0579087_f/crit_x
232 :
233 : jcoe = -0.148712_f + 0.00283508_f * temp &
234 : - 9.24619e-6_f * temp**2.0_f &
235 : + 5.00427e-9_f * temp**3.0_f &
236 0 : - 0.0127081_f/crit_x
237 :
238 : tmpa = ( &
239 : acoe &
240 : + bcoe * log (rh) &
241 : + ccoe * ( log (rh))**2.0_f &
242 : + dcoe * ( log (rh))**3.0_f &
243 : + ecoe * log (h2so4) &
244 : + fcoe * (log (rh)) * (log (h2so4)) &
245 : + gcoe * ((log (rh) ) **2.0_f) &
246 : * (log (h2so4)) &
247 : + hcoe * (log (h2so4)) **2.0_f &
248 : + icoe * log (rh) &
249 : * ((log (h2so4)) **2.0_f) &
250 : + jcoe * (log (h2so4)) **3.0_f &
251 0 : )
252 :
253 0 : tmpa = min( tmpa, log(1.0e38_f) )
254 0 : nucrate_cgs = exp ( tmpa )
255 :
256 : ! calc number of molecules in critical cluster
257 : acoe = -0.00295413_f - 0.0976834_f*temp &
258 : + 0.00102485_f * temp**2.0_f &
259 0 : - 2.18646e-6_f * temp**3.0_f - 0.101717_f/crit_x
260 :
261 : ! write(LUNOPRT,*)'291 acoe=',acoe
262 :
263 : bcoe = -0.00205064_f - 0.00758504_f*temp &
264 : + 0.000192654_f * temp**2.0_f &
265 0 : - 6.7043e-7_f * temp**3.0_f - 0.255774_f/crit_x
266 :
267 : ccoe = +0.00322308_f + 0.000852637_f * temp &
268 : - 0.0000154757_f * temp**2.0_f &
269 : + 5.66661e-8_f * temp**3.0_f &
270 0 : + 0.0338444_f/crit_x
271 :
272 : dcoe = +0.0474323_f - 0.000625104_f * temp &
273 : + 2.65066e-6_f * temp**2.0_f &
274 : - 3.67471e-9_f * temp**3.0_f &
275 0 : - 0.000267251_f/crit_x
276 :
277 : ecoe = -0.0125211_f + 0.00580655_f * temp &
278 : - 0.000101674_f * temp**2.0_f &
279 : + 2.88195e-7_f * temp**3.0_f &
280 0 : + 0.0942243_f/crit_x
281 :
282 : fcoe = -0.038546_f - 0.000672316_f * temp &
283 : + 2.60288e-6_f * temp**2.0_f &
284 : + 1.19416e-8_f * temp**3.0_f &
285 0 : - 0.00851515_f/crit_x
286 :
287 : gcoe = -0.0183749_f + 0.000172072_f * temp &
288 : - 3.71766e-7_f * temp**2.0_f &
289 : - 5.14875e-10_f * temp**3.0_f &
290 0 : + 0.00026866_f/crit_x
291 :
292 : hcoe = -0.0619974_f + 0.000906958_f * temp &
293 : - 9.11728e-7_f * temp**2.0_f &
294 : - 5.36796e-9_f * temp**3.0_f &
295 0 : - 0.00774234_f/crit_x
296 :
297 : icoe = +0.0121827_f - 0.00010665_f * temp &
298 : + 2.5346e-7_f * temp**2.0_f &
299 : - 3.63519e-10_f * temp**3.0_f &
300 0 : + 0.000610065_f/crit_x
301 :
302 : jcoe = +0.000320184_f - 0.0000174762_f * temp &
303 : + 6.06504e-8_f * temp**2.0_f &
304 : - 1.4177e-11_f * temp**3.0_f &
305 0 : + 0.000135751_f/crit_x
306 :
307 0 : cnum_tot = acoe + bcoe * log (rh)
308 :
309 : cnum_tot = exp ( &
310 : acoe &
311 : + bcoe * log (rh) &
312 : + ccoe * ( log (rh))**2.0_f &
313 : + dcoe * ( log (rh))**3.0_f &
314 : + ecoe * log (h2so4) &
315 : + fcoe * (log (rh)) * (log (h2so4)) &
316 : + gcoe * ((log (rh) ) **2.0_f) &
317 : * (log (h2so4)) &
318 : + hcoe * (log (h2so4)) **2.0_f &
319 : + icoe * log (rh) &
320 : * ((log (h2so4)) **2.0_f) &
321 : + jcoe * (log (h2so4)) **3.0_f &
322 0 : )
323 :
324 0 : cnum_h2so4 = cnum_tot * crit_x
325 :
326 : ! calc radius (nm) of critical cluster
327 : radius_cluster = exp( -1.6524245_f + 0.42316402_f*crit_x &
328 0 : + 0.3346648_f*log(cnum_tot) )
329 :
330 0 : radius_cluster = radius_cluster * 1e-7_f ! nm -> cm
331 :
332 0 : mass_cluster_dry = cnum_h2so4 * gwtmol(igash2so4) / AVG ! cluster dry mass in g
333 :
334 0 : return
335 0 : end subroutine binary_nuc_vehk2002
336 :
337 : !----------------------------------------------------------------------
338 : !-----------------------------------------------------------------------
339 1418703437 : subroutine binary_nuc_zhao1995( carma, cstate, temp, weight_percent, rh, h2so4, h2so4_cgs, h2o, h2o_cgs, beta1, &
340 : nucrate_cgs, mass_cluster_dry, radius_cluster, ftry, rc )
341 : !! Calculates particle production rates due to nucleation <rhompe>:
342 : !! binary homogeneous nucleation of sulfuric acid and water only
343 : !! Numerical method follows Zhao & Turco, JAS, V.26, No.5, 1995.
344 : use carma_precision_mod
345 : use carma_enums_mod
346 : use carma_constants_mod
347 : use carma_types_mod
348 : use carmastate_mod
349 : use carma_mod
350 : use sulfate_utils
351 :
352 : implicit none
353 : ! subr arguments (in)
354 : type(carma_type), intent(in) :: carma !! the carma object
355 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
356 : real(kind=f), intent(in) :: temp ! temperature (k)
357 : real(kind=f), intent(in) :: weight_percent ! weight percent H2SO4 (0-100)
358 : real(kind=f), intent(in) :: rh ! relative humidity of water wrt liquid water (0-1)
359 : real(kind=f), intent(in) :: h2so4 ! concentration of H2SO4 (molecules cm-3)
360 : real(kind=f), intent(in) :: h2o ! concentration of H2SO4 (molecules cm-3)
361 : real(kind=f), intent(in) :: h2so4_cgs ! H2SO4 densities in g/cm3
362 : real(kind=f), intent(in) :: h2o_cgs ! H2O densities in g/cm3
363 : real(kind=f), intent(in) :: beta1
364 :
365 : ! subr arguments (out)
366 : real(kind=f), intent(out) :: nucrate_cgs ! binary nucleation rate, j (# cm-3 s-1)
367 : real(kind=f), intent(out) :: radius_cluster ! the radius of cluster (cm)
368 : real(kind=f), intent(out) :: mass_cluster_dry ! dry mass of cluster (g)
369 : real(kind=f), intent(out) :: ftry
370 : integer, intent(inout) :: rc !! return code, negative indicates failure
371 :
372 : ! Local declarations
373 : integer :: i, ibin, ie
374 : real(kind=f) :: dens(46)
375 : real(kind=f) :: pa(46)
376 : real(kind=f) :: pb(46)
377 : real(kind=f) :: c1(46)
378 : real(kind=f) :: c2(46)
379 : real(kind=f) :: fct(46)
380 : real(kind=f) :: wtmolr ! molecular weight ration of H2SO4/H2O
381 : real(kind=f) :: h2oln ! H2O ambient vapor pressures [dynes/cm2]
382 : real(kind=f) :: h2so4ln ! H2SO4 ambient vapor pressures [dynes/cm2]
383 : real(kind=f) :: SA ! total surface area of pre-existing wet particles
384 : real(kind=f) :: SAbin ! bin surface area of pre-existing wet particles
385 : real(kind=f) :: cw
386 : real(kind=f) :: dw
387 : real(kind=f) :: wvp ! water eq.vp over solution
388 : real(kind=f) :: wvpln
389 : real(kind=f) :: t0_kulm
390 : real(kind=f) :: seqln
391 : real(kind=f) :: t_crit_kulm
392 : real(kind=f) :: factor_kulm
393 : real(kind=f) :: dw1, dw2
394 : real(kind=f) :: dens1
395 : real(kind=f) :: dens11
396 : real(kind=f) :: dens12
397 : real(kind=f) :: xfrac
398 : real(kind=f) :: wstar
399 : real(kind=f) :: dstar
400 : real(kind=f) :: rhln
401 : real(kind=f) :: raln
402 : real(kind=f) :: wfstar
403 : real(kind=f) :: sigma
404 : real(kind=f) :: ystar
405 : real(kind=f) :: r2
406 : real(kind=f) :: gstar
407 : real(kind=f) :: rpr
408 : real(kind=f) :: rpre
409 : real(kind=f) :: fracmol
410 : real(kind=f) :: zphi
411 : real(kind=f) :: zeld
412 : real(kind=f) :: cfac
413 : real(kind=f) :: ahom
414 : real(kind=f) :: exhom
415 : real(kind=f) :: frac_h2so4
416 : real(kind=f) :: rhomlim
417 : real(kind=f) :: dnpot(46), dnwf(46)
418 : real(kind=f) :: rho_H2SO4_wet
419 :
420 1418703437 : radius_cluster = -1._f
421 :
422 : ! Parameterized fit developed by Mike Mills in 1994 to the partial molal
423 : ! Gibbs energies (F2|o-F2) vs. weight percent H2SO4 table in Giauque et al.,
424 : ! J. Am. Chem. Soc, 82, 62-70, 1960. The parameterization gives excellent
425 : ! agreement. Ayers (GRL, 7, 433-436, 1980) refers to F2|o-F2 as mu - mu_0
426 : ! (chemical potential). This parameterization may be replaced by a lookup
427 : ! table, as was done ultimately in the Garcia-Solomon sulfate code.
428 66679061539 : do i = 1, 46
429 65260358102 : dnpot(i) = 4.184_f * (23624.8_f - 1.14208e8_f / ((dnwtp(i) - 105.318_f)**2 + 4798.69_f))
430 66679061539 : dnwf(i) = dnwtp(i) / 100._f
431 : end do
432 :
433 : ! Molecular weight ratio of H2SO4 / H2O:
434 1418703437 : wtmolr = gwtmol(igash2so4) / gwtmol(igash2o)
435 :
436 : ! Compute ln of H2O and H2SO4 ambient vapor pressures [dynes/cm2]
437 1418703437 : h2oln = log(h2o_cgs * (RGAS / gwtmol(igash2o)) * temp)
438 1418703437 : h2so4ln = log(h2so4_cgs * (RGAS / gwtmol(igash2so4)) * temp)
439 :
440 : ! loop through wt pcts and calculate vp/composition for each
441 66679061539 : do i = 1, 46
442 65260358102 : dens(i) = dnc0(i) + dnc1(i) * temp
443 :
444 : ! Calc. water eq.vp over solution using (Lin & Tabazadeh eqn 5, JGR, 2001)
445 65260358102 : cw = 22.7490_f + 0.0424817_f * dnwtp(i) - 0.0567432_f * dnwtp(i)**0.5_f - 0.000621533_f * dnwtp(i)**2
446 65260358102 : dw = -5850.24_f + 21.9744_f * dnwtp(i) - 44.5210_f * dnwtp(i)**0.5_f - 0.384362_f * dnwtp(i)**2
447 :
448 : ! pH20 | eq[mb]
449 65260358102 : wvp = exp(cw + dw / temp)
450 :
451 : ! Ln(pH2O | eq [dynes/cm2])
452 65260358102 : wvpln = log(wvp * 1013250._f / 1013.25_f)
453 :
454 : ! Save the water eq.vp over solution at each wt pct into this array:
455 : !
456 : ! Ln(pH2O/pH2O|eq) with both terms in dynes/cm2
457 65260358102 : pb(i) = h2oln - wvpln
458 :
459 : ! Calc. sulfuric acid eq.vp over solution using (Ayers et. al., GRL, V.7, No.6, June 1980)
460 : !
461 : ! T0 set in the low end of the Ayers measurement range (338-445K)
462 65260358102 : t0_kulm = 340._f
463 65260358102 : seqln = -10156._f / t0_kulm + 16.259_f
464 :
465 : ! Now calc. Kulmala correction (J. CHEM. PHYS. V.93, No.1, 1 July 1990)
466 : !
467 : ! Critical temperature = 1.5 * Boiling point
468 65260358102 : t_crit_kulm = 905._f
469 : factor_kulm = -1._f / temp + 1._f / t0_kulm + 0.38_f / (t_crit_kulm - t0_kulm) * &
470 65260358102 : (1.0_f + log(t0_kulm / temp) - t0_kulm / temp)
471 :
472 : ! For pure sulfuric acid
473 65260358102 : seqln = seqln + 10156._f * factor_kulm
474 :
475 : ! Now adjust vp based on weight % composition using parameterization of Giauque 1960
476 : !
477 : ! Adjust for weight percent composition
478 65260358102 : seqln = seqln - dnpot(i) / (8.3143_f * temp)
479 :
480 : ! Convert atmospheres => dynes/cm2
481 65260358102 : seqln = seqln + log(1013250._f)
482 :
483 : ! Save the sulfuric acid eq.vp over solution at each wt pct into this array:
484 : !
485 : ! Ln(pH2SO4/pH2SO4|eq) with both terms in dynes/cm2
486 65260358102 : pa(i) = h2so4ln - seqln
487 :
488 : ! Create 2-component solutions of varying composition c1 and c2
489 65260358102 : c1(i) = pa(i) - pb(i) * wtmolr
490 66679061539 : c2(i) = pa(i) * dnwf(i) + pb(i) * (1._f - dnwf(i)) * wtmolr
491 : end do ! end of loop through weight percents
492 :
493 : ! Now loop through until we find the c1+c2 combination with minimum Gibbs free energy
494 1418703437 : dw2 = dnwtp(46) - dnwtp(45)
495 1418703437 : dens1 = (dens(46) - dens(45)) / dw2
496 1418703437 : fct(46) = c1(46) + c2(46) * 100._f * dens1 / dens(46)
497 1418703437 : dens12 = dens1
498 :
499 49652936713 : do i = 45, 2, -1
500 49598666834 : dw1 = dw2
501 49598666834 : dens11 = dens12
502 49598666834 : dw2 = dnwtp(i) - dnwtp(i-1)
503 49598666834 : dens12 = (dens(i) - dens(i-1)) / dw2
504 49598666834 : dens1 = (dens11 * dw2 + dens12 * dw1) / (dw1 + dw2)
505 :
506 49598666834 : fct(i) = c1(i) + c2(i) * 100._f * dens1 / dens(i)
507 :
508 : ! Find saddle where fct(i)<0<fct(i+1)
509 49652936713 : if (fct(i) * fct(i+1) <= 0._f) exit
510 : end do
511 :
512 1418703437 : if (i == 1) then
513 54269879 : dens1 = (dens(2) - dens(1)) / (dnwtp(2) - dnwtp(1))
514 54269879 : fct(1) = c1(1) + c2(1) * 100._f * dens1 / dens(1)
515 : end if
516 :
517 : ! Possibility 1: loop finds no saddle, so no nucleation occurs:
518 1418703437 : if (fct(i) * fct(i+1) > 0._f) then
519 44570862 : nucrate_cgs = 0.0_f
520 44570862 : radius_cluster = 0.0_f
521 44570862 : mass_cluster_dry = 0.0_f
522 44570862 : ftry = 0.0_f
523 44570862 : return
524 :
525 : ! Possibility 2: loop crossed the saddle; interpolate to find exact value:
526 1374132575 : else if (fct(i) * fct(i+1) < 0._f) then
527 1374132575 : xfrac = fct(i+1) / (fct(i+1) - fct(i))
528 1374132575 : wstar = dnwtp(i+1) * (1.0_f - xfrac) + dnwtp(i) * xfrac ! critical weight percent
529 1374132575 : dstar = dens(i+1) * (1.0_f - xfrac) + dens(i) * xfrac
530 1374132575 : rhln = pb(i+1) * (1.0_f - xfrac) + pb(i) * xfrac
531 1374132575 : raln = pa(i+1) * (1.0_f - xfrac) + pa(i) * xfrac
532 :
533 : ! Possibility 3: loop found the saddle point exactly
534 : else
535 0 : dstar = dens(i)
536 :
537 : ! critical weight percent
538 0 : wstar = dnwtp(i)
539 0 : rhln = pb(i)
540 0 : raln = pa(i)
541 : end if
542 :
543 : ! Critical weight fraction
544 1374132575 : wfstar = wstar / 100._f
545 :
546 1374132575 : if ((wfstar < 0._f) .or. (wfstar > 1._f)) then
547 0 : write(LUNOPRT,*)'sulfnuc: wstar out of bounds!'
548 0 : rc = RC_ERROR
549 0 : return
550 : end if
551 :
552 : ! Critical surface tension [erg/cm2]
553 1374132575 : sigma = sulfate_surf_tens(carma, wstar, temp, rc)
554 :
555 : ! Critical Y (eqn 13 in Zhao & Turco 1993) [erg/cm3]
556 : ystar = dstar * RGAS * temp * (wfstar / gwtmol(igash2so4) &
557 1374132575 : * raln + (1._f - wfstar) / gwtmol(igash2o) * rhln)
558 1374132575 : if (ystar < 1.e-20_f) then
559 433249482 : nucrate_cgs = 0.0_f
560 433249482 : radius_cluster = 0.0_f
561 433249482 : mass_cluster_dry = 0.0_f
562 433249482 : ftry = 0.0_f
563 433249482 : return
564 : end if
565 :
566 : ! Critical cluster radius [cm]
567 940883093 : radius_cluster = 2._f * sigma / ystar
568 940883093 : radius_cluster = max(radius_cluster, 0.0_f)
569 940883093 : r2 = radius_cluster * radius_cluster
570 :
571 : ! Critical Gibbs free energy [erg]
572 940883093 : gstar = (4._f * PI / 3._f) * r2 * sigma
573 :
574 : ! RPR[molecules/s] = 4Pi * R2[cm2] * H2O[molecules/cm3] * Beta[cm/s]
575 940883093 : rpr = 4._f * PI * r2 * h2o * beta1
576 :
577 : ! RPRE[/cm3/s] = RPR[/s] * H2SO4[/cm3]; first part of Zhao & Turco eqn 16
578 940883093 : rpre = rpr * h2so4
579 :
580 : ! Zeldovitch non-equilibrium correction factor [unitless]
581 : ! Jaecker-Voirol & Mirabel, 1988 (not considered in Zhao & Turco)
582 940883093 : fracmol = 1._f /(1._f + wtmolr * (1._f - wfstar) / wfstar)
583 : zphi = atan(fracmol)
584 940883093 : zeld = 0.25_f / (sin(zphi))**2
585 :
586 : ! Empirical correction factor:
587 940883093 : cfac = 0.0_f
588 :
589 : ! Gstar exponential term in Zhao & Turco eqn 16 [unitless]
590 940883093 : ftry = (-gstar / BK / temp)
591 940883093 : ahom = ftry + cfac
592 940883093 : if (ahom .lt. -500._f) then
593 : exhom=0.0_f
594 : else
595 652215392 : exhom = exp(min(ahom, 28.0_f))
596 : endif
597 :
598 : ! Calculate mass of critical nucleus
599 940883093 : rho_H2SO4_wet = sulfate_density(carma, weight_percent, temp, rc)
600 940883093 : mass_cluster_dry = (4._f * PI / 3._f) * rho_H2SO4_wet * r2 * radius_cluster
601 :
602 : ! Calculate dry mass of critical nucleus
603 940883093 : mass_cluster_dry = mass_cluster_dry * wfstar
604 :
605 : ! Calculate the nucleation rate [#/cm3/s], Zhao & Turco eqn 16.
606 940883093 : nucrate_cgs = rpre * zeld * exhom
607 :
608 940883093 : return
609 1418703437 : end subroutine binary_nuc_zhao1995
|
