module_mosaic_initmixrats.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_initmixrats
10 
11 ! rce 2005-feb-18 - several fixes for indices of dlo/hi_sect, volumlo/hi_sect,
12 !     which are now (isize,itype)
13 
14 ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
15 !     variables in module_data_mosaic_asect
16 
17 	USE module_state_description
18 
19 	integer, parameter :: mosaic_init_wrf_mixrats_flagaa = 1
20                             ! turns subr mosaic_init_wrf_mixrats on/off
21 
22 	contains
23 
24 
25 !-----------------------------------------------------------------------
26 !wig 10-May-2004, added phb, ph, g
27 	subroutine mosaic_init_wrf_mixrats(  &
28 		iflagaa, config_flags,       &
29 		chem, alt, z_at_w, g,        &
30 		ids,ide, jds,jde, kds,kde,   &
31 		ims,ime, jms,jme, kms,kme,   &
32 		its,ite, jts,jte, kts,kte    )
33 
34 !
35 !   initializes the species and number mixing ratios for each section
36 !
37 !   this top level routine simply calls other routines depending on value
38 !	of config_flags%aer_ic_opt
39 !
40 	use module_configure, only:  grid_config_rec_type
41 	use module_state_description, only:  num_chem, param_first_scalar,   &
42 			aer_ic_pnnl
43 	use module_data_mosaic_asect
44 	use module_data_mosaic_other
45 	use module_peg_util, only:  peg_message, peg_error_fatal
46 
47 	implicit none
48 
49 
50 !   subr arguments
51 	type(grid_config_rec_type), intent(in) :: config_flags
52 
53 	integer, intent(in) ::   &
54 		iflagaa,   &
55 		ids, ide, jds, jde, kds, kde,   &
56 		ims, ime, jms, jme, kms, kme,   &
57 		its, ite, jts, jte, kts, kte
58 
59 	real, intent(inout),   &
60 		dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
61 		chem
62 
63 !wig 10-May-2004, added phb, ph, g
64 	real, intent(in),      &
65 		dimension( ims:ime, kms:kme, jms:jme ) :: &
66         alt, z_at_w
67     real, intent(in) :: g
68 
69 !   local variables
70         integer :: i, ic, j, k
71 
72 	if (config_flags%aer_ic_opt == aer_ic_pnnl) then
73 !wig 10-May-2004, added phb, ph, g
74 	    call mosaic_init_wrf_mixrats_opt2(   &
75 		iflagaa, config_flags,           &
76 		chem, z_at_w, g,                 &
77 		ids,ide, jds,jde, kds,kde,       &
78 		ims,ime, jms,jme, kms,kme,       &
79 		its,ite, jts,jte, kts,kte        )
80 
81 	else
82 	    call mosaic_init_wrf_mixrats_opt1(   &
83 		iflagaa, config_flags,   &
84 		chem,   &
85 		ids,ide, jds,jde, kds,kde,   &
86 		ims,ime, jms,jme, kms,kme,   &
87 		its,ite, jts,jte, kts,kte    )
88 
89 	end if
90 
91 ! Aerosol species are returned from above in concentration units. Convert
92 ! them to mixing ratio for use in advection, etc.
93     do ic = p_so4_a01,num_chem
94        do j = jts,jte
95           do k = kts,kte-1
96              do i = its,ite
97                 chem(i,k,j,ic) = chem(i,k,j,ic)*alt(i,k,j)
98              end do
99           end do
100        end do
101     end do
102 
103 ! Fill the top z-staggered location to prevent trouble during advection.
104     do ic = p_so4_a01,num_chem
105        do j = jts,jte
106           do i = its,ite
107              chem(i,kte,j,ic) = chem(i,kte-1,j,ic)
108           end do
109        end do
110     end do
111 
112 	return
113 	end subroutine mosaic_init_wrf_mixrats
114 
115 
116 !-----------------------------------------------------------------------
117 	subroutine mosaic_init_wrf_mixrats_opt1(   &
118 		iflagaa, config_flags,   &
119 		chem,   &
120 		ids,ide, jds,jde, kds,kde,   &
121 		ims,ime, jms,jme, kms,kme,   &
122 		its,ite, jts,jte, kts,kte    )
123 
124 !
125 !   initializes the species and number mixing ratios for each section
126 !   based on user-specified lognormal modes that span the size distribution
127 !
128 ! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
129 !     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now; 
130 !     added l1dum
131 !
132 	use module_configure, only:  grid_config_rec_type
133 	use module_state_description, only:  num_chem, param_first_scalar
134 	use module_data_mosaic_asect
135 	use module_data_mosaic_other
136 	use module_peg_util, only:  peg_message, peg_error_fatal
137 
138 	implicit none
139 
140 !   subr arguments
141 	type(grid_config_rec_type), intent(in) :: config_flags
142 
143 	integer, intent(in) ::   &
144 		iflagaa,   &
145 		ids, ide, jds, jde, kds, kde,   &
146 		ims, ime, jms, jme, kms, kme,   &
147 		its, ite, jts, jte, kts, kte
148 
149 	real, intent(inout),   &
150 		dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
151 		chem
152 
153 !   local variables
154 	integer i, j, k, l, ll, l1, l3, l1dum, m, mdum, nsm
155 	integer it, jt, kt
156 
157 	real dum, dumdp, dumrsfc, dumvol,   &
158       	  xlo, xhi,   &
159 	  dumvol1p,  &
160       	  pdummb, zdumkm, zscalekm, zfactor
161 
162 	real vtot_nsm_ofmode(maxd_asize)
163 	real dumarr(maxd_acomp+5)
164 
165 	real erfc
166 
167 !	double precision fracnum, fracvol, tlo, thi
168 	real fracvol, tlo, thi
169 
170 	integer nmaxd_nsm
171 	parameter (nmaxd_nsm = 4)
172 
173 	integer iphase, itype, ntot_nsm
174 	integer iiprof_nsm(nmaxd_nsm)
175 	integer lldum_so4, lldum_nh4, lldum_oc, lldum_bc,   &
176 		lldum_oin, lldum_na, lldum_cl, lldum_hysw
177 
178 	real sx_nsm(nmaxd_nsm), sxr2_nsm(nmaxd_nsm),   &
179       	  x0_nsm(nmaxd_nsm), x3_nsm(nmaxd_nsm),   &
180       	  rtot_nsm(maxd_acomp,nmaxd_nsm),   &
181       	  vtot_nsm(nmaxd_nsm), xntot_nsm(nmaxd_nsm)
182 
183 	real dgnum_nsm(nmaxd_nsm), sigmag_nsm(nmaxd_nsm)
184 	real aaprof_nsm(maxd_acomp+1,nmaxd_nsm)
185 	real bbprof_nsm(nmaxd_nsm)
186 
187 	character*80 msg
188 	character*10 dumname
189 
190 
191 !   check for on/off
192 	if (mosaic_init_wrf_mixrats_flagaa .le. 0) return
193 
194 
195 !   *** currently only works for ntype_aer = 1
196 	itype = 1
197 	iphase = ai_phase
198 	m = 1
199 
200 !   set values for initialization modes
201 	iiprof_nsm(:) = 1
202 	aaprof_nsm(:,:) = 0.0
203 	bbprof_nsm(:) = 0.0
204 
205 	ntot_nsm = 4
206 	ntot_nsm = min( ntot_nsm, nsize_aer(itype) )
207 
208 	lldum_so4 = 0
209 	lldum_nh4 = 0
210 	lldum_oc  = 0
211 	lldum_bc  = 0
212 	lldum_oin = 0
213 	lldum_na  = 0
214 	lldum_cl  = 0
215 	lldum_hysw = 0
216 
217 	do ll = 1, ncomp_plustracer_aer(itype)
218 	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
219 		lptr_so4_aer(m,itype,iphase)) lldum_so4 = ll
220 	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
221 		lptr_nh4_aer(m,itype,iphase)) lldum_nh4 = ll
222 	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
223 		lptr_oc_aer(m,itype,iphase))  lldum_oc  = ll
224 	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
225 		lptr_bc_aer(m,itype,iphase))  lldum_bc  = ll
226 	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
227 		lptr_oin_aer(m,itype,iphase)) lldum_oin = ll
228 	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
229 		lptr_na_aer(m,itype,iphase))  lldum_na  = ll
230 	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
231 		lptr_cl_aer(m,itype,iphase))  lldum_cl  = ll
232 	end do
233 	if (hyswptr_aer(m,itype) .gt. 0)   &
234 		lldum_hysw = ncomp_plustracer_aer(itype) + 1
235 
236 	msg = ' '
237 	if (lldum_so4 .le. 0)   &
238 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_so4 = 0'
239 	if (lldum_nh4 .le. 0)   &
240 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_nh4 = 0'
241 	if (lldum_oc .le. 0)   &
242 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_oc = 0'
243 	if (lldum_bc .le. 0)   &
244 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_bc = 0'
245 	if (lldum_oin .le. 0)   &
246 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_oin = 0'
247 	if (lldum_na .le. 0)   &
248 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_na = 0'
249 	if (lldum_cl .le. 0)   &
250 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_cl = 0'
251 	if (lldum_hysw .le. 0)   &
252 		msg = '*** subr mosaic_init_wrf_mixrats - lldum_hysw = 0'
253 	if (msg .ne. ' ') call peg_error_fatal( lunerr, msg )
254 
255 
256 	do nsm = 1, ntot_nsm
257 
258 	if (nsm .eq. 1) then
259 !   accum mode with so4, nh4, oc, bc
260 	    dgnum_nsm( nsm) = 0.15e-4
261 	    sigmag_nsm(nsm) = 2.0
262 	    aaprof_nsm(lldum_so4,nsm) = 0.50
263 	    aaprof_nsm(lldum_nh4,nsm) = aaprof_nsm(lldum_so4,nsm)   &
264 					* (mw_nh4_aer/mw_so4_aer)
265 	    aaprof_nsm(lldum_oc,nsm) = 0.25
266 	    aaprof_nsm(lldum_bc,nsm) = 0.05
267 	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_so4,nsm) * 0.2
268 
269 	else if (nsm .eq. 2) then
270 !   aitken mode with so4, nh4, oc, bc
271 	    dgnum_nsm( nsm) = 0.03e-4
272 	    sigmag_nsm(nsm) = 2.0
273 	    aaprof_nsm(lldum_so4,nsm) = 0.50 * 0.020
274 	    aaprof_nsm(lldum_nh4,nsm) = aaprof_nsm(lldum_so4,nsm)   &
275 					* (mw_nh4_aer/mw_so4_aer)
276 	    aaprof_nsm(lldum_oc,nsm) = 0.25 * 0.020
277 	    aaprof_nsm(lldum_bc,nsm) = 0.05 * 0.020
278 	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_so4,nsm) * 0.2
279 
280 	else if (nsm .eq. 3) then
281 !   coarse dust mode with oin
282 	    dgnum_nsm( nsm) = 1.0e-4
283 	    sigmag_nsm(nsm) = 2.0
284 	    aaprof_nsm(lldum_oin,nsm) = 0.5
285 	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm( 9,nsm) * 1.0e-3
286 
287 	else if (nsm .eq. 4) then
288 !   coarse seasalt mode with na, cl
289 	    dgnum_nsm( nsm) = 2.0e-4
290 	    sigmag_nsm(nsm) = 2.0
291 	    aaprof_nsm(lldum_cl,nsm) = 0.1
292 	    aaprof_nsm(lldum_na,nsm) = aaprof_nsm(lldum_cl,nsm)   &
293 					* (mw_na_aer/mw_cl_aer)
294 	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_cl,nsm) * 0.2
295 
296 	end if
297 
298 	end do
299 
300 !   when iflagaa = 1/2/3/4, use only the nsm-mode = iflagaa
301 	if (iflagaa .gt. 0) then
302 	    do nsm = 1, ntot_nsm
303 		if (nsm .ne. iflagaa) aaprof_nsm(:,nsm) = 0.0
304 	    end do
305 	end if
306 
307 
308 
309 !
310 !   do the initialization now
311 !
312 
313 !   calculate mode parameters
314 	do nsm = 1, ntot_nsm
315 	    sx_nsm(nsm) = alog( sigmag_nsm(nsm) )
316 	    sxr2_nsm(nsm) = sx_nsm(nsm) * sqrt(2.0)
317 	    x0_nsm(nsm) = alog( dgnum_nsm(nsm) )
318 	    x3_nsm(nsm) = x0_nsm(nsm) + 3.0*sx_nsm(nsm)*sx_nsm(nsm)
319 	end do
320 
321 !   initialize rclm array to zero
322 	rclm(:,:) = 0.
323 !	rclmsvaa(:,:) = 0.
324 
325 !
326 !   loop over all vertical levels
327 !
328 !	do 12900 k = 1, ktot
329 	do 12900 k = 1, 1
330 
331 !	pdummb = 1013.*scord(k)
332 !	zdumkm = ptoz( pdummb ) * 1.e-3
333 	zdumkm = 0.0
334 
335 
336 !   for each species and nsm mode, define total mixing ratio
337 !	(for all sizes) at level k
338 !
339 !   iiprof_nsm =  +1 gives constant mixing ratio
340 !	aaprof_nsm(l,nsm) = constant mixing ratio (ppb)
341 !   iiprof_nsm =  +2 gives exponential profile
342 !	aaprof_nsm(l,nsm) = surface mixing ratio (ppb)
343 !	bbprof_nsm(l) = scale height (km)
344 !   iiprof_nsm =  +3 gives linear profile (then zero at z > zscalekm)
345 !	aaprof_nsm(l,nsm) = surface mixing ratio (ppb)
346 !	bbprof_nsm(l) = height (km) at which mixing ratio = 0
347 
348 	do nsm = 1, ntot_nsm
349 
350 	if (iiprof_nsm(nsm) .eq. 2) then
351 	    zscalekm = bbprof_nsm(nsm)
352 	    zfactor = exp( -zdumkm/zscalekm )
353 	else if (iiprof_nsm(nsm) .eq. 3) then
354 	    zscalekm = bbprof_nsm(nsm)
355 	    zfactor = max( 0., (1. - zdumkm/zscalekm) )
356 	else
357 	    zfactor = 1.0
358 	end if
359 
360 	    do ll = 1, ncomp_plustracer_aer(itype) + 1
361 		rtot_nsm(ll,nsm) = aaprof_nsm(ll,nsm) * zfactor
362 	    end do
363 
364 	end do
365 
366 !   compute total volume and number mixing ratios for each nsm mode
367 !   rtot_nsm is ug/m3,  1.0e-6*rtot is g/m3,  vtot_nsm is cm3/m3
368 	do nsm = 1, ntot_nsm
369 	    dumvol = 0.
370 	    do ll = 1, ncomp_aer(itype)
371 		dum = 1.0e-6*rtot_nsm(ll,nsm)/dens_aer(ll,itype)
372 		dumvol = dumvol + max( 0., dum )
373 	    end do
374 	    vtot_nsm(nsm) = dumvol
375 	end do
376 
377 !   now compute species and number mixing ratios for each bin
378 	do 12700 m = 1, nsize_aer(itype)
379 
380 	vtot_nsm_ofmode(m) = 0.0
381 
382 	do 12500 nsm = 1, ntot_nsm
383 
384 !   for nsm_mode = n, compute fraction of number and volume
385 !   that is in bin m
386 	xlo = alog( dlo_sect(m,itype) )
387 	xhi = alog( dhi_sect(m,itype) )
388 
389 	tlo = (xlo - x3_nsm(nsm))/sxr2_nsm(nsm)
390 	thi = (xhi - x3_nsm(nsm))/sxr2_nsm(nsm)
391 	if (tlo .ge. 0.) then
392 !	    fracvol = 0.5*( erfc(tlo) - erfc(thi) )
393 	    fracvol = 0.5*( erfc_num_recipes(tlo) - erfc_num_recipes(thi) )
394 	else
395 !	    fracvol = 0.5*( erfc(-thi) - erfc(-tlo) )
396 	    fracvol = 0.5*( erfc_num_recipes(-thi) - erfc_num_recipes(-tlo) )
397 	end if
398 	fracvol = max( fracvol, 0.0 )
399 
400 	vtot_nsm_ofmode(m) = vtot_nsm_ofmode(m)  + vtot_nsm(nsm)*fracvol
401 
402 !   now load that fraction of species-mixing-ratio
403 !   into the appropriate rclm location
404 	do ll = 1, ncomp_plustracer_aer(itype)
405 	    rclm( k, massptr_aer(ll,m,itype,iphase) ) =   &
406       	    rclm( k, massptr_aer(ll,m,itype,iphase) ) +   &
407       		fracvol*rtot_nsm(ll,nsm)
408 	end do
409 
410 	if ((iphase .eq. ai_phase) .and.   &
411 	    (lldum_hysw .gt. 0) .and.   &
412 	    (hyswptr_aer(m,itype) .gt. 0)) then
413 
414 	    rclm( k, hyswptr_aer(m,itype) ) =   &
415       	    rclm( k, hyswptr_aer(m,itype) ) +   &
416       		fracvol*rtot_nsm(lldum_hysw,nsm)
417 	end if
418 
419 12500	continue
420 
421 !   now compute number from volume
422 	dum = sqrt( dlo_sect(m,itype)*dhi_sect(m,itype) )
423 	dumvol1p = (pi/6.0)*(dum**3)
424 	rclm( k, numptr_aer(m,itype,iphase) ) = vtot_nsm_ofmode(m)/dumvol1p
425 
426 !   set water = hyswatr
427 	if ((iphase .eq. ai_phase) .and.   &
428 	    (lldum_hysw .gt. 0) .and.   &
429 	    (hyswptr_aer(m,itype) .gt. 0) .and.   &
430 	    (waterptr_aer(m,itype) .gt. 0)) then
431 
432       	    rclm( k, waterptr_aer(m,itype) ) =    &
433 	    rclm( k, hyswptr_aer(m,itype) )
434 	end if
435 
436 12700	continue
437 
438 12900	continue
439 
440 
441 !
442 !   do diagnostic output
443 !
444 
445 ! temporary
446 ! temporary
447 
448 	dumarr(:) = 0.0
449 	msg = ' '
450 	call peg_message( lunout, msg )
451 	msg = '*** subr mosaic_init_wrf_mixrats_opt1 results'
452 	call peg_message( lunout, msg )
453 	msg = '    mass in ug/m3     number in #/m3     volume in cm3/m3'
454 	call peg_message( lunout, msg )
455 	msg = ' '
456 	call peg_message( lunout, msg )
457 	msg = ' mode  l  l1  species      conc'
458 	call peg_message( lunout, msg )
459 
460         do 14390 mdum = 1, nsize_aer(itype)+1
461 	m = min( mdum, nsize_aer(itype) )
462 	msg = ' '
463 	call peg_message( lunout, msg )
464         do 14350 l = 1, ncomp_plustracer_aer(itype)+4
465 
466             if (l .le. ncomp_plustracer_aer(itype)) then
467                 l1 = massptr_aer(l,m,itype,iphase)
468 		dumname = name_aer(l,itype)
469 		dum = rclm(1,l1)
470             else if (l .eq. ncomp_plustracer_aer(itype)+1) then
471                 l1 = hyswptr_aer(m,itype)
472 		dumname = 'hystwatr'
473 		dum = rclm(1,l1)
474             else if (l .eq. ncomp_plustracer_aer(itype)+2) then
475                 l1 = waterptr_aer(m,itype)
476 		dumname = 'water'
477 		dum = rclm(1,l1)
478             else if (l .eq. ncomp_plustracer_aer(itype)+3) then
479                 l1 = numptr_aer(m,itype,iphase)
480 		dumname = 'number'
481 		dum = rclm(1,l1)
482             else if (l .eq. ncomp_plustracer_aer(itype)+4) then
483                 l1 = 0
484 		dumname = 'volume'
485 		dum = vtot_nsm_ofmode(m)
486 	    else
487 		dumname = '????????'
488 		l1 = -1
489 		dum = -1.0
490             end if
491 
492 	    l1dum = l1
493 	    if (aboxtest_testmode .gt. 0) l1dum = max( l1-1, 0 )
494 
495 	    if (mdum .le. nsize_aer(itype)) then
496 		dumarr(l) = dumarr(l) + dum
497 		write(msg,9620) m, l, l1dum, dumname, dum
498 	    else
499 		write(msg,9625) l, dumname, dumarr(l)
500 	    end if
501 	    call peg_message( lunout, msg )
502 
503 14350   continue
504 14390   continue
505 
506 9620    format( 3i4, 2x, a, 3(1pe12.3) )
507 9625    format( ' sum', i4, '   -', 2x, a, 3(1pe12.3) )
508 
509 
510 !
511 !   load the chem array
512 !
513         do 16390 m = 1, nsize_aer(itype)
514         do 16350 l = 1, 15
515 
516 	    if (l .eq. 1) then
517 		l1 = lptr_so4_aer(m,itype,iphase)
518 	    else if (l .eq. 2) then
519 		l1 = lptr_no3_aer(m,itype,iphase)
520 	    else if (l .eq. 3) then
521 		l1 = lptr_cl_aer(m,itype,iphase)
522 	    else if (l .eq. 4) then
523 		l1 = lptr_msa_aer(m,itype,iphase)
524 	    else if (l .eq. 5) then
525 		l1 = lptr_co3_aer(m,itype,iphase)
526 	    else if (l .eq. 6) then
527 		l1 = lptr_nh4_aer(m,itype,iphase)
528 	    else if (l .eq. 7) then
529 		l1 = lptr_na_aer(m,itype,iphase)
530 	    else if (l .eq. 8) then
531 		l1 = lptr_ca_aer(m,itype,iphase)
532 	    else if (l .eq. 9) then
533 		l1 = lptr_oin_aer(m,itype,iphase)
534 	    else if (l .eq. 10) then
535 		l1 = lptr_oc_aer(m,itype,iphase)
536 	    else if (l .eq. 11) then
537 		l1 = lptr_bc_aer(m,itype,iphase)
538 	    else if (l .eq. 12) then
539 		l1 = hyswptr_aer(m,itype)
540 	    else if (l .eq. 13) then
541 		l1 = waterptr_aer(m,itype)
542 	    else if (l .eq. 14) then
543 		l1 = numptr_aer(m,itype,iphase)
544 	    else
545 		goto 16350
546 	    end if
547 	    l3 = l1
548 
549 	    if ((l1 .gt. 0) .and. (l1 .le. lmaxd) .and.   &
550 	        (l3 .ge. param_first_scalar)) then
551 		do it = its, ite
552 !   *** note:  not sure what the kt limits should be here 
553 		do kt = kts, kte-1
554 		do jt = jts, jte
555 		    chem(it,kt,jt,l3) = rclm(1,l1)
556 		end do
557 		end do
558 		end do
559 	    end if
560 
561 16350   continue
562 16390   continue
563 
564 
565 !   all done
566 	return
567 	end subroutine mosaic_init_wrf_mixrats_opt1
568 
569 
570 !-----------------------------------------------------------------------
571 !wig 10-May-2004, added phb, ph, g
572 	subroutine mosaic_init_wrf_mixrats_opt2(   &
573 		iflagaa, config_flags,             &
574 		chem, z_at_w, g,                   &
575 		ids,ide, jds,jde, kds,kde,         &
576 		ims,ime, jms,jme, kms,kme,         &
577 		its,ite, jts,jte, kts,kte          )
578 
579 !
580 !   provides initial values for mosaic aerosol species (mass and number
581 !	mixing ratio) for "Texas August 2000" simulations
582 !   currently all the initial values are uniform in x, y, AND z
583 !
584 ! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
585 !     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now
586 !
587 	use module_configure, only:  grid_config_rec_type
588 	use module_state_description, only:  num_chem, param_first_scalar
589 	use module_data_mosaic_asect
590 	use module_data_mosaic_other
591 	use module_peg_util, only:  peg_message, peg_error_fatal
592 
593 	implicit none
594 
595 !   subr arguments
596 	type(grid_config_rec_type), intent(in) :: config_flags
597 
598 	integer, intent(in) ::   &
599 		iflagaa,   &
600 		ids, ide, jds, jde, kds, kde,   &
601 		ims, ime, jms, jme, kms, kme,   &
602 		its, ite, jts, jte, kts, kte
603 
604 	real, intent(inout),     &
605 		dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
606 		chem
607 
608 !wig 10-May-2004, added phb, ph, g
609         real, intent(in),        &
610 		dimension( ims:ime, kms:kme, jms:jme ) :: &
611                 z_at_w
612         real :: g
613 
614 !   local variables
615 	integer l, l1, l3, m, mdum
616 	integer iphase, itype
617 	integer it, jt, kt
618 
619 !wig 10-May-2004, added mult
620 	real dum, dumvol1p, mult
621 	real qcoar, qfine, qval
622 
623 	real :: vtot_ofmode(maxd_asize)
624 	real :: dumarr(maxd_acomp+5)
625 	real :: fr_coar(8), fr_fine(8)
626 
627 !wig 01-Apr-2005, Updated fractional breakdown between bins. (See also
628 !                 bdy_chem_value_mosaic in this file and mosaic_addemiss in 
629 !                 module_mosaic_addemiss.F) Note that the values here no
630 !                 longer match those in mosaic_addemiss.
631 !rce 10-May-2005, changed fr8b_coar to values determined by jdf
632 	real, save :: fr8b_coar(8) =   &
633         (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.300, 0.700 /) ! 10-May-2005
634 !       (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.933, 0.067 /) ! 01-Apr-2005
635 !       (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.16,  0.84  /) ! "old"
636 
637 	real, save :: fr8b_fine(8) =   &
638         (/ 0.006, 0.024, 0.231, 0.341, 0.140, 0.258, 0., 0./) !5-Apr-2005 values
639 !       (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
640 !	(/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
641 !	(/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.00, 0.00 /)
642 
643 	real :: qfine_so4, qfine_no3, qfine_cl, qfine_msa,   &
644 		qfine_co3, qfine_nh4, qfine_na, qfine_ca, qfine_oin,   &
645 		qfine_oc, qfine_bc, qfine_hysw, qfine_watr, qfine_vol
646 	real :: qcoar_so4, qcoar_no3, qcoar_cl, qcoar_msa,   &
647 		qcoar_co3, qcoar_nh4, qcoar_na, qcoar_ca, qcoar_oin,   &
648 		qcoar_oc, qcoar_bc, qcoar_hysw, qcoar_watr, qcoar_vol
649 
650 !wig 10-May-2004, added z
651         real, dimension( ims:ime, kms:kme, jms:jme ) :: z
652 
653 	character*80 msg
654 	character*10 dumname
655 
656 
657 !   *** currently only works for ntype_aer = 1
658 	itype = 1
659 	iphase = ai_phase
660 	m = 1
661 
662 !wig 10-May-2004, block begin...
663 ! calculate the mid-level heights
664         do it = its, ite
665            do kt = kts, kte-1
666               do jt = jts, jte
667                  z(it,kt,jt) = (z_at_w(it,kt,jt)+z_at_w(it,kt+1,jt))*0.5
668               end do
669            end do
670         end do
671 !wig 10-May-2004, ...end block
672 
673 !   set the partitioning fractions for either 8 or 4 bins
674 	if (nsize_aer(itype) == 8) then
675 	    fr_coar(:) = fr8b_coar(:)
676 	    fr_fine(:) = fr8b_fine(:)
677 	else if (nsize_aer(itype) == 4) then
678 	    do m = 1, nsize_aer(itype)
679 		fr_coar(m) = fr8b_coar(2*m) + fr8b_coar(2*m-1)
680 		fr_fine(m) = fr8b_fine(2*m) + fr8b_fine(2*m-1)
681 	    end do
682 	else
683 	    write(msg,'(a,i5)')   &
684 		'subr mosaic_init_wrf_mixrats_opt2' //   &
685 		' - nsize_aer(itype) must be 4 or 8 but = ', nsize_aer(itype)
686 	    call peg_error_fatal( lunout, msg )
687 	end if
688 
689 !
690 !   compute initial values (currently uniform in x, y, AND z)
691 !   load them into the chem array
692 !   also load them into the rclm array for diagnostic output
693 !
694 	rclm(:,:) = 0.0
695 	vtot_ofmode(:) = 0.0
696 
697 !   Set values for fine and coarse mass concentrations, and
698 !   compute fine/coarse volume concentrations. These are also set in
699 !   bdy_chem_value_mosaic, below.
700 !   (these are latest values from 1-Apr-2005 discussions)
701 !
702 !   rce 10-may-2005 - changed qfine_cl, _na, _oin to values determined by jdf
703 	qfine_so4 = 2.14
704 	qcoar_so4 = 0.242
705 	qfine_no3 = 0.11
706 	qcoar_no3 = 0.03
707 !	qfine_cl  = 0.3
708 	qfine_cl  = 0.14     ! 10-May-2005
709 	qcoar_cl  = 0.139
710 	qfine_msa = 0.0
711 	qcoar_msa = 0.0
712 	qfine_co3 = 0.0
713 	qcoar_co3 = 0.0
714 	qfine_nh4 = 0.83
715 	qcoar_nh4 = 0.10
716 !	qfine_na  = 0.2
717 	qfine_na  = 0.1      ! 10-May-2005
718 	qcoar_na  = 0.09
719 	qfine_ca  = 0.0
720 	qcoar_ca  = 0.0
721 !	qfine_oin = 6.2
722 	qfine_oin = 3.48     ! 10-May-2005
723 	qcoar_oin = 0.35
724 	qfine_oc  = 1.00
725 	qcoar_oc  = 1.50
726 	qfine_bc  = 0.2
727 	qcoar_bc  = 0.075
728 	qfine_hysw = 0.0
729 	qcoar_hysw = 0.0
730 	qfine_watr = 0.0
731 	qcoar_watr = 0.0
732 
733 !!$! old values from 23-Apr-2004:
734 !!$	qfine_so4 = 2.554
735 !!$	qcoar_so4 = 0.242
736 !!$	qfine_no3 = 0.07
737 !!$	qcoar_no3 = 0.03
738 !!$	qfine_cl  = 0.324
739 !!$	qcoar_cl  = 0.139
740 !!$	qfine_msa = 0.0
741 !!$	qcoar_msa = 0.0
742 !!$	qfine_co3 = 0.0
743 !!$	qcoar_co3 = 0.0
744 !!$	qfine_nh4 = 0.98
745 !!$	qcoar_nh4 = 0.10
746 !!$	qfine_na  = 0.21
747 !!$	qcoar_na  = 0.09
748 !!$	qfine_ca  = 0.0
749 !!$	qcoar_ca  = 0.0
750 !!$	qfine_oin = 0.15
751 !!$	qcoar_oin = 0.35
752 !!$	qfine_oc  = 1.00
753 !!$	qcoar_oc  = 1.50
754 !!$	qfine_bc  = 0.175
755 !!$	qcoar_bc  = 0.075
756 !!$	qfine_hysw = 0.0
757 !!$	qcoar_hysw = 0.0
758 !!$	qfine_watr = 0.0
759 !!$	qcoar_watr = 0.0
760 
761 
762 ! qfine_so4 ... are ug/m3 so 1.0e-6 factor gives g/m3
763 ! dens_so4  ... are g/cm3;  qfine_vol ... are cm3/m3
764 	qfine_vol = 1.0e-6 * (   &
765 	    (qfine_so4/dens_so4_aer) + (qfine_no3/dens_no3_aer) +   &
766 	    (qfine_cl /dens_cl_aer ) + (qfine_msa/dens_msa_aer) +   &
767 	    (qfine_co3/dens_co3_aer) + (qfine_nh4/dens_nh4_aer) +   &
768 	    (qfine_na /dens_na_aer ) + (qfine_ca /dens_ca_aer ) +   &
769 	    (qfine_oin/dens_oin_aer) + (qfine_oc /dens_oc_aer ) +   &
770 	    (qfine_bc /dens_bc_aer ) )
771 	qcoar_vol =  1.0e-6 * (  &
772 	    (qcoar_so4/dens_so4_aer) + (qcoar_no3/dens_no3_aer) +   &
773 	    (qcoar_cl /dens_cl_aer ) + (qcoar_msa/dens_msa_aer) +   &
774 	    (qcoar_co3/dens_co3_aer) + (qcoar_nh4/dens_nh4_aer) +   &
775 	    (qcoar_na /dens_na_aer ) + (qcoar_ca /dens_ca_aer ) +   &
776 	    (qcoar_oin/dens_oin_aer) + (qcoar_oc /dens_oc_aer ) +   &
777 	    (qcoar_bc /dens_bc_aer ) )
778 
779         do 2900 m = 1, nsize_aer(itype)
780         do 2800 l = 1, 15
781 
782 	    if (l .eq. 1) then
783 		l1 = lptr_so4_aer(m,itype,iphase)
784 		qfine = qfine_so4
785 		qcoar = qcoar_so4
786 	    else if (l .eq. 2) then
787 		l1 = lptr_no3_aer(m,itype,iphase)
788 		qfine = qfine_no3
789 		qcoar = qcoar_no3
790 	    else if (l .eq. 3) then
791 		l1 = lptr_cl_aer(m,itype,iphase)
792 		qfine = qfine_cl
793 		qcoar = qcoar_cl
794 	    else if (l .eq. 4) then
795 		l1 = lptr_msa_aer(m,itype,iphase)
796 		qfine = qfine_msa
797 		qcoar = qcoar_msa
798 	    else if (l .eq. 5) then
799 		l1 = lptr_co3_aer(m,itype,iphase)
800 		qfine = qfine_co3
801 		qcoar = qcoar_co3
802 	    else if (l .eq. 6) then
803 		l1 = lptr_nh4_aer(m,itype,iphase)
804 		qfine = qfine_nh4
805 		qcoar = qcoar_nh4
806 	    else if (l .eq. 7) then
807 		l1 = lptr_na_aer(m,itype,iphase)
808 		qfine = qfine_na
809 		qcoar = qcoar_na
810 	    else if (l .eq. 8) then
811 		l1 = lptr_ca_aer(m,itype,iphase)
812 		qfine = qfine_ca
813 		qcoar = qcoar_ca
814 	    else if (l .eq. 9) then
815 		l1 = lptr_oin_aer(m,itype,iphase)
816 		qfine = qfine_oin
817 		qcoar = qcoar_oin
818 	    else if (l .eq. 10) then
819 		l1 = lptr_oc_aer(m,itype,iphase)
820 		qfine = qfine_oc
821 		qcoar = qcoar_oc
822 	    else if (l .eq. 11) then
823 		l1 = lptr_bc_aer(m,itype,iphase)
824 		qfine = qfine_bc
825 		qcoar = qcoar_bc
826 	    else if (l .eq. 12) then
827 		l1 = hyswptr_aer(m,itype)
828 		qfine = qfine_hysw
829 		qcoar = qcoar_hysw
830 	    else if (l .eq. 13) then
831 		l1 = waterptr_aer(m,itype)
832 		qfine = qfine_watr
833 		qcoar = qcoar_watr
834 	    else if (l .eq. 14) then
835 		l1 = numptr_aer(m,itype,iphase)
836 		dumvol1p = sqrt(volumlo_sect(m,itype)*volumhi_sect(m,itype))
837 		qfine = qfine_vol/dumvol1p
838 		qcoar = qcoar_vol/dumvol1p
839 		vtot_ofmode(m) =   &
840 			qfine_vol*fr_fine(m) + qcoar_vol*fr_coar(m)
841 	    else
842 		goto 2800
843 	    end if
844 	    l3 = l1
845 
846 	    if ((l1 .gt. 0) .and. (l1 .le. lmaxd) .and.   &
847 	        (l3 .ge. param_first_scalar)) then
848 		qval = qfine*fr_fine(m) + qcoar*fr_coar(m)
849 		rclm(1,l1) = qval
850 
851 		do it = its, ite
852 		do kt = kts, kte-1
853 		do jt = jts, jte
854 
855 !wig 28-May-2004, begin block...
856 ! Determine height multiplier...
857 ! This should mimic the calculation in sorgam_set_aer_bc_pnnl,
858 ! sorgam_init_aer_ic_pnnl, bdy_chem_value_mosaic
859 !!$!    Height(m)     Multiplier
860 !!$!    ---------     ----------
861 !!$!    <=2000        1.0
862 !!$!    2000<z<3000   linear transition zone to 0.5
863 !!$!    3000<z<5000   linear transision zone to 0.25
864 !!$!    >=5000        0.25
865 !!$!
866 !!$! which translates to:
867 !!$!    2000<z<3000   mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
868 !!$!    3000<z<5000   mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
869 !!$! or in reduced form:
870 !!$		   if( z(it,kt,jt) <= 2000. ) then
871 !!$		      mult = 1.0
872 !!$		   elseif( z(it,kt,jt) > 2000. &
873 !!$                        .and. z(it,kt,jt) <= 3000. ) then
874 !!$		      mult = 1.0 - 0.0005*(z(it,kt,jt)-2000.)
875 !!$		   elseif( z(it,kt,jt) > 3000. &
876 !!$                        .and. z(it,kt,jt) <= 5000. ) then
877 !!$		      mult = 0.5 - 1.25e-4*(z(it,kt,jt)-3000.)
878 !!$		   else
879 !!$		      mult = 0.25
880 !!$           end if
881 ! Updated 1-Apr-2005:
882 !    Height(m)     Multiplier
883 !    ---------     ----------
884 !    <=2000        1.0
885 !    2000<z<3000   linear transition zone to 0.25
886 !    3000<z<5000   linear transision zone to 0.125
887 !    >=5000        0.125
888 !
889 ! which translates to:
890 !    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
891 !    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
892 ! or in reduced form:
893 		   if( z(it,kt,jt) <= 2000. ) then
894 		      mult = 1.0
895 		   elseif( z(it,kt,jt) > 2000. &
896                         .and. z(it,kt,jt) <= 3000. ) then
897 		      mult = 1.0 - 0.00075*(z(it,kt,jt)-2000.)
898 		   elseif( z(it,kt,jt) > 3000. &
899                         .and. z(it,kt,jt) <= 5000. ) then
900 		      mult = 0.25 - 4.166666667e-5*(z(it,kt,jt)-3000.)
901 		   else
902 		      mult = 0.125
903            end if
904 
905 		    chem(it,kt,jt,l3) = mult*rclm(1,l1)
906 !wig 28-May-2004, ...end block
907 !		    chem(it,kt,jt,l3) = rclm(1,l1)
908 		end do
909 		end do
910 		end do
911 	    end if
912 
913 
914 2800   continue
915 2900   continue
916 
917 !
918 !   do diagnostic output
919 !
920 	dumarr(:) = 0.0
921 	msg = ' '
922 	call peg_message( lunout, msg )
923 	msg = '*** subr mosaic_init_wrf_mixrats_opt2 results'
924 	call peg_message( lunout, msg )
925 	msg = '    mass in ug/m3     number in #/m3     volume in cm3/m3'
926 	call peg_message( lunout, msg )
927 	msg = ' '
928 	call peg_message( lunout, msg )
929 	msg = ' mode  l  l1  species      conc'
930 	call peg_message( lunout, msg )
931 
932         do 3190 mdum = 1, nsize_aer(itype)+1
933 	m = min( mdum, nsize_aer(itype) )
934 	msg = ' '
935 	call peg_message( lunout, msg )
936         do 3150 l = 1, ncomp_plustracer_aer(itype)+4
937 
938             if (l .le. ncomp_plustracer_aer(itype)) then
939                 l1 = massptr_aer(l,m,itype,iphase)
940 		dumname = name_aer(l,itype)
941 		dum = rclm(1,l1)
942             else if (l .eq. ncomp_plustracer_aer(itype)+1) then
943                 l1 = hyswptr_aer(m,itype)
944 		dumname = 'hystwatr'
945 		dum = rclm(1,l1)
946             else if (l .eq. ncomp_plustracer_aer(itype)+2) then
947                 l1 = waterptr_aer(m,itype)
948 		dumname = 'water'
949 		dum = rclm(1,l1)
950             else if (l .eq. ncomp_plustracer_aer(itype)+3) then
951                 l1 = numptr_aer(m,itype,iphase)
952 		dumname = 'number'
953 		dum = rclm(1,l1)
954             else if (l .eq. ncomp_plustracer_aer(itype)+4) then
955                 l1 = 0
956 		dumname = 'volume'
957 		dum = vtot_ofmode(m)
958 	    else
959 		dumname = '????????'
960 		l1 = -1
961 		dum = -1.0
962             end if
963 
964 	    if (mdum .le. nsize_aer(itype)) then
965 		dumarr(l) = dumarr(l) + dum
966 		write(msg,9620) m, l, l1, dumname, dum
967 	    else
968 		write(msg,9625) l, dumname, dumarr(l)
969 	    end if
970 	    call peg_message( lunout, msg )
971 
972 3150   continue
973 3190   continue
974 
975 9620    format( 3i4, 2x, a, 3(1pe12.3) )
976 9625    format( ' sum', i4, '   -', 2x, a, 3(1pe12.3) )
977 
978 
979 !   all done
980 	return
981 	end subroutine mosaic_init_wrf_mixrats_opt2
982 
983 
984 !-----------------------------------------------------------------------
985 	real function erfc_num_recipes( x )
986 !
987 !   from press et al, numerical recipes, 1990, page 164
988 !
989 	implicit none
990 	real x
991 	double precision erfc_dbl, dum, t, zz
992 
993 	zz = abs(x)
994 	t = 1.0/(1.0 + 0.5*zz)
995 
996 !	erfc_num_recipes =
997 !     &	  t*exp( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
998 !     &	  t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
999 !     &	                                   t*(-1.13520398 +
1000 !     &	  t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
1001 
1002 	dum =  ( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +   &
1003       	  t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +   &
1004       	                                   t*(-1.13520398 +   &
1005       	  t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
1006 
1007 	erfc_dbl = t * exp(dum)
1008 	if (x .lt. 0.0) erfc_dbl = 2.0d0 - erfc_dbl
1009 
1010 	erfc_num_recipes = erfc_dbl
1011 
1012 	return
1013 	end function erfc_num_recipes     
1014 
1015 
1016 !-----------------------------------------------------------------------
1017 	end module module_mosaic_initmixrats
1018 
1019 
1020 
1021 
1022 !-----------------------------------------------------------------------
1023  	subroutine bdy_chem_value_mosaic ( chem_bv, alt, zz, nch, config_flags )
1024 !
1025 ! provides boundary values for the mosaic aerosol species
1026 !
1027 ! it is outside of the module declaration because of potential
1028 !    module1 --> module2 --> module1 use conflicts
1029 !
1030 ! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
1031 !     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now
1032 !
1033 	use module_configure, only:  grid_config_rec_type
1034 	use module_state_description, only:  param_first_scalar,   &
1035 			aer_bc_pnnl
1036 	use module_data_mosaic_asect
1037 	use module_data_mosaic_other
1038 	implicit none
1039 
1040 ! arguments
1041 	REAL,    intent(OUT)  :: chem_bv    ! boundary value for chem(-,-,-,nch)
1042 	REAL,    intent(IN)   :: alt        ! inverse density
1043 	REAL,    intent(IN)   :: zz         ! height
1044 	INTEGER, intent(IN)   :: nch        ! index number of chemical species
1045 	TYPE (grid_config_rec_type), intent(in) :: config_flags
1046 
1047 ! local variables
1048 	integer iphase, itype, m
1049 
1050 	real, parameter :: chem_bv_def = 1.0e-20
1051 !wig 10-May-2004, added mult
1052 	real dumvol1p, mult
1053 	real qcoar, qfine, qval
1054 
1055 	real :: fr_coar(8), fr_fine(8)
1056 
1057 !wig 1-Apr-2005,  Updated fractional breakdown between bins. (See also
1058 !                 mosaic_init_wrf_mixrats_opt2, above, and mosaic_addemiss
1059 !                 in module_mosaic_addemiss.F). Note that these values no
1060 !                 longer match those in mosaic_addemiss.
1061 	real, save :: fr8b_coar(8) =   &
1062         (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.933, 0.067 /)
1063 !		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.16, 0.84 /)
1064 	real, save :: fr8b_fine(8) =   &
1065         (/ 0.006, 0.024, 0.231, 0.341, 0.140, 0.258, 0., 0./) !5-Apr-2005 values
1066 !       (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
1067 !		(/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
1068 !		(/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.00, 0.00 /)
1069 
1070 	real :: qfine_so4, qfine_no3, qfine_cl, qfine_msa,   &
1071 		qfine_co3, qfine_nh4, qfine_na, qfine_ca, qfine_oin,   &
1072 		qfine_oc, qfine_bc, qfine_hysw, qfine_watr, qfine_vol
1073 	real :: qcoar_so4, qcoar_no3, qcoar_cl, qcoar_msa,   &
1074 		qcoar_co3, qcoar_nh4, qcoar_na, qcoar_ca, qcoar_oin,   &
1075 		qcoar_oc, qcoar_bc, qcoar_hysw, qcoar_watr, qcoar_vol
1076 
1077 	character*80 msg
1078 
1079 
1080 
1081 ! when aer_bc_opt /= aer_bc_pnnl,
1082 ! just chem_bv=near-zero for all species
1083 	chem_bv = chem_bv_def
1084 	if (config_flags%aer_bc_opt /= aer_bc_pnnl) return
1085 	if (nch < param_first_scalar) return
1086 
1087 
1088 !   *** currently only works for ntype_aer = 1
1089 	itype = 1
1090 	iphase = ai_phase
1091 	m = 1
1092 
1093 
1094 !
1095 !   following is for aer_bc_opt == aer_bc_pnnl
1096 !       when boundary values are set for Texas August 2000 simulations
1097 !
1098 !   set the partitioning fractions for either 8 or 4 bins
1099 	if (nsize_aer(itype) == 8) then
1100 	    fr_coar(:) = fr8b_coar(:)
1101 	    fr_fine(:) = fr8b_fine(:)
1102 	else if (nsize_aer(itype) == 4) then
1103 	    do m = 1, nsize_aer(itype)
1104 		fr_coar(m) = fr8b_coar(2*m) + fr8b_coar(2*m-1)
1105 		fr_fine(m) = fr8b_fine(2*m) + fr8b_fine(2*m-1)
1106 	    end do
1107 	else
1108 	    write(msg,'(a,i5)')   &
1109 		'subr bdy_chem_value_mosaic' //   &
1110 		' - nsize_aer(itype) must be 4 or 8 but = ', nsize_aer(itype)
1111 	    call wrf_error_fatal( msg )
1112 	end if
1113 
1114 ! Determine height multiplier...
1115 ! This should mimic the calculation in sorgam_set_aer_bc_pnnl,
1116 ! sorgam_init_aer_ic_pnnl, and mosaic_init_wrf_mixrats_opt2
1117 !!$!    Height(m)     Multiplier
1118 !!$!    ---------     ----------
1119 !!$!    <=2000        1.0
1120 !!$!    2000<z<3000   linear transition zone to 0.5
1121 !!$!    3000<z<5000   linear transision zone to 0.25
1122 !!$!    >=5000        0.25
1123 !!$!
1124 !!$! which translates to:
1125 !!$!    2000<z<3000   mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
1126 !!$!    3000<z<5000   mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
1127 !!$! or in reduced form:
1128 !!$        if( zz <= 2000. ) then
1129 !!$           mult = 1.0
1130 !!$        elseif( zz > 2000. &
1131 !!$             .and. zz <= 3000. ) then
1132 !!$           mult = 1.0 - 0.0005*(zz-2000.)
1133 !!$        elseif( zz > 3000. &
1134 !!$             .and. zz <= 5000. ) then
1135 !!$           mult = 0.5 - 1.25e-4*(zz-3000.)
1136 !!$        else
1137 !!$           mult = 0.25
1138 !!$        end if
1139 ! Updated aerosol profile multiplier 1-Apr-2005:
1140 !    Height(m)     Multiplier
1141 !    ---------     ----------
1142 !    <=2000        1.0
1143 !    2000<z<3000   linear transition zone to 0.25
1144 !    3000<z<5000   linear transision zone to 0.125
1145 !    >=5000        0.125
1146 !
1147 ! which translates to:
1148 !    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
1149 !    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
1150 ! or in reduced form:
1151         if( zz <= 2000. ) then
1152            mult = 1.0
1153         elseif( zz > 2000. &
1154              .and. zz <= 3000. ) then
1155            mult = 1.0 - 0.00075*(zz-2000.)
1156         elseif( zz > 3000. &
1157              .and. zz <= 5000. ) then
1158            mult = 0.25 - 4.166666667e-5*(zz-3000.)
1159         else
1160            mult = 0.125
1161         end if
1162 
1163 !   Set values for fine and coarse mass concentrations, and
1164 !   compute fine/coarse volume concentrations. These are also set in
1165 !   mosaic_init_wrf_mixrats_opt2.
1166 !   (these are latest values from 1-Apr-2005 discussions)
1167 !wig 10-May-2004, added height multiplier (mult*) to each species...
1168 	qfine_so4 = mult*2.14
1169 	qcoar_so4 = mult*0.242
1170 	qfine_no3 = mult*0.11
1171 	qcoar_no3 = mult*0.03
1172 	qfine_cl  = mult*0.3
1173 	qcoar_cl  = mult*0.139
1174 	qfine_msa = mult*0.0
1175 	qcoar_msa = mult*0.0
1176 	qfine_co3 = mult*0.0
1177 	qcoar_co3 = mult*0.0
1178 	qfine_nh4 = mult*0.83
1179 	qcoar_nh4 = mult*0.10
1180 	qfine_na  = mult*0.2
1181 	qcoar_na  = mult*0.09
1182 	qfine_ca  = mult*0.0
1183 	qcoar_ca  = mult*0.0
1184 	qfine_oin = mult*6.2
1185 	qcoar_oin = mult*0.35
1186 	qfine_oc  = mult*1.00
1187 	qcoar_oc  = mult*1.50
1188 	qfine_bc  = mult*0.2
1189 	qcoar_bc  = mult*0.075
1190 	qfine_hysw = mult*0.0
1191 	qcoar_hysw = mult*0.0
1192 	qfine_watr = mult*0.0
1193 	qcoar_watr = mult*0.0
1194 !!$! old values from 23-Apr-2004:
1195 !!$	qfine_so4 = mult*2.554
1196 !!$	qcoar_so4 = mult*0.242
1197 !!$	qfine_no3 = mult*0.07
1198 !!$	qcoar_no3 = mult*0.03
1199 !!$	qfine_cl  = mult*0.324
1200 !!$	qcoar_cl  = mult*0.139
1201 !!$	qfine_msa = mult*0.0
1202 !!$	qcoar_msa = mult*0.0
1203 !!$	qfine_co3 = mult*0.0
1204 !!$	qcoar_co3 = mult*0.0
1205 !!$	qfine_nh4 = mult*0.98
1206 !!$	qcoar_nh4 = mult*0.10
1207 !!$	qfine_na  = mult*0.21
1208 !!$	qcoar_na  = mult*0.09
1209 !!$	qfine_ca  = mult*0.0
1210 !!$	qcoar_ca  = mult*0.0
1211 !!$	qfine_oin = mult*0.15
1212 !!$	qcoar_oin = mult*0.35
1213 !!$	qfine_oc  = mult*1.00
1214 !!$	qcoar_oc  = mult*1.50
1215 !!$	qfine_bc  = mult*0.175
1216 !!$	qcoar_bc  = mult*0.075
1217 !!$	qfine_hysw = mult*0.0
1218 !!$	qcoar_hysw = mult*0.0
1219 !!$	qfine_watr = mult*0.0
1220 !!$	qcoar_watr = mult*0.0
1221 
1222 
1223 ! qfine_so4 ... are ug/m3 so 1.0e-6 factor gives g/m3
1224 ! dens_so4  ... are g/cm3;  qfine_vol ... are cm3/m3
1225 	qfine_vol = 1.0e-6 * (   &
1226 	    (qfine_so4/dens_so4_aer) + (qfine_no3/dens_no3_aer) +   &
1227 	    (qfine_cl /dens_cl_aer ) + (qfine_msa/dens_msa_aer) +   &
1228 	    (qfine_co3/dens_co3_aer) + (qfine_nh4/dens_nh4_aer) +   &
1229 	    (qfine_na /dens_na_aer ) + (qfine_ca /dens_ca_aer ) +   &
1230 	    (qfine_oin/dens_oin_aer) + (qfine_oc /dens_oc_aer ) +   &
1231 	    (qfine_bc /dens_bc_aer ) )
1232 	qcoar_vol = 1.0e-6 * (   &
1233 	    (qcoar_so4/dens_so4_aer) + (qcoar_no3/dens_no3_aer) +   &
1234 	    (qcoar_cl /dens_cl_aer ) + (qcoar_msa/dens_msa_aer) +   &
1235 	    (qcoar_co3/dens_co3_aer) + (qcoar_nh4/dens_nh4_aer) +   &
1236 	    (qcoar_na /dens_na_aer ) + (qcoar_ca /dens_ca_aer ) +   &
1237 	    (qcoar_oin/dens_oin_aer) + (qcoar_oc /dens_oc_aer ) +   &
1238 	    (qcoar_bc /dens_bc_aer ) )
1239 
1240 	qfine = -1.0e30
1241 	qcoar = -1.0e30
1242 
1243 !   identify the species by checking "nch" against the "lptr_xxx_a_amode(m)"
1244         do 2900 m = 1, nsize_aer(itype)
1245 
1246 	    if (nch .eq. lptr_so4_aer(m,itype,iphase)) then
1247 		qfine = qfine_so4
1248 		qcoar = qcoar_so4
1249 	    else if (nch .eq. lptr_no3_aer(m,itype,iphase)) then
1250 		qfine = qfine_no3
1251 		qcoar = qcoar_no3
1252 	    else if (nch .eq. lptr_cl_aer(m,itype,iphase)) then
1253 		qfine = qfine_cl
1254 		qcoar = qcoar_cl
1255 	    else if (nch .eq. lptr_msa_aer(m,itype,iphase)) then
1256 		qfine = qfine_msa
1257 		qcoar = qcoar_msa
1258 	    else if (nch .eq. lptr_co3_aer(m,itype,iphase)) then
1259 		qfine = qfine_co3
1260 		qcoar = qcoar_co3
1261 	    else if (nch .eq. lptr_nh4_aer(m,itype,iphase)) then
1262 		qfine = qfine_nh4
1263 		qcoar = qcoar_nh4
1264 	    else if (nch .eq. lptr_na_aer(m,itype,iphase)) then
1265 		qfine = qfine_na
1266 		qcoar = qcoar_na
1267 	    else if (nch .eq. lptr_ca_aer(m,itype,iphase)) then
1268 		qfine = qfine_ca
1269 		qcoar = qcoar_ca
1270 	    else if (nch .eq. lptr_oin_aer(m,itype,iphase)) then
1271 		qfine = qfine_oin
1272 		qcoar = qcoar_oin
1273 	    else if (nch .eq. lptr_oc_aer(m,itype,iphase)) then
1274 		qfine = qfine_oc
1275 		qcoar = qcoar_oc
1276 	    else if (nch .eq. lptr_bc_aer(m,itype,iphase)) then
1277 		qfine = qfine_bc
1278 		qcoar = qcoar_bc
1279 	    else if (nch .eq. hyswptr_aer(m,itype)) then
1280 		qfine = qfine_hysw
1281 		qcoar = qcoar_hysw
1282 	    else if (nch .eq. waterptr_aer(m,itype)) then
1283 		qfine = qfine_watr
1284 		qcoar = qcoar_watr
1285 	    else if (nch .eq. numptr_aer(m,itype,iphase)) then
1286 		dumvol1p = sqrt(volumlo_sect(m,itype)*volumhi_sect(m,itype))
1287 		qfine = qfine_vol/dumvol1p
1288 		qcoar = qcoar_vol/dumvol1p
1289 	    end if
1290 
1291 	    if ((qfine >= 0.0) .and. (qcoar >= 0.0)) then
1292 		qval = qfine*fr_fine(m) + qcoar*fr_coar(m)
1293 		chem_bv = qval*alt
1294 		goto 2910
1295 	    end if
1296 
1297 2900   continue
1298 2910   continue
1299 
1300 	return
1301  	end subroutine bdy_chem_value_mosaic
1302 
1303