From eaa93e73c7d77901884fb81c23107d715ff9ba1d Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Tue, 14 Nov 2023 17:31:28 -0600 Subject: [PATCH] add downward cape (#3120) * add downward cape --- docs/_templates/overrides/metpy.calc.rst | 1 + docs/api/references.rst | 2 + src/metpy/calc/thermo.py | 127 +++++++++++++++++++++++ tests/calc/test_thermo.py | 26 ++++- 4 files changed, 154 insertions(+), 2 deletions(-) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index c5920a81505..fc2836cfdd1 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -75,6 +75,7 @@ Soundings ccl critical_angle cross_totals + downdraft_cape el k_index lcl diff --git a/docs/api/references.rst b/docs/api/references.rst index c99e3d7d753..5a3796a5cbe 100644 --- a/docs/api/references.rst +++ b/docs/api/references.rst @@ -58,6 +58,8 @@ References Parameters in Forecasting Severe Storms `_. *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 `_. *Electronic J. Severe Storms Meteor.*, **3** (2), 1-50. diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 25357c2fa38..fa172978020 100644 --- a/src/metpy/calc/thermo.py +++ b/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) + (, , ) + + 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]') diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 05144fd408d..d53da112bfb 100644 --- a/tests/calc/test_thermo.py +++ b/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