diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py
index 2bb67054b6d..41b36bec3d8 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 fae1a08d1ce..4114ded9dd6 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,27 @@ 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
         else:
-            x = x[idx]
-            y = y[idx]
-            return x[0], y[0]
+            lfc_p_list, lfc_t_list = [], []
+            for i in np.where(x < this_lcl[0])[0]:
+                idp = np.where(pressure < x[i])[0][0]
+                env_lapse = (temperature[idp] - temperature[idp - 1]) / (pressure[idp]
+                                                                         - pressure[idp - 1])
+                prof_lapse = (parcel_temperature_profile[idp]
+                              - parcel_temperature_profile[idp - 1]) / (pressure[idp]
+                                                                        - pressure[idp - 1])
+                # An intersection is a LFC if the environmental lapse rate is greater than
+                # the profile lapse rate
+                if env_lapse > prof_lapse:
+                    lfc_p_list.append(x[i])
+                    lfc_t_list.append(y[i])
+            # Return the lowest-pressure LFC if return_all=False, or return all if True
+            if return_all is False:
+                x, y = lfc_p_list[-1], lfc_t_list[-1]
+            else:
+                x, y = lfc_p_list, lfc_t_list
+            return x, y
 
 
 @exporter.export