From 5ae1cf4e767112c0f1d958ea0827a920f3eead71 Mon Sep 17 00:00:00 2001
From: zbruick <zachary.bruick@valpo.edu>
Date: Thu, 13 Jun 2019 10:10:27 -0600
Subject: [PATCH] This update to the lfc() function allows for multiple LFCs to
 be identified within a sounding. This fixes #992 and begins to address

---
 metpy/calc/tests/test_thermo.py | 28 ++++++++++++++++++++++++++
 metpy/calc/thermo.py            | 35 +++++++++++++++++++++++++++------
 2 files changed, 57 insertions(+), 6 deletions(-)

diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py
index 2bb67054b6..41b36bec3d 100644
--- a/metpy/calc/tests/test_thermo.py
+++ b/metpy/calc/tests/test_thermo.py
@@ -1201,6 +1201,34 @@ def test_lfc_not_below_lcl():
     assert_almost_equal(lfc_temp, 6.4992871 * units.celsius, 6)
 
 
+def test_multiple_lfcs():
+    """Test sounding with multiple LFCs.
+
+    If return_all=False, return lowest-pressure LFC.
+    If return_all=True, return all LFCs
+
+    """
+    levels = np.array([966., 937.2, 925., 904.6, 872.6, 853., 850., 836., 821., 811.6, 782.3,
+                       754.2, 726.9, 700., 648.9, 624.6, 601.1, 595., 587., 576., 555.7,
+                       534.2, 524., 500., 473.3, 400., 384.5, 358., 343., 308.3, 300., 276.,
+                       273., 268.5, 250., 244.2, 233., 200.]) * units.mbar
+    temperatures = np.array([18.2, 16.8, 16.2, 15.1, 13.3, 12.2, 12.4, 14., 14.4,
+                             13.7, 11.4, 9.1, 6.8, 4.4, -1.4, -4.4, -7.3, -8.1,
+                             -7.9, -7.7, -8.7, -9.8, -10.3, -13.5, -17.1, -28.1, -30.7,
+                             -35.3, -37.1, -43.5, -45.1, -49.9, -50.4, -51.1, -54.1, -55.,
+                             -56.7, -57.5]) * units.degC
+    dewpoints = np.array([16.9, 15.9, 15.5, 14.2, 12.1, 10.8, 8.6, 0., -3.6, -4.4,
+                          -6.9, -9.5, -12., -14.6, -15.8, -16.4, -16.9, -17.1, -27.9, -42.7,
+                          -44.1, -45.6, -46.3, -45.5, -47.1, -52.1, -50.4, -47.3, -57.1,
+                          -57.9, -58.1, -60.9, -61.4, -62.1, -65.1, -65.6,
+                          -66.7, -70.5]) * units.degC
+    lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints)
+    lfc_pressure_all, _ = lfc(levels, temperatures, dewpoints, return_all=True)
+    assert_almost_equal(lfc_pressure, 705.5170051 * units.mbar, 6)
+    assert_almost_equal(lfc_temp, 4.8922235 * units.degC, 6)
+    assert_almost_equal(len(lfc_pressure_all), 2, 0)
+
+
 def test_cape_cin_custom_profile():
     """Test the CAPE and CIN calculation with a custom profile passed to LFC and EL."""
     p = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py
index fae1a08d1c..a98b669d42 100644
--- a/metpy/calc/thermo.py
+++ b/metpy/calc/thermo.py
@@ -362,7 +362,8 @@ 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, dewpt_start=None):
+def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_start=None,
+        return_all=False):
     r"""Calculate the level of free convection (LFC).
 
     This works by finding the first intersection of the ideal parcel path and
@@ -382,11 +383,13 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta
     dewpt_start: `pint.Quantity`, optional
         The dewpoint of the parcel for which to calculate the LFC. Defaults to the surface
         dewpoint.
+    return_all: bool, optional
+        Option to return all LFCs in profile. Defaults to the lowest-pressure LFC.
 
     Returns
     -------
     `pint.Quantity`
-        The LFC pressure and temperature
+        The LFC pressure and temperature, or lists of same if return_all=True
 
     See Also
     --------
@@ -439,11 +442,31 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta
         if not any(idx):
             x, y = this_lcl
             return x, y
-        # Otherwise, make select first candidate LFC above the LCL
+        # Otherwise, find all LFCs that exist above the LCL
+        # What is returned depends on return_all flag as described in the docstring
         else:
-            x = x[idx]
-            y = y[idx]
-            return x[0], y[0]
+            x, y = _multiple_lfcs(pressure, temperature, parcel_temperature_profile, x, y,
+                                  this_lcl[0], return_all)
+            return x, y
+
+
+def _multiple_lfcs(pressure, temperature, profile_temperatures, intersect_pressures,
+                   intersect_temperatures, lcl_pressure, return_all):
+    """Find multiple LFCs in a sounding."""
+    lfc_p_list, lfc_t_list = [], []
+    for i in np.where(intersect_pressures < lcl_pressure)[0]:
+        idp = np.where(pressure < intersect_pressures[i])[0][0]
+        env_lapse = (temperature[idp] - temperature[idp - 1]) / (pressure[idp]
+                                                                 - pressure[idp - 1])
+        prof_lapse = (profile_temperatures[idp] - profile_temperatures[idp - 1]) / \
+                     (pressure[idp] - pressure[idp - 1])
+        if env_lapse > prof_lapse:
+            lfc_p_list.append(intersect_pressures[i])
+            lfc_t_list.append(intersect_temperatures[i])
+    if return_all is True:
+        return lfc_p_list, lfc_t_list
+    else:
+        return lfc_p_list[-1], lfc_t_list[-1]
 
 
 @exporter.export