diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst
index f7a78e3ba58..57a438210f9 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
+      wet_bulb_potential_temperature
 
 
 Soundings
diff --git a/docs/api/references.rst b/docs/api/references.rst
index 7191be729fa..787cc6581da 100644
--- a/docs/api/references.rst
+++ b/docs/api/references.rst
@@ -60,6 +60,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.
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index 716886e9e15..eeba567cd6a 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -1644,6 +1644,57 @@ 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 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.
+    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`
+        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)
+    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
 @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 462a53786e8..118902b3f27 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_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,6 +802,15 @@ def test_equivalent_potential_temperature_masked():
     assert_array_almost_equal(ept, expected, 3)
 
 
+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 = wet_bulb_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