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