| blob | 1a503c168b77e0d0ea4de2ae70a66c536c88d048 |
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************
9 module module_mosaic_newnuc
13 use module_peg_util
17 implicit none
21 contains
25 !-----------------------------------------------------------------------
26 subroutine mosaic_newnuc_1clm( istat_newnuc, &
27 it, jt, kclm_calcbgn, kclm_calcend, &
28 idiagbb_in, &
29 dtchem, dtnuc_in, rsub0, &
30 id, ktau, ktauc, its, ite, jts, jte, kts, kte )
31 !
32 ! calculates new particle nucleation for grid points in
33 ! the i=it, j=jt vertical column over timestep dtnuc_in
34 ! works with mosaic sectional aerosol packages
35 !
36 ! when newnuc_method = 1 (internal parameter), uses
37 ! h2so4-nh3-h2o ternary nucleation from napari et al., 2002
38 ! *** this option is currently disabled because the ternary
39 ! parameterization was recently shown to be invalid
40 ! when newnuc_method = 2 (internal parameter), uses
41 ! h2so4-h2o binary nucleation from wexler et al., 1994
42 !
43 use module_data_mosaic_asect
44 use module_data_mosaic_other
45 use module_state_description, only: param_first_scalar
47 ! subr arguments
48 integer, intent(inout) :: istat_newnuc ! =0 if no problems
49 integer, intent(in) :: &
50 it, jt, kclm_calcbgn, kclm_calcend, &
51 idiagbb_in, &
52 id, ktau, ktauc, its, ite, jts, jte, kts, kte
53 real, intent(in) :: dtchem, dtnuc_in
54 real, intent(in) :: rsub0(l2maxd,kmaxd,nsubareamaxd)
55 ! rsub0 holds mixrat values before the aerchemistry calcs
57 ! NOTE - much information is passed via arrays in
58 ! module_data_mosaic_asect and module_data_mosaic_other
59 !
60 ! rsub (inout) - trace gas and aerosol mixing ratios
61 ! aqmassdry_sub, aqvoldry_sub (inout) -
62 ! aerosol dry-mass and dry-volume mixing ratios
63 ! adrydens_sub (inout) - aerosol dry density
64 ! rsub(ktemp,:,:), cairclm, relhumclm (in) -
65 ! air temperature, molar density, and relative humidity
68 ! local variables
69 integer, parameter :: newnuc_method = 2
71 integer :: k, l, ll, m
72 integer :: isize, itype, iphase
73 integer :: iconform_numb
74 integer :: idiagbb
75 integer :: nsize, ntau_nuc
76 integer :: ncount(10)
77 integer :: p1st
79 real, parameter :: densdefault = 2.0
80 real, parameter :: smallmassbb = 1.0e-30
82 ! nh4hso4 values for a_zsr and b_zsr
83 real, parameter :: a_zsr_xx1 = 1.15510
84 real, parameter :: a_zsr_xx2 = -3.20815
85 real, parameter :: a_zsr_xx3 = 2.71141
86 real, parameter :: a_zsr_xx4 = 2.01155
87 real, parameter :: a_zsr_xx5 = -4.71014
88 real, parameter :: a_zsr_xx6 = 2.04616
89 real, parameter :: b_zsr_xx = 29.4779
91 real :: aw
92 real :: cair_box
93 real :: dens_nh4so4a, dtnuc
94 real :: duma, dumb, dumc
95 real :: rh_box
96 real :: qh2so4_avg, qh2so4_cur, qh2so4_del
97 real :: qnh3_avg, qnh3_cur, qnh3_del
98 real :: qnuma_del, qso4a_del, qnh4a_del
99 real :: temp_box
100 real :: xxdens, xxmass, xxnumb, xxvolu
102 real,save :: dumveca(10), dumvecb(10), dumvecc(10), dumvecd(10), dumvece(10)
103 real :: volumlo_nuc(maxd_asize), volumhi_nuc(maxd_asize)
105 character(len=100) :: msg
109 ! check newnuc_method
110 istat_newnuc = 0
111 if (newnuc_method .ne. 2) then
112 if ((it .eq. its) .and. (jt .eq. jts)) &
113 call peg_message( lunerr, &
114 '*** mosaic_newnuc_1clm -- illegal newnuc_method' )
115 istat_newnuc = -1
116 return
117 end if
120 ! when dtnuc_in > dtchem, do not perform nucleation calcs
121 ! on every chemistry time step
122 ntau_nuc = nint( dtnuc_in/dtchem )
123 ntau_nuc = max( 1, ntau_nuc )
124 if (mod(ktau,ntau_nuc) .ne. 0) return
125 dtnuc = dtchem*ntau_nuc
128 ! set variables that do not change
129 idiagbb = idiagbb_in
131 itype = 1
132 iphase = ai_phase
133 nsize = nsize_aer(itype)
134 volumlo_nuc(1:nsize) = volumlo_sect(1:nsize,itype)
135 volumhi_nuc(1:nsize) = volumhi_sect(1:nsize,itype)
138 ! loop over subareas (currently only 1) and vertical levels
139 do 2900 m = 1, nsubareas
141 do 2800 k = kclm_calcbgn, kclm_calcend
144 ! initialize diagnostics
145 if ((it .eq. its) .and. &
146 (jt .eq. jts) .and. (k .eq. kclm_calcbgn)) then
147 dumveca(:) = 0.0 ! current grid param values
148 dumvecb(:) = +1.0e35 ! param minimums
149 dumvecc(:) = -1.0e35 ! param maximums
150 dumvecd(:) = 0.0 ! param averages
151 dumvece(:) = 0.0 ! param values for highest qnuma_del
152 ncount(:) = 0
153 end if
156 ncount(1) = ncount(1) + 1
157 if (afracsubarea(k,m) .lt. 1.e-4) goto 2700
159 cair_box = cairclm(k)
160 temp_box = rsub(ktemp,k,m)
161 rh_box = relhumclm(k)
163 qh2so4_cur = max(0.0,rsub(kh2so4,k,m))
164 qnh3_cur = max(0.0,rsub(knh3,k,m))
165 qh2so4_avg = 0.5*( qh2so4_cur + max(0.0,rsub0(kh2so4,k,m)) )
166 qnh3_avg = 0.5*( qnh3_cur + max(0.0,rsub0(knh3,k,m)) )
168 qh2so4_del = 0.0
169 qnh3_del = 0.0
170 qnuma_del = 0.0
171 qso4a_del = 0.0
172 qnh4a_del = 0.0
174 dens_nh4so4a = dens_so4_aer
176 isize = 0
178 ! make call to nucleation routine
179 if (newnuc_method .eq. 1) then
180 call ternary_nuc_mosaic_1box( &
181 dtnuc, temp_box, rh_box, cair_box, &
182 qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
183 nsize, maxd_asize, volumlo_nuc, volumhi_nuc, &
184 isize, qnuma_del, qso4a_del, qnh4a_del, &
185 qh2so4_del, qnh3_del, dens_nh4so4a )
186 else if (newnuc_method .eq. 2) then
187 call wexler_nuc_mosaic_1box( &
188 dtnuc, temp_box, rh_box, cair_box, &
189 qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
190 nsize, maxd_asize, volumlo_nuc, volumhi_nuc, &
191 isize, qnuma_del, qso4a_del, qnh4a_del, &
192 qh2so4_del, qnh3_del, dens_nh4so4a )
193 else
194 istat_newnuc = -1
195 return
196 end if
199 ! temporary diagnostics
200 dumveca(1) = temp_box
201 dumveca(2) = rh_box
202 dumveca(3) = rsub(kso2,k,m)
203 dumveca(4) = qh2so4_avg
204 dumveca(5) = qnh3_avg
205 dumveca(6) = qnuma_del
206 do l = 1, 6
207 dumvecb(l) = min( dumvecb(l), dumveca(l) )
208 dumvecc(l) = max( dumvecc(l), dumveca(l) )
209 dumvecd(l) = dumvecd(l) + dumveca(l)
210 if (qnuma_del .gt. dumvece(6)) dumvece(l) = dumveca(l)
211 end do
214 ! check for zero new particles
215 if (qnuma_del .le. 0.0) goto 2700
217 ! check for valid isize
218 if (isize .ne. 1) ncount(3) = ncount(3) + 1
219 if ((isize .lt. 1) .or. (isize .gt. nsize)) then
220 write(msg,93010) 'newnucxx bad isize_nuc' , it, jt, k, &
221 isize, nsize
222 call peg_message( lunerr, msg )
223 goto 2700
224 end if
225 93010 format( a, 3i3, 1p, 9e10.2 )
228 ncount(2) = ncount(2) + 1
230 ! update gas and aerosol so4 and nh3/nh4 mixing ratios
231 rsub(kh2so4,k,m) = max( 0.0, rsub(kh2so4,k,m) + qh2so4_del )
232 rsub(knh3, k,m) = max( 0.0, rsub(knh3, k,m) + qnh3_del )
234 l = lptr_so4_aer(isize,itype,iphase)
235 if (l .ge. p1st) then
236 rsub(l,k,m) = rsub(l,k,m) + qso4a_del
237 end if
238 l = lptr_nh4_aer(isize,itype,iphase)
239 if (l .ge. p1st) then
240 rsub(l,k,m) = rsub(l,k,m) + qnh4a_del
241 end if
242 l = numptr_aer(isize,itype,iphase)
243 rsub(l,k,m) = rsub(l,k,m) + qnuma_del
244 xxnumb = rsub(l,k,m)
246 ! update aerosol water, using mosaic parameterizations for nh4hso4
247 ! duma = (mole-salt)/(mole-salt+water)
248 ! dumb = (mole-salt)/(kg-water)
249 ! dumc = (mole-water)/(mole-salt)
250 l = waterptr_aer(isize,itype)
251 if ((rh_box .gt. 0.10) .and. (l .ge. p1st)) then
252 aw = min( rh_box, 0.98 )
253 if (aw .lt. 0.97) then
254 duma = a_zsr_xx1 + &
255 aw*( a_zsr_xx2 + &
256 aw*( a_zsr_xx3 + &
257 aw*( a_zsr_xx4 + &
258 aw*( a_zsr_xx5 + &
259 aw* a_zsr_xx6 ))))
260 else
261 dumb = -b_zsr_xx*log(aw)
262 dumb = max( dumb, 0.5 )
263 duma = 1.0/(1.0 + 55.509/dumb)
264 end if
265 duma = max( duma, 0.01 )
266 dumc = (1.0 - duma)/duma
267 rsub(l,k,m) = rsub(l,k,m) + qso4a_del*dumc
268 end if
271 !
272 ! update dry mass, density, and volume,
273 ! and check for mean dry-size within bounds
274 !
275 xxmass = aqmassdry_sub(isize,itype,k,m)
276 xxdens = adrydens_sub( isize,itype,k,m)
277 iconform_numb = 1
279 if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
280 ! (exception) case of drydensity not valid
281 continue
282 else
283 ! (normal) case of drydensity valid (which means drymass is valid also)
284 ! so increment mass and volume with the so4 & nh4 deltas, then calc density
285 xxvolu = xxmass/xxdens
286 duma = qso4a_del*mw_so4_aer + qnh4a_del*mw_nh4_aer
287 xxmass = xxmass + duma
288 xxvolu = xxvolu + duma/dens_nh4so4a
289 if (xxmass .le. smallmassbb) then
290 ! do this to force calc of dry mass, volume from rsub
291 xxdens = 0.001
292 else if (xxmass .gt. 1000.0*xxvolu) then
293 ! in this case, density is too large. setting density=1000 forces
294 ! next IF block while avoiding potential divide by zero or overflow
295 xxdens = 1000.0
296 else
297 xxdens = xxmass/xxvolu
298 end if
299 end if
301 if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
302 ! (exception) case of drydensity not valid (or drymass extremely small),
303 ! so compute from dry mass, volume from rsub
304 ncount(4) = ncount(4) + 1
305 xxmass = 0.0
306 xxvolu = 0.0
307 do ll = 1, ncomp_aer(itype)
308 l = massptr_aer(ll,isize,itype,iphase)
309 if (l .ge. p1st) then
310 duma = max( 0.0, rsub(l,k,m) )*mw_aer(ll,itype)
311 xxmass = xxmass + duma
312 xxvolu = xxvolu + duma/dens_aer(ll,itype)
313 end if
314 end do
315 end if
317 if (xxmass .le. smallmassbb) then
318 ! when drymass extremely small, use default density and bin center size,
319 ! and zero out water
320 ncount(5) = ncount(5) + 1
321 xxdens = densdefault
322 xxvolu = xxmass/xxdens
323 xxnumb = xxmass/(volumcen_sect(isize,itype)*xxdens)
324 iconform_numb = 0
325 l = waterptr_aer(isize,itype)
326 if (l .ge. p1st) rsub(l,k,m) = 0.0
327 l = hyswptr_aer(isize,itype)
328 if (l .ge. p1st) rsub(l,k,m) = 0.0
329 else
330 xxdens = xxmass/xxvolu
331 end if
333 if (iconform_numb .gt. 0) then
334 ! check for mean dry-size within bounds, and conform number if not
335 if (xxnumb .gt. xxvolu/volumlo_sect(isize,itype)) then
336 ncount(6) = ncount(6) + 1
337 xxnumb = xxvolu/volumlo_sect(isize,itype)
338 else if (xxnumb .lt. xxvolu/volumhi_sect(isize,itype)) then
339 ncount(7) = ncount(7) + 1
340 xxnumb = xxvolu/volumhi_sect(isize,itype)
341 end if
342 end if
344 ! load dry mass, density, volume, and (possibly conformed) number
345 l = numptr_aer(isize,itype,iphase)
346 rsub(l,k,m) = xxnumb
347 adrydens_sub( isize,itype,k,m) = xxdens
348 aqmassdry_sub(isize,itype,k,m) = xxmass
349 aqvoldry_sub( isize,itype,k,m) = xxvolu
352 2700 continue
354 ! temporary diagnostics
355 if ((idiagbb .ge. 100) .and. &
356 (it .eq. ite) .and. &
357 (jt .eq. jte) .and. (k .eq. kclm_calcend)) then
358 if (idiagbb .ge. 110) then
359 write(msg,93020) 'newnucbb mins ', dumvecb(1:6)
360 call peg_message( lunerr, msg )
361 write(msg,93020) 'newnucbb maxs ', dumvecc(1:6)
362 call peg_message( lunerr, msg )
363 duma = max( 1, ncount(1) )
364 write(msg,93020) 'newnucbb avgs ', dumvecd(1:6)/duma
365 call peg_message( lunerr, msg )
366 write(msg,93020) 'newnucbb hinuc', dumvece(1:6)
367 call peg_message( lunerr, msg )
368 write(msg,93020) 'newnucbb dtnuc', dtnuc
369 call peg_message( lunerr, msg )
370 end if
371 write(msg,93030) 'newnucbb ncnt ', ncount(1:7)
372 call peg_message( lunerr, msg )
373 end if
374 93020 format( a, 1p, 10e10.2 )
375 93030 format( a, 1p, 10i10 )
378 2800 continue ! k levels
380 2900 continue ! subareas
383 return
384 end subroutine mosaic_newnuc_1clm
387 !-----------------------------------------------------------------------
388 ! note - subrs ternary_nuc_mosaic_1box, ternary_nuc_napari, wexler_nuc_mosaic_1box
389 ! taken from ternucl04.f90 on 22-jul-2006
390 !-----------------------------------------------------------------------
391 subroutine ternary_nuc_mosaic_1box( &
392 dtnuc, temp_in, rh_in, cair, &
393 qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
394 nsize, maxd_asize, volumlo_sect, volumhi_sect, &
395 isize_nuc, qnuma_del, qso4a_del, qnh4a_del, &
396 qh2so4_del, qnh3_del, dens_nh4so4a )
397 !.......................................................................
398 !
399 ! calculates new particle production from h2so4-nh3-h2o ternary nucleation
400 ! over timestep dtnuc, using nucleation rates from the
401 ! napari et al. (2002) parameterization
402 !
403 ! the new particles are "grown" to the lower-bound size of the host code's
404 ! smallest size bin. (this "growth" is purely ad hoc, and would not be
405 ! necessary if the host code's size bins extend down to ~1 nm.)
406 ! if the h2so4 and nh3 mass mixing ratios (mixrats) of the grown new
407 ! particles exceed the current gas mixrats, the new particle production
408 ! is reduced so that the new particle mass mixrats match the gas mixrats.
409 !
410 ! revision history
411 ! coded by rc easter, pnnl, 20-mar-2006
412 !
413 ! key routines called: subr ternary_nuc_napari
414 !
415 ! references:
416 ! napari, i., m. noppel, h. vehkamaki, and m. kulmala,
417 ! parameterization of ternary nucleation rates for
418 ! h2so4-nh3-h2o vapors. j. geophys. res., 107, 4381, 2002.
419 !
420 !.......................................................................
421 implicit none
423 ! subr arguments (in)
424 real, intent(in) :: dtnuc ! nucleation time step (s)
425 real, intent(in) :: temp_in ! temperature, in k
426 real, intent(in) :: rh_in ! relative humidity, as fraction
427 real, intent(in) :: cair ! dry-air molar density (mole/cm3)
429 real, intent(in) :: qh2so4_avg, qh2so4_cur ! gas h2so4 mixing ratios (mole/mole-air)
430 real, intent(in) :: qnh3_avg, qnh3_cur ! gas nh3 mixing ratios (mole/mole-air)
431 ! qxxx_cur = current value (at end of condensation)
432 ! qxxx_avg = average value (from start to end of condensation)
434 integer, intent(in) :: nsize ! number of aerosol size bins
435 integer, intent(in) :: maxd_asize ! dimension for volumlo_sect, ...
436 real, intent(in) :: volumlo_sect(maxd_asize) ! dry volume at lower bnd of bin (cm3)
437 real, intent(in) :: volumhi_sect(maxd_asize) ! dry volume at upper bnd of bin (cm3)
439 ! subr arguments (out)
440 integer, intent(out) :: isize_nuc ! size bin into which new particles go
441 real, intent(out) :: qnuma_del ! change to aerosol number mixing ratio (#/mole-air)
442 real, intent(out) :: qso4a_del ! change to aerosol so4 mixing ratio (mole/mole-air)
443 real, intent(out) :: qnh4a_del ! change to aerosol nh4 mixing ratio (mole/mole-air)
444 real, intent(out) :: qh2so4_del ! change to gas h2so4 mixing ratio (mole/mole-air)
445 real, intent(out) :: qnh3_del ! change to gas nh3 mixing ratio (mole/mole-air)
446 ! aerosol changes are > 0; gas changes are < 0
448 ! subr arguments (inout)
449 real, intent(inout) :: dens_nh4so4a ! dry-density of the new nh4-so4 aerosol mass (g/cm3)
450 ! use 'in' value only if it is between 1.6-2.0 g/cm3
452 ! subr arguments (out) passed via common block
453 ! these are used to duplicate the outputs of yang zhang's original test driver
454 ! they are not really needed in wrf-chem, so just make them local
455 real :: ratenuclt ! j = ternary nucleation rate from napari param. (cm-3 s-1)
456 real :: rateloge ! ln (j)
457 real :: cnum_h2so4 ! number of h2so4 molecules in the critical nucleus
458 real :: cnum_nh3 ! number of nh3 molecules in the critical nucleus
459 real :: cnum_tot ! total number of molecules in the critical nucleus
460 real :: radius_cluster ! the radius of cluster (nm)
462 ! common / ternary_nuc_yzhang_cmn01 / &
463 ! ratenuclt, rateloge, &
464 ! cnum_h2so4, cnum_nh3, cnum_tot, radius_cluster
467 ! local variables
468 integer i
469 integer, save :: icase = 0, icase_reldiffmax = 0
471 real, parameter :: pi = 3.1415926536
472 real, parameter :: avogad = 6.022e23 ! avogadro number (molecules/mole)
473 real, parameter :: mw_air = 28.966 ! dry-air mean molecular weight (g/mole)
475 ! dry densities (g/cm3) molecular weights of aerosol
476 ! ammsulf, ammbisulf, and sulfacid (from mosaic dens_electrolyte values)
477 real, parameter :: dens_ammsulf = 1.769
478 real, parameter :: dens_ammbisulf = 1.78
479 real, parameter :: dens_sulfacid = 1.841
481 ! molecular weights (g/mole) of aerosol ammsulf, ammbisulf, and sulfacid
482 ! for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
483 ! because we don't keep track of aerosol hion mass
484 real, parameter :: mw_ammsulf = 132.0
485 real, parameter :: mw_ammbisulf = 114.0
486 real, parameter :: mw_sulfacid = 96.0
487 ! molecular weights of aerosol sulfate and ammonium
488 real, parameter :: mw_so4a = 96.0
489 real, parameter :: mw_nh4a = 18.0
491 real, save :: reldiffmax = 0.0
493 real dens_part ! "grown" single-particle dry density (g/cm3)
494 real duma, dumb, dumc, dume
495 real dum_m1, dum_m2, dum_m3, dum_n1, dum_n2, dum_n3
496 real fogas, foso4a, fonh4a, fonuma
497 real freduce ! reduction factor applied to nucleation rate
498 ! due to limited availability of h2so4 & nh3 gases
499 real freducea, freduceb
500 real gramaero_per_moleso4a ! (g dry aerosol)/(mole aerosol so4)
501 real mass_part ! "grown" single-particle mass (g)
502 real molenh4a_per_moleso4a ! (mole aerosol nh4)/(mole aerosol so4)
503 real nh3conc_in ! mixing ratio of nh3 for nucl. calc., pptv
504 real so4vol_in ! concentration of h2so4 for nucl. calc., molecules cm-3
505 real qmolnh4a_del_max ! max production of aerosol nh4 over dtnuc (mole/mole-air)
506 real qmolso4a_del_max ! max production of aerosol so4 over dtnuc (mole/mole-air)
507 real vol_cluster ! critical-cluster volume (cm3)
508 real vol_part ! "grown" single-particle volume (cm3)
511 !
512 ! if h2so4 vapor < 4.0e-16 mole/moleair ~= 1.0e4 molecules/cm3,
513 ! exit with new particle formation = 0
514 !
515 isize_nuc = 1
516 qnuma_del = 0.0
517 qso4a_del = 0.0
518 qnh4a_del = 0.0
519 qh2so4_del = 0.0
520 qnh3_del = 0.0
521 if (qh2so4_avg .le. 4.0e-16) return
522 if (qh2so4_cur .le. 4.0e-16) return
525 !
526 ! make call to napari parameterization routine
527 !
529 ! convert nh3 from mole/mole-air to ppt
530 ! qnh3_cur = nh3conc(m)*1.0e-12
531 nh3conc_in = qnh3_avg/1.0e-12
533 ! convert h2so4 from mole/mole-air to molecules/cm3
534 ! qh2so4_cur = (so4vol(k)/avogad) / cair
535 so4vol_in = (qh2so4_avg) * cair * avogad
537 call ternary_nuc_napari( &
538 temp_in, rh_in, nh3conc_in, so4vol_in, &
539 ratenuclt, rateloge, &
540 cnum_h2so4, cnum_nh3, cnum_tot, radius_cluster )
542 ! if nucleation rate is less than 0.1 particle/cm3/day,
543 ! exit with new particle formation = 0
544 if (ratenuclt .le. 1.0e-6) return
547 ! determine size bin into which the new particles go
548 ! (probably it will always be bin #1, but ...)
549 vol_cluster = (pi*4.0/3.0)* (radius_cluster**3) * 1.0e-21
550 isize_nuc = 1
551 vol_part = volumlo_sect(1)
552 if (vol_cluster .le. volumlo_sect(1)) then
553 continue
554 else if (vol_cluster .ge. volumhi_sect(nsize)) then
555 isize_nuc = nsize
556 vol_part = volumhi_sect(nsize)
557 else
558 do i = 1, nsize
559 if (vol_cluster .lt. volumhi_sect(i)) then
560 isize_nuc = i
561 vol_part = vol_cluster
562 vol_part = min( vol_part, volumhi_sect(i) )
563 vol_part = max( vol_part, volumlo_sect(i) )
564 exit
565 end if
566 end do
567 end if
570 !
571 ! determine composition and density of the "grown particles"
572 ! the grown particles are assumed to be liquid
573 ! (since critical clusters contain water)
574 ! so any (nh4/so4) molar ratio between 0 and 2 is allowed
575 ! assume that the grown particles will have
576 ! (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
577 !
578 if (qnh3_cur .ge. qh2so4_cur) then
579 ! combination of ammonium sulfate and ammonium bisulfate
580 ! dum_n1 & dum_n2 = mole fractions of the ammsulf & ammbisulf
581 dum_n1 = (qnh3_cur/qh2so4_cur) - 1.0
582 dum_n1 = max( 0.0, min( 1.0, dum_n1 ) )
583 dum_n2 = 1.0 - dum_n1
584 dum_n3 = 0.0
585 else
586 ! combination of ammonium bisulfate and sulfuric acid
587 ! dum_n2 & dum_n3 = mole fractions of the ammbisulf & sulfacid
588 dum_n1 = 0.0
589 dum_n2 = (qnh3_cur/qh2so4_cur)
590 dum_n2 = max( 0.0, min( 1.0, dum_n2 ) )
591 dum_n3 = 1.0 - dum_n2
592 end if
594 dum_m1 = dum_n1*mw_ammsulf
595 dum_m2 = dum_n2*mw_ammbisulf
596 dum_m3 = dum_n3*mw_sulfacid
597 dens_part = (dum_m1 + dum_m2 + dum_m3)/ &
598 ((dum_m1/dens_ammsulf) + (dum_m2/dens_ammbisulf) &
599 + (dum_m3/dens_sulfacid))
600 ! 25-jul-2006 - use 'in' value only if it is between 1.6-2.0 g/cm3
601 if (abs(dens_nh4so4a-1.8) .le. 0.2) then
602 dens_part = dens_nh4so4a
603 else
604 dens_nh4so4a = dens_part
605 end if
606 mass_part = vol_part*dens_part
607 molenh4a_per_moleso4a = 2.0*dum_n1 + dum_n2 ! (mole aerosol nh4)/(mole aerosol so4)
608 gramaero_per_moleso4a = dum_m1 + dum_m2 + dum_m3 ! (g dry aerosol)/(mole aerosol so4)
611 ! max production of aerosol dry mass (g-aero/cm3-air)
612 duma = max( 0.0, (ratenuclt*dtnuc*mass_part) )
613 ! max production of aerosol so4 (mole-so4a/cm3-air)
614 dumc = duma/gramaero_per_moleso4a
615 ! max production of aerosol so4 (mole-so4a/mole-air)
616 dume = dumc/cair
617 ! max production of aerosol so4 (mole/mole-air)
618 ! based on ratenuclt and mass_part
619 qmolso4a_del_max = dume
621 ! check if max production exceeds available h2so4 vapor
622 freducea = 1.0
623 if (qmolso4a_del_max .gt. qh2so4_cur) then
624 freducea = qh2so4_cur/qmolso4a_del_max
625 end if
627 ! check if max production exceeds available nh3 vapor
628 freduceb = 1.0
629 if (molenh4a_per_moleso4a .ge. 1.0e-10) then
630 ! max production of aerosol nh4 (ppm) based on ratenuclt and mass_part
631 qmolnh4a_del_max = qmolso4a_del_max*molenh4a_per_moleso4a
632 if (qmolnh4a_del_max .gt. qnh3_cur) then
633 freduceb = qnh3_cur/qmolnh4a_del_max
634 end if
635 end if
636 freduce = min( freducea, freduceb )
638 ! if adjusted nucleation rate is less than 0.1 particle/cm3/day,
639 ! exit with new particle formation = 0
640 if (freduce*ratenuclt .le. 1.0e-6) return
643 ! note: suppose that at this point, freduce = 1.0 (no gas-available
644 ! constraints) and molenh4a_per_moleso4a < 2.0
645 ! then it would be possible to condense "additional" nh3 and have
646 ! (nh3/h2so4 gas molar ratio) < (nh4/so4 aerosol molar ratio) <= 2
647 ! one could do some additional calculations of
648 ! dens_part & molenh4a_per_moleso4a to realize this
649 ! however, the particle "growing" is a crude approximate way to get
650 ! the new particles to the host code's minimum particle size,
651 ! so are such refinements worth the effort?
654 ! changes to h2so4 & nh3 gas (in mole/mole-air), limited by amounts available
655 duma = 0.9999
656 qh2so4_del = min( duma*qh2so4_cur, freduce*qmolso4a_del_max )
657 qnh3_del = min( duma*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
658 qh2so4_del = -qh2so4_del
659 qnh3_del = -qnh3_del
661 ! changes to so4 & nh4 aerosol (in mole/mole-air)
662 qso4a_del = -qh2so4_del
663 qnh4a_del = -qnh3_del
664 ! change to aerosol number (in #/mole-air)
665 qnuma_del = (qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
668 return
669 end subroutine ternary_nuc_mosaic_1box
673 !-----------------------------------------------------------------------
674 subroutine ternary_nuc_napari( &
675 temp_in, rh_in, nh3conc_in, so4vol_in, &
676 ratenuclt, rateloge, &
677 cnum_h2so4, cnum_nh3, cnum_tot, radius_cluster )
678 !*********************************************************************
679 ! purpose: calculate ternary nucleation rates based on *
680 ! napari et al. (2002) parameterization *
681 ! *
682 ! revision history: *
683 ! coded by yang zhang, ncsu, nov. 25, 2004 *
684 ! *
685 ! 14-mar-2006 rce - yang sent this on 10-mar-2005; *
686 ! converted to lowercase; *
687 ! moved the nucleation rate calcs *
688 ! from the main program unit to this subr; *
689 ! added temporary variables log_rh, log_nh3conc, log_so4vol; *
690 ! in the nuc subr, apply limits to the input parameters; *
691 ! converted to f90; *
692 ! *
693 ! reference: *
694 ! napari, i., m, noppel, h. vehkamaki, . and . m. kulmala, *
695 ! parameterization of ternary nucleation rates for *
696 ! h2so4-nh3-h2o vapors. j. geophys. res., 107, 4381, 2002b. *
697 ! *
698 ! note: *
699 ! advantage of this parameterization *
700 ! this parameterization reproduces the ternary nucleation *
701 ! rates obtained from the full model within the range of *
702 ! one order of magnitude, with a cpu saving by a factor *
703 ! of 10e5. *
704 ! *
705 ! limitations / assumptions of this parameterization *
706 ! 1. the limits of validity are *
707 ! t = 240-300 k, rh=0.05-0.95 *
708 ! [h2so4]=10e4-10e9 cm-3 (4e-4 - 40 ppt at stp) *
709 ! [nh3]=0.1-100 ppt, *
710 ! j=10e-5 - 10e6 cm-3 s-1 *
711 ! 2. it cannot be used to obtain binary h2o-h2so4 or h2o-nh3 *
712 ! limit, due to the logarithmic dependencies on rh, *
713 ! [h2so4], and [nh3] *
714 ! 3. the fit is the worst at high temperatures ( > 298 k) *
715 ! and low nucleation rates ( < 0.01 cm-3 s-1), but it is *
716 ! more accurate at significant nucleation rates (0.01-0.1 *
717 ! cm-3 s-1), thus adequate for simulating ternary *
718 ! nucleation in the atmosphere *
719 ! *
720 !*********************************************************************
722 implicit none
724 ! subr arguments (in)
725 real, intent(in) :: temp_in ! temperature, in k
726 real, intent(in) :: rh_in ! relative humidity, as fraction
727 real, intent(in) :: nh3conc_in ! mixing ratio of nh3, pptv
728 real, intent(in) :: so4vol_in ! concentration of h2so4, cm-3
730 ! subr arguments (out)
731 real, intent(out) :: ratenuclt ! ternary nucleation rate, j, cm-3 s-1
732 real, intent(out) :: rateloge ! ln (j)
734 real, intent(out) :: cnum_h2so4 ! number of h2so4 molecules
735 ! in the critical nucleus
736 real, intent(out) :: cnum_nh3 ! number of nh3 molecules
737 ! in the critical nucleus
738 real, intent(out) :: cnum_tot ! total number of molecules
739 ! in the critical nucleus
740 real, intent(out) :: radius_cluster ! the radius of cluster, nm
742 ! local variables
743 integer ncoeff
744 parameter ( ncoeff = 4 ) ! total fitting coefficients at each t
745 ! corresponds to 3rd order polynomial
746 integer npoly
747 parameter ( npoly = 20 ) ! total number of polynomial functions
749 integer n ! loop index for functions of polynomial
751 ! polynomial functions
752 real f ( npoly )
754 ! coefficients of polynomials
755 real a ( ncoeff, npoly )
757 real temp, rh, nh3conc, so4vol ! bounded values of the input args
758 real log_rh, log_nh3conc, log_so4vol
760 ! coefficients of polynomials fi(t), i = 1, 20
761 data a / -0.355297, -33.8449, 0.34536, -8.24007e-4, &
762 3.13735, -0.772861, 5.61204e-3, -9.74576e-6, &
763 19.0359, -0.170957, 4.79808e-4, -4.14699e-7, &
764 1.07605, 1.48932, -7.96052e-3, 7.61229e-6, &
765 6.0916, -1.25378, 9.39836e-3, -1.74927e-5, &
766 0.31176, 1.64009, -3.43852e-3, -1.09753e-5, &
767 -2.00738e-2, -0.752115, 5.25813e-3, -8.98038e-6, &
768 0.165536, 3.26623, -4.89703e-2, 1.46967e-4, &
769 6.52645, -0.258002, 1.43456e-3, -2.02036e-6, &
770 3.68024, -0.204098, 1.06259e-3, -1.2656e-6 , &
771 -6.6514e-2, -7.82382, 1.22938e-2, 6.18554e-5, &
772 0.65874, 0.190542, -1.65718e-3, 3.41744e-6, &
773 5.99321e-2, 5.96475, -3.62432e-2, 4.93337e-5, &
774 -0.732731, -1.84179e-2, 1.47186e-4, -2.37711e-7, &
775 0.728429, 3.64736, -2.7422e-2, 4.93478e-5, &
776 41.3016, -0.35752, 9.04383e-4, -5.73788e-7, &
777 -0.160336, 8.89881e-3, -5.39514e-5, 8.39522e-8, &
778 8.57868, -0.112358, 4.72626e-4, -6.48365e-7, &
779 5.30167e-2, -1.98815, 1.57827e-2, -2.93564e-5, &
780 -2.32736, 2.34646e-2, -7.6519e-5, 8.0459e-8 /
783 !
784 ! apply limits to input parameters
785 !
786 temp = max( 240.15, min (300.15, temp_in ) )
787 rh = max( 0.05, min (0.95, rh_in ) )
788 so4vol = max( 1.0e4, min (1.0e9, so4vol_in ) )
789 nh3conc = max( 0.1, min (100.0, nh3conc_in ) )
791 !
792 ! calculate the functions of polynomials, fi(t), i = 1, npoly
793 ! at temperature j
794 !
795 do n = 1, npoly
796 f ( n ) = a ( 1, n ) + a ( 2, n ) * temp &
797 + a ( 3, n ) * ( temp ) ** 2.0 &
798 + a ( 4, n ) * ( temp ) ** 3.0
799 end do
801 !
802 ! calculate ln (j)
803 !
804 log_rh = log ( rh )
805 log_nh3conc = log ( nh3conc )
806 log_so4vol = log ( so4vol )
807 rateloge = -84.7551 &
808 + f ( 1 ) / log_so4vol &
809 + f ( 2 ) * ( log_so4vol ) &
810 + f ( 3 ) * ( log_so4vol ) **2.0 &
811 + f ( 4 ) * ( log_nh3conc ) &
812 + f ( 5 ) * ( log_nh3conc ) **2.0 &
813 + f ( 6 ) * rh &
814 + f ( 7 ) * ( log_rh ) &
815 + f ( 8 ) * ( log_nh3conc / &
816 log_so4vol ) &
817 + f ( 9 ) * ( log_nh3conc * &
818 log_so4vol ) &
819 + f ( 10 ) * rh * &
820 ( log_so4vol ) &
821 + f ( 11 ) * ( rh / &
822 log_so4vol ) &
823 + f ( 12 ) * ( rh * &
824 log_nh3conc ) &
825 + f ( 13 ) * ( log_rh / &
826 log_so4vol ) &
827 + f ( 14 ) * ( log_rh * &
828 log_nh3conc ) &
829 + f ( 15 ) * (( log_nh3conc ) ** 2.0 &
830 / log_so4vol ) &
831 + f ( 16 ) * ( log_so4vol * &
832 ( log_nh3conc ) ** 2.0 ) &
833 + f ( 17 ) * (( log_so4vol ) ** 2.0 * &
834 log_nh3conc ) &
835 + f ( 18 ) * ( rh * &
836 ( log_nh3conc ) ** 2.0 ) &
837 + f ( 19 ) * ( rh * log_nh3conc &
838 / log_so4vol ) &
839 + f ( 20 ) * (( log_so4vol ) ** 2.0 * &
840 ( log_nh3conc ) ** 2.0 )
842 ratenuclt = exp ( rateloge )
843 !
844 ! calculate number of molecules and critical radius of the cluster
845 !
846 cnum_h2so4 = 38.1645 + 0.774106 * rateloge &
847 + 0.00298879 * ( rateloge ) ** 2.0 &
848 - 0.357605 * temp &
849 - 0.00366358 * temp * rateloge &
850 + 0.0008553 * ( temp ) ** 2.0
852 cnum_nh3 = 26.8982 + 0.682905 * rateloge &
853 + 0.00357521 * ( rateloge ) ** 2.0 &
854 - 0.265748 * temp &
855 - 0.00341895 * temp * rateloge &
856 + 0.000673454 * ( temp ) ** 2.0
858 cnum_tot = 79.3484 + 1.7384 * rateloge &
859 + 0.00711403 * ( rateloge ) ** 2.0 &
860 - 0.744993 * temp &
861 - 0.00820608 * temp * rateloge &
862 + 0.0017855 * ( temp ) ** 2.0
864 radius_cluster = 0.141027 - 0.00122625 * rateloge &
865 - 7.82211e-6 * ( rateloge ) ** 2.0 &
866 - 0.00156727 * temp &
867 - 0.00003076 * temp * rateloge &
868 + 0.0000108375 * ( temp ) ** 2.0
870 return
871 end subroutine ternary_nuc_napari
875 !-----------------------------------------------------------------------
876 subroutine wexler_nuc_mosaic_1box( &
877 dtnuc, temp_in, rh_in, cair, &
878 qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
879 nsize, maxd_asize, volumlo_sect, volumhi_sect, &
880 isize_nuc, qnuma_del, qso4a_del, qnh4a_del, &
881 qh2so4_del, qnh3_del, dens_nh4so4a )
882 !.......................................................................
883 !
884 ! calculates new particle production from h2so4-h2o binary nucleation
885 ! over timestep dtnuc, using the wexler et al. (1994) parameterization
886 !
887 ! the size of new particles is the lower-bound size of the host code's
888 ! smallest size bin. their composition is so4 and nh4, since the nuclei
889 ! would incorporate nh3 as they grow from ~1 nm to the lower-bound size.
890 ! (the new particle composition can be forced to pure so4 by setting
891 ! the qnh3_avg & qnh3_cur input arguments to 0.0).
892 !
893 ! revision history
894 ! coded by rc easter, pnnl, 20-mar-2006
895 !
896 ! key routines called: none
897 !
898 ! references:
899 ! wexler, a. s., f. w. lurmann, and j. h. seinfeld,
900 ! modelling urban and regional aerosols -- i. model development,
901 ! atmos. environ., 28, 531-546, 1994.
902 !
903 !.......................................................................
904 implicit none
906 ! subr arguments (in)
907 real, intent(in) :: dtnuc ! nucleation time step (s)
908 real, intent(in) :: temp_in ! temperature, in k
909 real, intent(in) :: rh_in ! relative humidity, as fraction
910 real, intent(in) :: cair ! dry-air molar density (mole-air/cm3)
912 real, intent(in) :: qh2so4_avg, qh2so4_cur ! gas h2so4 mixing ratios (ppm)
913 real, intent(in) :: qnh3_avg, qnh3_cur ! gas nh3 mixing ratios (ppm)
914 ! qxxx_cur = current value (at end of condensation)
915 ! qxxx_avg = average value (from start to end of condensation)
917 integer, intent(in) :: nsize ! number of aerosol size bins
918 integer, intent(in) :: maxd_asize ! dimension for volumlo_sect, ...
919 real, intent(in) :: volumlo_sect(maxd_asize) ! dry volume at lower bnd of bin (cm3)
920 real, intent(in) :: volumhi_sect(maxd_asize) ! dry volume at upper bnd of bin (cm3)
922 ! subr arguments (out)
923 integer, intent(out) :: isize_nuc ! size bin into which new particles go
924 real, intent(out) :: qnuma_del ! change to aerosol number mixing ratio (#/kg)
925 real, intent(out) :: qso4a_del ! change to aerosol so4 mixing ratio (ug/kg)
926 real, intent(out) :: qnh4a_del ! change to aerosol nh4 mixing ratio (ug/kg)
927 real, intent(out) :: qh2so4_del ! change to gas h2so4 mixing ratio (ppm)
928 real, intent(out) :: qnh3_del ! change to gas nh3 mixing ratio (ppm)
929 ! aerosol changes are > 0; gas changes are < 0
931 ! subr arguments (inout)
932 real, intent(inout) :: dens_nh4so4a ! dry-density of the new nh4-so4 aerosol mass (g/cm3)
933 ! use 'in' value only if it is between 1.6-2.0 g/cm3
935 ! local variables
936 integer i
937 integer, save :: icase = 0, icase_reldiffmax = 0
939 real, parameter :: pi = 3.1415926536
940 real, parameter :: avogad = 6.022e23 ! avogadro number (molecules/mole)
941 real, parameter :: mw_air = 28.966 ! dry-air mean molecular weight (g/mole)
943 ! dry densities (g/cm3) molecular weights of aerosol
944 ! ammsulf, ammbisulf, and sulfacid (from mosaic dens_electrolyte values)
945 real, parameter :: dens_ammsulf = 1.769
946 real, parameter :: dens_ammbisulf = 1.78
947 real, parameter :: dens_sulfacid = 1.841
949 ! molecular weights (g/mole) of aerosol ammsulf, ammbisulf, and sulfacid
950 ! for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
951 ! because we don't keep track of aerosol hion mass
952 real, parameter :: mw_ammsulf = 132.0
953 real, parameter :: mw_ammbisulf = 114.0
954 real, parameter :: mw_sulfacid = 96.0
955 ! molecular weights of aerosol sulfate and ammonium
956 real, parameter :: mw_so4a = 96.0
957 real, parameter :: mw_nh4a = 18.0
959 real, save :: reldiffmax = 0.0
961 real ch2so4_crit ! critical h2so4 conc (ug/m3)
962 real dens_part ! "grown" single-particle dry density (g/cm3)
963 real duma, dumb, dumc, dume
964 real dum_m1, dum_m2, dum_m3, dum_n1, dum_n2, dum_n3
965 real fogas, foso4a, fonh4a, fonuma
966 real mass_part ! "grown" single-particle mass (g)
967 real molenh4a_per_moleso4a ! (mole aerosol nh4)/(mole aerosol so4)
968 real qh2so4_crit ! critical h2so4 mixrat (ppm)
969 real qh2so4_avail ! amount of h2so4 available for new particles (ppm)
970 real vol_part ! "grown" single-particle volume (cm3)
973 !
974 ! initialization output arguments with "zero nucleation" values
975 !
976 isize_nuc = 1
977 qnuma_del = 0.0
978 qso4a_del = 0.0
979 qnh4a_del = 0.0
980 qh2so4_del = 0.0
981 qnh3_del = 0.0
983 !
984 ! calculate critical h2so4 concentration (ug/m3) and mixing ratio (mole/mole-air)
985 !
986 ch2so4_crit = 0.16 * exp( 0.1*temp_in - 3.5*rh_in - 27.7 )
987 ! ch2so4 = (ug-h2so4/m3-air)
988 ! ch2so4*1.0e-12/mwh2so4 = (mole-h2so4/cm3-air)
989 qh2so4_crit = (ch2so4_crit*1.0e-12/98.0)/cair
990 qh2so4_avail = qh2so4_cur - qh2so4_crit
992 ! if "available" h2so4 vapor < 4.0e-18 mole/mole-air ~= 1.0e2 molecules/cm3,
993 ! exit with new particle formation = 0
994 if (qh2so4_avail .le. 4.0e-18) then
995 return
996 end if
998 ! determine size bin into which the new particles go
999 isize_nuc = 1
1000 vol_part = volumlo_sect(1)
1002 !
1003 ! determine composition and density of the "grown particles"
1004 ! the grown particles are assumed to be liquid
1005 ! (since critical clusters contain water)
1006 ! so any (nh4/so4) molar ratio between 0 and 2 is allowed
1007 ! assume that the grown particles will have
1008 ! (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
1009 !
1010 if (qnh3_cur .ge. qh2so4_avail) then
1011 ! combination of ammonium sulfate and ammonium bisulfate
1012 ! dum_n1 & dum_n2 = mole fractions of the ammsulf & ammbisulf
1013 dum_n1 = (qnh3_cur/qh2so4_avail) - 1.0
1014 dum_n1 = max( 0.0, min( 1.0, dum_n1 ) )
1015 dum_n2 = 1.0 - dum_n1
1016 dum_n3 = 0.0
1017 else
1018 ! combination of ammonium bisulfate and sulfuric acid
1019 ! dum_n2 & dum_n3 = mole fractions of the ammbisulf & sulfacid
1020 dum_n1 = 0.0
1021 dum_n2 = (qnh3_cur/qh2so4_avail)
1022 dum_n2 = max( 0.0, min( 1.0, dum_n2 ) )
1023 dum_n3 = 1.0 - dum_n2
1024 end if
1026 dum_m1 = dum_n1*mw_ammsulf
1027 dum_m2 = dum_n2*mw_ammbisulf
1028 dum_m3 = dum_n3*mw_sulfacid
1029 dens_part = (dum_m1 + dum_m2 + dum_m3)/ &
1030 ((dum_m1/dens_ammsulf) + (dum_m2/dens_ammbisulf) &
1031 + (dum_m3/dens_sulfacid))
1032 ! 25-jul-2006 - use 'in' value only if it is between 1.6-2.0 g/cm3
1033 if (abs(dens_nh4so4a-1.8) .le. 0.2) then
1034 dens_part = dens_nh4so4a
1035 else
1036 dens_nh4so4a = dens_part
1037 end if
1038 mass_part = vol_part*dens_part
1039 molenh4a_per_moleso4a = 2.0*dum_n1 + dum_n2
1042 ! changes to h2so4 & nh3 gas (in mole/mole-air), limited by amounts available
1043 duma = 0.9999
1044 qh2so4_del = min( duma*qh2so4_cur, qh2so4_avail )
1045 qnh3_del = min( duma*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
1046 qh2so4_del = -qh2so4_del
1047 qnh3_del = -qnh3_del
1049 ! changes to so4 & nh4 aerosol (in mole/mole-air)
1050 qso4a_del = -qh2so4_del
1051 qnh4a_del = -qnh3_del
1052 ! change to aerosol number (in #/mole-air)
1053 qnuma_del = (qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
1056 return
1057 end subroutine wexler_nuc_mosaic_1box
1062 !-----------------------------------------------------------------------
1066 end module module_mosaic_newnuc