Skip to content

Compute the LCL from actual temperature in lfc/el used by cape_cin #4067

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

Open
gaoflow wants to merge 2 commits into
base: main
Choose a base branch
from
+101 −13
Open

Compute the LCL from actual temperature in lfc/el used by cape_cin #4067

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
c0dfd81
Compute LCL from actual temperatures in lfc/el for cape_cin (GH #3857)
gaoflow May 27, 2026
36034a6
Update cape_cin/most_unstable_cape_cin doctests for corrected CAPE
gaoflow May 27, 2026
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
60 changes: 48 additions & 12 deletions src/metpy/calc/thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -847,7 +847,7 @@ def ccl(pressure, temperature, dewpoint, height=None, mixed_layer_depth=None, wh
@preprocess_and_wrap()
@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]')
def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoint_start=None,
which='top'):
which='top', lcl_pressure=None, lcl_temperature=None):
r"""Calculate the level of free convection (LFC).

This works by finding the first intersection of the ideal parcel path and
Expand Down Expand Up @@ -882,6 +882,18 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi
'wide' returns the LFC whose corresponding EL is farthest away,
'most_cape' returns the LFC that results in the most CAPE in the profile.

lcl_pressure: `pint.Quantity`, optional
Pressure of the parcel's lifting condensation level (LCL). If not provided, it is
computed from ``parcel_temperature_profile`` and ``dewpoint_start``. Supplying it is
useful when ``temperature`` and ``parcel_temperature_profile`` are virtual
temperatures (as in :func:`cape_cin`), so that the LCL is computed from the actual
temperatures rather than the virtual ones. Must be given together with
``lcl_temperature``.

lcl_temperature: `pint.Quantity`, optional
Temperature of the parcel's lifting condensation level (LCL), used alongside
``lcl_pressure`` when the LFC coincides with the LCL.

Returns
-------
`pint.Quantity`
Expand Down Expand Up @@ -951,8 +963,13 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi
x, y = find_intersections(pressure, parcel_temperature_profile,
temperature, direction='increasing', log_x=True)

# Compute LCL for this parcel for future comparisons
this_lcl = lcl(pressure[0], parcel_temperature_profile[0], dewpoint_start)
# Compute LCL for this parcel for future comparisons. Allow the caller to supply a
# precomputed LCL so it can be based on the actual (not virtual) temperatures when the
# temperature profiles passed in are virtual temperatures (see cape_cin).
if lcl_pressure is None:
this_lcl = lcl(pressure[0], parcel_temperature_profile[0], dewpoint_start)
else:
this_lcl = (lcl_pressure, lcl_temperature)

# The LFC could:
# 1) Not exist
Expand Down Expand Up @@ -1067,7 +1084,8 @@ def _most_cape_option(intersect_type, p_list, t_list, pressure, temperature, dew
@exporter.export
@preprocess_and_wrap()
@check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]')
def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top'):
def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top',
lcl_pressure=None):
r"""Calculate the equilibrium level.

This works by finding the last intersection of the ideal parcel path and
Expand Down Expand Up @@ -1096,6 +1114,13 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
'wide' returns the EL whose corresponding LFC is farthest away.
'most_cape' returns the EL that results in the most CAPE in the profile.

lcl_pressure: `pint.Quantity`, optional
Pressure of the parcel's lifting condensation level (LCL). If not provided, it is
computed from ``temperature`` and ``dewpoint``. Supplying it is useful when
``temperature`` and ``parcel_temperature_profile`` are virtual temperatures (as in
:func:`cape_cin`), so that the LCL is computed from the actual temperatures rather
than the virtual ones.

Returns
-------
`pint.Quantity`
Expand Down Expand Up @@ -1163,7 +1188,12 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
# and reassigned to allow np.log() to function properly.
x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:],
direction='decreasing', log_x=True)
lcl_p, _ = lcl(pressure[0], temperature[0], dewpoint[0])
# Allow the caller to supply a precomputed LCL pressure so it can be based on the actual
# (not virtual) temperatures when the temperature profiles passed in are virtual
# temperatures (see cape_cin).
if lcl_pressure is None:
lcl_pressure, _ = lcl(pressure[0], temperature[0], dewpoint[0])
lcl_p = lcl_pressure
if len(x) > 0 and x[-1] < lcl_p:
idx = x < lcl_p
return _multiple_el_lfc_options(x, y, idx, which, pressure,
Expand Down Expand Up @@ -2786,7 +2816,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
>>> prof = parcel_profile(p, T[0], Td[0]).to('degC')
>>> # calculate surface based CAPE/CIN
>>> cape_cin(p, T, Td, prof)
(<Quantity(4830.74608, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
(<Quantity(4902.33809, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)

See Also
--------
Expand Down Expand Up @@ -2825,7 +2855,10 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature,
dewpoint, parcel_profile)

pressure_lcl, _ = lcl(pressure[0], temperature[0], dewpoint[0])
# Compute the LCL from the actual temperatures, before they are replaced below by
# virtual temperatures. This is passed to lfc/el so they do not recompute the LCL from
# the virtual temperatures (see GH #3857).
pressure_lcl, temperature_lcl = lcl(pressure[0], temperature[0], dewpoint[0])
below_lcl = pressure > pressure_lcl

# The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated
Expand All @@ -2840,19 +2873,22 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint)
parcel_profile = virtual_temperature(parcel_profile, parcel_mixing_ratio)

# Calculate LFC limit of integration
# Calculate LFC limit of integration. The LCL is supplied from the actual-temperature
# computation above so it is not recomputed from the virtual temperatures.
lfc_pressure, _ = lfc(pressure, temperature, dewpoint,
parcel_temperature_profile=parcel_profile, which=which_lfc)
parcel_temperature_profile=parcel_profile, which=which_lfc,
lcl_pressure=pressure_lcl, lcl_temperature=temperature_lcl)

# If there is no LFC, no need to proceed.
if np.isnan(lfc_pressure):
return units.Quantity(0, 'J/kg'), units.Quantity(0, 'J/kg')
else:
lfc_pressure = lfc_pressure.magnitude

# Calculate the EL limit of integration
# Calculate the EL limit of integration (LCL supplied as above).
el_pressure, _ = el(pressure, temperature, dewpoint,
parcel_temperature_profile=parcel_profile, which=which_el)
parcel_temperature_profile=parcel_profile, which=which_el,
lcl_pressure=pressure_lcl)

# No EL and we use the top reading of the sounding.
el_pressure = pressure[-1].magnitude if np.isnan(el_pressure) else el_pressure.magnitude
Expand Down Expand Up @@ -3427,7 +3463,7 @@ def most_unstable_cape_cin(pressure, temperature, dewpoint, **kwargs):
>>> Td = dewpoint_from_relative_humidity(T, rh)
>>> # calculate most unstbale CAPE/CIN
>>> most_unstable_cape_cin(p, T, Td)
(<Quantity(4830.74608, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
(<Quantity(4921.07107, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)

See Also
--------
Expand Down
54 changes: 53 additions & 1 deletion tests/calc/test_thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2379,11 +2379,63 @@ 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 = 2161.912443 * units('joules/kg'), 0.0 * units('joules/kg')
# The LFC for this profile is floored at the LCL; with the LCL now computed from the
# actual (not virtual) temperatures (GH #3857) the LFC moves to the correct, higher
# pressure, widening the integration range and increasing CAPE from the prior 2161.912.
expected_cape, expected_cin = 2206.730694 * units('joules/kg'), 0.0 * units('joules/kg')
assert_almost_equal(cape, expected_cape, 3)
assert_almost_equal(cin, expected_cin, 3)


def test_lfc_el_lcl_from_actual_temperature():
"""Test that lfc and el base the LCL on actual, not virtual, temperatures (GH #3857).

``cape_cin`` replaces the temperature profiles with virtual temperatures before calling
``lfc``/``el``. Since virtual temperature exceeds temperature, recomputing the LCL from
those profiles puts it at too low a pressure (too high an altitude), which in turn
biases the LFC and EL upward. The ``lcl_pressure``/``lcl_temperature`` keywords let the
caller supply the LCL computed from the actual temperatures instead.
"""
pressure = np.array([1000., 925., 850., 800., 700., 600., 500., 400., 300.]) * units.hPa
temperature = np.array([28., 22., 16., 13., 6., -2., -12., -25., -40.]) * units.degC
dewpoint = np.array([24., 19., 13., 9., 0., -8., -20., -35., -55.]) * units.degC
profile = parcel_profile(pressure, temperature[0], dewpoint[0])

# Reproduce the virtual-temperature profiles cape_cin builds internally.
lcl_pressure, lcl_temperature = lcl(pressure[0], temperature[0], dewpoint[0])
below_lcl = pressure > lcl_pressure
parcel_mixing_ratio = np.where(below_lcl,
saturation_mixing_ratio(pressure[0], dewpoint[0]),
saturation_mixing_ratio(pressure, profile))
temperature_v = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint)
profile_v = virtual_temperature(profile, parcel_mixing_ratio)

# The LCL recomputed from the virtual parcel temperature is spuriously low (high).
lcl_pressure_virtual, _ = lcl(pressure[0], profile_v[0], dewpoint[0])
assert lcl_pressure > lcl_pressure_virtual

# The LFC for this profile is floored at the LCL, so it tracks whichever LCL is used:
# without the override it floors at the wrong (virtual) LCL ...
lfc_pressure_virtual, _ = lfc(pressure, temperature_v, dewpoint,
parcel_temperature_profile=profile_v, which='bottom')
assert_almost_equal(lfc_pressure_virtual, lcl_pressure_virtual, 5)

# ... and with the actual-temperature LCL supplied it floors at the correct LCL.
lfc_pressure_actual, lfc_temperature_actual = lfc(
pressure, temperature_v, dewpoint, parcel_temperature_profile=profile_v,
which='bottom', lcl_pressure=lcl_pressure, lcl_temperature=lcl_temperature)
assert_almost_equal(lfc_pressure_actual, lcl_pressure, 5)
assert_almost_equal(lfc_temperature_actual, lcl_temperature, 5)
assert lfc_pressure_actual > lfc_pressure_virtual

# el likewise honors the supplied LCL when deciding which crossings lie above it.
el_default, _ = el(pressure, temperature_v, dewpoint,
parcel_temperature_profile=profile_v)
el_with_lcl, _ = el(pressure, temperature_v, dewpoint,
parcel_temperature_profile=profile_v, lcl_pressure=lcl_pressure)
assert_almost_equal(el_default, el_with_lcl, 5)


def test_lcl_grid_surface_lcls():
"""Test surface grid where some values have LCLs at the surface."""
pressure = np.array([1000, 990, 1010]) * units.hPa
Expand Down
Loading