WRF/Chem at NCAR

We use WRF/Chem 3.4.1 with modifications. Our setup options are given in the overview table here. Output variables are listed here.

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:

  1. 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"
  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) ;
  1. 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.

ftuv_noaawp3d_comparison

The flight covered most of the Californian coast between Sacramento and LA, and also a broad range of zenith angles:

calnex_wp3d_flighttrack_20100614
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
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License