diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index e5e7a4f787f..a9399b364b6 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -39,6 +39,8 @@ Moist Thermodynamics mixing_ratio mixing_ratio_from_relative_humidity mixing_ratio_from_specific_humidity + moist_air_gas_constant + moist_air_specific_heat_pressure moist_lapse moist_static_energy precipitable_water @@ -60,6 +62,9 @@ Moist Thermodynamics virtual_potential_temperature virtual_temperature virtual_temperature_from_dewpoint + water_latent_heat_melting + water_latent_heat_sublimation + water_latent_heat_vaporization wet_bulb_temperature wet_bulb_potential_temperature diff --git a/docs/api/references.rst b/docs/api/references.rst index dacd4189546..87af83af82f 100644 --- a/docs/api/references.rst +++ b/docs/api/references.rst @@ -2,6 +2,10 @@ References ========== +.. [Ambaum2020] Ambaum MHP, 2020: Accurate, simple equation for saturated vapour pressure over water and ice. + *QJR Meteorol Soc.*; **146**: 4252–4258, + doi: `10.1002/qj.3899 `_. + .. [Anderson2013] Anderson, G. B., M. L. Bell, and R. D. Peng, 2013: Methods to Calculate the Heat Index as an Exposure Metric in Environmental Health Research. *Environmental Health Perspectives*, **121**, 1111-1119, @@ -176,6 +180,10 @@ References .. [Rochette2006] Rochette, Scott M., and Patrick S. Market. "A primer on the ageostrophic wind." Natl. Weather Dig. 30 (2006): 17-28. +.. [Romps2017] Romps, D. M., 2017: Exact Expression for the Lifting Condensation Level. + *J. Atmos. Sci.*, **74**, 3891–3900, + doi: `10.1175/JAS-D-17-0102.1. `_. + .. [Rothfusz1990] Rothfusz, L.P.: *The Heat Index "Equation"*. Fort Worth, TX: Scientific Services Division, NWS Southern Region Headquarters, 1990. `SR90-23 <../_static/rothfusz-1990-heat-index-equation.pdf>`_, 2 pp. diff --git a/docs/conf.py b/docs/conf.py index f0d500c8f7f..baaf08a7645 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -455,6 +455,7 @@ # Couldn't fix these 403's with user agents r'https://doi\.org/10\.1029/2010GL045777', r'https://doi\.org/10\.1098/rspa\.2004\.1430', + r'https://doi\.org/10\.1002/qj\.3899', # Currently giving certificate errors on GitHub r'https://library.wmo.int/.*', # For some reason GHA gets a 403 from Stack Overflow diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 983e1315474..fad04411e81 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 @@ -29,9 +30,252 @@ exporter = Exporter(globals()) +@exporter.export +@preprocess_and_wrap(wrap_like='specific_humidity') +@process_units(input_dimensionalities={'specific_humidity': 'dimensionless'}, + output_dimensionalities='[specific_heat_capacity]', + output_to='J K**-1 kg**-1 ') +def moist_air_gas_constant(specific_humidity): + r"""Calculate R_m, the specific gas constant for a parcel of moist air. + + Parameters + ---------- + specific_humidity : `pint.Quantity` + + Returns + ------- + `pint.Quantity` + Specific gas constant + + Examples + -------- + >>> from metpy.calc import moist_air_gas_constant + >>> from metpy.units import units + >>> moist_air_gas_constant(11 * units('g/kg')) + + + See Also + -------- + moist_air_specific_heat_pressure, moist_air_poisson_exponent + + Notes + ----- + .. math:: R_m = (1 - q_v) R_a + q_v R_v + + Eq 16, [Romps2017]_ using MetPy-defined constants in place of cited values. + + """ + return ((1 - specific_humidity) * mpconsts.nounit.Rd + + specific_humidity * mpconsts.nounit.Rv) + + +@exporter.export +@preprocess_and_wrap(wrap_like='specific_humidity') +@process_units(input_dimensionalities={'specific_humidity': 'dimensionless'}, + output_dimensionalities='[specific_heat_capacity]', + output_to='J K**-1 kg**-1 ') +def moist_air_specific_heat_pressure(specific_humidity): + r"""Calculate C_pm, the specific heat at constant pressure for a moist air parcel. + + Parameters + ---------- + specific_humidity : `pint.Quantity` + + Returns + ------- + `pint.Quantity` + Specific heat capacity of air at constant pressure + + Examples + -------- + >>> from metpy.calc import moist_air_specific_heat_pressure + >>> from metpy.units import units + >>> moist_air_specific_heat_pressure(11 * units('g/kg')) + + + See Also + -------- + moist_air_gas_constant, moist_air_poisson_exponent + + Notes + ----- + .. math:: c_{pm} = (1 - q_v) c_{pa} + q_v c_{pv} + + Eq 17, [Romps2017]_ using MetPy-defined constants in place of cited values. + + """ + return ((1 - specific_humidity) * mpconsts.nounit.Cp_d + + specific_humidity * mpconsts.nounit.Cp_v) + + +@exporter.export +@preprocess_and_wrap(wrap_like='specific_humidity') +@process_units( + input_dimensionalities={'specific_humidity': 'dimensionless'}, + output_dimensionalities='[dimensionless]' +) +def moist_air_poisson_exponent(specific_humidity): + r"""Calculate kappa_m, the Poisson exponent for a moist air parcel. + + Parameters + ---------- + specific_humidity : `pint.Quantity` + + Returns + ------- + `pint.Quantity` + Poisson exponent of moist air parcel + + Examples + -------- + >>> from metpy.calc import moist_air_poisson_exponent + >>> from metpy.units import units + >>> moist_air_poisson_exponent(11 * units('g/kg')) + + + See Also + -------- + moist_air_gas_constant, moist_air_specific_heat_pressure + + """ + return (moist_air_gas_constant._nounit(specific_humidity) + / moist_air_specific_heat_pressure._nounit(specific_humidity)) + + +@exporter.export +@preprocess_and_wrap(wrap_like='temperature') +@process_units(input_dimensionalities={'temperature': '[temperature]'}, + output_dimensionalities='[specific_enthalpy]', + output_to='J kg**-1') +def water_latent_heat_vaporization(temperature): + r"""Calculate the latent heat of vaporization for water. + + Accounts for variations in latent heat across valid temperature range. + + Parameters + ---------- + temperature : `pint.Quantity` + + Returns + ------- + `pint.Quantity` + Latent heat of vaporization + + Examples + -------- + >>> from metpy.calc import water_latent_heat_vaporization + >>> from metpy.units import units + >>> water_latent_heat_vaporization(20 * units.degC) + + + See Also + -------- + water_latent_heat_sublimation, water_latent_heat_melting + + Notes + ----- + Assumption of constant :math:`C_pv` limits validity to :math:`0 \deg \en 100 \deg C` range. + + .. math:: L = L_0 - (c_{pl} - c_{pv}) (T - T_0) + + Eq 15, [Ambaum2020]_, using MetPy-defined constants in place of cited values. + + """ + return (mpconsts.nounit.Lv + - (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v) + * (temperature - mpconsts.nounit.T0)) + + +@exporter.export +@preprocess_and_wrap(wrap_like='temperature') +@process_units(input_dimensionalities={'temperature': '[temperature]'}, + output_dimensionalities='[specific_enthalpy]', + output_to='J kg**-1') +def water_latent_heat_sublimation(temperature): + r"""Calculate the latent heat of sublimation for water. + + Accounts for variations in latent heat across valid temperature range. + + Parameters + ---------- + temperature : `pint.Quantity` + + Returns + ------- + `pint.Quantity` + Latent heat of vaporization + + Examples + -------- + >>> from metpy.calc import water_latent_heat_sublimation + >>> from metpy.units import units + >>> water_latent_heat_sublimation(-15 * units.degC) + + + See Also + -------- + water_latent_heat_vaporization, water_latent_heat_melting + + Notes + ----- + .. math:: L_s = L_{s0} - (c_{pl} - c_{pv}) (T - T_0) + + Eq 18, [Ambaum2020]_, using MetPy-defined constants in place of cited values. + + """ + return (mpconsts.nounit.Ls + - (mpconsts.nounit.Cp_i - mpconsts.nounit.Cp_v) + * (temperature - mpconsts.nounit.T0)) + + +@exporter.export +@preprocess_and_wrap(wrap_like='temperature') +@process_units(input_dimensionalities={'temperature': '[temperature]'}, + output_dimensionalities='[specific_enthalpy]', + output_to='J kg**-1') +def water_latent_heat_melting(temperature): + r"""Calculate the latent heat of melting for water. + + Accounts for variations in latent heat across valid temperature range. + + Parameters + ---------- + temperature : `pint.Quantity` + + Returns + ------- + `pint.Quantity` + Latent heat of vaporization + + Examples + -------- + >>> from metpy.calc import water_latent_heat_melting + >>> from metpy.units import units + >>> water_latent_heat_melting(-15 * units.degC) + + + See Also + -------- + water_latent_heat_vaporization, water_latent_heat_sublimation + + Notes + ----- + .. math:: L_m = L_{m0} + (c_{pl} - c_{pi}) (T - T_0) + + Body text below Eq 20, [Ambaum2020]_, derived from Eq 15, Eq 18. + Uses MetPy-defined constants in place of cited values. + + """ + return (mpconsts.nounit.Lf + - (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_i) + * (temperature - mpconsts.nounit.T0)) + + @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. @@ -70,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 @@ -406,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`, @@ -432,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 -------- @@ -453,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) + if max_iters or eps: + _warnings.warn( + 'max_iters, eps arguments unused and will be deprecated in a future version.', + PendingDeprecationWarning) - # 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) + 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 - # 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) + 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 - return lcl_p, globals()['dewpoint']._nounit(vapor_pressure._nounit(lcl_p, w)) + 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 @@ -1300,7 +1538,7 @@ def saturation_vapor_pressure(temperature): >>> from metpy.calc import saturation_vapor_pressure >>> from metpy.units import units >>> saturation_vapor_pressure(25 * units.degC).to('hPa') - + See Also -------- @@ -1311,14 +1549,21 @@ def saturation_vapor_pressure(temperature): Instead of temperature, dewpoint may be used in order to calculate the actual (ambient) water vapor (partial) pressure. - The formula used is that from [Bolton1980]_ for T in degrees Celsius: - - .. math:: 6.112 e^\frac{17.67T}{T + 243.5} + Implemented solution from [Ambaum2020]_, Eq. 13, + .. math:: e = e_{s0} \frac{T_0}{T}^{(c_{pl} - c_{pv}) / R_v} \exp{ + \frac{L_0}{R_v T_0} - \frac{L}{R_v T}} """ - # Converted from original in terms of C to use kelvin. - return mpconsts.nounit.sat_pressure_0c * np.exp( - 17.67 * (temperature - 273.15) / (temperature - 29.65) + latent_heat = water_latent_heat_vaporization._nounit(temperature) + temp_ratio = mpconsts.nounit.T0 / temperature + heat_power = (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v) / mpconsts.nounit.Rv + exp_term = ((mpconsts.nounit.Lv / mpconsts.nounit.T0 - latent_heat / temperature) + / mpconsts.nounit.Rv) + + return ( + mpconsts.nounit.sat_pressure_0c + * temp_ratio ** heat_power + * np.exp(exp_term) ) diff --git a/src/metpy/constants/__init__.py b/src/metpy/constants/__init__.py index 9c2aab04b8a..ee8228be7ee 100644 --- a/src/metpy/constants/__init__.py +++ b/src/metpy/constants/__init__.py @@ -22,20 +22,22 @@ Water ----- -======================= ================ ========== ============================ ==================================================== -Name Symbol Short Name Units Description ------------------------ ---------------- ---------- ---------------------------- ---------------------------------------------------- -water_molecular_weight :math:`M_w` Mw :math:`\text{g mol}^{-1}` Molecular weight of water [5]_ -water_gas_constant :math:`R_v` Rv :math:`\text{J (K kg)}^{-1}` Gas constant for water vapor [2]_ [5]_ -density_water :math:`\rho_l` rho_l :math:`\text{kg m}^{-3}` Maximum recommended density of liquid water, 0-40C [5]_ -wv_specific_heat_press :math:`C_{pv}` Cp_v :math:`\text{J (K kg)}^{-1}` Specific heat at constant pressure for water vapor -wv_specific_heat_vol :math:`C_{vv}` Cv_v :math:`\text{J (K kg)}^{-1}` Specific heat at constant volume for water vapor -water_specific_heat :math:`Cp_l` Cp_l :math:`\text{J (K kg)}^{-1}` Specific heat of liquid water at 0C [6]_ -water_heat_vaporization :math:`L_v` Lv :math:`\text{J kg}^{-1}` Latent heat of vaporization for liquid water at 0C [7]_ -water_heat_fusion :math:`L_f` Lf :math:`\text{J kg}^{-1}` Latent heat of fusion for liquid water at 0C [7]_ -ice_specific_heat :math:`C_{pi}` Cp_i :math:`\text{J (K kg)}^{-1}` Specific heat of ice at 0C [7]_ -density_ice :math:`\rho_i` rho_i :math:`\text{kg m}^{-3}` Density of ice at 0C -======================= ================ ========== ============================ ==================================================== +============================== ================ ========== ============================== ========================================================== +Name Symbol Short Name Units Description +------------------------------ ---------------- ---------- ------------------------------ ---------------------------------------------------------- +water_molecular_weight :math:`M_w` Mw :math:`\text{g mol}^{-1}` Molecular weight of water [5]_ +water_gas_constant :math:`R_v` Rv :math:`\text{J (K kg)}^{-1}` Gas constant for water vapor [2]_ [5]_ +density_water :math:`\rho_l` rho_l :math:`\text{kg m}^{-3}` Maximum recommended density of liquid water, 0-40C [5]_ +wv_specific_heat_press :math:`C_{pv}` Cp_v :math:`\text{J (K kg)}^{-1}` Specific heat at constant pressure for water vapor +wv_specific_heat_vol :math:`C_{vv}` Cv_v :math:`\text{J (K kg)}^{-1}` Specific heat at constant volume for water vapor +water_specific_heat :math:`C_{pl}` Cp_l :math:`\text{J (K kg)}^{-1}` Specific heat of liquid water at 0C [6]_ +water_heat_vaporization :math:`L_v` Lv :math:`\text{J kg}^{-1}` Latent heat of vaporization for liquid water at 0C [7]_ +water_heat_fusion :math:`L_f` Lf :math:`\text{J kg}^{-1}` Latent heat of fusion for liquid water at 0C [7]_ +water_heat_sublimation :math:`L_s` Ls :math:`\text{J kg}^{-1}` Latent heat of sublimation for water, Lv + Lf +ice_specific_heat :math:`C_{pi}` Cp_i :math:`\text{J (K kg)}^{-1}` Specific heat of ice at 0C [7]_ +density_ice :math:`\rho_i` rho_i :math:`\text{kg m}^{-3}` Density of ice at 0C +water_triple_point_temperature :math:`T_0` T0 :math:`\text{K}` Triple-point temperature of water [2]_ +============================== ================ ========== ============================== ========================================================== Dry Air ------- diff --git a/src/metpy/constants/default.py b/src/metpy/constants/default.py index 28c9558b01f..c3067d85aac 100644 --- a/src/metpy/constants/default.py +++ b/src/metpy/constants/default.py @@ -37,9 +37,11 @@ Cp_l = water_specific_heat = units.Quantity(4.2194, 'kJ / kg / K').to('J / kg / K') Lv = water_heat_vaporization = units.Quantity(2.50084e6, 'J / kg') Lf = water_heat_fusion = units.Quantity(3.337e5, 'J / kg') + Ls = water_heat_sublimation = Lv + Lf Cp_i = ice_specific_heat = units.Quantity(2090, 'J / kg / K') rho_i = density_ice = units.Quantity(917, 'kg / m^3') sat_pressure_0c = units.Quantity(6.112, 'millibar') + T0 = water_triple_point_temperature = units.Quantity(273.16, 'K') # Dry air Md = dry_air_molecular_weight = units.Quantity(28.96546e-3, 'kg / mol') diff --git a/src/metpy/constants/nounit.py b/src/metpy/constants/nounit.py index 25d9f2cd228..2fd6e8a78fa 100644 --- a/src/metpy/constants/nounit.py +++ b/src/metpy/constants/nounit.py @@ -7,10 +7,17 @@ from ..units import units Rd = default.Rd.m_as('m**2 / K / s**2') +Rv = default.Rv.m_as('m**2 / K / s**2') Lv = default.Lv.m_as('m**2 / s**2') +Lf = default.Lf.m_as('m**2 / s**2') +Ls = default.Ls.m_as('m**2 / s**2') Cp_d = default.Cp_d.m_as('m**2 / K / s**2') +Cp_l = default.Cp_l.m_as('m**2 / K / s**2') +Cp_v = default.Cp_v.m_as('m**2 / K / s**2') +Cp_i = default.Cp_i.m_as('m**2 / K / s**2') zero_degc = units.Quantity(0., 'degC').m_as('K') sat_pressure_0c = default.sat_pressure_0c.m_as('Pa') epsilon = default.epsilon.m_as('') kappa = default.kappa.m_as('') g = default.g.m_as('m / s**2') +T0 = default.T0.m_as('K') diff --git a/src/metpy/units.py b/src/metpy/units.py index 1328c772251..cc39517ac02 100644 --- a/src/metpy/units.py +++ b/src/metpy/units.py @@ -39,7 +39,9 @@ '[temperature]': 'K', '[dimensionless]': '', '[length]': 'm', - '[speed]': 'm s**-1' + '[speed]': 'm s**-1', + '[specific_enthalpy]': 'm**2 s**-2', + '[specific_heat_capacity]': 'm**2 s**-2 K-1' } @@ -83,6 +85,8 @@ def setup_registry(reg): reg.define('degrees_east = degree = degrees_E = degreesE = degree_east = degree_E ' '= degreeE') reg.define('dBz = 1e-18 m^3; logbase: 10; logfactor: 10 = dBZ') + reg.define('[specific_enthalpy] = [energy] / [mass]') + reg.define('[specific_heat_capacity] = [specific_enthalpy] / [temperature]') # Alias geopotential meters (gpm) to just meters reg.define('@alias meter = gpm') diff --git a/tests/calc/test_indices.py b/tests/calc/test_indices.py index 0a2982a65f9..da5d6537752 100644 --- a/tests/calc/test_indices.py +++ b/tests/calc/test_indices.py @@ -21,7 +21,7 @@ def test_precipitable_water(): """Test precipitable water with observed sounding.""" data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC') pw = precipitable_water(data['pressure'], data['dewpoint'], top=400 * units.hPa) - truth = 22.60430651 * units.millimeters + truth = 22.5880919 * units.millimeters assert_array_almost_equal(pw, truth, 4) @@ -32,7 +32,7 @@ def test_precipitable_water_no_bounds(): pressure = data['pressure'] inds = pressure >= 400 * units.hPa pw = precipitable_water(pressure[inds], dewpoint[inds]) - truth = 22.60430651 * units.millimeters + truth = 22.5880919 * units.millimeters assert_array_almost_equal(pw, truth, 4) @@ -43,7 +43,7 @@ def test_precipitable_water_bound_error(): dewpoint = np.array([25.5, 24.1, 23.1, 21.2, 21.1, 19.4, 19.2, 19.2, -87.1, -86.5, -86.5, -86.5, -88.1]) * units.degC pw = precipitable_water(pressure, dewpoint) - truth = 89.86846252697836 * units('millimeters') + truth = 89.7845131 * units('millimeters') assert_almost_equal(pw, truth, 5) @@ -59,7 +59,7 @@ def test_precipitable_water_nans(): -28.3, np.nan, -32.6, np.nan, -33.8, -35., -35.1, -38.1, -40., -43.3, -44.6, -46.4, -47., -49.2, -50.7]) * units.degC pw = precipitable_water(pressure, dewpoint) - truth = 4.003660322395436 * units.mm + truth = 3.99738502 * units.mm assert_almost_equal(pw, truth, 5) @@ -161,7 +161,7 @@ def test_precipitable_water_xarray(): press = xr.DataArray(data['pressure'].m, attrs={'units': str(data['pressure'].units)}) dewp = xr.DataArray(data['dewpoint'], dims=('press',), coords=(press,)) pw = precipitable_water(press, dewp, top=400 * units.hPa) - truth = 22.60430651 * units.millimeters + truth = 22.5880919 * units.millimeters assert_almost_equal(pw, truth) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 05add1f17a2..3dde4c723c9 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -15,16 +15,18 @@ brunt_vaisala_period, cape_cin, ccl, cross_totals, density, dewpoint, dewpoint_from_relative_humidity, dewpoint_from_specific_humidity, downdraft_cape, dry_lapse, dry_static_energy, el, - equivalent_potential_temperature, exner_function, + equivalent_potential_temperature, exner_function, galvez_davison_index, gradient_richardson_number, InvalidSoundingError, isentropic_interpolation, isentropic_interpolation_as_dataset, k_index, lcl, lfc, lifted_index, mixed_layer, mixed_layer_cape_cin, mixed_parcel, mixing_ratio, mixing_ratio_from_relative_humidity, - mixing_ratio_from_specific_humidity, moist_lapse, moist_static_energy, - most_unstable_cape_cin, most_unstable_parcel, parcel_profile, - parcel_profile_with_lcl, parcel_profile_with_lcl_as_dataset, - potential_temperature, psychrometric_vapor_pressure_wet, - relative_humidity_from_dewpoint, relative_humidity_from_mixing_ratio, + mixing_ratio_from_specific_humidity, moist_air_gas_constant, + moist_air_poisson_exponent, moist_air_specific_heat_pressure, + moist_lapse, moist_static_energy, most_unstable_cape_cin, + most_unstable_parcel, parcel_profile, parcel_profile_with_lcl, + parcel_profile_with_lcl_as_dataset, potential_temperature, + psychrometric_vapor_pressure_wet, relative_humidity_from_dewpoint, + relative_humidity_from_mixing_ratio, relative_humidity_from_specific_humidity, relative_humidity_wet_psychrometric, saturation_equivalent_potential_temperature, saturation_mixing_ratio, @@ -36,30 +38,81 @@ vapor_pressure, vertical_totals, vertical_velocity, vertical_velocity_pressure, virtual_potential_temperature, virtual_temperature, virtual_temperature_from_dewpoint, - wet_bulb_potential_temperature, wet_bulb_temperature) -from metpy.calc.thermo import _find_append_zero_crossings, galvez_davison_index + water_latent_heat_melting, water_latent_heat_sublimation, + water_latent_heat_vaporization, wet_bulb_potential_temperature, + wet_bulb_temperature) +from metpy.calc.thermo import _find_append_zero_crossings +from metpy.constants import Cp_d, kappa, Lf, Ls, Lv, Rd, T0 from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_nan, version_check) from metpy.units import is_quantity, masked_array, units +def test_moist_air_gas_constant(): + """Test calculation of gas constant for moist air.""" + q = 9 * units('g/kg') + assert_almost_equal(moist_air_gas_constant(q), 288.62 * units('J / kg / K'), 2) + assert_almost_equal(moist_air_gas_constant(0), Rd) + + +def test_moist_air_specific_heat_pressure(): + """Test calculation of specific heat for moist air.""" + q = 9 * units('g/kg') + assert_almost_equal(moist_air_specific_heat_pressure(q), 1012.36 * units('J / kg /K'), 2) + assert_almost_equal(moist_air_specific_heat_pressure(0), Cp_d) + + +def test_moist_air_poisson_exponent(): + """Test calculation of kappa for moist air.""" + q = 9 * units('g/kg') + assert_almost_equal(moist_air_poisson_exponent(q), 0.2851, 3) + assert_almost_equal(moist_air_poisson_exponent(0), kappa) + + +def test_water_latent_heat_vaporization(): + """Test temperature-dependent calculation of latent heat of vaporization for water.""" + temperature = 300 * units.K + # Divide out sig figs in results for decimal comparison + assert_almost_equal(water_latent_heat_vaporization(temperature) / 10**6, + 2.4375 * units('J / kg'), 4) + assert_almost_equal(water_latent_heat_vaporization(T0), Lv) + + +def test_water_latent_heat_sublimation(): + """Test temperature-dependent calculation of latent heat of sublimation for water.""" + temperature = 233 * units.K + # Divide out sig figs in results for decimal comparison + assert_almost_equal(water_latent_heat_sublimation(temperature) / 10**6, + 2.8438 * units('J / kg'), 4) + assert_almost_equal(water_latent_heat_sublimation(T0), Ls) + + +def test_water_latent_heat_melting(): + """Test temperature-dependent calculation of latent heat of melting for water.""" + temperature = 233 * units.K + # Divide out sig figs in results for decimal comparison + assert_almost_equal(water_latent_heat_melting(temperature) / 10**6, + 0.4192 * units('J / kg'), 4) + assert_almost_equal(water_latent_heat_melting(T0), Lf) + + def test_relative_humidity_from_dewpoint(): """Test Relative Humidity calculation.""" assert_almost_equal(relative_humidity_from_dewpoint(25. * units.degC, 15. * units.degC), - 53.80 * units.percent, 2) + 53.86 * units.percent, 2) def test_relative_humidity_from_dewpoint_with_f(): """Test Relative Humidity accepts temperature in Fahrenheit.""" assert_almost_equal(relative_humidity_from_dewpoint(70. * units.degF, 55. * units.degF), - 58.935 * units.percent, 3) + 58.970 * units.percent, 3) def test_relative_humidity_from_dewpoint_xarray(): """Test Relative Humidity with xarray data arrays (quantified and unquantified).""" temp = xr.DataArray(25., attrs={'units': 'degC'}) dewp = xr.DataArray([15.] * units.degC) - assert_almost_equal(relative_humidity_from_dewpoint(temp, dewp), 53.80 * units.percent, 2) + assert_almost_equal(relative_humidity_from_dewpoint(temp, dewp), 53.86 * units.percent, 2) def test_exner_function(): @@ -126,7 +179,7 @@ def test_moist_lapse(temp_units): starting_temp = 19.85 * units.degC temp = moist_lapse(np.array([1000., 800., 600., 500., 400.]) * units.mbar, starting_temp.to(temp_units)) - true_temp = np.array([293, 284.64, 272.81, 264.42, 252.91]) * units.kelvin + true_temp = np.array([293, 284.64, 272.8, 264.41, 252.88]) * units.kelvin assert_array_almost_equal(temp, true_temp, 2) @@ -134,7 +187,7 @@ def test_moist_lapse_ref_pressure(): """Test moist_lapse with a reference pressure.""" temp = moist_lapse(np.array([1050., 800., 600., 500., 400.]) * units.mbar, 19.85 * units.degC, 1000. * units.mbar) - true_temp = np.array([294.76, 284.64, 272.81, 264.42, 252.91]) * units.kelvin + true_temp = np.array([294.76, 284.64, 272.8, 264.41, 252.88]) * units.kelvin assert_array_almost_equal(temp, true_temp, 2) @@ -142,9 +195,9 @@ def test_moist_lapse_multiple_temps(): """Test moist_lapse with multiple starting temperatures.""" temp = moist_lapse(np.array([1050., 800., 600., 500., 400.]) * units.mbar, np.array([19.85, 25.6, 19.85]) * units.degC, 1000. * units.mbar) - true_temp = np.array([[294.76, 284.64, 272.81, 264.42, 252.91], - [300.35, 291.27, 281.05, 274.05, 264.64], - [294.76, 284.64, 272.81, 264.42, 252.91]]) * units.kelvin + true_temp = np.array([[294.76, 284.64, 272.8, 264.41, 252.88], + [300.35, 291.27, 281.04, 274.04, 264.62], + [294.76, 284.64, 272.8, 264.41, 252.88]]) * units.kelvin assert_array_almost_equal(temp, true_temp, 2) @@ -184,14 +237,14 @@ def test_moist_lapse_nan_ref_press(): def test_moist_lapse_downwards(): """Test moist_lapse when integrating downwards (#2128).""" temp = moist_lapse(units.Quantity([600, 700], 'mbar'), units.Quantity(0, 'degC')) - assert_almost_equal(temp, units.Quantity([0, 6.47748353], units.degC), 4) + assert_almost_equal(temp, units.Quantity([0, 6.47879528], units.degC), 4) @pytest.mark.parametrize('direction', (1, -1)) @pytest.mark.parametrize('start', list(range(5))) def test_moist_lapse_starting_points(start, direction): """Test moist_lapse with a variety of reference points.""" - truth = units.Quantity([20.0804315, 17.2333509, 14.0752659, 6.4774835, 0.0], + truth = units.Quantity([20.0804315, 17.23238477, 14.07345674, 6.47393223, -0.00554555], 'degC')[::direction] pressure = units.Quantity([1000, 925, 850, 700, 600], 'hPa')[::direction] temp = moist_lapse(pressure, truth[start], pressure[start]) @@ -218,8 +271,8 @@ def test_moist_lapse_failure(): def test_parcel_profile(): """Test parcel profile calculation.""" levels = np.array([1000., 900., 800., 700., 600., 500., 400.]) * units.mbar - true_prof = np.array([303.15, 294.16, 288.026, 283.073, 277.058, 269.402, - 258.966]) * units.kelvin + true_prof = np.array([303.15, 294.16, 288.01, 283.06, 277.04, 269.37, + 258.92]) * units.kelvin prof = parcel_profile(levels, 30. * units.degC, 20. * units.degC) assert_array_almost_equal(prof, true_prof, 2) @@ -231,11 +284,11 @@ def test_parcel_profile_lcl(): t = np.array([24.2, 24., 20.2, 21.6, 21.4, 20.4, 20.2, 14.4, 13.2, 13.]) * units.degC td = np.array([21.9, 22.1, 19.2, 20.5, 20.4, 18.4, 17.4, 8.4, -2.8, -3.0]) * units.degC - true_prof = np.array([297.35, 297.01, 294.5, 293.48, 292.92, 292.81, 289.79, 289.32, - 285.15, 282.59, 282.53]) * units.kelvin - true_p = np.insert(p.m, 2, 970.711) * units.mbar - true_t = np.insert(t.m, 2, 22.047) * units.degC - true_td = np.insert(td.m, 2, 20.609) * units.degC + true_prof = np.array([297.35, 297.01, 294.48, 293.47, 292.9, 292.79, 289.77, 289.3, + 285.13, 282.57, 282.51]) * units.kelvin + true_p = np.insert(p.m, 2, 970.441) * units.mbar + true_t = np.insert(t.m, 2, 22.029) * units.degC + true_td = np.insert(td.m, 2, 20.596) * units.degC pressure, temp, dewp, prof = parcel_profile_with_lcl(p, t, td) assert_almost_equal(pressure, true_p, 3) @@ -266,19 +319,19 @@ def test_parcel_profile_with_lcl_as_dataset(): { 'ambient_temperature': ( ('isobaric',), - np.insert(t.m, 2, 22.047) * units.degC, + np.insert(t.m, 2, 22.029) * units.degC, {'standard_name': 'air_temperature'} ), 'ambient_dew_point': ( ('isobaric',), - np.insert(td.m, 2, 20.609) * units.degC, + np.insert(td.m, 2, 20.596) * units.degC, {'standard_name': 'dew_point_temperature'} ), 'parcel_temperature': ( ('isobaric',), [ - 297.35, 297.01, 294.5, 293.48, 292.92, 292.81, 289.79, 289.32, 285.15, - 282.59, 282.53 + 297.35, 297.01, 294.48, 293.47, 292.90, 292.79, 289.77, 289.30, 285.13, + 282.57, 282.51 ] * units.kelvin, {'long_name': 'air_temperature_of_lifted_parcel'} ) @@ -286,7 +339,7 @@ def test_parcel_profile_with_lcl_as_dataset(): coords={ 'isobaric': ( 'isobaric', - np.insert(p.m, 2, 970.699), + np.insert(p.m, 2, 970.441), {'units': 'hectopascal', 'standard_name': 'air_pressure'} ) } @@ -304,7 +357,7 @@ def test_parcel_profile_with_lcl_as_dataset(): def test_parcel_profile_saturated(): """Test parcel_profile works when LCL in levels (issue #232).""" levels = np.array([1000., 700., 500.]) * units.mbar - true_prof = np.array([296.95, 284.381, 271.123]) * units.kelvin + true_prof = np.array([296.95, 284.35, 271.08]) * units.kelvin prof = parcel_profile(levels, 23.8 * units.degC, 23.8 * units.degC) assert_array_almost_equal(prof, true_prof, 2) @@ -313,20 +366,20 @@ def test_parcel_profile_saturated(): def test_sat_vapor_pressure(): """Test saturation_vapor_pressure calculation.""" temp = np.array([5., 10., 18., 25.]) * units.degC - real_es = np.array([8.72, 12.27, 20.63, 31.67]) * units.mbar + real_es = np.array([8.72, 12.27, 20.61, 31.62]) * units.mbar assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 2) def test_sat_vapor_pressure_scalar(): """Test saturation_vapor_pressure handles scalar values.""" - es = saturation_vapor_pressure(0 * units.degC) + es = saturation_vapor_pressure(T0) assert_almost_equal(es, 6.112 * units.mbar, 3) def test_sat_vapor_pressure_fahrenheit(): """Test saturation_vapor_pressure handles temperature in Fahrenheit.""" temp = np.array([50., 68.]) * units.degF - real_es = np.array([12.2717, 23.3695]) * units.mbar + real_es = np.array([12.2666, 23.3475]) * units.mbar assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 4) @@ -388,8 +441,8 @@ def test_vapor_pressure(): def test_lcl(): """Test LCL calculation.""" lcl_pressure, lcl_temperature = lcl(1000. * units.mbar, 30. * units.degC, 20. * units.degC) - assert_almost_equal(lcl_pressure, 864.806 * units.mbar, 2) - assert_almost_equal(lcl_temperature, 17.676 * units.degC, 2) + assert_almost_equal(lcl_pressure, 864.614 * units.mbar, 2) + assert_almost_equal(lcl_temperature, 17.658 * units.degC, 2) def test_lcl_kelvin(): @@ -397,8 +450,8 @@ def test_lcl_kelvin(): temperature = 273.09723 * units.kelvin lcl_pressure, lcl_temperature = lcl(1017.16 * units.mbar, temperature, 264.5351 * units.kelvin) - assert_almost_equal(lcl_pressure, 889.459 * units.mbar, 2) - assert_almost_equal(lcl_temperature, 262.827 * units.kelvin, 2) + assert_almost_equal(lcl_pressure, 889.218 * units.mbar, 2) + assert_almost_equal(lcl_temperature, 262.802 * units.kelvin, 2) assert lcl_temperature.units == temperature.units @@ -415,10 +468,10 @@ def test_lcl_nans(): dewp = np.array([20., 20., np.nan, 20.]) * units.degC lcl_press, lcl_temp = lcl(press, temp, dewp) - assert_array_almost_equal(lcl_press, np.array([np.nan, 836.4098648012595, - np.nan, 836.4098648012595]) * units.hPa) - assert_array_almost_equal(lcl_temp, np.array([np.nan, 18.82281982535794, - np.nan, 18.82281982535794]) * units.degC) + assert_array_almost_equal(lcl_press, np.array([np.nan, 836.2231283, + np.nan, 836.2231283]) * units.hPa) + assert_array_almost_equal(lcl_temp, np.array([np.nan, 18.8041938, + np.nan, 18.8041938]) * units.degC) def test_ccl_basic(): @@ -433,9 +486,9 @@ def test_ccl_basic(): -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) * units.degC ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) - assert_almost_equal(ccl_p, 763.006048 * units.mbar, 5) - assert_almost_equal(ccl_t, 15.429946 * units.degC, 5) - assert_almost_equal(t_c, 37.991498 * units.degC, 5) + assert_almost_equal(ccl_p, 762.81930 * units.mbar, 5) + assert_almost_equal(ccl_t, 15.41210 * units.degC, 5) + assert_almost_equal(t_c, 37.99401 * units.degC, 5) def test_ccl_nans(): @@ -450,9 +503,9 @@ def test_ccl_nans(): -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) * units.degC ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) - assert_almost_equal(ccl_p, 763.006048 * units.mbar, 5) - assert_almost_equal(ccl_t, 15.429946 * units.degC, 5) - assert_almost_equal(t_c, 37.991498 * units.degC, 5) + assert_almost_equal(ccl_p, 762.81930 * units.mbar, 5) + assert_almost_equal(ccl_t, 15.41210 * units.degC, 5) + assert_almost_equal(t_c, 37.99401 * units.degC, 5) def test_ccl_unit(): @@ -468,9 +521,9 @@ def test_ccl_unit(): -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) + 273.15) * units.kelvin ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) - assert_almost_equal(ccl_p, (763.006048 * 100) * units.Pa, 3) - assert_almost_equal(ccl_t, (15.429946 + 273.15) * units.kelvin, 3) - assert_almost_equal(t_c, (37.991498 + 273.15) * units.kelvin, 3) + assert_almost_equal(ccl_p, (762.81930 * 100) * units.Pa, 3) + assert_almost_equal(ccl_t, (15.41210 + 273.15) * units.kelvin, 3) + assert_almost_equal(t_c, (37.99401 + 273.15) * units.kelvin, 3) assert ccl_p.units == pressure.units assert ccl_t.units == temperature.units @@ -498,19 +551,19 @@ def test_multiple_ccl(): -36.7, -28.9, -28.5, -22.5]) * units.degC ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) - assert_almost_equal(ccl_p, 680.191653 * units.mbar, 5) - assert_almost_equal(ccl_t, -0.204408 * units.degC, 5) - assert_almost_equal(t_c, 30.8678258 * units.degC, 5) + assert_almost_equal(ccl_p, 680.087037 * units.mbar, 5) + assert_almost_equal(ccl_t, -0.213239 * units.degC, 5) + assert_almost_equal(t_c, 30.8713506 * units.degC, 5) ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, which='bottom') - assert_almost_equal(ccl_p, 886.835325 * units.mbar, 5) - assert_almost_equal(ccl_t, 3.500840 * units.degC, 5) - assert_almost_equal(t_c, 12.5020423 * units.degC, 5) + assert_almost_equal(ccl_p, 886.724147 * units.mbar, 5) + assert_almost_equal(ccl_t, 3.492138 * units.degC, 5) + assert_almost_equal(t_c, 12.5032889 * units.degC, 5) ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, which='all') - assert_array_almost_equal(ccl_p, np.array([886.835325, 680.191653]) * units.mbar, 5) - assert_array_almost_equal(ccl_t, np.array([3.500840, -0.204408]) * units.degC, 5) - assert_array_almost_equal(t_c, np.array([12.5020423, 30.8678258]) * units.degC, 5) + assert_array_almost_equal(ccl_p, np.array([886.724147, 680.087037]) * units.mbar, 5) + assert_array_almost_equal(ccl_t, np.array([3.492138, -0.213239]) * units.degC, 5) + assert_array_almost_equal(t_c, np.array([12.5032889, 30.8713506]) * units.degC, 5) def test_ccl_with_ml(): @@ -537,11 +590,11 @@ def test_ccl_with_ml(): mixed_layer_depth=500 * units.m, which='all') assert_array_almost_equal(ccl_p, np.array( - [850.600930, 784.325312, 737.767377, 648.076147]) * units.mbar, 5) + [850.438506, 784.136599, 782.671594, 737.637719, 647.950662]) * units.mbar, 5) assert_array_almost_equal(ccl_t, np.array( - [0.840118, -0.280299, -1.118757, -2.875716]) * units.degC, 5) + [0.829278, -0.291727, -0.317423, -1.129214, -2.886271]) * units.degC, 5) assert_array_almost_equal(t_c, np.array( - [13.146845, 18.661621, 22.896152, 32.081388]) * units.degC, 5) + [13.151139, 18.669463, 18.797927, 22.899638, 32.086356]) * units.degC, 5) def test_lfc_basic(): @@ -550,8 +603,8 @@ def test_lfc_basic(): temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -49.]) * units.celsius dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints) - assert_almost_equal(lfc_pressure, 727.371 * units.mbar, 2) - assert_almost_equal(lfc_temp, 9.705 * units.celsius, 2) + assert_almost_equal(lfc_pressure, 727.079 * units.mbar, 2) + assert_almost_equal(lfc_temp, 9.672 * units.celsius, 2) def test_lfc_kelvin(): @@ -562,8 +615,8 @@ def test_lfc_kelvin(): dewpoint = (np.array([19., -11.2, -10.8, -10.4, -10., -53.2] ) + 273.15) * units.kelvin lfc_pressure, lfc_temp = lfc(pressure, temperature, dewpoint) - assert_almost_equal(lfc_pressure, 727.371 * units.mbar, 2) - assert_almost_equal(lfc_temp, 9.705 * units.degC, 2) + assert_almost_equal(lfc_pressure, 727.079 * units.mbar, 2) + assert_almost_equal(lfc_temp, 9.672 * units.degC, 2) assert lfc_temp.units == temperature.units @@ -575,8 +628,8 @@ def test_lfc_ml(): __, t_mixed, td_mixed = mixed_parcel(levels, temperatures, dewpoints) mixed_parcel_prof = parcel_profile(levels, t_mixed, td_mixed) lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints, mixed_parcel_prof) - assert_almost_equal(lfc_pressure, 601.225 * units.mbar, 2) - assert_almost_equal(lfc_temp, -1.90688 * units.degC, 2) + assert_almost_equal(lfc_pressure, 599.020 * units.mbar, 2) + assert_almost_equal(lfc_temp, -2.122 * units.degC, 2) def test_lfc_ml2(): @@ -616,8 +669,8 @@ def test_lfc_ml2(): __, t_mixed, td_mixed = mixed_parcel(levels, temperatures, dewpoints) mixed_parcel_prof = parcel_profile(levels, t_mixed, td_mixed) lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints, mixed_parcel_prof, td_mixed) - assert_almost_equal(lfc_pressure, 962.34 * units.mbar, 2) - assert_almost_equal(lfc_temp, 0.767 * units.degC, 2) + assert_almost_equal(lfc_pressure, 962.07 * units.mbar, 2) + assert_almost_equal(lfc_temp, 0.745 * units.degC, 2) def test_lfc_intersection(): @@ -629,7 +682,7 @@ def test_lfc_intersection(): _, mlt, mltd = mixed_parcel(p, t, td) ml_profile = parcel_profile(p, mlt, mltd) mllfc_p, mllfc_t = lfc(p, t, td, ml_profile, mltd) - assert_almost_equal(mllfc_p, 981.620 * units.hPa, 2) + assert_almost_equal(mllfc_p, 981.584 * units.hPa, 2) assert_almost_equal(mllfc_t, 272.045 * units.kelvin, 2) @@ -652,8 +705,8 @@ def test_lfc_inversion(): dewpoints = np.array([20.4, 0.4, -0.5, -4.3, -8., -8.2, -9., -23.9, -33.3, -54.1, -63.5]) * units.celsius lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints) - assert_almost_equal(lfc_pressure, 705.8806 * units.mbar, 2) - assert_almost_equal(lfc_temp, 10.6232 * units.celsius, 2) + assert_almost_equal(lfc_pressure, 705.5807 * units.mbar, 2) + assert_almost_equal(lfc_temp, 10.5875 * units.celsius, 2) def test_lfc_equals_lcl(): @@ -665,7 +718,7 @@ def test_lfc_equals_lcl(): dewpoints = np.array([18.4, 18.1, 16.6, 15.4, 13.2, 11.4, 9.6, 8.8, 0., -18.6, -22.9]) * units.celsius lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints) - assert_almost_equal(lfc_pressure, 777.0786 * units.mbar, 2) + assert_almost_equal(lfc_pressure, 776.935 * units.mbar, 2) assert_almost_equal(lfc_temp, 15.8714 * units.celsius, 2) @@ -675,8 +728,8 @@ def test_lfc_profile_nan(): temperatures = np.array([22.2, 14.6, np.nan, 9.4, 7., -38.]) * units.degC dewpoints = np.array([19., -11.2, -10.8, -10.4, np.nan, -53.2]) * units.degC lfc_pressure, lfc_temperature = lfc(levels, temperatures, dewpoints) - assert_almost_equal(lfc_pressure, 727.3365 * units.mbar, 3) - assert_almost_equal(lfc_temperature, 9.6977 * units.degC, 3) + assert_almost_equal(lfc_pressure, 727.0471 * units.mbar, 3) + assert_almost_equal(lfc_temperature, 9.6694 * units.degC, 3) def test_lfc_profile_nan_with_parcel_profile(): @@ -686,8 +739,8 @@ def test_lfc_profile_nan_with_parcel_profile(): dewpoints = np.array([19., -11.2, -10.8, -10.4, np.nan, -53.2]) * units.degC parcel_temps = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') lfc_pressure, lfc_temperature = lfc(levels, temperatures, dewpoints, parcel_temps) - assert_almost_equal(lfc_pressure, 727.3365 * units.mbar, 3) - assert_almost_equal(lfc_temperature, 9.6977 * units.degC, 3) + assert_almost_equal(lfc_pressure, 727.0471 * units.mbar, 3) + assert_almost_equal(lfc_temperature, 9.6694 * units.degC, 3) def test_sensitive_sounding(): @@ -704,12 +757,12 @@ def test_sensitive_sounding(): -20.3, -29.1, -27.7, -24.9, -39.5, -41.9, -51.9, -60.7, -62.7, -65.1, -71.9], 'degC') lfc_pressure, lfc_temp = lfc(p, t, td) - assert_almost_equal(lfc_pressure, 952.8445 * units.mbar, 2) - assert_almost_equal(lfc_temp, 20.94469 * units.degC, 2) + assert_almost_equal(lfc_pressure, 952.4893 * units.mbar, 2) + assert_almost_equal(lfc_temp, 20.9179 * units.degC, 2) pos, neg = surface_based_cape_cin(p, t, td) - assert_almost_equal(pos, 0.106791 * units('J/kg'), 3) - assert_almost_equal(neg, -282.620677 * units('J/kg'), 3) + assert_almost_equal(pos, 0.002948 * units('J/kg'), 3) + assert_almost_equal(neg, -286.683725 * units('J/kg'), 3) def test_lfc_sfc_precision(): @@ -782,7 +835,7 @@ def test_equivalent_potential_temperature(): t = 293. * units.kelvin td = 280. * units.kelvin ept = equivalent_potential_temperature(p, t, td) - assert_almost_equal(ept, 311.18586467284007 * units.kelvin, 3) + assert_almost_equal(ept, 311.17705883845116 * units.kelvin, 3) def test_equivalent_potential_temperature_masked(): @@ -809,7 +862,7 @@ def test_wet_bulb_potential_temperature(): t = 293. * units.kelvin td = 291. * units.kelvin wpt = wet_bulb_potential_temperature(p, t, td) - assert_almost_equal(wpt, 291.65839705486565 * units.kelvin, 3) + assert_almost_equal(wpt, 291.65007188488414 * units.kelvin, 3) def test_saturation_equivalent_potential_temperature(): @@ -819,7 +872,7 @@ def test_saturation_equivalent_potential_temperature(): s_ept = saturation_equivalent_potential_temperature(p, t) # 299.096584 comes from equivalent_potential_temperature(p,t,t) # where dewpoint and temperature are equal, which means saturations. - assert_almost_equal(s_ept, 299.10542 * units.kelvin, 3) + assert_almost_equal(s_ept, 299.09407 * units.kelvin, 3) def test_saturation_equivalent_potential_temperature_masked(): @@ -828,7 +881,7 @@ def test_saturation_equivalent_potential_temperature_masked(): t = units.Quantity(np.ma.array([293., 294., 295.]), units.kelvin) s_ept = saturation_equivalent_potential_temperature(p, t) expected = units.Quantity( - np.ma.array([335.02750, 338.95813, 343.08740]), + np.ma.array([334.98378, 338.90573, 343.02461]), units.kelvin ) assert is_quantity(s_ept) @@ -850,7 +903,7 @@ def test_virtual_temperature_from_dewpoint(): t = 30 * units.degC td = 25 * units.degC tv = virtual_temperature_from_dewpoint(p, t, td) - assert_almost_equal(tv, 33.6740 * units.degC, 3) + assert_almost_equal(tv, 33.6680 * units.degC, 3) def test_virtual_potential_temperature(): @@ -877,8 +930,8 @@ def test_el(): temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius el_pressure, el_temperature = el(levels, temperatures, dewpoints) - assert_almost_equal(el_pressure, 471.83286 * units.mbar, 3) - assert_almost_equal(el_temperature, -11.5603 * units.degC, 3) + assert_almost_equal(el_pressure, 476.07144 * units.mbar, 3) + assert_almost_equal(el_temperature, -11.1395 * units.degC, 3) def test_el_kelvin(): @@ -887,8 +940,8 @@ def test_el_kelvin(): temperatures = (np.array([22.2, 14.6, 12., 9.4, 7., -38.]) + 273.15) * units.kelvin dewpoints = (np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) + 273.15) * units.kelvin el_pressure, el_temp = el(levels, temperatures, dewpoints) - assert_almost_equal(el_pressure, 471.8329 * units.mbar, 3) - assert_almost_equal(el_temp, -11.5603 * units.degC, 3) + assert_almost_equal(el_pressure, 476.07144 * units.mbar, 3) + assert_almost_equal(el_temp, -11.139 * units.degC, 3) assert el_temp.units == temperatures.units @@ -900,8 +953,8 @@ def test_el_ml(): __, t_mixed, td_mixed = mixed_parcel(levels, temperatures, dewpoints) mixed_parcel_prof = parcel_profile(levels, t_mixed, td_mixed) el_pressure, el_temperature = el(levels, temperatures, dewpoints, mixed_parcel_prof) - assert_almost_equal(el_pressure, 350.0561 * units.mbar, 3) - assert_almost_equal(el_temperature, -28.36156 * units.degC, 3) + assert_almost_equal(el_pressure, 350.5072 * units.mbar, 3) + assert_almost_equal(el_temperature, -28.32911 * units.degC, 3) def test_no_el(): @@ -959,8 +1012,8 @@ def test_el_lfc_equals_lcl(): -54.4, -68., -70.1, -70., -70., -70., -70., -70., -70., -70., -70., -70., -70., -70., -70., -70., -70.]) * units.celsius el_pressure, el_temperature = el(levels, temperatures, dewpoints) - assert_almost_equal(el_pressure, 175.7663 * units.mbar, 3) - assert_almost_equal(el_temperature, -57.03994 * units.degC, 3) + assert_almost_equal(el_pressure, 175.8445 * units.mbar, 3) + assert_almost_equal(el_temperature, -57.05514 * units.degC, 3) def test_el_small_surface_instability(): @@ -1029,8 +1082,8 @@ def test_el_profile_nan(): temperatures = np.array([22.2, 14.6, np.nan, 9.4, 7., -38.]) * units.degC dewpoints = np.array([19., -11.2, -10.8, -10.4, np.nan, -53.2]) * units.degC el_pressure, el_temperature = el(levels, temperatures, dewpoints) - assert_almost_equal(el_pressure, 673.0104 * units.mbar, 3) - assert_almost_equal(el_temperature, 5.8853 * units.degC, 3) + assert_almost_equal(el_pressure, 678.5703 * units.mbar, 3) + assert_almost_equal(el_temperature, 6.2790 * units.degC, 3) def test_el_profile_nan_with_parcel_profile(): @@ -1040,8 +1093,8 @@ def test_el_profile_nan_with_parcel_profile(): dewpoints = np.array([19., -11.2, -10.8, -10.4, np.nan, -53.2]) * units.degC parcel_temps = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') el_pressure, el_temperature = el(levels, temperatures, dewpoints, parcel_temps) - assert_almost_equal(el_pressure, 673.0104 * units.mbar, 3) - assert_almost_equal(el_temperature, 5.8853 * units.degC, 3) + assert_almost_equal(el_pressure, 678.5703 * units.mbar, 3) + assert_almost_equal(el_temperature, 6.2790 * units.degC, 3) def test_wet_psychrometric_vapor_pressure(): @@ -1051,7 +1104,7 @@ def test_wet_psychrometric_vapor_pressure(): wet_bulb_temperature = 18. * units.degC psychrometric_vapor_pressure = psychrometric_vapor_pressure_wet(p, dry_bulb_temperature, wet_bulb_temperature) - assert_almost_equal(psychrometric_vapor_pressure, 19.3673 * units.mbar, 3) + assert_almost_equal(psychrometric_vapor_pressure, 19.3518 * units.mbar, 3) def test_wet_psychrometric_rh(): @@ -1061,7 +1114,7 @@ def test_wet_psychrometric_rh(): wet_bulb_temperature = 18. * units.degC psychrometric_rh = relative_humidity_wet_psychrometric(p, dry_bulb_temperature, wet_bulb_temperature) - assert_almost_equal(psychrometric_rh, 82.8747 * units.percent, 3) + assert_almost_equal(psychrometric_rh, 82.8861 * units.percent, 3) def test_wet_psychrometric_rh_kwargs(): @@ -1073,7 +1126,7 @@ def test_wet_psychrometric_rh_kwargs(): psychrometric_rh = relative_humidity_wet_psychrometric(p, dry_bulb_temperature, wet_bulb_temperature, psychrometer_coefficient=coeff) - assert_almost_equal(psychrometric_rh, 82.9701 * units.percent, 3) + assert_almost_equal(psychrometric_rh, 82.9816 * units.percent, 3) def test_mixing_ratio_from_relative_humidity(): @@ -1091,7 +1144,7 @@ def test_rh_mixing_ratio(): temperature = 20. * units.degC w = 0.012 * units.dimensionless rh = relative_humidity_from_mixing_ratio(p, temperature, w) - assert_almost_equal(rh, 82.0709069 * units.percent, 3) + assert_almost_equal(rh, 82.1482060 * units.percent, 3) def test_mixing_ratio_from_specific_humidity(): @@ -1128,7 +1181,7 @@ def test_rh_specific_humidity(): temperature = 20. * units.degC q = 0.012 * units.dimensionless rh = relative_humidity_from_specific_humidity(p, temperature, q) - assert_almost_equal(rh, 83.0486264 * units.percent, 3) + assert_almost_equal(rh, 83.1268463 * units.percent, 3) def test_cape_cin(): @@ -1138,8 +1191,8 @@ def test_cape_cin(): dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 215.056976 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -9.94798721 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 210.952968 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -10.37750324 * units('joule / kilogram'), 2) def test_cape_cin_no_el(): @@ -1149,8 +1202,8 @@ def test_cape_cin_no_el(): dewpoint = np.array([19., -11.2, -10.8, -10.4]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC') cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 12.74623773 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -9.947987213 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 12.52109932 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -10.377503242 * units('joule / kilogram'), 2) def test_cape_cin_no_lfc(): @@ -1547,8 +1600,8 @@ def test_surface_based_cape_cin(array_class): temperature = array_class([22.2, 14.6, 12., 9.4, 7., -38.], units.celsius) dewpoint = array_class([19., -11.2, -10.8, -10.4, -10., -53.2], units.celsius) cape, cin = surface_based_cape_cin(p, temperature, dewpoint) - assert_almost_equal(cape, 215.05697634 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -33.0633599455 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 210.95296816 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -33.6620801701 * units('joule / kilogram'), 2) def test_surface_based_cape_cin_with_xarray(): @@ -1571,8 +1624,8 @@ def test_surface_based_cape_cin_with_xarray(): data['temperature'], data['dewpoint'] ) - assert_almost_equal(cape, 215.056976346 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -33.0633599455 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 210.95296816 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -33.6620801701 * units('joule / kilogram'), 2) def test_profile_with_nans(): @@ -1612,8 +1665,8 @@ def test_most_unstable_cape_cin_surface(): temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 215.056976346 * units('joule / kilogram'), 2) - assert_almost_equal(mucin, -33.0633599455 * units('joule / kilogram'), 2) + assert_almost_equal(mucape, 210.95296816 * units('joule / kilogram'), 2) + assert_almost_equal(mucin, -33.6620801701 * units('joule / kilogram'), 2) def test_most_unstable_cape_cin(): @@ -1622,8 +1675,8 @@ def test_most_unstable_cape_cin(): temperature = np.array([18.2, 22.2, 17.4, 10., 0., 15]) * units.celsius dewpoint = np.array([19., 19., 14.3, 0., -10., 0.]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 173.749389796 * units('joule / kilogram'), 4) - assert_almost_equal(mucin, -20.968278741 * units('joule / kilogram'), 4) + assert_almost_equal(mucape, 173.492800262 * units('joule / kilogram'), 4) + assert_almost_equal(mucin, -21.159922894 * units('joule / kilogram'), 4) def test_mixed_parcel(): @@ -1643,8 +1696,8 @@ def test_mixed_layer_cape_cin(multiple_intersections): """Test the calculation of mixed layer cape/cin.""" pressure, temperature, dewpoint = multiple_intersections mlcape, mlcin = mixed_layer_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mlcape, 1132.706800436 * units('joule / kilogram'), 2) - assert_almost_equal(mlcin, -13.4809966289 * units('joule / kilogram'), 2) + assert_almost_equal(mlcape, 1120.723057388 * units('joule / kilogram'), 2) + assert_almost_equal(mlcin, -13.6568786318 * units('joule / kilogram'), 2) def test_mixed_layer_cape_cin_bottom_pressure(multiple_intersections): @@ -1652,8 +1705,8 @@ def test_mixed_layer_cape_cin_bottom_pressure(multiple_intersections): pressure, temperature, dewpoint = multiple_intersections mlcape_middle, mlcin_middle = mixed_layer_cape_cin(pressure, temperature, dewpoint, parcel_start_pressure=903 * units.hPa) - assert_almost_equal(mlcape_middle, 1177.86 * units('joule / kilogram'), 2) - assert_almost_equal(mlcin_middle, -37. * units('joule / kilogram'), 2) + assert_almost_equal(mlcape_middle, 1166.88 * units('joule / kilogram'), 2) + assert_almost_equal(mlcin_middle, -37.71 * units('joule / kilogram'), 2) def test_dcape(): @@ -1813,14 +1866,14 @@ def test_wet_bulb_temperature(temp_units): temp = 25 * units.degC dewp = 15 * units.degC val = wet_bulb_temperature(1000 * units.hPa, temp.to(temp_units), dewp.to(temp_units)) - truth = 18.3432116 * units.degC # 18.59 from NWS calculator + truth = 18.3395910 * units.degC # 18.59 from NWS calculator assert_almost_equal(val, truth, 5) def test_wet_bulb_temperature_saturated(): """Test wet bulb calculation works properly with saturated conditions.""" val = wet_bulb_temperature(850. * units.hPa, 17.6 * units.degC, 17.6 * units.degC) - assert_almost_equal(val, 17.6 * units.degC, 7) + assert_almost_equal(val, 17.5918788 * units.degC, 7) def test_wet_bulb_temperature_numpy_scalars(): @@ -1829,7 +1882,7 @@ def test_wet_bulb_temperature_numpy_scalars(): temperature = units.Quantity(np.float32(25), 'degC') dewpoint = units.Quantity(np.float32(15), 'degC') val = wet_bulb_temperature(pressure, temperature, dewpoint) - truth = 18.3432116 * units.degC + truth = 18.3395909 * units.degC assert_almost_equal(val, truth, 5) @@ -1839,7 +1892,7 @@ def test_wet_bulb_temperature_1d(): temperatures = [25, 20, 15] * units.degC dewpoints = [20, 15, 10] * units.degC val = wet_bulb_temperature(pressures, temperatures, dewpoints) - truth = [21.44487, 16.73673, 12.06554] * units.degC + truth = [21.43577, 16.73217, 12.06271] * units.degC # 21.58, 16.86, 12.18 from NWS Calculator assert_array_almost_equal(val, truth, 5) @@ -1853,8 +1906,8 @@ def test_wet_bulb_temperature_2d(): dewpoints = [[20, 15, 10], [19, 14, 9]] * units.degC val = wet_bulb_temperature(pressures, temperatures, dewpoints) - truth = [[21.44487, 16.73673, 12.06554], - [20.50205, 15.80108, 11.13603]] * units.degC + truth = [[21.43577, 16.73217, 12.06271], + [20.49416, 15.79706, 11.13332]] * units.degC # 21.58, 16.86, 12.18 # 20.6, 15.9, 11.2 from NWS Calculator assert_array_almost_equal(val, truth, 5) @@ -1884,13 +1937,13 @@ def test_wet_bulb_nan(temp_units): -14.55002441, -11.74678955, -25.58999634] * units.degC val = wet_bulb_temperature(pressure, temperature.to(temp_units), dewpoint.to(temp_units)) - truth = [19.238071735308814, 18.033294139060633, 16.65179610640866, 14.341431051131467, - 10.713278865013098, 7.2703265785039, 4.63234087372236, 1.8324379627895773, - -1.0154897545814394, -4.337334561885717, -8.23596994210175, -11.902727397896111, - -15.669313544076992, -20.78875735056887, -27.342629368898884, -33.42946313024462, - -41.8026422159221, -46.17279371976847, -48.99179569697857, -60.13602538549741, - -74.63516192059605, -79.80028104362006, -74.59295016865613, np.nan, - -59.20059644897026, -55.76040402608365, np.nan, -49.692666440433335, np.nan, + truth = [19.23045181, 18.02610206, 16.64537114, 14.33739256, + 10.71130716, 7.26884942, 4.63084162, 1.83105242, + -1.01674304, -4.33748074, -8.23285723, -11.89759095, + -15.66482285, -20.78620567, -27.34145216, -33.42872025, + -41.8065101, -46.17511911, -49.00112528, -60.1408904, + -74.63591645, -79.79987632, -74.59567764, np.nan, + -59.23708051, -55.82144907, np.nan, -49.8085793, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan] * units.degC assert_array_almost_equal(val, truth, 5) @@ -1993,7 +2046,7 @@ def test_dewpoint_specific_humidity_xarray(index_xarray_data): p = index_xarray_data.isobaric q = specific_humidity_from_dewpoint(p, index_xarray_data.dewpoint) td = dewpoint_from_specific_humidity(p, specific_humidity=q) - assert_array_almost_equal(td, index_xarray_data.dewpoint) + assert_array_almost_equal(td, index_xarray_data.dewpoint, decimal=1) def test_lfc_not_below_lcl(): @@ -2009,8 +2062,8 @@ def test_lfc_not_below_lcl(): 8.6, 8.1, 7.6, 7., 6.5, 6., 5.4]) * units.degC lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints) # Before patch, LFC pressure would show 1000.5912165339967 hPa - assert_almost_equal(lfc_pressure, 811.618879 * units.mbar, 3) - assert_almost_equal(lfc_temp, 6.48644650 * units.celsius, 3) + assert_almost_equal(lfc_pressure, 810.908808 * units.mbar, 3) + assert_almost_equal(lfc_temp, 6.44247208 * units.celsius, 3) @pytest.fixture @@ -2046,10 +2099,10 @@ def test_multiple_lfcs_simple(multiple_intersections): lfc_pressure_bottom, lfc_temp_bottom = lfc(levels, temperatures, dewpoints, which='bottom') lfc_pressure_all, _ = lfc(levels, temperatures, dewpoints, which='all') - assert_almost_equal(lfc_pressure_top, 705.3534497 * units.mbar, 3) - assert_almost_equal(lfc_temp_top, 4.884899 * units.degC, 3) - assert_almost_equal(lfc_pressure_bottom, 884.1478935 * units.mbar, 3) - assert_almost_equal(lfc_temp_bottom, 13.95706975 * units.degC, 3) + assert_almost_equal(lfc_pressure_top, 704.9853394 * units.mbar, 3) + assert_almost_equal(lfc_temp_top, 4.851675 * units.degC, 3) + assert_almost_equal(lfc_pressure_bottom, 883.5277233 * units.mbar, 3) + assert_almost_equal(lfc_temp_bottom, 13.92200100 * units.degC, 3) assert_almost_equal(len(lfc_pressure_all), 2, 0) @@ -2057,8 +2110,8 @@ def test_multiple_lfs_wide(multiple_intersections): """Test 'wide' LFC for sounding with multiple LFCs.""" levels, temperatures, dewpoints = multiple_intersections lfc_pressure_wide, lfc_temp_wide = lfc(levels, temperatures, dewpoints, which='wide') - assert_almost_equal(lfc_pressure_wide, 705.3534497 * units.hPa, 3) - assert_almost_equal(lfc_temp_wide, 4.88489902 * units.degC, 3) + assert_almost_equal(lfc_pressure_wide, 704.9853394 * units.hPa, 3) + assert_almost_equal(lfc_temp_wide, 4.8516747 * units.degC, 3) def test_invalid_which(multiple_intersections): @@ -2082,10 +2135,10 @@ def test_multiple_els_simple(multiple_intersections): el_pressure_top, el_temp_top = el(levels, temperatures, dewpoints) el_pressure_bottom, el_temp_bottom = el(levels, temperatures, dewpoints, which='bottom') el_pressure_all, _ = el(levels, temperatures, dewpoints, which='all') - assert_almost_equal(el_pressure_top, 228.151524 * units.mbar, 3) - assert_almost_equal(el_temp_top, -56.810153566 * units.degC, 3) - assert_almost_equal(el_pressure_bottom, 849.7998957 * units.mbar, 3) - assert_almost_equal(el_temp_bottom, 12.42268288 * units.degC, 3) + assert_almost_equal(el_pressure_top, 228.290837 * units.mbar, 3) + assert_almost_equal(el_temp_top, -56.806956 * units.degC, 3) + assert_almost_equal(el_pressure_bottom, 849.85394 * units.mbar, 3) + assert_almost_equal(el_temp_bottom, 12.4165561 * units.degC, 3) assert_almost_equal(len(el_pressure_all), 2, 0) @@ -2093,24 +2146,24 @@ def test_multiple_el_wide(multiple_intersections): """Test 'wide' EL for sounding with multiple ELs.""" levels, temperatures, dewpoints = multiple_intersections el_pressure_wide, el_temp_wide = el(levels, temperatures, dewpoints, which='wide') - assert_almost_equal(el_pressure_wide, 228.151524 * units.hPa, 3) - assert_almost_equal(el_temp_wide, -56.81015357 * units.degC, 3) + assert_almost_equal(el_pressure_wide, 228.290837 * units.hPa, 3) + assert_almost_equal(el_temp_wide, -56.8069560 * units.degC, 3) def test_muliple_el_most_cape(multiple_intersections): """Test 'most_cape' EL for sounding with multiple ELs.""" levels, temperatures, dewpoints = multiple_intersections el_pressure_wide, el_temp_wide = el(levels, temperatures, dewpoints, which='most_cape') - assert_almost_equal(el_pressure_wide, 228.151524 * units.hPa, 3) - assert_almost_equal(el_temp_wide, -56.81015356 * units.degC, 3) + assert_almost_equal(el_pressure_wide, 228.290837 * units.hPa, 3) + assert_almost_equal(el_temp_wide, -56.8069560 * units.degC, 3) def test_muliple_lfc_most_cape(multiple_intersections): """Test 'most_cape' LFC for sounding with multiple LFCs.""" levels, temperatures, dewpoints = multiple_intersections lfc_pressure_wide, lfc_temp_wide = lfc(levels, temperatures, dewpoints, which='most_cape') - assert_almost_equal(lfc_pressure_wide, 705.35344968 * units.hPa, 3) - assert_almost_equal(lfc_temp_wide, 4.88489902 * units.degC, 3) + assert_almost_equal(lfc_pressure_wide, 704.98533942 * units.hPa, 3) + assert_almost_equal(lfc_temp_wide, 4.85167474 * units.degC, 3) def test_el_lfc_most_cape_bottom(): @@ -2123,10 +2176,10 @@ def test_el_lfc_most_cape_bottom(): -6.9, -9.5, -12., -14.6, -15.8]) * units.degC lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints, which='most_cape') el_pressure, el_temp = el(levels, temperatures, dewpoints, which='most_cape') - assert_almost_equal(lfc_pressure, 714.358842 * units.hPa, 3) - assert_almost_equal(lfc_temp, 5.415421178 * units.degC, 3) - assert_almost_equal(el_pressure, 658.61917328 * units.hPa, 3) - assert_almost_equal(el_temp, 1.951095056 * units.degC, 3) + assert_almost_equal(lfc_pressure, 714.147342 * units.hPa, 3) + assert_almost_equal(lfc_temp, 5.39186322 * units.degC, 3) + assert_almost_equal(el_pressure, 659.284276 * units.hPa, 3) + assert_almost_equal(el_temp, 1.98172061 * units.degC, 3) def test_cape_cin_top_el_lfc(multiple_intersections): @@ -2134,8 +2187,8 @@ def test_cape_cin_top_el_lfc(multiple_intersections): levels, temperatures, dewpoints = multiple_intersections parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='top') - assert_almost_equal(cape, 1345.18884959 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -35.179268355 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 1335.95888464 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -35.816974306 * units('joule / kilogram'), 3) def test_cape_cin_bottom_el_lfc(multiple_intersections): @@ -2143,8 +2196,8 @@ def test_cape_cin_bottom_el_lfc(multiple_intersections): levels, temperatures, dewpoints = multiple_intersections parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_el='bottom') - assert_almost_equal(cape, 4.57262449 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -5.9471237534 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 4.43749259 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -6.08659801232 * units('joule / kilogram'), 3) def test_cape_cin_wide_el_lfc(multiple_intersections): @@ -2153,8 +2206,8 @@ def test_cape_cin_wide_el_lfc(multiple_intersections): parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='wide', which_el='wide') - assert_almost_equal(cape, 1345.1888496 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -35.179268355 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 1335.9588846 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -35.816974306 * units('joule / kilogram'), 3) def test_cape_cin_custom_profile(): @@ -2164,7 +2217,7 @@ def test_cape_cin_custom_profile(): dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) + 5 * units.delta_degC cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 1650.61208729 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 1641.74927485 * units('joule / kilogram'), 2) assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 2) @@ -2258,7 +2311,7 @@ def test_cape_cin_value_error(): -35.9, -26.7, -37.7, -43.1, -33.9, -40.9, -46.1, -34.9, -33.9, -33.7, -33.3, -42.5, -50.3, -49.7, -49.5, -58.3, -61.3]) * units.degC cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint) - expected_cape, expected_cin = 2098.688061 * units('joules/kg'), 0.0 * units('joules/kg') + expected_cape, expected_cin = 2088.286477 * units('joules/kg'), 0.0 * units('joules/kg') assert_almost_equal(cape, expected_cape, 3) assert_almost_equal(cin, expected_cin, 3) @@ -2269,8 +2322,8 @@ def test_lcl_grid_surface_lcls(): temperature = np.array([15, 14, 13]) * units.degC dewpoint = np.array([15, 10, 13]) * units.degC lcl_pressure, lcl_temperature = lcl(pressure, temperature, dewpoint) - pres_truth = np.array([1000, 932.1719, 1010]) * units.hPa - temp_truth = np.array([15, 9.10424, 13]) * units.degC + pres_truth = np.array([999.87072, 932.08387, 1009.89001]) * units.hPa + temp_truth = np.array([14.98936, 9.09663, 12.99110]) * units.degC assert_array_almost_equal(lcl_pressure, pres_truth, 4) assert_array_almost_equal(lcl_temperature, temp_truth, 4) @@ -2375,7 +2428,7 @@ def test_lifted_index(): -57.5]) * units.degC parcel_prof = parcel_profile(pressure, temperature[0], dewpoint[0]) li = lifted_index(pressure, temperature, parcel_prof) - assert_almost_equal(li, -7.9115691 * units.delta_degree_Celsius, 2) + assert_almost_equal(li, -7.8752751 * units.delta_degree_Celsius, 2) def test_lifted_index_500hpa_missing(): @@ -2488,8 +2541,8 @@ def test_gdi_xarray(index_xarray_data_expanded): assert_array_almost_equal( result, - np.array([[[189.5890429, 157.4307982, 129.9739099], - [106.6763526, 87.0637477, 70.7202505]]]) + np.array([[[188.8392281, 156.8406814, 129.5099932], + [106.3119035, 86.7774799, 70.4952583]]]) ) @@ -2508,8 +2561,8 @@ def test_gdi_arrays(index_xarray_data_expanded): assert_array_almost_equal( result, - np.array([[189.5890429, 157.4307982, 129.9739099], - [106.6763526, 87.0637477, 70.7202505]]) + np.array([[188.8392281, 156.8406814, 129.5099932], + [106.3119035, 86.7774799, 70.4952583]]) ) @@ -2523,7 +2576,7 @@ def test_gdi_profile(index_xarray_data_expanded): pressure, temperature, relative_humidity_from_dewpoint(temperature, dewpoint)) assert_almost_equal(galvez_davison_index(pressure, temperature, mixing_ratio, pressure[0]), - 189.5890429, 4) + 188.8392281, 4) def test_gdi_no_950_raises_valueerror(index_xarray_data): @@ -2617,7 +2670,7 @@ def test_showalter_index(): -48.9, -50.2, -51.5, -53.3, -55.5, -55.9]), 'degC') result = showalter_index(pressure, temps, dewp) - assert_almost_equal(result, units.Quantity(7.6024, 'delta_degC'), 4) + assert_almost_equal(result, units.Quantity(7.6145, 'delta_degC'), 4) def test_total_totals_index(): @@ -2718,11 +2771,11 @@ def test_parcel_profile_drop_duplicates(): dewpoint = units.Quantity(18.6, 'degC') - truth = np.array([292.75, 291.78965331, 291.12778784, 290.61996294, - 289.93681828, 289.84313902, 289.36183185, 288.5626898, - 135.46280886, 134.99220142, 131.27369084, 131.07055878, - 125.93977169, 123.63877507, 120.85291224, 120.85291224, - 119.20448296]) * units.kelvin + truth = np.array([292.75, 291.78965, 291.11868, 290.61062, + 289.92716, 289.83344, 289.35192, 288.55244, + 135.44123, 134.9707, 131.25278, 131.04968, + 125.91971, 123.61908, 120.83366, 120.83366, + 119.18549]) * units.kelvin with pytest.warns(UserWarning, match='Duplicate pressure'): profile = parcel_profile(pressure, temperature, dewpoint) @@ -2747,24 +2800,24 @@ def test_parcel_profile_with_lcl_as_dataset_duplicates(): { 'ambient_temperature': ( ('isobaric',), - np.insert(temperature.m, 2, 19.679237747615478) * units.degC + np.insert(temperature.m, 2, 19.67170422) * units.degC ), 'ambient_dew_point': ( ('isobaric',), - np.insert(dewpoint.m, 2, 19.143390198092384) * units.degC + np.insert(dewpoint.m, 2, 19.13736338) * units.degC ), 'parcel_temperature': ( ('isobaric',), [ - 293.15, 293.15, 292.40749167, 292.22841462, 291.73069653, 291.06139433, - 125.22698955, 122.40534065, 122.40534065, 120.73573642, 119.0063293 + 293.15, 293.15, 292.38999616, 292.2183331, 291.72035678, 291.05071707, + 125.20436659, 122.38322754, 122.38322754, 120.71392503, 118.98483107 ] * units.kelvin ) }, coords={ 'isobaric': ( 'isobaric', - np.insert(pressure.m, 2, 942.6) + np.insert(pressure.m, 2, 942.4) ) } )