Offline CMAQ Dust Emulator

README.md

Overview

---

title: CMAQ Dust Emulator
author: Barron H. Henderson
organization: US EPA
last updated: 2025-11-03

---

github.com/USEPA/CMAQ/emis/emis/DUST_EMIS.F routines and functions converted from Fortran to Python. All modules are replaced by inline code and direct access to MCIP files. The original DUST_EMIS.F was developed by Foroutan et al. (https://dx.doi.org/10.1002/2016MS0008230).

How To

Steps

1. Create a configuration file (e.g., mp2022_dust.json)
2. Run cmaqdust
3. Plot

Each step is detailed below.

Create a Configuration File

The example below calculates hourly emissions for 2022 and saves them out in month-long files. The inputs are mcip using the strftime templates. For example, when the date is 2022-07-31 then %Y=2022, %y=22, %m=07, %d=31.

{
    "start": "2022-01-01T00", "end": "2022-12-31T23", "freq": "1h",
    "luftmpl": "inputs/mcip/LUFRAC_CRO.12US1.35L.%y%m%d",
    "g2dtmpl": "inputs/mcip/GRIDCRO2D.12US1.35L.%y%m%d",
    "m2dtmpl": "inputs/mcip/METCRO2D.12US1.35L.%y%m%d",
    "m3dtmpl": "inputs/mcip/METCRO3D.12US1.35L.%y%m%d",
    "vegtmpl": "inputs/mcip/VEG.12US1.35L.%y%m%d",
    "outtmpl": "outputs/dust_emis.12US1.%Y-%m.nc",
    "n_jobs": 3
}

Run cmaqdust.py

You must install prerequisites the first time (see below) and then run. Every other time, you can just run.

python cmaqdust.py mp2022_dust.json

Plot

python plot_dust.py

Prerequisites:

• Python version>=3.9
• Python libraries
  • pandas
  • xarray
  • numpy
  • joblib
  • (otional) pycno for adding maps
• MCIP meteorological outputs
  • required: GRIDCRO2D, METCRO2D, LUFRAC_CRO
  • optional: METCRO3D, VEG

To install prereqs:

All prerequisites are easily installed with pip as shown below.

python -m pip install numpy pandas xarray joblib

To download example inputs:

The 2022 modeling platform is avaialable on AWS. You can download it using wget or your favorite download tool.

mkdir -p inputs/mcip
wget -c -P inputs/mcip https://epa-2022-modeling-platform.s3.amazonaws.com/MCIP/12US1/GRIDCRO2D.12US1.35L.220701
wget -c -P inputs/mcip https://epa-2022-modeling-platform.s3.amazonaws.com/MCIP/12US1/METCRO2D.12US1.35L.220701
wget -c -P inputs/mcip https://epa-2022-modeling-platform.s3.amazonaws.com/MCIP/12US1/LUFRAC_CRO.12US1.35L.220701
# skip METCRO3D because it is big; use DENS = PRSFC / TEMP2 / R * (1 - Q2 * MWAIR / MWWAT) * MWAIR
# wget -c -P inputs/mcip https://epa-2022-modeling-platform.s3.amazonaws.com/MCIP/12US1/METCRO3D.12US1.35L.220701

Change Log

2025-11-03: some basic file renaming and restructuring. 2025-10-30: updated to match CMAQ v5.5+

Annotated Directory

.
|-- README.md  # this file
|-- cmaqdust.py  # The CMAQ dust logic
|-- plot_dust.py  # Example plotting routine
|-- inputs  # not distributed
|   `-- mcip/
|       |-- GRIDCRO2D.12US1.35L.%y%m%d
|       |-- LUFRAC_CRO.12US1.35L.%y%m%d
|       |-- METCRO2D.12US1.35L.%y%m%d
|       |-- METCRO3D.12US1.35L.%y%m%d  (optional, used in example)
|       `-- VEG.12US1.35L.%y%m%d       (optional, not used in example)
`-- outputs/
    |-- dust_emis.12US1.%Y-%m.nc
    `-- dust_emis.12US1.%Y-%m_DUST_EMIS.png

cmaqdust.py

__doc__ = """
Offline CMAQ Python WBD approximator
====================================

github.com/USEPA/CMAQ/emis/emis/DUST_EMIS.F routines and functions converted
from Fortran to Python. All modules are replaced by inline code and direct
access to MCIP files. The original DUST_EMIS.F was developed by Foroutan et
al. (https://dx.doi.org/10.1002/2016MS0008230).

updated: 2025-10-30 to match CMAQ v5.5+

To-Do:
1. Currently only operational for NLCD40 and MODIS_NOAH.
  a. Added MODIS/NLCD, but not tested.
2. Currently overestimates compared to comparable run
  a. Joey's 2022 run had 28 Tg in April.
  b. 2022 run has 29.6 Tg in April
  c. May be within uncertainty of eyeballing the plot to get 28Tg


Contents
--------

process : function
    Apply and output results for a desired time range.

get_dust_emis : function
    Calculate dust emissions

dust_hflux : function
    Used by get_dust_emis

Example
-------

### Process dates

```python
import cmaqdust
config = {
    'start': '2022-07-01', 'end': '2022-07-01T23', 'freq': '1h',
    'g2dtmpl': 'inputs/mcip/GRIDCRO2D.12US1.35L.%y%m%d',
    'm2dtmpl': 'inputs/mcip/METCRO2D.12US1.35L.%y%m%d',
    'm3dtmpl': 'inputs/mcip/METCRO3D.12US1.35L.%y%m%d',
    'luftmpl': 'inputs/mcip/LUFRAC_CRO.12US1.35L.%y%m%d',
    'outtmpl': 'outputs/dust_emis.12US1.%Y-%m-%d.nc'
}
cmaqdust.process(**config)
```

"""
__version__ = '1.1.0'
__cmaqversion__ = 'v5.5+'
import warnings
import pandas as pd
import xarray as xr
import numpy as np
import os

warnings.simplefilter('once')
_dwmsg = 'No METCRO3D; set DENS = PRSFC / TEMP2 / R'
_dwmsg += ' * (1 - Q2 * MWAIR / MWWAT) * MWAIR'
_vwmsg = 'No VEG; using VEG from METCRO2D'
warnings.filterwarnings('once', message=_dwmsg)
warnings.filterwarnings('once', message=_vwmsg)


def get_dust_emis(
    jdate, jtime, rjacm, cellhgt, lufpath, g2dpath, m2dpath, m3dpath,
    vegpath=None
):
    # Met_Data%
    lud = xr.open_dataset(lufpath)
    g2d = xr.open_dataset(g2dpath)
    m2d = xr.open_dataset(m2dpath)
    if m3dpath is None:
        # create an approximation of cell center density
        # constants copied on 2025-11-03 from:
        # https://github.com/USEPA/CMAQ/blob/main/CCTM/src/ICL/fixed/const/CONST.EXT
        # updated g/mol to kg/mol
        MWAIR = .0289628  # kg/mol
        RGASUNIV = 8.314510  # m3 Pa K-1 mol-1
        MWWAT = .0180153  # kg/mol
        m3d = m2d[[]]
        X2 = m2d['Q2'] / MWWAT * MWAIR  # kg_w / kw_a -> mol_w / mol_a
        m3d['DENS'] = (
            m2d['PRSFC'] / m2d['TEMP2'] / RGASUNIV  # mol_air/m3
        ) * (1 - X2) * MWAIR  # kg_dryair / m3
    else:
        m3d = xr.open_dataset(m3dpath)
    if vegpath is None:
        vegd = m2d
    else:
        vegd = xr.open_dataset(vegpath)
    ti = jtime // 10000
    sltyp = m2d['SLTYP'][ti, 0].astype('i')
    lufrac = lud['LUFRAC'][0]
    soim1 = m2d['SOIM1'][ti, 0]
    dens1 = m3d['DENS'][ti, 0]
    rc = m2d['RC'][ti, 0]
    rn = m2d['RN'][ti, 0]
    rain = rc + rn
    snocov = m2d['SNOCOV'][ti, 0]
    WSPD10 = m2d['WSPD10'][ti, 0]
    lwmask = g2d['LWMASK'][0, 0]
    clay = m2d['CLAY_PX'][ti, 0]
    csand = m2d['CSAND_PX'][ti, 0]
    fmsand = m2d['FMSAND_PX'][ti, 0]
    veg = vegd['VEG'][ti, 0]

    nr = rc.shape[0]
    nc = rc.shape[1]

    ndp = 4     # number of soil texture type particle sizes:
    #  1  Coarse sand
    #  2  Fine-medium sand
    #  3  Silt
    #  4  Clay
    karman = 0.40      # von Karman constant
    mv = 0.16
    sigv = 1.45
    betav = 202.0
    sigv_mv = sigv * mv   # = 0.232
    betav_mv = betav * mv  # = 32.32
    mb = 0.5
    sigb = 1.0
    betab = 90.0
    sigb_mb = sigb * mb   # = 0.5
    betab_mb = betab * mb  # = 45.0
    vnmld_nlcd_modis = [
        ('LUFRAC_37', 'Closed Shrublands           ', 37),
        ('LUFRAC_38', 'Open Shrublands             ', 38),
        ('LUFRAC_07', 'Barren Land (Rock-Sand-Clay)', 7),
        ('LUFRAC_12', 'Dwarf Scrub                 ', 12),
        ('LUFRAC_13', 'Shrub-Scrub                 ', 13),
        ('LUFRAC_18', 'Tundra                      ', 18),
        ('LUFRAC_47', 'Barren or Sparsely Vegetated', 47)
    ]
    # all decremented by 1 because 0-based
    # updated to match https://github.com/USEPA/CMAQ_Dev/pull/1181
    dmap_nlcd_modis = np.array([  # land use type desert map to BELD3
        0,       # shrubland
        0,       # shrubland
        2,       # barrenland
        0,       # updated from barrenland (2) to shrubland (0) per PR 1181
        0,       # updated from barrenland (2) to shrubland (0) per PR 1181
        2,       # barrenland
        2,       # barrenland
        2,       # ag landuse surrogate
    ])
    vnmld_nlcd40 = [
        ('LUFRAC_06', 'Closed Shrublands           ', 6),
        ('LUFRAC_07', 'Open Shrublands             ', 7),
        ('LUFRAC_27', 'Barren Land (Rock-Sand-Clay)', 27),
        ('LUFRAC_31', 'Dwarf Scrub                 ', 31),
        ('LUFRAC_32', 'Shrub-Scrub                 ', 32),
        ('LUFRAC_16', 'Barren or Sparsely Vegetated', 16),
    ]
    # all decremented by 1 because 0-based
    # updated to match https://github.com/USEPA/CMAQ_Dev/pull/1181
    dmap_nlcd40 = np.array([  # land use type desert map to BELD3
        0,       # shrubland
        0,       # shrubland
        2,       # barrenland
        2,       # barrenland
        0,       # updated from barrenland (2) to shrubland (0) per PR 1181
        0,       # updated from barrenland (2) to shrubland (0) per PR 1181
        2,       # ag landuse surrogate
    ])
    vnmld_modis_noah = [
        ('LUFRAC_06', 'Closed Shrublands           ', 6),
        ('LUFRAC_07', 'Open Shrublands             ', 7),
        ('LUFRAC_16', 'Barren or Sparsely Vegetated', 16),
        ('LUFRAC_20', 'Barren Tundra               ', 20)
    ]
    dmap_modis_noah = np.array([  # land use type desert map to BELD3
        0,       # shrubland
        0,       # shrubland
        2,       # barrenland
        2,       # barrenland
        2        # ag landuse surrogate
    ])
    vnmld_beld3 = [
        ('LUFRAC_??', 'USGS_shrubland ', 1),
        ('LUFRAC_??', 'USGS_shrubgrass', 2),
        ('LUFRAC_??', 'USGS_sprsbarren', 3),
    ]
    dmap_beld3 = np.array([  # land use type desert map to BELD3
        0,       # shrubland
        1,       # shrubgras
        2,       # barrenland
        2        # ag landuse surrogate
    ])
    if lufrac.shape[0] == 40:
        vnmld = vnmld_nlcd40
        dmap = dmap_nlcd40
    elif lufrac.shape[0] == 50:
        vnmld = vnmld_nlcd_modis
        dmap = dmap_nlcd_modis
    elif lufrac.shape[0] == 20:
        vnmld = vnmld_modis_noah
        dmap = dmap_modis_noah
    elif lufrac.shape[0] == 3:
        vnmld = vnmld_beld3
        dmap = dmap_beld3
    else:
        print('WARNING: Unknown number of LUFRAC; must be 40 or 20')
        vnmld = vnmld_nlcd40
        dmap = dmap_nlcd40
    n_dlcat = dmap.shape[0] - 1
    emap = np.zeros((n_dlcat + 1), dtype='i')
    # ---Roughness density for solid elements
    # from Darmenova et al. [JGR,2009] and Xi and Sokolik [JGR,2015]
    # shrubland, shrubgrass, barrenland, cropland
    lambdab = np.array([0.03, 0.04, 0.0001, 0.15])
    # ---Compound for computational efficiency
    hb_lambdab = np.array([6.0e-04, 8.0e-04, 2.0e-06, 3.0e-03])
    eropot = np.array([    # erodible potential of soil components
        0.08,   # clay
        1.00,   # silt
        0.12    # sand
    ])
    # converted to gravimetric [kg/kg]
    soilml1 = np.array([
        0.242,     # Sand
        0.257,     # Loamy Sand
        0.286,     # Sandy Loam
        0.350,     # Silt Loam
        0.350,     # Silt
        0.307,     # Loam
        0.277,     # Sandy Clay Loam
        0.350,     # Silty Clay Loam
        0.332,     # Clay Loam
        0.284,     # Sandy Clay
        0.357,     # Silty Clay
        0.344,     # Clay
        0.329,     # Organic Material
        0.000,     # Water
        0.170,     # BedRock
        0.280      # Other
    ])

    # ---Soil texture: the amount of
    #   1: Coarse sand, 2: Fine-medium sand, 3: Silt, 4: Clay
    #   in each soil type [Kg/Kg]. from Menut et al. [JGR,2013]
    # soiltxt = np.array([
    #     [0.46, 0.46, 0.05, 0.03],     # Sand
    #     [0.41, 0.41, 0.18, 0.00],     # Loamy Sand
    #     [0.29, 0.29, 0.32, 0.10],     # Sandy Loam
    #     [0.00, 0.17, 0.70, 0.13],     # Silt Loam
    #     [0.00, 0.10, 0.85, 0.05],     # Silt
    #     [0.00, 0.43, 0.39, 0.18],     # Loam
    #     [0.29, 0.29, 0.15, 0.27],     # Sandy Clay Loam
    #     [0.00, 0.10, 0.56, 0.34],     # Silty Clay Loam
    #     [0.00, 0.32, 0.34, 0.34],     # Clay Loam
    #     [0.00, 0.52, 0.06, 0.42],     # Sandy Clay
    #     [0.00, 0.06, 0.47, 0.47],     # Silty Clay
    #     [0.00, 0.22, 0.20, 0.58],     # Clay
    #     [0.00, 0.00, 0.00, 0.00],     # Organic Material
    #     [0.00, 0.00, 0.00, 0.00],     # Water
    #     [0.00, 0.00, 0.00, 0.00],     # BedRock
    #     [0.00, 0.00, 0.00, 0.00]
    # ])

    # ---Mean mass median particle diameter (m) for each soil texture type
    #   Chatenet et al. [Sedimentology,1996] and Menut et al. [JGR,2013]
    dp = np.array([
        690.0E-6,     # Coarse sand
        210.0E-6,     # Fine-medium sand
        125.0E-6,     # Silt
        2.0E-6        # Clay
    ])

    # ---Get Julian day number in year
    jday = float(jdate % 1000)

    # ---Vegetation height dynamically changed based on the month of the year
    #   Veg. heights in [m] for 1: Shrubland 2: shrubgrass 3: barrenland
    #   4: Cropland following the idea of Xi and Sokolik [JGR,2015]
    if (jday > 59 and jday <= 90):         # Mar
        hv = np.array([0.15, 0.05, 0.10, 0.05])
    elif (jday > 90 and jday <= 120):   # Apr
        hv = np.array([0.15, 0.10, 0.10, 0.05])
    elif (jday > 120 and jday <= 151):  # May
        hv = np.array([0.12, 0.20, 0.10, 0.10])
    elif (jday > 151 and jday <= 181):  # Jun
        hv = np.array([0.12, 0.15, 0.10, 0.30])
    elif (jday > 181 and jday <= 212):  # Jul
        hv = np.array([0.10, 0.12, 0.10, 0.50])
    elif (jday > 212 and jday <= 243):  # Aug
        hv = np.array([0.10, 0.12, 0.10, 0.50])
    elif (jday > 243 and jday <= 273):  # Sep
        hv = np.array([0.10, 0.10, 0.10, 0.30])
    elif (jday > 273 and jday <= 304):  # Oct
        hv = np.array([0.05, 0.08, 0.10, 0.10])
    else:                                 # Nov-Feb
        hv = np.array([0.05, 0.05, 0.05, 0.05])

    # set erodible landuse map
    for m in range(n_dlcat):
        emap[m] = dmap[m]  # dmap maps to one of the 3 BELD3 desert types

    emap[m + 1] = 4

    # --------- ###### Start Main Loop ###### ---------
    dust_em = np.zeros((nr, nc), dtype='f')
    soimt = np.zeros_like(dust_em)
    fmoit = np.zeros_like(dust_em)
    vegfrac = np.zeros_like(dust_em)
    ustr = np.zeros((n_dlcat, nr, nc), dtype='f')
    qam = np.zeros_like(ustr)
    elus = np.zeros_like(ustr)
    fruf = np.zeros_like(ustr)
    kvh = np.zeros_like(ustr)

    rlay1hgt = rjacm / cellhgt

# --- Set Clay, coarse and fine/medium sand fractions.
# --- If value from WRF is missing (-9999.) use old table values
# --- If value from WRF is from WRFV4.1 PX LSM csand_px, etc use those
    sandf = csand + fmsand
    siltf = 1.0 - clay - sandf
    soiltxt_gcell = np.stack([
        csand, fmsand, siltf, clay
    ])

# -- Vegetation fraction based on the WRF/MCIP VEG variable. In WRF that would
# -- be VEGF_PX for the case of PX and VEGFRA in the case of other LSMs. In
# -- more recent WRFv4+ versions high resolution MODIS veg data is availiable
# -- and can be used in PX with pxlsm_modis_veg = 1
    # read from METCRO2D
    vegfrac = np.maximum(np.minimum(veg, 0.95), 0.005)
    vegfree = 1.0 - vegfrac
    lambdav = -0.35 * np.log(vegfree)  # Shao et al. [Aus. J. Soil Res.,1996]

# -- Dust possiblity only if 1. not over water
#                           2. rain < 1/100 in. (1 in. = 2.540 cm)
#                           3. not snow-covered
#                           4. if soimt <= limit
#                           5. desert type or ag landuse
#                           6. erodible landuse
#                           7. friction velocity > threshold

# Calculate maximum amount of the adsorbed water
#    w` = 0.0014(%clay)**2 + 0.17(%clay) - w` in %
#    Fecan et al. [1999,Annales Geophys.,17,144-157]
    wmax = (14.0 * clay + 17.0) * clay   # [%]

# Change soil moisture units from volumetric (m**3/m**3) to gravimetric (Kg/Kg)
    soimt = soim1 * 1000.0 / (2650.0 * (0.511 + 0.126 * sandf))
# -- Dust possiblity 1,2,3 4
    iszero = (
        (lwmask == 0.0) | (rain > 0.0254) | (snocov >= 0.001)
        | (soilml1[sltyp] < soimt)
    )
    fmoit = np.where(
        soimt <= (0.01 * wmax),
        1,
        np.sqrt(1.0 + 1.21 * (100.0 * soimt - wmax)**0.68)
    )

# -- Erodibility potential of soil component
    sd_ep = clay * eropot[0] + siltf * eropot[1] + sandf * eropot[2]

# -- Lu and Shao [JGR,1999] and Kang et al. [JGR,2011]
#   First, mapping soil types into 4 main soil types following Kang et al.
#   [JGR,2011]
    flxcondlist = [
        np.isin(sltyp, [1, 2]),
        np.isin(sltyp, [3, 4, 6, 8, 9]),
        np.isin(sltyp, [7]),
        np.isin(sltyp, [5, 10, 11, 12]),
    ]
    flxfac1list = [
        5.886e-05,
        5.2974e-05,
        9.4176e-05,
        2.3544e-05,
    ]
    flxfac2list = [
        1.5215430,
        1.0758933,
        1.0758933,
        0.1964303,
    ]
    flxfac1 = np.select(flxcondlist, flxfac1list, default=0)
    flxfac2 = np.select(flxcondlist, flxfac2list, default=0.3402273)

    # ------- Start Loop Over Erodible Landuse ----
    for ei, (lvkey, lname, luidx) in enumerate(vnmld):
        m = ei
        n = emap[ei]
        # lufrac is in units of 0-1, so added 100 multiplier
        elus[m] = lufrac[luidx - 1] * vegfree * 100
        rlambda = lambdab[n] + lambdav
        vegheight = (hb_lambdab[n] + hv[n] * lambdav) / rlambda

# -- New parametrization for surface roughness by H. Foroutan - Oct. 2015
        z0 = np.where(
            rlambda <= 0.2,
            0.96 * (rlambda ** 1.07) * vegheight,
            0.083 * (rlambda ** (-0.46)) * vegheight,
        )

# -- Calculate friction velocity (U*) at the surafce applicable to dust
# -- emission
        ustr[m] = karman * WSPD10 / np.log(10.0 / z0)

# -- Roughness effect on U*t (Drag partitioning)
#   Xi and Sokolik [JGR,2015]
        fruf2 = (
            (1.0 - sigv_mv * lambdav)
            * (1.0 + betav_mv * lambdav)
            * (1.0 - sigb_mb * lambdab[n] / vegfree)
            * (1.0 + betab_mb * lambdab[n] / vegfree)
        )
        fruf[m] = np.where(fruf2 > 1.0, np.sqrt(fruf2), 10)

# -- Vert-to-Horiz dust flux ratio : Kang et al. [JGR, 2011] : Eq. (12)
        kvh[m] = flxfac1 * (0.24 + flxfac2 * ustr[m])

        hflux = dust_hflux(
            ndp, dp, soiltxt_gcell,
            fmoit, fruf[m], ustr[m], sd_ep, dens1
        )
        vflux = hflux * kvh[m]                         # [g/m**2/s]
        qam[m] += vflux * rlay1hgt * (elus[m] * 0.01)  # [g/m**3/s]
        dust_em[:] += qam[m]
    dust_em = np.where(iszero, 0, dust_em)
    return dust_em


def dust_hflux(ndp, dp, soiltxt, fmoit, fruf, ustr, sd_ep, dens):
    """
    Arguments
    ---------
    ndp : int
        Number of particle size bins for csand, fmsand, siltf, clay
    dp : array
        Array of particle diameters for csand, fmsand, siltf, clay
    soiltxt : array
        soil texture fraction (ndp, nr, nc) where data has fraction
    fmoit : array
        Moisture factor

    C usage: hflux = dust_flux( ndp, dp,
    C                           soiltxt( : ),
    C                           fmoit( c,r ),
    C                           fruf( c,r,m ),
    C                           ustr( c,r,m ),
    C                           sd_ep( c,r ),
    C                           dens( c,r ) )
    """
    grav = 9.80622
    amen = 1.0          # Marticorena and Bergametti [JGR,1997]
    cfac = 1000.0 * amen / grav
    A = 260.60061       # 0.0123 * 2650.0 * 9.81 / 1.227
    B = 1.6540342e-06   # 0.0123 * 0.000165 / 1.227
    # real utstar                            ! threshold U* [m/s]
    # real utem                              ! U term [(m/s)**3]
    # real fac
    # integer n

    fac = cfac * dens * sd_ep
    utem = 0.0
    utstar = 0.0
    hflux = 0.0
    for n in range(ndp):   # loop over dust particle size
        # Shao and Lu [JGR,2000]
        utstar = np.sqrt(A * dp[n] + B / dp[n]) * fmoit * fruf
        # ---Horiz. Flux from White (1979)
        utem = (ustr + utstar) * (ustr * ustr - utstar * utstar)
        # ---Horiz. Flux from Owen (1964)
        hfluxtmp = fac * utem * soiltxt[n]   # [g/m/s]
        # wind erosion occurs only if U* > U*t
        hfluxtmp = np.where(ustr > utstar, hfluxtmp, 0)
        hflux += hfluxtmp
    return hflux


def process(
    start, end, freq='1h',
    luftmpl='inputs/mcip/LUFRAC_CRO.12US1.35L.%y%m%d',
    g2dtmpl='inputs/mcip/GRIDCRO2D.12US1.35L.%y%m%d',
    m2dtmpl='inputs/mcip/METCRO2D.12US1.35L.%y%m%d',
    m3dtmpl='inputs/mcip/METCRO3D.12US1.35L.%y%m%d',
    vegtmpl='inputs/mcip/VEG.12US1.35L.%y%m%d',
    outtmpl='outputs/dust_emis.12US1.%Y-%m-%d.nc',
    n_jobs=1
):
    """
    Process CMAQ dust for dates from start to end (inclusive) in freq steps.
    Outputs are collected in files whose output frequency depends on the
    strftime frequency of the output template.

    Arguments
    ---------
    start : str or date-like
        Starting date in UTC
    end : str or date-like
        Ending date in UTC
    freq : str
       Defaults to 1h, but could be 2h, 3h, 4h, 6h
    luftmpl : str
        Path to landuse fraction on grid centers (CRO) in strftime format.
        Must have LUFRAC variable.
    g2dtmpl : str
        Path to GRIDCRO2D file in strftime format. Must have MSFX2 and LWMASK.
    m2dtmpl : str
        Path to METCRO2D file in strftime format. Must have SLTYP, SOIM1, RC,
        RN, SNOCOV, WSPD10, CLAY_PX, CSAND_PX, FMSAND_PX variables.
    m3dtmpl : str
        Optional, path to METCRO3D file in strftime format. If provided, must
        have DENS variable. If not provided, DENS is calculated from METCRO2D:
        DENS = PRSFC / TEMP2 / R * (1 - Q2 * MWAIR / MWWAT) * MWAIR
    vegtmpl : str
        Optional, path to VEG file in strftime format. If provided, must have
        VEG variable.
    outtmpl : str
        Path strftime template for outputs. Must include a directory (./ if
        working directory)

    Returns
    -------
    outputs : list
        List of outputs.
    """
    from joblib import Parallel, delayed
    from os.path import basename
    dates = pd.date_range(start, end, freq=freq)
    dateoutpaths = dates.strftime(outtmpl)
    outpaths = []

    for outpath, mydates in dates.groupby(dateoutpaths).items():
        if os.path.exists(outpath):
            print(f'WARN:: Keeping cached {outpath}; delete to remake')
            outpaths.append(outpath)
            continue
        os.makedirs(os.path.dirname(outpath), exist_ok=True)
        jobs = []
        dust_ems = []
        for date in mydates:
            jdate = int(date.strftime('%Y%j'))
            jtime = int(date.strftime('%H%M%S'))
            lufpath = date.strftime(luftmpl)
            g2dpath = date.strftime(g2dtmpl)
            m2dpath = date.strftime(m2dtmpl)
            if m3dtmpl is None:
                m3dpath = None
            else:
                m3dpath = date.strftime(m3dtmpl)
                if not os.path.exists(m3dpath):
                    m3dpath = None
            if m3dpath is None:
                warnings.warn(_dwmsg)

            if vegtmpl is None:
                vegpath = None
            else:
                vegpath = date.strftime(vegtmpl)
                if not os.path.exists(vegpath):
                    vegpath = None
            if vegpath is None:
                warnings.warn(_vwmsg)
            # layerdepth is used in get_dust_emis to return concentration units
            # This assumes 1 / 20 as conversion from g/cm2/ to g/cm3
            layerdepth = 20.
            jobs.append(
                delayed(get_dust_emis)(
                    jdate, jtime, 1, layerdepth,
                    lufpath=lufpath, g2dpath=g2dpath, m2dpath=m2dpath,
                    m3dpath=m3dpath, vegpath=vegpath
                )
            )
        dust_ems = Parallel(n_jobs=n_jobs, verbose=1)(jobs)

        dust_ems = np.stack(dust_ems)
        nr, nc = dust_ems[0].shape
        dustvar = xr.DataArray(
            dust_ems[:, None, :, :] / 1e3 * layerdepth, name='dust_emis',
            dims=('TSTEP', 'LAY', 'ROW', 'COL'),
            attrs=dict(
                long_name='dust_emis'.ljust(16), units='kg/m**2/s'.ljust(16),
                description='Calculated assuming 20m layer depth',
                var_desc='Calculated assuming 20m layer depth'.ljust(80)
            ),
            coords={
                'TSTEP': mydates, 'LAY': np.array([0.9985]),
                'ROW': np.arange(nr) + 0.5, 'COL': np.arange(nc) + 0.5
            }
        )
        g2d = xr.open_dataset(g2dpath)
        msfx2 = g2d['MSFX2']
        area = g2d.XCELL * g2d.YCELL / msfx2
        dustf = xr.Dataset(
            {'dust_emis': dustvar}
        )
        dustf['AREA'] = (
            ('ROW', 'COL'), area[0, 0].data,
            {'units': 'm**2'.ljust(16), 'long_name': 'area'.ljust(16),
             'var_desc': 'area'.ljust(80)}
        )
        outf = dustf[[]]
        jd = mydates.strftime('%Y%j').astype('i')
        jh = mydates.strftime('%H%M%S').astype('i')
        if freq.endswith('h'):
            tstep = int(freq[:-1]) * 10000
        elif freq.endswith('d'):
            tstep = int(freq[:-1]) * 24 * 10000
        else:
            raise Exception('oops, h or d')
        tf = np.array([jd, jh]).T[:, None, :].repeat(2, 1)
        tattrs = dict(
            long_name='TFLAG'.ljust(16), units='<YYYYJJJ,HHMMSS>',
            var_desc='time'.ljust(80)
        )
        outf['TFLAG'] = xr.DataArray(
            tf, dims=('TSTEP', 'VAR', 'DATE-TIME'), attrs=tattrs
        )
        outf['DUST_EMIS'] = dustf['dust_emis']
        outf['AREA'] = dustf['AREA']
        outf.attrs.update(g2d.attrs)
        bn = basename(__file__)
        hist = f'Made by {bn} v{__version__} process'
        flines = [
            hist + ' with arguments:',
            (str(start) + ' ' + str(end) + f' by {freq}'),
            luftmpl, g2dtmpl, m2dtmpl, str(m3dtmpl), str(vegtmpl)
        ]
        fdesc = ''.join([_l.ljust(80)[:80] for _l in flines])
        fdesc = fdesc.ljust(4800)[:4800]
        outf.attrs.update(
            NVARS=2, TSTEP=tstep, SDATE=tf[0, 0, 0], STIME=tf[0, 0, 1],
            FILEDESC=fdesc, HISTORY=hist.ljust(4800)[:4800]
        )
        outf.attrs['VAR-LIST'] = 'DUST_EMIS'.ljust(16) + 'AREA'.ljust(16)
        idx = list(outf.indexes)
        outf = outf.drop_indexes(idx).reset_coords(drop=True)
        outf.to_netcdf(outpath, format='NETCDF4_CLASSIC')
        outpaths.append(outpath)
    return outpaths


def speciate(exprcfg, inpath, outpath):
    """
    Arguments
    ---------
    exprcfg : dict
        Keys are output variable names, values are dict that include expression
        (to be evaluated), units that must match the expression (note:
        DUST_EMIS is total dust in kg/m2/s and AREA is m2), and other
        attributes to be added to the variable as properites.
    inpath : str
        Path to input
    outpath : str
        Path to output

    Returns
    -------
    None
    """
    import json
    import xarray as xr
    import pprint
    from os.path import basename
    bn = basename(__file__)
    hlines = [
        f'Made by {bn} v{__version__} speciate',
    ] + pprint.pformat(exprcfg).split('\n')

    if isinstance(exprcfg, str):
        with open(exprcfg, 'r') as jf:
            exprcfg = json.load(jf)

    varnames = sorted(exprcfg)
    nvars = len(varnames)
    ef = xr.open_dataset(inpath).load()
    tf = ef['TFLAG'][:, 0].expand_dims(VAR=nvars).transpose(
        'TSTEP', 'VAR', 'DATE-TIME'
    )
    sf = ef[list(ef.coords)]
    sf['TFLAG'] = tf
    gbls = {'DUST_EMIS': ef['DUST_EMIS'], 'AREA': ef['AREA']}
    for k in varnames:
        opts = exprcfg[k]
        exprstr = opts['expression']
        sf[k] = eval(exprstr, None, gbls)
        sf[k].attrs.update(opts)
        sf[k].attrs.update(
            long_name=k.ljust(16),
            var_desc=exprstr.ljust(80)[:80]
        )

    sf.attrs.update(ef.attrs)
    hist = ef.attrs.get('HISTORY', '').strip()
    nhist = len(hist) // 80 + 1
    hist = hist.ljust(80 * nhist)
    hist += ''.join([_l.ljust(80)[:80] for _l in hlines])
    sf.attrs['VAR-LIST'] = ''.join([k.ljust(16) for k in varnames])
    sf.attrs['NVARS'] = nvars
    sf.attrs['HISTORY'] = hist.ljust(4800)[:4800]
    sf.to_netcdf(outpath, format='NETCDF4_CLASSIC')


if __name__ == '__main__':
    import json
    import argparse

    prsr = argparse.ArgumentParser()
    prsr.add_argument('config')
    args = prsr.parse_args()
    with open(args.config, 'r') as cf:
        cfg = json.load(cf)
    process(**cfg)

mp2022_dust.json

{
    "start": "2022-01-01T00", "end": "2022-12-31T23", "freq": "1h",
    "luftmpl": "inputs/mcip/LUFRAC_CRO.12US1.35L.%y%m%d",
    "g2dtmpl": "inputs/mcip/GRIDCRO2D.12US1.35L.%y%m%d",
    "m2dtmpl": "inputs/mcip/METCRO2D.12US1.35L.%y%m%d",
    "m3dtmpl": "inputs/mcip/METCRO3D.12US1.35L.%y%m%d",
    "vegtmpl": "inputs/mcip/VEG.12US1.35L.%y%m%d",
    "outtmpl": "outputs/dust_emis.12US1.%Y-%m.nc",
    "n_jobs": 3
}

plot_dust.py

try:
    import pycno
    addmaps = True
except Exception:
    addmaps = False

import xarray as xr
import matplotlib.pyplot as plt
import glob

projtmpl = (
    '+proj=lcc +lat_1={P_ALP} +lat_2={P_BET} +lat_0={YCENT} +lon_0={P_GAM}'
    + ' +x_0={x_0} +y_0={y_0} +R=6370000. +to_meter={XCELL} +no_defs'
)

outpaths = sorted(glob.glob('outputs/dust_emis.12US1.2022-??.nc'))

for outpath in outpaths:
    figpath = outpath.replace('.nc', '_DUST_EMIS.png')
    print(outpath)
    print(figpath)
    dustf = xr.open_dataset(outpath)
    attrs = {k: v for k, v in dustf.attrs.items()}
    attrs['x_0'] = -attrs['XORIG']
    attrs['y_0'] = -attrs['YORIG']
    attrs['y_0'] = -attrs['YORIG']
    dust_em = dustf['DUST_EMIS']
    fig, ax = plt.subplots()
    vmin, vmax = dust_em.where(lambda x: x > 0).quantile([.05, .95])
    Z = dust_em.mean(('TSTEP', 'LAY'))
    qm = ax.pcolormesh(Z, norm=plt.matplotlib.colors.LogNorm(vmin, vmax))
    ax.set(title=outpath)
    fig.colorbar(
        qm,
        label='Mean dust emissions [{units}]'.format(**dust_em.attrs)
    )
    if addmaps:
        cno = pycno.cno(projtmpl.format(**attrs))
        cno.drawstates(ax=ax)
    fig.savefig(figpath)

speciate_dust.json

{
    "ASO4J": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.0225", "units": "g/s"},
    "ANO3J": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.0002", "units": "g/s"},
    "ACLJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.00945", "units": "g/s"},
    "ANH4J": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 5e-05", "units": "g/s"},
    "ANAJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.03935", "units": "g/s"},
    "ACAJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.0794", "units": "g/s"},
    "AMGJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.019", "units": "g/s"},
    "AKJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.0377", "units": "g/s"},
    "APOCJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.01075", "units": "g/s"},
    "APNCOMJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.0043", "units": "g/s"},
    "ALVPO1J": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.01501", "units": "g/s"},
    "ALVOO1J": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 2.23e-05", "units": "g/s"},
    "AFEJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.03355", "units": "g/s"},
    "AALJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.05695", "units": "g/s"},
    "ASIJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.19425", "units": "g/s"},
    "ATIJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.0028", "units": "g/s"},
    "AMNJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.00115", "units": "g/s"},
    "AH2OJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.00541", "units": "g/s"},
    "AOTHRJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.48319", "units": "g/s"},
    "AMN_HAPSJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 0.00041", "units": "g/s"},
    "ANIJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 2e-05", "units": "g/s"},
    "ACR_IIIJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 5.6e-05", "units": "g/s"},
    "AASJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 5.5e-06", "units": "g/s"},
    "APBJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 1.6e-05", "units": "g/s"},
    "ACDJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 9.7e-08", "units": "g/s"},
    "APHGJ": {"expression": "DUST_EMIS * AREA * 1000 * 0.07 * 1.4e-08", "units": "g/s"},
    "ASO4K": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.0225", "units": "g/s"},
    "ANO3K": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.0002", "units": "g/s"},
    "ACLK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.00945", "units": "g/s"},
    "ANH4K": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 5e-05", "units": "g/s"},
    "ANAK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.03935", "units": "g/s"},
    "ACAK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.0794", "units": "g/s"},
    "AMGK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.019", "units": "g/s"},
    "AKK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.0377", "units": "g/s"},
    "APOCK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.01075", "units": "g/s"},
    "APNCOMK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.0043", "units": "g/s"},
    "ALVPO1K": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.01501", "units": "g/s"},
    "ALVOO1K": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 2.23e-05", "units": "g/s"},
    "AFEK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.03355", "units": "g/s"},
    "AALK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.05695", "units": "g/s"},
    "ASIK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.19425", "units": "g/s"},
    "ATIK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.0028", "units": "g/s"},
    "AMNK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.00115", "units": "g/s"},
    "AH2OK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.00541", "units": "g/s"},
    "AOTHRK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.48319", "units": "g/s"},
    "AMN_HAPSK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 0.00041", "units": "g/s"},
    "ANIK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 2e-05", "units": "g/s"},
    "ACR_IIIK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 5.6e-05", "units": "g/s"},
    "AASK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 5.5e-06", "units": "g/s"},
    "APBK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 1.6e-05", "units": "g/s"},
    "ACDK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 9.7e-08", "units": "g/s"},
    "APHGK": {"expression": "DUST_EMIS * AREA * 1000 * 0.93 * 1.4e-08", "units": "g/s"}
}

Needs update USEPA/CMAQ#222

Updated to match 1181. Includes improve shrubland representation and nlcd_modis combo option.

Updated to be json configuration file driven and separate plotting from operation.

MCIP files for a 2022 12US1 (not 12US2) are available on AWS.

https://epa-2022-modeling-platform.s3.amazonaws.com/index.html

Speciation is easily applied:

PMFINE = DUST_EMIS * 0.07
PMCOARSE = DUST_EMIS * 0.93
ASO4J = PMFINE * 0.0225
ANO3J = PMFINE * 0.0002
ACLJ = PMFINE * 0.00945
ANH4J = PMFINE * 5e-05
ANAJ = PMFINE * 0.03935
ACAJ = PMFINE * 0.0794
AMGJ = PMFINE * 0.019
AKJ = PMFINE * 0.0377
APOCJ = PMFINE * 0.01075
APNCOMJ = PMFINE * 0.0043
ALVPO1J = PMFINE * 0.01501
ALVOO1J = PMFINE * 2.23e-05
AFEJ = PMFINE * 0.03355
AALJ = PMFINE * 0.05695
ASIJ = PMFINE * 0.19425
ATIJ = PMFINE * 0.0028
AMNJ = PMFINE * 0.00115
AH2OJ = PMFINE * 0.00541
AOTHRJ = PMFINE * 0.48319
AMN_HAPSJ = PMFINE * 0.00041
ANIJ = PMFINE * 2e-05
ACR_IIIJ = PMFINE * 5.6e-05
AASJ = PMFINE * 5.5e-06
APBJ = PMFINE * 1.6e-05
ACDJ = PMFINE * 9.7e-08
APHGJ = PMFINE * 1.4e-08
ASO4K = PMCOARSE * 0.0225
ANO3K = PMCOARSE * 0.0002
ACLK = PMCOARSE * 0.00945
ANH4K = PMCOARSE * 5e-05
ANAK = PMCOARSE * 0.03935
ACAK = PMCOARSE * 0.0794
AMGK = PMCOARSE * 0.019
AKK = PMCOARSE * 0.0377
APOCK = PMCOARSE * 0.01075
APNCOMK = PMCOARSE * 0.0043
ALVPO1K = PMCOARSE * 0.01501
ALVOO1K = PMCOARSE * 2.23e-05
AFEK = PMCOARSE * 0.03355
AALK = PMCOARSE * 0.05695
ASIK = PMCOARSE * 0.19425
ATIK = PMCOARSE * 0.0028
AMNK = PMCOARSE * 0.00115
AH2OK = PMCOARSE * 0.00541
AOTHRK = PMCOARSE * 0.48319
AMN_HAPSK = PMCOARSE * 0.00041
ANIK = PMCOARSE * 2e-05
ACR_IIIK = PMCOARSE * 5.6e-05
AASK = PMCOARSE * 5.5e-06
APBK = PMCOARSE * 1.6e-05
ACDK = PMCOARSE * 9.7e-08
APHGK = PMCOARSE * 1.4e-08

p.s., these are based on CMAQ’s code lines below.

https://github.com/USEPA/CMAQ/blob/main/CCTM/src/emis/emis/DUST_EMIS.F#L82C1-L86C1
C updated values of mass fraction for "freshly emitted dust"
C based on Kok [PNAS, 2011] and Nabat et al. [ACP, 2012]
real, parameter :: fracmj = 0.07 ! mass fraction assigned to accum mode
real, parameter :: fracmk = 0.93 ! mass fraction assigned to coarse mode

https://github.com/USEPA/CMAQ/blob/main/CCTM/src/aero/aero6/AERO_DATA.F#L502-L532
Type( emis_table ) :: dust_spc( ndust_spc ) = (/
C --description-- -- Species -- --MW-- ------ spcfac ------
C fine coarse
& emis_table( 'Sulfate ', 'SO4 ', 96.0, (/ 0.02250, 0.02655/) ), ! Sulfate
& emis_table( 'Nitrate ', 'NO3 ', 62.0, (/ 0.00020, 0.00160/) ), ! Nitrate
& emis_table( 'Chlorine ', 'CL ', 35.5, (/ 0.00945, 0.01190/) ), ! Chlorine
& emis_table( 'Ammonium ', 'NH4 ', 18.0, (/ 0.00005, 0.0 /) ), ! Ammonium
& emis_table( 'Sodium ', 'NA ', 23.0, (/ 0.03935, 0.0 /) ), ! Sodium
& emis_table( 'Calcium ', 'CA ', 40.1, (/ 0.07940, 0.0 /) ), ! Calcium
& emis_table( 'Magnesium ', 'MG ', 24.3, (/ 0.01900, 0.0 /) ), ! Magnesium
& emis_table( 'Potassium ', 'K ', 39.1, (/ 0.03770, 0.0 /) ), ! Potassium
& emis_table( 'Org. Carbon ', 'POC ', 220.0, (/ 0.01075, 0.0 /) ), ! Organic Carbon
& emis_table( 'NonCarbon Org.', 'PNCOM ', 220.0, (/ 0.00430, 0.0 /) ), ! Non-Carbon Organic Matter
& emis_table( 'Low-Vol. POA ', 'LVPO1 ', 218.0, (/ 0.01501, 0.0 /) ), ! Non-Carbon Organic Matter
& emis_table( 'Low-Vol. OOA ', 'LVOO1 ', 136.0, (/ 2.23E-5, 0.0 /) ), ! Non-Carbon Organic Matter
& emis_table( 'Black Carbon ', 'EC ', 12.0, (/ 0.0, 0.0 /) ), ! Black or Elemental Carbon
& emis_table( 'Iron ', 'FE ', 55.8, (/ 0.03355, 0.0 /) ), ! Iron
& emis_table( 'Aluminum ', 'AL ', 27.0, (/ 0.05695, 0.0 /) ), ! Aluminum
& emis_table( 'Silicon ', 'SI ', 28.1, (/ 0.19425, 0.0 /) ), ! Silicon
& emis_table( 'Titanium ', 'TI ', 47.9, (/ 0.00280, 0.0 /) ), ! Titanium
& emis_table( 'Manganese ', 'MN ', 54.9, (/ 0.00115, 0.0 /) ), ! Manganese
& emis_table( 'Water ', 'H2O ', 18.0, (/ 0.00541, 0.00637/) ), ! Water
& emis_table( 'Undefined Mass', 'OTHR ', 200.0, (/ 0.48319, 0.0 /) ), ! Other
& emis_table( 'Non-Anion Dust', 'SOIL ', 100.0, (/ 0.0, 0.95358/) ), ! Non-Anion Dust
& emis_table( 'Air Toxics Mn ', 'MN_HAPS', 54.9, (/ 0.00041, 0.00041/) ), ! Air toxics Manganese J and K mode
& emis_table( 'Nickel ', 'NI ', 58.7, (/ 2.0E-05, 2.0E-05/) ), ! Nickel J and K mode
& emis_table( 'Chromium III ', 'CR_III ', 52.0, (/ 5.6E-05, 5.6E-05/) ), ! Trivalent Chromium J and K mode
& emis_table( 'Arsenic ', 'AS ', 74.92,(/ 5.5E-06, 5.5E-06/) ), ! Arsenic J and K mode
& emis_table( 'Lead ', 'PB ', 207.2, (/ 1.6E-05, 1.6E-05/) ), ! Lead J and K mode
& emis_table( 'Cadmium ', 'CD ', 112.4, (/ 9.7E-08, 9.7E-08/) ), ! Cadmium J and K mode
& emis_table( 'Mercury ', 'PHG ', 200.5, (/ 1.4E-08, 1.4E-08/) ) /) ! Mercury J and K mode