Merge branch 'master' into devel
[wrffire.git] / wrfv2_fire / chem / module_mosaic_initmixrats.F
blob167e86e20bc1b97a0d8e7368acdaac2a13aa1d95
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.
7 ! MOSAIC module: see module_mosaic_driver.F for references 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 !             2   Mexico City
15 !             3   Mexico City for Testbed, to be consistent with MADE/SORGAM too
16 !             4   bc,ic values from wrfinput and wrfbdy files
17 #define CASENAME 0
19         module module_mosaic_initmixrats
21 ! rce 2005-feb-18 - several fixes for indices of dlo/hi_sect, volumlo/hi_sect,
22 !     which are now (isize,itype)
24 ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
25 !     variables in module_data_mosaic_asect
27         USE module_state_description
29         integer, parameter :: mosaic_init_wrf_mixrats_flagaa = 1
30                             ! turns subr mosaic_init_wrf_mixrats on/off
32         contains
35 !-----------------------------------------------------------------------
36         subroutine mosaic_init_wrf_mixrats(  &
37                 iflagaa, config_flags,       &
38                 chem, alt, z_at_w, g,        &
39                 ids,ide, jds,jde, kds,kde,   &
40                 ims,ime, jms,jme, kms,kme,   &
41                 its,ite, jts,jte, kts,kte    )
44 !   initializes the species and number mixing ratios for each section
46 !   this top level routine simply calls other routines depending on value
47 !       of config_flags%aer_ic_opt
49         use module_configure, only:  grid_config_rec_type
50         use module_state_description, only:  num_chem, param_first_scalar,   &
51                         aer_ic_pnnl
52         use module_data_mosaic_asect
53         use module_data_mosaic_other
54         use module_peg_util, only:  peg_message, peg_error_fatal
56         implicit none
59 !   subr arguments
60         type(grid_config_rec_type), intent(in) :: config_flags
62         integer, intent(in) ::   &
63                 iflagaa,   &
64                 ids, ide, jds, jde, kds, kde,   &
65                 ims, ime, jms, jme, kms, kme,   &
66                 its, ite, jts, jte, kts, kte
68         real, intent(inout),   &
69                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
70                 chem
72         real, intent(in),      &
73                 dimension( ims:ime, kms:kme, jms:jme ) :: &
74         alt, z_at_w
75     real, intent(in) :: g
77 !   local variables
78         integer :: i, ic, j, k
80 #if (CASENAME == 0) || (CASENAME == 1) || (CASENAME == 2) || (CASENAME == 3)
81         if (config_flags%aer_ic_opt == aer_ic_pnnl) then
82             call mosaic_init_wrf_mixrats_opt2(   &
83                 iflagaa, config_flags,           &
84                 chem, z_at_w, g,                 &
85                 ids,ide, jds,jde, kds,kde,       &
86                 ims,ime, jms,jme, kms,kme,       &
87                 its,ite, jts,jte, kts,kte        )
89         else
90             call mosaic_init_wrf_mixrats_opt1(   &
91                 iflagaa, config_flags,   &
92                 chem,   &
93                 ids,ide, jds,jde, kds,kde,   &
94                 ims,ime, jms,jme, kms,kme,   &
95                 its,ite, jts,jte, kts,kte    )
97         end if
99 ! Aerosol species are returned from above in concentration units. Convert
100 ! them to mixing ratio for use in advection, etc.
101     do ic = p_so4_a01,num_chem
102        do j = jts,jte
103           do k = kts,kte-1
104              do i = its,ite
105                 chem(i,k,j,ic) = chem(i,k,j,ic)*alt(i,k,j)
106              end do
107           end do
108        end do
109     end do
111 ! Fill the top z-staggered location to prevent trouble during advection.
112     do ic = p_so4_a01,num_chem
113        do j = jts,jte
114           do i = its,ite
115              chem(i,kte,j,ic) = chem(i,kte-1,j,ic)
116           end do
117        end do
118     end do
119 #endif
121         return
122         end subroutine mosaic_init_wrf_mixrats
125 !-----------------------------------------------------------------------
126         subroutine mosaic_init_wrf_mixrats_opt1(   &
127                 iflagaa, config_flags,   &
128                 chem,   &
129                 ids,ide, jds,jde, kds,kde,   &
130                 ims,ime, jms,jme, kms,kme,   &
131                 its,ite, jts,jte, kts,kte    )
134 !   initializes the species and number mixing ratios for each section
135 !   based on user-specified lognormal modes that span the size distribution
137 ! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
138 !     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now; 
139 !     added l1dum
141         use module_configure, only:  grid_config_rec_type
142         use module_state_description, only:  num_chem, param_first_scalar
143         use module_data_mosaic_asect
144         use module_data_mosaic_other
145         use module_peg_util, only:  peg_message, peg_error_fatal
147         implicit none
149 !   subr arguments
150         type(grid_config_rec_type), intent(in) :: config_flags
152         integer, intent(in) ::   &
153                 iflagaa,   &
154                 ids, ide, jds, jde, kds, kde,   &
155                 ims, ime, jms, jme, kms, kme,   &
156                 its, ite, jts, jte, kts, kte
158         real, intent(inout),   &
159                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
160                 chem
162 !   local variables
163         integer i, j, k, l, ll, l1, l3, l1dum, m, mdum, nsm
164         integer it, jt, kt
166         real dum, dumdp, dumrsfc, dumvol,   &
167           xlo, xhi,   &
168           dumvol1p,  &
169           pdummb, zdumkm, zscalekm, zfactor
171         real vtot_nsm_ofmode(maxd_asize)
172         real dumarr(maxd_acomp+5)
174         real erfc
176 !       double precision fracnum, fracvol, tlo, thi
177         real fracvol, tlo, thi
179         integer nmaxd_nsm
180         parameter (nmaxd_nsm = 4)
182         integer iphase, itype, ntot_nsm
183         integer iiprof_nsm(nmaxd_nsm)
184         integer lldum_so4, lldum_nh4, lldum_oc, lldum_bc,   &
185                 lldum_oin, lldum_na, lldum_cl, lldum_hysw
187         real sx_nsm(nmaxd_nsm), sxr2_nsm(nmaxd_nsm),   &
188           x0_nsm(nmaxd_nsm), x3_nsm(nmaxd_nsm),   &
189           rtot_nsm(maxd_acomp,nmaxd_nsm),   &
190           vtot_nsm(nmaxd_nsm), xntot_nsm(nmaxd_nsm)
192         real dgnum_nsm(nmaxd_nsm), sigmag_nsm(nmaxd_nsm)
193         real aaprof_nsm(maxd_acomp+1,nmaxd_nsm)
194         real bbprof_nsm(nmaxd_nsm)
196         character*80 msg
197         character*10 dumname
200 !   check for on/off
201         if (mosaic_init_wrf_mixrats_flagaa .le. 0) return
204 !   *** currently only works for ntype_aer = 1
205         itype = 1
206         iphase = ai_phase
207         m = 1
209 !   set values for initialization modes
210         iiprof_nsm(:) = 1
211         aaprof_nsm(:,:) = 0.0
212         bbprof_nsm(:) = 0.0
214         ntot_nsm = 4
215         ntot_nsm = min( ntot_nsm, nsize_aer(itype) )
217         lldum_so4 = 0
218         lldum_nh4 = 0
219         lldum_oc  = 0
220         lldum_bc  = 0
221         lldum_oin = 0
222         lldum_na  = 0
223         lldum_cl  = 0
224         lldum_hysw = 0
226         do ll = 1, ncomp_plustracer_aer(itype)
227             if (massptr_aer(ll,m,itype,iphase) .eq.   &
228                 lptr_so4_aer(m,itype,iphase)) lldum_so4 = ll
229             if (massptr_aer(ll,m,itype,iphase) .eq.   &
230                 lptr_nh4_aer(m,itype,iphase)) lldum_nh4 = ll
231             if (massptr_aer(ll,m,itype,iphase) .eq.   &
232                 lptr_oc_aer(m,itype,iphase))  lldum_oc  = ll
233             if (massptr_aer(ll,m,itype,iphase) .eq.   &
234                 lptr_bc_aer(m,itype,iphase))  lldum_bc  = ll
235             if (massptr_aer(ll,m,itype,iphase) .eq.   &
236                 lptr_oin_aer(m,itype,iphase)) lldum_oin = ll
237             if (massptr_aer(ll,m,itype,iphase) .eq.   &
238                 lptr_na_aer(m,itype,iphase))  lldum_na  = ll
239             if (massptr_aer(ll,m,itype,iphase) .eq.   &
240                 lptr_cl_aer(m,itype,iphase))  lldum_cl  = ll
241         end do
242         if (hyswptr_aer(m,itype) .gt. 0)   &
243                 lldum_hysw = ncomp_plustracer_aer(itype) + 1
245         msg = ' '
246         if (lldum_so4 .le. 0)   &
247                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_so4 = 0'
248         if (lldum_nh4 .le. 0)   &
249                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_nh4 = 0'
250         if (lldum_oc .le. 0)   &
251                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_oc = 0'
252         if (lldum_bc .le. 0)   &
253                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_bc = 0'
254         if (lldum_oin .le. 0)   &
255                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_oin = 0'
256         if (lldum_na .le. 0)   &
257                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_na = 0'
258         if (lldum_cl .le. 0)   &
259                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_cl = 0'
260         if (lldum_hysw .le. 0)   &
261                 msg = '*** subr mosaic_init_wrf_mixrats - lldum_hysw = 0'
262         if (msg .ne. ' ') call peg_error_fatal( lunerr, msg )
265         do nsm = 1, ntot_nsm
267         if (nsm .eq. 1) then
268 !   accum mode with so4, nh4, oc, bc
269             dgnum_nsm( nsm) = 0.15e-4
270             sigmag_nsm(nsm) = 2.0
271             aaprof_nsm(lldum_so4,nsm) = 0.50
272             aaprof_nsm(lldum_nh4,nsm) = aaprof_nsm(lldum_so4,nsm)   &
273                                         * (mw_nh4_aer/mw_so4_aer)
274             aaprof_nsm(lldum_oc,nsm) = 0.25
275             aaprof_nsm(lldum_bc,nsm) = 0.05
276             aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_so4,nsm) * 0.2
278         else if (nsm .eq. 2) then
279 !   aitken mode with so4, nh4, oc, bc
280             dgnum_nsm( nsm) = 0.03e-4
281             sigmag_nsm(nsm) = 2.0
282             aaprof_nsm(lldum_so4,nsm) = 0.50 * 0.020
283             aaprof_nsm(lldum_nh4,nsm) = aaprof_nsm(lldum_so4,nsm)   &
284                                         * (mw_nh4_aer/mw_so4_aer)
285             aaprof_nsm(lldum_oc,nsm) = 0.25 * 0.020
286             aaprof_nsm(lldum_bc,nsm) = 0.05 * 0.020
287             aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_so4,nsm) * 0.2
289         else if (nsm .eq. 3) then
290 !   coarse dust mode with oin
291             dgnum_nsm( nsm) = 1.0e-4
292             sigmag_nsm(nsm) = 2.0
293             aaprof_nsm(lldum_oin,nsm) = 0.5
294             aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm( 9,nsm) * 1.0e-3
296         else if (nsm .eq. 4) then
297 !   coarse seasalt mode with na, cl
298             dgnum_nsm( nsm) = 2.0e-4
299             sigmag_nsm(nsm) = 2.0
300             aaprof_nsm(lldum_cl,nsm) = 0.1
301             aaprof_nsm(lldum_na,nsm) = aaprof_nsm(lldum_cl,nsm)   &
302                                         * (mw_na_aer/mw_cl_aer)
303             aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_cl,nsm) * 0.2
305         end if
307         end do
309 !   when iflagaa = 1/2/3/4, use only the nsm-mode = iflagaa
310         if (iflagaa .gt. 0) then
311             do nsm = 1, ntot_nsm
312                 if (nsm .ne. iflagaa) aaprof_nsm(:,nsm) = 0.0
313             end do
314         end if
319 !   do the initialization now
322 !   calculate mode parameters
323         do nsm = 1, ntot_nsm
324             sx_nsm(nsm) = alog( sigmag_nsm(nsm) )
325             sxr2_nsm(nsm) = sx_nsm(nsm) * sqrt(2.0)
326             x0_nsm(nsm) = alog( dgnum_nsm(nsm) )
327             x3_nsm(nsm) = x0_nsm(nsm) + 3.0*sx_nsm(nsm)*sx_nsm(nsm)
328         end do
330 !   initialize rclm array to zero
331         rclm(:,:) = 0.
332 !       rclmsvaa(:,:) = 0.
335 !   loop over all vertical levels
337 !       do 12900 k = 1, ktot
338         do 12900 k = 1, 1
340 !       pdummb = 1013.*scord(k)
341 !       zdumkm = ptoz( pdummb ) * 1.e-3
342         zdumkm = 0.0
345 !   for each species and nsm mode, define total mixing ratio
346 !       (for all sizes) at level k
348 !   iiprof_nsm =  +1 gives constant mixing ratio
349 !       aaprof_nsm(l,nsm) = constant mixing ratio (ppb)
350 !   iiprof_nsm =  +2 gives exponential profile
351 !       aaprof_nsm(l,nsm) = surface mixing ratio (ppb)
352 !       bbprof_nsm(l) = scale height (km)
353 !   iiprof_nsm =  +3 gives linear profile (then zero at z > zscalekm)
354 !       aaprof_nsm(l,nsm) = surface mixing ratio (ppb)
355 !       bbprof_nsm(l) = height (km) at which mixing ratio = 0
357         do nsm = 1, ntot_nsm
359         if (iiprof_nsm(nsm) .eq. 2) then
360             zscalekm = bbprof_nsm(nsm)
361             zfactor = exp( -zdumkm/zscalekm )
362         else if (iiprof_nsm(nsm) .eq. 3) then
363             zscalekm = bbprof_nsm(nsm)
364             zfactor = max( 0., (1. - zdumkm/zscalekm) )
365         else
366             zfactor = 1.0
367         end if
369             do ll = 1, ncomp_plustracer_aer(itype) + 1
370                 rtot_nsm(ll,nsm) = aaprof_nsm(ll,nsm) * zfactor
371             end do
373         end do
375 !   compute total volume and number mixing ratios for each nsm mode
376 !   rtot_nsm is ug/m3,  1.0e-6*rtot is g/m3,  vtot_nsm is cm3/m3
377         do nsm = 1, ntot_nsm
378             dumvol = 0.
379             do ll = 1, ncomp_aer(itype)
380                 dum = 1.0e-6*rtot_nsm(ll,nsm)/dens_aer(ll,itype)
381                 dumvol = dumvol + max( 0., dum )
382             end do
383             vtot_nsm(nsm) = dumvol
384         end do
386 !   now compute species and number mixing ratios for each bin
387         do 12700 m = 1, nsize_aer(itype)
389         vtot_nsm_ofmode(m) = 0.0
391         do 12500 nsm = 1, ntot_nsm
393 !   for nsm_mode = n, compute fraction of number and volume
394 !   that is in bin m
395         xlo = alog( dlo_sect(m,itype) )
396         xhi = alog( dhi_sect(m,itype) )
398         tlo = (xlo - x3_nsm(nsm))/sxr2_nsm(nsm)
399         thi = (xhi - x3_nsm(nsm))/sxr2_nsm(nsm)
400         if (tlo .ge. 0.) then
401 !           fracvol = 0.5*( erfc(tlo) - erfc(thi) )
402             fracvol = 0.5*( erfc_num_recipes(tlo) - erfc_num_recipes(thi) )
403         else
404 !           fracvol = 0.5*( erfc(-thi) - erfc(-tlo) )
405             fracvol = 0.5*( erfc_num_recipes(-thi) - erfc_num_recipes(-tlo) )
406         end if
407         fracvol = max( fracvol, 0.0 )
409         vtot_nsm_ofmode(m) = vtot_nsm_ofmode(m)  + vtot_nsm(nsm)*fracvol
411 !   now load that fraction of species-mixing-ratio
412 !   into the appropriate rclm location
413         do ll = 1, ncomp_plustracer_aer(itype)
414             rclm( k, massptr_aer(ll,m,itype,iphase) ) =   &
415             rclm( k, massptr_aer(ll,m,itype,iphase) ) +   &
416                 fracvol*rtot_nsm(ll,nsm)
417         end do
419         if ((iphase .eq. ai_phase) .and.   &
420             (lldum_hysw .gt. 0) .and.   &
421             (hyswptr_aer(m,itype) .gt. 0)) then
423             rclm( k, hyswptr_aer(m,itype) ) =   &
424             rclm( k, hyswptr_aer(m,itype) ) +   &
425                 fracvol*rtot_nsm(lldum_hysw,nsm)
426         end if
428 12500   continue
430 !   now compute number from volume
431         dum = sqrt( dlo_sect(m,itype)*dhi_sect(m,itype) )
432         dumvol1p = (pi/6.0)*(dum**3)
433         rclm( k, numptr_aer(m,itype,iphase) ) = vtot_nsm_ofmode(m)/dumvol1p
435 !   set water = hyswatr
436         if ((iphase .eq. ai_phase) .and.   &
437             (lldum_hysw .gt. 0) .and.   &
438             (hyswptr_aer(m,itype) .gt. 0) .and.   &
439             (waterptr_aer(m,itype) .gt. 0)) then
441             rclm( k, waterptr_aer(m,itype) ) =    &
442             rclm( k, hyswptr_aer(m,itype) )
443         end if
445 12700   continue
447 12900   continue
451 !   do diagnostic output
454 ! temporary
455 ! temporary
457         dumarr(:) = 0.0
458         msg = ' '
459         call peg_message( lunout, msg )
460         msg = '*** subr mosaic_init_wrf_mixrats_opt1 results'
461         call peg_message( lunout, msg )
462         msg = '    mass in ug/m3     number in #/m3     volume in cm3/m3'
463         call peg_message( lunout, msg )
464         msg = ' '
465         call peg_message( lunout, msg )
466         msg = ' mode  l  l1  species      conc'
467         call peg_message( lunout, msg )
469         do 14390 mdum = 1, nsize_aer(itype)+1
470         m = min( mdum, nsize_aer(itype) )
471         msg = ' '
472         call peg_message( lunout, msg )
473         do 14350 l = 1, ncomp_plustracer_aer(itype)+4
475             if (l .le. ncomp_plustracer_aer(itype)) then
476                 l1 = massptr_aer(l,m,itype,iphase)
477                 dumname = name_aer(l,itype)
478                 dum = rclm(1,l1)
479             else if (l .eq. ncomp_plustracer_aer(itype)+1) then
480                 l1 = hyswptr_aer(m,itype)
481                 dumname = 'hystwatr'
482                 dum = rclm(1,l1)
483             else if (l .eq. ncomp_plustracer_aer(itype)+2) then
484                 l1 = waterptr_aer(m,itype)
485                 dumname = 'water'
486                 dum = rclm(1,l1)
487             else if (l .eq. ncomp_plustracer_aer(itype)+3) then
488                 l1 = numptr_aer(m,itype,iphase)
489                 dumname = 'number'
490                 dum = rclm(1,l1)
491             else if (l .eq. ncomp_plustracer_aer(itype)+4) then
492                 l1 = 0
493                 dumname = 'volume'
494                 dum = vtot_nsm_ofmode(m)
495             else
496                 dumname = '=BADBAD='
497                 l1 = -1
498                 dum = -1.0
499             end if
501             l1dum = l1
502             if (aboxtest_testmode .gt. 0) l1dum = max( l1-1, 0 )
504             if (mdum .le. nsize_aer(itype)) then
505                 dumarr(l) = dumarr(l) + dum
506                 write(msg,9620) m, l, l1dum, dumname, dum
507             else
508                 write(msg,9625) l, dumname, dumarr(l)
509             end if
510             call peg_message( lunout, msg )
512 14350   continue
513 14390   continue
515 9620    format( 3i4, 2x, a, 3(1pe12.3) )
516 9625    format( ' sum', i4, '   -', 2x, a, 3(1pe12.3) )
520 !   load the chem array
522         do 16390 m = 1, nsize_aer(itype)
523         do 16350 l = 1, 15
525             if (l .eq. 1) then
526                 l1 = lptr_so4_aer(m,itype,iphase)
527             else if (l .eq. 2) then
528                 l1 = lptr_no3_aer(m,itype,iphase)
529             else if (l .eq. 3) then
530                 l1 = lptr_cl_aer(m,itype,iphase)
531             else if (l .eq. 4) then
532                 l1 = lptr_msa_aer(m,itype,iphase)
533             else if (l .eq. 5) then
534                 l1 = lptr_co3_aer(m,itype,iphase)
535             else if (l .eq. 6) then
536                 l1 = lptr_nh4_aer(m,itype,iphase)
537             else if (l .eq. 7) then
538                 l1 = lptr_na_aer(m,itype,iphase)
539             else if (l .eq. 8) then
540                 l1 = lptr_ca_aer(m,itype,iphase)
541             else if (l .eq. 9) then
542                 l1 = lptr_oin_aer(m,itype,iphase)
543             else if (l .eq. 10) then
544                 l1 = lptr_oc_aer(m,itype,iphase)
545             else if (l .eq. 11) then
546                 l1 = lptr_bc_aer(m,itype,iphase)
547             else if (l .eq. 12) then
548                 l1 = hyswptr_aer(m,itype)
549             else if (l .eq. 13) then
550                 l1 = waterptr_aer(m,itype)
551             else if (l .eq. 14) then
552                 l1 = numptr_aer(m,itype,iphase)
553             else
554                 goto 16350
555             end if
556             l3 = l1
558             if ((l1 .gt. 0) .and. (l1 .le. lmaxd) .and.   &
559                 (l3 .ge. param_first_scalar)) then
560                 do it = its, ite
561 !   *** note:  not sure what the kt limits should be here 
562                 do kt = kts, kte-1
563                 do jt = jts, jte
564                     chem(it,kt,jt,l3) = rclm(1,l1)
565                 end do
566                 end do
567                 end do
568             end if
570 16350   continue
571 16390   continue
574 !   all done
575         return
576         end subroutine mosaic_init_wrf_mixrats_opt1
579 !-----------------------------------------------------------------------
580 !wig 10-May-2004, added phb, ph, g
581         subroutine mosaic_init_wrf_mixrats_opt2(   &
582                 iflagaa, config_flags,             &
583                 chem, z_at_w, g,                   &
584                 ids,ide, jds,jde, kds,kde,         &
585                 ims,ime, jms,jme, kms,kme,         &
586                 its,ite, jts,jte, kts,kte          )
589 !   provides initial values for mosaic aerosol species (mass and number
590 !       mixing ratio) for "Texas August 2000" simulations
591 !   modified to use different profiles for different aerosols for NEAQS case, egc 7/2005
592 !   currently all the initial values are uniform in x, y, AND z
594 ! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
595 !     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now
597         use module_configure, only:  grid_config_rec_type
598         use module_state_description, only:  num_chem, param_first_scalar
599         use module_data_mosaic_asect
600         use module_data_mosaic_other
601         use module_peg_util, only:  peg_message, peg_error_fatal
603         implicit none
605 !   subr arguments
606         type(grid_config_rec_type), intent(in) :: config_flags
608         integer, intent(in) ::   &
609                 iflagaa,   &
610                 ids, ide, jds, jde, kds, kde,   &
611                 ims, ime, jms, jme, kms, kme,   &
612                 its, ite, jts, jte, kts, kte
614         real, intent(inout),     &
615                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
616                 chem
618         real, intent(in),        &
619                 dimension( ims:ime, kms:kme, jms:jme ) :: &
620                 z_at_w
621         real :: g
623 !   local variables
624         integer l, l1, l3, m, mdum
625         integer iphase, itype
626         integer it, jt, kt
628 !wig 10-May-2004, added mult
629         real dum, dumvol1p, mult
630         real qcoar, qfine, qval
632         real :: vtot_ofmode(maxd_asize)
633         real :: dumarr(maxd_acomp+5)
634         real :: fr_coar(8), fr_fine(8)
636 !wig 01-Apr-2005, Updated fractional breakdown between bins. (See also
637 !                 bdy_chem_value_mosaic in this file and mosaic_addemiss in 
638 !                 module_mosaic_addemiss.F) Note that the values here no
639 !                 longer match those in mosaic_addemiss.
640 !rce 10-May-2005, changed fr8b_coar to values determined by jdf
641         real, save :: fr8b_coar(8) =   &
642         (/ 0.0,  0.0,  0.0,  0.019,  0.0212,  0.1101,  0.2751, 0.7018 /) ! jdf Testbed
643 !jdf    (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.300, 0.700 /) ! 10-May-2005
644 !       (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.933, 0.067 /) ! 01-Apr-2005
645 !       (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.16,  0.84  /) ! "old"
647         real, save :: fr8b_fine(8) =   &
648         (/ 0.0245, 0.1472, 0.3501, 0.3317, 0.1251, 0.0186, 0.0011, 0.0/) !jdf Testbed
649 !jdf    (/ 0.006, 0.024, 0.231, 0.341, 0.140, 0.258, 0., 0./) !5-Apr-2005 values
650 !       (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
651 !       (/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
652 !       (/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.00, 0.00 /)
654         real :: qfine_so4, qfine_no3, qfine_cl, qfine_msa,   &
655                 qfine_co3, qfine_nh4, qfine_na, qfine_ca, qfine_oin,   &
656                 qfine_oc, qfine_bc, qfine_hysw, qfine_watr, qfine_vol
657         real :: qcoar_so4, qcoar_no3, qcoar_cl, qcoar_msa,   &
658                 qcoar_co3, qcoar_nh4, qcoar_na, qcoar_ca, qcoar_oin,   &
659                 qcoar_oc, qcoar_bc, qcoar_hysw, qcoar_watr, qcoar_vol
661 !wig 10-May-2004, added z
662         real, dimension( ims:ime, kms:kme, jms:jme ) :: z
664         character*80 msg
665         character*10 dumname
668 !   *** currently only works for ntype_aer = 1
669         itype = 1
670         iphase = ai_phase
671         m = 1
673 !wig 10-May-2004, block begin...
674 ! calculate the mid-level heights
675     do jt = jts, min(jte,jde-1)
676        do kt = kts, kte-1
677           do it = its, min(ite,ide-1)
678              z(it,kt,jt) = (z_at_w(it,kt,jt)+z_at_w(it,kt+1,jt))*0.5
679           end do
680        end do
681     end do
682 !wig 10-May-2004, ...end block
684 !   set the partitioning fractions for either 8 or 4 bins
685         if (nsize_aer(itype) == 8) then
686             fr_coar(:) = fr8b_coar(:)
687             fr_fine(:) = fr8b_fine(:)
688         else if (nsize_aer(itype) == 4) then
689             do m = 1, nsize_aer(itype)
690                 fr_coar(m) = fr8b_coar(2*m) + fr8b_coar(2*m-1)
691                 fr_fine(m) = fr8b_fine(2*m) + fr8b_fine(2*m-1)
692             end do
693         else
694             write(msg,'(a,i5)')   &
695                 'subr mosaic_init_wrf_mixrats_opt2' //   &
696                 ' - nsize_aer(itype) must be 4 or 8 but = ', nsize_aer(itype)
697             call peg_error_fatal( lunout, msg )
698         end if
701 !   compute initial values (currently uniform in x, y, AND z)
702 !   load them into the chem array
703 !   also load them into the rclm array for diagnostic output
705         rclm(:,:) = 0.0
706         vtot_ofmode(:) = 0.0
708 !   Set values for fine and coarse mass concentrations, and
709 !   compute fine/coarse volume concentrations. These are also set in
710 !   bdy_chem_value_mosaic, below.
711 !   (these are latest values from 1-Apr-2005 discussions)
713 !   rce 10-may-2005 - changed qfine_cl, _na, _oin to values determined by jdf
714         qfine_so4 = 2.14
715         qcoar_so4 = 0.242
716         qfine_no3 = 0.11
717         qcoar_no3 = 0.03
718 !       qfine_cl  = 0.3
719         qfine_cl  = 0.14     ! 10-May-2005
720         qcoar_cl  = 0.139
721         qfine_msa = 0.0
722         qcoar_msa = 0.0
723         qfine_co3 = 0.0
724         qcoar_co3 = 0.0
725         qfine_nh4 = 0.83
726         qcoar_nh4 = 0.10
727 !       qfine_na  = 0.2
728         qfine_na  = 0.1      ! 10-May-2005
729         qcoar_na  = 0.09
730         qfine_ca  = 0.0
731         qcoar_ca  = 0.0
732 !       qfine_oin = 6.2
733         qfine_oin = 3.48     ! 10-May-2005
734         qcoar_oin = 0.35
735         qfine_oc  = 1.00
736         qcoar_oc  = 1.50
737         qfine_bc  = 0.2
738         qcoar_bc  = 0.075
739         qfine_hysw = 0.0
740         qcoar_hysw = 0.0
741         qfine_watr = 0.0
742         qcoar_watr = 0.0
743 #if (CASENAME == 2)
744 !       qfine_so4 = 0.50
745 !       qcoar_so4 = 0.50/2.0
746 !       qfine_no3 = 0.75
747 !       qcoar_no3 = 0.75/2.0
748 !       qfine_cl  = 0.10
749 !       qcoar_cl  = 0.10/2.0
750 !       qfine_nh4 = 0.30
751 !       qcoar_nh4 = 0.30/2.0
752 !       qfine_na  = 0.20
753 !       qcoar_na  = 0.20/2.0
754 !       qfine_oin = 4.00
755 !       qcoar_oin = 4.00/2.0
756 !       qfine_oc  = 3.00
757 !       qcoar_oc  = 3.00/2.0
758 !       qfine_bc  = 0.50
759 !       qcoar_bc  = 0.50/2.0
760         qfine_so4 = 0.300
761         qcoar_so4 = 0.300/2.0
762         qfine_no3 = 0.001
763         qcoar_no3 = 0.001/2.0
764         qfine_cl  = 1.050
765         qcoar_cl  = 1.05/02.0
766         qfine_nh4 = 0.094
767         qcoar_nh4 = 0.094/2.0
768         qfine_na  = 0.700
769         qcoar_na  = 0.700/2.0
770         qfine_oin = 4.500
771         qcoar_oin = 4.500/2.0
772         qfine_oc  = 0.088
773         qcoar_oc  = 0.088/2.0
774         qfine_bc  = 0.013
775         qcoar_bc  = 0.013/2.0
776 #endif
777 #if (CASENAME == 3)
778         qfine_so4 = 0.300
779         qcoar_so4 = 0.0
780         qfine_nh4 = 0.094
781         qcoar_nh4 = 0.0
782         qfine_no3 = 0.001
783         qcoar_no3 = 0.0
784         qfine_cl  = 0.0
785         qcoar_cl  = 1.050
786         qfine_na  = 0.000
787         qcoar_na  = 0.700
788         qfine_oin = 4.500
789         qcoar_oin = 4.500/2.0
790         qfine_oc  = 0.088
791         qcoar_oc  = 0.0
792         qfine_bc  = 0.013
793         qcoar_bc  = 0.0
794 #endif
796 !!$! old values from 23-Apr-2004:
797 !!$     qfine_so4 = 2.554
798 !!$     qcoar_so4 = 0.242
799 !!$     qfine_no3 = 0.07
800 !!$     qcoar_no3 = 0.03
801 !!$     qfine_cl  = 0.324
802 !!$     qcoar_cl  = 0.139
803 !!$     qfine_msa = 0.0
804 !!$     qcoar_msa = 0.0
805 !!$     qfine_co3 = 0.0
806 !!$     qcoar_co3 = 0.0
807 !!$     qfine_nh4 = 0.98
808 !!$     qcoar_nh4 = 0.10
809 !!$     qfine_na  = 0.21
810 !!$     qcoar_na  = 0.09
811 !!$     qfine_ca  = 0.0
812 !!$     qcoar_ca  = 0.0
813 !!$     qfine_oin = 0.15
814 !!$     qcoar_oin = 0.35
815 !!$     qfine_oc  = 1.00
816 !!$     qcoar_oc  = 1.50
817 !!$     qfine_bc  = 0.175
818 !!$     qcoar_bc  = 0.075
819 !!$     qfine_hysw = 0.0
820 !!$     qcoar_hysw = 0.0
821 !!$     qfine_watr = 0.0
822 !!$     qcoar_watr = 0.0
825 ! qfine_so4 ... are ug/m3 so 1.0e-6 factor gives g/m3
826 ! dens_so4  ... are g/cm3;  qfine_vol ... are cm3/m3
827         qfine_vol = 1.0e-6 * (   &
828             (qfine_so4/dens_so4_aer) + (qfine_no3/dens_no3_aer) +   &
829             (qfine_cl /dens_cl_aer ) + (qfine_msa/dens_msa_aer) +   &
830             (qfine_co3/dens_co3_aer) + (qfine_nh4/dens_nh4_aer) +   &
831             (qfine_na /dens_na_aer ) + (qfine_ca /dens_ca_aer ) +   &
832             (qfine_oin/dens_oin_aer) + (qfine_oc /dens_oc_aer ) +   &
833             (qfine_bc /dens_bc_aer ) )
834         qcoar_vol =  1.0e-6 * (  &
835             (qcoar_so4/dens_so4_aer) + (qcoar_no3/dens_no3_aer) +   &
836             (qcoar_cl /dens_cl_aer ) + (qcoar_msa/dens_msa_aer) +   &
837             (qcoar_co3/dens_co3_aer) + (qcoar_nh4/dens_nh4_aer) +   &
838             (qcoar_na /dens_na_aer ) + (qcoar_ca /dens_ca_aer ) +   &
839             (qcoar_oin/dens_oin_aer) + (qcoar_oc /dens_oc_aer ) +   &
840             (qcoar_bc /dens_bc_aer ) )
842         do 2900 m = 1, nsize_aer(itype)
843         do 2800 l = 1, 15
845             if (l .eq. 1) then
846                 l1 = lptr_so4_aer(m,itype,iphase)
847                 qfine = qfine_so4
848                 qcoar = qcoar_so4
849             else if (l .eq. 2) then
850                 l1 = lptr_no3_aer(m,itype,iphase)
851                 qfine = qfine_no3
852                 qcoar = qcoar_no3
853             else if (l .eq. 3) then
854                 l1 = lptr_cl_aer(m,itype,iphase)
855                 qfine = qfine_cl
856                 qcoar = qcoar_cl
857             else if (l .eq. 4) then
858                 l1 = lptr_msa_aer(m,itype,iphase)
859                 qfine = qfine_msa
860                 qcoar = qcoar_msa
861             else if (l .eq. 5) then
862                 l1 = lptr_co3_aer(m,itype,iphase)
863                 qfine = qfine_co3
864                 qcoar = qcoar_co3
865             else if (l .eq. 6) then
866                 l1 = lptr_nh4_aer(m,itype,iphase)
867                 qfine = qfine_nh4
868                 qcoar = qcoar_nh4
869             else if (l .eq. 7) then
870                 l1 = lptr_na_aer(m,itype,iphase)
871                 qfine = qfine_na
872                 qcoar = qcoar_na
873             else if (l .eq. 8) then
874                 l1 = lptr_ca_aer(m,itype,iphase)
875                 qfine = qfine_ca
876                 qcoar = qcoar_ca
877             else if (l .eq. 9) then
878                 l1 = lptr_oin_aer(m,itype,iphase)
879                 qfine = qfine_oin
880                 qcoar = qcoar_oin
881             else if (l .eq. 10) then
882                 l1 = lptr_oc_aer(m,itype,iphase)
883                 qfine = qfine_oc
884                 qcoar = qcoar_oc
885             else if (l .eq. 11) then
886                 l1 = lptr_bc_aer(m,itype,iphase)
887                 qfine = qfine_bc
888                 qcoar = qcoar_bc
889             else if (l .eq. 12) then
890                 l1 = hyswptr_aer(m,itype)
891                 qfine = qfine_hysw
892                 qcoar = qcoar_hysw
893             else if (l .eq. 13) then
894                 l1 = waterptr_aer(m,itype)
895                 qfine = qfine_watr
896                 qcoar = qcoar_watr
897             else if (l .eq. 14) then
898                 l1 = numptr_aer(m,itype,iphase)
899                 dumvol1p = sqrt(volumlo_sect(m,itype)*volumhi_sect(m,itype))
900                 qfine = qfine_vol/dumvol1p
901                 qcoar = qcoar_vol/dumvol1p
902                 vtot_ofmode(m) =   &
903                         qfine_vol*fr_fine(m) + qcoar_vol*fr_coar(m)
904             else
905                 goto 2800
906             end if
907             l3 = l1
909             if ((l1 .gt. 0) .and. (l1 .le. lmaxd) .and.   &
910                 (l3 .ge. param_first_scalar)) then
911                 qval = qfine*fr_fine(m) + qcoar*fr_coar(m)
912                 rclm(1,l1) = qval
914                 do jt = jts, min(jte,jde-1)
915                 do kt = kts, kte-1
916                 do it = its, min(ite,ide-1)
918 !wig 28-May-2004, begin block...
919 ! Determine height multiplier...
920 ! This should mimic the calculation in sorgam_set_aer_bc_pnnl,
921 ! sorgam_init_aer_ic_pnnl, bdy_chem_value_mosaic
922 !!$!    Height(m)     Multiplier
923 !!$!    ---------     ----------
924 !!$!    <=2000        1.0
925 !!$!    2000<z<3000   linear transition zone to 0.5
926 !!$!    3000<z<5000   linear transition zone to 0.25
927 !!$!    >=5000        0.25
928 !!$!
929 !!$! which translates to:
930 !!$!    2000<z<3000   mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
931 !!$!    3000<z<5000   mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
932 !!$! or in reduced form:
933 !!$                if( z(it,kt,jt) <= 2000. ) then
934 !!$                   mult = 1.0
935 !!$                elseif( z(it,kt,jt) > 2000. &
936 !!$                        .and. z(it,kt,jt) <= 3000. ) then
937 !!$                   mult = 1.0 - 0.0005*(z(it,kt,jt)-2000.)
938 !!$                elseif( z(it,kt,jt) > 3000. &
939 !!$                        .and. z(it,kt,jt) <= 5000. ) then
940 !!$                   mult = 0.5 - 1.25e-4*(z(it,kt,jt)-3000.)
941 !!$                else
942 !!$                   mult = 0.25
943 !!$           end if
944 ! Updated 1-Apr-2005:
945 #if (CASENAME == 1)
946 ! further updated 7-20-05 egc:   species specific profiles based on MIRAG2 output
947          mult = 1.0
948          if ( (l3 == 1) .or. (l3 == 2) .or. (l3 == 6) ) then
949 ! so4, no3 and nh4 aerosol profiles
950              if ( z(it,kt,jt) <= 1000. ) then
951                       mult = 1.0
952                   elseif( z(it,kt,jt) > 1000. &
953                         .and. z(it,kt,jt) <= 4000. ) then
954                       mult = 1.0 - 2.367e-4*(z(it,kt,jt)-1000.)
955                    elseif( z(it,kt,jt) > 4000. &
956                         .and. z(it,kt,jt) <= 5500. ) then
957                       mult = 0.29 - 4.0e-5*(z(it,kt,jt)-4000.)
958                    else
959                       mult = 0.23
960                end if
961            else if ( ( l3 == 3) .or. (l3 ==7) ) then
962 ! na and cl aerosol profiles
963              if ( z(it,kt,jt) <= 100. ) then
964                        mult = 1.0
965                   elseif( z(it,kt,jt) > 100. &
966                         .and. z(it,kt,jt) <= 265. ) then
967                       mult = 1.0 - 2.9e-3*(z(it,kt,jt)-100.)
968                    elseif( z(it,kt,jt) > 265. &
969                         .and. z(it,kt,jt) <= 2000. ) then
970                       mult = 0.52 - 2.94e-4*(z(it,kt,jt)-265.)
971                    elseif( z(it,kt,jt) > 2000. &
972                         .and. z(it,kt,jt) <= 7000. ) then
973                        mult = 0.01
974                    else
975                       mult = 1.e-10
976                end if
977           else if  ( l3 == 10)  then
978 ! oc aerosol profile
979              if ( z(it,kt,jt) <= 600. ) then
980                        mult = 1.0
981                   elseif( z(it,kt,jt) > 600. &
982                         .and. z(it,kt,jt) <= 3389. ) then
983                       mult = 1.0 - 2.04e-4*(z(it,kt,jt)-600.)
984                    elseif( z(it,kt,jt) > 3389. &
985                         .and. z(it,kt,jt) <= 8840. ) then
986                       mult = 0.43 - 7.045e-5*(z(it,kt,jt)-3389.)
987                    else
988                       mult = 0.046
989                end if
990           else if  ( l3 == 11)  then
991 ! bc aerosol profile
992              if ( z(it,kt,jt) <= 100. ) then
993                        mult = 1.0
994                   elseif( z(it,kt,jt) > 100. &
995                         .and. z(it,kt,jt) <= 3400. ) then
996                       mult = 1.0 - 2.51e-4*(z(it,kt,jt)-100.)
997                    elseif( z(it,kt,jt) > 3400. &
998                         .and. z(it,kt,jt) <= 8400. ) then
999                       mult = 0.172 -2.64e-5*(z(it,kt,jt)-3400.)
1000                    else
1001                       mult = 0.04
1002                end if
1003          else if  ( l3 == 9)  then
1004 ! oin aerosol profile
1005              if ( z(it,kt,jt) <= 1580. ) then
1006                        mult = 1.0 + 1.77e-4 *z(it,kt,jt)
1007                   elseif( z(it,kt,jt) > 1580. &
1008                         .and. z(it,kt,jt) <= 5280. ) then
1009                       mult = 1.28 - 2.65e-4*(z(it,kt,jt)-1580.)
1010                    else
1011                       mult = 0.30
1012                end if
1013         else
1014 ! generic profile for other other species (which should have groundlevel values=0)
1015 #endif
1016 !    Height(m)     Multiplier
1017 !    ---------     ----------
1018 !    <=2000        1.0
1019 !    2000<z<3000   linear transition zone to 0.25
1020 !    3000<z<5000   linear transision zone to 0.125
1021 !    >=5000        0.125
1023 ! which translates to:
1024 !    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
1025 !    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
1026 ! or in reduced form:
1027                    if( z(it,kt,jt) <= 2000. ) then
1028                       mult = 1.0
1029                    elseif( z(it,kt,jt) > 2000. &
1030                         .and. z(it,kt,jt) <= 3000. ) then
1031                       mult = 1.0 - 0.00075*(z(it,kt,jt)-2000.)
1032                    elseif( z(it,kt,jt) > 3000. &
1033                         .and. z(it,kt,jt) <= 5000. ) then
1034                       mult = 0.25 - 4.166666667e-5*(z(it,kt,jt)-3000.)
1035                    else
1036                       mult = 0.125
1037            end if
1038 #if (CASENAME == 1)
1039         end if                          !close l3 comparison block
1040 #endif
1041 #if (CASENAME == 2)
1042 !     if( z(it,kt,jt) <= 3000. ) then
1043 !       mult = 1.0
1044 !     elseif( z(it,kt,jt) > 3000. &
1045 !       .and. z(it,kt,jt) <= 5000. ) then
1046 !       mult = 1.0 - 0.0004375*(z(it,kt,jt)-3000.)
1047 !     else
1048 !        mult = 0.125
1049 !     end if
1051 !     if( z(it,kt,jt) <= 5000. ) then
1052 !       mult = 1.0
1053 !     elseif( z(it,kt,jt) > 5000. &
1054 !       .and. z(it,kt,jt) <= 7000. ) then
1055 !       mult = 1.0 - 0.0004375*(z(it,kt,jt)-5000.)
1056 !     else
1057 !        mult = 0.125
1058 !     end if
1059       mult=1.0
1060       if ( (l3 == lptr_so4_aer(m,itype,iphase)) .or. &
1061            (l3 == lptr_no3_aer(m,itype,iphase)) .or. &
1062            (l3 == lptr_nh4_aer(m,itype,iphase)) ) then
1063 ! so4, no3 and nh4 aerosol profiles
1064         if( z(it,kt,jt) <= 500. ) then
1065           mult = 1.000
1066         elseif( z(it,kt,jt) > 500. &
1067           .and. z(it,kt,jt) <= 1000. ) then
1068           mult = 1.000 - 0.0010740*(z(it,kt,jt)-500.)
1069         elseif( z(it,kt,jt) > 1000. &
1070           .and. z(it,kt,jt) <= 5000. ) then
1071           mult = 0.463 - 0.0001110*(z(it,kt,jt)-1000.)
1072         elseif( z(it,kt,jt) > 5000. &
1073           .and. z(it,kt,jt) <= 20000. ) then
1074           mult = 0.019
1075         else
1076           mult = 0.019
1077         end if
1078       else if ( (l3 == lptr_na_aer(m,itype,iphase)) .or. &
1079                 (l3 == lptr_cl_aer(m,itype,iphase)) ) then
1080 ! na and cl aerosol profiles
1081         if( z(it,kt,jt) <= 500. ) then
1082           mult = 1.000
1083         elseif( z(it,kt,jt) > 500. &
1084           .and. z(it,kt,jt) <= 1000. ) then
1085           mult = 1.000 - 0.0014260*(z(it,kt,jt)-500.)
1086         elseif( z(it,kt,jt) > 1000. &
1087           .and. z(it,kt,jt) <= 5000. ) then
1088           mult = 0.287 - 0.0000643*(z(it,kt,jt)-1000.)
1089         elseif( z(it,kt,jt) > 5000. &
1090           .and. z(it,kt,jt) <= 10000. ) then
1091           mult = 0.030
1092         elseif( z(it,kt,jt) > 10000. &
1093           .and. z(it,kt,jt) <= 20000. ) then
1094           mult = 0.030 - 0.0000028*(z(it,kt,jt)-10000.)
1095         else
1096           mult = 0.002
1097         end if
1098       else if ( (l3 == lptr_oc_aer(m,itype,iphase)) .or. &
1099                 (l3 == lptr_bc_aer(m,itype,iphase)) ) then
1100 ! oc and bc aerosol profile
1101         if( z(it,kt,jt) <= 500. ) then
1102           mult = 1.000
1103         elseif( z(it,kt,jt) > 500. &
1104           .and. z(it,kt,jt) <= 1000. ) then
1105           mult = 1.000 - 0.0002500*(z(it,kt,jt)-500.)
1106         elseif( z(it,kt,jt) > 1000. &
1107           .and. z(it,kt,jt) <= 5000. ) then
1108           mult = 0.875 - 0.0001843*(z(it,kt,jt)-1000.)
1109         elseif( z(it,kt,jt) > 5000. &
1110           .and. z(it,kt,jt) <= 20000. ) then
1111           mult = 0.138
1112         else
1113           mult = 0.138
1114         end if
1115       else if ( l3 == lptr_oin_aer(m,itype,iphase) ) then
1116 ! oin aerosol profile
1117         if( z(it,kt,jt) <= 500. ) then
1118           mult = 1.000
1119         elseif( z(it,kt,jt) > 500. &
1120           .and. z(it,kt,jt) <= 1000. ) then
1121           mult = 1.000 - 0.0007340*(z(it,kt,jt)-500.)
1122         elseif( z(it,kt,jt) > 1000. &
1123           .and. z(it,kt,jt) <= 5000. ) then
1124           mult = 0.633 - 0.0001103*(z(it,kt,jt)-1000.)
1125         elseif( z(it,kt,jt) > 5000. &
1126           .and. z(it,kt,jt) <= 10000. ) then
1127           mult = 0.192 - 0.0000080*(z(it,kt,jt)-5000.)
1128         elseif( z(it,kt,jt) > 10000. &
1129           .and. z(it,kt,jt) <= 20000. ) then
1130           mult = 0.152 - 0.0000137*(z(it,kt,jt)-10000.)
1131         else
1132           mult = 0.015
1133         end if
1134       else
1135         if( z(it,kt,jt) <= 5000. ) then
1136           mult = 1.0
1137         elseif( z(it,kt,jt) > 5000. &
1138           .and. z(it,kt,jt) <= 7000. ) then
1139           mult = 1.0 - 0.0004375*(z(it,kt,jt)-5000.)
1140         else
1141            mult = 0.125
1142         end if
1143       end if
1144 #endif
1145 #if (CASENAME == 3)
1146       mult = 1.000
1147       if( z(it,kt,jt) <= 500. ) then
1148         mult = 1.000
1149       elseif( z(it,kt,jt) > 500. &
1150         .and. z(it,kt,jt) <= 1000. ) then
1151         mult = 1.000 - 0.0010740*(z(it,kt,jt)-500.)
1152       elseif( z(it,kt,jt) > 1000. &
1153         .and. z(it,kt,jt) <= 5000. ) then
1154         mult = 0.463 - 0.0001110*(z(it,kt,jt)-1000.)
1155       else
1156         mult = 0.019
1157       end if
1158 #endif
1160                     chem(it,kt,jt,l3) = mult*rclm(1,l1)
1161 !wig 28-May-2004, ...end block
1162 !                   chem(it,kt,jt,l3) = rclm(1,l1)
1163                 end do
1164                 end do
1165                 end do
1166             end if
1169 2800   continue
1170 2900   continue
1173 !   do diagnostic output
1175         dumarr(:) = 0.0
1176         msg = ' '
1177         call peg_message( lunout, msg )
1178         msg = '*** subr mosaic_init_wrf_mixrats_opt2 results'
1179         call peg_message( lunout, msg )
1180         msg = '    mass in ug/m3     number in #/m3     volume in cm3/m3'
1181         call peg_message( lunout, msg )
1182         msg = ' '
1183         call peg_message( lunout, msg )
1184         msg = ' mode  l  l1  species      conc'
1185         call peg_message( lunout, msg )
1187         do 3190 mdum = 1, nsize_aer(itype)+1
1188         m = min( mdum, nsize_aer(itype) )
1189         msg = ' '
1190         call peg_message( lunout, msg )
1191         do 3150 l = 1, ncomp_plustracer_aer(itype)+4
1193             if (l .le. ncomp_plustracer_aer(itype)) then
1194                 l1 = massptr_aer(l,m,itype,iphase)
1195                 dumname = name_aer(l,itype)
1196                 dum = rclm(1,l1)
1197             else if (l .eq. ncomp_plustracer_aer(itype)+1) then
1198                 l1 = hyswptr_aer(m,itype)
1199                 dumname = 'hystwatr'
1200                 dum = rclm(1,l1)
1201             else if (l .eq. ncomp_plustracer_aer(itype)+2) then
1202                 l1 = waterptr_aer(m,itype)
1203                 dumname = 'water'
1204                 dum = rclm(1,l1)
1205             else if (l .eq. ncomp_plustracer_aer(itype)+3) then
1206                 l1 = numptr_aer(m,itype,iphase)
1207                 dumname = 'number'
1208                 dum = rclm(1,l1)
1209             else if (l .eq. ncomp_plustracer_aer(itype)+4) then
1210                 l1 = 0
1211                 dumname = 'volume'
1212                 dum = vtot_ofmode(m)
1213             else
1214                 dumname = '=BADBAD='
1215                 l1 = -1
1216                 dum = -1.0
1217             end if
1219             if (mdum .le. nsize_aer(itype)) then
1220                 dumarr(l) = dumarr(l) + dum
1221                 write(msg,9620) m, l, l1, dumname, dum
1222             else
1223                 write(msg,9625) l, dumname, dumarr(l)
1224             end if
1225             call peg_message( lunout, msg )
1227 3150   continue
1228 3190   continue
1230 9620    format( 3i4, 2x, a, 3(1pe12.3) )
1231 9625    format( ' sum', i4, '   -', 2x, a, 3(1pe12.3) )
1234 !   all done
1235         return
1236         end subroutine mosaic_init_wrf_mixrats_opt2
1239 !-----------------------------------------------------------------------
1240         real function erfc_num_recipes( x )
1242 !   from press et al, numerical recipes, 1990, page 164
1244         implicit none
1245         real x
1246         double precision erfc_dbl, dum, t, zz
1248         zz = abs(x)
1249         t = 1.0/(1.0 + 0.5*zz)
1251 !       erfc_num_recipes =
1252 !     &   t*exp( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
1253 !     &   t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
1254 !     &                                    t*(-1.13520398 +
1255 !     &   t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
1257         dum =  ( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +   &
1258           t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +   &
1259                                            t*(-1.13520398 +   &
1260           t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
1262         erfc_dbl = t * exp(dum)
1263         if (x .lt. 0.0) erfc_dbl = 2.0d0 - erfc_dbl
1265         erfc_num_recipes = erfc_dbl
1267         return
1268         end function erfc_num_recipes     
1271 !-----------------------------------------------------------------------
1272         end module module_mosaic_initmixrats
1277 !-----------------------------------------------------------------------
1278         subroutine bdy_chem_value_mosaic ( chem_bv, alt, zz, nch, config_flags )
1280 ! provides boundary values for the mosaic aerosol species
1282 ! it is outside of the module declaration because of potential
1283 !    module1 --> module2 --> module1 use conflicts
1285 ! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
1286 !     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now
1288         use module_configure, only:  grid_config_rec_type
1289         use module_state_description, only:  param_first_scalar,   &
1290                         aer_bc_pnnl
1291         use module_data_mosaic_asect
1292         use module_data_mosaic_other
1293         implicit none
1295 ! arguments
1296         REAL,    intent(OUT)  :: chem_bv    ! boundary value for chem(-,-,-,nch)
1297         REAL,    intent(IN)   :: alt        ! inverse density
1298         REAL,    intent(IN)   :: zz         ! height
1299         INTEGER, intent(IN)   :: nch        ! index number of chemical species
1300         TYPE (grid_config_rec_type), intent(in) :: config_flags
1302 ! local variables
1303         integer :: iphase, itype, m
1304     logical :: foundit
1306         real, parameter :: chem_bv_def = 1.0e-20
1307 !wig 10-May-2004, added mult
1308         real :: dumvol1p, mult
1309         integer :: iflg
1310         real :: qcoar, qfine, qval
1312         real :: fr_coar(8), fr_fine(8)
1314 !wig 1-Apr-2005,  Updated fractional breakdown between bins. (See also
1315 !                 mosaic_init_wrf_mixrats_opt2, above, and mosaic_addemiss
1316 !                 in module_mosaic_addemiss.F). Note that these values no
1317 !                 longer match those in mosaic_addemiss.
1318         real, save :: fr8b_coar(8) =   &
1319         (/ 0.0,  0.0,  0.0,  0.019,  0.0212,  0.1101,  0.2751, 0.7018 /) ! jdf Testbed
1320 !jdf    (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.300, 0.700 /) ! 10-May-2005
1321 !       (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.933, 0.067 /) ! 01-Apr-2005
1322 !               (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.16, 0.84 /)
1323         real, save :: fr8b_fine(8) =   &
1324         (/ 0.0245, 0.1472, 0.3501, 0.3317, 0.1251, 0.0186, 0.0011, 0.0/) !jdf Testbed
1325 !jdf    (/ 0.006, 0.024, 0.231, 0.341, 0.140, 0.258, 0., 0./) !5-Apr-2005 values
1326 !       (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
1327 !               (/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
1328 !               (/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.00, 0.00 /)
1330         real :: qfine_so4, qfine_no3, qfine_cl, qfine_msa,   &
1331                 qfine_co3, qfine_nh4, qfine_na, qfine_ca, qfine_oin,   &
1332                 qfine_oc, qfine_bc, qfine_hysw, qfine_watr, qfine_vol
1333         real :: qcoar_so4, qcoar_no3, qcoar_cl, qcoar_msa,   &
1334                 qcoar_co3, qcoar_nh4, qcoar_na, qcoar_ca, qcoar_oin,   &
1335                 qcoar_oc, qcoar_bc, qcoar_hysw, qcoar_watr, qcoar_vol
1337         character*80 msg
1341 ! when aer_bc_opt /= aer_bc_pnnl,
1342 ! just chem_bv=near-zero for all species
1343         chem_bv = chem_bv_def
1344         if (config_flags%aer_bc_opt /= aer_bc_pnnl) return
1345         if (nch < param_first_scalar) return
1348 !   *** currently only works for ntype_aer = 1
1349         itype = 1
1350         iphase = ai_phase
1351         m = 1 !This is here for size, but is overridden by loop below.
1355 !   following is for aer_bc_opt == aer_bc_pnnl
1356 !       when boundary values are set for Texas August 2000 simulations
1358 !   set the partitioning fractions for either 8 or 4 bins
1359         if (nsize_aer(itype) == 8) then
1360             fr_coar(:) = fr8b_coar(:)
1361             fr_fine(:) = fr8b_fine(:)
1362         else if (nsize_aer(itype) == 4) then
1363             do m = 1, nsize_aer(itype)
1364                 fr_coar(m) = fr8b_coar(2*m) + fr8b_coar(2*m-1)
1365                 fr_fine(m) = fr8b_fine(2*m) + fr8b_fine(2*m-1)
1366             end do
1367         else
1368             write(msg,'(a,i5)')   &
1369                 'subr bdy_chem_value_mosaic' //   &
1370                 ' - nsize_aer(itype) must be 4 or 8 but = ', nsize_aer(itype)
1371             call wrf_error_fatal( msg )
1372         end if
1374 ! Determine height multiplier...
1375 ! This should mimic the calculation in sorgam_set_aer_bc_pnnl,
1376 ! sorgam_init_aer_ic_pnnl, and mosaic_init_wrf_mixrats_opt2
1377 !!$!    Height(m)     Multiplier
1378 !!$!    ---------     ----------
1379 !!$!    <=2000        1.0
1380 !!$!    2000<z<3000   linear transition zone to 0.5
1381 !!$!    3000<z<5000   linear transision zone to 0.25
1382 !!$!    >=5000        0.25
1383 !!$!
1384 !!$! which translates to:
1385 !!$!    2000<z<3000   mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
1386 !!$!    3000<z<5000   mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
1387 !!$! or in reduced form:
1388 !!$        if( zz <= 2000. ) then
1389 !!$           mult = 1.0
1390 !!$        elseif( zz > 2000. &
1391 !!$             .and. zz <= 3000. ) then
1392 !!$           mult = 1.0 - 0.0005*(zz-2000.)
1393 !!$        elseif( zz > 3000. &
1394 !!$             .and. zz <= 5000. ) then
1395 !!$           mult = 0.5 - 1.25e-4*(zz-3000.)
1396 !!$        else
1397 !!$           mult = 0.25
1398 !!$        end if
1399 #if (CASENAME == 1)
1400     mult = 1.0
1401     SIZE_LOOP : do m=1,nsize_aer(itype)
1402        if( ( nch .eq. lptr_so4_aer(m,itype,iphase) ) .or.  &
1403             (nch .eq. lptr_no3_aer(m,itype,iphase) ) .or.  &
1404             (nch .eq. lptr_nh4_aer(m,itype,iphase) )  )then
1405 ! so4, no3 and nh4 aerosol profiles
1406           if ( zz <= 1000. ) then
1407              mult = 1.0
1408                   elseif( zz > 1000. &
1409                  .and. zz <= 4000. ) then
1410              mult = 1.0 - 2.367e-4*(zz-1000.)
1411           elseif( zz > 4000. &
1412                  .and. zz <= 5500. ) then
1413              mult = 0.29 - 4.0e-5*(zz-4000.)
1414           else
1415              mult = 0.23
1416           end if
1417           exit SIZE_LOOP
1418        else if ( (nch .eq. lptr_na_aer(m,itype,iphase) ) .or.   &
1419                      (nch .eq. lptr_cl_aer(m,itype,iphase) ) ) then
1420 ! na and cl aerosol profiles
1421           if ( zz <= 100. ) then
1422              mult = 1.0
1423                   elseif( zz > 100. &
1424                   .and. zz <= 265. ) then
1425              mult = 1.0 - 2.9e-3*(zz-100.)
1426           elseif( zz > 265. &
1427                   .and. zz <= 2000. ) then
1428              mult = 0.52 - 2.94e-4*(zz-265.)
1429           elseif( zz > 2000. &
1430                   .and. zz <= 7000. ) then
1431              mult = 0.01
1432           else
1433              mult = 1.e-10
1434           end if
1435           exit SIZE_LOOP
1436        else if (nch .eq. lptr_oc_aer(m,itype,iphase) )  then
1437 ! oc aerosol profile
1438           if ( zz <= 600. ) then
1439              mult = 1.0
1440                   elseif( zz > 600. &
1441                   .and. zz <= 3389. ) then
1442              mult = 1.0 - 2.04e-4*(zz-600.)
1443           elseif( zz > 3389. &
1444                   .and. zz <= 8840. ) then
1445              mult = 0.43 - 7.045e-5*(zz-3389.)
1446           else
1447              mult = 0.046
1448           end if
1449           exit SIZE_LOOP
1450        else if  (nch .eq. lptr_bc_aer(m,itype,iphase) )   then
1451 ! bc aerosol profile
1452           if ( zz <= 100. ) then
1453              mult = 1.0
1454                   elseif( zz > 100. &
1455                   .and. zz <= 3400. ) then
1456              mult = 1.0 - 2.51e-4*(zz-100.)
1457           elseif( zz > 3400. &
1458                   .and. zz <= 8400. ) then
1459              mult = 0.172 -2.64e-5*(zz-3400.)
1460           else
1461              mult = 0.04
1462           end if
1463           exit SIZE_LOOP
1464        else if  (nch .eq. lptr_oin_aer(m,itype,iphase))  then
1465 ! oin aerosol profile
1466           if ( zz <= 1580. ) then
1467              mult = 1.0 + 1.77e-4 *zz
1468                   elseif( zz > 1580. &
1469                   .and. zz <= 5280. ) then
1470              mult = 1.28 - 2.65e-4*(zz-1580.)
1471           else
1472              mult = 0.30
1473           end if
1474           exit SIZE_LOOP
1475        else
1476 ! generic profile
1477 #endif
1478 ! Updated aerosol profile multiplier 1-Apr-2005:
1479 !    Height(m)     Multiplier
1480 !    ---------     ----------
1481 !    <=2000        1.0
1482 !    2000<z<3000   linear transition zone to 0.25
1483 !    3000<z<5000   linear transision zone to 0.125
1484 !    >=5000        0.125
1486 ! which translates to:
1487 !    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
1488 !    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
1489 ! or in reduced form:
1490         if( zz <= 2000. ) then
1491            mult = 1.0
1492         elseif( zz > 2000. &
1493              .and. zz <= 3000. ) then
1494            mult = 1.0 - 0.00075*(zz-2000.)
1495         elseif( zz > 3000. &
1496              .and. zz <= 5000. ) then
1497            mult = 0.25 - 4.166666667e-5*(zz-3000.)
1498         else
1499            mult = 0.125
1500         end if
1501 #if (CASENAME == 1)
1502       end if                    ! end nc block comparison
1503     end do SIZE_LOOP
1504 #endif
1505 #if (CASENAME == 2)
1506 !     if( zz <= 3000. ) then
1507 !       mult = 1.0
1508 !     elseif( zz > 3000. &
1509 !       .and. zz <= 5000. ) then
1510 !       mult = 1.0 - 0.0004375*(zz-3000.)
1511 !     else
1512 !        mult = 0.125
1513 !     end if
1514 !     if( zz <= 5000. ) then
1515 !       mult = 1.0
1516 !     elseif( zz > 5000. &
1517 !       .and. zz <= 7000. ) then
1518 !       mult = 1.0 - 0.0004375*(zz-5000.)
1519 !     else
1520 !        mult = 0.125
1521 !     end if
1522       mult=1.0
1523       iflg=0
1524       do 2901 m = 1, nsize_aer(itype)
1525       if(iflg.eq.1) goto 2902
1526       if( ( nch .eq. lptr_so4_aer(m,itype,iphase) ) .or.  &
1527         (nch .eq. lptr_no3_aer(m,itype,iphase) ) .or.  &
1528         (nch .eq. lptr_nh4_aer(m,itype,iphase) )  )then
1529 ! so4, no3 and nh4 aerosol profiles
1530         if( zz <= 500. ) then
1531           mult = 1.000
1532         elseif( zz > 500. &
1533           .and. zz <= 1000. ) then
1534           mult = 1.000 - 0.0010740*(zz-500.)
1535         elseif( zz > 1000. &
1536           .and. zz <= 5000. ) then
1537           mult = 0.463 - 0.0001110*(zz-1000.)
1538         elseif( zz > 5000. &
1539           .and. zz <= 20000. ) then
1540           mult = 0.019
1541         else
1542           mult = 0.019
1543         end if
1544         iflg=1
1545       else if ( (nch .eq. lptr_na_aer(m,itype,iphase) ) .or.   &
1546                 (nch .eq. lptr_cl_aer(m,itype,iphase) ) ) then
1547 ! na and cl aerosol profiles
1548         if( zz <= 500. ) then
1549           mult = 1.000
1550         elseif( zz > 500. &
1551           .and. zz <= 1000. ) then
1552           mult = 1.000 - 0.0014260*(zz-500.)
1553         elseif( zz > 1000. &
1554           .and. zz <= 5000. ) then
1555           mult = 0.287 - 0.0000643*(zz-1000.)
1556         elseif( zz > 5000. &
1557           .and. zz <= 10000. ) then
1558           mult = 0.030
1559         elseif( zz > 10000. &
1560           .and. zz <= 20000. ) then
1561           mult = 0.030 - 0.0000028*(zz-10000.)
1562         else
1563           mult = 0.002
1564         end if
1565         iflg=1
1566       else if ( (nch .eq. lptr_oc_aer(m,itype,iphase) ) .or.   &
1567                 (nch .eq. lptr_bc_aer(m,itype,iphase) ) ) then
1568 ! oc and bc aerosol profile
1569         if( zz <= 500. ) then
1570           mult = 1.000
1571         elseif( zz > 500. &
1572           .and. zz <= 1000. ) then
1573           mult = 1.000 - 0.0002500*(zz-500.)
1574         elseif( zz > 1000. &
1575           .and. zz <= 5000. ) then
1576           mult = 0.875 - 0.0001843*(zz-1000.)
1577         elseif( zz > 5000. &
1578           .and. zz <= 20000. ) then
1579           mult = 0.138
1580         else
1581           mult = 0.138
1582         end if
1583         iflg=1
1584       else if  (nch .eq. lptr_oin_aer(m,itype,iphase))  then
1585 ! oin aerosol profile
1586         if( zz <= 500. ) then
1587           mult = 1.000
1588         elseif( zz > 500. &
1589           .and. zz <= 1000. ) then
1590           mult = 1.000 - 0.0007340*(zz-500.)
1591         elseif( zz > 1000. &
1592           .and. zz <= 5000. ) then
1593           mult = 0.633 - 0.0001103*(zz-1000.)
1594         elseif( zz > 5000. &
1595           .and. zz <= 10000. ) then
1596           mult = 0.192 - 0.0000080*(zz-5000.)
1597         elseif( zz > 10000. &
1598           .and. zz <= 20000. ) then
1599           mult = 0.152 - 0.0000137*(zz-10000.)
1600         else
1601           mult = 0.015
1602         end if
1603         iflg=1
1604       else
1605         if( zz <= 5000. ) then
1606           mult = 1.0
1607         elseif( zz > 5000. &
1608           .and. zz <= 7000. ) then
1609           mult = 1.0 - 0.0004375*(zz-5000.)
1610         else
1611            mult = 0.125
1612         end if
1613       end if
1614 2901  continue
1615 2902  continue
1616 #endif
1617 #if (CASENAME == 3)
1618       if( zz <= 500. ) then
1619         mult = 1.000
1620       elseif( zz > 500. &
1621         .and. zz <= 1000. ) then
1622         mult = 1.000 - 0.0010740*(zz-500.)
1623       elseif( zz > 1000. &
1624         .and. zz <= 5000. ) then
1625         mult = 0.463 - 0.0001110*(zz-1000.)
1626       else
1627         mult = 0.019
1628       end if
1629 #endif
1630 !   Set values for fine and coarse mass concentrations, and
1631 !   compute fine/coarse volume concentrations. These are also set in
1632 !   mosaic_init_wrf_mixrats_opt2.
1633 !   (these are latest values from 1-Apr-2005 discussions)
1634 !wig 10-May-2004, added height multiplier (mult*) to each species...
1635         qfine_so4 = mult*2.14
1636         qcoar_so4 = mult*0.242
1637         qfine_no3 = mult*0.11
1638         qcoar_no3 = mult*0.03
1639 !       qfine_cl  = mult*0.3
1640         qfine_cl  = mult*0.14     ! 10-May-2005
1641         qcoar_cl  = mult*0.139
1642         qfine_msa = mult*0.0
1643         qcoar_msa = mult*0.0
1644         qfine_co3 = mult*0.0
1645         qcoar_co3 = mult*0.0
1646         qfine_nh4 = mult*0.83
1647         qcoar_nh4 = mult*0.10
1648 !       qfine_na  = mult*0.2
1649         qfine_na  = mult*0.1      ! 10-May-2005
1650         qcoar_na  = mult*0.09
1651         qfine_ca  = mult*0.0
1652         qcoar_ca  = mult*0.0
1653 !       qfine_oin = mult*6.2
1654 !       qfine_oin = 3.48     ! 10-May-2005
1655         qfine_oin = mult*3.48     ! not sure why mult was dropped from qfine_oin
1656         qcoar_oin = mult*0.35
1657         qfine_oc  = mult*1.00
1658         qcoar_oc  = mult*1.50
1659         qfine_bc  = mult*0.2
1660         qcoar_bc  = mult*0.075
1661         qfine_hysw = mult*0.0
1662         qcoar_hysw = mult*0.0
1663         qfine_watr = mult*0.0
1664         qcoar_watr = mult*0.0
1665 #if (CASENAME == 2)
1666 !       qfine_so4 = mult*0.50
1667 !       qcoar_so4 = mult*0.50/2.0
1668 !       qfine_no3 = mult*0.75
1669 !       qcoar_no3 = mult*0.75/2.0
1670 !       qfine_cl  = mult*0.10
1671 !       qcoar_cl  = mult*0.10/2.0
1672 !       qfine_nh4 = mult*0.30
1673 !       qcoar_nh4 = mult*0.30/2.0
1674 !       qfine_na  = mult*0.20
1675 !       qcoar_na  = mult*0.20/2.0
1676 !       qfine_oin = mult*4.00
1677 !       qcoar_oin = mult*4.00/2.0
1678 !       qfine_oc  = mult*3.00
1679 !       qcoar_oc  = mult*3.00/2.0
1680 !       qfine_bc  = mult*0.50
1681 !       qcoar_bc  = mult*0.50/2.0
1682         qfine_so4 = mult*0.300
1683         qcoar_so4 = mult*0.300/2.0
1684         qfine_no3 = mult*0.001
1685         qcoar_no3 = mult*0.001/2.0
1686         qfine_cl  = mult*1.050
1687         qcoar_cl  = mult*1.05/02.0
1688         qfine_nh4 = mult*0.094
1689         qcoar_nh4 = mult*0.094/2.0
1690         qfine_na  = mult*0.700
1691         qcoar_na  = mult*0.700/2.0
1692         qfine_oin = mult*4.500
1693         qcoar_oin = mult*4.500/2.0
1694         qfine_oc  = mult*0.088
1695         qcoar_oc  = mult*0.088/2.0
1696         qfine_bc  = mult*0.013
1697         qcoar_bc  = mult*0.013/2.0
1698 #endif
1699 #if (CASENAME == 3)
1700         qfine_so4 = mult*0.300
1701         qcoar_so4 = mult*0.0
1702         qfine_nh4 = mult*0.094
1703         qcoar_nh4 = mult*0.0
1704         qfine_no3 = mult*0.001
1705         qcoar_no3 = mult*0.0
1706         qfine_cl  = mult*0.0
1707         qcoar_cl  = mult*1.050
1708         qfine_na  = mult*0.000
1709         qcoar_na  = mult*0.700
1710         qfine_oin = mult*4.500
1711         qcoar_oin = mult*4.500/2.0
1712         qfine_oc  = mult*0.088
1713         qcoar_oc  = mult*0.0
1714         qfine_bc  = mult*0.013
1715         qcoar_bc  = mult*0.0
1716 #endif
1717 !!$! old values from 23-Apr-2004:
1718 !!$     qfine_so4 = mult*2.554
1719 !!$     qcoar_so4 = mult*0.242
1720 !!$     qfine_no3 = mult*0.07
1721 !!$     qcoar_no3 = mult*0.03
1722 !!$     qfine_cl  = mult*0.324
1723 !!$     qcoar_cl  = mult*0.139
1724 !!$     qfine_msa = mult*0.0
1725 !!$     qcoar_msa = mult*0.0
1726 !!$     qfine_co3 = mult*0.0
1727 !!$     qcoar_co3 = mult*0.0
1728 !!$     qfine_nh4 = mult*0.98
1729 !!$     qcoar_nh4 = mult*0.10
1730 !!$     qfine_na  = mult*0.21
1731 !!$     qcoar_na  = mult*0.09
1732 !!$     qfine_ca  = mult*0.0
1733 !!$     qcoar_ca  = mult*0.0
1734 !!$     qfine_oin = mult*0.15
1735 !!$     qcoar_oin = mult*0.35
1736 !!$     qfine_oc  = mult*1.00
1737 !!$     qcoar_oc  = mult*1.50
1738 !!$     qfine_bc  = mult*0.175
1739 !!$     qcoar_bc  = mult*0.075
1740 !!$     qfine_hysw = mult*0.0
1741 !!$     qcoar_hysw = mult*0.0
1742 !!$     qfine_watr = mult*0.0
1743 !!$     qcoar_watr = mult*0.0
1746 ! qfine_so4 ... are ug/m3 so 1.0e-6 factor gives g/m3
1747 ! dens_so4  ... are g/cm3;  qfine_vol ... are cm3/m3
1748         qfine_vol = 1.0e-6 * (   &
1749             (qfine_so4/dens_so4_aer) + (qfine_no3/dens_no3_aer) +   &
1750             (qfine_cl /dens_cl_aer ) + (qfine_msa/dens_msa_aer) +   &
1751             (qfine_co3/dens_co3_aer) + (qfine_nh4/dens_nh4_aer) +   &
1752             (qfine_na /dens_na_aer ) + (qfine_ca /dens_ca_aer ) +   &
1753             (qfine_oin/dens_oin_aer) + (qfine_oc /dens_oc_aer ) +   &
1754             (qfine_bc /dens_bc_aer ) )
1755         qcoar_vol = 1.0e-6 * (   &
1756             (qcoar_so4/dens_so4_aer) + (qcoar_no3/dens_no3_aer) +   &
1757             (qcoar_cl /dens_cl_aer ) + (qcoar_msa/dens_msa_aer) +   &
1758             (qcoar_co3/dens_co3_aer) + (qcoar_nh4/dens_nh4_aer) +   &
1759             (qcoar_na /dens_na_aer ) + (qcoar_ca /dens_ca_aer ) +   &
1760             (qcoar_oin/dens_oin_aer) + (qcoar_oc /dens_oc_aer ) +   &
1761             (qcoar_bc /dens_bc_aer ) )
1763         qfine = -1.0e30
1764         qcoar = -1.0e30
1766 !   identify the species by checking "nch" against the "lptr_xxx_a_amode(m)"
1767         do 2900 m = 1, nsize_aer(itype)
1769             if (nch .eq. lptr_so4_aer(m,itype,iphase)) then
1770                 qfine = qfine_so4
1771                 qcoar = qcoar_so4
1772             else if (nch .eq. lptr_no3_aer(m,itype,iphase)) then
1773                 qfine = qfine_no3
1774                 qcoar = qcoar_no3
1775             else if (nch .eq. lptr_cl_aer(m,itype,iphase)) then
1776                 qfine = qfine_cl
1777                 qcoar = qcoar_cl
1778             else if (nch .eq. lptr_msa_aer(m,itype,iphase)) then
1779                 qfine = qfine_msa
1780                 qcoar = qcoar_msa
1781             else if (nch .eq. lptr_co3_aer(m,itype,iphase)) then
1782                 qfine = qfine_co3
1783                 qcoar = qcoar_co3
1784             else if (nch .eq. lptr_nh4_aer(m,itype,iphase)) then
1785                 qfine = qfine_nh4
1786                 qcoar = qcoar_nh4
1787             else if (nch .eq. lptr_na_aer(m,itype,iphase)) then
1788                 qfine = qfine_na
1789                 qcoar = qcoar_na
1790             else if (nch .eq. lptr_ca_aer(m,itype,iphase)) then
1791                 qfine = qfine_ca
1792                 qcoar = qcoar_ca
1793             else if (nch .eq. lptr_oin_aer(m,itype,iphase)) then
1794                 qfine = qfine_oin
1795                 qcoar = qcoar_oin
1796             else if (nch .eq. lptr_oc_aer(m,itype,iphase)) then
1797                 qfine = qfine_oc
1798                 qcoar = qcoar_oc
1799             else if (nch .eq. lptr_bc_aer(m,itype,iphase)) then
1800                 qfine = qfine_bc
1801                 qcoar = qcoar_bc
1802             else if (nch .eq. hyswptr_aer(m,itype)) then
1803                 qfine = qfine_hysw
1804                 qcoar = qcoar_hysw
1805             else if (nch .eq. waterptr_aer(m,itype)) then
1806                 qfine = qfine_watr
1807                 qcoar = qcoar_watr
1808             else if (nch .eq. numptr_aer(m,itype,iphase)) then
1809                 dumvol1p = sqrt(volumlo_sect(m,itype)*volumhi_sect(m,itype))
1810                 qfine = qfine_vol/dumvol1p
1811                 qcoar = qcoar_vol/dumvol1p
1812             end if
1814             if ((qfine >= 0.0) .and. (qcoar >= 0.0)) then
1815                 qval = qfine*fr_fine(m) + qcoar*fr_coar(m)
1816                 chem_bv = qval*alt
1817                 goto 2910
1818             end if
1820 2900   continue
1821 2910   continue
1823         return
1824         end subroutine bdy_chem_value_mosaic