From 882f84f74de8aef1a74161ffcfb7663087a222e9 Mon Sep 17 00:00:00 2001 From: ahijevyc Date: Wed, 23 Apr 2025 16:28:13 -0600 Subject: [PATCH 1/4] Use parcel temperature not env temperature when getting parcel mixing ratio --- src/metpy/calc/thermo.py | 7 ++++-- tests/calc/test_thermo.py | 48 +++++++++++++++++++-------------------- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 10c221deef5..46b235a3c80 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -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) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 01e547c2135..52455545721 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -1200,8 +1200,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 +1211,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 +1609,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 +1633,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 +1674,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 +1684,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,8 +1705,8 @@ 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): @@ -1714,8 +1714,8 @@ def test_mixed_layer_cape_cin_bottom_pressure(multiple_intersections): 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,8 +2220,8 @@ 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): @@ -2229,8 +2229,8 @@ def test_cape_cin_bottom_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_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 +2239,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 +2250,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 +2344,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) From e98afcf67540526f1c283f88c113d60b5a5a4b35 Mon Sep 17 00:00:00 2001 From: ahijevyc Date: Thu, 24 Apr 2025 14:47:09 -0600 Subject: [PATCH 2/4] skip test for sounding sensitive to LCL inclusion --- tests/calc/test_thermo.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 52455545721..65162befa7c 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -744,6 +744,7 @@ def test_lfc_profile_nan_with_parcel_profile(): assert_almost_equal(lfc_temperature, 9.6977 * units.degC, 3) +@pytest.mark.skip(reason="Skipping because of difficulty constructing a sounding sensitive to inclusion of LCL after bug fix (#3751)") 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 From 269f7d5524cb9ac4ce73315e86433c0adc35fa88 Mon Sep 17 00:00:00 2001 From: ahijevyc Date: Thu, 24 Apr 2025 15:02:54 -0600 Subject: [PATCH 3/4] fixed cape/cin in function docstrings --- src/metpy/calc/thermo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 46b235a3c80..971be0e1b38 100644 --- a/src/metpy/calc/thermo.py +++ b/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) - (, ) + (, ) See Also -------- @@ -3283,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) - (, ) + (, ) See Also -------- @@ -3359,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) - (, ) + (, ) See Also -------- From ce6ba87c0539ccf14a9df922034f55ed9126eefe Mon Sep 17 00:00:00 2001 From: ahijevyc Date: Thu, 24 Apr 2025 15:29:42 -0600 Subject: [PATCH 4/4] remove sensitive_sounding test instead of skipping --- tests/calc/test_thermo.py | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 65162befa7c..76d10465f54 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -744,29 +744,6 @@ def test_lfc_profile_nan_with_parcel_profile(): assert_almost_equal(lfc_temperature, 9.6977 * units.degC, 3) -@pytest.mark.skip(reason="Skipping because of difficulty constructing a sounding sensitive to inclusion of LCL after bug fix (#3751)") -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,