add downward cape (#3120)

* add downward cape

docs/_templates/overrides/metpy.calc.rst

@@ -75,6 +75,7 @@ Soundings
       ccl
       critical_angle
       cross_totals
+      downdraft_cape
       el
       k_index
       lcl

docs/api/references.rst

@@ -58,6 +58,8 @@ References
            Parameters in Forecasting Severe Storms <https://ejssm.org/archives/2006/vol-1-3-2006/>`_.
            *Electronic J. Severe Storms Meteor.*, **1** (3), 1-22.
 
+.. [Emanuel1994] Emanuel, K. A., 1994: Atmospheric Convection. Oxford University Press, 592 pp.
+
 .. [Esterheld2008] Esterheld, J. M. and D. J. Giuliano, 2008: `Discriminating between Tornadic and
            Non-Tornadic Supercells: A New Hodograph Technique <https://ejssm.org/archives/2008/vol-3-2-2008/>`_.
            *Electronic J. Severe Storms Meteor.*, **3** (2), 1-50.

src/metpy/calc/thermo.py

@@ -3042,6 +3042,133 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
     return cape_cin(p, t, td, ml_profile)
 
 
+@exporter.export
+@preprocess_and_wrap()
+def downdraft_cape(pressure, temperature, dewpoint):
+    r"""Calculate downward CAPE (DCAPE).
+
+    Calculate the downward convective available potential energy (DCAPE) of a given upper air
+    profile. Downward CAPE is the maximum negative buoyancy energy available to a descending
+    parcel. Parcel descent is assumed to begin from the lowest equivalent potential temperature
+    between 700 and 500 hPa. This parcel is lowered moist adiabatically from the environmental
+    wet bulb temperature to the surface.  This assumes the parcel remains saturated
+    throughout the descent.
+
+    Parameters
+    ----------
+    pressure : `pint.Quantity`
+        Pressure profile
+
+    temperature : `pint.Quantity`
+        Temperature profile
+
+    dewpoint : `pint.Quantity`
+        Dewpoint profile
+
+    Returns
+    -------
+    dcape: `pint.Quantity`
+        Downward Convective Available Potential Energy (DCAPE)
+    down_pressure: `pint.Quantity`
+        Pressure levels of the descending parcel
+    down_parcel_trace: `pint.Quantity`
+        Temperatures of the descending parcel
+
+    Examples
+    --------
+    >>> from metpy.calc import dewpoint_from_relative_humidity, downdraft_cape
+    >>> from metpy.units import units
+    >>> # pressure
+    >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600.,
+    ...      550., 500., 450., 400., 350., 300., 250., 200.,
+    ...      175., 150., 125., 100., 80., 70., 60., 50.,
+    ...      40., 30., 25., 20.] * units.hPa
+    >>> # temperature
+    >>> T = [29.3, 28.1, 25.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1,
+    ...      -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4,
+    ...      -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3,
+    ...      -56.3, -51.7, -50.7, -47.5] * units.degC
+    >>> # relative humidity
+    >>> rh = [.85, .75, .56, .39, .82, .72, .75, .86, .65, .22, .52,
+    ...       .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88,
+    ...       .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless
+    >>> # calculate dewpoint
+    >>> Td = dewpoint_from_relative_humidity(T, rh)
+    >>> downdraft_cape(p, T, Td)
+    (<Quantity(1222.67968, 'joule / kilogram')>, <Quantity([1008. 1000.  950.
+    900.  850.  800.  750.  700.  650.  600.], 'hectopascal')>, <Quantity([17.50959548
+    17.20643425 15.237249 13.12607097 10.85045704 8.38243809 5.68671014 2.71808363
+    -0.58203825 -4.29053485], 'degree_Celsius')>)
+
+    See Also
+    --------
+    cape_cin, surface_based_cape_cin, most_unstable_cape_cin, mixed_layer_cape_cin
+
+    Notes
+    -----
+    Formula adopted from [Emanuel1994]_.
+
+    .. math:: \text{DCAPE} = -R_d \int_{SFC}^{p_\text{top}}
+            (T_{{v}_{env}} - T_{{v}_{parcel}}) d\text{ln}(p)
+
+
+    * :math:`DCAPE` is downward convective available potential energy
+    * :math:`SFC` is the level of the surface or beginning of parcel path
+    * :math:`p_\text{top}` is pressure of the start of descent path
+    * :math:`R_d` is the gas constant
+    * :math:`T_{{v}_{env}}` is environment virtual temperature
+    * :math:`T_{{v}_{parcel}}` is the parcel virtual temperature
+    * :math:`p` is atmospheric pressure
+
+    Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
+    Since this function returns scalar values when given a profile, this will return Pint
+    Quantities even when given xarray DataArray profiles.
+
+
+
+    """
+    pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint)
+    if not len(pressure) == len(temperature) == len(dewpoint):
+        raise ValueError('Provided pressure, temperature,'
+                         'and dewpoint must be the same length')
+
+    # Get layer between 500 and 700 hPa
+    p_layer, t_layer, td_layer = get_layer(pressure, temperature, dewpoint,
+                                           bottom=700 * units.hPa,
+                                           depth=200 * units.hPa, interpolate=True)
+    theta_e = equivalent_potential_temperature(p_layer, t_layer, td_layer)
+
+    # Find parcel with minimum thetae in the layer
+    min_idx = np.argmin(theta_e)
+    parcel_start_p = p_layer[min_idx]
+
+    parcel_start_td = td_layer[min_idx]
+    parcel_start_wb = wet_bulb_temperature(parcel_start_p, t_layer[min_idx], parcel_start_td)
+
+    # Descend parcel moist adiabatically to surface
+    down_pressure = pressure[pressure >= parcel_start_p].to(units.hPa)
+    down_parcel_trace = moist_lapse(down_pressure, parcel_start_wb,
+                                    reference_pressure=parcel_start_p)
+
+    # Find virtual temperature of parcel and environment
+    parcel_virt_temp = virtual_temperature_from_dewpoint(down_pressure, down_parcel_trace,
+                                                         down_parcel_trace)
+    env_virt_temp = virtual_temperature_from_dewpoint(down_pressure,
+                                                      temperature[pressure >= parcel_start_p],
+                                                      dewpoint[pressure >= parcel_start_p])
+
+    # calculate differences (remove units for NumPy)
+    diff = (env_virt_temp - parcel_virt_temp).to(units.degK).magnitude
+    lnp = np.log(down_pressure.magnitude)
+
+    # Find DCAPE
+    dcape = -(mpconsts.Rd
+              * units.Quantity(np.trapz(diff, lnp), 'K')
+              ).to(units('J/kg'))
+
+    return dcape, down_pressure, down_parcel_trace
+
+
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]', '[temperature]', '[temperature]')

tests/calc/test_thermo.py

@@ -14,8 +14,9 @@
 from metpy.calc import (brunt_vaisala_frequency, brunt_vaisala_frequency_squared,
                         brunt_vaisala_period, cape_cin, ccl, cross_totals, density, dewpoint,
                         dewpoint_from_relative_humidity, dewpoint_from_specific_humidity,
-                        dry_lapse, dry_static_energy, el, equivalent_potential_temperature,
-                        exner_function, gradient_richardson_number, InvalidSoundingError,
+                        downdraft_cape, dry_lapse, dry_static_energy, el,
+                        equivalent_potential_temperature, exner_function,
+                        gradient_richardson_number, InvalidSoundingError,
                         isentropic_interpolation, isentropic_interpolation_as_dataset, k_index,
                         lcl, lfc, lifted_index, mixed_layer, mixed_layer_cape_cin,
                         mixed_parcel, mixing_ratio, mixing_ratio_from_relative_humidity,
@@ -1555,6 +1556,27 @@ def test_mixed_layer_cape_cin(multiple_intersections):
     assert_almost_equal(mlcin, -13.4809966289 * units('joule / kilogram'), 2)
 
 
+def test_dcape():
+    """Test the calculation of DCAPE."""
+    pressure = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600.,
+                550., 500., 450., 400., 350., 300., 250., 200.,
+                175., 150., 125., 100., 80., 70., 60., 50.,
+                40., 30., 25., 20.] * units.hPa
+    temperature = [29.3, 28.1, 25.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1,
+                   -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4,
+                   -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3,
+                   -56.3, -51.7, -50.7, -47.5] * units.degC
+    dewpoint = [26.5, 23.3, 16.1, 6.4, 15.3, 10.9, 8.8, 7.9, 0.6,
+                -16.6, -9.2, -9.9, -14.6, -32.8, -51.2, -32.7, -42.6, -58.9,
+                -69.5, -71.7, -75.9, -79.3, -79.7, -72.5, -73.3, -64.3, -70.6,
+                -75.8, -51.2, -56.4] * units.degC
+    dcape, down_press, down_t = downdraft_cape(pressure, temperature, dewpoint)
+    assert_almost_equal(dcape, 1222 * units('joule / kilogram'), 0)
+    assert_array_almost_equal(down_press, pressure[:10], 0)
+    assert_almost_equal(down_t, [17.5, 17.2, 15.2, 13.1, 10.9, 8.4,
+                                 5.7, 2.7, -0.6, -4.3] * units.degC, 1)
+
+
 def test_mixed_layer():
     """Test the mixed layer calculation."""
     pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa