diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py index 1ae7e0fd828..2b79fbc85b7 100644 --- a/metpy/calc/tests/test_thermo.py +++ b/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) diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py index f9b217a25ce..aa2886b8885 100644 --- a/metpy/calc/thermo.py +++ b/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]