Skip to content

ENH: air density correction factor to V_D #142

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
isilber merged 2 commits into from
May 16, 2025
Merged

ENH: air density correction factor to V_D #142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d6f91e9
ENH: air density correction factor to V_D
isilber May 16, 2025
1be780e
ENH: add `Model` class method to calculate air density during init
isilber May 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 57 additions & 6 deletions emc2/core/model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,11 @@
This module contains the Model class and example Models for your use.

"""
import warnings

import xarray as xr
import numpy as np
from netCDF4 import Dataset

from scipy.special import gamma
from scipy.interpolate import RegularGridInterpolator
Expand All @@ -18,7 +21,6 @@
print('Using act-atmos v1.5.3 or earlier. Please update to v2.0.0 or newer')
from act.io.armfiles import read_netcdf
from .instrument import ureg, quantity
from netCDF4 import Dataset
from ..scattering import brandes

try:
Expand Down Expand Up @@ -85,6 +87,13 @@ class Model():
hydrometeor class
asp_ratio_func: dict
A dictionary that returns hydrometeor aspect ratios as a function of maximum dimension in mm.
Vd_rhoa_scaling_ref: dict
dict keys are hydrometeor classes with dict values being the reference air density [kg m-3] used
for the scaling of reflectivity-weighted fall velocity (Mean Doppler velocity minus air motion)
by the ratio (rhoa_0 / rhoa) ** 0.54, where rho_a_0 is the reference air density and rhoa is
the air density. This correction, implemented in many microphysics schemes
(e.g., MG2, P3), followes Heymsfield et al. (JAS, 2007).
Set key values to None (or exclude hydrometeor class from keys) to skip this correction.
hyd_types: list
list of hydrometeor classes to include in calcuations. By default set to be consistent
with the model represented by the Model subclass.
Expand Down Expand Up @@ -160,6 +169,7 @@ def __init__(self):
"Avogadro_c": 6.022140857e23,
"R": 8.3144598} # J K^-1 mol^-1
self.asp_ratio_func = {}
self.Vd_rhoa_scaling_ref = {}
self.ice_hyd_types = ["ci", "pi"]

def _add_vel_units(self):
Expand All @@ -176,6 +186,19 @@ def _prepare_variables(self):
continue
self.ds[variable].attrs = attrs

def _calc_air_density(self):
"""
Calculates the air density and stores it in the dataset under the key "rho_a".
The calculation is based on the ideal gas law using temperature [K] and pressure [hPa] fields.

"""
if np.all([self.T_field in self.ds.keys(), self.T_field in self.ds.keys()]):
self.ds["rho_a"] = self.ds[self.p_field] * 1e2 / (self.consts["R_d"] * self.ds[self.T_field])
self.ds["rho_a"].attrs["units"] = "kg / m ** 3"
else:
warnings.warn(f"Temperature and pressure field not provided in `Model` init; "
f"this could introduce errors if the microphysics approach is applied", RuntimeWarning)

def _crop_bounding_box(self, bounding_box):
"""
Crop the input region to a given bounding box for a regional model.
Expand Down Expand Up @@ -638,6 +661,10 @@ def __init__(self, file_path, time_range=None, load_processed=False):
self.strat_re_fields = {'cl': 're_sscl', 'ci': 're_ssci', 'pi': 're_sspi', 'pl': 're_sspl'}
self.q_names_convective = {'cl': 'QCLmc', 'ci': 'QCImc', 'pl': 'QPLmc', 'pi': 'QPImc'}
self.q_names_stratiform = {'cl': 'qcl', 'ci': 'qci', 'pl': 'qpl', 'pi': 'qpi'}
self.Vd_rhoa_scaling_ref = {"cl": 1.0, # Ikawa and Saito (1990, p. 85)
"ci": 1.0, # Ikawa and Saito (1990, p. 85)
"pl": 1.1, # Liu and Orville (1969, p. 1269),
"pi": 1.15} # Localleti and Hobbs (1974) - est. from 750-1500 alt + T < 273.15
self.hyd_types = ["cl", "ci", "pl", "pi"]
self.mcphys_scheme = "MG2" # per GM2015 (J. Clim.)
self.rad_scheme_family = "ModelE"
Expand All @@ -664,6 +691,7 @@ def __init__(self, file_path, time_range=None, load_processed=False):

# ModelE has pressure units in mb, but pint only supports hPa
self.ds["p_3d"].attrs["units"] = "hPa"
self._calc_air_density() # calculate air density with input T and p fields

self.model_name = "ModelE3"

Expand Down Expand Up @@ -741,6 +769,11 @@ def __init__(self, file_path, time_range=None, load_processed=False, time_dim="t
self.mu_field = {'cl': 'mu_cloud', 'ci': None, 'pl': None, 'pi': None}
self.lambda_field = {'cl': 'lambda_cloud', 'ci': None, 'pl': None, 'pi': None}

self.Vd_rhoa_scaling_ref = {"cl": 1.0, # Ikawa and Saito (1990, p. 85)
"ci": 1.0, # Ikawa and Saito (1990, p. 85)
"pl": 1.1, # Liu and Orville (1969, p. 1269),
"pi": 1.15} # Localleti and Hobbs (1974) - est. from 750-1500 alt + T < 273.15

if include_rain_in_rt:
self.hyd_types = ["cl", "ci", "pl"]
else:
Expand Down Expand Up @@ -813,8 +846,7 @@ def __init__(self, file_path, time_range=None, load_processed=False, time_dim="t
self.ds["zeros_cf"] = xr.DataArray(np.zeros_like(self.ds[self.p_field].values),
dims=self.ds[self.p_field].dims)
self.ds["zeros_cf"].attrs["long_name"] = "An array of zeros as only strat output is used for this model"
self.ds["rho_a"] = self.ds[self.p_field] * 1e2 / (self.consts["R_d"] * self.ds[self.T_field])
self.ds["rho_a"].attrs["units"] = "kg / m ** 3"
self._calc_air_density() # calculate air density with input T and p fields
precip_classes = [x for x in self.hyd_types if x[0] == "p"]
for hyd in self.hyd_types:

Expand Down Expand Up @@ -917,6 +949,10 @@ def __init__(self, file_path, time_range=None, load_processed=False, time_dim="t
if mcphys_scheme == "P3":
self.rad_scheme_family = "P3" # E3SMv3 uses CESM bulk LUTs but with a P3 twist (crp and Fr dependence)

self.Vd_rhoa_scaling_ref = {"cl": None, # N/A (see Morrison and Milbrandt (2015, p. 294)
"ci": 0.8254, # p=600 hPa, T=253.15 K (used in P3 LUT generator code)
"pl": 1.1} # Liu and Orville (1969, p. 1269)

# if Instrument object was input, we use the automatically-loaded scattering LUT to process PSD params
# (instrument type (HSRL, etc.) has no effect on PSD params).
if (instrument is not None):
Expand Down Expand Up @@ -1136,12 +1172,17 @@ def __init__(self, file_path, z_range=None, time_range=None,
'pl': 're_plc', 'pi': 're_pic'}
self.strat_re_fields = {'cl': 're_cls', 'ci': 're_cis',
'pl': 're_pls', 'pi': 're_pis'}
self.Vd_rhoa_scaling_ref = {"cl": 1.0, # Ikawa and Saito (1990, p. 85)
"ci": 1.0, # Ikawa and Saito (1990, p. 85)
"pl": 1.1, # Liu and Orville (1969, p. 1269),
"pi": 1.15} # Localleti and Hobbs (1974) - est. based on paper info
if include_gr_class:
self.hyd_types += ['gr']
if gr_class_is_hail:
self.Rho_hyd.update({'gr': 900.0 * ureg.kg / (ureg.m**3)})
self.vel_param_a.update({'gr': 114.5}) # MATSUN AND HUGGINS 1980
self.vel_param_b.update({'gr': 0.5 * ureg.dimensionless})
self.Vd_rhoa_scaling_ref["gr"] = 1.23 # Matsun and Huggins (1980, p. 1116)
else:
self.Rho_hyd.update({'gr': 400.0 * ureg.kg / (ureg.m**3)})
self.vel_param_a.update({'gr': 19.3})
Expand Down Expand Up @@ -1243,9 +1284,9 @@ def __init__(self, file_path, z_range=None, time_range=None,
self.ds["pressure"].attrs["units"] = "hPa"
self.ds["T"].attrs["units"] = "K"
self.ds["Z"].attrs["units"] = "m"
rho = (self.ds["pressure"] * 1e2) / (287.058 * self.ds["T"])
self._calc_air_density() # calculate air density with input T and p fields
# Qn in kg-1 --> cm-3 * rho to get m-3 * 1e-6 for cm-3
qn_conversion = rho.values * 1e-6
qn_conversion = self.ds["rho_a"].values * 1e-6
W = getvar(wrfin, "wa", units="m s-1", timeidx=ALL_TIMES, squeeze=False)
if NUWRF:
cldfrac = ds["CLDFRA"].values
Expand Down Expand Up @@ -1403,6 +1444,10 @@ def __init__(self, file_path, time_range=None, time_dim="dom_col", single_pi_cla
'pi': 're_strat_pi', 'pl': 're_strat_pl'}
self.q_names_convective = {'cl': 'conv_dat', 'ci': 'conv_dat', 'pl': 'conv_dat', 'pi': 'conv_dat'}
self.q_names_stratiform = {'cl': 'qcl', 'ci': 'qci', 'pl': 'qpl', 'pi': 'qpi'}
self.Vd_rhoa_scaling_ref = {"cl": 1.0, # Ikawa and Saito (1990, p. 85)
"ci": 1.0, # Ikawa and Saito (1990, p. 85)
"pl": 1.1, # Liu and Orville (1969, p. 1269),
"pi": 1.15} # Localleti and Hobbs (1974) - est. from 750-1500 alt + T < 273.15
self.hyd_types = ["cl", "ci", "pl", "pi"]
if not single_pi_class:
self.N_field.update({'pir': 'npir'})
Expand All @@ -1414,6 +1459,7 @@ def __init__(self, file_path, time_range=None, time_dim="dom_col", single_pi_cla
self.strat_re_fields.update({'pir': 'strat_pir_frac'})
self.q_names_convective.update({'pir': 'conv_dat'})
self.q_names_stratiform.update({'pir': 'qpir'})
self.Vd_rhoa_scaling_ref['pir'] = 1.15 # Localleti and Hobbs (1974) - est. based on paper info
self.hyd_types.append("pir")
self.mcphys_scheme="MG" # MG2008
self.rad_scheme_family = "ModelE" # Similar implementation to ModelE3
Expand All @@ -1427,7 +1473,8 @@ def __init__(self, file_path, time_range=None, time_dim="dom_col", single_pi_cla
my_attrs = self.ds[variable].attrs
self.ds[variable] = self.ds[variable].astype('float64')
self.ds[variable].attrs = my_attrs

self._calc_air_density() # calculate air density with input T and p fields

if bounding_box is not None:
super()._crop_bounding_box_x_y(bounding_box)

Expand Down Expand Up @@ -1531,6 +1578,7 @@ def __init__(self):
self.ds = my_ds
self.height_dim = "height"
self.time_dim = "time"
self._calc_air_density() # calculate air density with input T and p fields
self.hyd_types = ["cl", "ci", "pl", "pi"]


Expand Down Expand Up @@ -1670,6 +1718,7 @@ def __init__(self):
self.mcphys_scheme = "MG2" # required to prevent errors from being raised.
self.rad_scheme_family = "ModelE" # required to prevent errors from being raised.
self.ds = my_ds
self._calc_air_density() # calculate air density with input T and p fields
self.height_dim = "height"
self.time_dim = "time"

Expand Down Expand Up @@ -1812,6 +1861,7 @@ def __init__(self):
self.mcphys_scheme = "MG2" # required to prevent errors from being raised.
self.rad_scheme_family = "ModelE" # required to prevent errors from being raised.
self.ds = my_ds
self._calc_air_density() # calculate air density with input T and p fields
self.height_dim = "height"
self.time_dim = "time"

Expand Down Expand Up @@ -1934,5 +1984,6 @@ def __init__(self):
self.mcphys_scheme = "MG2" # required to prevent errors from being raised.
self.rad_scheme_family = "ModelE" # required to prevent errors from being raised.
self.ds = my_ds
self._calc_air_density() # calculate air density with input T and p fields
self.height_dim = "height"
self.time_dim = "time"
Loading