From 9d95c5b6a4d36a06ee2b36f6689b665b7e1dd722 Mon Sep 17 00:00:00 2001 From: Connor Cozad <23ccozad@gmail.com> Date: Tue, 15 Jun 2021 13:20:05 -0400 Subject: [PATCH] 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. --- src/metpy/calc/thermo.py | 30 ++++++++++++++++++++++++------ tests/calc/test_thermo.py | 22 ++++++++++++++++++++++ 2 files changed, 46 insertions(+), 6 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index f69c4d1efa9..153e466dcc7 100644 --- a/src/metpy/calc/thermo.py +++ b/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]) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 1b71c169be9..ad8d1d6b9b5 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -1323,6 +1323,28 @@ 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