From b5fcf6a4a02edc332fa1a8f484d875507df43a50 Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Sun, 16 Jul 2023 13:08:15 -0500 Subject: [PATCH 1/8] add downward cape --- docs/_templates/overrides/metpy.calc.rst | 1 + docs/api/references.rst | 2 + src/metpy/calc/thermo.py | 126 +++++++++++++++++++++++ tests/calc/test_thermo.py | 26 ++++- 4 files changed, 153 insertions(+), 2 deletions(-) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index c5920a81505..fe5c8effd1a 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -93,6 +93,7 @@ Soundings storm_relative_helicity supercell_composite surface_based_cape_cin + down_cape sweat_index total_totals_index vertical_totals diff --git a/docs/api/references.rst b/docs/api/references.rst index e4df24baa3e..8b58f6cc770 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 327aff534c1..7e6cb742402 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -3028,6 +3028,132 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs): return cape_cin(p, t, td, ml_profile) +@exporter.export +@preprocess_and_wrap() +def down_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 bouyancy 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, down_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) + >>> down_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.Pa) + 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]) + + diff = (env_virt_temp - parcel_virt_temp).to(units.degK) + + # Find DCAPE + dcape = -(mpconsts.Rd + * units.Quantity( + np.trapz(diff.magnitude, np.log(down_pressure.magnitude)), '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 3119ed1da40..f94da633add 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -15,8 +15,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, + down_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, @@ -1551,6 +1552,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 = down_cape(pressure, temperature, dewpoint) + assert_almost_equal(dcape, 1222 * units('joule / kilogram'), 0) + assert np.all(down_press == pressure[:10]) + 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 From 1870a40d0b118c210605e8d4c1bccc043e6237b9 Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Sun, 16 Jul 2023 13:14:28 -0500 Subject: [PATCH 2/8] that's a hard word --- src/metpy/calc/thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 7e6cb742402..66bb64119a1 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -3034,7 +3034,7 @@ def down_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 bouyancy energy available to a descending + 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 From f93a481adfe6ce7199cd239317ae0c33d35fdb2b Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Sun, 16 Jul 2023 13:36:52 -0500 Subject: [PATCH 3/8] use hPa for pressure --- src/metpy/calc/thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 66bb64119a1..d97fbcb1f96 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -3132,7 +3132,7 @@ def down_cape(pressure, temperature, dewpoint): 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.Pa) + 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) From 06461273e9e919031e5a3adcbaf1111344837b3c Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Sun, 16 Jul 2023 14:11:09 -0500 Subject: [PATCH 4/8] fix spacing --- src/metpy/calc/thermo.py | 13 +++++++------ tests/calc/test_thermo.py | 20 ++++++++++---------- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index d97fbcb1f96..ed6cb2c1ddf 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -3120,8 +3120,8 @@ def down_cape(pressure, temperature, dewpoint): # 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) + 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 @@ -3143,13 +3143,14 @@ def down_cape(pressure, temperature, dewpoint): temperature[pressure >= parcel_start_p], dewpoint[pressure >= parcel_start_p]) - diff = (env_virt_temp - parcel_virt_temp).to(units.degK) + # 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.magnitude, np.log(down_pressure.magnitude)), 'K') - ).to(units('J/kg')) + * units.Quantity(np.trapz(diff, lnp), 'K') + ).to(units('J/kg')) return dcape, down_pressure, down_parcel_trace diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index f94da633add..bab4105069b 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -1555,22 +1555,22 @@ def test_mixed_layer_cape_cin(multiple_intersections): 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 + 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, + -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 = down_cape(pressure, temperature, dewpoint) assert_almost_equal(dcape, 1222 * units('joule / kilogram'), 0) assert np.all(down_press == pressure[:10]) - 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) + 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(): From b9c33009c916d4c5a14fc1e5a584f98363d0afbd Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Sun, 16 Jul 2023 17:31:55 -0500 Subject: [PATCH 5/8] check for almost instead of exact --- tests/calc/test_thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index bab4105069b..97c374e1520 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -1568,7 +1568,7 @@ def test_dcape(): -75.8, -51.2, -56.4] * units.degC dcape, down_press, down_t = down_cape(pressure, temperature, dewpoint) assert_almost_equal(dcape, 1222 * units('joule / kilogram'), 0) - assert np.all(down_press == pressure[:10]) + 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) From d701b3ef3524e749eb012b6b6bd7f89a84addb02 Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Wed, 9 Aug 2023 21:35:36 +0000 Subject: [PATCH 6/8] I'm just moving whitespace around at this point --- src/metpy/calc/thermo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index ed6cb2c1ddf..ac8bf5e7e74 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -3082,9 +3082,9 @@ def down_cape(pressure, temperature, dewpoint): >>> Td = dewpoint_from_relative_humidity(T, rh) >>> down_cape(p, T, Td) (, , ) + 900. 850. 800. 750. 700. 650. 600.], 'hectopascal')>, ) See Also -------- From 56a9e4c6814ca1ce13c50da4490630bf5a094e9a Mon Sep 17 00:00:00 2001 From: Sam Gardner Date: Thu, 19 Oct 2023 15:36:50 -0500 Subject: [PATCH 7/8] rename to downdraft_cape --- docs/_templates/overrides/metpy.calc.rst | 2 +- src/metpy/calc/thermo.py | 6 +++--- tests/calc/test_thermo.py | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index fe5c8effd1a..11d1310c745 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -93,7 +93,7 @@ Soundings storm_relative_helicity supercell_composite surface_based_cape_cin - down_cape + downdraft_cape sweat_index total_totals_index vertical_totals diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index ac8bf5e7e74..6f1867c8520 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -3030,7 +3030,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs): @exporter.export @preprocess_and_wrap() -def down_cape(pressure, temperature, dewpoint): +def downdraft_cape(pressure, temperature, dewpoint): r"""Calculate downward CAPE (DCAPE). Calculate the downward convective available potential energy (DCAPE) of a given upper air @@ -3062,7 +3062,7 @@ def down_cape(pressure, temperature, dewpoint): Examples -------- - >>> from metpy.calc import dewpoint_from_relative_humidity, down_cape + >>> 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., @@ -3080,7 +3080,7 @@ def down_cape(pressure, temperature, dewpoint): ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless >>> # calculate dewpoint >>> Td = dewpoint_from_relative_humidity(T, rh) - >>> down_cape(p, T, Td) + >>> downdraft_cape(p, T, Td) (, , Date: Tue, 14 Nov 2023 16:26:10 -0600 Subject: [PATCH 8/8] alphabetize --- docs/_templates/overrides/metpy.calc.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index 11d1310c745..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 @@ -93,7 +94,6 @@ Soundings storm_relative_helicity supercell_composite surface_based_cape_cin - downdraft_cape sweat_index total_totals_index vertical_totals