Skip to content

Mixed-layer CAPE/CIN accepts depth in units of height #1912

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
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Mixed-layer CAPE/CIN accepts depth in units of height #1912

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
30 changes: 24 additions & 6 deletions src/metpy/calc/thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
import scipy.optimize as so
import xarray as xr

from .basic import height_to_pressure_std
from .tools import (_greater_or_close, _less_or_close, _remove_nans, find_bounding_indices,
find_intersections, first_derivative, get_layer)
from .. import constants as mpconsts
Expand Down Expand Up @@ -2337,8 +2338,8 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
of a given upper air profile and mixed-layer parcel path. CIN is integrated between the
surface and LFC, CAPE is integrated between the LFC and EL (or top of sounding).
Intersection points of the measured temperature profile and parcel profile are
logarithmically interpolated. Kwargs for `mixed_parcel` can be provided, such as `depth`.
Default mixed-layer depth is 100 hPa.
logarithmically interpolated. Kwargs for `mixed_parcel` can be provided, such as `depth`
and `height`. Default mixed-layer depth is 100 hPa.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From the name of the height kwarg it is not clear that this represents an array of height values associated with a vertical profile. Would a variable name of height_profile be more appropriate and descriptive for our users? I think this would also help alleviate some confusion between that kwarg and depth as they are clearly distinct.


Parameters
----------
Expand All @@ -2352,7 +2353,9 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
Dewpoint profile

kwargs
Additional keyword arguments to pass to `mixed_parcel`
Additional keyword arguments to pass to `mixed_parcel`. If `depth` is passed with units
of height without also passing a `height` profile, `depth` will be converted to a
pressure with `height_to_pressure_std`.

Returns
-------
Expand All @@ -2372,14 +2375,29 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
Quantities even when given xarray DataArray profiles.

"""
height = kwargs.get('height')
depth = kwargs.get('depth', units.Quantity(100, 'hPa'))

# Convert depth from a height to a presure if needed
if depth.check('[length]') and height is None:
depth = height_to_pressure_std(units.Quantity(0, 'm')) - height_to_pressure_std(depth)
kwargs['depth'] = depth
warnings.warn('Depth of the mixed layer was given as a height, but no height profile '
'was provided. Depth of mixed layer was converted to '
+ str(round(depth, 2)) + ' using US standard atmosphere.')

parcel_pressure, parcel_temp, parcel_dewpoint = mixed_parcel(pressure, temperature,
dewpoint, **kwargs)

# Remove values below top of mixed layer and add in the mixed layer values
pressure_prof = pressure[pressure < (pressure[0] - depth)]
temp_prof = temperature[pressure < (pressure[0] - depth)]
dew_prof = dewpoint[pressure < (pressure[0] - depth)]
if depth.check('[length]'): # Depth is given as a height
pressure_prof = pressure[height > (depth - height[0])]
temp_prof = temperature[height > (depth - height[0])]
dew_prof = dewpoint[height > (depth - height[0])]
else: # Depth is given as a pressure
pressure_prof = pressure[pressure < (pressure[0] - depth)]
temp_prof = temperature[pressure < (pressure[0] - depth)]
dew_prof = dewpoint[pressure < (pressure[0] - depth)]
pressure_prof = concatenate([parcel_pressure, pressure_prof])
temp_prof = concatenate([parcel_temp, temp_prof])
dew_prof = concatenate([parcel_dewpoint, dew_prof])
Expand Down
21 changes: 21 additions & 0 deletions tests/calc/test_thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1323,6 +1323,27 @@ def test_mixed_layer_cape_cin(multiple_intersections):
assert_almost_equal(mlcin, -20.6727628 * units('joule / kilogram'), 2)


def test_mixed_layer_cape_cin_with_depth_as_height():
"""Test passing depth in units of height to mixed layer cape/cin calculation."""
pressure = np.array([1003., 1000., 972., 925., 901., 851., 850., 727., 708.]) * units.hPa
temperature = np.array([24.8, 24.4, 23., 20.6, 19.3, 16.4, 16.4, 11.4, 11.]) * units.degC
dewpoint = np.array([21.1, 21.5, 21., 20.2, 19., 16.3, 16.3, 8.1, 6.8]) * units.degC
height = np.array([140, 110, 358, 792, 1018, 1510, 1520, 2843, 3065]) * units.meter

# Test passing depth as height without a height profile
with pytest.warns(UserWarning):
mlcape, mlcin = mixed_layer_cape_cin(pressure, temperature, dewpoint,
depth=units.Quantity(500, 'm'))
assert_almost_equal(mlcape, 6.6591 * units('joule / kilogram'), 2)
assert_almost_equal(mlcin, -24.0292 * units('joule / kilogram'), 2)

# Test passing depth as height with a height profile
mlcape, mlcin = mixed_layer_cape_cin(pressure, temperature, dewpoint, height=height,
depth=units.Quantity(500, 'm'))
assert_almost_equal(mlcape, 7.2832 * units('joule / kilogram'), 2)
assert_almost_equal(mlcin, -23.3026 * units('joule / kilogram'), 2)


def test_mixed_layer():
"""Test the mixed layer calculation."""
pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa
Expand Down