Merge pull request #1022 from sgdecker/fix_iss945

Fix crash with mixed-layer LFC calculation
Unidata · Apr 9, 2019 · 1b31cda · 1b31cda
2 parents f24fa7c + 1f49602
commit 1b31cda
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 5 deletions.
diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py
@@ -286,6 +286,47 @@ def test_lfc_ml():
     assert_almost_equal(lfc_temp, -1.862 * units.degC, 2)
 
 
+def test_lfc_ml2():
+    """Test a mixed-layer LFC calculation that previously crashed."""
+    levels = np.array([1024.95703125, 1016.61474609, 1005.33056641, 991.08544922, 973.4163208,
+                       951.3381958, 924.82836914, 898.25482178, 873.46124268, 848.69830322,
+                       823.92553711, 788.49304199, 743.44580078, 700.50970459, 659.62017822,
+                       620.70861816, 583.69421387, 548.49719238, 515.03826904, 483.24401855,
+                       453.0418396, 424.36477661, 397.1505127, 371.33441162, 346.85922241,
+                       323.66995239, 301.70935059, 280.92651367, 261.27053833, 242.69168091,
+                       225.14237976, 208.57781982, 192.95333862, 178.22599792, 164.39630127,
+                       151.54336548, 139.68635559, 128.74923706, 118.6588974, 109.35111237,
+                       100.76405334, 92.84288025, 85.53556824, 78.79430389, 72.57549286,
+                       66.83885193, 61.54678726, 56.66480637, 52.16108322]) * units.mbar
+    temperatures = np.array([6.00750732, 5.14892578, 4.177948, 3.00268555, 1.55535889,
+                             -0.25527954, -1.93988037, -3.57766724, -4.40600586, -4.19238281,
+                             -3.71185303, -4.47943115, -6.81280518, -8.08685303, -8.41287231,
+                             -10.79302979, -14.13262939, -16.85784912, -19.51675415,
+                             -22.28689575, -24.99938965, -27.79664612, -30.90414429,
+                             -34.49435425, -38.438797, -42.27981567, -45.99230957,
+                             -49.75340271, -53.58230591, -57.30686951, -60.76026917,
+                             -63.92070007, -66.72470093, -68.97846985, -70.4264679,
+                             -71.16407776, -71.53797913, -71.64375305, -71.52735901,
+                             -71.53523254, -71.61097717, -71.92687988, -72.68682861,
+                             -74.129776, -76.02471924, -76.88977051, -76.26008606,
+                             -75.90351868, -76.15809631]) * units.celsius
+    dewpoints = np.array([4.50012302, 3.42483997, 2.78102994, 2.24474645, 1.593485, 0.94408154,
+                          -3.8044982, -3.55629468, -9.7376976, -10.2950449, -9.67498302,
+                          -10.30486488, -8.70559597, -8.71669006, -12.66509628, -18.6697197,
+                          -23.00351334, -29.46240425, -36.82178497, -41.68824768, -44.50320816,
+                          -48.54426575, -52.50753403, -51.09564209, -48.92690659, -49.97380829,
+                          -51.57516098, -52.62096405, -54.24332809, -57.09109879, -60.5596199,
+                          -63.93486404, -67.07530212, -70.01263428, -72.9258728, -76.12271881,
+                          -79.49847412, -82.2350769, -83.91127014, -84.95665741, -85.61238861,
+                          -86.16391754, -86.7653656, -87.34436035, -87.87495422, -88.34281921,
+                          -88.74453735, -89.04680634, -89.26436615]) * units.celsius
+    __, t_mixed, td_mixed = mixed_parcel(levels, temperatures, dewpoints)
+    mixed_parcel_prof = parcel_profile(levels, t_mixed, td_mixed)
+    lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints, mixed_parcel_prof)
+    assert_almost_equal(lfc_pressure, 1001.532 * units.mbar, 2)
+    assert_almost_equal(lfc_temp, 4.17 * units.degC, 2)
+
+
 def test_no_lfc():
     """Test LFC calculation when there is no LFC in the data."""
     levels = np.array([959., 867.9, 779.2, 647.5, 472.5, 321.9, 251.]) * units.mbar

diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py
@@ -421,12 +421,18 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
             x, y = this_lcl
             return x, y
 
-    # LFC exists and is not LCL. Make sure it is above the LCL.
+    # LFC exists. Make sure it is no lower than the LCL
     else:
-        idx = x < lcl(pressure[0], temperature[0], dewpt[0])[0]
-        x = x[idx]
-        y = y[idx]
-        return x[0], y[0]
+        idx = x < this_lcl[0]
+        # LFC height < LCL height, so set LFC = LCL
+        if not any(idx):
+            x, y = this_lcl
+            return x, y
+        # Otherwise, make select first candidate LFC above the LCL
+        else:
+            x = x[idx]
+            y = y[idx]
+            return x[0], y[0]
 
 
 @exporter.export