!********************************************************************************** ! 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=5000 0.25 !!$! !!$! which translates to: !!$! 2000 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=5000 0.125 ! ! which translates to: ! 2000 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=5000 0.25 !!$! !!$! which translates to: !!$! 2000 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=5000 0.125 ! ! which translates to: ! 2000 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