Merge pull request #1023 from sgdecker/fix_lfc_assumptions

Fix poor assumptions in lfc() logic

metpy/calc/tests/test_thermo.py

@@ -310,7 +310,7 @@ def test_lfc_ml2():
                              -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,
+    dewpoints = np.array([4.50012302, 3.42483997, 2.78102994, 2.24474645, 1.593485, -0.9440815,
                           -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,
@@ -322,9 +322,22 @@ def test_lfc_ml2():
                           -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)
+    lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints, mixed_parcel_prof, td_mixed)
+    assert_almost_equal(lfc_pressure, 962.34 * units.mbar, 2)
+    assert_almost_equal(lfc_temp, 0.767 * units.degC, 2)
+
+
+def test_lfc_intersection():
+    """Test LFC calculation when LFC is below a tricky intersection."""
+    p = np.array([1024.957, 930., 924.828, 898.255, 873.461, 848.698, 823.926,
+                  788.493]) * units('hPa')
+    t = np.array([6.008, -10., -6.94, -8.58, -4.41, -4.19, -3.71, -4.48]) * units('degC')
+    td = np.array([5., -10., -7., -9., -4.5, -4.2, -3.8, -4.5]) * units('degC')
+    _, mlt, mltd = mixed_parcel(p, t, td)
+    ml_profile = parcel_profile(p, mlt, mltd)
+    mllfc_p, mllfc_t = lfc(p, t, td, ml_profile, mltd)
+    assert_almost_equal(mllfc_p, 982.762 * units.hPa, 2)
+    assert_almost_equal(mllfc_t, 272.045 * units.kelvin, 2)
 
 
 def test_no_lfc():

metpy/calc/thermo.py

@@ -362,7 +362,7 @@ def _lcl_iter(p, p0, w, t):
 @exporter.export
 @preprocess_xarray
 @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]')
-def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
+def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_start=None):
     r"""Calculate the level of free convection (LFC).
 
     This works by finding the first intersection of the ideal parcel path and
@@ -379,6 +379,9 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
     parcel_temperature_profile: `pint.Quantity`, optional
         The parcel temperature profile from which to calculate the LFC. Defaults to the
         surface parcel profile.
+    dewpt_start: `pint.Quantity`, optional
+        The dewpoint of the parcel for which to calculate the LFC. Defaults to the surface
+        dewpoint.
 
     Returns
     -------
@@ -397,13 +400,21 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
         temperature = temperature.to('degC')
         parcel_temperature_profile = parcel_temperature_profile.to('degC')
 
-    # The parcel profile and data have the same first data point, so we ignore
-    # that point to get the real first intersection for the LFC calculation.
-    x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:],
-                              temperature[1:], direction='increasing')
+    if dewpt_start is None:
+        dewpt_start = dewpt[0]
+
+    # The parcel profile and data may have the same first data point.
+    # If that is the case, ignore that point to get the real first
+    # intersection for the LFC calculation.
+    if np.isclose(parcel_temperature_profile[0].m, temperature[0].m):
+        x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:],
+                                  temperature[1:], direction='increasing')
+    else:
+        x, y = find_intersections(pressure, parcel_temperature_profile,
+                                  temperature, direction='increasing')
 
     # Compute LCL for this parcel for future comparisons
-    this_lcl = lcl(pressure[0], temperature[0], dewpt[0])
+    this_lcl = lcl(pressure[0], parcel_temperature_profile[0], dewpt_start)
 
     # The LFC could:
     # 1) Not exist