Fix up Showalter Index calculation

src/metpy/calc/thermo.py

@@ -3314,27 +3314,26 @@ def gradient_richardson_number(height, potential_temperature, u, v, vertical_dim
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]', '[temperature]', '[temperature]')
-def showalter_index(pressure, temperature, dewpt):
-    """Calculate Showalter Index from pressure temperature and 850 hPa lcl.
+def showalter_index(pressure, temperature, dewpoint):
+    """Calculate Showalter Index.
 
     Showalter Index derived from [Galway1956]_:
     SI = T500 - Tp500
 
     where:
     T500 is the measured temperature at 500 hPa
-    Tp500 is the temperature of the lifted parcel at 500 hPa
+    Tp500 is the temperature of the parcel at 500 hPa when lifted from 850 hPa
 
     Parameters
     ----------
         pressure : `pint.Quantity`
-            Atmospheric pressure level(s) of interest, in order from highest to
-            lowest pressure
+            Atmospheric pressure, in order from highest to lowest pressure
 
         temperature : `pint.Quantity`
-            Parcel temperature for corresponding pressure
+            Ambient temperature corresponding to ``pressure``
 
-        dewpt : `pint.Quantity`
-            Parcel dew point temperatures for corresponding pressure
+        dewpoint : `pint.Quantity`
+            Ambient dew point temperatures corresponding to ``pressure``
 
     Returns
     -------
@@ -3343,22 +3342,13 @@ def showalter_index(pressure, temperature, dewpt):
 
     """
     # find the measured temperature and dew point temperature at 850 hPa.
-    t850, td850 = interpolate_1d(units.Quantity(850, 'hPa'), pressure, temperature, dewpt)
+    t850, td850 = interpolate_1d(units.Quantity(850, 'hPa'), pressure, temperature, dewpoint)
 
     # find the parcel profile temperature at 500 hPa.
     tp500 = interpolate_1d(units.Quantity(500, 'hPa'), pressure, temperature)
 
-    # Calculate lcl at the 850 hPa level
-    lcl_calc, _ = lcl(units.Quantity(850, 'hPa'), t850[0], td850[0])
-
-    # Define end height for moist lapse rate calculation
-    p_end = units.Quantity(500, 'hPa')
-
-    # Calculate parcel temp when raised dry adiabatically from surface to lcl
-    dl = dry_lapse(lcl_calc, temperature[0], pressure[0])
-
-    # Calculate parcel temp when raised moist adiabatically from lcl to 500mb
-    ml = moist_lapse(p_end, dl, lcl_calc)
+    # Lift parcel from 850 to 500, handling any needed dry vs. saturated adiabatic processes
+    prof = parcel_profile(units.Quantity([850., 500.], 'hPa'), t850, td850)
 
     # Calculate the Showalter index
-    return tp500 - ml
+    return tp500 - prof[-1]

tests/calc/test_thermo.py

@@ -1912,12 +1912,19 @@ def test_gradient_richardson_number_with_xarray():
 
 def test_showalter_index():
     """Test the Showalter index calculation."""
-    p_upper = np.arange(1000, 200, -50) * units.hPa
-    p_lower = np.arange(175, 0, -25) * units.hPa
-    p = np.append(p_upper, p_lower,)
-    tc = np.linspace(30, -30, 25) * units.degC
-    tdc = np.linspace(10, -30, 25) * units.degC
-
-    result = showalter_index(p, tc, tdc)
-    expected = 23.73036498600014 * units.delta_degree_Celsius
-    assert_array_almost_equal(result, expected, 4)
+    pressure = units.Quantity(np.array([931.0, 925.0, 911.0, 891.0, 886.9, 855.0, 850.0, 825.6,
+                                        796.3, 783.0, 768.0, 759.0, 745.0, 740.4, 733.0, 715.0,
+                                        700.0, 695.0, 687.2, 684.0, 681.0, 677.0, 674.0, 661.9,
+                                        657.0, 639.0, 637.6, 614.0, 592.0, 568.9, 547.4, 526.8,
+                                        500.0, 487.5, 485.0]), 'hPa')
+    temps = units.Quantity(np.array([18.4, 19.8, 20.0, 19.6, 19.3, 16.8, 16.4, 15.1, 13.4,
+                                     12.6, 11.2, 10.4, 8.6, 8.3, 7.8, 5.8, 4.6, 4.2, 3.4, 3.0,
+                                     3.0, 4.4, 5.0, 5.1, 5.2, 3.4, 3.3, 2.4, 1.4, -0.4, -2.2,
+                                     -3.9, -6.3, -7.6, -7.9]), 'degC')
+    dewp = units.Quantity(np.array([9.4, 8.8, 6.0, 8.6, 8.4, 6.8, 6.4, 4.0, 1.0, -0.4, -1.1,
+                                    -1.6, 1.6, -0.2, -3.2, -3.2, -4.4, -2.8, -3.6, -4.0, -6.0,
+                                    -17.6, -25.0, -31.2, -33.8, -29.6, -30.1, -39.0, -47.6,
+                                    -48.9, -50.2, -51.5, -53.3, -55.5, -55.9]), 'degC')
+
+    result = showalter_index(pressure, temps, dewp)
+    assert_almost_equal(result, units.Quantity(7.6024, 'delta_degC'), 4)