Repository: Unidata/MetPy

Issue:
moist_lapse fails for increasing pressure arrays whose first element is the reference pressure
I'm using Mac, Python 3.8.5 and MetPy 1.1.0, and the following code
```python
import metpy.calc as mpcalc
from metpy.units import units

mpcalc.moist_lapse([600, 700]*units.mbar, 0*units.celsius)
```
produces a `ValueError` with the following traceback:
<details>

```python-traceback
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-65-8d4bc4174249> in <module>
----> 1 mpcalc.moist_lapse([600, 700]*units.mbar, 0*units.celsius)

~/anaconda3/lib/python3.8/site-packages/metpy/xarray.py in wrapper(*args, **kwargs)
   1214 
   1215             # Evaluate inner calculation
-> 1216             result = func(*bound_args.args, **bound_args.kwargs)
   1217 
   1218             # Wrap output based on match and match_unit

~/anaconda3/lib/python3.8/site-packages/metpy/units.py in wrapper(*args, **kwargs)
    244                                'that the function is being called properly.\n') + msg
    245                 raise ValueError(msg)
--> 246             return func(*args, **kwargs)
    247 
    248         return wrapper

~/anaconda3/lib/python3.8/site-packages/metpy/calc/thermo.py in moist_lapse(pressure, temperature, reference_pressure)
    332         # Integrate downward in pressure
    333         pres_down = np.append(reference_pressure.m, pressure[(ref_pres_idx - 1)::-1].m)
--> 334         trace_down = si.odeint(dt, temperature.m.squeeze(), pres_down.squeeze())
    335         ret_temperatures = np.concatenate((ret_temperatures, trace_down[:0:-1]))
    336 

~/anaconda3/lib/python3.8/site-packages/scipy/integrate/odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
    233     dt = np.diff(t)
    234     if not((dt >= 0).all() or (dt <= 0).all()):
--> 235         raise ValueError("The values in t must be monotonically increasing "
    236                          "or monotonically decreasing; repeated values are "
    237                          "allowed.")

ValueError: The values in t must be monotonically increasing or monotonically decreasing; repeated values are allowed.
```

</details>

when I would expect to get `[0.0 6.47748352956819] degree_Celsius`. My current workaround is quite cumbersome:
```python
import metpy.calc as mpcalc
from metpy.units import units
from metpy.units import concatenate

concatenate([0*units.celsius, mpcalc.moist_lapse(700*units.mbar, 0*units.celsius, reference_pressure=600*units.mbar)])
```

Looking at [lines 307-335](https://github.com/Unidata/MetPy/blob/3e084846414ffe63162951cacd4e60a7d6881fe1/src/metpy/calc/thermo.py#L307-L335) of thermo.py, 
<details>

```python
pressure = np.atleast_1d(pressure)
if reference_pressure is None:
    reference_pressure = pressure[0]

if np.isnan(reference_pressure):
    return units.Quantity(np.full(pressure.shape, np.nan), temperature.units)

pressure = pressure.to('mbar')
reference_pressure = reference_pressure.to('mbar')
org_units = temperature.units
temperature = np.atleast_1d(temperature).to('kelvin')

side = 'left'

pres_decreasing = (pressure[0] > pressure[-1])
if pres_decreasing:
    # Everything is easier if pressures are in increasing order
    pressure = pressure[::-1]
    side = 'right'

ref_pres_idx = np.searchsorted(pressure.m, reference_pressure.m, side=side)

ret_temperatures = np.empty((0, temperature.shape[0]))

if _greater_or_close(reference_pressure, pressure.min()):
    # Integrate downward in pressure
    pres_down = np.append(reference_pressure.m, pressure[(ref_pres_idx - 1)::-1].m)
    trace_down = si.odeint(dt, temperature.m.squeeze(), pres_down.squeeze())
    ret_temperatures = np.concatenate((ret_temperatures, trace_down[:0:-1]))
```
</details>

this is what I think is happening in my above example:
- `reference_pressure` is set to 600 mbar,
- `pressure` is increasing to start with, so `ref_pres_idx` is set to `0`,
- `pressure[(ref_pres_idx - 1)::-1].m` evaluates to `array([700, 600])`,
- `pres_down` is set to `array([600, 700, 600])`,
- which when passed to `si.odeint` causes the `ValueError` because it's not monotonic.

Discussion:

Patch:
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index c433f63a3b..bbc7a34e84 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -316,16 +316,12 @@ def moist_lapse(pressure, temperature, reference_pressure=None):
     org_units = temperature.units
     temperature = np.atleast_1d(temperature).to('kelvin')
 
-    side = 'left'
-
     pres_decreasing = (pressure[0] > pressure[-1])
     if pres_decreasing:
         # Everything is easier if pressures are in increasing order
         pressure = pressure[::-1]
-        side = 'right'
-
-    ref_pres_idx = np.searchsorted(pressure.m, reference_pressure.m, side=side)
 
+    ref_pres_idx = np.searchsorted(pressure.m, reference_pressure.m, side='right')
     ret_temperatures = np.empty((0, temperature.shape[0]))
 
     if _greater_or_close(reference_pressure, pressure.min()):


Test patch:
diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py
index 687f15df12..36b50bb9f1 100644
--- a/tests/calc/test_thermo.py
+++ b/tests/calc/test_thermo.py
@@ -155,6 +155,23 @@ def test_moist_lapse_nan_ref_press():
     assert_nan(temp, units.degC)
 
 
+def test_moist_lapse_downwards():
+    """Test moist_lapse when integrating downwards (#2128)."""
+    temp = moist_lapse(units.Quantity([600, 700], 'mbar'), units.Quantity(0, 'degC'))
+    assert_almost_equal(temp, units.Quantity([0, 6.47748353], units.degC))
+
+
+@pytest.mark.parametrize('direction', (1, -1))
+@pytest.mark.parametrize('start', list(range(5)))
+def test_moist_lapse_starting_points(start, direction):
+    """Test moist_lapse with a variety of reference points."""
+    truth = units.Quantity([20.0804315, 17.2333509, 14.0752659, 6.4774835, 0.0],
+                           'degC')[::direction]
+    pressure = units.Quantity([1000, 925, 850, 700, 600], 'hPa')[::direction]
+    temp = moist_lapse(pressure, truth[start], pressure[start])
+    assert_almost_equal(temp, truth, 4)
+
+
 def test_parcel_profile():
     """Test parcel profile calculation."""
     levels = np.array([1000., 900., 800., 700., 600., 500., 400.]) * units.mbar