diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 25c9109f79..fad04411e8 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -14,6 +14,7 @@ import scipy.integrate as si import scipy.optimize as so +from scipy.special import lambertw import xarray as xr from .exceptions import InvalidSoundingError @@ -272,7 +273,9 @@ def water_latent_heat_melting(temperature): @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'dewpoint')) -@check_units('[temperature]', '[temperature]') +@process_units(input_dimensionalities={'temperature': '[temperature]', + 'dewpoint': '[temperature]'}, + output_dimensionalities='[dimensionless]') def relative_humidity_from_dewpoint(temperature, dewpoint): r"""Calculate the relative humidity. @@ -311,8 +314,8 @@ def relative_humidity_from_dewpoint(temperature, dewpoint): Renamed ``dewpt`` parameter to ``dewpoint`` """ - e = saturation_vapor_pressure(dewpoint) - e_s = saturation_vapor_pressure(temperature) + e = saturation_vapor_pressure._nounit(dewpoint) + e_s = saturation_vapor_pressure._nounit(temperature) return e / e_s @@ -647,7 +650,7 @@ def dt(p, t): {'pressure': '[pressure]', 'temperature': '[temperature]', 'dewpoint': '[temperature]'}, ('[pressure]', '[temperature]') ) -def lcl(pressure, temperature, dewpoint, max_iters=50, eps=1e-5): +def lcl(pressure, temperature, dewpoint, max_iters=None, eps=None): r"""Calculate the lifted condensation level (LCL) from the starting point. The starting state for the parcel is defined by `temperature`, `dewpoint`, @@ -673,20 +676,12 @@ def lcl(pressure, temperature, dewpoint, max_iters=50, eps=1e-5): `pint.Quantity` LCL temperature - Other Parameters - ---------------- - max_iters : int, optional - The maximum number of iterations to use in calculation, defaults to 50. - - eps : float, optional - The desired relative error in the calculated value, defaults to 1e-5. - Examples -------- >>> from metpy.calc import lcl >>> from metpy.units import units >>> lcl(943 * units.hPa, 33 * units.degC, 28 * units.degC) - (, ) + (, ) See Also -------- @@ -694,39 +689,41 @@ def lcl(pressure, temperature, dewpoint, max_iters=50, eps=1e-5): Notes ----- - This function is implemented using an iterative approach to solve for the - LCL. The basic algorithm is: - - 1. Find the dewpoint from the LCL pressure and starting mixing ratio - 2. Find the LCL pressure from the starting temperature and dewpoint - 3. Iterate until convergence - - The function is guaranteed to finish by virtue of the `max_iters` counter. + From [Romps2017]_, this directly solves for the temperature at the LCL, Eq 22a, + .. math:: T_{LCL} = c [W_{-1}(RH_l^{1/a} c \exp{c})]^{-1} T + and similarly for pressure at the LCL, Eq 22b, + .. math:: p_{LCL} = p \left( \frac{T_{LCL}}{T} \right)^{c_{pm} / R_m} + where :math:`a` (Eq 22d), :math:`b` (Eq 22e), and :math:`c` (Eq 22f) are derived constants, + and :math:`W_{-1}` is the :math:`k=-1` branch of the Lambert :math:`W` function. .. versionchanged:: 1.0 Renamed ``dewpt`` parameter to ``dewpoint`` """ - def _lcl_iter(p, p0, w, t): - nonlocal nan_mask - td = globals()['dewpoint']._nounit(vapor_pressure._nounit(p, w)) - p_new = (p0 * (td / t) ** (1. / mpconsts.nounit.kappa)) - nan_mask = nan_mask | np.isnan(p_new) - return np.where(np.isnan(p_new), p, p_new) - - # Handle nans by creating a mask that gets set by our _lcl_iter function if it - # ever encounters a nan, at which point pressure is set to p, stopping iteration. - nan_mask = False - w = mixing_ratio._nounit(saturation_vapor_pressure._nounit(dewpoint), pressure) - lcl_p = so.fixed_point(_lcl_iter, pressure, args=(pressure, w, temperature), - xtol=eps, maxiter=max_iters) - lcl_p = np.where(nan_mask, np.nan, lcl_p) - - # np.isclose needed if surface is LCL due to precision error with np.log in dewpoint. - # Causes issues with parcel_profile_with_lcl if removed. Issue #1187 - lcl_p = np.where(np.isclose(lcl_p, pressure), pressure, lcl_p) - - return lcl_p, globals()['dewpoint']._nounit(vapor_pressure._nounit(lcl_p, w)) + if max_iters or eps: + _warnings.warn( + 'max_iters, eps arguments unused and will be deprecated in a future version.', + PendingDeprecationWarning) + + q = specific_humidity_from_dewpoint._nounit(pressure, dewpoint) + moist_heat_ratio = (moist_air_specific_heat_pressure._nounit(q) + / moist_air_gas_constant._nounit(q)) + spec_heat_diff = mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v + + a = moist_heat_ratio + spec_heat_diff / mpconsts.nounit.Rv + b = (-(mpconsts.nounit.Lv + spec_heat_diff * mpconsts.nounit.T0) + / (mpconsts.nounit.Rv * temperature)) + c = b / a + + w_minus1 = lambertw( + (relative_humidity_from_dewpoint._nounit(temperature, dewpoint)**(1 / a) + * c * np.exp(c)), + k=-1).real + + t_lcl = c / w_minus1 * temperature + p_lcl = pressure * (t_lcl / temperature) ** moist_heat_ratio + + return p_lcl, t_lcl @exporter.export