Update thickness_hydrostatic to have an option to return full profile

Author: Z-Richard

src/metpy/calc/thermo.py

@@ -3338,16 +3338,9 @@ def moist_static_energy(height, temperature, specific_humidity):
 )
 def thickness_hydrostatic(pressure, temperature, mixing_ratio=None,
                           molecular_weight_ratio=mpconsts.nounit.epsilon, bottom=None,
-                          depth=None):
+                          depth=None, full_profile=False):
     r"""Calculate the thickness of a layer via the hypsometric equation.
 
-    This thickness calculation uses the pressure and temperature profiles (and optionally
-    mixing ratio) via the hypsometric equation with virtual temperature adjustment.
-
-    .. math:: Z_2 - Z_1 = -\frac{R_d}{g} \int_{p_1}^{p_2} T_v d\ln p,
-
-    Which is based off of Equation 3.24 in [Hobbs2006]_.
-
     This assumes a hydrostatic atmosphere. Layer bottom and depth specified in pressure.
 
     Parameters
@@ -3374,6 +3367,12 @@ def thickness_hydrostatic(pressure, temperature, mixing_ratio=None,
         The depth of the layer in hPa. Defaults to the full profile if bottom is not given,
         and 100 hPa if bottom is given.
 
+    full_profile : `bool`, optional
+        Controls whether or not to return the full thickness profile, i.e., the thickness
+        between individual layers as a `Pint.Quantity` array or the thickness of
+        the whole atmospheric profile as a `Pint.Quantity`. Defaults to False, which returns
+        the thickness of the atmosphere as a single value.
+
     Returns
     -------
     `pint.Quantity`
@@ -3427,6 +3426,13 @@ def thickness_hydrostatic(pressure, temperature, mixing_ratio=None,
     Since this function returns scalar values when given a profile, this will return Pint
     Quantities even when given xarray DataArray profiles.
 
+    This thickness calculation uses the pressure and temperature profiles (and optionally
+    mixing ratio) via the hypsometric equation with virtual temperature adjustment.
+
+    .. math:: Z_2 - Z_1 = -\frac{R_d}{g} \int_{p_1}^{p_2} T_v d\ln p,
+
+    Which is based off of Equation 3.24 in [Hobbs2006]_.
+
     .. versionchanged:: 1.0
        Renamed ``mixing`` parameter to ``mixing_ratio``
 
@@ -3470,31 +3476,20 @@ def thickness_hydrostatic(pressure, temperature, mixing_ratio=None,
             )
 
     # Take the integral
-    return (
-        -mpconsts.nounit.Rd / mpconsts.nounit.g
-        * np.trapz(layer_virttemp, np.log(layer_p))
-    )
+    diff_logp = np.diff(np.log(layer_p))
+    thickness = (-mpconsts.nounit.Rd / mpconsts.nounit.g *
+                diff_logp * (layer_virttemp[1:] + layer_virttemp[:-1]) / 2.)
+    return thickness if full_profile else np.sum(thickness)
 
 
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]', '[temperature]')
 def thickness_hydrostatic_from_relative_humidity(pressure, temperature, relative_humidity,
-                                                 bottom=None, depth=None):
+                                                 bottom=None, depth=None, full_profile=False):
     r"""Calculate the thickness of a layer given pressure, temperature and relative humidity.
 
-    Similar to ``thickness_hydrostatic``, this thickness calculation uses the pressure,
-    temperature, and relative humidity profiles via the hypsometric equation with virtual
-    temperature adjustment
-
-    .. math:: Z_2 - Z_1 = -\frac{R_d}{g} \int_{p_1}^{p_2} T_v d\ln p,
-
-    which is based off of Equation 3.24 in [Hobbs2006]_. Virtual temperature is calculated
-    from the profiles of temperature and relative humidity.
-
-    This assumes a hydrostatic atmosphere.
-
-    Layer bottom and depth specified in pressure.
+    This assumes a hydrostatic atmosphere. Layer bottom and depth specified in pressure.
 
     Parameters
     ----------
@@ -3516,6 +3511,12 @@ def thickness_hydrostatic_from_relative_humidity(pressure, temperature, relative
         The depth of the layer in hPa. Defaults to the full profile if bottom is not given,
         and 100 hPa if bottom is given.
 
+    full_profile : `bool`, optional
+        Controls whether or not to return the full thickness profile, i.e., the thickness
+        between individual layers as a `Pint.Quantity` array or the thickness of
+        the whole atmospheric profile as a `Pint.Quantity`. Defaults to False, which returns
+        the thickness of the atmosphere as a single value.
+
     Returns
     -------
     `pint.Quantity`
@@ -3557,11 +3558,19 @@ def thickness_hydrostatic_from_relative_humidity(pressure, temperature, relative
     Since this function returns scalar values when given a profile, this will return Pint
     Quantities even when given xarray DataArray profiles.
 
+    Similar to ``thickness_hydrostatic``, this thickness calculation uses the pressure,
+    temperature, and relative humidity profiles via the hypsometric equation with virtual
+    temperature adjustment
+
+    .. math:: Z_2 - Z_1 = -\frac{R_d}{g} \int_{p_1}^{p_2} T_v d\ln p,
+
+    which is based off of Equation 3.24 in [Hobbs2006]_. Virtual temperature is calculated
+    from the profiles of temperature and relative humidity.
     """
     mixing = mixing_ratio_from_relative_humidity(pressure, temperature, relative_humidity)
 
     return thickness_hydrostatic(pressure, temperature, mixing_ratio=mixing, bottom=bottom,
-                                 depth=depth)
+                                 depth=depth, full_profile=full_profile)
 
 
 @exporter.export

tests/calc/test_thermo.py

@@ -1607,6 +1607,17 @@ def test_thickness_hydrostatic_isothermal_subset():
     assert_almost_equal(thickness, 4242.527 * units.m, 2)
 
 
+def test_thickness_hydrostatic_full_profile():
+    """Test the thickness calculation for a moist layer."""
+    pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa
+    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.degC
+    mixing = np.array([0.01458, 0.00209, 0.00224, 0.00240, 0.00256, 0.00010])
+    thickness = thickness_hydrostatic(pressure, temperature, mixing_ratio=mixing, 
+                                      full_profile=True)
+    result = [1780.778, 306.126, 304.513, 281.455, 7218.833] * units.m
+    assert_almost_equal(thickness, result, 2)
+
+
 def test_thickness_hydrostatic_from_relative_humidity():
     """Test the thickness calculation for a moist layer using RH data."""
     pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa