From 0f673d14208229457e765cd856b3dbf4cd55ca64 Mon Sep 17 00:00:00 2001 From: sec147 Date: Sun, 19 Jun 2022 23:43:28 +0100 Subject: [PATCH 1/3] add wet bulb potential temperature function, based on existing equivalent potential temperature function and formula given in DaviesJones2009 --- docs/_templates/overrides/metpy.calc.rst | 1 + docs/api/references.rst | 4 +++ src/metpy/calc/thermo.py | 39 ++++++++++++++++++++++++ tests/calc/test_thermo.py | 11 ++++++- 4 files changed, 54 insertions(+), 1 deletion(-) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index fc2836cfdd1..e24a5f559c5 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -61,6 +61,7 @@ Moist Thermodynamics virtual_temperature virtual_temperature_from_dewpoint wet_bulb_temperature + wetbulb_potential_temperature Soundings diff --git a/docs/api/references.rst b/docs/api/references.rst index 5a3796a5cbe..3b4215c13c8 100644 --- a/docs/api/references.rst +++ b/docs/api/references.rst @@ -54,6 +54,10 @@ References *Mon. Wea. Rev.*, **137**, 3137-3148, doi: `10.1175/2009mwr2774.1 `_. +.. [DaviesJones2008] Davies-Jones, R., 2008: An Efficient and Accurate Method for Computing the Wet-Bulb Temperature along Pseudoadiabats. + *Mon. Wea. Rev.*, **136**, 2764-2785, + doi: `10.1175/2007mwr2224.1 `_. + .. [DoswellSchultz2006] Doswell, C. A. III, and D. M. Schultz, 2006: `On the Use of Indices and Parameters in Forecasting Severe Storms `_. *Electronic J. Severe Storms Meteor.*, **1** (3), 1-22. diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index fa172978020..94039aa6c6c 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1644,6 +1644,45 @@ def saturation_equivalent_potential_temperature(pressure, temperature): return units.Quantity(th_es, units.kelvin) +@exporter.export +@preprocess_and_wrap( + wrap_like='temperature', + broadcast=('pressure', 'temperature', 'dewpoint') +) +@check_units('[pressure]', '[temperature]', '[temperature]') +def wetbulb_potential_temperature(pressure, temperature, dewpoint): + r"""Calculate wet-bulb potential temperature. + + This calculation must be given an air parcel's pressure, temperature, and dewpoint. + The implementation uses the formula outlined in [DaviesJones2008]_: + + First, theta-e is calculated + then use the formula from [DaviesJones2008]_ + + Parameters + ---------- + pressure: `pint.Quantity` + Total atmospheric pressure + + temperature: `pint.Quantity` + Temperature of parcel + + dewpoint: `pint.Quantity` + Dewpoint of parcel + + Returns + ------- + `pint.Quantity` + wet-bulb potential temperature of the parcel + + """ + theta_e = equivalent_potential_temperature(pressure, temperature, dewpoint).magnitude + x = theta_e / 273.15 + a = 7.101574 - 20.68208 * x + 16.11182 * x ** 2 + 2.574631 * x ** 3 - 5.205688 * x ** 4 + b = 1 - 3.552497 * x + 3.781782 * x ** 2 - 0.6899655 * x ** 3 - 0.5929340 * x ** 4 + return units.Quantity(theta_e - np.exp((a) / (b)), 'kelvin') + + @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'mixing_ratio')) @process_units( diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index d53da112bfb..8af5c439a8d 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -36,7 +36,7 @@ vapor_pressure, vertical_totals, vertical_velocity, vertical_velocity_pressure, virtual_potential_temperature, virtual_temperature, virtual_temperature_from_dewpoint, - wet_bulb_temperature) + wet_bulb_temperature, wetbulb_potential_temperature) from metpy.calc.thermo import _find_append_zero_crossings from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_nan, version_check) @@ -802,6 +802,15 @@ def test_equivalent_potential_temperature_masked(): assert_array_almost_equal(ept, expected, 3) +def test_wetbulb_potential_temperature(): + """Test wetbulb potential temperature calculation.""" + p = 1000 * units.mbar + t = 293. * units.kelvin + td = 291. * units.kelvin + wpt = wetbulb_potential_temperature(p, t, td) + assert_almost_equal(wpt, 291.65839705486565 * units.kelvin, 3) + + def test_saturation_equivalent_potential_temperature(): """Test saturation equivalent potential temperature calculation.""" p = 700 * units.mbar From ecb5136fe2a1261717e59e6af98120c08fa2d240 Mon Sep 17 00:00:00 2001 From: kgoebber Date: Thu, 5 Oct 2023 14:59:16 -0500 Subject: [PATCH 2/3] make requested changes --- src/metpy/calc/thermo.py | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 94039aa6c6c..94227ed46e5 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1654,11 +1654,18 @@ def wetbulb_potential_temperature(pressure, temperature, dewpoint): r"""Calculate wet-bulb potential temperature. This calculation must be given an air parcel's pressure, temperature, and dewpoint. - The implementation uses the formula outlined in [DaviesJones2008]_: - + The implementation uses the formula outlined in [DaviesJones2008]_. First, theta-e is calculated then use the formula from [DaviesJones2008]_ + .. math:: \theta_w = \theta_e - + exp(\frac{a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4} + {1 + b_1 x + b_2 x^2 + b_3 x^3 + b_4 x^4}) + + where :math:`x = \theta_e / 273.15 K`. + + When :math:`\theta_e <= -173.15 K` then :math:`\theta_w = \theta_e`. + Parameters ---------- pressure: `pint.Quantity` @@ -1676,11 +1683,16 @@ def wetbulb_potential_temperature(pressure, temperature, dewpoint): wet-bulb potential temperature of the parcel """ - theta_e = equivalent_potential_temperature(pressure, temperature, dewpoint).magnitude - x = theta_e / 273.15 - a = 7.101574 - 20.68208 * x + 16.11182 * x ** 2 + 2.574631 * x ** 3 - 5.205688 * x ** 4 - b = 1 - 3.552497 * x + 3.781782 * x ** 2 - 0.6899655 * x ** 3 - 0.5929340 * x ** 4 - return units.Quantity(theta_e - np.exp((a) / (b)), 'kelvin') + theta_e = equivalent_potential_temperature(pressure, temperature, dewpoint) + x = theta_e.m_as('kelvin') / 273.15 + x2 = x * x + x3 = x2 * x + x4 = x2 * x2 + a = 7.101574 - 20.68208 * x + 16.11182 * x2 + 2.574631 * x3 - 5.205688 * x4 + b = 1 - 3.552497 * x + 3.781782 * x2 - 0.6899655 * x3 - 0.5929340 * x4 + + theta_w = units.Quantity(theta_e.m_as('kelvin') - np.exp(a / b), 'kelvin') + return np.where(theta_e <= units.Quantity(173.15, 'kelvin'), theta_e, theta_w) @exporter.export From 457898ea955bf9db503da2692e3ad61370f4d2dc Mon Sep 17 00:00:00 2001 From: Ryan May Date: Tue, 14 Nov 2023 17:11:12 -0700 Subject: [PATCH 3/3] MNT: wetbulb -> wet_bulb --- docs/_templates/overrides/metpy.calc.rst | 2 +- src/metpy/calc/thermo.py | 2 +- tests/calc/test_thermo.py | 8 ++++---- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index e24a5f559c5..19f05d47fc8 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -61,7 +61,7 @@ Moist Thermodynamics virtual_temperature virtual_temperature_from_dewpoint wet_bulb_temperature - wetbulb_potential_temperature + wet_bulb_potential_temperature Soundings diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 94227ed46e5..dae5401f890 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1650,7 +1650,7 @@ def saturation_equivalent_potential_temperature(pressure, temperature): broadcast=('pressure', 'temperature', 'dewpoint') ) @check_units('[pressure]', '[temperature]', '[temperature]') -def wetbulb_potential_temperature(pressure, temperature, dewpoint): +def wet_bulb_potential_temperature(pressure, temperature, dewpoint): r"""Calculate wet-bulb potential temperature. This calculation must be given an air parcel's pressure, temperature, and dewpoint. diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 8af5c439a8d..39faa743cc9 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -36,7 +36,7 @@ vapor_pressure, vertical_totals, vertical_velocity, vertical_velocity_pressure, virtual_potential_temperature, virtual_temperature, virtual_temperature_from_dewpoint, - wet_bulb_temperature, wetbulb_potential_temperature) + wet_bulb_potential_temperature, wet_bulb_temperature) from metpy.calc.thermo import _find_append_zero_crossings from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_nan, version_check) @@ -802,12 +802,12 @@ def test_equivalent_potential_temperature_masked(): assert_array_almost_equal(ept, expected, 3) -def test_wetbulb_potential_temperature(): - """Test wetbulb potential temperature calculation.""" +def test_wet_bulb_potential_temperature(): + """Test wet_bulb potential temperature calculation.""" p = 1000 * units.mbar t = 293. * units.kelvin td = 291. * units.kelvin - wpt = wetbulb_potential_temperature(p, t, td) + wpt = wet_bulb_potential_temperature(p, t, td) assert_almost_equal(wpt, 291.65839705486565 * units.kelvin, 3)