Add downdraft CAPE calculation

metpy/calc/tests/test_thermo.py

@@ -9,7 +9,7 @@
 from metpy.calc import (brunt_vaisala_frequency, brunt_vaisala_frequency_squared,
                         brunt_vaisala_period, cape_cin, density, dewpoint,
                         dewpoint_from_specific_humidity, dewpoint_rh,
-                        dry_lapse, dry_static_energy, el,
+                        downdraft_cape, dry_lapse, dry_static_energy, el,
                         equivalent_potential_temperature,
                         exner_function, isentropic_interpolation, lcl, lfc, mixed_layer,
                         mixed_parcel, mixing_ratio, mixing_ratio_from_relative_humidity,
@@ -1020,3 +1020,48 @@ def test_lfc_not_below_lcl():
     # Before patch, LFC pressure would show 1000.5912165339967 hPa
     assert_almost_equal(lfc_pressure, 811.8456357 * units.mbar, 6)
     assert_almost_equal(lfc_temp, 6.4992871 * units.celsius, 6)
+
+
+def test_dcape_defaults():
+    """Test DCAPE with the default behavior."""
+    pressure = np.array([973, 943.5, 925, 910.6, 878.4, 865, 850, 848, 847.2, 816.9,
+                        793, 787.8, 759.6, 759, 732.2, 700, 654.5, 607.1,
+                        607, 580, 562, 549.]) * units.hPa
+    temperature = np.array([23.4, 20.4, 18.4, 17.1, 14.1, 12.8, 13.2, 13.4, 13.4,
+                           13.1, 12.8, 12.7, 12.4, 12.4, 10, 7, 3.4, -0.7,
+                           -0.7, -3.3, -5.5, -7.1]) * units.degC
+    dewpoint = np.array([5.4, 3.6, 2.4, 1.8, 0.4, -0.2, 0.2, 0.4, 0.2, -5.5, -10.2,
+                         -11.9, -21.4, -21.6, -21.8, -22, -21.4, -20.7, -20.7,
+                         -38.3, -27.8, -20.1]) * units.degC
+    dcape, dcape_pressure, dcape_temperature = downdraft_cape(pressure, temperature, dewpoint)
+    dcape_truth = 74.11506371089433 * units.joule / units.kilogram
+    dcape_pressure_truth = np.array([973, 943.5, 925, 910.6, 878.4]) * units.hPa
+    dcape_temperature_truth = np.array([17.95718657, 16.80836487,
+                                        16.06406211, 15.47121721, 14.1]) * units.degC
+    assert_almost_equal(dcape, dcape_truth, 6)
+    assert_almost_equal(dcape_pressure, dcape_pressure_truth, 6)
+    assert_almost_equal(dcape_temperature, dcape_temperature_truth, 6)
+
+
+def test_dcape_custom_parcel_start():
+    """Test DCAPE with a custom parcel starting point."""
+    pressure = np.array([973, 943.5, 925, 910.6, 878.4, 865, 850, 848, 847.2, 816.9,
+                        793, 787.8, 759.6, 759, 732.2, 700, 654.5]) * units.hPa
+    temperature = np.array([23.4, 20.4, 18.4, 17.1, 14.1, 12.8, 13.2, 13.4, 13.4,
+                           13.1, 12.8, 12.7, 12.4, 12.4, 10, 7, 3.4]) * units.degC
+    dewpoint = np.array([5.4, 3.6, 2.4, 1.8, 0.4, -0.2, 0.2, 0.4, 0.2, -5.5, -10.2,
+                         -11.9, -21.4, -21.6, -21.8, -22, -21.4]) * units.degC
+    custom_parcel = (670 * units.hPa, 3.5 * units.degC)
+    dcape, dcape_pressure, dcape_temperature = downdraft_cape(pressure, temperature, dewpoint,
+                                                              parcel=custom_parcel)
+    dcape_truth = 101.56717405359117 * units.joule / units.kilogram
+    dcape_pressure_truth = np.array([973, 943.5, 925, 910.6, 878.4, 865, 850, 848, 847.2,
+                                     816.9, 793, 787.8, 759.6, 759, 732.2, 700, 670]) * units.hPa
+    dcape_temperature_truth = np.array([19.0633538, 17.93772489, 17.20885559, 16.62853588,
+                                        15.2870874, 14.70992377, 14.04982066, 13.96065181,
+                                        13.92490675, 12.53726347, 11.39336396, 11.13832389,
+                                        9.71426456, 9.68318326, 8.25933429, 6.44940391,
+                                        4.65373642]) * units.degC
+    assert_almost_equal(dcape, dcape_truth, 6)
+    assert_almost_equal(dcape_pressure, dcape_pressure_truth, 6)
+    assert_almost_equal(dcape_temperature, dcape_temperature_truth, 6)

metpy/calc/thermo.py

@@ -2105,3 +2105,103 @@ def dewpoint_from_specific_humidity(specific_humidity, temperature, pressure):
     return dewpoint_rh(temperature, relative_humidity_from_specific_humidity(specific_humidity,
                                                                              temperature,
                                                                              pressure))
+
+
+@exporter.export
+@check_units('[pressure]', '[temperature]', '[temperature]')
+def downdraft_cape(pressure, temperature, dewpoint, parcel=None, bottom=None,
+                   heights=None, depth=400 * units.hPa):
+    r"""Calculate downdraft CAPE.
+
+    Calculate the downdraft convective available potential energy (CAPE).The minimum theta-e
+    value between the surface and top_pressure is used as the parcel starting point. Parcels
+    descend along a moist adiabat. Area between the parcel path and environmental temperature
+    is integrated to calculate DCAPE.
+
+    Parameters
+    ----------
+    pressure : `pint.Quantity`
+        The atmospheric pressure level(s) of interest. The first entry should be the starting
+        point pressure.
+    temperature : `pint.Quantity`
+        The atmospheric temperature corresponding to pressure.
+    dewpoint : `pint.Quantity`
+        The atmospheric dew point corresponding to pressure.
+    parcel : tuple
+        Tuple of pressure and temperature for the starting parcel.
+    bottom : `pint.Quantity`
+        Lowest point in the parcel path to consider in integration in height or pressure.
+        If no heights are given, a standard atmosphere will be assumed. Defaults to None.
+    top_pressure : `pint.Quantity`
+        The lowest pressure value to be considered in the calculation. Defaults to 400 hPa.
+
+    Returns
+    -------
+    `pint.Quantity`
+        Downdraft convective available potential energy (CAPE).
+
+    Notes
+    -----
+    .. math:: \text{DCAPE} = R_d \int_{P_i}^{P_{sfc}} (T_{env} - T_{parcel}) d\text{ln}(p)
+
+    * :math:`DCAPE` Downdraft convective available potential energy
+    * :math:`R_d` Gas constant
+    * :math:`g` Gravitational acceleration
+    * :math:`T_{parcel}` Parcel temperature
+    * :math:`T_{env}` Environment temperature
+    * :math:`p` Atmospheric pressure
+
+    """
+
+    # If the user specified a parcel starting point, we'll use that instead of the default
+    if parcel:
+        if depth:
+            ValueError('Cannot specify a depth when using a custom parcel.')
+
+        parcel_starting_pressure, parcel_starting_temperature = parcel
+
+        # If no bottom is specified, we'll use the surface (default of get_layer),but
+        # done explicitly here to depth calculation.
+        if not bottom:
+            bottom = np.nanmax(pressure) * pressure.units
+
+        # Calculate the depth of the sounding to use (bottom - the parcel starting point)
+        depth = bottom - parcel_starting_pressure
+
+        # Trim data to the bottom and depth above it
+        pressure, temperature = get_layer(pressure, temperature, bottom=bottom,
+                                          depth=depth, heights=heights)
+
+    # The user did not give us a parcel, so we'll calculate a sensible default.
+    else:
+
+        # Trim data to only top height and below as well as bottom height and above
+        pressure, temperature, dewpoint = get_layer(pressure, temperature, dewpoint,
+                                                    depth=depth,
+                                                    bottom=bottom,
+                                                    heights=heights)
+
+        # Find minimum theta-e
+        theta_e = equivalent_potential_temperature(pressure, temperature, dewpoint)
+        minimum_theta_e_index = np.argmin(theta_e)
+
+        # Trim data to be minimum theta-e down (dewpoint is no longer required)
+        pressure = pressure[:minimum_theta_e_index]
+        temperature = temperature[:minimum_theta_e_index]
+
+    # Sort for monotonically increasing pressure that trapz needs to work.
+    # Also guarantees the logic for calculating a descending parcel is sound.
+    sort_inds = np.argsort(pressure)
+    pressure = pressure[sort_inds]
+    temperature = temperature[sort_inds]
+
+    # Create the parcel profile for decent along a moist adiabat
+    profile_temperatures = moist_lapse(pressure, temperature[0])
+
+    # Calculate the difference in profile and environmental temperature to integrate
+    y_vals = temperature - profile_temperatures
+
+    # Calculate DCAPE
+    dcape = Rd * (np.trapz(y_vals,
+                           x=np.log(pressure.m)) * units.kelvin)
+    return dcape.to('J/kg'), pressure[::-1], profile_temperatures[::-1]