Using Virtual Temperature in CAPE/CIN calculations

src/metpy/calc/thermo.py

@@ -715,10 +715,19 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
                 units.Quantity(np.nan, temperature.units))
 
 
+class ParcelPathAssumptions():
+    """Holds assumptions made about the parcel path during calculations."""
+
+    def __init__(self):
+        self.use_virtual_temperature = True
+        self.moist_adiabat = 'pseudoadiabatic'
+
+
 @exporter.export
 @preprocess_and_wrap(wrap_like='pressure')
 @check_units('[pressure]', '[temperature]', '[temperature]')
-def parcel_profile(pressure, temperature, dewpoint):
+def parcel_profile(pressure, temperature, dewpoint,
+                   assumptions=ParcelPathAssumptions()):
     r"""Calculate the profile a parcel takes through the atmosphere.
 
     The parcel starts at `temperature`, and `dewpoint`, lifted up
@@ -760,7 +769,8 @@ def parcel_profile(pressure, temperature, dewpoint):
        Renamed ``dewpt`` parameter to ``dewpoint``
 
     """
-    _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint)
+    _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint,
+                                                  assumptions=assumptions)
     return concatenate((t_l, t_u))
 
 
@@ -911,7 +921,8 @@ def _check_pressure(pressure):
     return np.all(pressure[:-1] >= pressure[1:])
 
 
-def _parcel_profile_helper(pressure, temperature, dewpoint):
+def _parcel_profile_helper(pressure, temperature, dewpoint,
+                           assumptions=ParcelPathAssumptions()):
     """Help calculate parcel profiles.
 
     Returns the temperature and pressure, above, below, and including the LCL. The
@@ -924,6 +935,11 @@ def _parcel_profile_helper(pressure, temperature, dewpoint):
         Pressure increases between at least two points in your sounding.
         Using scipy.signal.medfilt may fix this."""
         raise InvalidSoundingError(msg)
+
+    if assumptions.use_virtual_temperature:
+        specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint)
+        mixing_ratio = mixing_ratio_from_specific_humidity(specific_humidity)
+        temperature = virtual_temperature(temperature, mixing_ratio)
 
     # Find the LCL
     press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpoint)
@@ -955,6 +971,7 @@ def _parcel_profile_helper(pressure, temperature, dewpoint):
     temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units)
     temp_upper = temp_upper[::-1][indices]
 
+
     # Return profile pieces
     return (press_lower[:-1], press_lcl, press_upper[1:],
             temp_lower[:-1], temp_lcl, temp_upper[1:])

tests/calc/test_thermo.py

@@ -1261,8 +1261,8 @@ def test_surface_based_cape_cin(array_class):
     temperature = array_class([22.2, 14.6, 12., 9.4, 7., -38.], units.celsius)
     dewpoint = array_class([19., -11.2, -10.8, -10.4, -10., -53.2], units.celsius)
     cape, cin = surface_based_cape_cin(p, temperature, dewpoint)
-    assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2)
-    assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2)
+    assert_almost_equal(cape, 327.1348028 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -40.0654774 * units('joule / kilogram'), 2)
 
 
 def test_surface_based_cape_cin_with_xarray():
@@ -1285,8 +1285,8 @@ def test_surface_based_cape_cin_with_xarray():
         data['temperature'],
         data['dewpoint']
     )
-    assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2)
-    assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2)
+    assert_almost_equal(cape, 327.1348028 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -40.0654774 * units('joule / kilogram'), 2)
 
 
 def test_profile_with_nans():
@@ -1326,8 +1326,8 @@ def test_most_unstable_cape_cin_surface():
     temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint)
-    assert_almost_equal(mucape, 75.0535446 * units('joule / kilogram'), 2)
-    assert_almost_equal(mucin, -136.685967 * units('joule / kilogram'), 2)
+    assert_almost_equal(mucape, 327.1348028 * units('joule / kilogram'), 2)
+    assert_almost_equal(mucin, -40.0654774 * units('joule / kilogram'), 2)
 
 
 def test_most_unstable_cape_cin():
@@ -1336,8 +1336,8 @@ def test_most_unstable_cape_cin():
     temperature = np.array([18.2, 22.2, 17.4, 10., 0., 15]) * units.celsius
     dewpoint = np.array([19., 19., 14.3, 0., -10., 0.]) * units.celsius
     mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint)
-    assert_almost_equal(mucape, 157.11404 * units('joule / kilogram'), 4)
-    assert_almost_equal(mucin, -31.8406578 * units('joule / kilogram'), 4)
+    assert_almost_equal(mucape, 174.66467 * units('joule / kilogram'), 4)
+    assert_almost_equal(mucin, 0 * units('joule / kilogram'), 4)
 
 
 def test_mixed_parcel():