Skip to content
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

Return input temperature units for LCL/LFC/EL temperature, disallow LFC if EL below LCL #1172

Merged
merged 2 commits into from
Sep 26, 2019
Merged

Return input temperature units for LCL/LFC/EL temperature, disallow LFC if EL below LCL #1172

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
10bbda7
Return input temperature units for LCL/LFC/EL.
zbruick Sep 18, 2019
60ca790
Reduce return statements in `lfc`
zbruick Sep 19, 2019
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
50 changes: 49 additions & 1 deletion metpy/calc/tests/test_thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# Copyright (c) 2008,2015,2016,2017,2018 MetPy Developers.
# Copyright (c) 2008,2015,2016,2017,2018,2019 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Test the `thermo` module."""
Expand Down Expand Up @@ -260,6 +260,16 @@ def test_lcl():
assert_almost_equal(lcl_temperature, 17.676 * units.degC, 2)


def test_lcl_kelvin():
"""Test LCL temperature is returned as Kelvin, if temperature is 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.416 * units.mbar, 2)
assert_almost_equal(lcl_temperature, 262.827 * units.kelvin, 2)
assert lcl_temperature.units == temperature.units


def test_lcl_convergence():
"""Test LCL calculation convergence failure."""
with pytest.raises(RuntimeError):
Expand All @@ -276,6 +286,19 @@ def test_lfc_basic():
assert_almost_equal(lfc_temp, 9.705 * units.celsius, 2)


def test_lfc_kelvin():
"""Test that LFC temperature returns Kelvin if Kelvin is provided."""
pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
temperature = (np.array([22.2, 14.6, 12., 9.4, 7., -49.]
) + 273.15) * units.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.415 * units.mbar, 2)
assert_almost_equal(lfc_temp, 9.705 * units.degC, 2)
assert lfc_temp.units == temperature.units


def test_lfc_ml():
"""Test Mixed-Layer LFC calculation."""
levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
Expand Down Expand Up @@ -499,6 +522,17 @@ def test_el():
assert_almost_equal(el_temperature, -11.7027 * units.degC, 3)


def test_el_kelvin():
"""Test that EL temperature returns Kelvin if Kelvin is provided."""
levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
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, 470.4075 * units.mbar, 3)
assert_almost_equal(el_temp, -11.7027 * units.degC, 3)
assert el_temp.units == temperatures.units


def test_el_ml():
"""Test equilibrium layer calculation for a mixed parcel."""
levels = np.array([959., 779.2, 751.3, 724.3, 700., 400., 269.]) * units.mbar
Expand Down Expand Up @@ -535,6 +569,20 @@ def test_no_el_multi_crossing():
assert_nan(el_temperature, temperatures.units)


def test_lfc_and_el_below_lcl():
"""Test that LFC and EL are returned as NaN if both are below LCL."""
dewpoint = [264.5351, 261.13443, 259.0122, 252.30063, 248.58017, 242.66582] * units.kelvin
temperature = [273.09723, 268.40173, 263.56207, 260.257, 256.63538,
252.91345] * units.kelvin
pressure = [1017.16, 950, 900, 850, 800, 750] * units.hPa
el_pressure, el_temperature = el(pressure, temperature, dewpoint)
lfc_pressure, lfc_temperature = lfc(pressure, temperature, dewpoint)
assert_nan(lfc_pressure, pressure.units)
assert_nan(lfc_temperature, temperature.units)
assert_nan(el_pressure, pressure.units)
assert_nan(el_temperature, temperature.units)


def test_el_lfc_equals_lcl():
"""Test equilibrium layer calculation when the lfc equals the lcl."""
levels = np.array([912., 905.3, 874.4, 850., 815.1, 786.6, 759.1, 748.,
Expand Down
22 changes: 13 additions & 9 deletions metpy/calc/thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# Copyright (c) 2008,2015,2016,2017,2018 MetPy Developers.
# Copyright (c) 2008,2015,2016,2017,2018,2019 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains a collection of thermodynamic calculations."""
Expand Down Expand Up @@ -359,7 +359,7 @@ def _lcl_iter(p, p0, w, t):
fp = so.fixed_point(_lcl_iter, pressure.m, args=(pressure.m, w, temperature),
xtol=eps, maxiter=max_iters)
lcl_p = fp * pressure.units
return lcl_p, dewpoint(vapor_pressure(lcl_p, w))
return lcl_p, dewpoint(vapor_pressure(lcl_p, w)).to(temperature.units)


@exporter.export
Expand Down Expand Up @@ -406,8 +406,7 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta
if parcel_temperature_profile is None:
new_stuff = parcel_profile_with_lcl(pressure, temperature, dewpt)
pressure, temperature, _, parcel_temperature_profile = new_stuff
temperature = temperature.to('degC')
parcel_temperature_profile = parcel_temperature_profile.to('degC')
parcel_temperature_profile = parcel_temperature_profile.to(temperature.units)

if dewpt_start is None:
dewpt_start = dewpt[0]
Expand Down Expand Up @@ -436,17 +435,23 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta
mask = pressure < this_lcl[0]
if np.all(_less_or_close(parcel_temperature_profile[mask], temperature[mask])):
# LFC doesn't exist
return np.nan * pressure.units, np.nan * temperature.units
x, y = np.nan * pressure.units, np.nan * temperature.units
else: # LFC = LCL
x, y = this_lcl
return x, y
return x, y

# LFC exists. Make sure it is no lower than the LCL
else:
idx = x < this_lcl[0]
# LFC height < LCL height, so set LFC = LCL
if not any(idx):
x, y = this_lcl
el_pres, _ = find_intersections(pressure[1:], parcel_temperature_profile[1:],
temperature[1:], direction='decreasing',
log_x=True)
if np.min(el_pres) > this_lcl[0]:
x, y = np.nan * pressure.units, np.nan * temperature.units
else:
x, y = this_lcl
return x, y
# Otherwise, find all LFCs that exist above the LCL
# What is returned depends on which flag as described in the docstring
Expand Down Expand Up @@ -508,8 +513,7 @@ def el(pressure, temperature, dewpt, parcel_temperature_profile=None, which='top
if parcel_temperature_profile is None:
new_stuff = parcel_profile_with_lcl(pressure, temperature, dewpt)
pressure, temperature, _, parcel_temperature_profile = new_stuff
temperature = temperature.to('degC')
parcel_temperature_profile = parcel_temperature_profile.to('degC')
parcel_temperature_profile = parcel_temperature_profile.to(temperature.units)

# If the top of the sounding parcel is warmer than the environment, there is no EL
if parcel_temperature_profile[-1] > temperature[-1]:
Expand Down