We use WRF/Chem 3.4.1 with modifications. Our setup options are given in the overview table here. Output variables are listed here.
Table of Contents
|
Aerosol module
MOSAIC with 4 size bins, simple SOA scheme.
Bin size ranges:
(1) 3 nm - 156 nm
(2) 156 nm - 625 nm
(3) 625 nm - 2.5 μm
(4) 2.5 μm - 10 μm
Mass contributions per bin: bc, biog1_o (biogenic SOA from isoprene), biog1_c (biogenic SOA from pinenes), ca, cl, co3, glysoa (SOA from glyoxal), na, nh4, no3, oc, oin, smpa (simplified anthropogenic SOA), smpbb(simplified biomass burning SOA), so4, water
Conversion of MACC aerosol boundary conditions
The following script creates files that can be read by the mozbc script on the ACD website:
; NAME:
; MACC2_MOZART_AQMEII.ncl
;
; PURPOSE:
; Convert IFS-MOZART output into the same format as provided by
; NCAR MOZART files.
; This conversion is the first step of a two-step pre-processing
; of IFS-MOZART files to be compatible with the format required
; by the COSMO-ART preprocessor int2lm.
;
; The three input files PS_, GRG_, AER_ containing surface pressure,
; Global Reactive Gases and AERosols are merged into one, adds a variable
; time with units of days since 0000-01-01 00:00:00 as "gregorian"
; calender, a variable date (6 digit integer) a variable datesec
; (seconds of day) and a reference pressure P0 which is set to 1.0 here
; for reasons mentioned in a comment below.
;
; Note: although the date is given as "gregorian" calender, it is actually
; a "proleptic_gregorian" calender. This wrong annotation is made to be
; compatible with the MOZART-NCEP output.
;
; CALLING:
; ncl fileid="fileid" indir="indir" outdir="outdir" MACC2_MOZART_AQMEII.ncl
;
; INPUTS:
; fileid: the date to process in format YYYYMMDD, e.g. 20100213
; indir: directory of IFS-MOZART files (default is /nas/input/IFS-MOZART)
; outdir: directory of processed file (default is /nas/input/IFS-MOZART/processed)
;
; AUTHOR:
; Christoph Knote
; v1: CK, Jan 2013: first implementation
; v2: DB, 15 Jan 2013: added possibility to provide command line input
;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
if (.not. isvar("fileid")) then ; is fileid on command line?
print("parameter fileid missing")
exit
end if
if (.not. isvar("indir")) then
indir = "/glade/p/acd/knote/AQMEII/MACC2_MOZART"
end if
if (.not. isvar("outdir")) then
outdir = "/glade/p/acd/knote/AQMEII/MACC2_MOZART/processed"
end if
; print("fileid = "+fileid)
; print("indir = "+indir)
; print("outdir = "+outdir)
; input files
psfile = indir+"/PS_"+fileid+"_NA_AQ.nc"
grgfile = indir+"/GRG_"+fileid+"_NA_AQ.nc"
aerfile = indir+"/AER_"+fileid+"_NA_AQ.nc"
if ( isfilepresent(psfile) .and. isfilepresent(grgfile) .and. isfilepresent(aerfile) ) then
ps = addfile(psfile,"r")
grg = addfile(grgfile,"r")
aer = addfile(aerfile,"r")
; output file
system ("rm "+outdir+"/AQMEII_"+fileid+".nc") ; remove pre-existing file
out = addfile(outdir+"/AQMEII_"+fileid+".nc","c")
; coordinate variables
; levels are cut out from level 7 to 60 (from top)
hyam = ps->hyam(6:59)
hybm = ps->hybm(6:59)
hyam!0 = "lev"
hybm!0 = "lev"
hyam&lev = ispan(0,53,1)
hybm&lev = ispan(0,53,1)
; compatible with MOZART output
hyam = hyam / 100000.0
; date and time
date = toint(cd_calendar(ps->time, -2))
date@long_name = "current date as 6 digit integer (YYMMDD)"
date!0 = "time"
date&time = ps->time&time
datesec1 = cd_calendar(ps->time, 0)
datesec = toint(datesec1(:,3) * 3600 + datesec1(:,4) * 60 + datesec1(:,5))
datesec@units = "s"
datesec@long_name = "seconds to complete current date"
datesec!0 = "time"
datesec&time = ps->time&time
time = cd_convert(ps->time, "days since 0000-01-01 00:00:00")
time@long_name = "simulation time"
time@calendar = "gregorian"
time!0 = "time"
time&time = ps->time&time
out->time = time
out->hyam = hyam
out->hybm = hybm
out->lat = ps->lat
out->lon = ps->lon
out->date = date
out->datesec = datesec
; surface pressure
pres1 = ps->var152
pres = pres1(:,0,:,:)
pres@units = "PA"
out->PS = pres
p0 = 100000
p0@long_name = "reference pressure"
p0@units = "Pa"
out->P0 = p0
; gases
gases = getfilevarnames(grg)
ngases = dimsizes(gases)
do n=0,ngases-1
gas = grg->$gases(n)$
dimg = dimsizes(gas)
rank = dimsizes(dimg)
if (rank .eq. 4) then
gas@units = "VMR"
gas&lev = (/ ispan(0,53,1) /)
out->$gases(n)$ = gas
end if
delete(dimg)
delete(rank)
delete(gas)
end do
; aerosols
;aeros = (/ "var1", "var2", "var3", "var4", "var5", "var6", "var8","var10","var11" /)
;aerosNice = (/ "SS01", "SS02", "SS03", "DU01", "DU02", "DU03", "OC", "BC", "SO4" /)
; DB: sea salt ignored due to strong biases
aeros = (/ "var4", "var5", "var6", "var8","var10","var11" /)
aerosNice = (/ "DU01", "DU02", "DU03", "OC", "BC", "SO4" /)
naeros = dimsizes(aeros)
do n=0,naeros-1
aero = aer->$aeros(n)$
aero@units = "UG M-3"
aero&lev = ispan(0,53,1)
out->$aerosNice(n)$ = aero
delete(aero)
end do
else
print("Not all input files available for fileid "+fileid)
end if
Emissions preprocessing
We use model ready files as provided by EPA on their FTP site. Here is the NCL script we use to convert emissions to wrfchemi* files:
;
; AQMEII project emissions converter
; Christoph Knote, NCAR, ACD, 03/28/2013
;
; Description:
; Conversion of EPA provided model ready emissions files to wrfchemi* input.
; Area emissions are aggregated to model resolution - only works for
; model resolutions coarser than ~15 km. Plume rise for point sources
; based on constant meteorology - no consideration of actual conditions!
; Use at your own risk. Please report bugs found to knote@ucar.edu.
;
; Call:
; ncl date=YYYYmmdd emiss_preproc.ncl
; creates emissions for a whole day (i.e. 24 wrfchemi files in ./output)
;
; Directory structure expected:
; anc - ancillary files
; * put a wrfinput of your target domain in there - called "wrfdomain.nc"
; * needs to contain the NEI_grid.nc file that describes the location of the input data
; inln - inln files from EPA ftp site
; layer1 - layer1 files from EPA ftp site
; output - well, it's what it says.
; stack_groups - stack_groups files from EPA ftp site
; tmp - temporary files will be created in here.
;
; Notes:
; First call will take a long time, as two kinds of cached data are created:
; tmp/NEI_grid_rel.nc - correspondence between input grid point and output grid point
; tmp/stacklocations_*.nc - locations of all stacks in output grid
; Subsequent calls will be much faster...
;
; License:
; use and modify at your own discretion, just keep header and add revisions.
;
; Dependencies:
; ncl >= 6.1
;
; version date creator
; --------------------------------
; 1.0 20130328 Christoph Knote
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
; date is YYYYmmdd
if (.not. isvar("date")) then ; is date on command line?
print("parameter date (YYYYmmdd) missing")
exit
end if
print(" ")
print("---")
print(tostring(date))
print(" ")
vars = (/ "ACROLEIN","ALD2_PRIMARY","ALD2","ALDX","BENZENE","BUTADIENE13","CH4","CL2","CO","ETHA","ETH","ETOH","FORM_PRIMARY","FORM","HCL","HONO","IOLE","ISOP","MEOH","NH3_FERT","NH3","NO2","NO","NVOL","OLE","SO2","SULF","TERP","TOL","UNK","UNR","VOC_INV","XYL","PAL","PAR","PCA","PCL","PEC","PFE","PH2O","PK","PMC","PMFINE","PMG","PMN","PMOTHR","PNA","PNCOM","PNH4","PNO3","POC","PSI","PSO4","PTI","CO_A","CO_BB" /)
ntime = 25
nvar = dimsizes(vars)
layer1FileRaw = "layer1/emis_mole_all_"+date+"_12US1_cmaq_cb05_soa_E21_2007ed_10.ncf.gz"
layer1File = "tmp/layer1_"+date+"_tmp.nc"
; fires (*ptfire* files) have time varying position - fireidx is the index in the inlnFiles / stackFiles array
fireidx=0
inlnFilesRaw = (/ "inln/inln_mole_ptfire_"+date+"_12US1_cmaq_cb05_soa_E21_2007ed_10.ncf.gz", \
"inln/inln_mole_c3marine_"+date+"_12US1_cmaq_cb05_soa_E21_2007ed_10.ncf.gz", \
"inln/inln_mole_canpt_"+date+"_12US1_cmaq_cb05_soa_E21_2007ed_10.ncf.gz", \
"inln/inln_mole_mexpt_"+date+"_12US1_cmaq_cb05_soa_E21_2007ed_10.ncf.gz", \
"inln/inln_mole_ptipm_"+date+"_12US1_cmaq_cb05_soa_E21_2007ed_10.ncf.gz", \
"inln/inln_mole_ptnonipm_"+date+"_12US1_cmaq_cb05_soa_E21_2007ed_10.ncf.gz" /)
; inlnFilesRaw and inlnFiles have to have the same ordering!
inlnFiles = (/ "tmp/inln_fire_"+date+"_tmp.nc", \
"tmp/inln_c3marine_"+date+"_tmp.nc", \
"tmp/inln_canpt_"+date+"_tmp.nc", \
"tmp/inln_mexpt_"+date+"_tmp.nc", \
"tmp/inln_ptipm_"+date+"_tmp.nc", \
"tmp/inln_ptnonipm_"+date+"_tmp.nc" /)
; note that all stack files are time invariant, except for fires (because stacks don't move...)
stackFilesRaw = (/ "stack_groups/stack_groups_ptfire_"+date+"_12US1_E21_2007ed_10.ncf.gz", \
"stack_groups/stack_groups_c3marine_12US1_E21_2007ed_10.ncf.gz", \
"stack_groups/stack_groups_canpt_12US1_E21_2007ed_10.ncf.gz", \
"stack_groups/stack_groups_mexpt_12US1_E21_2007ed_10.ncf.gz", \
"stack_groups/stack_groups_ptipm_12US1_E21_2007ed_10.ncf.gz", \
"stack_groups/stack_groups_ptnonipm_12US1_E21_2007ed_10.ncf.gz" /)
; inlnFilesRaw and inlnFiles have to have the same ordering!
stackFiles = (/ "tmp/stack_groups_ptfire_"+date+"_tmp.nc", \
"tmp/stack_groups_c3marine_"+date+"_tmp.nc", \
"tmp/stack_groups_canpt_"+date+"_tmp.nc", \
"tmp/stack_groups_mexpt_"+date+"_tmp.nc", \
"tmp/stack_groups_ptipm_"+date+"_tmp.nc", \
"tmp/stack_groups_ptnonipm_"+date+"_tmp.nc" /)
; safe guard
if (.not. isfilepresent(layer1FileRaw)) then
print("Input file "+layer1FileRaw+" missing")
exit
end if
; -----------------------------------------------------------------------------
; input data unzipping
; -----------------------------------------------------------------------------
print("unzipping input data")
system("gunzip -c "+layer1FileRaw+" > "+layer1File)
npointFiles = dimsizes(inlnFiles)
do i = 0, npointFiles - 1
if (.not. isfilepresent(inlnFiles(i)))
system("gunzip -c "+inlnFilesRaw(i)+" > "+inlnFiles(i))
end if
if (.not. isfilepresent(stackFiles(i)))
system("gunzip -c "+stackFilesRaw(i)+" > "+stackFiles(i))
end if
end do
wrfgrid = addfile("anc/wrfdomain.nc","r")
wrflat = wrfgrid->XLAT(0,:,:)
wrflon = wrfgrid->XLONG(0,:,:)
wrflonmin = min(wrflon)
wrflonmax = max(wrflon)
wrflatmin = min(wrflat)
wrflatmax = max(wrflat)
wrfdim = dimsizes(wrflat)
nrow_wrf = wrfdim(0)
ncol_wrf = wrfdim(1)
WRF_grid_area = (wrfgrid@DX * wrfgrid@DY) / 1e6 ; km^2
; -----------------------------------------------------------------------------
; area emissions
; -----------------------------------------------------------------------------
print("area emissions interpolation")
nrow_nei = 299
ncol_nei = 459
if (.not. isfilepresent("tmp/NEI_grid_rel.nc"))
; make a frame of points around the target grid to catch contributions from
; surrounding grid cells
ext = 5
wrflonext = new( (/ nrow_wrf+2*ext, ncol_wrf+2*ext /), "float", -999.0)
wrflatext = new( (/ nrow_wrf+2*ext, ncol_wrf+2*ext /), "float", -999.0)
wrflonext(ext:(nrow_wrf+ext-1),ext:(ncol_wrf+ext-1)) = (/ wrflon /)
wrflatext(ext:(nrow_wrf+ext-1),ext:(ncol_wrf+ext-1)) = (/ wrflat /)
wrfdimext = dimsizes(wrflonext)
; just a guess, doesn't have to be accurate
adx = wrfgrid@DX/40000000*360*0.6
ady = wrfgrid@DY/40000000*360
do i = 1, ext
wrflonext(i-1,:) = wrflonext(ext,:)
wrflonext(nrow_wrf-1+ext+i,:) = wrflonext(nrow_wrf-1+ext,:)
wrflonext(:,ext-i) = wrflonext(:,ext) - i * adx
wrflonext(:,ncol_wrf+ext-1+i) = wrflonext(:,ncol_wrf+ext-1) + i * adx
wrflatext(:,i-1) = wrflatext(:,ext)
wrflatext(:,ncol_wrf-1+ext+i) = wrflatext(:,ncol_wrf-1+ext)
wrflatext(ext-i,:) = wrflatext(ext,:) - i * adx
wrflatext(nrow_wrf+ext-1+i,:) = wrflatext(nrow_wrf+ext-1,:) + i * adx
end do
neigrid = addfile("anc/NEI_grid.nc","r")
neicelat = neigrid->CElat
neicelon = neigrid->CElon
grid_rel = new( (/ nrow_nei, ncol_nei, 2 /), "integer", "No_FillValue")
grid_rel(:,:,:) = -1
do i = 0, nrow_nei - 1
do j = 0, ncol_nei - 1
if ( (neicelon(i,j) .ge. min(wrflonext) .and. neicelon(i,j) .le. max(wrflonext)) .and. (neicelat(i,j) .ge. min(wrflatext) .and. neicelat(i,j) .le. max(wrflatext)))
dlon = abs(neicelon(i,j) - wrflonext)
dlat = abs(neicelat(i,j) - wrflatext)
dtot = dlon + dlat
dtot1D = ndtooned(dtot)
indices = ind_resolve(minind(dtot1D),wrfdimext)
if ( (indices(0,0) .ge. ext .and. indices(0,0) .lt. nrow_wrf+ext-1) .and. (indices(0,1) .ge. ext .and. indices(0,1) .lt. ncol_wrf+ext-1) )
print(i+"x"+j+" is at "+(indices(0,:)-ext)+" is good")
grid_rel(i,j,:) = indices(0,:)-ext
else
print(i+"x"+j+" is at "+(indices(0,:)-ext)+" is bad")
end if
end if
end do
end do
grid_rel_file = addfile("tmp/NEI_grid_rel.nc","c")
grid_rel_file->grid_rel = grid_rel
delete(grid_rel_file)
delete(grid_rel)
end if
grid_rel_file = addfile("tmp/NEI_grid_rel.nc","r")
grid_rel = grid_rel_file->grid_rel
delete(grid_rel_file)
; interpolation:
area = new( (/ ntime, nrow_wrf, ncol_wrf, nvar /), "float", "No_FillValue" )
area!0 = "TSTEP"
area!1 = "ROW"
area!2 = "COL"
area!3 = "VAR"
in = addfile(layer1File, "r")
do l = 0, nvar - 1
if (isfilevar(in, vars(l))) then
dat = in->$vars(l)$(:,0,:,:)
do j = 0, ncol_nei - 1
do i = 0, nrow_nei - 1
if ( all(grid_rel(i,j,:) .gt. -1) )
ii = grid_rel(i,j,0)
jj = grid_rel(i,j,1)
area(:,ii,jj,l) = area(:,ii,jj,l) + dat(:,i,j)
end if
end do
end do
delete(dat)
print("done with "+vars(l))
end if
end do
; all area CO is anthropogenic CO
area(:,:,:,ind(vars .eq. "CO_A")) = area(:,:,:,ind(vars .eq. "CO"))
delete(ii)
delete(jj)
; fill missing values with zeros
area = where(ismissing(area), 0.0, area)
; -----------------------------------------------------------------------------
; point emissions
; -----------------------------------------------------------------------------
print("point emissions gridding")
; calculate mean WRF level heights
; maximum emission height considered
emis_height_max = 1100.
nlay = wrfgrid@$"BOTTOM-TOP_GRID_DIMENSION"$
ph = wrfgrid->PH(0,:,:,:)
phb = wrfgrid->PHB(0,:,:,:)
hgt = wrfgrid->HGT(0,:,:)
alt = (ph+phb)/9.81
do i = 0, nlay - 1
alt(i,:,:)=alt(i,:,:)-hgt
end do
mean_alt = new( nlay, "float", "No_FillValue" )
do i = 0, nlay - 1
mean_alt(i) = avg(alt(i,:,:)-alt(0,:,:))
if (mean_alt(i) .lt. emis_height_max) then
nlay_emis = min( (/ i + 1, nlay - 1 /) )
end if
end do
wrflevheights = mean_alt(0:(nlay_emis-1))
; point emission files - plume rise calculations
npointFiles = dimsizes(inlnFiles)
points = new( (/ ntime, nlay_emis, nrow_wrf, ncol_wrf, nvar /), "float", "No_FillValue" )
points!0 = "TSTEP"
points!1 = "LAY"
points!2 = "ROW"
points!3 = "COL"
points!4 = "VAR"
do i = 0, npointFiles - 1
pointData = addfile(inlnFiles(i), "r")
pointProp = addfile(stackFiles(i), "r")
pointDims = getfiledimsizes(pointProp)
npoints = pointDims(4) ; we "just know" that the fourth dimension of the point properties file is "ROW", which represents the number of point sources in file
pointDims = getfiledimsizes(pointData)
nvarfile = pointDims(3) ; we "just know" that the third dimension of the point data file is "VAR", which represents the number of species
print("Working on "+inlnFiles(i))
lons = pointProp->LONGITUDE(0,0,:,0)
lats = pointProp->LATITUDE(0,0,:,0)
idxFile = "tmp/stacklocations_inlnfile"+i+".nc"
; fires move - we need to recalculate the positions for each date individually
; whereas stacks can be cached...
if ((i .eq. fireidx) .or. (.not. isfilepresent(idxFile)))
print("no locations file for "+stackFiles(i)+" present - recreating")
; 0 is False, 1 is True
good = new(npoints, integer)
good(:) = 0
ii = new(npoints, integer)
ii(:) = -1
jj = new(npoints, integer)
jj(:) = -1
do j = 0, npoints - 1
; position in grid
if ( (lons(j) .gt. wrflonmin) .and. \
(lons(j) .lt. wrflonmax) .and. \
(lats(j) .gt. wrflatmin) .and. \
(lats(j) .lt. wrflatmax)) then
pos = getind_latlon2d(wrflat, wrflon, lats(j), lons(j))
ii(j) = pos(0,0)
jj(j) = pos(0,1)
; valid point?
if (ii(j) .ge. 0 .and. \
ii(j) .lt. nrow_wrf .and. \
jj(j) .ge. 0 .and. \
jj(j) .lt. ncol_wrf)
good(j) = 1
end if
end if
print(j+" of "+(npoints-1))
end do
if (i .ne. fireidx)
idxData = addfile(idxFile, "c")
idxData->good = good
idxData->ii = ii
idxData->jj = jj
end if
else
idxData = addfile(idxFile, "r")
good = idxData->good
ii = idxData->ii
jj = idxData->jj
end if
dia = pointProp->STKDM(0,0,:,0) ; stack diameter (m)
height = pointProp->STKHT(0,0,:,0) ; stack height (m)
temp = pointProp->STKTK(0,0,:,0) ; stack exit temperature (K)
vel = pointProp->STKVE(0,0,:,0) ; stack exit velocity (m/s)
flow = pointProp->STKFLW(0,0,:,0) ; stack exit flow rate (m3/s)
ngroup = pointProp->STKCNT(0,0,:,0) ; number of stacks in group
pres_e = 1013.25 ; hPa - fixed
u_e = 2.5 ; m/s - fixed
temp_e = 289.15 ; K - fixed
Rgas = 8.314 ; J K-1 mol-1 - univ. gas constant
mw = 0.02896 ; kg/mol - stack gas molweight (air?!)
cp = 29.19 ; J K-1 mol-1 - stack gas specific heat at const. pressure (air?!)
pi = acos(-1.0)
varData = new( (/ nvar, ntime, npoints /), "float", "No_FillValue" )
varData(:,:,:) = 0.0
do k = 0, nvar - 1
if (isfilevar(pointData, vars(k))) then
varData(k,:,:) = pointData->$vars(k)$(:,0,:,0)
end if
end do
print(npoints+" point sources to process")
do j = 0, npoints - 1
; good point? Calculated above...
if ( good(j) .gt. 0 ) then
; plume rise as in eqn. 17 of
; James E. Carson & Harry Moses (1969):
; The Validity of Several Plume Rise Formulas, Journal of the
; Air Pollution Control Association, 19:11, 862-866
; stack gas mass flow rate (kg/s)
m = (pi * flow(j) * pres_e * mw) / (4.0 * Rgas * temp(j))
; heat emission rate (kJ/s)
qh = m * cp * (temp(j) - temp_e)
; plume rise (m)
dz = -0.029 * (vel(j) * dia(j)) / u_e + 2.62 * ( sqrt( qh ) / u_e )
levdiff = abs( (height(j) + dz) - alt(:,ii(j),jj(j)) )
levidx = ind( levdiff .eq. min(levdiff) )
levidx = min( (/ levidx, (nlay_emis-1) /) )
do x = 0, nvar - 1
points(:,levidx,ii(j),jj(j),x) = points(:,levidx,ii(j),jj(j),x) + varData(x,:,j) * ngroup(j)
end do
if(i .eq. fireidx) then
points(:,levidx,ii(j),jj(j),ind(vars .eq. "CO_BB")) = points(:,levidx,ii(j),jj(j),ind(vars .eq. "CO_BB")) + varData(ind(vars .eq. "CO"),:,j) * ngroup(j)
else
points(:,levidx,ii(j),jj(j),ind(vars .eq. "CO_A")) = points(:,levidx,ii(j),jj(j),ind(vars .eq. "CO_A")) + varData(ind(vars .eq. "CO"),:,j) * ngroup(j)
end if
end if
end do
delete(varData)
delete(lons)
delete(lats)
delete(ii)
delete(jj)
delete(good)
delete(dia)
delete(height)
delete(temp)
delete(vel)
delete(flow)
delete(ngroup)
end do
; -----------------------------------------------------------------------------
; delete tmp input files
; -----------------------------------------------------------------------------
system("rm "+layer1File)
npointFiles = dimsizes(inlnFiles)
do i = 0, npointFiles - 1
system("rm "+inlnFiles(i))
system("rm "+stackFiles(i))
end do
; -----------------------------------------------------------------------------
; merge point and area sources
; -----------------------------------------------------------------------------
print("merging")
; "area" (ntime, nrow_wrf, ncol_wrf, nvar) contains area emissions
; "points" (ntime, nlay_emis, nrow_wrf, ncol_wrf, nvar) contains gridded point emissions
; simplest way to do: add area emissions to lowest level of point sources
do j = 0, nvar - 1
points(:,0,:,:,j) = points(:,0,:,:,j) + area(:,:,:,j)
end do
delete(area)
; -----------------------------------------------------------------------------
; conversion to MOZART-MOSAIC speciation
; -----------------------------------------------------------------------------
print("Conversion to MOZART-MOSAIC species")
nvarout = 46
nvargas = 29
outvars = new( nvarout, "string", "No_FillValue" )
out = new( (/ ntime, nlay_emis, nrow_wrf, ncol_wrf, nvarout /), "float", "No_FillValue" )
out!0 = "TSTEP"
out!1 = "LAY"
out!2 = "ROW"
out!3 = "COL"
out!4 = "VAR"
; keep ordering - gases first! (units in output are different)
out(:,:,:,:,:) = 0.0
i = 0
; JO: assume acrolein is methacrolein
; JO: don't have 1,3-butadiene, and nobody cares anyway about it
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "ACROLEIN")) + \
points(:,:,:,:,ind(vars .eq. "BUTADIENE13"))
outvars(i) = "E_MACR"
i = i + 1
; JO: we don't represent higher aldehydes, put ALDX into acetaldehyde
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "ALD2")) + \
points(:,:,:,:,ind(vars .eq. "ALDX"))
outvars(i) = "E_CH3CHO"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "BENZENE"))
outvars(i) = "E_BENZENE"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CH4"))
outvars(i) = "E_CH4"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CO"))
outvars(i) = "E_CO"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CO_A"))
outvars(i) = "E_CO_A"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CO_BB"))
outvars(i) = "E_CO_BB"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "ETHA"))
outvars(i) = "E_C2H6"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "ETH"))
outvars(i) = "E_C2H4"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "ETOH"))
outvars(i) = "E_C2H5OH"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "FORM"))
outvars(i) = "E_CH2O"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "ISOP"))
outvars(i) = "E_ISOP"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "MEOH"))
outvars(i) = "E_CH3OH"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "NH3"))
outvars(i) = "E_NH3"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "NO2"))
outvars(i) = "E_NO2"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "NO"))
outvars(i) = "E_NO"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "HONO"))
outvars(i) = "E_HONO"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "SO2"))
outvars(i) = "E_SO2"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "SULF"))
outvars(i) = "E_SULF"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "TERP"))
outvars(i) = "E_C10H16"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "TOL"))
outvars(i) = "E_TOLUENE"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "XYL"))
outvars(i) = "E_XYLENE"
; now the tricky ones
; CBO5 is a carbon-bond mechanism (lumped structure),
; while MOZART uses surrogate species (lumped molecules)
; need to translate that (i.e. make an educated guess):
;
; IOLE (internal olefine carbon bond)
; OLE (terminal olefine carbon bond)
; PAR (paraffins)
;
; OLE + PAR = C3H6
; OLE + 2 PAR = BIGENE
; IOLE + 2 PAR = BIGENE
; 5 PAR = BIGALK
;
; hence (mole-wise)
;
; C3H6 + BIGENE + BIGALK = a(OLE + PAR)/2 + b(OLE + 2 PAR)/3 + c(IOLE + 2 PAR)/3 + d(5 PAR)/5
; from previously speciated data (CalNex, based on scaling to CO by Borbon et al. 2013):
a = 0.1433
bc = 0.0196 ; that's b + c
d = 0.8371
;
i = i + 1
out(:,:,:,:,i) = 0.0
outvars(i) = "E_C3H6"
c3h6_idx = i
i = i + 1
out(:,:,:,:,i) = 0.0
outvars(i) = "E_BIGENE"
bigene_idx = i
i = i + 1
out(:,:,:,:,i) = 0.0
outvars(i) = "E_BIGALK"
bigalk_idx = i
do x = 0, ncol_wrf - 1
do y = 0, nrow_wrf - 1
do z = 0, nlay_emis - 1
do t = 0, ntime - 1
par = points(t,z,y,x,ind(vars .eq. "PAR"))
ole = points(t,z,y,x,ind(vars .eq. "OLE"))
iole = points(t,z,y,x,ind(vars .eq. "IOLE"))
c3h6 = 0.0
bigene = 0.0
bigalk = 0.0
; first to create: C3H6 (~14%)
moles = min( (/ a * ole, par /) )
c3h6 = moles / 2.0
par = par - moles
ole = ole - moles
; then BIGENE (~2%)
moles = min( (/ ole, par / 2.0 /) )
bigene = moles / 3.0
ole = ole - moles
par = par - moles
moles = min( (/ iole, par / 2.0 /) )
bigene = bigene + moles / 3.0
iole = iole - moles
par = par - moles
; the rest is BIGALK (~83%)
bigalk = par / 5.0
out(t,z,y,x,c3h6_idx) = c3h6
out(t,z,y,x,bigene_idx) = bigene
out(t,z,y,x,bigalk_idx) = bigalk
end do
end do
end do
end do
; some VOC emissions are unknown - we scale them based on Borbon et al., 2012
; to CO
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CO_A")) * 1.12E-02
outvars(i) = "E_C3H8"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CO_A")) * 1.18E-02
outvars(i) = "E_CH3COCH3"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CO_A")) * 2.40E-04
outvars(i) = "E_MVK"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "CO_A")) * 5.87E-03
outvars(i) = "E_C2H2"
; --- aerosols ---
; split between Aitken and accumulation mode emitted mass:
aitkenFrac = 0.1
i = i + 1
out(:,:,:,:,i) = aitkenFrac * \
(points(:,:,:,:,ind(vars .eq. "PAL")) + \
points(:,:,:,:,ind(vars .eq. "PCA")) + \
points(:,:,:,:,ind(vars .eq. "PFE")) + \
points(:,:,:,:,ind(vars .eq. "PK")) + \
points(:,:,:,:,ind(vars .eq. "PMFINE")) + \
points(:,:,:,:,ind(vars .eq. "PMG")) + \
points(:,:,:,:,ind(vars .eq. "PMN")) + \
points(:,:,:,:,ind(vars .eq. "PMOTHR")) + \
points(:,:,:,:,ind(vars .eq. "PSI")) + \
points(:,:,:,:,ind(vars .eq. "PTI")))
outvars(i) = "E_PM25I"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * \
(points(:,:,:,:,ind(vars .eq. "PAL")) + \
points(:,:,:,:,ind(vars .eq. "PCA")) + \
points(:,:,:,:,ind(vars .eq. "PFE")) + \
points(:,:,:,:,ind(vars .eq. "PK")) + \
points(:,:,:,:,ind(vars .eq. "PMFINE")) + \
points(:,:,:,:,ind(vars .eq. "PMG")) + \
points(:,:,:,:,ind(vars .eq. "PMN")) + \
points(:,:,:,:,ind(vars .eq. "PMOTHR")) + \
points(:,:,:,:,ind(vars .eq. "PSI")) + \
points(:,:,:,:,ind(vars .eq. "PTI")))
outvars(i) = "E_PM25J"
i = i + 1
out(:,:,:,:,i) = aitkenFrac * points(:,:,:,:,ind(vars .eq. "PNA"))
outvars(i) = "E_NAI"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * points(:,:,:,:,ind(vars .eq. "PNA"))
outvars(i) = "E_NAJ"
i = i + 1
out(:,:,:,:,i) = aitkenFrac * points(:,:,:,:,ind(vars .eq. "PCL"))
outvars(i) = "E_CLI"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * points(:,:,:,:,ind(vars .eq. "PCL"))
outvars(i) = "E_CLJ"
i = i + 1
out(:,:,:,:,i) = aitkenFrac * points(:,:,:,:,ind(vars .eq. "PEC"))
outvars(i) = "E_ECI"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * points(:,:,:,:,ind(vars .eq. "PEC"))
outvars(i) = "E_ECJ"
i = i + 1
out(:,:,:,:,i) = aitkenFrac * \
(points(:,:,:,:,ind(vars .eq. "PNCOM")) + \
points(:,:,:,:,ind(vars .eq. "POC")))
outvars(i) = "E_ORGI"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * \
(points(:,:,:,:,ind(vars .eq. "PNCOM")) + \
points(:,:,:,:,ind(vars .eq. "POC")))
outvars(i) = "E_ORGJ"
i = i + 1
out(:,:,:,:,i) = aitkenFrac * points(:,:,:,:,ind(vars .eq. "PSO4"))
outvars(i) = "E_SO4I"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * points(:,:,:,:,ind(vars .eq. "PSO4"))
outvars(i) = "E_SO4J"
i = i + 1
out(:,:,:,:,i) = aitkenFrac * points(:,:,:,:,ind(vars .eq. "PNO3"))
outvars(i) = "E_NO3I"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * points(:,:,:,:,ind(vars .eq. "PNO3"))
outvars(i) = "E_NO3J"
i = i + 1
out(:,:,:,:,i) = aitkenFrac * points(:,:,:,:,ind(vars .eq. "PNH4"))
outvars(i) = "E_NH4I"
i = i + 1
out(:,:,:,:,i) = (1.0 - aitkenFrac) * points(:,:,:,:,ind(vars .eq. "PNH4"))
outvars(i) = "E_NH4J"
i = i + 1
out(:,:,:,:,i) = points(:,:,:,:,ind(vars .eq. "PMC"))
outvars(i) = "E_PM_10"
; -----------------------------------------------------------------------------
; create wrfchemi files
; -----------------------------------------------------------------------------
print("output")
; dump raw data (in CB05 speciation as given by EPA)
dmp = addfile("output/NEI_raw_"+date+".nc", "c")
do j = 0, nvar - 1
dmp->$vars(j)$ = points(:,:,:,:,j)
end do
delete(dmp)
delete(points)
; copy global attributes from WRF domain file
attsToCopy = (/ "SIMULATION_START_DATE", "WEST-EAST_GRID_DIMENSION", \
"SOUTH-NORTH_GRID_DIMENSION", "BOTTOM-TOP_GRID_DIMENSION", \
"DX", "DY", "GRIDTYPE", "DIFF_OPT", "KM_OPT", "DAMP_OPT", \
"DAMPCOEF", "KHDIF", "KVDIF", "MP_PHYSICS", "RA_LW_PHYSICS", \
"RA_SW_PHYSICS", "SF_SFCLAY_PHYSICS", "SF_SURFACE_PHYSICS", \
"BL_PBL_PHYSICS", "CU_PHYSICS", "SURFACE_INPUT_SOURCE", \
"SST_UPDATE", "GRID_FDDA", "GFDDA_INTERVAL_M", "GFDDA_END_H", \
"GRID_SFDDA", "SGFDDA_INTERVAL_M", "SGFDDA_END_H", \
"HYPSOMETRIC_OPT", "WEST-EAST_PATCH_START_UNSTAG", \
"WEST-EAST_PATCH_END_UNSTAG", "WEST-EAST_PATCH_START_STAG", \
"WEST-EAST_PATCH_END_STAG", "SOUTH-NORTH_PATCH_START_UNSTAG", \
"SOUTH-NORTH_PATCH_END_UNSTAG", \
"SOUTH-NORTH_PATCH_START_STAG", \
"SOUTH-NORTH_PATCH_END_STAG", \
"BOTTOM-TOP_PATCH_START_UNSTAG", \
"BOTTOM-TOP_PATCH_END_UNSTAG", \
"BOTTOM-TOP_PATCH_START_STAG", \
"BOTTOM-TOP_PATCH_END_STAG", "GRID_ID", "PARENT_ID", \
"I_PARENT_START", "J_PARENT_START", "PARENT_GRID_RATIO", \
"DT", "CEN_LAT", "CEN_LON", "TRUELAT1", "TRUELAT2", \
"MOAD_CEN_LAT", "STAND_LON", "POLE_LAT", "POLE_LON", "GMT", \
"JULYR", "JULDAY", "MAP_PROJ", "MMINLU", "NUM_LAND_CAT", \
"ISWATER", "ISLAKE", "ISICE", "ISURBAN", "ISOILWATER" /)
natts = dimsizes(attsToCopy)
att = True
do i = 0, natts - 1
att@$attsToCopy(i)$ = wrfgrid@$attsToCopy(i)$
end do
; dimensions
dimNames = (/ "Time", "emissions_zdim_stag", "south_north", "west_east", "DateStrLen" /)
dimSizes = (/ -1 , nlay_emis , nrow_wrf , ncol_wrf , 19 /)
dimUnlim = (/ True , False , False , False , False /)
; unit conversion
; gases: moles s-1 -> moles km-2 hr-1
conv_gas = 3600.0 / WRF_grid_area
; area interpolation (better: aggregation) should have conserved mass,
; point sources were just added to each cell...
; aerosols: g s-1 -> ug m-2 s-1
conv_aer = 1.0 / WRF_grid_area
do t = 0, ntime - 1
year = str_get_cols(date, 0, 3)
month = str_get_cols(date, 4, 5)
day = str_get_cols(date, 6, 7)
hour = "0"+t
if (t .gt. 9)
hour = tostring(t)
end if
niceDate = year + "-" + month + "-" + day + "_" + hour + ":00:00"
outfileName = "output/wrfchemi_d01_"+niceDate+".nc"
if (isfilepresent(outfileName)) then
system("rm "+outfileName)
end if
outfile = addfile(outfileName, "c")
setfileoption(outfile, "DefineMode", True)
; global attributes
fileattdef( outfile, att )
; dimensions
filedimdef(outfile, dimNames, dimSizes, dimUnlim)
setfileoption(outfile, "DefineMode", False)
do j = 0, nvarout - 1
var = new( (/ 1, nlay_emis, nrow_wrf, ncol_wrf /), "float", "No_FillValue" )
var!0 = "Time"
var!1 = "emissions_zdim_stag"
var!2 = "south_north"
var!3 = "west_east"
var@FieldType = 104
var@MemoryOrder = "XYZ"
var@description = "EMISSIONS"
conv_var = 1.0
if (j .lt. nvargas)
; gases
var@units = "mole km-2 hr-1"
conv_var = conv_gas
else
; aerosols
var@units = "ug m-2 s-1"
conv_var = conv_aer
end if
var = (/ out(t,:,:,:,j) * conv_var /)
outfile->$outvars(j)$ = var
delete(var)
end do
times = new((/1,19/), character)
times!0 = "Time"
times!1 = "DateStrLen"
outfile->Times = times
niceDateC = stringtochar(niceDate)
outfile->Times(0,:) = (/ niceDateC(:18) /)
delete(times)
delete(outfile)
print("Written wrfchemi file for "+niceDate)
end do
You will need this small IDL script to create the NEI grid file:
; purpose:
; create a NetCDF file with the lat-lon locations of grid cells
; in a given Lambert Conformal Conical grid (NEI emissions).
; writes corner and center positions as 2D arrays to "NEI_grid.nc".
NCOLS = 459L
NROWS = 299L
XORIG = -2556000.
YORIG = -1728000.
XCELL = 12000.
YCELL = 12000.
xidx = fltarr(NCOLS, NROWS)
yidx = fltarr(NCOLS, NROWS)
; add & $ for multiline execution
;FOR i = 0L, ncols - 1 DO BEGIN & $
; xidx[i,*] = i & $
; FOR j = 0L, nrows - 1 DO BEGIN & $
; yidx[i,j] = j & $
; ENDFOR & $
;ENDFOR
FOR i = 0L, ncols - 1 DO BEGIN
xidx[i,*] = i
FOR j = 0L, nrows - 1 DO BEGIN
yidx[i,j] = j
ENDFOR
ENDFOR
llx = xorig + xidx * xcell
lly = yorig + yidx * ycell
lrx = xorig + (xidx + 1.) * xcell
lry = yorig + yidx * ycell
ulx = xorig + xidx * xcell
uly = yorig + (yidx + 1.) * ycell
urx = xorig + (xidx + 1.) * xcell
ury = yorig + (yidx + 1.) * ycell
cex = xorig + (xidx + 0.5) * xcell
cey = yorig + (yidx + 0.5) * ycell
mapStruct = MAP_PROJ_INIT(104, /GCTP, CENTER_LONGITUDE=-97., CENTER_LATITUDE=40., $
STANDARD_PAR1=33., STANDARD_PAR2=45., DATUM=19)
ijnei = MAP_PROJ_INVERSE(llx, lly, MAP_STRUCTURE=mapStruct)
lllon = reform(ijnei[0,*], ncols, nrows)
lllat = reform(ijnei[1,*], ncols, nrows)
ijnei = MAP_PROJ_INVERSE(lrx, lry, MAP_STRUCTURE=mapStruct)
lrlon = reform(ijnei[0,*], ncols, nrows)
lrlat = reform(ijnei[1,*], ncols, nrows)
ijnei = MAP_PROJ_INVERSE(ulx, uly, MAP_STRUCTURE=mapStruct)
ullon = reform(ijnei[0,*], ncols, nrows)
ullat = reform(ijnei[1,*], ncols, nrows)
ijnei = MAP_PROJ_INVERSE(urx, ury, MAP_STRUCTURE=mapStruct)
urlon = reform(ijnei[0,*], ncols, nrows)
urlat = reform(ijnei[1,*], ncols, nrows)
ijnei = MAP_PROJ_INVERSE(cex, cey, MAP_STRUCTURE=mapStruct)
celon = reform(ijnei[0,*], ncols, nrows)
celat = reform(ijnei[1,*], ncols, nrows)
ncid = ncdf_create("NEI_grid.nc", /clobber)
; create dimensions
xid = ncdf_dimdef(ncid, 'X',NCOLS)
yid = ncdf_dimdef(ncid, 'Y',NROWS)
; define variables
llxid = ncdf_vardef(ncid,'LLx',[xid, yid],/float)
llyid = ncdf_vardef(ncid,'LLy',[xid, yid],/float)
lllonid = ncdf_vardef(ncid,'LLlon',[xid, yid],/float)
lllatid = ncdf_vardef(ncid,'LLlat',[xid, yid],/float)
lrlonid = ncdf_vardef(ncid,'LRlon',[xid, yid],/float)
lrlatid = ncdf_vardef(ncid,'LRlat',[xid, yid],/float)
ullonid = ncdf_vardef(ncid,'ULlon',[xid, yid],/float)
ullatid = ncdf_vardef(ncid,'ULlat',[xid, yid],/float)
urlonid = ncdf_vardef(ncid,'URlon',[xid, yid],/float)
urlatid = ncdf_vardef(ncid,'URlat',[xid, yid],/float)
celonid = ncdf_vardef(ncid,'CElon',[xid, yid],/float)
celatid = ncdf_vardef(ncid,'CElat',[xid, yid],/float)
ncdf_control,ncid, /endef
ncdf_varput,ncid,'LLx',llx
ncdf_varput,ncid,'LLy',lly
ncdf_varput,ncid,'LLlon',lllon
ncdf_varput,ncid,'LLlat',lllat
ncdf_varput,ncid,'LRlon',lrlon
ncdf_varput,ncid,'LRlat',lrlat
ncdf_varput,ncid,'ULlon',ullon
ncdf_varput,ncid,'ULlat',ullat
ncdf_varput,ncid,'URlon',urlon
ncdf_varput,ncid,'URlat',urlat
ncdf_varput,ncid,'CElon',celon
ncdf_varput,ncid,'CElat',celat
ncdf_close,ncid
Operational setup
Runscripts for chained runs with meteo spinup will be linked here soon…
Major modifications / extensions:
Accumulated dry deposition amounts as output fields
dry_dep_driver.F, dry_dep_driver:
END SELECT mix_select
! CK 20130116 write out dry deposition amounts >>
! chem is in ppmv
! dry deposition is mangled with vertical mixing, but column independent.
! Hence, all molecules lost per column must be dry deposited.
! old and new column totals (mol/m2 or ug/m2)
old = 0.0
new = 0.0
do k=kts,kte-1
fac = 1.0
if (nv <= numgas) then
! from ppmv to mol/m2
! fac = 1e-6 * rho * 1/mw_air * dz
! kg/m3 mol/kg m
! CK 20131028 bugfix - mwdry is in g/mol, to get kg/mol we need to use 1e-3, not 1e3
! fac = 1e-6 * dryrho_1d(k) * 1./(mwdry*1.e3) * dz8w(i,k,j)
fac = 1e-6 * dryrho_1d(k) * 1./(mwdry*1.e-3) * dz8w(i,k,j)
else
! from ug/kg to ug/m2
! fac = rho * dz
! kg/m3 m
fac = dryrho_1d(k) * dz8w(i,k,j)
endif
old = old + max(epsilc,chem(i,k,j,nv)) * fac
new = new + max(epsilc,pblst(k)) * fac
enddo
! we add new dry deposition to existing field (accumulated deposition!)
ddmassn(nv) = max( 0.0, (old - new) )
! CK <<
do k=kts,kte-1
chem(i,k,j,nv)=max(epsilc,pblst(k))
enddo
enddo
! CK 20130116 write out dry deposition amounts >>
if( config_flags%diagnostic_chem == DEPVEL1 .and. &
(config_flags%chem_opt == MOZCART_KPP .or. &
config_flags%chem_opt == MOZART_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_AQ_KPP) ) then
dvel(i,1,j,p_ddmass_o3) = dvel(i,1,j,p_ddmass_o3) + ddmassn(p_o3)
dvel(i,1,j,p_ddmass_no) = dvel(i,1,j,p_ddmass_no) + ddmassn(p_no)
dvel(i,1,j,p_ddmass_no2) = dvel(i,1,j,p_ddmass_no2) + ddmassn(p_no2)
dvel(i,1,j,p_ddmass_nh3) = dvel(i,1,j,p_ddmass_nh3) + ddmassn(p_nh3)
dvel(i,1,j,p_ddmass_hno3) = dvel(i,1,j,p_ddmass_hno3) + ddmassn(p_hno3)
dvel(i,1,j,p_ddmass_hno4) = dvel(i,1,j,p_ddmass_hno4) + ddmassn(p_hno4)
dvel(i,1,j,p_ddmass_h2o2) = dvel(i,1,j,p_ddmass_h2o2) + ddmassn(p_h2o2)
dvel(i,1,j,p_ddmass_co) = dvel(i,1,j,p_ddmass_co) + ddmassn(p_co)
dvel(i,1,j,p_ddmass_ch3ooh) = dvel(i,1,j,p_ddmass_ch3ooh) + ddmassn(p_ch3ooh)
dvel(i,1,j,p_ddmass_hcho) = dvel(i,1,j,p_ddmass_hcho) + ddmassn(p_hcho)
dvel(i,1,j,p_ddmass_ch3oh) = dvel(i,1,j,p_ddmass_ch3oh) + ddmassn(p_ch3oh)
dvel(i,1,j,p_ddmass_eo2) = dvel(i,1,j,p_ddmass_eo2) + ddmassn(p_eo2)
dvel(i,1,j,p_ddmass_ald) = dvel(i,1,j,p_ddmass_ald) + ddmassn(p_ald)
dvel(i,1,j,p_ddmass_ch3cooh) = dvel(i,1,j,p_ddmass_ch3cooh) + ddmassn(p_ch3cooh)
dvel(i,1,j,p_ddmass_acet) = dvel(i,1,j,p_ddmass_acet) + ddmassn(p_acet)
dvel(i,1,j,p_ddmass_mgly) = dvel(i,1,j,p_ddmass_mgly) + ddmassn(p_mgly)
dvel(i,1,j,p_ddmass_gly) = dvel(i,1,j,p_ddmass_gly) + ddmassn(p_gly)
dvel(i,1,j,p_ddmass_paa) = dvel(i,1,j,p_ddmass_paa) + ddmassn(p_paa)
dvel(i,1,j,p_ddmass_pooh) = dvel(i,1,j,p_ddmass_pooh) + ddmassn(p_c3h6ooh)
dvel(i,1,j,p_ddmass_mpan) = dvel(i,1,j,p_ddmass_mpan) + ddmassn(p_mpan)
dvel(i,1,j,p_ddmass_mco3) = dvel(i,1,j,p_ddmass_mco3) + ddmassn(p_mco3)
dvel(i,1,j,p_ddmass_mvkooh) = dvel(i,1,j,p_ddmass_mvkooh) + ddmassn(p_mvkooh)
dvel(i,1,j,p_ddmass_c2h5oh) = dvel(i,1,j,p_ddmass_c2h5oh) + ddmassn(p_c2h5oh)
dvel(i,1,j,p_ddmass_etooh) = dvel(i,1,j,p_ddmass_etooh) + ddmassn(p_etooh)
dvel(i,1,j,p_ddmass_prooh) = dvel(i,1,j,p_ddmass_prooh) + ddmassn(p_prooh)
dvel(i,1,j,p_ddmass_acetp) = dvel(i,1,j,p_ddmass_acetp) + ddmassn(p_acetp)
dvel(i,1,j,p_ddmass_onit) = dvel(i,1,j,p_ddmass_onit) + ddmassn(p_onit)
dvel(i,1,j,p_ddmass_onitr) = dvel(i,1,j,p_ddmass_onitr) + ddmassn(p_onitr)
dvel(i,1,j,p_ddmass_isooh) = dvel(i,1,j,p_ddmass_isooh) + ddmassn(p_isooh)
dvel(i,1,j,p_ddmass_acetol) = dvel(i,1,j,p_ddmass_acetol) + ddmassn(p_acetol)
dvel(i,1,j,p_ddmass_glyald) = dvel(i,1,j,p_ddmass_glyald) + ddmassn(p_glyald)
dvel(i,1,j,p_ddmass_hydrald) = dvel(i,1,j,p_ddmass_hydrald) + ddmassn(p_hydrald)
dvel(i,1,j,p_ddmass_alkooh) = dvel(i,1,j,p_ddmass_alkooh) + ddmassn(p_alkooh)
dvel(i,1,j,p_ddmass_mekooh) = dvel(i,1,j,p_ddmass_mekooh) + ddmassn(p_mekooh)
dvel(i,1,j,p_ddmass_tolooh) = dvel(i,1,j,p_ddmass_tolooh) + ddmassn(p_tolooh)
dvel(i,1,j,p_ddmass_xooh) = dvel(i,1,j,p_ddmass_xooh) + ddmassn(p_xooh)
dvel(i,1,j,p_ddmass_so2) = dvel(i,1,j,p_ddmass_so2) + ddmassn(p_so2)
dvel(i,1,j,p_ddmass_so4) = dvel(i,1,j,p_ddmass_so4) + ddmassn(p_sulf)
dvel(i,1,j,p_ddmass_pan) = dvel(i,1,j,p_ddmass_pan) + ddmassn(p_pan)
dvel(i,1,j,p_ddmass_terpooh) = dvel(i,1,j,p_ddmass_terpooh) + ddmassn(p_terpooh)
if (config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_AQ_KPP) then
dvel(i,1,j,p_ddmass_so4_a01) = dvel(i,1,j,p_ddmass_so4_a01) + ddmassn(p_so4_a01)
dvel(i,1,j,p_ddmass_no3_a01) = dvel(i,1,j,p_ddmass_no3_a01) + ddmassn(p_no3_a01)
dvel(i,1,j,p_ddmass_cl_a01) = dvel(i,1,j,p_ddmass_cl_a01) + ddmassn(p_cl_a01)
dvel(i,1,j,p_ddmass_nh4_a01) = dvel(i,1,j,p_ddmass_nh4_a01) + ddmassn(p_nh4_a01)
dvel(i,1,j,p_ddmass_na_a01) = dvel(i,1,j,p_ddmass_na_a01) + ddmassn(p_na_a01)
dvel(i,1,j,p_ddmass_oin_a01) = dvel(i,1,j,p_ddmass_oin_a01) + ddmassn(p_oin_a01)
dvel(i,1,j,p_ddmass_oc_a01) = dvel(i,1,j,p_ddmass_oc_a01) + ddmassn(p_oc_a01)
dvel(i,1,j,p_ddmass_bc_a01) = dvel(i,1,j,p_ddmass_bc_a01) + ddmassn(p_bc_a01)
dvel(i,1,j,p_ddmass_smpa_a01) = dvel(i,1,j,p_ddmass_smpa_a01) + ddmassn(p_smpa_a01)
dvel(i,1,j,p_ddmass_smpbb_a01) = dvel(i,1,j,p_ddmass_smpbb_a01) + ddmassn(p_smpbb_a01)
dvel(i,1,j,p_ddmass_glysoa_a01) = dvel(i,1,j,p_ddmass_glysoa_a01) + &
ddmassn(p_glysoa_r1_a01) + &
ddmassn(p_glysoa_r2_a01) + &
ddmassn(p_glysoa_pht_a01) + &
ddmassn(p_glysoa_drk_a01)
dvel(i,1,j,p_ddmass_biog1_c_a01) = dvel(i,1,j,p_ddmass_biog1_c_a01) + ddmassn(p_biog1_c_a01)
dvel(i,1,j,p_ddmass_biog1_o_a01) = dvel(i,1,j,p_ddmass_biog1_o_a01) + ddmassn(p_biog1_o_a01)
dvel(i,1,j,p_ddmass_so4_a02) = dvel(i,1,j,p_ddmass_so4_a02) + ddmassn(p_so4_a02)
dvel(i,1,j,p_ddmass_no3_a02) = dvel(i,1,j,p_ddmass_no3_a02) + ddmassn(p_no3_a02)
dvel(i,1,j,p_ddmass_cl_a02) = dvel(i,1,j,p_ddmass_cl_a02) + ddmassn(p_cl_a02)
dvel(i,1,j,p_ddmass_nh4_a02) = dvel(i,1,j,p_ddmass_nh4_a02) + ddmassn(p_nh4_a02)
dvel(i,1,j,p_ddmass_na_a02) = dvel(i,1,j,p_ddmass_na_a02) + ddmassn(p_na_a02)
dvel(i,1,j,p_ddmass_oin_a02) = dvel(i,1,j,p_ddmass_oin_a02) + ddmassn(p_oin_a02)
dvel(i,1,j,p_ddmass_oc_a02) = dvel(i,1,j,p_ddmass_oc_a02) + ddmassn(p_oc_a02)
dvel(i,1,j,p_ddmass_bc_a02) = dvel(i,1,j,p_ddmass_bc_a02) + ddmassn(p_bc_a02)
dvel(i,1,j,p_ddmass_smpa_a02) = dvel(i,1,j,p_ddmass_smpa_a02) + ddmassn(p_smpa_a02)
dvel(i,1,j,p_ddmass_smpbb_a02) = dvel(i,1,j,p_ddmass_smpbb_a02) + ddmassn(p_smpbb_a02)
dvel(i,1,j,p_ddmass_glysoa_a02) = dvel(i,1,j,p_ddmass_glysoa_a02) + &
ddmassn(p_glysoa_r1_a02) + &
ddmassn(p_glysoa_r2_a02) + &
ddmassn(p_glysoa_pht_a02) + &
ddmassn(p_glysoa_drk_a02)
dvel(i,1,j,p_ddmass_biog1_c_a02) = dvel(i,1,j,p_ddmass_biog1_c_a02) + ddmassn(p_biog1_c_a02)
dvel(i,1,j,p_ddmass_biog1_o_a02) = dvel(i,1,j,p_ddmass_biog1_o_a02) + ddmassn(p_biog1_o_a02)
dvel(i,1,j,p_ddmass_so4_a03) = dvel(i,1,j,p_ddmass_so4_a03) + ddmassn(p_so4_a03)
dvel(i,1,j,p_ddmass_no3_a03) = dvel(i,1,j,p_ddmass_no3_a03) + ddmassn(p_no3_a03)
dvel(i,1,j,p_ddmass_cl_a03) = dvel(i,1,j,p_ddmass_cl_a03) + ddmassn(p_cl_a03)
dvel(i,1,j,p_ddmass_nh4_a03) = dvel(i,1,j,p_ddmass_nh4_a03) + ddmassn(p_nh4_a03)
dvel(i,1,j,p_ddmass_na_a03) = dvel(i,1,j,p_ddmass_na_a03) + ddmassn(p_na_a03)
dvel(i,1,j,p_ddmass_oin_a03) = dvel(i,1,j,p_ddmass_oin_a03) + ddmassn(p_oin_a03)
dvel(i,1,j,p_ddmass_oc_a03) = dvel(i,1,j,p_ddmass_oc_a03) + ddmassn(p_oc_a03)
dvel(i,1,j,p_ddmass_bc_a03) = dvel(i,1,j,p_ddmass_bc_a03) + ddmassn(p_bc_a03)
dvel(i,1,j,p_ddmass_smpa_a03) = dvel(i,1,j,p_ddmass_smpa_a03) + ddmassn(p_smpa_a03)
dvel(i,1,j,p_ddmass_smpbb_a03) = dvel(i,1,j,p_ddmass_smpbb_a03) + ddmassn(p_smpbb_a03)
dvel(i,1,j,p_ddmass_glysoa_a03) = dvel(i,1,j,p_ddmass_glysoa_a03) + &
ddmassn(p_glysoa_r1_a03) + &
ddmassn(p_glysoa_r2_a03) + &
ddmassn(p_glysoa_pht_a03) + &
ddmassn(p_glysoa_drk_a03)
dvel(i,1,j,p_ddmass_biog1_c_a03) = dvel(i,1,j,p_ddmass_biog1_c_a03) + ddmassn(p_biog1_c_a03)
dvel(i,1,j,p_ddmass_biog1_o_a03) = dvel(i,1,j,p_ddmass_biog1_o_a03) + ddmassn(p_biog1_o_a03)
dvel(i,1,j,p_ddmass_so4_a04) = dvel(i,1,j,p_ddmass_so4_a04) + ddmassn(p_so4_a04)
dvel(i,1,j,p_ddmass_no3_a04) = dvel(i,1,j,p_ddmass_no3_a04) + ddmassn(p_no3_a04)
dvel(i,1,j,p_ddmass_cl_a04) = dvel(i,1,j,p_ddmass_cl_a04) + ddmassn(p_cl_a04)
dvel(i,1,j,p_ddmass_nh4_a04) = dvel(i,1,j,p_ddmass_nh4_a04) + ddmassn(p_nh4_a04)
dvel(i,1,j,p_ddmass_na_a04) = dvel(i,1,j,p_ddmass_na_a04) + ddmassn(p_na_a04)
dvel(i,1,j,p_ddmass_oin_a04) = dvel(i,1,j,p_ddmass_oin_a04) + ddmassn(p_oin_a04)
dvel(i,1,j,p_ddmass_oc_a04) = dvel(i,1,j,p_ddmass_oc_a04) + ddmassn(p_oc_a04)
dvel(i,1,j,p_ddmass_bc_a04) = dvel(i,1,j,p_ddmass_bc_a04) + ddmassn(p_bc_a04)
dvel(i,1,j,p_ddmass_smpa_a04) = dvel(i,1,j,p_ddmass_smpa_a04) + ddmassn(p_smpa_a04)
dvel(i,1,j,p_ddmass_smpbb_a04) = dvel(i,1,j,p_ddmass_smpbb_a04) + ddmassn(p_smpbb_a04)
dvel(i,1,j,p_ddmass_glysoa_a04) = dvel(i,1,j,p_ddmass_glysoa_a04) + &
ddmassn(p_glysoa_r1_a04) + &
ddmassn(p_glysoa_r2_a04) + &
ddmassn(p_glysoa_pht_a04) + &
ddmassn(p_glysoa_drk_a04)
dvel(i,1,j,p_ddmass_biog1_c_a04) = dvel(i,1,j,p_ddmass_biog1_c_a04) + ddmassn(p_biog1_c_a04)
dvel(i,1,j,p_ddmass_biog1_o_a04) = dvel(i,1,j,p_ddmass_biog1_o_a04) + ddmassn(p_biog1_o_a04)
dvel(i,1,j,p_ddmass_ca_a01) = dvel(i,1,j,p_ddmass_ca_a01) + ddmassn(p_ca_a01)
dvel(i,1,j,p_ddmass_ca_a02) = dvel(i,1,j,p_ddmass_ca_a02) + ddmassn(p_ca_a02)
dvel(i,1,j,p_ddmass_ca_a03) = dvel(i,1,j,p_ddmass_ca_a03) + ddmassn(p_ca_a03)
dvel(i,1,j,p_ddmass_ca_a04) = dvel(i,1,j,p_ddmass_ca_a04) + ddmassn(p_ca_a04)
dvel(i,1,j,p_ddmass_co3_a01) = dvel(i,1,j,p_ddmass_co3_a01) + ddmassn(p_co3_a01)
dvel(i,1,j,p_ddmass_co3_a02) = dvel(i,1,j,p_ddmass_co3_a02) + ddmassn(p_co3_a02)
dvel(i,1,j,p_ddmass_co3_a03) = dvel(i,1,j,p_ddmass_co3_a03) + ddmassn(p_co3_a03)
dvel(i,1,j,p_ddmass_co3_a04) = dvel(i,1,j,p_ddmass_co3_a04) + ddmassn(p_co3_a04)
endif
endif
! CK 20130116 write out dry deposition amounts <<
Aerosols feed back on photolysis rates:
module_ftuv_driver.F, ftuv_driver:
elseif ( config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_AQ_KPP ) then
atm_mass_den = rho_phy(i,k,j)
! smaller bins are probably not mass relevant, added anyway
acb1(kp1) = acb1(kp1) + &
(chem(i,k,j,p_bc_a01) + &
chem(i,k,j,p_bc_a02) + &
chem(i,k,j,p_bc_a03) + &
chem(i,k,j,p_bc_a04)) * atm_mass_den
acb2(kp1) = 0._dp
aoc1(kp1) = aoc1(kp1) + &
(chem(i,k,j,p_oc_a01) + &
chem(i,k,j,p_oc_a02) + &
chem(i,k,j,p_oc_a03) + &
chem(i,k,j,p_oc_a04)) * atm_mass_den
aoc2(kp1) = aoc2(kp1) + &
(chem(i,k,j,p_smpa_a01) + &
chem(i,k,j,p_smpa_a02) + &
chem(i,k,j,p_smpa_a03) + &
chem(i,k,j,p_smpa_a04)) * atm_mass_den + &
(chem(i,k,j,p_smpbb_a01) + &
chem(i,k,j,p_smpbb_a02) + &
chem(i,k,j,p_smpbb_a03) + &
chem(i,k,j,p_smpbb_a04)) * atm_mass_den + &
(chem(i,k,j,p_biog1_c_a01) + &
chem(i,k,j,p_biog1_c_a02) + &
chem(i,k,j,p_biog1_c_a03) + &
chem(i,k,j,p_biog1_c_a04)) * atm_mass_den + &
(chem(i,k,j,p_biog1_o_a01) + &
chem(i,k,j,p_biog1_o_a02) + &
chem(i,k,j,p_biog1_o_a03) + &
chem(i,k,j,p_biog1_o_a04)) * atm_mass_den + &
(chem(i,k,j,p_glysoa_r1_a01) + &
chem(i,k,j,p_glysoa_r1_a02) + &
chem(i,k,j,p_glysoa_r1_a03) + &
chem(i,k,j,p_glysoa_r1_a04)) * atm_mass_den + &
(chem(i,k,j,p_glysoa_r2_a01) + &
chem(i,k,j,p_glysoa_r2_a02) + &
chem(i,k,j,p_glysoa_r2_a03) + &
chem(i,k,j,p_glysoa_r2_a04)) * atm_mass_den + &
(chem(i,k,j,p_glysoa_pht_a01) + &
chem(i,k,j,p_glysoa_pht_a02) + &
chem(i,k,j,p_glysoa_pht_a03) + &
chem(i,k,j,p_glysoa_pht_a04)) * atm_mass_den + &
(chem(i,k,j,p_glysoa_drk_a01) + &
chem(i,k,j,p_glysoa_drk_a02) + &
chem(i,k,j,p_glysoa_drk_a03) + &
chem(i,k,j,p_glysoa_drk_a04)) * atm_mass_den
aant(kp1) = aant(kp1) + &
(chem(i,k,j,p_nh4_a01) + &
chem(i,k,j,p_nh4_a02) + &
chem(i,k,j,p_nh4_a03) + &
chem(i,k,j,p_nh4_a04)) * atm_mass_den + &
(chem(i,k,j,p_no3_a01) + &
chem(i,k,j,p_no3_a02) + &
chem(i,k,j,p_no3_a03) + &
chem(i,k,j,p_no3_a04)) * atm_mass_den
aso4(kp1) = aso4(kp1) + &
(chem(i,k,j,p_so4_a01) + &
chem(i,k,j,p_so4_a02) + &
chem(i,k,j,p_so4_a03) + &
chem(i,k,j,p_so4_a04)) * atm_mass_den
asal(kp1) = asal(kp1) + &
(chem(i,k,j,p_na_a01) + &
chem(i,k,j,p_na_a02) + &
chem(i,k,j,p_na_a03) + &
chem(i,k,j,p_na_a04)) * atm_mass_den + &
(chem(i,k,j,p_cl_a01) + &
chem(i,k,j,p_cl_a02) + &
chem(i,k,j,p_cl_a03) + &
chem(i,k,j,p_cl_a04)) * atm_mass_den
Extended MOZART-4 chemistry
- explicit description of aromatics (instead of lumped "TOLUENE")
- ethyne as explicit species
- updated description of glyoxal chemistry
see upcoming paper (Knote et al., in preparation).
Glyoxal SOA description
see upcoming paper (Knote et al., in preparation).
Emissions of particulate NH4, Na, Cl as given in AQMEII dataset into MOSAIC
module_mosaic_addemiss.F, mosaic_addemiss:
! CK 20130206 emit more stuff (AQMEII)
chem_select_2 : SELECT CASE( config_flags%chem_opt )
CASE(MOZART_MOSAIC_4BIN_VBS0_KPP, MOZART_MOSAIC_4BIN_VBS0_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
Updates to dust_opt == 3 for MOSAIC aerosols with 4 bins for correct size distribution
module_mosaic_addemiss.F, mosaic_dust_gocartemis
! CK 20130304 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
if (n <= 5) densdust=2.5
if (n > 5 ) densdust=2.65
! CK 20130304 added option for 4bin WRF
if (nsize_aer(itype) == 4) then
if (n <= 2) densdust=2.5
if (n > 2 ) densdust=2.65
endif
Wet scavenging of aerosols
The implementation of convective wet scavenging in module_ctrans_grell.F is used, adapted for the 4-bin MOSAIC aerosol scheme. Aqueous-phase chemistry during convective transport is considered according to the CMAQ (Walcek and Taylor) method (module_ctrans_aqchem.F). Large-scale wet scavenging uses the version of MOSAIC with an explicit aqueous phase (module_mosaic_wetscav.F), again adapted for 4-bin MOSAIC and additional SOA species. Aqueous-phase chemistry according to the Carnegie Mellon University (S. Pandis et al.) scheme is active for species scavenged in large-scale clouds (module_mosaic_cloudchem.F).
(— adaptions in the code to follow —)
Reading in emissions at first time step after restart (courtesy of Alessandra Balzarini, RSE)
#ifdef WRF_CHEM
! - Get chemistry data
ELSE IF( ialarm .EQ. AUXINPUT5_ALARM .AND. config_flags%chem_opt > 0 ) THEN
IF( config_flags%emiss_inpt_opt /= 0 ) THEN
! CK 20130315 fix for reading in emissions immediately after restart
IF( WRFU_AlarmIsRinging( grid%alarms( ialarm ), rc=rc ) .OR. ( &
(config_flags%restart) .AND. ( currTime .EQ. startTime ))) THEN
call wrf_debug(15,' CALL med_read_wrf_chem_emiss ')
CALL med_read_wrf_chem_emiss ( grid , config_flags )
VOC reactivity
Our definition of VOC reactivity is found here. We track CO separately. It is calculated online in KPP as follows:
- 2 new species in the chem array are created, called "vocr_nmvoc" and "vocr_co"
state real vocr_nmvoc ikjftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "vocr_nmvoc" "total NMVOC reactivity against OH" "s-1"
state real vocr_co ikjftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "vocr_co" "CO reactivity against OH" "s-1"
- all reactions of VOC + OH create 1 vocr_nmvoc on the right hand side, here the example for CO+OH
{ 81:042 } CO+OH=HO2{+CO2}+VOCR_CO : usr8(temp, c_m) ;
- in chem_driver.f90, vocr_nmvoc is zeroed out before the call to kpp_mechanism_driver, and divided by ([OH] * dt) afterwards
! CK 20130318 VOC-OH reactivity
if( config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_AQ_KPP ) then
do k=kts,kte
do i=its,ite
do j=jts,jte
! zero out VOC reactivity, and save old OH
oh_old(i,k,j) = chem(i,k,j,p_ho)
chem(i,k,j,p_vocr_nmvoc) = 0.0
chem(i,k,j,p_vocr_co) = 0.0
enddo
enddo
enddo
endif
CALL kpp_mechanism_driver (chem, &
grid%id,dtstepc,config_flags, &
p_phy,t_phy,rho, &
#if (NMM_CORE == 1)
moist_trans, &
#endif
#if (EM_CORE == 1)
moist, &
#endif
vdrog3, ldrog, vdrog3_vbs, ldrog_vbs, &
!
#include <call_to_kpp_mech_drive.inc>
!
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite,jts,jte,kts,kte)
! CK 20130318 VOC-OH reactivity
if( config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_AQ_KPP ) then
do k=kts,kte
do i=its,ite
do j=jts,jte
! vocr_* is calculated in KPP as product from all bimolecular reactions of VOC + OH:
! e.g. CH2O+OH=CO+HO2+H2O+vocr_nmvoc : 9.e-12_dp ;
! that is k * [CH2O] * [OH], and vocr_* is in "ppmv" (created per timestep)
! we calculate reactivity (s-1) back by
! vocr_* / ([OH] * dt) with [OH] being the concentration of OH before KPP
! ignore [OH] < 1e-3 pptv
if (oh_old(i,k,j) > 1.0e-9) then
chem(i,k,j,p_vocr_nmvoc) = chem(i,k,j,p_vocr_nmvoc) / (oh_old(i,k,j) * dtstepc)
chem(i,k,j,p_vocr_co) = chem(i,k,j,p_vocr_co) / (oh_old(i,k,j) * dtstepc)
endif
enddo
enddo
enddo
endif
VOC reactivity in general is the sum of all VOCs (A, B, C) times their reaction rates (a, b, c) of with OH:
a *A *OH + b * B * OH + c * C * OH
equivalent to
OH * ( a * A + b * B + c * C)
KPP calculates vocr_nmvoc as OH * ( a * A + b * B + c * C) * dt, hence we divide by (OH * dt) afterwards.
fTUV problems with some photolysis rates
fTUV has been found to overestimate photolysis rates of glyoxal and methylglyoxal, and some other species. Due to the lack of time, and in the light of the possibility that this will be a problem with the way fTUV has been developed rather than a bug, we scale some rates to measurements made during CalNex on the NOAA WP-3D aircraft.
This is a comparison of photolysis rates between WRF/Chem 3.4.1 (MOZART, MOSAIC 4-bin, RRTMG long and shortwave, fTUV, 4km resolution) and measurements by Harald Stark onboard the NOAA WP-3D during its flight on 06/14/10 (CalNex). Each plot shows photolysis rates in min-1 on the y-, and zenith angle on the x-axis. Note that while scaling of MEK would require a factor of 25 to match measurements, we only scale by a factor of 10, as measurements seem overestimated, and a reference TUV run indicated lower values as well.
The flight covered most of the Californian coast between Sacramento and LA, and also a broad range of zenith angles:
adapted from http://www.esrl.noaa.gov/csd/groups/csd7/measurements/2010calnex/P3/flighttrack/# (last accessed Feb 2013)
This is the resulting code in module_ftuv_driver.f90:
! CK 20130304 added scaling factors based on measurements
! on NOAA WP-3D flights over CalNex
! !! these are specific for a certain model configuration and
! !! should not be used for other runs
do k = kts, kte
ph_no2(i,k,j) = max( 0._dp,prate(k,pid_no2)*m2s )
ph_o31d(i,k,j) = max( 0._dp,prate(k,pid_o31d)*m2s )
ph_o33p(i,k,j) = max( 0._dp,prate(k,pid_o33p)*m2s )
ph_hono(i,k,j) = max( 0._dp,prate(k,pid_hono)*m2s )
ph_hno3(i,k,j) = max( 0._dp,prate(k,pid_hno3)*m2s )
ph_hno4(i,k,j) = max( 0._dp,prate(k,pid_hno4)*m2s )
ph_no3o2(i,k,j) = max( 0._dp,prate(k,pid_no3o2)*m2s )
ph_no3o(i,k,j) = max( 0._dp,prate(k,pid_no3o)*m2s*ph_no3o_factor )
! CK 20130304 rate scaling
ph_h2o2(i,k,j) = max( 0._dp,prate(k,pid_h2o2)*m2s ) * 1.4
ph_ch2om(i,k,j) = max( 0._dp,prate(k,pid_ch2om)*m2s )
ph_ch2or(i,k,j) = max( 0._dp,prate(k,pid_ch2or)*m2s )
ph_ald(i,k,j) = max( 0._dp,prate(k,pid_ald)*m2s )
ph_op1(i,k,j) = max( 0._dp,prate(k,pid_op1)*m2s )
ph_op2(i,k,j) = max( 0._dp,prate(k,pid_op2)*m2s )
ph_paa(i,k,j) = max( 0._dp,prate(k,pid_paa)*m2s )
ph_ket(i,k,j) = max( 0._dp,prate(k,pid_ket)*m2s )
ph_glya(i,k,j) = max( 0._dp,prate(k,pid_glya)*m2s )
ph_glyb(i,k,j) = max( 0._dp,prate(k,pid_glyb)*m2s )
! CK 20130304 rate scaling
ph_mgly(i,k,j) = max( 0._dp,prate(k,pid_mgly)*m2s ) * 0.2
ph_dcb(i,k,j) = max( 0._dp,prate(k,pid_dcb)*m2s )
ph_onit(i,k,j) = max( 0._dp,prate(k,pid_onit)*m2s )
ph_macr(i,k,j) = max( 0._dp,prate(k,pid_macr) *m2s )
ph_ch3coc2h5(i,k,j) = max( 0._dp,prate(k,23)*m2s )
enddo
do m = 1,nw-1
do k = 1,nz
ph_radfld(i,k,j,m) = radfld(k,m)
end do
end do
do m = 1,tuv_jmax
do k = 1,nz
ph_adjcoe(i,k,j,m) = adjcoe(k,m)
ph_prate (i,k,j,m) = prate0(k,m)
end do
end do
wc_x(1:nw-1) = wc(1:nw-1)
zref_x(1:nref) = zref(1:nref)
! CK 20130118 added aqueous-phase MOZART-MOSAIC-VBS0 version
if( config_flags%chem_opt == MOZART_KPP .or. &
config_flags%chem_opt == MOZCART_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_KPP .or. &
config_flags%chem_opt == MOZART_MOSAIC_4BIN_VBS0_AQ_KPP ) then
do k = kts, kte
ph_n2o(i,k,j) = max( 0._dp,prate(k,pid_n2o)*m2s )
ph_n2o5(i,k,j) = max( 0._dp,prate(k,pid_n2o5)*m2s )
ph_mvk(i,k,j) = max( 0._dp,prate(k,pid_mvk)*m2s )
ph_pan(i,k,j) = max( 0._dp,prate(k,pid_pan)*m2s )
ph_ald(i,k,j) = max( 0._dp,(prate(k,pid_ch3choa)+prate(k,pid_ch3chob)+prate(k,pid_ch3choc))*m2s )
! CK 20130304 rate scaling
ph_gly(i,k,j) = max( 0._dp,prate(k,pid_mgly)*m2s ) * 0.14
! CK 20130304 rate scaling
!ph_mek(i,k,j) = ph_ket(i,k,j)
ph_mek(i,k,j) = max( 0._dp,prate(k,pid_ket)*m2s ) * 10.0
ph_acetol(i,k,j) = ph_op1(i,k,j)
ph_pooh(i,k,j) = ph_op1(i,k,j)
ph_glyald(i,k,j) = max( 0._dp,prate(k,pid_glyald) *m2s )
ph_hyac(i,k,j) = max( 0._dp,2._dp*prate(k,pid_hyac) *m2s )
! CK 20130304 rate scaling
!ph_hno4(i,k,j) = ph_hno4(i,k,j) + 1.e-5_dp*m2s
ph_hno4(i,k,j) = max( 0._dp,prate(k,pid_hno4)*m2s )
enddo
end if