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.

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

Fix poor assumptions in lfc() logic #1023

Merged
merged 3 commits into from Apr 10, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
21 changes: 17 additions & 4 deletions metpy/calc/tests/test_thermo.py
Expand Up @@ -310,7 +310,7 @@ def test_lfc_ml2():
-71.53523254, -71.61097717, -71.92687988, -72.68682861,
-74.129776, -76.02471924, -76.88977051, -76.26008606,
-75.90351868, -76.15809631]) * units.celsius
dewpoints = np.array([4.50012302, 3.42483997, 2.78102994, 2.24474645, 1.593485, 0.94408154,
dewpoints = np.array([4.50012302, 3.42483997, 2.78102994, 2.24474645, 1.593485, -0.9440815,
-3.8044982, -3.55629468, -9.7376976, -10.2950449, -9.67498302,
-10.30486488, -8.70559597, -8.71669006, -12.66509628, -18.6697197,
-23.00351334, -29.46240425, -36.82178497, -41.68824768, -44.50320816,
Expand All @@ -322,9 +322,22 @@ def test_lfc_ml2():
-88.74453735, -89.04680634, -89.26436615]) * units.celsius
__, 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, 1001.532 * units.mbar, 2)
assert_almost_equal(lfc_temp, 4.17 * units.degC, 2)
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)


def test_lfc_intersection():
"""Test LFC calculation when LFC is below a tricky intersection."""
p = np.array([1024.957, 930., 924.828, 898.255, 873.461, 848.698, 823.926,
788.493]) * units('hPa')
t = np.array([6.008, -10., -6.94, -8.58, -4.41, -4.19, -3.71, -4.48]) * units('degC')
td = np.array([5., -10., -7., -9., -4.5, -4.2, -3.8, -4.5]) * units('degC')
_, 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, 982.762 * units.hPa, 2)
assert_almost_equal(mllfc_t, 272.045 * units.kelvin, 2)


def test_no_lfc():
Expand Down
23 changes: 17 additions & 6 deletions metpy/calc/thermo.py
Expand Up @@ -362,7 +362,7 @@ def _lcl_iter(p, p0, w, t):
@exporter.export
@preprocess_xarray
@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]')
def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_start=None):
r"""Calculate the level of free convection (LFC).

This works by finding the first intersection of the ideal parcel path and
Expand All @@ -379,6 +379,9 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
parcel_temperature_profile: `pint.Quantity`, optional
The parcel temperature profile from which to calculate the LFC. Defaults to the
surface parcel profile.
dewpt_start: `pint.Quantity`, optional
The dewpoint of the parcel for which to calculate the LFC. Defaults to the surface
dewpoint.

Returns
-------
Expand All @@ -397,13 +400,21 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
temperature = temperature.to('degC')
parcel_temperature_profile = parcel_temperature_profile.to('degC')

# The parcel profile and data have the same first data point, so we ignore
# that point to get the real first intersection for the LFC calculation.
x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:],
temperature[1:], direction='increasing')
if dewpt_start is None:
dewpt_start = dewpt[0]

# The parcel profile and data may have the same first data point.
# If that is the case, ignore that point to get the real first
# intersection for the LFC calculation.
if np.isclose(parcel_temperature_profile[0].m, temperature[0].m):
x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:],
temperature[1:], direction='increasing')
else:
x, y = find_intersections(pressure, parcel_temperature_profile,
temperature, direction='increasing')

# Compute LCL for this parcel for future comparisons
this_lcl = lcl(pressure[0], temperature[0], dewpt[0])
this_lcl = lcl(pressure[0], parcel_temperature_profile[0], dewpt_start)

# The LFC could:
# 1) Not exist
Expand Down