!**********************************************************************************  
! This computer software was prepared by Battelle Memorial Institute, hereinafter
! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of 
! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! MOSAIC module: see module_mosaic_driver.F for references and terms of use
!**********************************************************************************  
!CPP directives to control ic/bc conditions...
!(The directive in module_input_chem_data also needs to be set.)
!  CASENAME = 0   Uses Texas August 2004 case values and profiles
!             1   Uses same concentrations as TX, but uses different
!                 profiles depending on the species. (NEAQS2004 case)
!             2   Mexico City
!             3   Mexico City for Testbed, to be consistent with MADE/SORGAM too
!             4   bc,ic values from wrfinput and wrfbdy files
#define CASENAME 4

	module module_mosaic_initmixrats

! rce 2005-feb-18 - several fixes for indices of dlo/hi_sect, volumlo/hi_sect,
!     which are now (isize,itype)

! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
!     variables in module_data_mosaic_asect

	USE module_state_description

	integer, parameter :: mosaic_init_wrf_mixrats_flagaa = 1
                            ! turns subr mosaic_init_wrf_mixrats on/off

	contains


!-----------------------------------------------------------------------
	subroutine mosaic_init_wrf_mixrats(  &
		iflagaa, config_flags,       &
		chem, alt, z_at_w, g,        &
		ids,ide, jds,jde, kds,kde,   &
		ims,ime, jms,jme, kms,kme,   &
		its,ite, jts,jte, kts,kte    )

!
!   initializes the species and number mixing ratios for each section
!
!   this top level routine simply calls other routines depending on value
!	of config_flags%aer_ic_opt
!
	use module_configure, only:  grid_config_rec_type
	use module_state_description, only:  num_chem, param_first_scalar,   &
			aer_ic_pnnl
	use module_data_mosaic_asect
	use module_data_mosaic_other
	use module_peg_util, only:  peg_message, peg_error_fatal

	implicit none


!   subr arguments
	type(grid_config_rec_type), intent(in) :: config_flags

	integer, intent(in) ::   &
		iflagaa,   &
		ids, ide, jds, jde, kds, kde,   &
		ims, ime, jms, jme, kms, kme,   &
		its, ite, jts, jte, kts, kte

	real, intent(inout),   &
		dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
		chem

	real, intent(in),      &
		dimension( ims:ime, kms:kme, jms:jme ) :: &
        alt, z_at_w
    real, intent(in) :: g

!   local variables
        integer :: i, ic, j, k

#if (CASENAME == 0) || (CASENAME == 1) || (CASENAME == 2) || (CASENAME == 3)
	if (config_flags%aer_ic_opt == aer_ic_pnnl) then
	    call mosaic_init_wrf_mixrats_opt2(   &
		iflagaa, config_flags,           &
		chem, z_at_w, g,                 &
		ids,ide, jds,jde, kds,kde,       &
		ims,ime, jms,jme, kms,kme,       &
		its,ite, jts,jte, kts,kte        )

	else
	    call mosaic_init_wrf_mixrats_opt1(   &
		iflagaa, config_flags,   &
		chem,   &
		ids,ide, jds,jde, kds,kde,   &
		ims,ime, jms,jme, kms,kme,   &
		its,ite, jts,jte, kts,kte    )

	end if

! Aerosol species are returned from above in concentration units. Convert
! them to mixing ratio for use in advection, etc.
    do ic = p_so4_a01,num_chem
       do j = jts,jte
          do k = kts,kte-1
             do i = its,ite
                chem(i,k,j,ic) = chem(i,k,j,ic)*alt(i,k,j)
             end do
          end do
       end do
    end do

! Fill the top z-staggered location to prevent trouble during advection.
    do ic = p_so4_a01,num_chem
       do j = jts,jte
          do i = its,ite
             chem(i,kte,j,ic) = chem(i,kte-1,j,ic)
          end do
       end do
    end do
#endif

	return
	end subroutine mosaic_init_wrf_mixrats


!-----------------------------------------------------------------------
	subroutine mosaic_init_wrf_mixrats_opt1(   &
		iflagaa, config_flags,   &
		chem,   &
		ids,ide, jds,jde, kds,kde,   &
		ims,ime, jms,jme, kms,kme,   &
		its,ite, jts,jte, kts,kte    )

!
!   initializes the species and number mixing ratios for each section
!   based on user-specified lognormal modes that span the size distribution
!
! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
!     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now; 
!     added l1dum
!
	use module_configure, only:  grid_config_rec_type
	use module_state_description, only:  num_chem, param_first_scalar
	use module_data_mosaic_asect
	use module_data_mosaic_other
	use module_peg_util, only:  peg_message, peg_error_fatal

	implicit none

!   subr arguments
	type(grid_config_rec_type), intent(in) :: config_flags

	integer, intent(in) ::   &
		iflagaa,   &
		ids, ide, jds, jde, kds, kde,   &
		ims, ime, jms, jme, kms, kme,   &
		its, ite, jts, jte, kts, kte

	real, intent(inout),   &
		dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
		chem

!   local variables
	integer i, j, k, l, ll, l1, l3, l1dum, m, mdum, nsm
	integer it, jt, kt

	real dum, dumdp, dumrsfc, dumvol,   &
      	  xlo, xhi,   &
	  dumvol1p,  &
      	  pdummb, zdumkm, zscalekm, zfactor

	real vtot_nsm_ofmode(maxd_asize)
	real dumarr(maxd_acomp+5)

	real erfc

!	double precision fracnum, fracvol, tlo, thi
	real fracvol, tlo, thi

	integer nmaxd_nsm
	parameter (nmaxd_nsm = 4)

	integer iphase, itype, ntot_nsm
	integer iiprof_nsm(nmaxd_nsm)
	integer lldum_so4, lldum_nh4, lldum_oc, lldum_bc,   &
		lldum_oin, lldum_na, lldum_cl, lldum_hysw

	real sx_nsm(nmaxd_nsm), sxr2_nsm(nmaxd_nsm),   &
      	  x0_nsm(nmaxd_nsm), x3_nsm(nmaxd_nsm),   &
      	  rtot_nsm(maxd_acomp,nmaxd_nsm),   &
      	  vtot_nsm(nmaxd_nsm), xntot_nsm(nmaxd_nsm)

	real dgnum_nsm(nmaxd_nsm), sigmag_nsm(nmaxd_nsm)
	real aaprof_nsm(maxd_acomp+1,nmaxd_nsm)
	real bbprof_nsm(nmaxd_nsm)

	character*80 msg
	character*10 dumname


!   check for on/off
	if (mosaic_init_wrf_mixrats_flagaa .le. 0) return


!   *** currently only works for ntype_aer = 1
	itype = 1
	iphase = ai_phase
	m = 1

!   set values for initialization modes
	iiprof_nsm(:) = 1
	aaprof_nsm(:,:) = 0.0
	bbprof_nsm(:) = 0.0

	ntot_nsm = 4
	ntot_nsm = min( ntot_nsm, nsize_aer(itype) )

	lldum_so4 = 0
	lldum_nh4 = 0
	lldum_oc  = 0
	lldum_bc  = 0
	lldum_oin = 0
	lldum_na  = 0
	lldum_cl  = 0
	lldum_hysw = 0

	do ll = 1, ncomp_plustracer_aer(itype)
	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
		lptr_so4_aer(m,itype,iphase)) lldum_so4 = ll
	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
		lptr_nh4_aer(m,itype,iphase)) lldum_nh4 = ll
	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
		lptr_oc_aer(m,itype,iphase))  lldum_oc  = ll
	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
		lptr_bc_aer(m,itype,iphase))  lldum_bc  = ll
	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
		lptr_oin_aer(m,itype,iphase)) lldum_oin = ll
	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
		lptr_na_aer(m,itype,iphase))  lldum_na  = ll
	    if (massptr_aer(ll,m,itype,iphase) .eq.   &
		lptr_cl_aer(m,itype,iphase))  lldum_cl  = ll
	end do
	if (hyswptr_aer(m,itype) .gt. 0)   &
		lldum_hysw = ncomp_plustracer_aer(itype) + 1

	msg = ' '
	if (lldum_so4 .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_so4 = 0'
	if (lldum_nh4 .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_nh4 = 0'
	if (lldum_oc .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_oc = 0'
	if (lldum_bc .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_bc = 0'
	if (lldum_oin .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_oin = 0'
	if (lldum_na .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_na = 0'
	if (lldum_cl .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_cl = 0'
	if (lldum_hysw .le. 0)   &
		msg = '*** subr mosaic_init_wrf_mixrats - lldum_hysw = 0'
	if (msg .ne. ' ') call peg_error_fatal( lunerr, msg )


	do nsm = 1, ntot_nsm

	if (nsm .eq. 1) then
!   accum mode with so4, nh4, oc, bc
	    dgnum_nsm( nsm) = 0.15e-4
	    sigmag_nsm(nsm) = 2.0
	    aaprof_nsm(lldum_so4,nsm) = 0.50
	    aaprof_nsm(lldum_nh4,nsm) = aaprof_nsm(lldum_so4,nsm)   &
					* (mw_nh4_aer/mw_so4_aer)
	    aaprof_nsm(lldum_oc,nsm) = 0.25
	    aaprof_nsm(lldum_bc,nsm) = 0.05
	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_so4,nsm) * 0.2

	else if (nsm .eq. 2) then
!   aitken mode with so4, nh4, oc, bc
	    dgnum_nsm( nsm) = 0.03e-4
	    sigmag_nsm(nsm) = 2.0
	    aaprof_nsm(lldum_so4,nsm) = 0.50 * 0.020
	    aaprof_nsm(lldum_nh4,nsm) = aaprof_nsm(lldum_so4,nsm)   &
					* (mw_nh4_aer/mw_so4_aer)
	    aaprof_nsm(lldum_oc,nsm) = 0.25 * 0.020
	    aaprof_nsm(lldum_bc,nsm) = 0.05 * 0.020
	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_so4,nsm) * 0.2

	else if (nsm .eq. 3) then
!   coarse dust mode with oin
	    dgnum_nsm( nsm) = 1.0e-4
	    sigmag_nsm(nsm) = 2.0
	    aaprof_nsm(lldum_oin,nsm) = 0.5
	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm( 9,nsm) * 1.0e-3

	else if (nsm .eq. 4) then
!   coarse seasalt mode with na, cl
	    dgnum_nsm( nsm) = 2.0e-4
	    sigmag_nsm(nsm) = 2.0
	    aaprof_nsm(lldum_cl,nsm) = 0.1
	    aaprof_nsm(lldum_na,nsm) = aaprof_nsm(lldum_cl,nsm)   &
					* (mw_na_aer/mw_cl_aer)
	    aaprof_nsm(lldum_hysw,nsm) = aaprof_nsm(lldum_cl,nsm) * 0.2

	end if

	end do

!   when iflagaa = 1/2/3/4, use only the nsm-mode = iflagaa
	if (iflagaa .gt. 0) then
	    do nsm = 1, ntot_nsm
		if (nsm .ne. iflagaa) aaprof_nsm(:,nsm) = 0.0
	    end do
	end if



!
!   do the initialization now
!

!   calculate mode parameters
	do nsm = 1, ntot_nsm
	    sx_nsm(nsm) = alog( sigmag_nsm(nsm) )
	    sxr2_nsm(nsm) = sx_nsm(nsm) * sqrt(2.0)
	    x0_nsm(nsm) = alog( dgnum_nsm(nsm) )
	    x3_nsm(nsm) = x0_nsm(nsm) + 3.0*sx_nsm(nsm)*sx_nsm(nsm)
	end do

!   initialize rclm array to zero
	rclm(:,:) = 0.
!	rclmsvaa(:,:) = 0.

!
!   loop over all vertical levels
!
!	do 12900 k = 1, ktot
	do 12900 k = 1, 1

!	pdummb = 1013.*scord(k)
!	zdumkm = ptoz( pdummb ) * 1.e-3
	zdumkm = 0.0


!   for each species and nsm mode, define total mixing ratio
!	(for all sizes) at level k
!
!   iiprof_nsm =  +1 gives constant mixing ratio
!	aaprof_nsm(l,nsm) = constant mixing ratio (ppb)
!   iiprof_nsm =  +2 gives exponential profile
!	aaprof_nsm(l,nsm) = surface mixing ratio (ppb)
!	bbprof_nsm(l) = scale height (km)
!   iiprof_nsm =  +3 gives linear profile (then zero at z > zscalekm)
!	aaprof_nsm(l,nsm) = surface mixing ratio (ppb)
!	bbprof_nsm(l) = height (km) at which mixing ratio = 0

	do nsm = 1, ntot_nsm

	if (iiprof_nsm(nsm) .eq. 2) then
	    zscalekm = bbprof_nsm(nsm)
	    zfactor = exp( -zdumkm/zscalekm )
	else if (iiprof_nsm(nsm) .eq. 3) then
	    zscalekm = bbprof_nsm(nsm)
	    zfactor = max( 0., (1. - zdumkm/zscalekm) )
	else
	    zfactor = 1.0
	end if

	    do ll = 1, ncomp_plustracer_aer(itype) + 1
		rtot_nsm(ll,nsm) = aaprof_nsm(ll,nsm) * zfactor
	    end do

	end do

!   compute total volume and number mixing ratios for each nsm mode
!   rtot_nsm is ug/m3,  1.0e-6*rtot is g/m3,  vtot_nsm is cm3/m3
	do nsm = 1, ntot_nsm
	    dumvol = 0.
	    do ll = 1, ncomp_aer(itype)
		dum = 1.0e-6*rtot_nsm(ll,nsm)/dens_aer(ll,itype)
		dumvol = dumvol + max( 0., dum )
	    end do
	    vtot_nsm(nsm) = dumvol
	end do

!   now compute species and number mixing ratios for each bin
	do 12700 m = 1, nsize_aer(itype)

	vtot_nsm_ofmode(m) = 0.0

	do 12500 nsm = 1, ntot_nsm

!   for nsm_mode = n, compute fraction of number and volume
!   that is in bin m
	xlo = alog( dlo_sect(m,itype) )
	xhi = alog( dhi_sect(m,itype) )

	tlo = (xlo - x3_nsm(nsm))/sxr2_nsm(nsm)
	thi = (xhi - x3_nsm(nsm))/sxr2_nsm(nsm)
	if (tlo .ge. 0.) then
!	    fracvol = 0.5*( erfc(tlo) - erfc(thi) )
	    fracvol = 0.5*( erfc_num_recipes(tlo) - erfc_num_recipes(thi) )
	else
!	    fracvol = 0.5*( erfc(-thi) - erfc(-tlo) )
	    fracvol = 0.5*( erfc_num_recipes(-thi) - erfc_num_recipes(-tlo) )
	end if
	fracvol = max( fracvol, 0.0 )

	vtot_nsm_ofmode(m) = vtot_nsm_ofmode(m)  + vtot_nsm(nsm)*fracvol

!   now load that fraction of species-mixing-ratio
!   into the appropriate rclm location
	do ll = 1, ncomp_plustracer_aer(itype)
	    rclm( k, massptr_aer(ll,m,itype,iphase) ) =   &
      	    rclm( k, massptr_aer(ll,m,itype,iphase) ) +   &
      		fracvol*rtot_nsm(ll,nsm)
	end do

	if ((iphase .eq. ai_phase) .and.   &
	    (lldum_hysw .gt. 0) .and.   &
	    (hyswptr_aer(m,itype) .gt. 0)) then

	    rclm( k, hyswptr_aer(m,itype) ) =   &
      	    rclm( k, hyswptr_aer(m,itype) ) +   &
      		fracvol*rtot_nsm(lldum_hysw,nsm)
	end if

12500	continue

!   now compute number from volume
	dum = sqrt( dlo_sect(m,itype)*dhi_sect(m,itype) )
	dumvol1p = (pi/6.0)*(dum**3)
	rclm( k, numptr_aer(m,itype,iphase) ) = vtot_nsm_ofmode(m)/dumvol1p

!   set water = hyswatr
	if ((iphase .eq. ai_phase) .and.   &
	    (lldum_hysw .gt. 0) .and.   &
	    (hyswptr_aer(m,itype) .gt. 0) .and.   &
	    (waterptr_aer(m,itype) .gt. 0)) then

      	    rclm( k, waterptr_aer(m,itype) ) =    &
	    rclm( k, hyswptr_aer(m,itype) )
	end if

12700	continue

12900	continue


!
!   do diagnostic output
!

! temporary
! temporary

	dumarr(:) = 0.0
	msg = ' '
	call peg_message( lunout, msg )
	msg = '*** subr mosaic_init_wrf_mixrats_opt1 results'
	call peg_message( lunout, msg )
	msg = '    mass in ug/m3     number in #/m3     volume in cm3/m3'
	call peg_message( lunout, msg )
	msg = ' '
	call peg_message( lunout, msg )
	msg = ' mode  l  l1  species      conc'
	call peg_message( lunout, msg )

        do 14390 mdum = 1, nsize_aer(itype)+1
	m = min( mdum, nsize_aer(itype) )
	msg = ' '
	call peg_message( lunout, msg )
        do 14350 l = 1, ncomp_plustracer_aer(itype)+4

            if (l .le. ncomp_plustracer_aer(itype)) then
                l1 = massptr_aer(l,m,itype,iphase)
		dumname = name_aer(l,itype)
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+1) then
                l1 = hyswptr_aer(m,itype)
		dumname = 'hystwatr'
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+2) then
                l1 = waterptr_aer(m,itype)
		dumname = 'water'
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+3) then
                l1 = numptr_aer(m,itype,iphase)
		dumname = 'number'
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+4) then
                l1 = 0
		dumname = 'volume'
		dum = vtot_nsm_ofmode(m)
	    else
		dumname = '=BADBAD='
		l1 = -1
		dum = -1.0
            end if

	    l1dum = l1
	    if (aboxtest_testmode .gt. 0) l1dum = max( l1-1, 0 )

	    if (mdum .le. nsize_aer(itype)) then
		dumarr(l) = dumarr(l) + dum
		write(msg,9620) m, l, l1dum, dumname, dum
	    else
		write(msg,9625) l, dumname, dumarr(l)
	    end if
	    call peg_message( lunout, msg )

14350   continue
14390   continue

9620    format( 3i4, 2x, a, 3(1pe12.3) )
9625    format( ' sum', i4, '   -', 2x, a, 3(1pe12.3) )


!
!   load the chem array
!
        do 16390 m = 1, nsize_aer(itype)
        do 16350 l = 1, 15

	    if (l .eq. 1) then
		l1 = lptr_so4_aer(m,itype,iphase)
	    else if (l .eq. 2) then
		l1 = lptr_no3_aer(m,itype,iphase)
	    else if (l .eq. 3) then
		l1 = lptr_cl_aer(m,itype,iphase)
	    else if (l .eq. 4) then
		l1 = lptr_msa_aer(m,itype,iphase)
	    else if (l .eq. 5) then
		l1 = lptr_co3_aer(m,itype,iphase)
	    else if (l .eq. 6) then
		l1 = lptr_nh4_aer(m,itype,iphase)
	    else if (l .eq. 7) then
		l1 = lptr_na_aer(m,itype,iphase)
	    else if (l .eq. 8) then
		l1 = lptr_ca_aer(m,itype,iphase)
	    else if (l .eq. 9) then
		l1 = lptr_oin_aer(m,itype,iphase)
	    else if (l .eq. 10) then
		l1 = lptr_oc_aer(m,itype,iphase)
	    else if (l .eq. 11) then
		l1 = lptr_bc_aer(m,itype,iphase)
	    else if (l .eq. 12) then
		l1 = hyswptr_aer(m,itype)
	    else if (l .eq. 13) then
		l1 = waterptr_aer(m,itype)
	    else if (l .eq. 14) then
		l1 = numptr_aer(m,itype,iphase)
	    else
		goto 16350
	    end if
	    l3 = l1

	    if ((l1 .gt. 0) .and. (l1 .le. lmaxd) .and.   &
	        (l3 .ge. param_first_scalar)) then
		do it = its, ite
!   *** note:  not sure what the kt limits should be here 
		do kt = kts, kte-1
		do jt = jts, jte
		    chem(it,kt,jt,l3) = rclm(1,l1)
		end do
		end do
		end do
	    end if

16350   continue
16390   continue


!   all done
	return
	end subroutine mosaic_init_wrf_mixrats_opt1


!-----------------------------------------------------------------------
!wig 10-May-2004, added phb, ph, g
	subroutine mosaic_init_wrf_mixrats_opt2(   &
		iflagaa, config_flags,             &
		chem, z_at_w, g,                   &
		ids,ide, jds,jde, kds,kde,         &
		ims,ime, jms,jme, kms,kme,         &
		its,ite, jts,jte, kts,kte          )

!
!   provides initial values for mosaic aerosol species (mass and number
!	mixing ratio) for "Texas August 2000" simulations
!   modified to use different profiles for different aerosols for NEAQS case, egc 7/2005
!   currently all the initial values are uniform in x, y, AND z
!
! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
!     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now
!
	use module_configure, only:  grid_config_rec_type
	use module_state_description, only:  num_chem, param_first_scalar
	use module_data_mosaic_asect
	use module_data_mosaic_other
	use module_peg_util, only:  peg_message, peg_error_fatal

	implicit none

!   subr arguments
	type(grid_config_rec_type), intent(in) :: config_flags

	integer, intent(in) ::   &
		iflagaa,   &
		ids, ide, jds, jde, kds, kde,   &
		ims, ime, jms, jme, kms, kme,   &
		its, ite, jts, jte, kts, kte

	real, intent(inout),     &
		dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
		chem

        real, intent(in),        &
		dimension( ims:ime, kms:kme, jms:jme ) :: &
                z_at_w
        real :: g

!   local variables
	integer l, l1, l3, m, mdum
	integer iphase, itype
	integer it, jt, kt

!wig 10-May-2004, added mult
	real dum, dumvol1p, mult
	real qcoar, qfine, qval

	real :: vtot_ofmode(maxd_asize)
	real :: dumarr(maxd_acomp+5)
	real :: fr_coar(8), fr_fine(8)

!wig 01-Apr-2005, Updated fractional breakdown between bins. (See also
!                 bdy_chem_value_mosaic in this file and mosaic_addemiss in 
!                 module_mosaic_addemiss.F) Note that the values here no
!                 longer match those in mosaic_addemiss.
!rce 10-May-2005, changed fr8b_coar to values determined by jdf
	real, save :: fr8b_coar(8) =   &
        (/ 0.0,  0.0,  0.0,  0.019,  0.0212,  0.1101,  0.2751, 0.7018 /) ! jdf Testbed
!jdf    (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.300, 0.700 /) ! 10-May-2005
!       (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.933, 0.067 /) ! 01-Apr-2005
!       (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.16,  0.84  /) ! "old"

	real, save :: fr8b_fine(8) =   &
        (/ 0.0245, 0.1472, 0.3501, 0.3317, 0.1251, 0.0186, 0.0011, 0.0/) !jdf Testbed
!jdf    (/ 0.006, 0.024, 0.231, 0.341, 0.140, 0.258, 0., 0./) !5-Apr-2005 values
!       (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
!	(/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
!	(/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.00, 0.00 /)

	real :: qfine_so4, qfine_no3, qfine_cl, qfine_msa,   &
		qfine_co3, qfine_nh4, qfine_na, qfine_ca, qfine_oin,   &
		qfine_oc, qfine_bc, qfine_hysw, qfine_watr, qfine_vol
	real :: qcoar_so4, qcoar_no3, qcoar_cl, qcoar_msa,   &
		qcoar_co3, qcoar_nh4, qcoar_na, qcoar_ca, qcoar_oin,   &
		qcoar_oc, qcoar_bc, qcoar_hysw, qcoar_watr, qcoar_vol

!wig 10-May-2004, added z
        real, dimension( ims:ime, kms:kme, jms:jme ) :: z

	character*80 msg
	character*10 dumname


!   *** currently only works for ntype_aer = 1
	itype = 1
	iphase = ai_phase
	m = 1

!wig 10-May-2004, block begin...
! calculate the mid-level heights
    do jt = jts, min(jte,jde-1)
       do kt = kts, kte-1
          do it = its, min(ite,ide-1)
             z(it,kt,jt) = (z_at_w(it,kt,jt)+z_at_w(it,kt+1,jt))*0.5
          end do
       end do
    end do
!wig 10-May-2004, ...end block

!   set the partitioning fractions for either 8 or 4 bins
	if (nsize_aer(itype) == 8) then
	    fr_coar(:) = fr8b_coar(:)
	    fr_fine(:) = fr8b_fine(:)
	else if (nsize_aer(itype) == 4) then
	    do m = 1, nsize_aer(itype)
		fr_coar(m) = fr8b_coar(2*m) + fr8b_coar(2*m-1)
		fr_fine(m) = fr8b_fine(2*m) + fr8b_fine(2*m-1)
	    end do
	else
	    write(msg,'(a,i5)')   &
		'subr mosaic_init_wrf_mixrats_opt2' //   &
		' - nsize_aer(itype) must be 4 or 8 but = ', nsize_aer(itype)
	    call peg_error_fatal( lunout, msg )
	end if

!
!   compute initial values (currently uniform in x, y, AND z)
!   load them into the chem array
!   also load them into the rclm array for diagnostic output
!
	rclm(:,:) = 0.0
	vtot_ofmode(:) = 0.0

!   Set values for fine and coarse mass concentrations, and
!   compute fine/coarse volume concentrations. These are also set in
!   bdy_chem_value_mosaic, below.
!   (these are latest values from 1-Apr-2005 discussions)
!
!   rce 10-may-2005 - changed qfine_cl, _na, _oin to values determined by jdf
	qfine_so4 = 2.14
	qcoar_so4 = 0.242
	qfine_no3 = 0.11
	qcoar_no3 = 0.03
!	qfine_cl  = 0.3
	qfine_cl  = 0.14     ! 10-May-2005
	qcoar_cl  = 0.139
	qfine_msa = 0.0
	qcoar_msa = 0.0
	qfine_co3 = 0.0
	qcoar_co3 = 0.0
	qfine_nh4 = 0.83
	qcoar_nh4 = 0.10
!	qfine_na  = 0.2
	qfine_na  = 0.1      ! 10-May-2005
	qcoar_na  = 0.09
	qfine_ca  = 0.0
	qcoar_ca  = 0.0
!	qfine_oin = 6.2
	qfine_oin = 3.48     ! 10-May-2005
	qcoar_oin = 0.35
	qfine_oc  = 1.00
	qcoar_oc  = 1.50
	qfine_bc  = 0.2
	qcoar_bc  = 0.075
	qfine_hysw = 0.0
	qcoar_hysw = 0.0
	qfine_watr = 0.0
	qcoar_watr = 0.0
#if (CASENAME == 2)
!       qfine_so4 = 0.50
!       qcoar_so4 = 0.50/2.0
!       qfine_no3 = 0.75
!       qcoar_no3 = 0.75/2.0
!       qfine_cl  = 0.10
!       qcoar_cl  = 0.10/2.0
!       qfine_nh4 = 0.30
!       qcoar_nh4 = 0.30/2.0
!       qfine_na  = 0.20
!       qcoar_na  = 0.20/2.0
!       qfine_oin = 4.00
!       qcoar_oin = 4.00/2.0
!       qfine_oc  = 3.00
!       qcoar_oc  = 3.00/2.0
!       qfine_bc  = 0.50
!       qcoar_bc  = 0.50/2.0
        qfine_so4 = 0.300
        qcoar_so4 = 0.300/2.0
        qfine_no3 = 0.001
        qcoar_no3 = 0.001/2.0
        qfine_cl  = 1.050
        qcoar_cl  = 1.05/02.0
        qfine_nh4 = 0.094
        qcoar_nh4 = 0.094/2.0
        qfine_na  = 0.700
        qcoar_na  = 0.700/2.0
        qfine_oin = 4.500
        qcoar_oin = 4.500/2.0
        qfine_oc  = 0.088
        qcoar_oc  = 0.088/2.0
        qfine_bc  = 0.013
        qcoar_bc  = 0.013/2.0
#endif
#if (CASENAME == 3)
        qfine_so4 = 0.300
        qcoar_so4 = 0.0
        qfine_nh4 = 0.094
        qcoar_nh4 = 0.0
        qfine_no3 = 0.001
        qcoar_no3 = 0.0
        qfine_cl  = 0.0
        qcoar_cl  = 1.050
        qfine_na  = 0.000
        qcoar_na  = 0.700
        qfine_oin = 4.500
        qcoar_oin = 4.500/2.0
        qfine_oc  = 0.088
        qcoar_oc  = 0.0
        qfine_bc  = 0.013
        qcoar_bc  = 0.0
#endif

!!$! old values from 23-Apr-2004:
!!$	qfine_so4 = 2.554
!!$	qcoar_so4 = 0.242
!!$	qfine_no3 = 0.07
!!$	qcoar_no3 = 0.03
!!$	qfine_cl  = 0.324
!!$	qcoar_cl  = 0.139
!!$	qfine_msa = 0.0
!!$	qcoar_msa = 0.0
!!$	qfine_co3 = 0.0
!!$	qcoar_co3 = 0.0
!!$	qfine_nh4 = 0.98
!!$	qcoar_nh4 = 0.10
!!$	qfine_na  = 0.21
!!$	qcoar_na  = 0.09
!!$	qfine_ca  = 0.0
!!$	qcoar_ca  = 0.0
!!$	qfine_oin = 0.15
!!$	qcoar_oin = 0.35
!!$	qfine_oc  = 1.00
!!$	qcoar_oc  = 1.50
!!$	qfine_bc  = 0.175
!!$	qcoar_bc  = 0.075
!!$	qfine_hysw = 0.0
!!$	qcoar_hysw = 0.0
!!$	qfine_watr = 0.0
!!$	qcoar_watr = 0.0


! qfine_so4 ... are ug/m3 so 1.0e-6 factor gives g/m3
! dens_so4  ... are g/cm3;  qfine_vol ... are cm3/m3
	qfine_vol = 1.0e-6 * (   &
	    (qfine_so4/dens_so4_aer) + (qfine_no3/dens_no3_aer) +   &
	    (qfine_cl /dens_cl_aer ) + (qfine_msa/dens_msa_aer) +   &
	    (qfine_co3/dens_co3_aer) + (qfine_nh4/dens_nh4_aer) +   &
	    (qfine_na /dens_na_aer ) + (qfine_ca /dens_ca_aer ) +   &
	    (qfine_oin/dens_oin_aer) + (qfine_oc /dens_oc_aer ) +   &
	    (qfine_bc /dens_bc_aer ) )
	qcoar_vol =  1.0e-6 * (  &
	    (qcoar_so4/dens_so4_aer) + (qcoar_no3/dens_no3_aer) +   &
	    (qcoar_cl /dens_cl_aer ) + (qcoar_msa/dens_msa_aer) +   &
	    (qcoar_co3/dens_co3_aer) + (qcoar_nh4/dens_nh4_aer) +   &
	    (qcoar_na /dens_na_aer ) + (qcoar_ca /dens_ca_aer ) +   &
	    (qcoar_oin/dens_oin_aer) + (qcoar_oc /dens_oc_aer ) +   &
	    (qcoar_bc /dens_bc_aer ) )

        do 2900 m = 1, nsize_aer(itype)
        do 2800 l = 1, 15

	    if (l .eq. 1) then
		l1 = lptr_so4_aer(m,itype,iphase)
		qfine = qfine_so4
		qcoar = qcoar_so4
	    else if (l .eq. 2) then
		l1 = lptr_no3_aer(m,itype,iphase)
		qfine = qfine_no3
		qcoar = qcoar_no3
	    else if (l .eq. 3) then
		l1 = lptr_cl_aer(m,itype,iphase)
		qfine = qfine_cl
		qcoar = qcoar_cl
	    else if (l .eq. 4) then
		l1 = lptr_msa_aer(m,itype,iphase)
		qfine = qfine_msa
		qcoar = qcoar_msa
	    else if (l .eq. 5) then
		l1 = lptr_co3_aer(m,itype,iphase)
		qfine = qfine_co3
		qcoar = qcoar_co3
	    else if (l .eq. 6) then
		l1 = lptr_nh4_aer(m,itype,iphase)
		qfine = qfine_nh4
		qcoar = qcoar_nh4
	    else if (l .eq. 7) then
		l1 = lptr_na_aer(m,itype,iphase)
		qfine = qfine_na
		qcoar = qcoar_na
	    else if (l .eq. 8) then
		l1 = lptr_ca_aer(m,itype,iphase)
		qfine = qfine_ca
		qcoar = qcoar_ca
	    else if (l .eq. 9) then
		l1 = lptr_oin_aer(m,itype,iphase)
		qfine = qfine_oin
		qcoar = qcoar_oin
	    else if (l .eq. 10) then
		l1 = lptr_oc_aer(m,itype,iphase)
		qfine = qfine_oc
		qcoar = qcoar_oc
	    else if (l .eq. 11) then
		l1 = lptr_bc_aer(m,itype,iphase)
		qfine = qfine_bc
		qcoar = qcoar_bc
	    else if (l .eq. 12) then
		l1 = hyswptr_aer(m,itype)
		qfine = qfine_hysw
		qcoar = qcoar_hysw
	    else if (l .eq. 13) then
		l1 = waterptr_aer(m,itype)
		qfine = qfine_watr
		qcoar = qcoar_watr
	    else if (l .eq. 14) then
		l1 = numptr_aer(m,itype,iphase)
		dumvol1p = sqrt(volumlo_sect(m,itype)*volumhi_sect(m,itype))
		qfine = qfine_vol/dumvol1p
		qcoar = qcoar_vol/dumvol1p
		vtot_ofmode(m) =   &
			qfine_vol*fr_fine(m) + qcoar_vol*fr_coar(m)
	    else
		goto 2800
	    end if
	    l3 = l1

	    if ((l1 .gt. 0) .and. (l1 .le. lmaxd) .and.   &
	        (l3 .ge. param_first_scalar)) then
		qval = qfine*fr_fine(m) + qcoar*fr_coar(m)
		rclm(1,l1) = qval

		do jt = jts, min(jte,jde-1)
		do kt = kts, kte-1
		do it = its, min(ite,ide-1)

!wig 28-May-2004, begin block...
! Determine height multiplier...
! This should mimic the calculation in sorgam_set_aer_bc_pnnl,
! sorgam_init_aer_ic_pnnl, bdy_chem_value_mosaic
!!$!    Height(m)     Multiplier
!!$!    ---------     ----------
!!$!    <=2000        1.0
!!$!    2000<z<3000   linear transition zone to 0.5
!!$!    3000<z<5000   linear transition zone to 0.25
!!$!    >=5000        0.25
!!$!
!!$! which translates to:
!!$!    2000<z<3000   mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
!!$!    3000<z<5000   mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
!!$! or in reduced form:
!!$		   if( z(it,kt,jt) <= 2000. ) then
!!$		      mult = 1.0
!!$		   elseif( z(it,kt,jt) > 2000. &
!!$                        .and. z(it,kt,jt) <= 3000. ) then
!!$		      mult = 1.0 - 0.0005*(z(it,kt,jt)-2000.)
!!$		   elseif( z(it,kt,jt) > 3000. &
!!$                        .and. z(it,kt,jt) <= 5000. ) then
!!$		      mult = 0.5 - 1.25e-4*(z(it,kt,jt)-3000.)
!!$		   else
!!$		      mult = 0.25
!!$           end if
! Updated 1-Apr-2005:
#if (CASENAME == 1)
! further updated 7-20-05 egc:   species specific profiles based on MIRAG2 output
         mult = 1.0
         if ( (l3 == 1) .or. (l3 == 2) .or. (l3 == 6) ) then
! so4, no3 and nh4 aerosol profiles
	     if ( z(it,kt,jt) <= 1000. ) then
	              mult = 1.0
		  elseif( z(it,kt,jt) > 1000. &
                        .and. z(it,kt,jt) <= 4000. ) then
		      mult = 1.0 - 2.367e-4*(z(it,kt,jt)-1000.)
		   elseif( z(it,kt,jt) > 4000. &
                        .and. z(it,kt,jt) <= 5500. ) then
		      mult = 0.29 - 4.0e-5*(z(it,kt,jt)-4000.)
		   else
		      mult = 0.23
               end if
           else if ( ( l3 == 3) .or. (l3 ==7) ) then
! na and cl aerosol profiles
	     if ( z(it,kt,jt) <= 100. ) then
	               mult = 1.0
		  elseif( z(it,kt,jt) > 100. &
                        .and. z(it,kt,jt) <= 265. ) then
		      mult = 1.0 - 2.9e-3*(z(it,kt,jt)-100.)
		   elseif( z(it,kt,jt) > 265. &
                        .and. z(it,kt,jt) <= 2000. ) then
		      mult = 0.52 - 2.94e-4*(z(it,kt,jt)-265.)
                   elseif( z(it,kt,jt) > 2000. &
                        .and. z(it,kt,jt) <= 7000. ) then
		       mult = 0.01
		   else
		      mult = 1.e-10
               end if
          else if  ( l3 == 10)  then
! oc aerosol profile
	     if ( z(it,kt,jt) <= 600. ) then
	               mult = 1.0
		  elseif( z(it,kt,jt) > 600. &
                        .and. z(it,kt,jt) <= 3389. ) then
		      mult = 1.0 - 2.04e-4*(z(it,kt,jt)-600.)
		   elseif( z(it,kt,jt) > 3389. &
                        .and. z(it,kt,jt) <= 8840. ) then
		      mult = 0.43 - 7.045e-5*(z(it,kt,jt)-3389.)
		   else
		      mult = 0.046
               end if
          else if  ( l3 == 11)  then
! bc aerosol profile
	     if ( z(it,kt,jt) <= 100. ) then
	               mult = 1.0
		  elseif( z(it,kt,jt) > 100. &
                        .and. z(it,kt,jt) <= 3400. ) then
		      mult = 1.0 - 2.51e-4*(z(it,kt,jt)-100.)
		   elseif( z(it,kt,jt) > 3400. &
                        .and. z(it,kt,jt) <= 8400. ) then
		      mult = 0.172 -2.64e-5*(z(it,kt,jt)-3400.)
		   else
		      mult = 0.04
               end if
         else if  ( l3 == 9)  then
! oin aerosol profile
	     if ( z(it,kt,jt) <= 1580. ) then
	               mult = 1.0 + 1.77e-4 *z(it,kt,jt)
		  elseif( z(it,kt,jt) > 1580. &
                        .and. z(it,kt,jt) <= 5280. ) then
		      mult = 1.28 - 2.65e-4*(z(it,kt,jt)-1580.)
		   else
		      mult = 0.30
               end if
	else
! generic profile for other other species (which should have groundlevel values=0)
#endif
!    Height(m)     Multiplier
!    ---------     ----------
!    <=2000        1.0
!    2000<z<3000   linear transition zone to 0.25
!    3000<z<5000   linear transision zone to 0.125
!    >=5000        0.125
!
! which translates to:
!    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
!    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
! or in reduced form:
		   if( z(it,kt,jt) <= 2000. ) then
		      mult = 1.0
		   elseif( z(it,kt,jt) > 2000. &
                        .and. z(it,kt,jt) <= 3000. ) then
		      mult = 1.0 - 0.00075*(z(it,kt,jt)-2000.)
		   elseif( z(it,kt,jt) > 3000. &
                        .and. z(it,kt,jt) <= 5000. ) then
		      mult = 0.25 - 4.166666667e-5*(z(it,kt,jt)-3000.)
		   else
		      mult = 0.125
           end if
#if (CASENAME == 1)
        end if				!close l3 comparison block
#endif
#if (CASENAME == 2)
!     if( z(it,kt,jt) <= 3000. ) then
!       mult = 1.0
!     elseif( z(it,kt,jt) > 3000. &
!       .and. z(it,kt,jt) <= 5000. ) then
!       mult = 1.0 - 0.0004375*(z(it,kt,jt)-3000.)
!     else
!        mult = 0.125
!     end if
!
!     if( z(it,kt,jt) <= 5000. ) then
!       mult = 1.0
!     elseif( z(it,kt,jt) > 5000. &
!       .and. z(it,kt,jt) <= 7000. ) then
!       mult = 1.0 - 0.0004375*(z(it,kt,jt)-5000.)
!     else
!        mult = 0.125
!     end if
      mult=1.0
      if ( (l3 == lptr_so4_aer(m,itype,iphase)) .or. &
           (l3 == lptr_no3_aer(m,itype,iphase)) .or. &
           (l3 == lptr_nh4_aer(m,itype,iphase)) ) then
! so4, no3 and nh4 aerosol profiles
        if( z(it,kt,jt) <= 500. ) then
          mult = 1.000
        elseif( z(it,kt,jt) > 500. &
          .and. z(it,kt,jt) <= 1000. ) then
          mult = 1.000 - 0.0010740*(z(it,kt,jt)-500.)
        elseif( z(it,kt,jt) > 1000. &
          .and. z(it,kt,jt) <= 5000. ) then
          mult = 0.463 - 0.0001110*(z(it,kt,jt)-1000.)
        elseif( z(it,kt,jt) > 5000. &
          .and. z(it,kt,jt) <= 20000. ) then
          mult = 0.019
        else
          mult = 0.019
        end if
      else if ( (l3 == lptr_na_aer(m,itype,iphase)) .or. &
                (l3 == lptr_cl_aer(m,itype,iphase)) ) then
! na and cl aerosol profiles
        if( z(it,kt,jt) <= 500. ) then
          mult = 1.000
        elseif( z(it,kt,jt) > 500. &
          .and. z(it,kt,jt) <= 1000. ) then
          mult = 1.000 - 0.0014260*(z(it,kt,jt)-500.)
        elseif( z(it,kt,jt) > 1000. &
          .and. z(it,kt,jt) <= 5000. ) then
          mult = 0.287 - 0.0000643*(z(it,kt,jt)-1000.)
        elseif( z(it,kt,jt) > 5000. &
          .and. z(it,kt,jt) <= 10000. ) then
          mult = 0.030
        elseif( z(it,kt,jt) > 10000. &
          .and. z(it,kt,jt) <= 20000. ) then
          mult = 0.030 - 0.0000028*(z(it,kt,jt)-10000.)
        else
          mult = 0.002
        end if
      else if ( (l3 == lptr_oc_aer(m,itype,iphase)) .or. &
                (l3 == lptr_bc_aer(m,itype,iphase)) ) then
! oc and bc aerosol profile
        if( z(it,kt,jt) <= 500. ) then
          mult = 1.000
        elseif( z(it,kt,jt) > 500. &
          .and. z(it,kt,jt) <= 1000. ) then
          mult = 1.000 - 0.0002500*(z(it,kt,jt)-500.)
        elseif( z(it,kt,jt) > 1000. &
          .and. z(it,kt,jt) <= 5000. ) then
          mult = 0.875 - 0.0001843*(z(it,kt,jt)-1000.)
        elseif( z(it,kt,jt) > 5000. &
          .and. z(it,kt,jt) <= 20000. ) then
          mult = 0.138
        else
          mult = 0.138
        end if
      else if ( l3 == lptr_oin_aer(m,itype,iphase) ) then
! oin aerosol profile
        if( z(it,kt,jt) <= 500. ) then
          mult = 1.000
        elseif( z(it,kt,jt) > 500. &
          .and. z(it,kt,jt) <= 1000. ) then
          mult = 1.000 - 0.0007340*(z(it,kt,jt)-500.)
        elseif( z(it,kt,jt) > 1000. &
          .and. z(it,kt,jt) <= 5000. ) then
          mult = 0.633 - 0.0001103*(z(it,kt,jt)-1000.)
        elseif( z(it,kt,jt) > 5000. &
          .and. z(it,kt,jt) <= 10000. ) then
          mult = 0.192 - 0.0000080*(z(it,kt,jt)-5000.)
        elseif( z(it,kt,jt) > 10000. &
          .and. z(it,kt,jt) <= 20000. ) then
          mult = 0.152 - 0.0000137*(z(it,kt,jt)-10000.)
        else
          mult = 0.015
        end if
      else
        if( z(it,kt,jt) <= 5000. ) then
          mult = 1.0
        elseif( z(it,kt,jt) > 5000. &
          .and. z(it,kt,jt) <= 7000. ) then
          mult = 1.0 - 0.0004375*(z(it,kt,jt)-5000.)
        else
           mult = 0.125
        end if
      end if
#endif
#if (CASENAME == 3)
      mult = 1.000
      if( z(it,kt,jt) <= 500. ) then
        mult = 1.000
      elseif( z(it,kt,jt) > 500. &
        .and. z(it,kt,jt) <= 1000. ) then
        mult = 1.000 - 0.0010740*(z(it,kt,jt)-500.)
      elseif( z(it,kt,jt) > 1000. &
        .and. z(it,kt,jt) <= 5000. ) then
        mult = 0.463 - 0.0001110*(z(it,kt,jt)-1000.)
      else
        mult = 0.019
      end if
#endif

		    chem(it,kt,jt,l3) = mult*rclm(1,l1)
!wig 28-May-2004, ...end block
!		    chem(it,kt,jt,l3) = rclm(1,l1)
		end do
		end do
		end do
	    end if


2800   continue
2900   continue

!
!   do diagnostic output
!
	dumarr(:) = 0.0
	msg = ' '
	call peg_message( lunout, msg )
	msg = '*** subr mosaic_init_wrf_mixrats_opt2 results'
	call peg_message( lunout, msg )
	msg = '    mass in ug/m3     number in #/m3     volume in cm3/m3'
	call peg_message( lunout, msg )
	msg = ' '
	call peg_message( lunout, msg )
	msg = ' mode  l  l1  species      conc'
	call peg_message( lunout, msg )

        do 3190 mdum = 1, nsize_aer(itype)+1
	m = min( mdum, nsize_aer(itype) )
	msg = ' '
	call peg_message( lunout, msg )
        do 3150 l = 1, ncomp_plustracer_aer(itype)+4

            if (l .le. ncomp_plustracer_aer(itype)) then
                l1 = massptr_aer(l,m,itype,iphase)
		dumname = name_aer(l,itype)
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+1) then
                l1 = hyswptr_aer(m,itype)
		dumname = 'hystwatr'
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+2) then
                l1 = waterptr_aer(m,itype)
		dumname = 'water'
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+3) then
                l1 = numptr_aer(m,itype,iphase)
		dumname = 'number'
		dum = rclm(1,l1)
            else if (l .eq. ncomp_plustracer_aer(itype)+4) then
                l1 = 0
		dumname = 'volume'
		dum = vtot_ofmode(m)
	    else
		dumname = '=BADBAD='
		l1 = -1
		dum = -1.0
            end if

	    if (mdum .le. nsize_aer(itype)) then
		dumarr(l) = dumarr(l) + dum
		write(msg,9620) m, l, l1, dumname, dum
	    else
		write(msg,9625) l, dumname, dumarr(l)
	    end if
	    call peg_message( lunout, msg )

3150   continue
3190   continue

9620    format( 3i4, 2x, a, 3(1pe12.3) )
9625    format( ' sum', i4, '   -', 2x, a, 3(1pe12.3) )


!   all done
	return
	end subroutine mosaic_init_wrf_mixrats_opt2


!-----------------------------------------------------------------------
	real function erfc_num_recipes( x )
!
!   from press et al, numerical recipes, 1990, page 164
!
	implicit none
	real x
	double precision erfc_dbl, dum, t, zz

	zz = abs(x)
	t = 1.0/(1.0 + 0.5*zz)

!	erfc_num_recipes =
!     &	  t*exp( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
!     &	  t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
!     &	                                   t*(-1.13520398 +
!     &	  t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))

	dum =  ( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +   &
      	  t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +   &
      	                                   t*(-1.13520398 +   &
      	  t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))

	erfc_dbl = t * exp(dum)
	if (x .lt. 0.0) erfc_dbl = 2.0d0 - erfc_dbl

	erfc_num_recipes = erfc_dbl

	return
	end function erfc_num_recipes     


!-----------------------------------------------------------------------
	end module module_mosaic_initmixrats




!-----------------------------------------------------------------------
 	subroutine bdy_chem_value_mosaic ( chem_bv, alt, zz, nch, config_flags )
!
! provides boundary values for the mosaic aerosol species
!
! it is outside of the module declaration because of potential
!    module1 --> module2 --> module1 use conflicts
!
! rce 11-sep-2004 - eliminated use of the _wrfch pointers (lptr_xxx_a_wrfch, 
!     lwaterptr_wrfch, numptr_wrfch); use only the _amode pointers now
!
	use module_configure, only:  grid_config_rec_type
	use module_state_description, only:  param_first_scalar,   &
			aer_bc_pnnl
	use module_data_mosaic_asect
	use module_data_mosaic_other
	implicit none

! arguments
	REAL,    intent(OUT)  :: chem_bv    ! boundary value for chem(-,-,-,nch)
	REAL,    intent(IN)   :: alt        ! inverse density
	REAL,    intent(IN)   :: zz         ! height
	INTEGER, intent(IN)   :: nch        ! index number of chemical species
	TYPE (grid_config_rec_type), intent(in) :: config_flags

! local variables
	integer :: iphase, itype, m
    logical :: foundit

	real, parameter :: chem_bv_def = 1.0e-20
!wig 10-May-2004, added mult
	real :: dumvol1p, mult
        integer :: iflg
	real :: qcoar, qfine, qval

	real :: fr_coar(8), fr_fine(8)

!wig 1-Apr-2005,  Updated fractional breakdown between bins. (See also
!                 mosaic_init_wrf_mixrats_opt2, above, and mosaic_addemiss
!                 in module_mosaic_addemiss.F). Note that these values no
!                 longer match those in mosaic_addemiss.
	real, save :: fr8b_coar(8) =   &
        (/ 0.0,  0.0,  0.0,  0.019,  0.0212,  0.1101,  0.2751, 0.7018 /) ! jdf Testbed
!jdf    (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.300, 0.700 /) ! 10-May-2005
!       (/ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.933, 0.067 /) ! 01-Apr-2005
!		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.16, 0.84 /)
	real, save :: fr8b_fine(8) =   &
        (/ 0.0245, 0.1472, 0.3501, 0.3317, 0.1251, 0.0186, 0.0011, 0.0/) !jdf Testbed
!jdf    (/ 0.006, 0.024, 0.231, 0.341, 0.140, 0.258, 0., 0./) !5-Apr-2005 values
!       (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
!		(/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
!		(/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.00, 0.00 /)

	real :: qfine_so4, qfine_no3, qfine_cl, qfine_msa,   &
		qfine_co3, qfine_nh4, qfine_na, qfine_ca, qfine_oin,   &
		qfine_oc, qfine_bc, qfine_hysw, qfine_watr, qfine_vol
	real :: qcoar_so4, qcoar_no3, qcoar_cl, qcoar_msa,   &
		qcoar_co3, qcoar_nh4, qcoar_na, qcoar_ca, qcoar_oin,   &
		qcoar_oc, qcoar_bc, qcoar_hysw, qcoar_watr, qcoar_vol

	character*80 msg



! when aer_bc_opt /= aer_bc_pnnl,
! just chem_bv=near-zero for all species
	chem_bv = chem_bv_def
	if (config_flags%aer_bc_opt /= aer_bc_pnnl) return
	if (nch < param_first_scalar) return


!   *** currently only works for ntype_aer = 1
	itype = 1
	iphase = ai_phase
	m = 1 !This is here for size, but is overridden by loop below.


!
!   following is for aer_bc_opt == aer_bc_pnnl
!       when boundary values are set for Texas August 2000 simulations
!
!   set the partitioning fractions for either 8 or 4 bins
	if (nsize_aer(itype) == 8) then
	    fr_coar(:) = fr8b_coar(:)
	    fr_fine(:) = fr8b_fine(:)
	else if (nsize_aer(itype) == 4) then
	    do m = 1, nsize_aer(itype)
		fr_coar(m) = fr8b_coar(2*m) + fr8b_coar(2*m-1)
		fr_fine(m) = fr8b_fine(2*m) + fr8b_fine(2*m-1)
	    end do
	else
	    write(msg,'(a,i5)')   &
		'subr bdy_chem_value_mosaic' //   &
		' - nsize_aer(itype) must be 4 or 8 but = ', nsize_aer(itype)
	    call wrf_error_fatal( msg )
	end if

! Determine height multiplier...
! This should mimic the calculation in sorgam_set_aer_bc_pnnl,
! sorgam_init_aer_ic_pnnl, and mosaic_init_wrf_mixrats_opt2
!!$!    Height(m)     Multiplier
!!$!    ---------     ----------
!!$!    <=2000        1.0
!!$!    2000<z<3000   linear transition zone to 0.5
!!$!    3000<z<5000   linear transision zone to 0.25
!!$!    >=5000        0.25
!!$!
!!$! which translates to:
!!$!    2000<z<3000   mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
!!$!    3000<z<5000   mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
!!$! or in reduced form:
!!$        if( zz <= 2000. ) then
!!$           mult = 1.0
!!$        elseif( zz > 2000. &
!!$             .and. zz <= 3000. ) then
!!$           mult = 1.0 - 0.0005*(zz-2000.)
!!$        elseif( zz > 3000. &
!!$             .and. zz <= 5000. ) then
!!$           mult = 0.5 - 1.25e-4*(zz-3000.)
!!$        else
!!$           mult = 0.25
!!$        end if
#if (CASENAME == 1)
    mult = 1.0
    SIZE_LOOP : do m=1,nsize_aer(itype)
       if( ( nch .eq. lptr_so4_aer(m,itype,iphase) ) .or.  &
            (nch .eq. lptr_no3_aer(m,itype,iphase) ) .or.  &
            (nch .eq. lptr_nh4_aer(m,itype,iphase) )  )then
! so4, no3 and nh4 aerosol profiles
          if ( zz <= 1000. ) then
             mult = 1.0
		  elseif( zz > 1000. &
                 .and. zz <= 4000. ) then
             mult = 1.0 - 2.367e-4*(zz-1000.)
          elseif( zz > 4000. &
                 .and. zz <= 5500. ) then
             mult = 0.29 - 4.0e-5*(zz-4000.)
          else
             mult = 0.23
          end if
          exit SIZE_LOOP
       else if ( (nch .eq. lptr_na_aer(m,itype,iphase) ) .or.   &
	             (nch .eq. lptr_cl_aer(m,itype,iphase) ) ) then
! na and cl aerosol profiles
          if ( zz <= 100. ) then
             mult = 1.0
		  elseif( zz > 100. &
                  .and. zz <= 265. ) then
             mult = 1.0 - 2.9e-3*(zz-100.)
          elseif( zz > 265. &
                  .and. zz <= 2000. ) then
             mult = 0.52 - 2.94e-4*(zz-265.)
          elseif( zz > 2000. &
                  .and. zz <= 7000. ) then
             mult = 0.01
          else
             mult = 1.e-10
          end if
          exit SIZE_LOOP
       else if (nch .eq. lptr_oc_aer(m,itype,iphase) )  then
! oc aerosol profile
          if ( zz <= 600. ) then
             mult = 1.0
		  elseif( zz > 600. &
                  .and. zz <= 3389. ) then
             mult = 1.0 - 2.04e-4*(zz-600.)
          elseif( zz > 3389. &
                  .and. zz <= 8840. ) then
             mult = 0.43 - 7.045e-5*(zz-3389.)
          else
             mult = 0.046
          end if
          exit SIZE_LOOP
       else if  (nch .eq. lptr_bc_aer(m,itype,iphase) )   then
! bc aerosol profile
          if ( zz <= 100. ) then
             mult = 1.0
		  elseif( zz > 100. &
                  .and. zz <= 3400. ) then
             mult = 1.0 - 2.51e-4*(zz-100.)
          elseif( zz > 3400. &
                  .and. zz <= 8400. ) then
             mult = 0.172 -2.64e-5*(zz-3400.)
          else
             mult = 0.04
          end if
          exit SIZE_LOOP
       else if  (nch .eq. lptr_oin_aer(m,itype,iphase))  then
! oin aerosol profile
          if ( zz <= 1580. ) then
             mult = 1.0 + 1.77e-4 *zz
		  elseif( zz > 1580. &
                  .and. zz <= 5280. ) then
             mult = 1.28 - 2.65e-4*(zz-1580.)
          else
             mult = 0.30
          end if
          exit SIZE_LOOP
       else
! generic profile
#endif
! Updated aerosol profile multiplier 1-Apr-2005:
!    Height(m)     Multiplier
!    ---------     ----------
!    <=2000        1.0
!    2000<z<3000   linear transition zone to 0.25
!    3000<z<5000   linear transision zone to 0.125
!    >=5000        0.125
!
! which translates to:
!    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
!    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
! or in reduced form:
        if( zz <= 2000. ) then
           mult = 1.0
        elseif( zz > 2000. &
             .and. zz <= 3000. ) then
           mult = 1.0 - 0.00075*(zz-2000.)
        elseif( zz > 3000. &
             .and. zz <= 5000. ) then
           mult = 0.25 - 4.166666667e-5*(zz-3000.)
        else
           mult = 0.125
        end if
#if (CASENAME == 1)
      end if			! end nc block comparison
    end do SIZE_LOOP
#endif
#if (CASENAME == 2)
!     if( zz <= 3000. ) then
!       mult = 1.0
!     elseif( zz > 3000. &
!       .and. zz <= 5000. ) then
!       mult = 1.0 - 0.0004375*(zz-3000.)
!     else
!        mult = 0.125
!     end if
!     if( zz <= 5000. ) then
!       mult = 1.0
!     elseif( zz > 5000. &
!       .and. zz <= 7000. ) then
!       mult = 1.0 - 0.0004375*(zz-5000.)
!     else
!        mult = 0.125
!     end if
      mult=1.0
      iflg=0
      do 2901 m = 1, nsize_aer(itype)
      if(iflg.eq.1) goto 2902
      if( ( nch .eq. lptr_so4_aer(m,itype,iphase) ) .or.  &
        (nch .eq. lptr_no3_aer(m,itype,iphase) ) .or.  &
        (nch .eq. lptr_nh4_aer(m,itype,iphase) )  )then
! so4, no3 and nh4 aerosol profiles
        if( zz <= 500. ) then
          mult = 1.000
        elseif( zz > 500. &
          .and. zz <= 1000. ) then
          mult = 1.000 - 0.0010740*(zz-500.)
        elseif( zz > 1000. &
          .and. zz <= 5000. ) then
          mult = 0.463 - 0.0001110*(zz-1000.)
        elseif( zz > 5000. &
          .and. zz <= 20000. ) then
          mult = 0.019
        else
          mult = 0.019
        end if
        iflg=1
      else if ( (nch .eq. lptr_na_aer(m,itype,iphase) ) .or.   &
                (nch .eq. lptr_cl_aer(m,itype,iphase) ) ) then
! na and cl aerosol profiles
        if( zz <= 500. ) then
          mult = 1.000
        elseif( zz > 500. &
          .and. zz <= 1000. ) then
          mult = 1.000 - 0.0014260*(zz-500.)
        elseif( zz > 1000. &
          .and. zz <= 5000. ) then
          mult = 0.287 - 0.0000643*(zz-1000.)
        elseif( zz > 5000. &
          .and. zz <= 10000. ) then
          mult = 0.030
        elseif( zz > 10000. &
          .and. zz <= 20000. ) then
          mult = 0.030 - 0.0000028*(zz-10000.)
        else
          mult = 0.002
        end if
        iflg=1
      else if ( (nch .eq. lptr_oc_aer(m,itype,iphase) ) .or.   &
                (nch .eq. lptr_bc_aer(m,itype,iphase) ) ) then
! oc and bc aerosol profile
        if( zz <= 500. ) then
          mult = 1.000
        elseif( zz > 500. &
          .and. zz <= 1000. ) then
          mult = 1.000 - 0.0002500*(zz-500.)
        elseif( zz > 1000. &
          .and. zz <= 5000. ) then
          mult = 0.875 - 0.0001843*(zz-1000.)
        elseif( zz > 5000. &
          .and. zz <= 20000. ) then
          mult = 0.138
        else
          mult = 0.138
        end if
        iflg=1
      else if  (nch .eq. lptr_oin_aer(m,itype,iphase))  then
! oin aerosol profile
        if( zz <= 500. ) then
          mult = 1.000
        elseif( zz > 500. &
          .and. zz <= 1000. ) then
          mult = 1.000 - 0.0007340*(zz-500.)
        elseif( zz > 1000. &
          .and. zz <= 5000. ) then
          mult = 0.633 - 0.0001103*(zz-1000.)
        elseif( zz > 5000. &
          .and. zz <= 10000. ) then
          mult = 0.192 - 0.0000080*(zz-5000.)
        elseif( zz > 10000. &
          .and. zz <= 20000. ) then
          mult = 0.152 - 0.0000137*(zz-10000.)
        else
          mult = 0.015
        end if
        iflg=1
      else
        if( zz <= 5000. ) then
          mult = 1.0
        elseif( zz > 5000. &
          .and. zz <= 7000. ) then
          mult = 1.0 - 0.0004375*(zz-5000.)
        else
           mult = 0.125
        end if
      end if
2901  continue
2902  continue
#endif
#if (CASENAME == 3)
      if( zz <= 500. ) then
        mult = 1.000
      elseif( zz > 500. &
        .and. zz <= 1000. ) then
        mult = 1.000 - 0.0010740*(zz-500.)
      elseif( zz > 1000. &
        .and. zz <= 5000. ) then
        mult = 0.463 - 0.0001110*(zz-1000.)
      else
        mult = 0.019
      end if
#endif
!   Set values for fine and coarse mass concentrations, and
!   compute fine/coarse volume concentrations. These are also set in
!   mosaic_init_wrf_mixrats_opt2.
!   (these are latest values from 1-Apr-2005 discussions)
!wig 10-May-2004, added height multiplier (mult*) to each species...
	qfine_so4 = mult*2.14
	qcoar_so4 = mult*0.242
	qfine_no3 = mult*0.11
	qcoar_no3 = mult*0.03
!	qfine_cl  = mult*0.3
	qfine_cl  = mult*0.14     ! 10-May-2005
	qcoar_cl  = mult*0.139
	qfine_msa = mult*0.0
	qcoar_msa = mult*0.0
	qfine_co3 = mult*0.0
	qcoar_co3 = mult*0.0
	qfine_nh4 = mult*0.83
	qcoar_nh4 = mult*0.10
!	qfine_na  = mult*0.2
	qfine_na  = mult*0.1      ! 10-May-2005
	qcoar_na  = mult*0.09
	qfine_ca  = mult*0.0
	qcoar_ca  = mult*0.0
!	qfine_oin = mult*6.2
!	qfine_oin = 3.48     ! 10-May-2005
	qfine_oin = mult*3.48     ! not sure why mult was dropped from qfine_oin
	qcoar_oin = mult*0.35
	qfine_oc  = mult*1.00
	qcoar_oc  = mult*1.50
	qfine_bc  = mult*0.2
	qcoar_bc  = mult*0.075
	qfine_hysw = mult*0.0
	qcoar_hysw = mult*0.0
	qfine_watr = mult*0.0
	qcoar_watr = mult*0.0
#if (CASENAME == 2)
!       qfine_so4 = mult*0.50
!       qcoar_so4 = mult*0.50/2.0
!       qfine_no3 = mult*0.75
!       qcoar_no3 = mult*0.75/2.0
!       qfine_cl  = mult*0.10
!       qcoar_cl  = mult*0.10/2.0
!       qfine_nh4 = mult*0.30
!       qcoar_nh4 = mult*0.30/2.0
!       qfine_na  = mult*0.20
!       qcoar_na  = mult*0.20/2.0
!       qfine_oin = mult*4.00
!       qcoar_oin = mult*4.00/2.0
!       qfine_oc  = mult*3.00
!       qcoar_oc  = mult*3.00/2.0
!       qfine_bc  = mult*0.50
!       qcoar_bc  = mult*0.50/2.0
        qfine_so4 = mult*0.300
        qcoar_so4 = mult*0.300/2.0
        qfine_no3 = mult*0.001
        qcoar_no3 = mult*0.001/2.0
        qfine_cl  = mult*1.050
        qcoar_cl  = mult*1.05/02.0
        qfine_nh4 = mult*0.094
        qcoar_nh4 = mult*0.094/2.0
        qfine_na  = mult*0.700
        qcoar_na  = mult*0.700/2.0
        qfine_oin = mult*4.500
        qcoar_oin = mult*4.500/2.0
        qfine_oc  = mult*0.088
        qcoar_oc  = mult*0.088/2.0
        qfine_bc  = mult*0.013
        qcoar_bc  = mult*0.013/2.0
#endif
#if (CASENAME == 3)
        qfine_so4 = mult*0.300
        qcoar_so4 = mult*0.0
        qfine_nh4 = mult*0.094
        qcoar_nh4 = mult*0.0
        qfine_no3 = mult*0.001
        qcoar_no3 = mult*0.0
        qfine_cl  = mult*0.0
        qcoar_cl  = mult*1.050
        qfine_na  = mult*0.000
        qcoar_na  = mult*0.700
        qfine_oin = mult*4.500
        qcoar_oin = mult*4.500/2.0
        qfine_oc  = mult*0.088
        qcoar_oc  = mult*0.0
        qfine_bc  = mult*0.013
        qcoar_bc  = mult*0.0
#endif
!!$! old values from 23-Apr-2004:
!!$	qfine_so4 = mult*2.554
!!$	qcoar_so4 = mult*0.242
!!$	qfine_no3 = mult*0.07
!!$	qcoar_no3 = mult*0.03
!!$	qfine_cl  = mult*0.324
!!$	qcoar_cl  = mult*0.139
!!$	qfine_msa = mult*0.0
!!$	qcoar_msa = mult*0.0
!!$	qfine_co3 = mult*0.0
!!$	qcoar_co3 = mult*0.0
!!$	qfine_nh4 = mult*0.98
!!$	qcoar_nh4 = mult*0.10
!!$	qfine_na  = mult*0.21
!!$	qcoar_na  = mult*0.09
!!$	qfine_ca  = mult*0.0
!!$	qcoar_ca  = mult*0.0
!!$	qfine_oin = mult*0.15
!!$	qcoar_oin = mult*0.35
!!$	qfine_oc  = mult*1.00
!!$	qcoar_oc  = mult*1.50
!!$	qfine_bc  = mult*0.175
!!$	qcoar_bc  = mult*0.075
!!$	qfine_hysw = mult*0.0
!!$	qcoar_hysw = mult*0.0
!!$	qfine_watr = mult*0.0
!!$	qcoar_watr = mult*0.0


! qfine_so4 ... are ug/m3 so 1.0e-6 factor gives g/m3
! dens_so4  ... are g/cm3;  qfine_vol ... are cm3/m3
	qfine_vol = 1.0e-6 * (   &
	    (qfine_so4/dens_so4_aer) + (qfine_no3/dens_no3_aer) +   &
	    (qfine_cl /dens_cl_aer ) + (qfine_msa/dens_msa_aer) +   &
	    (qfine_co3/dens_co3_aer) + (qfine_nh4/dens_nh4_aer) +   &
	    (qfine_na /dens_na_aer ) + (qfine_ca /dens_ca_aer ) +   &
	    (qfine_oin/dens_oin_aer) + (qfine_oc /dens_oc_aer ) +   &
	    (qfine_bc /dens_bc_aer ) )
	qcoar_vol = 1.0e-6 * (   &
	    (qcoar_so4/dens_so4_aer) + (qcoar_no3/dens_no3_aer) +   &
	    (qcoar_cl /dens_cl_aer ) + (qcoar_msa/dens_msa_aer) +   &
	    (qcoar_co3/dens_co3_aer) + (qcoar_nh4/dens_nh4_aer) +   &
	    (qcoar_na /dens_na_aer ) + (qcoar_ca /dens_ca_aer ) +   &
	    (qcoar_oin/dens_oin_aer) + (qcoar_oc /dens_oc_aer ) +   &
	    (qcoar_bc /dens_bc_aer ) )

	qfine = -1.0e30
	qcoar = -1.0e30

!   identify the species by checking "nch" against the "lptr_xxx_a_amode(m)"
        do 2900 m = 1, nsize_aer(itype)

	    if (nch .eq. lptr_so4_aer(m,itype,iphase)) then
		qfine = qfine_so4
		qcoar = qcoar_so4
	    else if (nch .eq. lptr_no3_aer(m,itype,iphase)) then
		qfine = qfine_no3
		qcoar = qcoar_no3
	    else if (nch .eq. lptr_cl_aer(m,itype,iphase)) then
		qfine = qfine_cl
		qcoar = qcoar_cl
	    else if (nch .eq. lptr_msa_aer(m,itype,iphase)) then
		qfine = qfine_msa
		qcoar = qcoar_msa
	    else if (nch .eq. lptr_co3_aer(m,itype,iphase)) then
		qfine = qfine_co3
		qcoar = qcoar_co3
	    else if (nch .eq. lptr_nh4_aer(m,itype,iphase)) then
		qfine = qfine_nh4
		qcoar = qcoar_nh4
	    else if (nch .eq. lptr_na_aer(m,itype,iphase)) then
		qfine = qfine_na
		qcoar = qcoar_na
	    else if (nch .eq. lptr_ca_aer(m,itype,iphase)) then
		qfine = qfine_ca
		qcoar = qcoar_ca
	    else if (nch .eq. lptr_oin_aer(m,itype,iphase)) then
		qfine = qfine_oin
		qcoar = qcoar_oin
	    else if (nch .eq. lptr_oc_aer(m,itype,iphase)) then
		qfine = qfine_oc
		qcoar = qcoar_oc
	    else if (nch .eq. lptr_bc_aer(m,itype,iphase)) then
		qfine = qfine_bc
		qcoar = qcoar_bc
	    else if (nch .eq. hyswptr_aer(m,itype)) then
		qfine = qfine_hysw
		qcoar = qcoar_hysw
	    else if (nch .eq. waterptr_aer(m,itype)) then
		qfine = qfine_watr
		qcoar = qcoar_watr
	    else if (nch .eq. numptr_aer(m,itype,iphase)) then
		dumvol1p = sqrt(volumlo_sect(m,itype)*volumhi_sect(m,itype))
		qfine = qfine_vol/dumvol1p
		qcoar = qcoar_vol/dumvol1p
	    end if

	    if ((qfine >= 0.0) .and. (qcoar >= 0.0)) then
		qval = qfine*fr_fine(m) + qcoar*fr_coar(m)
		chem_bv = qval*alt
		goto 2910
	    end if

2900   continue
2910   continue

	return
 	end subroutine bdy_chem_value_mosaic