Skip to content

Instantly share code, notes, and snippets.

@barronh

barronh/cmaq_dust.py

Last active April 23, 2024 16:07
Show Gist options
  • Save barronh/a34a50c8d8a5982361cc3bbc52d06e30 to your computer and use it in GitHub Desktop.
Save barronh/a34a50c8d8a5982361cc3bbc52d06e30 to your computer and use it in GitHub Desktop.
Download ZIP
Offline CMAQ Dust Emulator
Raw
cmaq_dust.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.
To-Do:
1. Currently only operational for NLCD40 and MODIS_NOAH
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
outpath = process('2022-04-01', '2022-05-01')
### Run for a specific date/time
dust_em = get_dust_emis(
2022091, 0, 1, 20., metroot='/work/ROMO/met/12US_2022/mcip'
)
"""
import warnings
import pandas as pd
import xarray as xr
import numpy as np
import os
warnings.simplefilter('ignore')
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)
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_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
dmap_nlcd40 = np.array([ # land use type desert map to BELD3
0, # shrubland
0, # shrubland
2, # barrenland
2, # barrenland
2, # barrenland
2, # barrenland
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] == 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(sdate, edate, metroot, datefmt='.12US1.35L.%y%m%d', outdir='.'):
"""
# Required
lufpath = date.strftime('{metroot}/LUFRAC_CRO{datefmt}')
g2dpath = date.strftime('{metroot}/GRIDCRO2D{datefmt}')
m2dpath = date.strftime('{metroot}/METCRO2D{datefmt}')
m3dpath = date.strftime('{metroot}/METCRO3D{datefmt}')
# Optional, if not available, it will be skipped
vegpath = date.strftime('{metroot}/VEG{appl}{datefmt}')
"""
dust_ems = []
dates = pd.date_range(sdate, edate, freq='h')
bdate = dates.min()
edate = dates.max()
outpath = f'{outdir}/dust_emis_{bdate:%FT%H%M%S}_{edate:%FT%H%M%S}.nc'
if os.path.exists(outpath):
print(f'WARNING: keeping cached {outpath}')
return outpath
for date in dates:
print(f'\r{date:%F}', end='')
jdate = int(date.strftime('%Y%j'))
jtime = int(date.strftime('%H%M%S'))
lufpath = date.strftime(f'{metroot}/LUFRAC_CRO{datefmt}')
g2dpath = date.strftime(f'{metroot}/GRIDCRO2D{datefmt}')
m2dpath = date.strftime(f'{metroot}/METCRO2D{datefmt}')
m3dpath = date.strftime(f'{metroot}/METCRO3D{datefmt}')
vegpath = date.strftime(f'{metroot}/VEG{datefmt}')
if os.path.exists(vegpath):
print('.v', end='')
else:
vegpath = None
dust_em = get_dust_emis(
jdate, jtime, 1, 20., # assume 1 / 20 as conversion from g/cm2/
lufpath=lufpath, g2dpath=g2dpath, m2dpath=m2dpath,
m3dpath=m3dpath, vegpath=vegpath
)
dust_ems.append(dust_em)
print()
dust_ems = np.stack(dust_ems)
dust_em.mean(0)
nr, nc = dust_em.shape
dustvar = xr.DataArray(
dust_ems[:, None, :, :] / 1e3, name='dust_emis',
dims=('TSTEP', 'LAY', 'ROW', 'COL'),
attrs=dict(
long_name='dust_emis', units='kg/m**3/s',
description='Assuming 20m layer depth'
),
coords={
'TSTEP': dates, 'LAY': np.array([0.9985]),
'ROW': np.arange(nr) + 0.5, 'COL': np.arange(nc) + 0.5
}
)
g2d = xr.open_dataset(f'{metroot}/GRIDCRO2D.12US1.35L.220101')
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', 'long_name': 'area'}
)
dustf.to_netcdf(outpath)
return outpath
if __name__ == '__main__':
import matplotlib.pyplot as plt
metroot = input('Enter the path to your met: ')
# metroot = '/home/bhenders/temp/dust/modis_noah/mcip'
print('metroot', metroot)
outpath = process('2022-01-01', '2022-01-31T23', metroot=metroot)
outpath = process('2022-02-01', '2022-02-28T23', metroot=metroot)
outpath = process('2022-03-01', '2022-03-31T23', metroot=metroot)
outpath = process('2022-04-01', '2022-04-30T23', metroot=metroot)
outpath = process('2022-05-01', '2022-05-31T23', metroot=metroot)
outpath = process('2022-06-01', '2022-06-30T23', metroot=metroot)
outpath = process('2022-07-01', '2022-07-31T23', metroot=metroot)
outpath = process('2022-08-01', '2022-08-31T23', metroot=metroot)
outpath = process('2022-09-01', '2022-09-30T23', metroot=metroot)
outpath = process('2022-10-01', '2022-10-31T23', metroot=metroot)
outpath = process('2022-11-01', '2022-11-30T23', metroot=metroot)
outpath = process('2022-12-01', '2022-12-31T23', metroot=metroot)
dustf = xr.open_dataset(outpath)
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='dust emissions [{units}]'.format(**dust_em.attrs))
fig.savefig('test.png')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment