module_mosaic_newnuc.F

References to this file elsewhere.
1 !**********************************************************************************  
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of 
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************  
9 	module module_mosaic_newnuc
10 
11 
12 
13 	use module_peg_util
14 
15 
16 
17 	implicit none
18 
19 
20 
21 	contains
22 
23 
24 
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
46 
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
56 
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
66 
67 
68 !   local variables
69 	integer, parameter :: newnuc_method = 2
70 
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, save :: ncount(10)
77 	integer :: p1st
78 
79 	real, parameter :: densdefault = 2.0
80 	real, parameter :: smallmassbb = 1.0e-30
81 
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
90 
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
101 
102 	real,save :: dumveca(10), dumvecb(10), dumvecc(10), dumvecd(10), dumvece(10)
103 	real :: volumlo_nuc(maxd_asize), volumhi_nuc(maxd_asize)
104 
105 	character(len=100) :: msg
106 
107     p1st = PARAM_FIRST_SCALAR
108 
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
118 
119 
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
126 
127 
128 !   set variables that do not change
129 	idiagbb = idiagbb_in
130 
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)
136 
137 
138 !   loop over subareas (currently only 1) and vertical levels
139 	do 2900 m = 1, nsubareas
140 
141 	do 2800 k = kclm_calcbgn, kclm_calcend
142 
143 
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
154 
155 
156 	ncount(1) = ncount(1) + 1
157 	if (afracsubarea(k,m) .lt. 1.e-4) goto 2700
158 
159 	cair_box = cairclm(k)
160 	temp_box = rsub(ktemp,k,m)
161 	rh_box = relhumclm(k)
162 
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)) )
167 
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
173 
174 	dens_nh4so4a = dens_so4_aer
175 
176 	isize = 0
177 
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
197 
198 
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
212 
213 
214 !   check for zero new particles
215 	if (qnuma_del .le. 0.0) goto 2700
216 
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 )
226 
227 
228 	ncount(2) = ncount(2) + 1
229 
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 )
233 
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)
245 
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
269 
270 
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
278 
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
300 
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
316 
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
332 
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
343 
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
350 
351 
352 2700	continue
353 
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 )
376 
377 
378 2800	continue	! k levels
379 
380 2900	continue	! subareas
381 
382 
383 	return
384 	end subroutine mosaic_newnuc_1clm
385 
386 
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
422 
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)
428 
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)
433 
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)
438 
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
447 
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
451 
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)
461 
462 !       common / ternary_nuc_yzhang_cmn01 /   &
463 !         ratenuclt, rateloge,   &
464 !         cnum_h2so4, cnum_nh3, cnum_tot, radius_cluster
465 
466 
467 ! local variables
468         integer i
469         integer, save :: icase = 0, icase_reldiffmax = 0
470 
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)
474 
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
480 
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
490 
491         real, save :: reldiffmax = 0.0
492 
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)
509 
510 
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
523 
524 
525 !
526 ! make call to napari parameterization routine
527 !
528 
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
532 
533 ! convert h2so4 from mole/mole-air to molecules/cm3
534 !       qh2so4_cur = (so4vol(k)/avogad) / cair
535         so4vol_in  = (qh2so4_avg) * cair * avogad
536 
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 )
541 
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
545 
546 
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
568 
569 
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
593 
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)
609 
610 
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
620 
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
626 
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 )
637 
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
641 
642 
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?
652 
653 
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
660 
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
666 
667 
668         return
669         end subroutine ternary_nuc_mosaic_1box
670 
671 
672 
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 !*********************************************************************
721 
722       implicit none
723 
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
729 
730 ! subr arguments (out)
731         real, intent(out) :: ratenuclt        ! ternary nucleation rate, j, cm-3 s-1
732         real, intent(out) :: rateloge         ! ln (j)
733 
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
741 
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
748 
749         integer n                     ! loop index for functions of polynomial
750 
751 ! polynomial functions
752         real f ( npoly )
753 
754 ! coefficients of polynomials
755         real a  ( ncoeff, npoly )
756 
757         real temp, rh, nh3conc, so4vol     ! bounded values of the input args
758         real log_rh, log_nh3conc, log_so4vol
759 
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 /
781 
782 
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 ) )
790 
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
800 
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 )
841 
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
851 
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
857 
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
863 
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
869 
870         return
871         end subroutine ternary_nuc_napari
872 
873 
874 
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
905 
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)
911 
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)
916 
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)
921 
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
930 
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
934 
935 ! local variables
936         integer i
937         integer, save :: icase = 0, icase_reldiffmax = 0
938 
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)
942 
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
948 
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
958 
959         real, save :: reldiffmax = 0.0
960 
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)
971 
972 
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
982 
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
991 
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
997 
998 ! determine size bin into which the new particles go
999         isize_nuc = 1
1000         vol_part = volumlo_sect(1)
1001 
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
1025 
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
1040 
1041 
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
1048 
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
1054 
1055 
1056         return
1057         end subroutine wexler_nuc_mosaic_1box
1058 
1059 
1060 
1061 
1062 !-----------------------------------------------------------------------
1063 
1064 
1065 
1066 	end module module_mosaic_newnuc