MODULE module_cam_mam_addemiss
!WRF:MODEL_LAYER:CHEMICS
!----------------------------------------------------------------------
! module_cam_mam_addemiss.F
!
! adapted from module_mosaic_addemiss.F in nov, 2010 by r.c.easter
!
!----------------------------------------------------------------------
#include "MODAL_AERO_CPP_DEFINES.h"
! 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
private
public :: cam_mam_addemiss
CONTAINS
!----------------------------------------------------------------------
subroutine cam_mam_addemiss( id, dtstep, u10, v10, alt, dz8w, xland, &
config_flags, chem, slai, ust, smois, ivgtyp, isltyp, &
emis_ant,ebio_iso,ebio_olt,ebio_oli,rho_phy, &
dust_emiss_active, seasalt_emiss_active, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!
! adds emissions for cam-mam aerosol species
! (i.e., emissions tendencies over time dtstep are applied
! to the aerosol mass and number mixing ratios)
!
! currently this routine
! - applies emissions held in the emis_ant array
! - does online sea-salt emissions using a modified version of the
! mosaic sea-salt emissions routine
! - does online dust emissions using a modified version of the
! mosaic dust emissions routine
!
! still to do
! - for emis_ant emissions, implement emitted sizes similar to those used in cam
! - implement the cam sea-salt and dust emissions routines
!
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_cam_mam_asect, only: ai_phase, &
dens_so4_aer, dens_nh4_aer, dens_pom_aer, &
dens_bc_aer, dens_dust_aer, dens_seas_aer, &
lptr_so4_aer, lptr_nh4_aer, lptr_pom_aer, &
lptr_bc_aer, lptr_dust_aer, lptr_seas_aer, &
maxd_asize, maxd_atype, nsize_aer, ntype_aer, numptr_aer
USE modal_aero_data, only: &
modeptr_accum, modeptr_aitken, modeptr_coarse, &
modeptr_coardust, modeptr_finedust, modeptr_pcarbon, &
modeptr_fineseas, modeptr_coarseas
USE module_cam_mam_init, only: &
pom_emit_1p4_factor, so4_emit_1p2_factor
USE module_data_sorgam, only: &
dgvem_i, dgvem_j, dgvem_c, sgem_i, sgem_j, sgem_c
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
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
!jdf
REAL, DIMENSION( ims:ime , jms:jme ) , &
INTENT(IN ) :: ebio_iso, ebio_olt, ebio_oli
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: rho_phy
!jdf
! 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(IN ) :: smois
! local variables
! anth_emiss_size_opt = 1 --> sizes follow those used in cam5
! anth_emiss_size_opt = 2 --> sizes follow those used in sorgam
integer, parameter :: anth_emiss_size_opt = 1 ! user can change this value
! dust_emiss_size_opt = 1 --> sizes follow those used in cam5
! dust_emiss_size_opt = 2 --> sizes follow those used in sorgam
integer, parameter :: dust_emiss_size_opt = 1 ! user can change this value
! seas_emiss_size_opt = 1 --> sizes follow those used in cam5
! seas_emiss_size_opt = 2 --> sizes follow those used in sorgam
integer, parameter :: seas_emiss_size_opt = 1 ! user can change this value
integer :: i, iphase, isize, itype
integer :: itype_accum, itype_aitken, itype_pcarbon, &
itype_finedust, itype_coardust
integer :: j, k
integer :: l, lemit, lemit_type, lnumb
integer :: m, n
integer :: p1st
real, parameter :: pi = 3.1415926536
real :: dp_anth_emiss_aitken, dp_anth_emiss_accum, dp_anth_emiss_coarse
real :: dlo_seas_emiss( maxd_asize, maxd_atype ), &
dhi_seas_emiss( maxd_asize, maxd_atype )
real :: fact_mass, factbb_mass, factbb_numb
real :: tmp_dens, tmp_dp, tmp_vol1p
real :: total_soag, fact_soag
character(len=100) :: msg
iphase = ai_phase
!
! code applies emissions in the emis_ant array,
! then does "online" sea-salt and dust emissions
!
! currently the code only works with
! config_flags%emiss_inpt_opt = emiss_inpt_default OR emiss_inpt_pnnl_rs
!
emiss_inpt_select_1: SELECT CASE( config_flags%emiss_inpt_opt )
!jdf
CASE( emiss_inpt_default, emiss_inpt_pnnl_rs, emiss_inpt_pnnl_mam )
write(*,*) 'cam_mam_addemiss DOING emiss - emiss_inpt_opt =', &
config_flags%emiss_inpt_opt
!jdf
CASE DEFAULT
write(*,*) 'cam_mam_addemiss SKIPPING emiss - emiss_inpt_opt =', &
config_flags%emiss_inpt_opt
return
END SELECT emiss_inpt_select_1
p1st = param_first_scalar
itype_accum = modeptr_accum
itype_aitken = modeptr_aitken
#if ( defined MODAL_AERO_3MODE )
itype_pcarbon = modeptr_accum
itype_finedust = modeptr_accum
itype_coardust = modeptr_coarse
#elif ( defined MODAL_AERO_7MODE )
if ( (config_flags%chem_opt /= CBMZ_CAM_MAM7_NOAQ) .and. &
(config_flags%chem_opt /= CBMZ_CAM_MAM7_AQ ) ) then
call wrf_error_fatal( 'emission package for MAM7 is not implemented yet' )
end if
#else
itype_pcarbon = modeptr_pcarbon
itype_finedust = modeptr_finedust
itype_coardust = modeptr_coardust
#endif
! set emitted sizes
if (anth_emiss_size_opt == 1) then
! volume-mean diameter for emitted particles (cm)
dp_anth_emiss_aitken = 0.0504e-4
dp_anth_emiss_accum = 0.134e-4
dp_anth_emiss_coarse = 2.06e-4
else if (anth_emiss_size_opt == 2) then
! volume-mean diameter for emitted particles (cm)
dp_anth_emiss_aitken = 1.0e2 * dgvem_i / exp( 1.5 * (log(sgem_i)**2) )
dp_anth_emiss_accum = 1.0e2 * dgvem_j / exp( 1.5 * (log(sgem_j)**2) )
dp_anth_emiss_coarse = 1.0e2 * dgvem_c / exp( 1.5 * (log(sgem_c)**2) )
else
write(msg,'(2a,i7)') 'subr cam_mam_addemiss', &
' - illegal anth_emiss_size_opt = ', anth_emiss_size_opt
call wrf_error_fatal( msg )
end if
!jdf add case statements here to switch between 2 emission file options
emiss_inpt_select_2: SELECT CASE( config_flags%emiss_inpt_opt )
CASE( emiss_inpt_pnnl_rs )
!jdf
do 1900 lemit_type = 1, 11
iphase = 1
isize = 1
! following if/then/else blocks determine
! which chem species will receive a particular emis_ant(...,p_e_...)
! additional mass factor for consistency with cam_mam assumptions
! particle size and density for converting mass to number emissions
!
! oc (org), ec, and pm25 - both the "i" and "j" go to the same
! cam_mam mode
! pm25 and pm_10 - go to the cam_mam dust species
! no3 - ignore cam_mam does not treat no3
!
if (lemit_type == 1) then
lemit = p_e_so4i
itype = itype_aitken
l = lptr_so4_aer(isize,itype,iphase)
tmp_dens = dens_so4_aer
factbb_mass = so4_emit_1p2_factor
tmp_dp = dp_anth_emiss_aitken
else if (lemit_type == 2) then
lemit = p_e_so4j
itype = itype_accum
l = lptr_so4_aer(isize,itype,iphase)
tmp_dens = dens_so4_aer
! so4_emit_1p2_factor is 115/96 if so4 (in 3-mode) is being treated
! as ammonium bisulfate, or 1.0 if treated as pure sulfate
factbb_mass = so4_emit_1p2_factor
tmp_dp = dp_anth_emiss_accum
else if (lemit_type == 3) then
lemit = p_e_no3i
cycle ! cam_mam does not treat no3
else if (lemit_type == 4) then
lemit = p_e_no3j
cycle ! cam_mam does not treat no3
else if (lemit_type == 5) then
lemit = p_e_orgi
itype = itype_pcarbon
l = lptr_pom_aer(isize,itype,iphase)
tmp_dens = dens_pom_aer
! pom_emit_1p4_factor is 1.4 if pom emissions are carbon-only emissions,
! or 1.0 if they are organic matter emissions
factbb_mass = pom_emit_1p4_factor
tmp_dp = dp_anth_emiss_accum
else if (lemit_type == 6) then
lemit = p_e_orgj
itype = itype_pcarbon
l = lptr_pom_aer(isize,itype,iphase)
tmp_dens = dens_pom_aer
factbb_mass = pom_emit_1p4_factor
tmp_dp = dp_anth_emiss_accum
else if (lemit_type == 7) then
lemit = p_e_eci
itype = itype_pcarbon
l = lptr_bc_aer(isize,itype,iphase)
tmp_dens = dens_bc_aer
factbb_mass = 1.0
tmp_dp = dp_anth_emiss_accum
else if (lemit_type == 8) then
lemit = p_e_ecj
itype = itype_pcarbon
l = lptr_bc_aer(isize,itype,iphase)
tmp_dens = dens_bc_aer
factbb_mass = 1.0
tmp_dp = dp_anth_emiss_accum
else if (lemit_type == 9) then
lemit = p_e_pm25i
itype = itype_finedust
l = lptr_dust_aer(isize,itype,iphase)
tmp_dens = dens_dust_aer
factbb_mass = 1.0
tmp_dp = dp_anth_emiss_accum
else if (lemit_type == 10) then
lemit = p_e_pm25j
itype = itype_finedust
l = lptr_dust_aer(isize,itype,iphase)
tmp_dens = dens_dust_aer
factbb_mass = 1.0
tmp_dp = dp_anth_emiss_accum
else if (lemit_type == 11) then
lemit = p_e_pm_10
itype = itype_coardust
l = lptr_dust_aer(isize,itype,iphase)
tmp_dens = dens_dust_aer
factbb_mass = 1.0
tmp_dp = dp_anth_emiss_coarse
else
cycle
end if
if ((l < p1st) .or. (l > num_chem)) cycle
lnumb = numptr_aer(isize,itype,iphase)
if ((lnumb < p1st) .or. (lnumb > num_chem)) lnumb = -999888777
tmp_vol1p = (pi/6.0)*(tmp_dp**3) ! mean volume of emitted particles (cm3)
! factbb_numb convert mass emissions (ug/m2/s) to number emissions (#/m2/s)
factbb_numb = 1.0e-6/(tmp_dens*tmp_vol1p)
do j = jts, jte
do k = 1, min(config_flags%kemit,kte)
do i = its, ite
! emissions fluxes are ug/m2/s
! mult by tmp_fact_mass to get mixing ratio change (ug/kg) over dtstep
fact_mass = (dtstep/dz8w(i,k,j))*alt(i,k,j)*factbb_mass
chem(i,k,j,l) = chem(i,k,j,l) + emis_ant(i,k,j,lemit)*fact_mass
if (lnumb > 0) then
chem(i,k,j,lnumb) = chem(i,k,j,lnumb) + &
emis_ant(i,k,j,lemit)*fact_mass*factbb_numb
end if
end do ! i
end do ! k
end do ! j
1900 continue
!
! now do emissions for e_soag
! add biogenic isoprene for soag_isoprene
! no monoterpenes in CBMZ at present for soag_terpene
! lump anthropogenic ald, hc3, hc5, hc8, ket, oli, olt, and ora2 which is par in CBMZ for soag_bigalk
! lump biogenic olet and olei for soag_bigene
! lump anthropogenic tol and xyl for soag_toluene
!
do j = jts, jte
do k = 1, 1
do i = its, ite
fact_soag = 4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
total_soag = (emis_ant(i,k,j,p_e_tol) * (92.14*0.15/12.0)) + &
(emis_ant(i,k,j,p_e_xyl) * (106.16*0.15/12.0)) + &
(emis_ant(i,k,j,p_e_hc5) * (77.6*0.05/12.0)) + &
(emis_ant(i,k,j,p_e_olt) * (72.34*0.05/12.0)) + &
(emis_ant(i,k,j,p_e_oli) * (75.78*0.05/12.0)) + &
(ebio_iso(i,j) * (68.11*0.04/12.0))
! 0.4*emis_ant(i,k,j,p_e_ald) + 2.9*emis_ant(i,k,j,p_e_hc3) + &
! 4.8*emis_ant(i,k,j,p_e_hc5) + 7.9*emis_ant(i,k,j,p_e_hc8) + &
! 0.9*emis_ant(i,k,j,p_e_ket) + 2.8*emis_ant(i,k,j,p_e_oli) + &
! 1.8*emis_ant(i,k,j,p_e_olt) + 1.0*emis_ant(i,k,j,p_e_ora2) + &
! (ebio_iso(i,j) * (68.11*0.04/12.0)) + &
! ebio_olt(i,j) + &
! ebio_oli(i,j)
chem(i,k,j,p_soag) = chem(i,k,j,p_soag) + total_soag*fact_soag
end do ! i
end do ! k
end do ! j
do j = jts, jte
do k = 2, min(config_flags%kemit,kte)
do i = its, ite
fact_soag = 4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
total_soag = (emis_ant(i,k,j,p_e_tol) * (92.14*0.15/12.0)) + &
(emis_ant(i,k,j,p_e_xyl) * (106.16*0.15/12.0)) + &
(emis_ant(i,k,j,p_e_hc5) * (77.6*0.05/12.0)) + &
(emis_ant(i,k,j,p_e_olt) * (72.34*0.05/12.0)) + &
(emis_ant(i,k,j,p_e_oli) * (75.78*0.05/12.0))
! emis_ant(i,k,j,p_e_xyl) + &
! 0.4*emis_ant(i,k,j,p_e_ald) + 2.9*emis_ant(i,k,j,p_e_hc3) + &
! 4.8*emis_ant(i,k,j,p_e_hc5) + 7.9*emis_ant(i,k,j,p_e_hc8) + &
! 0.9*emis_ant(i,k,j,p_e_ket) + 2.8*emis_ant(i,k,j,p_e_oli) + &
! 1.8*emis_ant(i,k,j,p_e_olt) + 1.0*emis_ant(i,k,j,p_e_ora2)
chem(i,k,j,p_soag) = chem(i,k,j,p_soag) + total_soag*fact_soag
end do ! i
end do ! k
end do ! j
CASE( emiss_inpt_pnnl_mam )
do 1910 lemit_type = 1, 14
iphase = 1
isize = 1
if (lemit_type == 1) then
lemit = p_e_so4i
itype = itype_aitken
l = lptr_so4_aer(isize,itype,iphase)
factbb_mass = so4_emit_1p2_factor
else if (lemit_type == 2) then
lemit = p_e_so4j
itype = itype_accum
l = lptr_so4_aer(isize,itype,iphase)
! so4_emit_1p2_factor is 115/96 if so4 (in 3-mode) is being treated
! as ammonium bisulfate, or 1.0 if treated as pure sulfate
factbb_mass = so4_emit_1p2_factor
else if (lemit_type == 3) then
lemit = p_e_orgj
itype = itype_pcarbon
l = lptr_pom_aer(isize,itype,iphase)
factbb_mass = pom_emit_1p4_factor
else if (lemit_type == 4) then
lemit = p_e_ecj
itype = itype_pcarbon
l = lptr_bc_aer(isize,itype,iphase)
factbb_mass = 1.0
else if (lemit_type == 5) then
lemit = p_e_dust_a1
itype = itype_finedust
l = lptr_dust_aer(isize,itype,iphase)
factbb_mass = 1.0
else if (lemit_type == 6) then
lemit = p_e_dust_a3
itype = itype_coardust
l = lptr_dust_aer(isize,itype,iphase)
factbb_mass = 1.0
else if (lemit_type == 7) then
lemit = p_e_ncl_a1
itype = itype_accum
l = lptr_seas_aer(isize,itype,iphase)
factbb_mass = 1.0
else if (lemit_type == 8) then
lemit = p_e_ncl_a2
itype = itype_aitken
l = lptr_seas_aer(isize,itype,iphase)
factbb_mass = 1.0
else if (lemit_type == 9) then
lemit = p_e_ncl_a3
itype = itype_coardust
l = lptr_seas_aer(isize,itype,iphase)
factbb_mass = 1.0
!PMA
!Currently my pre-processing script put num_a1 and num_a2 emissions from
!sea salt and dust into so4j_num and so4i_num
!num_a3 emissions are solely from dust and sea salt
!
else if (lemit_type == 10) then
lemit = p_e_so4i_num
itype = itype_aitken
l = numptr_aer(isize,itype,iphase)
else if (lemit_type == 11) then
lemit = p_e_so4j_num
itype = itype_accum
l = numptr_aer(isize,itype,iphase)
else if (lemit_type == 12) then
lemit = p_e_orgj_num
itype = itype_accum
l = numptr_aer(isize,itype,iphase)
else if (lemit_type == 13) then
lemit = p_e_ecj_num
itype = itype_accum
l = numptr_aer(isize,itype,iphase)
else if (lemit_type == 14) then
lemit = p_e_num_a3
itype = itype_coardust
l = numptr_aer(isize,itype,iphase)
else
cycle
end if
if ((l < p1st) .or. (l > num_chem)) cycle
do j = jts, jte
do k = 1, min(config_flags%kemit,kte)
do i = its, ite
! emissions fluxes are ug/m2/s
! mult by tmp_fact_mass to get mixing ratio change (ug/kg) over dtstep
fact_mass = (dtstep/dz8w(i,k,j))*alt(i,k,j)*factbb_mass
!PMA
factbb_numb = (dtstep/dz8w(i,k,j))*alt(i,k,j)
if (lemit_type < 10) then
chem(i,k,j,l) = chem(i,k,j,l) + emis_ant(i,k,j,lemit)*fact_mass
else
chem(i,k,j,l) = chem(i,k,j,l) + emis_ant(i,k,j,lemit)*factbb_numb
end if
end do ! i
end do ! k
end do ! j
1910 continue
!
! units for e_soag* assumed to be ug/m2/s for now
!
do j = jts, jte
do k = 1, min(config_flags%kemit,kte)
do i = its, ite
fact_soag = (dtstep/dz8w(i,k,j))*alt(i,k,j)/1000.0
total_soag = emis_ant(i,k,j,p_e_soag_bigalk) + &
emis_ant(i,k,j,p_e_soag_bigene) + &
emis_ant(i,k,j,p_e_soag_isoprene) + &
emis_ant(i,k,j,p_e_soag_terpene) + &
emis_ant(i,k,j,p_e_soag_toluene)
chem(i,k,j,p_soag) = chem(i,k,j,p_soag) + total_soag*fact_soag
end do ! i
end do ! k
end do ! j
END SELECT emiss_inpt_select_2
!jdf end of add case statements here to switch between 2 emission file options
!
! seasalt emissions
!
! set emitted sizes
dlo_seas_emiss(:,:) = 0.0
dhi_seas_emiss(:,:) = 0.0
if (seas_emiss_size_opt == 1) then
#if ( defined MODAL_AERO_3MODE )
! cut sizes are 0.02, 0.08, 1.0, 10.0 um
dlo_seas_emiss(1,modeptr_aitken ) = 0.02e-4
dhi_seas_emiss(1,modeptr_aitken ) = 0.08e-4
dlo_seas_emiss(1,modeptr_accum ) = 0.08e-4
dhi_seas_emiss(1,modeptr_accum ) = 1.00e-4
dlo_seas_emiss(1,modeptr_coarse ) = 1.00e-4
dhi_seas_emiss(1,modeptr_coarse ) = 10.0e-4
#else
! cut sizes are 0.02, 0.08, 0.3, 1.0, 10.0 um
dlo_seas_emiss(1,modeptr_aitken ) = 0.02e-4
dhi_seas_emiss(1,modeptr_aitken ) = 0.08e-4
dlo_seas_emiss(1,modeptr_accum ) = 0.08e-4
dhi_seas_emiss(1,modeptr_accum ) = 0.30e-4
dlo_seas_emiss(1,modeptr_fineseas) = 0.30e-4
dhi_seas_emiss(1,modeptr_fineseas) = 1.00e-4
dlo_seas_emiss(1,modeptr_coarseas) = 1.00e-4
dhi_seas_emiss(1,modeptr_coarseas) = 10.0e-4
#endif
else if (seas_emiss_size_opt == 2) then
#if ( defined MODAL_AERO_3MODE )
! cut sizes are 0.1, 1.25, 10.0 um and no aitken
dlo_seas_emiss(1,modeptr_accum ) = 0.10e-4
dhi_seas_emiss(1,modeptr_accum ) = 1.25e-4
dlo_seas_emiss(1,modeptr_coarse ) = 1.25e-4
dhi_seas_emiss(1,modeptr_coarse ) = 10.0e-4
#else
! cut sizes are 0.1, 0.3, 1.25, 10.0 um and no aitken
dlo_seas_emiss(1,modeptr_accum ) = 0.10e-4
dhi_seas_emiss(1,modeptr_accum ) = 0.30e-4
dlo_seas_emiss(1,modeptr_fineseas) = 0.30e-4
dhi_seas_emiss(1,modeptr_fineseas) = 1.25e-4
dlo_seas_emiss(1,modeptr_coarseas) = 1.25e-4
dhi_seas_emiss(1,modeptr_coarseas) = 10.0e-4
#endif
else
write(msg,'(2a,i7)') 'subr cam_mam_addemiss', &
' - illegal seas_emiss_size_opt = ', seas_emiss_size_opt
call wrf_error_fatal( msg )
end if
! now do the sea-salt emissions
if (config_flags%seas_opt == 2) then
call cam_mam_mosaic_seasalt_emiss( &
id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
dlo_seas_emiss, dhi_seas_emiss, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
end if
!
! dust emissions
!
! jdf: preliminary version that has not been made generic for situation
if (config_flags%dust_opt == 2) then
call cam_mam_mosaic_dust_emiss( &
slai, ust, smois, ivgtyp, isltyp, &
id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
dust_emiss_size_opt, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
end if
return
END subroutine cam_mam_addemiss
!----------------------------------------------------------------------
subroutine cam_mam_mosaic_seasalt_emiss( &
id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
dlo_seas_emiss, dhi_seas_emiss, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!
! 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_cam_mam_asect, only: ai_phase, &
lptr_seas_aer, &
maxd_asize, maxd_atype, nsize_aer, ntype_aer, numptr_aer
! USE module_mosaic_addemiss, only: seasalt_emitfactors_1bin
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
! 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
! emissions cut sizes (cm)
REAL, DIMENSION( maxd_asize, maxd_atype ), &
INTENT(IN ) :: dlo_seas_emiss, dhi_seas_emiss
! local variables
integer i, j, k, l, l_seas, n
integer iphase, itype
integer p1st
real dum, dumdlo, dumdhi, dumoceanfrac, dumspd10
real factaa, factbb
real :: ssemfact_numb( maxd_asize, maxd_atype )
real :: ssemfact_mass( maxd_asize, maxd_atype )
write(*,*) 'in subr cam_mam_mosaic_seasalt_emiss'
p1st = PARAM_FIRST_SCALAR
! for now just do itype=1
iphase = ai_phase
! compute emissions factors for each size bin
! (limit emissions to dp > 0.1 micrometer)
ssemfact_mass(:,:) = 0.0
ssemfact_numb(:,:) = 0.0
do itype = 1, ntype_aer
do n = 1, nsize_aer(itype)
dumdlo = max( dlo_seas_emiss(n,itype), 0.1e-4 )
dumdhi = max( dhi_seas_emiss(n,itype), 0.1e-4 )
if (dumdlo >= dumdhi) cycle
call seasalt_emitfactors_1bin( 1, dumdlo, dumdhi, &
ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype) )
! convert mass emissions factor from (g/m2/s) to (ug/m2/s)
ssemfact_mass(n,itype) = ssemfact_mass(n,itype)*1.0e6
end do
end do
! loop over i,j and apply seasalt emissions
k = kts
do 1830 j = jts, jte
do 1820 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
do 1815 itype = 1, ntype_aer
do 1810 n = 1, nsize_aer(itype)
if (ssemfact_mass(n,itype) <= 0.0) cycle
! only apply emissions if bin has both na and cl species
l_seas = lptr_seas_aer(n,itype,iphase)
if (l_seas < p1st) cycle
chem(i,k,j,l_seas) = chem(i,k,j,l_seas) + &
factbb * ssemfact_mass(n,itype)
l = numptr_aer(n,itype,iphase)
if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
factbb * ssemfact_numb(n,itype)
1810 continue
1815 continue
1820 continue
1830 continue
return
END subroutine cam_mam_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
!----------------------------------------------------------------------
subroutine cam_mam_mosaic_dust_emiss( slai,ust, smois, ivgtyp, isltyp, &
id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
dust_emiss_size_opt, &
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).
!
USE module_configure, only: grid_config_rec_type
USE module_state_description, only: num_chem, param_first_scalar
USE module_data_cam_mam_asect, only: ai_phase, &
lptr_dust_aer, ntype_aer, numptr_aer
USE modal_aero_data, only: &
modeptr_accum, modeptr_coarse, modeptr_coardust, modeptr_finedust
IMPLICIT NONE
TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
INTEGER, INTENT(IN ) :: id, &
dust_emiss_size_opt, &
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(IN ) :: smois
! local variables
integer, parameter :: max_dust_mode = 2
integer :: i, ii, j, k, l, l_dust, l_numb, n
integer :: imode, iphase, isize, itype, izob
integer :: p1st
integer :: isize_dust_mode(max_dust_mode)
integer :: itype_dust_mode(max_dust_mode)
integer :: n_dust_mode
real :: dum, dumdlo, dumdhi, dumlandfrac, dumspd10
real :: factaa, factbb, factcc, 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
real :: mass1part_mos8bin(8)
real :: sz_wght_dust_mode(8,max_dust_mode)
real :: dcen_mos8bin(8)
character(len=100) :: msg
write(*,*) 'in subr cam_mam_mosaic_dust_emiss'
p1st = param_first_scalar
! the mass emission fractions [sz(1:8)] were derived for the
! mosaic 8-bin size bins
! the emissions for each mosaic bin are assigned to the cam-mam fine or coarse mode
! (e.g, the fine mode may get bins 1-5 and the coarse mode bins 6-8)
! for calculating number emissions from mass emissions, use the mosaic bin
! sizes (and corresponding 1-particle masses) in order to get number emissions
! consistent with mosaic
!
! note: the sz() values were recommend by Natalie Mahowold, so probably follow the
! assumptions used in the "dead" dust module.
! applying those assumptions along with the cut sizes might be a better approach
! than adapting the mosaic 8-bin sz() values
#if (defined MODAL_AERO)
! 3 modes -- dust is in accum and coarse modes
n_dust_mode = 2
itype_dust_mode(1) = modeptr_accum
itype_dust_mode(2) = modeptr_coarse
isize_dust_mode(:) = 1
if (dust_emiss_size_opt == 1) then
! cam5 cut sizes for dust emiss are 0.1, 1.0, 10.0 um
! fine mode gets bins 1-4 and 60% of bin 5
sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 0.6, 0., 0., 0. /)
sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0.4, 1., 1., 1. /)
else if (dust_emiss_size_opt == 2) then
! sorgam cut sizes for dust emiss are 0.1, 2.5, 10.0 um
! fine mode gets bins 1-6
sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 1., 1., 0., 0. /)
sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0., 0., 1., 1. /)
else
write(msg,'(2a,i7)') 'subr cam_mam_mosaic_dust_emiss', &
' - illegal dust_emiss_size_opt = ', dust_emiss_size_opt
call wrf_error_fatal( msg )
end if
#else
! 7 modes -- dust is in fine-dust and coarse-dust modes
n_dust_mode = 2
itype_dust_mode(1) = modeptr_finedust
itype_dust_mode(2) = modeptr_coardust
isize_dust_mode(:) = 1
if (dust_emiss_size_opt == 1) then
! cam5 cut sizes for dust emiss are 0.1, 2.0, 10.0 um
! fine mode gets bins 1-5 and 60% of bin 6
! (for 7-mode cam, the fine-dust mode is larger than the accumulation mode,
! so the cut size is larger than with 3-mode cam)
sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 1., 0.6, 0., 0. /)
sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0., 0.4, 1., 1. /)
else if (dust_emiss_size_opt == 2) then
! sorgam cut sizes for dust emiss are 0.1, 2.5, 10.0 um
! fine mode gets bins 1-6
sz_wght_dust_mode(1:8,1) = (/ 1., 1., 1., 1., 1., 1., 0., 0. /)
sz_wght_dust_mode(1:8,2) = (/ 0., 0., 0., 0., 0., 0., 1., 1. /)
else
write(msg,'(2a,i7)') 'subr cam_mam_mosaic_dust_emiss', &
'illegal dust_emiss_size_opt = ', dust_emiss_size_opt
call wrf_error_fatal( msg )
end if
#endif
! bin center diameters for mosaic-8bin (cm)
dcen_mos8bin(8) = 1.0e-4*(10.0/sqrt(2.0))
do n = 7, 1, -1
dcen_mos8bin(n) = dcen_mos8bin(n+1)*0.5
end do
! mass1part_mos8bin is mass of a single particle in ug,
! assuming density of dust ~2.5 g cm-3
densdust=2.5
do n = 1, 8
mass1part_mos8bin(n)=0.523598*(dcen_mos8bin(n)**3)*densdust*1.0e06
end do
!
! 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
!jdf sz(1)=0.0
!jdf sz(2)=0.0
!jdf sz(3)=0.0005
!jdf sz(4)=0.0095
!jdf sz(5)=0.03
!jdf sz(6)=0.10
!jdf sz(7)=0.18
!jdf sz(8)=0.68
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
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 1820
!
! 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
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=1.0e-14*(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
fractot = 1.0
do imode = 1, n_dust_mode
! loop over the cam-mam modes that have dust
itype = itype_dust_mode(imode)
isize = isize_dust_mode(imode)
l_dust = lptr_dust_aer(isize,itype,iphase)
if ((l_dust < p1st) .or. (l_dust > num_chem)) cycle
l_numb = numptr_aer(isize,itype,iphase)
if ((l_numb < p1st) .or. (l_numb > num_chem)) l_numb = -1
! if (ivgtyp(i,j) .eq. 8) print*,'jdf',i,j,ustar1,ustart0,factaa,ftot
! do 1810 n = 1, nsize_aer(itype)
do n = 1, 8
! calculate contribution of each mosaic-8bin bin to the current cam-mam mode
if (sz_wght_dust_mode(n,imode) <= 0.0) cycle
factcc = factbb * sz(n) * fractot * sz_wght_dust_mode(n,imode)
chem(i,k,j,l_dust) = chem(i,k,j,l_dust) + factcc
if (l_numb >= p1st) &
chem(i,k,j,l_numb) = chem(i,k,j,l_numb) + factcc/mass1part_mos8bin(n)
end do ! n
end do ! imode
1820 continue
1830 continue
return
END subroutine cam_mam_mosaic_dust_emiss
!----------------------------------------------------------------------
END MODULE module_cam_mam_addemiss