Merge pull request #1200 from zbruick/lcl_probs

Add support for grids when surface is LCL

src/metpy/calc/thermo.py

@@ -310,7 +310,8 @@ def lcl(pressure, temperature, dewpt, max_iters=50, eps=1e-5):
     r"""Calculate the lifted condensation level (LCL) using from the starting point.
 
     The starting state for the parcel is defined by `temperature`, `dewpt`,
-    and `pressure`.
+    and `pressure`. If these are arrays, this function will return a LCL
+    for every index. This function does work with surface grids as a result.
 
     Parameters
     ----------
@@ -356,14 +357,12 @@ def _lcl_iter(p, p0, w, t):
         return (p0 * (td / t) ** (1. / mpconsts.kappa)).m
 
     w = mixing_ratio(saturation_vapor_pressure(dewpt), pressure)
-    fp = so.fixed_point(_lcl_iter, pressure.m, args=(pressure.m, w, temperature),
-                        xtol=eps, maxiter=max_iters)
-    lcl_p = fp * pressure.units
+    lcl_p = so.fixed_point(_lcl_iter, pressure.m, args=(pressure.m, w, temperature),
+                           xtol=eps, maxiter=max_iters)
 
-    # Conditional needed due to precision error with np.log in dewpoint.
+    # np.isclose needed if surface is LCL due to precision error with np.log in dewpoint.
     # Causes issues with parcel_profile_with_lcl if removed. Issue #1187
-    if np.isclose(lcl_p, pressure):
-        lcl_p = pressure
+    lcl_p = np.where(np.isclose(lcl_p, pressure), pressure, lcl_p) * pressure.units
 
     return lcl_p, dewpoint(vapor_pressure(lcl_p, w)).to(temperature.units)
 

tests/calc/test_thermo.py

@@ -1415,3 +1415,15 @@ def test_lcl_convergence_issue():
     dewpoint = np.array([14.4, 11.7, 8.2, 7.8, 7.6]) * units.degC
     lcl_pressure, _ = lcl(pressure[0], temperature[0], dewpoint[0])
     assert_almost_equal(lcl_pressure, 990 * units.hPa, 0)
+
+
+def test_lcl_grid_surface_LCLs():
+    """Test surface grid where some values have LCLs at the surface."""
+    pressure = np.array([1000, 990, 1010]) * units.hPa
+    temperature = np.array([15, 14, 13]) * units.degC
+    dewpoint = np.array([15, 10, 13]) * units.degC
+    lcl_pressure, lcl_temperature = lcl(pressure, temperature, dewpoint)
+    pres_truth = np.array([1000, 932.1515324, 1010]) * units.hPa
+    temp_truth = np.array([15, 9.10391763, 13]) * units.degC
+    assert_array_almost_equal(lcl_pressure, pres_truth, 7)
+    assert_array_almost_equal(lcl_temperature, temp_truth, 7)