diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py index d528c505d71..0efe0a4038f 100644 --- a/metpy/calc/tests/test_thermo.py +++ b/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(): diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py index 662820cd336..48e4536038e 100644 --- a/metpy/calc/thermo.py +++ b/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