add wet bulb potential temperature function, based on existing equiva…

…lent potential

temperature function and formula given in DaviesJones2009

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

docs/api/references.rst

@@ -54,6 +54,10 @@ References
            *Mon. Wea. Rev.*, **137**, 3137-3148,
            doi: `10.1175/2009mwr2774.1 <https://doi.org/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 <https://doi.org/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 <https://ejssm.org/archives/2006/vol-1-3-2006/>`_.
            *Electronic J. Severe Storms Meteor.*, **1** (3), 1-22.

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(

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