Merge pull request #3751 from ahijevyc/use_parcel_prof_for_parcel_mix_ratio

Use parcel prof for parcel mix ratio

Author: dopplershift

src/metpy/calc/thermo.py

@@ -2642,7 +2642,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
     >>> prof = parcel_profile(p, T[0], Td[0]).to('degC')
     >>> # calculate surface based CAPE/CIN
     >>> cape_cin(p, T, Td, prof)
-    (<Quantity(4703.77308, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
+    (<Quantity(4910.59044, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
 
     See Also
     --------
@@ -2686,8 +2686,11 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
 
     # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated
     # based on the temperature above the LCL
-    parcel_mixing_ratio = np.where(below_lcl, saturation_mixing_ratio(pressure, dewpoint),
-                                   saturation_mixing_ratio(pressure, temperature))
+    parcel_mixing_ratio = np.where(
+        below_lcl,
+        saturation_mixing_ratio(pressure[0], dewpoint[0]),
+        saturation_mixing_ratio(pressure, parcel_profile)
+    )
 
     # Convert the temperature/parcel profile to virtual temperature
     temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint)
@@ -3280,7 +3283,7 @@ def most_unstable_cape_cin(pressure, temperature, dewpoint, **kwargs):
     >>> Td = dewpoint_from_relative_humidity(T, rh)
     >>> # calculate most unstbale CAPE/CIN
     >>> most_unstable_cape_cin(p, T, Td)
-    (<Quantity(4703.77308, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
+    (<Quantity(4910.59044, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
 
     See Also
     --------
@@ -3356,7 +3359,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
     >>> # calculate dewpoint
     >>> Td = dewpoint_from_relative_humidity(T, rh)
     >>> mixed_layer_cape_cin(p, T, Td, depth=50 * units.hPa)
-    (<Quantity(711.239032, 'joule / kilogram')>, <Quantity(-5.48053989, 'joule / kilogram')>)
+    (<Quantity(740.140593, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
 
     See Also
     --------

tests/calc/test_thermo.py

@@ -744,28 +744,6 @@ def test_lfc_profile_nan_with_parcel_profile():
     assert_almost_equal(lfc_temperature, 9.6977 * units.degC, 3)
 
 
-def test_sensitive_sounding():
-    """Test quantities for a sensitive sounding (#902)."""
-    # This sounding has a very small positive area in the low level. It's only captured
-    # properly if the parcel profile includes the LCL, otherwise it breaks LFC and CAPE
-    p = units.Quantity([1004., 1000., 943., 928., 925., 850., 839., 749., 700., 699.,
-                        603., 500., 404., 400., 363., 306., 300., 250., 213., 200.,
-                        176., 150.], 'hectopascal')
-    t = units.Quantity([25.1, 24.5, 20.2, 21.6, 21.4, 20.4, 20.2, 14.4, 13.2, 13., 6.8, -3.3,
-                        -13.1, -13.7, -17.9, -25.5, -26.9, -37.9, -46.7, -48.7, -52.1, -58.9],
-                       'degC')
-    td = units.Quantity([21.9, 22.1, 19.2, 20.5, 20.4, 18.4, 17.4, 8.4, -2.8, -3.0, -15.2,
-                         -20.3, -29.1, -27.7, -24.9, -39.5, -41.9, -51.9, -60.7, -62.7, -65.1,
-                         -71.9], 'degC')
-    lfc_pressure, lfc_temp = lfc(p, t, td)
-    assert_almost_equal(lfc_pressure, 952.8445 * units.mbar, 2)
-    assert_almost_equal(lfc_temp, 20.94469 * units.degC, 2)
-
-    pos, neg = surface_based_cape_cin(p, t, td)
-    assert_almost_equal(pos, 0.106791 * units('J/kg'), 3)
-    assert_almost_equal(neg, -282.620677 * units('J/kg'), 3)
-
-
 def test_lfc_sfc_precision():
     """Test LFC when there are precision issues with the parcel path."""
     levels = np.array([839., 819.4, 816., 807., 790.7, 763., 736.2,
@@ -1200,8 +1178,8 @@ def test_cape_cin():
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     parcel_prof = parcel_profile(p, temperature[0], dewpoint[0])
     cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
-    assert_almost_equal(cape, 215.056976 * units('joule / kilogram'), 2)
-    assert_almost_equal(cin, -9.94798721 * units('joule / kilogram'), 2)
+    assert_almost_equal(cape, 228.61081997000744 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -20.8938 * units('joule / kilogram'), 2)
 
 
 def test_cape_cin_no_el():
@@ -1211,8 +1189,8 @@ def test_cape_cin_no_el():
     dewpoint = np.array([19., -11.2, -10.8, -10.4]) * units.celsius
     parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC')
     cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
-    assert_almost_equal(cape, 12.74623773 * units('joule / kilogram'), 2)
-    assert_almost_equal(cin, -9.947987213 * units('joule / kilogram'), 2)
+    assert_almost_equal(cape, 11.10149 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -20.8938 * units('joule / kilogram'), 2)
 
 
 def test_cape_cin_no_lfc():
@@ -1609,8 +1587,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, 215.05697634 * units('joule / kilogram'), 2)
-    assert_almost_equal(cin, -33.0633599455 * units('joule / kilogram'), 2)
+    assert_almost_equal(cape, 228.61081997000744 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -52.46449098033761 * units('joule / kilogram'), 2)
 
 
 def test_surface_based_cape_cin_with_xarray():
@@ -1633,8 +1611,8 @@ def test_surface_based_cape_cin_with_xarray():
         data['temperature'],
         data['dewpoint']
     )
-    assert_almost_equal(cape, 215.056976346 * units('joule / kilogram'), 2)
-    assert_almost_equal(cin, -33.0633599455 * units('joule / kilogram'), 2)
+    assert_almost_equal(cape, 228.61081997000744 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -52.46449098033761 * units('joule / kilogram'), 2)
 
 
 def test_profile_with_nans():
@@ -1674,8 +1652,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, 215.056976346 * units('joule / kilogram'), 2)
-    assert_almost_equal(mucin, -33.0633599455 * units('joule / kilogram'), 2)
+    assert_almost_equal(mucape, 228.61081997000744 * units('joule / kilogram'), 2)
+    assert_almost_equal(mucin, -52.46449098033761 * units('joule / kilogram'), 2)
 
 
 def test_most_unstable_cape_cin():
@@ -1684,8 +1662,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, 173.749389796 * units('joule / kilogram'), 4)
-    assert_almost_equal(mucin, -20.968278741 * units('joule / kilogram'), 4)
+    assert_almost_equal(mucape, 189.41067504060692 * units('joule / kilogram'), 4)
+    assert_almost_equal(mucin, -26.225902591840672 * units('joule / kilogram'), 4)
 
 
 def test_mixed_parcel():
@@ -1705,17 +1683,17 @@ def test_mixed_layer_cape_cin(multiple_intersections):
     """Test the calculation of mixed layer cape/cin."""
     pressure, temperature, dewpoint = multiple_intersections
     mlcape, mlcin = mixed_layer_cape_cin(pressure, temperature, dewpoint)
-    assert_almost_equal(mlcape, 1132.706800436 * units('joule / kilogram'), 2)
-    assert_almost_equal(mlcin, -13.4809966289 * units('joule / kilogram'), 2)
+    assert_almost_equal(mlcape, 1143.3981 * units('joule / kilogram'), 2)
+    assert_almost_equal(mlcin, -16.240379524041845 * units('joule / kilogram'), 2)
 
 
 def test_mixed_layer_cape_cin_bottom_pressure(multiple_intersections):
     """Test the calculation of mixed layer cape/cin with a specified bottom pressure."""
     pressure, temperature, dewpoint = multiple_intersections
     mlcape_middle, mlcin_middle = mixed_layer_cape_cin(pressure, temperature, dewpoint,
                                                        parcel_start_pressure=903 * units.hPa)
-    assert_almost_equal(mlcape_middle, 1177.86 * units('joule / kilogram'), 2)
-    assert_almost_equal(mlcin_middle, -37. * units('joule / kilogram'), 2)
+    assert_almost_equal(mlcape_middle, 1200.528254 * units('joule / kilogram'), 2)
+    assert_almost_equal(mlcin_middle, -46.99243161905505 * units('joule / kilogram'), 2)
 
 
 def test_dcape():
@@ -2220,17 +2198,17 @@ def test_cape_cin_top_el_lfc(multiple_intersections):
     levels, temperatures, dewpoints = multiple_intersections
     parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC')
     cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='top')
-    assert_almost_equal(cape, 1345.18884959 * units('joule / kilogram'), 3)
-    assert_almost_equal(cin, -35.179268355 * units('joule / kilogram'), 3)
+    assert_almost_equal(cape, 1371.747661 * units('joule / kilogram'), 3)
+    assert_almost_equal(cin, -46.084968610767184 * units('joule / kilogram'), 3)
 
 
 def test_cape_cin_bottom_el_lfc(multiple_intersections):
     """Test using LFC/EL options for CAPE/CIN."""
     levels, temperatures, dewpoints = multiple_intersections
     parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC')
     cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_el='bottom')
-    assert_almost_equal(cape, 4.57262449 * units('joule / kilogram'), 3)
-    assert_almost_equal(cin, -5.9471237534 * units('joule / kilogram'), 3)
+    assert_almost_equal(cape, 4.76113 * units('joule / kilogram'), 3)
+    assert_almost_equal(cin, -7.058249615394496 * units('joule / kilogram'), 3)
 
 
 def test_cape_cin_wide_el_lfc(multiple_intersections):
@@ -2239,8 +2217,8 @@ def test_cape_cin_wide_el_lfc(multiple_intersections):
     parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC')
     cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='wide',
                          which_el='wide')
-    assert_almost_equal(cape, 1345.1888496 * units('joule / kilogram'), 3)
-    assert_almost_equal(cin, -35.179268355 * units('joule / kilogram'), 3)
+    assert_almost_equal(cape, 1371.747661 * units('joule / kilogram'), 3)
+    assert_almost_equal(cin, -46.084968610767184 * units('joule / kilogram'), 3)
 
 
 def test_cape_cin_custom_profile():
@@ -2250,7 +2228,7 @@ def test_cape_cin_custom_profile():
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) + 5 * units.delta_degC
     cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
-    assert_almost_equal(cape, 1650.61208729 * units('joule / kilogram'), 2)
+    assert_almost_equal(cape, 1785.7865489238366 * units('joule / kilogram'), 2)
     assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 2)
 
 
@@ -2344,7 +2322,7 @@ def test_cape_cin_value_error():
                          -35.9, -26.7, -37.7, -43.1, -33.9, -40.9, -46.1, -34.9, -33.9,
                          -33.7, -33.3, -42.5, -50.3, -49.7, -49.5, -58.3, -61.3]) * units.degC
     cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint)
-    expected_cape, expected_cin = 2098.688061 * units('joules/kg'), 0.0 * units('joules/kg')
+    expected_cape, expected_cin = 2182.6035338 * units('joules/kg'), 0.0 * units('joules/kg')
     assert_almost_equal(cape, expected_cape, 3)
     assert_almost_equal(cin, expected_cin, 3)