!********************************************************************************** ! 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 !********************************************************************************** MODULE module_mosaic_addemiss !WRF:MODEL_LAYER:CHEMICS ! rce 2005-feb-18 - one fix for indices of volumcen_sect, [now (isize,itype)] ! rce 2005-jan-14 - added subr mosaic_seasalt_emiss (and a call to it) ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer" ! variables in module_data_mosaic_asect integer, parameter :: mosaic_addemiss_active = 1 ! only do emissions when this is positive ! (when it is negative, emissions tendencies are zero) integer, parameter :: mosaic_addemiss_masscheck = -1 ! only do emissions masscheck calcs when this is positive CONTAINS !---------------------------------------------------------------------- subroutine mosaic_addemiss( id, dtstep, u10, v10, alt, dz8w, xland, & config_flags, chem, slai, ust, smois, ivgtyp, isltyp, & emis_ant,ebu,biom_active,dust_opt, & !czhao ktau,p8w, u_phy,v_phy,rho_phy,g,dx,erod, & ! GOCART DUST dust_emiss_active, seasalt_emiss_active, & dust_flux, seas_flux, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ! ! adds emissions for mosaic aerosol species ! (i.e., emissions tendencies over time dtstep are applied ! to the aerosol concentrations) ! USE module_configure, only: grid_config_rec_type USE module_state_description ! USE module_state_description, only: num_chem, param_first_scalar, & ! emiss_inpt_default, emiss_inpt_pnnl_rs, emiss_inpt_pnnl_cm,num_emis_ant USE module_data_mosaic_asect IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: id, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN) :: dust_emiss_active, seasalt_emiss_active, biom_active, dust_opt REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: dust_flux, seas_flux !czhao INTEGER, INTENT(IN) :: ktau REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN) :: p8w, u_phy,v_phy,rho_phy REAL, INTENT(IN ) :: dx, g REAL, DIMENSION( ims:ime, jms:jme,3),& INTENT(IN ) :: erod REAL, INTENT(IN ) :: dtstep ! 10-m wind speed components (m/s) REAL, DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) :: u10, v10, xland, slai, ust INTEGER, DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) :: ivgtyp, isltyp ! trace species mixing ratios (aerosol mass = ug/kg-air; number = #/kg-air) REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem ! ! aerosol emissions arrays ((ug/m3)*m/s) ! ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & REAL, DIMENSION( ims:ime, 1:config_flags%kemit, jms:jme,num_emis_ant ), & INTENT(IN ) :: & emis_ant REAL, DIMENSION( ims:ime, kms:kme, jms:jme,num_ebu ), & INTENT(IN ) :: & ebu ! 1/(dry air density) and layer thickness (m) REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: alt, dz8w REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , & INTENT(INOUT) :: smois ! local variables integer i, j, k, l, n integer iphase, itype integer p1st real, parameter :: efact1 = 1.0 real :: aem_so4, aem_no3, aem_cl, aem_msa, aem_co3, aem_nh4, & aem_na, aem_ca, aem_oin, aem_oc, aem_bc, aem_num real dum, fact ! fraction of sorgam i/aitken mode emissions that go to each ! of the mosaic 8 "standard" sections !now the anthropogenic and biomass burning emissiosn use same distribution !it can be improve later --czhao ! not 100%, because most mass falling down to the size less than the lowest boundary real, save :: fr8b_aem_sorgam_i(8) = & (/ 0.965, 0.035, 0.000, 0.000, & 0.000, 0.000, 0.000, 0.000 /) ! fraction of sorgam j/accum mode emissions that go to each ! of the mosaic 8 "standard" sections real, save :: fr8b_aem_sorgam_j(8) = & (/ 0.026, 0.147, 0.350, 0.332, & 0.125, 0.019, 0.001, 0.000/) ! fraction of sorgam coarse mode emissions that go to each ! of the mosaic 8 "standard" sections ! not 100%, because most mass falling up to the size larger than the highest boundary real, save :: fr8b_aem_sorgam_c(8) = & (/ 0.000, 0.000, 0.000, 0.002, & 0.021, 0.110, 0.275, 0.592 /) ! fraction of mosaic fine (< 2.5 um) emissions that go to each ! of the mosaic 8 "standard" sections !wig 1-Apr-2005, Updated fractional breakdown between bins. (See also ! bdy_chem_value_mosaic and mosaic_init_wrf_mixrats_opt2 ! in module_mosaic_initmixrats.F.) Note that the values ! here no longer match the other two subroutines. !rce 10-may-2005, changed fr8b_aem_mosaic_f & fr8b_aem_mosaic_c ! to values determined by jdf real, save :: fr8b_aem_mosaic_f(8) = & (/ 0.060, 0.045, 0.245, 0.400, 0.100, 0.150, 0., 0./) !10-may-2005 ! (/ 0.100, 0.045, 0.230, 0.375, 0.100, 0.150, 0., 0./) !1-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.0, 0.0 /) ! fraction of mosaic coarse (> 2.5 um) emissions that go to each ! of the mosaic 8 "standard" sections real, save :: fr8b_aem_mosaic_c(8) = & (/ 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 /) ! as of apr-2005 ! (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.16, 0.84 /) ! "old" !++ SAN Fractionation for PREP-Chem Fire emissions ++! ! New fractions for biomass burning emissions. Based on Janhall et. al. (2010) ! Based on average fresh smoke emissions from 20 data points ! mean_Dp = .117um, geometric_std = 1.7 real, save :: fr8b_bburn_mosaic(8) = & (/ 0.0092, 0.1385, 0.4548, 0.3388, 0.0567, 0.002, 0.0, 0.0/) !--SAN ! fraction of TNO black carbon <1um emissions that go into each of ! the mosaic 8 "standard" sections - Doug (4/3/2011) real, save :: fr8b_tno_bc1(8) = & (/ 0.0494, 0.3795, 0.4714, 0.0967, 0.003, 0.0, 0.0, 0.0 /) ! fraction of TNO elemental carbon 1um-2.5um emissions that go into each of ! the mosaic 8 "standard" sections - Doug (4/3/2011) real, save :: fr8b_tno_ec25(8) = & (/ 0.0, 0.0, 0.0, 0.0, 0.40, 0.60, 0.0, 0.0 /) ! fraction of TNO organic carbon <2.5um domestic combustion emissions that go into each of ! the mosaic 8 "standard" sections - Doug (4/3/2011) real, save :: fr8b_tno_oc_dom(8) = & (/ 0.0358, 0.1325, 0.2704, 0.3502, 0.1904, 0.0657, 0.0, 0.0 /) ! fraction of TNO organic carbon <2.5um traffic (and other) emissions that go into each of ! the mosaic 8 "standard" sections - Doug (4/3/2011) real, save :: fr8b_tno_oc_tra(8) = & (/ 0.0063, 0.0877, 0.3496, 0.4054, 0.1376, 0.0134, 0.0, 0.0 /) ! fraction of TNO "OIN" [PM2.5 minus sum(carbon<2.5)] emissions that go into each of ! the mosaic 8 "standard" sections - Doug (4/3/2011) --- based on mosaic fine mode real, save :: fr8b_tno_mosaic_f(8) = & (/ 0.060, 0.045, 0.245, 0.400, 0.100, 0.150, 0.0, 0.0/) ! fraction of TNO 2.5-10um emissions that go into each of ! the mosaic 8 "standard" sections - Doug (4/3/2011) --- based on mosaic coarse mode real, save :: fr8b_tno_mosaic_c(8) = & (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.300, 0.700 /) ! following 5 arrays correspond to the above "fr8b_" arrays ! but are set at run time base on input (namelist) parameters: ! only the sorgam or mosaic arrays are non-zero, depending on ! emiss_inpt_opt ! when nsize_aer=4 (=number of sections), the values are ! calculated by adding two of the 8-section values real :: fr_aem_sorgam_i(8) real :: fr_aem_sorgam_j(8) real :: fr_aem_sorgam_c(8) real :: fr_aem_mosaic_f(8) real :: fr_aem_mosaic_c(8) real :: fr_tno_bc1(8) real :: fr_tno_ec25(8) real :: fr_tno_oc_dom(8) real :: fr_tno_oc_tra(8) real :: fr_tno_mosaic_f(8) real :: fr_tno_mosaic_c(8) !++ SAN real :: fr_aem_gc2mosaic_f(8) ! extra arrays for GOCART->MOSAIC mapping real :: fr_aem_gc2mosaic_c(8) real :: bburn_mosaic_f(8) ! Arrays to hold biomass-burning size real :: bburn_mosaic_c(8) ! Arrays to hold biomass-burning size !-- SAN double precision :: chem_sum(num_chem) character(len=80) :: msg ! *** currently only works with ntype_aer = 1 itype = 1 iphase = ai_phase ! if (num_ebu.le.0 ) then ! call wrf_error_fatal( 'mosaic_addemiss: no biomass burning species' ) ! print*,'no biomass burning species',num_ebu ! stop ! endif ! compute factors used for apportioning either ! the MADE-SORGAM emissions (i=aitken, j=accum, c=coarse modes) OR ! the MOSAIC emission (f=fine (< 2.5 um), c=coarse (> 2.5 um)) ! to each size section ! ! note: the fr8b_aer_xxxxxx_y values are specific to the mosaic 8 bin ! structure with dlo_sect(1)=0.039 um and dhi_sect(8)=10.0 um, ! also, the fr8b_aem_sorgam_y are specific for the assumed ! dgvem_i/j/c used in the module_data_sorgam.F code ! also, the fr8b_aem_mosaic_y values are specific for the assumed (by us) ! size distribution for fine and coarse primary emissions ! ! when there are 4 bins (nsize_aer=4), each of these "wider" bins ! corresponds to 2 of the "narrower" bins of the 8 bin structure ! ! note: if fr_aem_sorgam_y > 0, then fr_aem_mosaic_y = 0, and vice-versa ! if ((nsize_aer(itype) .ne. 4) .and. (nsize_aer(itype) .ne. 8)) then write(msg,'(a,i5)') & 'subr mosaic_addemiss - nsize_aer(itype) must be ' // & '4 or 8 but = ', nsize_aer(itype) call wrf_error_fatal( msg ) end if fr_aem_sorgam_i(:) = 0.0 fr_aem_sorgam_j(:) = 0.0 fr_aem_sorgam_c(:) = 0.0 fr_aem_mosaic_f(:) = 0.0 fr_aem_mosaic_c(:) = 0.0 fr_tno_bc1(:) = 0.0 fr_tno_ec25(:) = 0.0 fr_tno_oc_dom(:) = 0.0 fr_tno_oc_tra(:) = 0.0 fr_tno_mosaic_f(:) = 0.0 fr_tno_mosaic_c(:) = 0.0 fr_aem_gc2mosaic_f(:) = 0.0 fr_aem_gc2mosaic_c(:) = 0.0 emiss_inpt_select_1: SELECT CASE( config_flags%emiss_inpt_opt ) CASE( emiss_inpt_default, emiss_inpt_pnnl_rs ) if (nsize_aer(itype) .eq. 8) then fr_aem_sorgam_i(:) = fr8b_aem_sorgam_i(:) fr_aem_sorgam_j(:) = fr8b_aem_sorgam_j(:) fr_aem_sorgam_c(:) = fr8b_aem_sorgam_c(:) else if (nsize_aer(itype) .eq. 4) then do n = 1, nsize_aer(itype) fr_aem_sorgam_i(n) = fr8b_aem_sorgam_i(2*n-1) & + fr8b_aem_sorgam_i(2*n) fr_aem_sorgam_j(n) = fr8b_aem_sorgam_j(2*n-1) & + fr8b_aem_sorgam_j(2*n) fr_aem_sorgam_c(n) = fr8b_aem_sorgam_c(2*n-1) & + fr8b_aem_sorgam_c(2*n) end do end if !++ SAN, 2015-04-09 ! Added mapping routine from GOCART aerosol emission variables to MOSAIC ! This will be the case if emissions were prepared using PREP-Chem in RAMD2-GOCART mode. ! (Use with emiss_opt = 5/6, ECPTEC/GOCART_ECPTEC, emiss_inpt_opt == emiss_inpt_default) !!! Maybe more consistant to add a new emiss_inpt_opt, but this requires more extensive changes !!! and testing in other parts of WRF-Chem, may cause confusion... Thoughts? if( config_flags%emiss_opt == 5 .or. config_flags%emiss_opt == 6 ) then CALL wrf_debug(15,'mosaic_addemiss: emiss_opt = eptec') CALL wrf_debug(15,'mosaic_addemiss: gocart speciation being mapped to mosaic') fr_aem_sorgam_i(:) = 0.0 fr_aem_sorgam_j(:) = 0.0 fr_aem_sorgam_c(:) = 0.0 ! Use FM MOSAIC mapping for SO4, OC, BC and PM2.5, CM mosaic for PM10 if (nsize_aer(itype) .eq. 8) then fr_aem_gc2mosaic_f(:) = fr8b_aem_mosaic_f(:) fr_aem_gc2mosaic_c(:) = fr8b_aem_mosaic_c(:) elseif (nsize_aer(itype) .eq. 4) then do n = 1, nsize_aer(itype) fr_aem_gc2mosaic_f(n) = fr8b_aem_mosaic_f(2*n-1) + fr8b_aem_mosaic_f(2*n) fr_aem_gc2mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) + fr8b_aem_mosaic_c(2*n) end do endif endif !-- SAN CASE( emiss_inpt_pnnl_cm ) if (nsize_aer(itype) .eq. 8) then fr_aem_mosaic_f(:) = fr8b_aem_mosaic_f(:) fr_aem_mosaic_c(:) = fr8b_aem_mosaic_c(:) else if (nsize_aer(itype) .eq. 4) then do n = 1, nsize_aer(itype) fr_aem_mosaic_f(n) = fr8b_aem_mosaic_f(2*n-1) & + fr8b_aem_mosaic_f(2*n) fr_aem_mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) & + fr8b_aem_mosaic_c(2*n) end do end if CASE( emiss_inpt_tno ) ! Doug - added 4/3/2011 if (nsize_aer(itype) .eq. 8) then fr_tno_bc1(:) = fr8b_tno_bc1(:) fr_tno_ec25(:) = fr8b_tno_ec25(:) fr_tno_oc_dom(:) = fr8b_tno_oc_dom(:) fr_tno_oc_tra(:) = fr8b_tno_oc_tra(:) fr_tno_mosaic_c(:) = fr8b_tno_mosaic_c(:) fr_tno_mosaic_f(:) = fr8b_tno_mosaic_f(:) else if (nsize_aer(itype) .eq. 4) then do n = 1, nsize_aer(itype) fr_tno_bc1(n) = fr8b_tno_bc1(2*n-1) & + fr8b_tno_bc1(2*n) fr_tno_ec25(n) = fr8b_tno_ec25(2*n-1) & + fr8b_tno_ec25(2*n) fr_tno_oc_dom(n) = fr8b_tno_oc_dom(2*n-1) & + fr8b_tno_oc_dom(2*n) fr_tno_oc_tra(n) = fr8b_tno_oc_tra(2*n-1) & + fr8b_tno_oc_tra(2*n) fr_tno_mosaic_c(n) = fr8b_tno_mosaic_c(2*n-1) & + fr8b_tno_mosaic_c(2*n) fr_tno_mosaic_f(n) = fr8b_tno_mosaic_f(2*n-1) & + fr8b_tno_mosaic_f(2*n) end do end if CASE DEFAULT return END SELECT emiss_inpt_select_1 !++ SAN 2015-04-09: Mapping fire emissions to MOSAIC ! get arrays for apportioning mass into MOSAIC bins for fire emissions fire_inpt_select: SELECT CASE (config_flags%biomass_burn_opt) CASE (BIOMASSB,BIOMASSB_MOZC) if (nsize_aer(itype) .eq. 8) then bburn_mosaic_f(:) = fr8b_bburn_mosaic(:) !++SAN bburn_mosaic_c(:) = fr8b_aem_mosaic_c(:) !--SAN else if (nsize_aer(itype) .eq. 4) then do n = 1, nsize_aer(itype) bburn_mosaic_f(n) = fr8b_bburn_mosaic(2*n-1) + fr8b_bburn_mosaic(2*n) !++SAN bburn_mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) + fr8b_aem_mosaic_c(2*n) !--SAN end do end if END SELECT fire_inpt_select ! when mosaic_addemiss_active <= 0, set fr's to zero, ! which causes the changes to chem(...) to be zero if (mosaic_addemiss_active <= 0 .and. biom_active <= 0) then fr_aem_sorgam_i(:) = 0.0 fr_aem_sorgam_j(:) = 0.0 fr_aem_sorgam_c(:) = 0.0 fr_aem_mosaic_f(:) = 0.0 fr_aem_mosaic_c(:) = 0.0 fr_tno_bc1(:) = 0.0 fr_tno_ec25(:) = 0.0 fr_tno_oc_dom(:) = 0.0 fr_tno_oc_tra(:) = 0.0 fr_tno_mosaic_f(:) = 0.0 fr_tno_mosaic_c(:) = 0.0 end if ! do mass check initial calc if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( & id, config_flags, 1, 'mosaic_ademiss', & dtstep, efact1, dz8w, chem, chem_sum, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & 14, & emis_ant(ims, 1,jms,p_e_pm_10),emis_ant(ims, 1,jms,p_e_pm_25), & emis_ant(ims, 1,jms,p_e_pm25i),emis_ant(ims, 1,jms,p_e_pm25j), & emis_ant(ims, 1,jms,p_e_eci),emis_ant(ims, 1,jms,p_e_ecj), & emis_ant(ims, 1,jms,p_e_orgi),emis_ant(ims, 1,jms,p_e_orgj), & emis_ant(ims, 1,jms,p_e_so4j),emis_ant(ims, 1,jms,p_e_so4c), & emis_ant(ims, 1,jms,p_e_no3j),emis_ant(ims, 1,jms,p_e_no3c), & emis_ant(ims, 1,jms,p_e_orgc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc)) p1st = param_first_scalar ! ! apply emissions at each section and grid point ! do 1900 n = 1, nsize_aer(itype) do 1830 j = jts, jte do 1820 k = 1, min(config_flags%kemit,kte) do 1810 i = its, ite ! compute anthropogenic mass emissions [(ug/m3)*m/s] for each species ! using the apportioning fractions aem_so4 = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_so4j) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_so4c) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_so4i) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_so4j) aem_no3 = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_no3j) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_no3c) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_no3i) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_no3j) aem_oc = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_orgj) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_orgc) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_orgi) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_orgj) & + fr_tno_oc_dom(n)*emis_ant(i,k,j,p_e_oc_dom) & + fr_tno_oc_tra(n)*emis_ant(i,k,j,p_e_oc_tra) & + fr_tno_mosaic_c(n)*emis_ant(i,k,j,p_e_oc_25_10) & !++SAN + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_oc) !--SAN chem_select_1 : SELECT CASE( config_flags%chem_opt ) CASE(SAPRC99_MOSAIC_4BIN_VBS2_KPP, SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP, & SAPRC99_MOSAIC_8BIN_VBS2_KPP)!BSINGH(12/04/2013): Added SAPRC 8 bin and non-aq on 04/03/2014! Set the oc to zero for VBS aem_oc = 0.0 END SELECT chem_select_1 aem_bc = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_ecj) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_ecc) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_eci) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_ecj) & + fr_tno_bc1(n)*emis_ant(i,k,j,p_e_bc_1) & + fr_tno_ec25(n)*emis_ant(i,k,j,p_e_ec_1_25) & + fr_tno_mosaic_c(n)*emis_ant(i,k,j,p_e_ec_25_10) & !++SAN + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_bc) !--SAN aem_oin = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_pm25j) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_pm_10) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_pm25i) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_pm25j) & !++SAN + fr_aem_sorgam_c(n)*emis_ant(i,k,j,p_e_pm_10) & !--SAN + fr_tno_mosaic_f(n)*emis_ant(i,k,j,p_e_oin_25) & + fr_tno_mosaic_c(n)*emis_ant(i,k,j,p_e_oin_10) & !++SAN + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_pm25) & + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_pm10) !--SAN ! emissions for these species are currently zero aem_nh4 = 0.0 aem_na = 0.0 aem_cl = 0.0 aem_ca = 0.0 aem_co3 = 0.0 aem_msa = 0.0 chem_select_2 : SELECT CASE( config_flags%chem_opt ) CASE(MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP) aem_nh4 = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_nh4j) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_nh4c) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_nh4i) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_nh4j) aem_na = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_naj) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_nac) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_nai) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_naj) aem_cl = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_clj) & + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_clc) & + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_cli) & + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_clj) END SELECT chem_select_2 ! compute number emissions ! first sum the mass-emissions/density aem_num = & (aem_so4/dens_so4_aer) + (aem_no3/dens_no3_aer) + & (aem_cl /dens_cl_aer ) + (aem_msa/dens_msa_aer) + & (aem_co3/dens_co3_aer) + (aem_nh4/dens_nh4_aer) + & (aem_na /dens_na_aer ) + (aem_ca /dens_ca_aer ) + & (aem_oin/dens_oin_aer) + (aem_oc /dens_oc_aer ) + & (aem_bc /dens_bc_aer ) ! then multiply by 1.0e-6 to convert ug to g ! and divide by particle volume at center of section (cm3) aem_num = aem_num*1.0e-6/volumcen_sect(n,itype) ! apply the emissions and convert from flux to mixing ratio ! fact = (dtstep/dz8w(i,k,j))*(28.966/1000.) fact = (dtstep/dz8w(i,k,j))*alt(i,k,j) ! rce 22-nov-2004 - change to using the "..._aer" species pointers l = lptr_so4_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_so4*fact l = lptr_no3_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_no3*fact l = lptr_cl_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_cl*fact l = lptr_msa_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_msa*fact l = lptr_co3_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_co3*fact l = lptr_nh4_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_nh4*fact l = lptr_na_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_na*fact l = lptr_ca_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_ca*fact l = lptr_oin_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oin*fact l = lptr_oc_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oc*fact l = lptr_bc_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_bc*fact l = numptr_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_num*fact 1810 continue 1820 continue 1830 continue 1900 continue if (num_ebu <= 0 ) then return endif !++ SAN, 2015-04-09. ! ---------------------- ADD BIOMASS BURNING EMISSIONS ----------------------------! ! - new case select for adding biomass burning emissions. ! First reset all emissions to make sure we don't double count anthropogenic emissions aem_so4 = 0. aem_no3 = 0. aem_oc = 0. aem_bc = 0. aem_oin = 0. aem_num = 0. fire_select: SELECT CASE(config_flags%biomass_burn_opt) CASE (BIOMASSB,BIOMASSB_MOZC) CALL wrf_debug(15,'mosaic_addemiss: adding fire emissions to MOSAIC') ! CALL wrf_debug(15,'e_oin_fm = PM2.5 - OC - BC') ! CALL wrf_debug(15,'e_oin_cm = PM10 - PM2.5') size_loop: & do n = 1, nsize_aer(itype) do j = jts, jte do k = kts, kte ! Loop up to kte for fire emissions (with plumerise) do i = its, ite ! compute mass biomass burning emissions [(ug/m3)*m/s] for each species ! using the apportioning fractions aem_so4 = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_sulf) aem_oc = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_oc) aem_bc = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_bc) ! Option to calculate OIN fraction of total PM for fire emissions ! Assume all OC and BC is in FM. aem_oin = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_pm25) & + bburn_mosaic_c(n)*ebu(i,k,j,p_ebu_pm10) chem_select_3 : SELECT CASE( config_flags%chem_opt ) CASE(SAPRC99_MOSAIC_4BIN_VBS2_KPP) ! Set the oc to zero for VBS aem_oc = 0.0 END SELECT chem_select_3 ! emissions for these species are currently zero aem_nh4 = 0.0 aem_na = 0.0 aem_cl = 0.0 aem_ca = 0.0 aem_co3 = 0.0 aem_msa = 0.0 ! compute number emissions ! first sum the mass-emissions/density aem_num = & (aem_so4/dens_so4_aer) + (aem_no3/dens_no3_aer) + & (aem_cl /dens_cl_aer ) + (aem_msa/dens_msa_aer) + & (aem_co3/dens_co3_aer) + (aem_nh4/dens_nh4_aer) + & (aem_na /dens_na_aer ) + (aem_ca /dens_ca_aer ) + & (aem_oin/dens_oin_aer) + (aem_oc /dens_oc_aer ) + & (aem_bc /dens_bc_aer ) ! then multiply by 1.0e-6 to convert ug to g ! and divide by particle volume at center of section (cm3) aem_num = aem_num*1.0e-6/volumcen_sect(n,itype) ! apply the emissions and convert from flux to mixing ratio ! fact = (dtstep/dz8w(i,k,j))*(28.966/1000.) fact = (dtstep/dz8w(i,k,j))*alt(i,k,j) ! rce 22-nov-2004 - change to using the "..._aer" species pointers l = lptr_so4_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_so4*fact l = lptr_no3_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_no3*fact l = lptr_cl_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_cl*fact l = lptr_msa_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_msa*fact l = lptr_co3_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_co3*fact l = lptr_nh4_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_nh4*fact l = lptr_na_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_na*fact l = lptr_ca_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_ca*fact l = lptr_oin_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oin*fact l = lptr_oc_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oc*fact l = lptr_bc_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_bc*fact l = numptr_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_num*fact end do end do end do end do size_loop END SELECT fire_select !-- SAN - end adding fire emissions. ! do mass check final calc if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( & id, config_flags, 2, 'mosaic_ademiss', & dtstep, efact1, dz8w, chem, chem_sum, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & 14, & emis_ant(ims, 1,jms,p_e_pm_10),emis_ant(ims, 1,jms,p_e_pm_25), & emis_ant(ims, 1,jms,p_e_pm25i),emis_ant(ims, 1,jms,p_e_pm25j), & emis_ant(ims, 1,jms,p_e_eci),emis_ant(ims, 1,jms,p_e_ecj), & emis_ant(ims, 1,jms,p_e_orgi),emis_ant(ims, 1,jms,p_e_orgj), & emis_ant(ims, 1,jms,p_e_so4j),emis_ant(ims, 1,jms,p_e_so4c), & emis_ant(ims, 1,jms,p_e_no3j),emis_ant(ims, 1,jms,p_e_no3c), & emis_ant(ims, 1,jms,p_e_orgc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), & emis_ant(ims, 1,jms,p_e_ecc)) ! do seasalt emissions if (seasalt_emiss_active > 0) & call mosaic_seasalt_emiss( & id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, & seas_flux, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, seasalt_emiss_active ) ! if (seasalt_emiss_active == 2) & ! call mosaic_seasalt_emiss( & ! id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, & ! ids,ide, jds,jde, kds,kde, & ! ims,ime, jms,jme, kms,kme, & ! its,ite, jts,jte, kts,kte ) ! ! do dust emissions ! jdf: preliminary version that has not been made generic for situation if (dust_opt == 2) then !czhao+++++++++++++++++++++++++++ call wrf_message("WARNING: You are calling DUSTRAN dust emission scheme with MOSAIC, which is highly experimental and not recommended for use. Please use dust_opt==13") !czhao--------------------------- call mosaic_dust_emiss( slai, ust, smois, ivgtyp, isltyp, & id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) endif !czhao+++++++++++++++++++++++++++ if (dust_opt == 13) then !czhao--------------------------- call mosaic_dust_gocartemis (ktau,dtstep,config_flags%num_soil_layers,alt,u_phy, & v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod, & ivgtyp,isltyp,xland,dx,g, & dust_flux, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) endif return END subroutine mosaic_addemiss !---------------------------------------------------------------------- subroutine mosaic_seasalt_emiss( & id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, & seas_flux, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, seasalt_emiss_active ) ! ! adds seasalt emissions for mosaic aerosol species ! (i.e., seasalt emissions tendencies over time dtstep are applied ! to the aerosol mixing ratios) ! USE module_configure, only: grid_config_rec_type USE module_state_description, only: num_chem, param_first_scalar USE module_data_mosaic_asect IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: id, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL, INTENT(IN ) :: dtstep INTEGER, INTENT(IN) :: seasalt_emiss_active ! 10-m wind speed components (m/s) REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: u10, v10, xland ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg) REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem ! alt = 1.0/(dry air density) in (m3/kg) ! dz8w = layer thickness in (m) REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: alt, dz8w REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: seas_flux ! local variables integer i, j, k, l, l_na, l_cl, n, l_oc integer iphase, itype integer p1st real dum, dumdlo, dumdhi, dumoceanfrac, dumspd10 real factaa, factbb, fracna, fraccl, fracorg real :: ssemfact_numb( maxd_asize, maxd_atype ) real :: ssemfact_mass( maxd_asize, maxd_atype ) ! Factors for Fuentes et al (ACP, 2011, doi:10.5194/acp-11-2585-2011) ! seasalt emission scheme. ! These are for calculating the seawater_OC content dependent factors. real, save :: alpha_f1(4) = & (/ 12.328, 38.077, 102.31, 281.65 /) real, save :: alpha_f2(4) = & (/ 2.2958, 8.0935, 25.251, 46.80 /) real, save :: alpha_f3(4) = & (/ 0.00452, 0.00705, 0.00080, 0.000761 /) real, save :: beta_f1(4) = & (/ 0.0311, -0.031633, 0.013154, -0.0017762 /) real, save :: beta_f2(4) = & (/ -13.916, 35.73, -9.7651, 1.1665 /) real, save :: beta_f3(4) = & (/ 4747.8, 12920.0, 7313.4, 6610.0 /) real :: nti(4), dp0gi(4) ! seawater OC<0.2um concentration (in uM) - first is for low activity, the 2nd for high activity ! High activity conc is from the average for RHaMBLe cruise in high activity regions real, save :: oc02um(2) = (/ 0.0, 280.0 /) ! Organic fraction for seasalt emissions (by size bin). ! These are estimated from Figure 5 of Fuentes et al (ACP, 2011), assuming a value of ! 2 ug/l for Chl-a for the high activity real, save :: org_frac_low_activity(8) = & (/ 0.05, 0.05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /) real, save :: org_frac_high_activity(8) = & (/ 0.2, 0.2, 0.1, 0.01, 0.01, 0.01, 0.01, 0.01 /) p1st = PARAM_FIRST_SCALAR ! for now just do itype=1 itype = 1 iphase = ai_phase ! if using Feuntes et al (2011) then calculate the seawater OC dependent factors if(seasalt_emiss_active .eq. 3 .or. seasalt_emiss_active .eq. 4)then do i=1,4 nti(i) = beta_f1(i) * oc02um(seasalt_emiss_active-2)**2.0 + beta_f2(i) & * oc02um(seasalt_emiss_active-2) + beta_f3(i) dp0gi(i) = alpha_f1(i) + alpha_f2(i) & * exp(-alpha_f3(i)*oc02um(seasalt_emiss_active-2)) end do end if ! compute emissions factors for each size bin do n = 1, nsize_aer(itype) ! changed the lower bound of the emission size limit to match lowest section edge, Qing.Yang@pnl.gov ! DL - 30/3/2012 - select emission scheme (1=original, 3&4=Feuntes et al, 2011) if(seasalt_emiss_active == 1)then ! changed the lower bound of the emission size limit to match lowest section edge, Qing.Yang@pnl.gov ! dumdlo = max( dlo_sect(n,itype), 0.1e-4 ) ! dumdhi = max( dhi_sect(n,itype), 0.1e-4 ) dumdlo = max( dlo_sect(n,itype), 0.02e-4 ) dumdhi = max( dhi_sect(n,itype), 0.02e-4 ) call seasalt_emitfactors_1bin( 1, dumdlo, dumdhi, & ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype) ) elseif(seasalt_emiss_active .eq. 3 .or. seasalt_emiss_active .eq. 4)then call seasalt_emit_feuntes_1bin( dlo_sect(n,itype), dhi_sect(n,itype), & ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype), nti, dp0gi, oc02um(seasalt_emiss_active-2) ) endif ! convert mass emissions factor from (g/m2/s) to (ug/m2/s) ssemfact_mass(n,itype) = ssemfact_mass(n,itype)*1.0e6 end do seas_flux(:,:) = 0.0 ! loop over i,j and apply seasalt emissions k = kts do j = jts, jte do i = its, ite ! Skip this point if over land. xland=1 for land and 2 for water. ! Also, there is no way to differentiate fresh from salt water. ! Currently, this assumes all water is salty. if( xland(i,j) < 1.5 ) cycle !wig: As far as I can tell, only real.exe knows the fractional breakdown ! of land use. So, in wrf.exe, dumoceanfrac will always be 1. dumoceanfrac = 1. !fraction of grid i,j that is salt water dumspd10 = dumoceanfrac* & ( (u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5*3.41) ) ! factaa is (s*m2/kg-air) ! factaa*ssemfact_mass(n) is (s*m2/kg-air)*(ug/m2/s) = ug/kg-air ! factaa*ssemfact_numb(n) is (s*m2/kg-air)*( #/m2/s) = #/kg-air factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j) factbb = factaa * dumspd10 if(seasalt_emiss_active == 1)then seas_flux(i,j) = dumspd10*SUM(ssemfact_mass(1:nsize_aer(itype),itype)) ! apportion seasalt mass emissions assumming that seasalt is pure nacl fracna = mw_na_aer / (mw_na_aer + mw_cl_aer) fraccl = 1.0 - fracna do n = 1, nsize_aer(itype) ! only apply emissions if bin has both na and cl species l_na = lptr_na_aer(n,itype,iphase) l_cl = lptr_cl_aer(n,itype,iphase) if ((l_na >= p1st) .and. (l_cl >= p1st)) then chem(i,k,j,l_na) = chem(i,k,j,l_na) + & factbb * ssemfact_mass(n,itype) * fracna chem(i,k,j,l_cl) = chem(i,k,j,l_cl) + & factbb * ssemfact_mass(n,itype) * fraccl l = numptr_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + & factbb * ssemfact_numb(n,itype) end if end do ! n = 1, nsize_aer(itype) elseif(seasalt_emiss_active.eq.3 .or. seasalt_emiss_active.eq.4)then do n = 1, nsize_aer(itype) ! apportion seasalt mass emissions assumming that seasalt is a ! simple mix of pure nacl and a single primary organic if(seasalt_emiss_active.eq.3)then fracorg = org_frac_low_activity(n) elseif(seasalt_emiss_active.eq.4)then fracorg = org_frac_high_activity(n) endif fracna = mw_na_aer / (mw_na_aer + mw_cl_aer) fraccl = 1.0 - fracna fracna = (1.0-fracorg)*fracna fraccl = (1.0-fracorg)*fraccl ! only apply emissions if bin has both na and cl species l_na = lptr_na_aer(n,itype,iphase) l_cl = lptr_cl_aer(n,itype,iphase) l_oc = lptr_oc_aer(n,itype,iphase) if ((l_na >= p1st) .and. (l_cl >= p1st) .and. (l_oc >= p1st)) then chem(i,k,j,l_na) = chem(i,k,j,l_na) + & factbb * ssemfact_mass(n,itype) * fracna chem(i,k,j,l_cl) = chem(i,k,j,l_cl) + & factbb * ssemfact_mass(n,itype) * fraccl chem(i,k,j,l_oc) = chem(i,k,j,l_oc) + & factbb * ssemfact_mass(n,itype) * fracorg l = numptr_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + & factbb * ssemfact_numb(n,itype) end if end do ! n = 1, nsize_aer(itype) endif end do ! i = its, ite end do ! j = jts, jte return END subroutine mosaic_seasalt_emiss !c---------------------------------------------------------------------- !c following is from gong06b.f in !c /net/cirrus/files1/home/rce/oldfiles1/box/seasaltg !c---------------------------------------------------------------------- subroutine seasalt_emitfactors_1bin( ireduce_smallr_emit, & dpdrylo_cm, dpdryhi_cm, & emitfact_numb, emitfact_surf, emitfact_mass ) !c !c computes seasalt emissions factors for a specifed !c dry particle size range !c dpdrylo_cm = lower dry diameter (cm) !c dpdryhi_cm = upper dry diameter (cm) !c !c number and mass emissions are then computed as !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41) !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41) !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41) !c !c where spd10 = 10 m windspeed in m/s !c !c uses bubble emissions formula (eqn 5a) from !c Gong et al. [JGR, 1997, p 3805-3818] !c !c *** for rdry < rdry_star, this formula overpredicts emissions. !c A strictly ad hoc correction is applied to the formula, !c based on sea-salt size measurements of !c O'Dowd et al. [Atmos Environ, 1997, p 73-80] !c !c *** the correction is only applied when ireduce_smallr_emit > 0 !c implicit none !c subr arguments integer ireduce_smallr_emit real dpdrylo_cm, dpdryhi_cm, & emitfact_numb, emitfact_surf, emitfact_mass !c local variables integer isub_bin, nsub_bin real alnrdrylo real drydens, drydens_43pi_em12, x_4pi_em8 real dum, dumadjust, dumb, dumexpb real dumsum_na, dumsum_ma, dumsum_sa real drwet, dlnrdry real df0drwet, df0dlnrdry, df0dlnrdry_star real relhum real rdry, rdrylo, rdryhi, rdryaa, rdrybb real rdrylowermost, rdryuppermost, rdry_star real rwet, rwetaa, rwetbb real rdry_cm, rwet_cm real sigmag_star real xmdry, xsdry real pi parameter (pi = 3.1415936536) !c c1-c4 are constants for seasalt hygroscopic growth parameterization !c in Eqn 3 and Table 2 of Gong et al. [1997] real c1, c2, c3, c4, onethird parameter (c1 = 0.7674) parameter (c2 = 3.079) parameter (c3 = 2.573e-11) parameter (c4 = -1.424) parameter (onethird = 1.0/3.0) !c dry particle density (g/cm3) drydens = 2.165 !c factor for radius (micrometers) to mass (g) drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12 !c factor for radius (micrometers) to surface (cm2) x_4pi_em8 = 4.0*pi*1.0e-8 !c bubble emissions formula assume 80% RH relhum = 0.80 !c rdry_star = dry radius (micrometers) below which the !c dF0/dr emission formula is adjusted downwards rdry_star = 0.1 if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20 !c sigmag_star = geometric standard deviation used for !c rdry < rdry_star sigmag_star = 1.9 !c initialize sums dumsum_na = 0.0 dumsum_sa = 0.0 dumsum_ma = 0.0 !c rdrylowermost, rdryuppermost = lower and upper !c dry radii (micrometers) for overall integration rdrylowermost = dpdrylo_cm*0.5e4 rdryuppermost = dpdryhi_cm*0.5e4 !c !c "section 1" !c integrate over rdry > rdry_star, where the dF0/dr emissions !c formula is applicable !c (when ireduce_smallr_emit <= 0, rdry_star = -1.0e20, !c and the entire integration is done here) !c if (rdryuppermost .le. rdry_star) goto 2000 !c rdrylo, rdryhi = lower and upper dry radii (micrometers) !c for this part of the integration rdrylo = max( rdrylowermost, rdry_star ) rdryhi = rdryuppermost nsub_bin = 1000 alnrdrylo = log( rdrylo ) dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin !c compute rdry, rwet (micrometers) at lowest size rdrybb = exp( alnrdrylo ) rdry_cm = rdrybb*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetbb = rwet_cm*1.0e4 do 1900 isub_bin = 1, nsub_bin !c rdry, rwet at sub_bin lower boundary are those !c at upper boundary of previous sub_bin rdryaa = rdrybb rwetaa = rwetbb !c compute rdry, rwet (micrometers) at sub_bin upper boundary dum = alnrdrylo + isub_bin*dlnrdry rdrybb = exp( dum ) rdry_cm = rdrybb*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetbb = rwet_cm*1.0e4 !c geometric mean rdry, rwet (micrometers) for sub_bin rdry = sqrt(rdryaa * rdrybb) rwet = sqrt(rwetaa * rwetbb) drwet = rwetbb - rwetaa !c xmdry is dry mass in g xmdry = drydens_43pi_em12 * (rdry**3.0) !c xsdry is dry surface in cm2 xsdry = x_4pi_em8 * (rdry**2.0) !c dumb is "B" in Gong's Eqn 5a !c df0drwet is "dF0/dr" in Gong's Eqn 5a dumb = ( 0.380 - log10(rwet) ) / 0.650 dumexpb = exp( -dumb*dumb) df0drwet = 1.373 * (rwet**(-3.0)) * & (1.0 + 0.057*(rwet**1.05)) * & (10.0**(1.19*dumexpb)) dumsum_na = dumsum_na + drwet*df0drwet dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry 1900 continue !c !c "section 2" !c integrate over rdry < rdry_star, where the dF0/dr emissions !c formula is just an extrapolation and predicts too many emissions !c !c 1. compute dF0/dln(rdry) = (dF0/drwet)*(drwet/dlnrdry) !c at rdry_star !c 2. for rdry < rdry_star, assume dF0/dln(rdry) is lognormal, !c with the same lognormal parameters observed in !c O'Dowd et al. [1997] !c 2000 if (rdrylowermost .ge. rdry_star) goto 3000 !c compute dF0/dln(rdry) at rdry_star rdryaa = 0.99*rdry_star rdry_cm = rdryaa*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetaa = rwet_cm*1.0e4 rdrybb = 1.01*rdry_star rdry_cm = rdrybb*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetbb = rwet_cm*1.0e4 rwet = 0.5*(rwetaa + rwetbb) dumb = ( 0.380 - log10(rwet) ) / 0.650 dumexpb = exp( -dumb*dumb) df0drwet = 1.373 * (rwet**(-3.0)) * & (1.0 + 0.057*(rwet**1.05)) * & (10.0**(1.19*dumexpb)) drwet = rwetbb - rwetaa dlnrdry = log( rdrybb/rdryaa ) df0dlnrdry_star = df0drwet * (drwet/dlnrdry) !c rdrylo, rdryhi = lower and upper dry radii (micrometers) !c for this part of the integration rdrylo = rdrylowermost rdryhi = min( rdryuppermost, rdry_star ) nsub_bin = 1000 alnrdrylo = log( rdrylo ) dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin do 2900 isub_bin = 1, nsub_bin !c geometric mean rdry (micrometers) for sub_bin dum = alnrdrylo + (isub_bin-0.5)*dlnrdry rdry = exp( dum ) !c xmdry is dry mass in g xmdry = drydens_43pi_em12 * (rdry**3.0) !c xsdry is dry surface in cm2 xsdry = x_4pi_em8 * (rdry**2.0) !c dumadjust is adjustment factor to reduce dF0/dr dum = log( rdry/rdry_star ) / log( sigmag_star ) dumadjust = exp( -0.5*dum*dum ) df0dlnrdry = df0dlnrdry_star * dumadjust dumsum_na = dumsum_na + dlnrdry*df0dlnrdry dumsum_ma = dumsum_ma + dlnrdry*df0dlnrdry*xmdry dumsum_sa = dumsum_sa + dlnrdry*df0dlnrdry*xsdry 2900 continue !c !c all done !c 3000 emitfact_numb = dumsum_na emitfact_mass = dumsum_ma emitfact_surf = dumsum_sa return end subroutine seasalt_emitfactors_1bin !c---------------------------------------------------------------------- !c Following is an adaption of the subroutine above. !c It has been modified to include flux from Fuentes et al, ACP, 2010 !c for dealing with smaller particle emissions. !c---------------------------------------------------------------------- subroutine seasalt_emit_feuntes_1bin( & dpdrylo_cm, dpdryhi_cm, & emitfact_numb, emitfact_surf, emitfact_mass, nti, dp0gi, oc02um ) !c !c computes seasalt emissions factors for a specifed !c dry particle size range !c dpdrylo_cm = lower dry diameter (cm) !c dpdryhi_cm = upper dry diameter (cm) !c !c number and mass emissions are then computed as !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41) !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41) !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41) !c !c where spd10 = 10 m windspeed in m/s !c !c Uses bubble emissions formula (eqn 5a) from !c Gong et al. [JGR, 1997, p 3805-3818] for particles !c with dry diameters of 0.45 um or more !c !c For smaller particles we use the parameterisation of !c Fuentes et al, ACP, 2010. !c implicit none ! subr arguments !integer ireduce_smallr_emit real dpdrylo_cm, dpdryhi_cm, & emitfact_numb, emitfact_surf, emitfact_mass real, intent(in) :: nti(4), dp0gi(4), oc02um ! local variables integer isub_bin, nsub_bin, jd, nsub_lower_bin, nsub_upper_bin real alnrdrylo real drydens, drydens_43pi_em12, x_4pi_em8, drydens_16pi_em12, x_pi_em8 real dum, dumadjust, dumb, dumexpb real dumsum_na, dumsum_ma, dumsum_sa real drwet, dlnrdry, ddwet, ddry, dwet real df0drwet, df0dlnrdry, df0dlnrdry_star real relhum real rdry, rdrylo, rdryhi, rdryaa, rdrybb real rdrylowermost, rdryuppermost, rdry_star real rwet, rwetaa, rwetbb real rdry_cm, rwet_cm real sigmag_star real xmdry, xsdry real ddrylo, ddryhi, alogddrylo real ddrybb, ddry_cm, dwet_cm, dwetbb real ddryaa, dwetaa, dlogddry, df0dlogddry real pi parameter (pi = 3.1415936536) ! c1-c4 are constants for seasalt hygroscopic growth parameterization ! in Eqn 3 and Table 2 of Gong et al. [1997] real c1, c2, c3, c4, onethird parameter (c1 = 0.7674) parameter (c2 = 3.079) parameter (c3 = 2.573e-11) parameter (c4 = -1.424) parameter (onethird = 1.0/3.0) ! constants for seasalt/organic aerosol distribution parameterisation ! from Fuentes et al, ACP, 2010 real, save :: width_ssinorg_f(4) = & (/ 1.55, 1.7, 1.5, 1.7 /) real, save :: width_ssorg_f(4) = & (/ 1.55, 1.9, 1.5, 1.7 /) real :: width_ss_f(4), frac ! scale_factor is a combination of: ! 1) (Q) sweep air flow (58.3 cm3/s) ! 2) (Ab) total surface area covered by bubbles (0.0146 m2) ! 4) scaling factor of 3.84e-6 from the whitecap coverage formulation of ! Monahan and O'Muircheartaigh (1980) - as the Gong et al formulation ! includes a whitecap coverage factor too ! this is for converting from dNt to dFp for Fuentes et al (ACP, 2010) real, parameter :: scale_factor = 58.3/(0.0146)*3.84e-6 !0.01533 ! select which distribution width to use for Fuentes et al, if(oc02um.gt.0e0)then width_ss_f = width_ssorg_f else width_ss_f = width_ssinorg_f end if ! dry particle density (g/cm3) drydens = 2.165 ! factor for radius (micrometers) to mass (g) drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12 ! factor for diameter (micrometers) to mass (g) drydens_16pi_em12 = drydens*(1.0/6.0)*pi*1.0e-12 ! factor for radius (micrometers) to surface (cm2) x_4pi_em8 = 4.0*pi*1.0e-8 ! factor for diameter (micrometers) to surface (cm2) x_pi_em8 = pi*1.0e-8 ! bubble emissions formula assume 80% RH relhum = 0.80 ! rdry_star = dry radius (micrometers) below which the ! dF0/dr emission formula from Fuentes et al, ACP, 2010 ! will be used. rdry_star = 0.45 / 2.0 !if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20 ! sigmag_star = geometric standard deviation used for ! rdry < rdry_star sigmag_star = 1.9 ! initialize sums dumsum_na = 0.0 dumsum_sa = 0.0 dumsum_ma = 0.0 ! rdrylowermost, rdryuppermost = lower and upper ! dry radii (micrometers) for overall integration rdrylowermost = dpdrylo_cm*0.5e4 rdryuppermost = dpdryhi_cm*0.5e4 ! "section 1" ! integrate over rdry > rdry_star, where the dF0/dr emissions ! formula from the original subroutine is applicable ! (so using Gong et al for all emission profile) if (rdrylowermost .ge. rdry_star) then ! rdrylo, rdryhi = lower and upper dry radii (micrometers) ! for this part of the integration rdrylo = rdrylowermost rdryhi = rdryuppermost nsub_bin = 1000 alnrdrylo = log( rdrylo ) dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin ! compute rdry, rwet (micrometers) at lowest size rdrybb = exp( alnrdrylo ) rdry_cm = rdrybb*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetbb = rwet_cm*1.0e4 do isub_bin = 1, nsub_bin ! rdry, rwet at sub_bin lower boundary are those ! at upper boundary of previous sub_bin rdryaa = rdrybb rwetaa = rwetbb ! compute rdry, rwet (micrometers) at sub_bin upper boundary dum = alnrdrylo + isub_bin*dlnrdry rdrybb = exp( dum ) rdry_cm = rdrybb*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetbb = rwet_cm*1.0e4 ! geometric mean rdry, rwet (micrometers) for sub_bin rdry = sqrt(rdryaa * rdrybb) rwet = sqrt(rwetaa * rwetbb) drwet = rwetbb - rwetaa ! xmdry is dry mass in g xmdry = drydens_43pi_em12 * (rdry**3.0) ! xsdry is dry surface in cm2 xsdry = x_4pi_em8 * (rdry**2.0) ! dumb is "B" in Gong's Eqn 5a ! df0drwet is "dF0/dr" in Gong's Eqn 5a dumb = ( 0.380 - log10(rwet) ) / 0.650 dumexpb = exp( -dumb*dumb) df0drwet = 1.373 * (rwet**(-3.0)) * & (1.0 + 0.057*(rwet**1.05)) * & (10.0**(1.19*dumexpb)) dumsum_na = dumsum_na + drwet*df0drwet dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry end do ! "section 2" ! integrate over rdry < rdry_star for old parameterisation, and use ! Fuentes et al below that ! ! 1. for rdry < rdry_star, use Fuentes et al parameterisation ! 2. for rdry > rdry_star, use Gong et al parameterisation, as above elseif (rdryuppermost .gt. rdry_star) then ! determine the fraction of size bin below rdry_star frac = (log(rdry_star)-log(rdrylowermost)) / (log(rdryuppermost)-log(rdrylowermost)) nsub_lower_bin = floor(frac*1000.0) ! calc number of size bins to use for section below rdry_star nsub_upper_bin = 1000-nsub_lower_bin ! use remaining size bins for section above rdry_star ! 1................. ! ddrylo, ddryhi = lower and upper dry diameter (micrometers) ! for this part of the integration ddrylo = rdrylowermost*2.0 ddryhi = rdry_star*2.0 nsub_bin = nsub_lower_bin alogddrylo = log10( ddrylo ) dlogddry = (log10( ddryhi ) - alogddrylo)/nsub_bin ! compute ddry, dwet (micrometers) at lowest size ddrybb = 10.0**( alogddrylo ) ddry_cm = ddrybb*1.0e-4 dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ & ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird dwetbb = dwet_cm*1.0e4 do isub_bin = 1, nsub_bin ! ddry, dwet at sub_bin lower boundary are those ! at upper boundary of previous sub_bin ddryaa = ddrybb !dwetaa = dwetbb ! compute ddry, dwet (micrometers) at sub_bin upper boundary dum = alogddrylo + isub_bin*dlogddry ddrybb = 10.0**( dum ) ddry_cm = ddrybb*1.0e-4 !dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ & ! ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird !dwetbb = dwet_cm*1.0e4 ! geometric mean rdry, rwet (micrometers) for sub_bin ddry = sqrt(ddryaa * ddrybb) !dwet = sqrt(dwetaa * dwetbb) ! xmdry is dry mass in g xmdry = drydens_16pi_em12 * (ddry**3.0) ! xsdry is dry surface in cm2 xsdry = x_pi_em8 * (ddry**2.0) ! Equation 2 from Fuentes et al (ACP, 2011). Wet diameter is scaled by a factor ! of 1e3 to convert from um to nm for this calculation df0dlogddry = 0.0 do jd = 1, 4 df0dlogddry = df0dlogddry + nti(jd)/( (2.0*pi)**0.5 * log10(width_ss_f(jd)) ) & * exp(-1.0/2.0 * (log10(ddry*1e3/dp0gi(jd))/log10(width_ss_f(jd)))**2.0 ) end do dumsum_na = dumsum_na + dlogddry*df0dlogddry*scale_factor dumsum_ma = dumsum_ma + dlogddry*df0dlogddry*scale_factor*xmdry dumsum_sa = dumsum_sa + dlogddry*df0dlogddry*scale_factor*xsdry end do ! 2................ ! rdrylo, rdryhi = lower and upper dry radii (micrometers) ! for this part of the integration rdrylo = rdry_star rdryhi = rdryuppermost nsub_bin = nsub_upper_bin alnrdrylo = log( rdrylo ) dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin ! compute rdry, rwet (micrometers) at lowest size rdrybb = exp( alnrdrylo ) rdry_cm = rdrybb*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetbb = rwet_cm*1.0e4 do isub_bin = 1, nsub_bin ! rdry, rwet at sub_bin lower boundary are those ! at upper boundary of previous sub_bin rdryaa = rdrybb rwetaa = rwetbb ! compute rdry, rwet (micrometers) at sub_bin upper boundary dum = alnrdrylo + isub_bin*dlnrdry rdrybb = exp( dum ) rdry_cm = rdrybb*1.0e-4 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ & ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird rwetbb = rwet_cm*1.0e4 ! geometric mean rdry, rwet (micrometers) for sub_bin rdry = sqrt(rdryaa * rdrybb) rwet = sqrt(rwetaa * rwetbb) drwet = rwetbb - rwetaa ! xmdry is dry mass in g xmdry = drydens_43pi_em12 * (rdry**3.0) ! xsdry is dry surface in cm2 xsdry = x_4pi_em8 * (rdry**2.0) ! dumb is "B" in Gong's Eqn 5a ! df0drwet is "dF0/dr" in Gong's Eqn 5a dumb = ( 0.380 - log10(rwet) ) / 0.650 dumexpb = exp( -dumb*dumb) df0drwet = 1.373 * (rwet**(-3.0)) * & (1.0 + 0.057*(rwet**1.05)) * & (10.0**(1.19*dumexpb)) dumsum_na = dumsum_na + drwet*df0drwet dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry end do ! "section 3" ! where rdry < rdry_star for the whole bin we use Fuentes et al ! for all emissions ! else ! ddrylo, ddryhi = lower and upper dry diameter (micrometers) ! for this part of the integration ddrylo = rdrylowermost*2.0 ddryhi = rdryuppermost*2.0 nsub_bin = 1000 alogddrylo = log10( ddrylo ) dlogddry = (log10( ddryhi ) - alogddrylo)/nsub_bin ! compute ddry, dwet (micrometers) at lowest size ddrybb = 10.0**( alogddrylo ) ddry_cm = ddrybb*1.0e-4 dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ & ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird dwetbb = dwet_cm*1.0e4 do isub_bin = 1, nsub_bin ! ddry, dwet at sub_bin lower boundary are those ! at upper boundary of previous sub_bin ddryaa = ddrybb !dwetaa = dwetbb ! compute ddry, dwet (micrometers) at sub_bin upper boundary dum = alogddrylo + isub_bin*dlogddry ddrybb = 10.0**( dum ) ddry_cm = ddrybb*1.0e-4 !dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ & ! ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird !dwetbb = dwet_cm*1.0e4 ! geometric mean rdry, rwet (micrometers) for sub_bin ddry = sqrt(ddryaa * ddrybb) !dwet = sqrt(dwetaa * dwetbb) ! xmdry is dry mass in g xmdry = drydens_16pi_em12 * (ddry**3.0) ! xsdry is dry surface in cm2 xsdry = x_pi_em8 * (ddry**2.0) ! Equation 2 from Fuentes et al (ACP, 2011). Wet diameter is scaled by a factor ! of 1e3 to convert from um to nm for this calculation df0dlogddry = 0.0 do jd = 1, 4 df0dlogddry = df0dlogddry + nti(jd)/( (2.0*pi)**0.5 * log10(width_ss_f(jd)) ) & * exp(-1.0/2.0 * (log10(ddry*1e3/dp0gi(jd))/log10(width_ss_f(jd)))**2.0 ) end do dumsum_na = dumsum_na + dlogddry*df0dlogddry*scale_factor dumsum_ma = dumsum_ma + dlogddry*df0dlogddry*scale_factor*xmdry dumsum_sa = dumsum_sa + dlogddry*df0dlogddry*scale_factor*xsdry end do end if !c !c all done !c emitfact_numb = dumsum_na emitfact_mass = dumsum_ma emitfact_surf = dumsum_sa return end subroutine seasalt_emit_feuntes_1bin END MODULE module_mosaic_addemiss !---------------------------------------------------------------------- subroutine mosaic_dust_emiss( slai,ust, smois, ivgtyp, isltyp, & id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ! ! adds dust emissions for mosaic aerosol species (i.e. emission tendencies ! over time dtstep are applied to the aerosol mixing ratios) ! ! This is a simple dust scheme based on Shaw et al. (2008) to appear in ! Atmospheric Environment, recoded by Jerome Fast ! ! NOTE: ! 1) This version only works with the 8-bin version of MOSAIC. ! 2) Dust added to MOSAIC's other inorganic specie, OIN. If Ca and CO3 are ! activated in the Registry, a small fraction also added to Ca and CO3. ! 3) The main departure from Shaw et al., is now alphamask is computed since ! the land-use categories in that paper and in WRF differ. WRF currently ! does not have that many land-use categories and adhoc assumptions had to ! be made. This version was tested for Mexico in the dry season. The main ! land-use categories in WRF that are likely dust sources are grass, shrub, ! and savannna (that WRF has in the desert regions of NW Mexico). Having ! dust emitted from these types for other locations and other times of the ! year is not likely to be valid. ! 4) An upper bound on ustar was placed because the surface parameterizations ! in WRF can produce unrealistically high values that lead to very high ! dust emission rates. ! 5) Other departures' from Shaw et al. noted below, but are probably not as ! important as 2) and 3). ! 6) Now the dust added into dust species USE module_configure, only: grid_config_rec_type USE module_state_description, only: num_chem, param_first_scalar USE module_data_mosaic_asect IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: id, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL, INTENT(IN ) :: dtstep ! 10-m wind speed components (m/s) REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: u10, v10, xland, slai, ust INTEGER, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: ivgtyp, isltyp ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg) REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem ! alt = 1.0/(dry air density) in (m3/kg) ! dz8w = layer thickness in (m) REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: alt, dz8w REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , & INTENT(INOUT) :: smois ! local variables integer i, j, k, l, l_oin, l_ca, l_co3, n, ii integer iphase, itype, izob integer p1st real dum, dumdlo, dumdhi, dumlandfrac, dumspd10 real factaa, factbb, fracoin, fracca, fracco3, fractot real ustart, ustar1, ustart0 real alphamask, f8, f50, f51, f52, wetfactor, sumdelta, ftot real smois_grav, wp, pclay real :: beta(4,7) real :: gamma(4), delta(4) real :: sz(8) real :: dustflux, densdust, mass1part real :: c_const p1st = param_first_scalar ! ! from Nickovic et al., JGR, 2001 and Shaw et al. 2007 ! beta: fraction of clay, small silt, large silt, and sand correcsponding to Zobler class (7) ! beta (1,*) for 0.5-1 um ! beta (2,*) for 1-10 um ! beta (3,*) for 10-25 um ! beta (4,*) for 25-50 um ! beta(1,1)=0.12 beta(2,1)=0.04 beta(3,1)=0.04 beta(4,1)=0.80 beta(1,2)=0.34 beta(2,2)=0.28 beta(3,2)=0.28 beta(4,2)=0.10 beta(1,3)=0.45 beta(2,3)=0.15 beta(3,3)=0.15 beta(4,3)=0.25 beta(1,4)=0.12 beta(2,4)=0.09 beta(3,4)=0.09 beta(4,4)=0.70 beta(1,5)=0.40 beta(2,5)=0.05 beta(3,5)=0.05 beta(4,5)=0.50 beta(1,6)=0.34 beta(2,6)=0.18 beta(3,6)=0.18 beta(4,6)=0.30 beta(1,7)=0.22 beta(2,7)=0.09 beta(3,7)=0.09 beta(4,7)=0.60 gamma(1)=0.08 gamma(2)=1.00 gamma(3)=1.00 gamma(4)=0.12 ! ! * Mass fractions for each size bin. These values were recommended by ! Natalie Mahowold, with bins 7 and 8 the same as bins 3 and 4 from CAM. ! * Changed slightly since Natelie's estimates do not add up to 1.0 ! * This would need to be made more generic for other bin sizes. ! sz(1)=0 ! sz(2)=1.78751e-06 ! sz(3)=0.000273786 ! sz(4)=0.00847978 ! sz(5)=0.056055 ! sz(6)=0.0951896 ! sz(7)=0.17 ! sz(8)=0.67 ! sz(1)=1.78751e-06 ! sz(2)=0.00875357 ! sz(3)=0.1512446 ! sz(4)=0.84 itype=1 if (nsize_aer(itype) .eq. 8) then !BSINGH(12/11/2013): Based on the suggestions by Manish Shrivastva, sz variable is modified below. !Original values are commented out !sz(1)=0.0 !sz(2)=0.0 !sz(3)=0.0005 !sz(4)=0.0095 !sz(5)=0.01 !sz(6)=0.06 !sz(7)=0.20 !sz(8)=0.72 sz(1)=0 sz(2)=1.78751e-06 sz(3)=0.000273786 sz(4)=0.00847978 sz(5)=0.056055 sz(6)=0.0951896 sz(7)=0.17 sz(8)=0.67 else if (nsize_aer(itype) .eq. 4) then !BSINGH(12/11/2013): Based on the suggestions by Manish Shrivastva, sz variable is modified below. !Original values are commented out !sz(1)=0.0 !sz(2)=0.01 !sz(3)=0.07 !sz(4)=0.92 sz(1)=1.78751e-06 sz(2)=0.00875357 sz(3)=0.1512446 sz(4)=0.84 sz(5)=0.0 sz(6)=0.0 sz(7)=0.0 sz(8)=0.0 endif ! for now just do itype=1 itype = 1 iphase = ai_phase ! loop over i,j and apply dust emissions k = kts do 1830 j = jts, jte do 1820 i = its, ite if( xland(i,j) > 1.5 ) cycle ! compute wind speed anyway, even though ustar is used below dumlandfrac = 1. dumspd10=(u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5) if(dumspd10 >= 5.0) then dumspd10 = dumlandfrac* & ( dumspd10*dumspd10*(dumspd10-5.0)) else dumspd10=0. endif ! part1 - compute vegetation mask ! ! * f8, f50, f51, f52 refer to vegetation classes from the Olsen categories ! for desert, sand desert, grass semi-desert, and shrub semi-desert ! * in WRF, vegetation types 7, 8 and 10 are grassland, shrubland, and savanna ! that are dominate types in Mexico and probably have some erodable surface ! during the dry season ! * currently modified these values so that only a small fraction of cell ! area is erodable ! * these values are highly tuneable! alphamask=0.001 if (ivgtyp(i,j) .eq. 7) then f8=0.005 f50=0.00 f51=0.10 f51=0.066 f52=0.00 alphamask=(f8+f50)*1.0+(f51+f52)*0.5 endif if (ivgtyp(i,j) .eq. 8) then f8=0.010 f50=0.00 f51=0.00 f52=0.15 f52=0.10 alphamask=(f8+f50)*1.0+(f51+f52)*0.5 endif if (ivgtyp(i,j) .eq. 10) then f8=0.00 f50=0.00 f51=0.01 f52=0.00 alphamask=(f8+f50)*1.0+(f51+f52)*0.5 endif ! part2 - zobler ! ! * in Shaw's paper, dust is computed for 4 size ranges: ! 0.5-1 um ! 1-10 um ! 10-25 um ! 25-50 um ! * Shaw's paper also accounts for sub-grid variability in soil ! texture, but here we just assume the same soil texture for each ! grid cell ! * since MOSAIC is currently has a maximum size range up to 10 um, ! neglect upper 2 size ranges and lowest size range (assume small) ! * map WRF soil classes arbitrarily to Zolber soil textural classes ! * skip dust computations for WRF soil classes greater than 13, i.e. ! do not compute dust over water, bedrock, and other surfaces ! * should be skipping for water surface at this point anyway ! izob=0 if(isltyp(i,j).eq.1) izob=1 if(isltyp(i,j).eq.2) izob=1 if(isltyp(i,j).eq.3) izob=4 if(isltyp(i,j).eq.4) izob=2 if(isltyp(i,j).eq.5) izob=2 if(isltyp(i,j).eq.6) izob=2 if(isltyp(i,j).eq.7) izob=7 if(isltyp(i,j).eq.8) izob=2 if(isltyp(i,j).eq.9) izob=6 if(isltyp(i,j).eq.10) izob=5 if(isltyp(i,j).eq.11) izob=2 if(isltyp(i,j).eq.12) izob=3 if(isltyp(i,j).ge.13) izob=0 if(izob.eq.0) goto 1840 ! ! part3 - dustprod ! do ii=1,4 delta(ii)=0.0 enddo sumdelta=0.0 do ii=1,4 delta(ii)=beta(ii,izob)*gamma(ii) if(ii.lt.4) then sumdelta=sumdelta+delta(ii) endif enddo do ii=1,4 delta(ii)=delta(ii)/sumdelta enddo ! part4 - wetness ! ! * assume dry for now, have passed in soil moisture to this routine ! but needs to be included here ! * wetfactor less than 1 would reduce dustflux ! * convert model soil moisture (m3/m3) to gravimetric soil moisture ! (mass of water / mass of soil in %) assuming a constant density ! for soil pclay=beta(1,izob)*100. wp=0.0014*pclay*pclay+0.17*pclay smois_grav=(smois(i,1,j)/2.6)*100. if(smois_grav.gt.wp) then wetfactor=sqrt(1.0+1.21*(smois_grav-wp)**0.68) else wetfactor=1.0 endif ! wetfactor=1.0 ! part5 - dustflux ! lower bound on ustar = 20 cm/s as in Shaw et al, but set upper ! bound to 100 cm/s c_const=1.e-14 ! default ustar1=ust(i,j)*100.0 if(ustar1.gt.100.0) ustar1=100.0 ustart0=20.0 ustart=ustart0*wetfactor if(ustar1.le.ustart) then dustflux=0.0 else dustflux=c_const*(ustar1**4)*(1.0-(ustart/ustar1)) endif dustflux=dustflux*10.0 ! units kg m-2 s-1 ftot=0.0 do ii=1,2 ftot=ftot+dustflux*alphamask*delta(ii) enddo ! convert to ug m-2 s-1 ftot=ftot*1.0e+09 ! apportion other inorganics only factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j) factbb = factaa * ftot fracoin = 0.97 fracca = 0.03*0.4 fracco3 = 0.03*0.6 fractot = fracoin + fracca + fracco3 ! if (ivgtyp(i,j) .eq. 8) print*,'jdf',i,j,ustar1,ustart0,factaa,ftot do 1810 n = 1, nsize_aer(itype) l_oin = lptr_oin_aer(n,itype,iphase) l_ca = lptr_ca_aer(n,itype,iphase) l_co3 = lptr_co3_aer(n,itype,iphase) if (l_oin >= p1st) then chem(i,k,j,l_oin) = chem(i,k,j,l_oin) + & factbb * sz(n) * fracoin if (l_ca >= p1st) & chem(i,k,j,l_ca) = chem(i,k,j,l_ca) + & factbb * sz(n) * fracca if (l_co3 >= p1st) & chem(i,k,j,l_co3) = chem(i,k,j,l_co3) + & factbb * sz(n) * fracco3 ! mass1part is mass of a single particle in ug, density of dust ~2.5 g cm-3 densdust=2.5 mass1part=0.523598*(dcen_sect(n,itype)**3)*densdust*1.0e06 l = numptr_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + & factbb * sz(n) * fractot / mass1part end if 1810 continue 1840 continue 1820 continue 1830 continue return END subroutine mosaic_dust_emiss !==================================================================================== !add another dust emission scheme following GOCART mechanism --czhao 09/17/2009 !==================================================================================== subroutine mosaic_dust_gocartemis (ktau,dt,num_soil_layers,alt,u_phy, & v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod, & ivgtyp,isltyp,xland,dx,g, & dust_flux, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) USE module_data_gocart_dust USE module_configure USE module_state_description USE module_model_constants, ONLY: mwdry USE module_data_mosaic_asect IMPLICIT NONE INTEGER, INTENT(IN ) :: ktau, num_soil_layers, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER,DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) :: & ivgtyp, & isltyp REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , & INTENT(INOUT) :: smois REAL, DIMENSION( ims:ime , jms:jme, 3 ) , & INTENT(IN ) :: erod REAL, DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) :: & u10, & v10, & xland REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: & alt, & dz8w,p8w, & u_phy,v_phy,rho_phy REAL, INTENT(IN ) :: dt,dx,g REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: dust_flux ! ! local variables ! integer :: nmx,i,j,k,ndt,imx,jmx,lmx integer ilwi, start_month real*8, DIMENSION (3) :: erodin real*8, DIMENSION (5) :: bems real*8 w10m,gwet,airden,airmas real*8 cdustemis,jdustemis,cdustcon,jdustcon real*8 cdustdens,jdustdens,mass1part,jdustdiam,cdustdiam,dp_meanvol_tmp real factaa, factbb, fracoin, fracca, fracco3, fractot real*8 dxy real*8 conver,converi real dttt real soilfacj,rhosoilj,rhosoilc integer p1st,l_oin,l,n real :: densdust real sz(8),ftot,frac_10um real totalemis,accfrac,corfrac,rscale integer iphase, itype rscale=1.01 ! to account for the emitted dust with radius larger than 10 um accfrac=0.07 ! assign 7% to accumulation mode corfrac=0.93 ! assign 93% to coarse mode accfrac=accfrac*rscale corfrac=corfrac*rscale p1st = param_first_scalar !czhao give an example of MOSAIC 8 size bins here (unit: um) from Jerome 2005 paper ! 0.039-0.078; 0.078-0.156; 0.156-0.312; 0.312-0.625; 0.625-1.25; 1.25-2.5; 2.5-5.0; 5.0-10.0 !in the model the size bound for 8 bins is calculated in the code using dlo_sec and dhi_sec !size3 ! sz(1)=0.0 ! sz(2)=0.0 ! sz(3)=0.0003 ! sz(4)=0.0048 ! sz(5)=0.0130 ! sz(6)=0.0250 ! sz(7)=0.1400 ! sz(8)=0.8169 !the size is following the best match2 of modal distribution !we have ~3.6% within 1 um radius and ~7.4% within 1.25 um radius !versus ~6.0% within 1 um radius following GOCART original and ~9.3% within 1.25 um radius ! totally ~82% are within the MOSAIC size range 0.03-10 um diameter ! now the dust distribution follows the parameter defined in module_data_sorgam.F sz(1)=1.0e-8 sz(2)=1.0e-6 sz(3)=3.0e-4 sz(4)=3.5e-3 sz(5)=0.018 sz(6)=0.070 sz(7)=0.2595 sz(8)=0.42 ! frac_10um=sum(sz(1:8)) ! sz(:)=sz(:)/sum(sz(1:8)) ! relative ratio in diameter range of 0-10 um ! for now just do itype=1 itype = 1 iphase = ai_phase ! added option for 4bin WRF IF (nsize_aer(itype) == 4) THEN sz(1) = sz(1) + sz(2) sz(2) = sz(3) + sz(4) sz(3) = sz(5) + sz(6) sz(4) = sz(7) + sz(8) ENDIF conver=1.e-9 converi=1.e9 dust_flux(:,:) = 0.0 ! ! number of dust bins nmx=5 k=kts do j=jts,jte do i=its,ite ! ! don't do dust over water!!! if(xland(i,j).lt.1.5)then ilwi=1 start_month = 3 ! it doesn't matter, ch_dust is not a month dependent now, a constant w10m=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j)) airmas=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g ! kg ! we don't trust the u10,v10 values, if model layers are very thin near surface if(dz8w(i,kts,j).lt.12.)w10m=sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j)) !erodin(1)=erod(i,j,1)/dx/dx ! czhao erod shouldn't be scaled to the area, because it's a fraction !erodin(2)=erod(i,j,2)/dx/dx !erodin(3)=erod(i,j,3)/dx/dx erodin(1)=erod(i,j,1) erodin(2)=erod(i,j,2) erodin(3)=erod(i,j,3) ! ! volumetric soil moisture over porosity gwet=smois(i,1,j)/porosity(isltyp(i,j)) ndt=ifix(dt) airden=rho_phy(i,kts,j) dxy=dx*dx call mosaic_source_du( nmx, dt, & erodin, ilwi, dxy, w10m, gwet, airden, airmas, & bems,start_month,g) !bems: kg/timestep/cell !sum up the dust emission from 0.1-10 um in radius ! unit change from kg/timestep/cell to ug/m2/s totalemis=(sum(bems(1:5))/dt)*converi/dxy dust_flux(i,j) = totalemis !totalemis=totalemis*rscale !to account for the particles larger than 10 um ! based on assumed size distribution jdustemis = totalemis*accfrac ! accumulation mode cdustemis = totalemis*corfrac ! coarse mode ftot=jdustemis+cdustemis ! total dust emission in ug/m2/s ! ftot=ftot*frac_10um ! only in diameter range of 0-10um ! apportion other inorganics only factaa = (dt/dz8w(i,k,j))*alt(i,k,j) factbb = factaa * ftot fracoin = 0.97 fracca = 0.03*0.4 fracco3 = 0.03*0.6 fractot = fracoin do 1810 n = 1, nsize_aer(itype) l_oin = lptr_oin_aer(n,itype,iphase) if (l_oin >= p1st) then chem(i,k,j,l_oin) = chem(i,k,j,l_oin) + & factbb * sz(n) * fracoin ! mass1part is mass of a single particle in ug, density of dust ~2.5 g cm-3 if (n <= 5) densdust=2.5 if (n > 5 ) densdust=2.65 ! added option for 4bin WRF if (nsize_aer(itype) == 4) then if (n <= 2) densdust=2.5 if (n > 2 ) densdust=2.65 endif mass1part=0.523598*(dcen_sect(n,itype)**3)*densdust*1.0e06 l = numptr_aer(n,itype,iphase) if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + & factbb * sz(n) * fractot / mass1part end if 1810 continue endif ! xland enddo ! i enddo ! j end subroutine mosaic_dust_gocartemis SUBROUTINE mosaic_source_du( nmx, dt1, & erod, ilwi, dxy, w10m, gwet, airden, airmas, & bems,month,g0) ! **************************************************************************** ! * Evaluate the source of each dust particles size classes (kg/m3) ! * by soil emission. ! * Input: ! * EROD Fraction of erodible grid cell (-) ! * for 1: Sand, 2: Silt, 3: Clay ! * DUSTDEN Dust density (kg/m3) ! * DXY Surface of each grid cell (m2) ! * AIRVOL Volume occupy by each grid boxes (m3) ! * NDT1 Time step (s) ! * W10m Velocity at the anemometer level (10meters) (m/s) ! * u_tresh Threshold velocity for particule uplifting (m/s) ! * CH_dust Constant to fudge the total emission of dust (s2/m2) ! * ! * Output: ! * DSRC Source of each dust type (kg/timestep/cell) ! * ! * Working: ! * SRC Potential source (kg/m/timestep/cell) ! * ! **************************************************************************** USE module_data_gocart_dust INTEGER, INTENT(IN) :: nmx REAL*8, INTENT(IN) :: erod(ndcls) INTEGER, INTENT(IN) :: ilwi,month REAL*8, INTENT(IN) :: w10m, gwet REAL*8, INTENT(IN) :: dxy REAL*8, INTENT(IN) :: airden, airmas REAL*8, INTENT(OUT) :: bems(nmx) REAL*8 :: den(nmx), diam(nmx) REAL*8 :: tsrc, u_ts0, cw, u_ts, dsrc, srce REAL, intent(in) :: g0 REAL :: rhoa, g,dt1 INTEGER :: i, j, n, m, k ! default is 1 ug s2 m-5 == 1.e-9 kg s2 m-5 !ch_dust(:,:)=0.8D-9 ! ch_dust is defined here instead of in the chemics_ini.F if with SORGAM -czhao ch_dust(:,:)=1.0D-9 ! default !ch_dust(:,:)=0.65D-9 ! ch_dust is scaled for Sahara ! executable statemenst DO n = 1, nmx ! Threshold velocity as a function of the dust density and the diameter from Bagnold (1941) den(n) = den_dust(n)*1.0D-3 diam(n) = 2.0*reff_dust(n)*1.0D2 g = g0*1.0E2 ! Pointer to the 3 classes considered in the source data files m = ipoint(n) tsrc = 0.0 rhoa = airden*1.0D-3 u_ts0 = 0.13*1.0D-2*SQRT(den(n)*g*diam(n)/rhoa)* & SQRT(1.0+0.006/den(n)/g/(diam(n))**2.5)/ & SQRT(1.928*(1331.0*(diam(n))**1.56+0.38)**0.092-1.0) ! Case of surface dry enough to erode IF (gwet < 0.5) THEN ! Pete's modified value ! IF (gwet < 0.2) THEN u_ts = MAX(0.0D+0,u_ts0*(1.2D+0+2.0D-1*LOG10(MAX(1.0D-3, gwet)))) ELSE ! Case of wet surface, no erosion u_ts = 100.0 END IF srce = frac_s(n)*erod(m)*dxy ! (m2) IF (ilwi == 1 ) THEN dsrc = ch_dust(n,month)*srce*w10m**2 & * (w10m - u_ts)*dt1 ! (kg) ELSE dsrc = 0.0 END IF IF (dsrc < 0.0) dsrc = 0.0 ! Update dust mixing ratio at first model level. !tc(n) = tc(n) + dsrc / airmas !kg/kg-dryair -czhao bems(n) = dsrc ! kg/timestep/cell ENDDO END SUBROUTINE mosaic_source_du !=========================================================================== !---------------------------------------------------------------------- subroutine addemiss_masscheck( id, config_flags, iflagaa, fromwhere, & dtstep, efact1, dz8w, chem, chem_sum, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & nemit, & e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, & e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21 ) ! ! produces test diagnostics for "addemiss" routines ! ! 1. computes {sum over i,j,k ( chem * dz8w )} before and after ! emissions tendencies are added to chem, ! then prints (sum_after - sum_before)/(dtstep*efact1) ! 2. computes {sum over i,j,k ( e_xxx )}, then prints them ! the two should be equal ! USE module_configure, only: grid_config_rec_type USE module_state_description, only: num_chem IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: id, iflagaa, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & nemit REAL, INTENT(IN ) :: dtstep, efact1 ! trace species mixing ratios REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(IN ) :: chem ! trace species integrals DOUBLE PRECISION, DIMENSION( num_chem ), & INTENT(INOUT ) :: chem_sum ! layer thickness (m) REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: dz8w ! emissions ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & REAL, DIMENSION( ims:ime , kms:config_flags%kemit , jms:jme ), & INTENT(IN ) :: & e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, & e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21 character(len=*), intent(in) :: fromwhere ! local variables integer, parameter :: nemit_maxd = 21 integer :: i, j, k, l double precision :: chem_sum_prev real :: fact real :: emit_sum(nemit_maxd) ! compute column integral, summed over i-j grids ! compute {sum over i,j,k ( chem * dz8w ) } ! on second pass (iflagaa==2), subtract the pass-one sum do 1900 l = 1, num_chem chem_sum_prev = chem_sum(l) chem_sum(l) = 0.0 do j = jts, jte do k = kts, kte do i = its, ite chem_sum(l) = chem_sum(l) + dble( chem(i,k,j,l)*dz8w(i,k,j) ) end do end do end do if (iflagaa == 2) chem_sum(l) = (chem_sum(l) - chem_sum_prev) 1900 continue if (iflagaa /= 2) return ! compute {sum over i,j,k ( e_xxx ) } emit_sum(:) = 0.0 do 2900 l = 1, min(nemit,nemit_maxd) do j = jts, jte do k = kts, min(config_flags%kemit,kte) do i = its, ite if (l== 1) emit_sum(l) = emit_sum(l) + e01(i,k,j) if (l== 2) emit_sum(l) = emit_sum(l) + e02(i,k,j) if (l== 3) emit_sum(l) = emit_sum(l) + e03(i,k,j) if (l== 4) emit_sum(l) = emit_sum(l) + e04(i,k,j) if (l== 5) emit_sum(l) = emit_sum(l) + e05(i,k,j) if (l== 6) emit_sum(l) = emit_sum(l) + e06(i,k,j) if (l== 7) emit_sum(l) = emit_sum(l) + e07(i,k,j) if (l== 8) emit_sum(l) = emit_sum(l) + e08(i,k,j) if (l== 9) emit_sum(l) = emit_sum(l) + e09(i,k,j) if (l==10) emit_sum(l) = emit_sum(l) + e10(i,k,j) if (l==11) emit_sum(l) = emit_sum(l) + e11(i,k,j) if (l==12) emit_sum(l) = emit_sum(l) + e12(i,k,j) if (l==13) emit_sum(l) = emit_sum(l) + e13(i,k,j) if (l==14) emit_sum(l) = emit_sum(l) + e14(i,k,j) if (l==15) emit_sum(l) = emit_sum(l) + e15(i,k,j) if (l==16) emit_sum(l) = emit_sum(l) + e16(i,k,j) if (l==17) emit_sum(l) = emit_sum(l) + e17(i,k,j) if (l==18) emit_sum(l) = emit_sum(l) + e18(i,k,j) if (l==19) emit_sum(l) = emit_sum(l) + e19(i,k,j) if (l==20) emit_sum(l) = emit_sum(l) + e20(i,k,j) if (l==21) emit_sum(l) = emit_sum(l) + e21(i,k,j) end do end do end do 2900 continue ! output the chem_sum and emit_sum print 9110, fromwhere, its, ite, jts, jte print 9100, 'chem_sum' fact = 1.0/(dtstep*efact1) print 9120, (l, fact*chem_sum(l), l=1,num_chem) print 9100, 'emit_sum' print 9120, (l, emit_sum(l), l=1,min(nemit,nemit_maxd)) 9100 format( a ) 9110 format( / 'addemiss_masscheck output, fromwhere = ', a / & 'its, ite, jts, jte =', 4i5 ) 9120 format( 5( i5, 1pe11.3 ) ) return END subroutine addemiss_masscheck