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

mixed_layer_cape_cin() can handle `depth` in units of height passed through kwargs. This works regardless of whether a `height` profile is provided.

Author: 23ccozad

src/metpy/calc/thermo.py

@@ -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
@@ -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.
 
     Parameters
     ----------
@@ -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
     -------
@@ -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])

tests/calc/test_thermo.py

@@ -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