Repository: NCAR/geocat-comp

Issue:
geocat.comp.meteorology.delta_pressure is different with dpres_plev
###this is the ncl  code
plev = (/ 1000.,950.,900.,850.,800.,750.,700.,650.,600., \
    550.,500.,450.,400.,350.,300.,250.,200., \
    175.,150.,125.,100., 80., 70., 60., 50., \
     40., 30., 25., 20. /)


plev       = plev*100.        
plev@units = "Pa"        ; Pa   =  kg/(m s2) (Pascal)


ptop       = min(plev)    

psfc       = 850*100.           
psfc@units = "Pa" 

dp   = dpres_plevel(plev, psfc, ptop, 0)  ; dp(30)
print(dp)





##this is the python code:
import geocat.comp as gc
import numpy as np 
plev = [ 1000.,950.,900.,850.,800.,750.,700.,650.,600., \
    550.,500.,450.,400.,350.,300.,250.,200., \
    175.,150.,125.,100., 80., 70., 60., 50., \
     40., 30., 25., 20. ]
plev       = np.array(plev)*100.        
ptop       = min(plev)    
psfc       = 850*100.           
dp   = gc.meteorology.delta_pressure(plev, psfc)  
print(dp)

but the result is different if the psfc is lower than the first value of plev
![屏幕截图 2024-01-05 162143](https://github.com/NCAR/geocat-comp/assets/98377186/649ccf69-99cc-4a99-9b33-8d182814cc21)
![屏幕截图 2024-01-05 162521](https://github.com/NCAR/geocat-comp/assets/98377186/fe282ef1-360f-49d1-b9b4-01c36dd0339a)


Discussion:
kafitzgerald: Rolling this back to draft status while I add in some better test coverage.

Patch:
diff --git a/docs/release-notes.rst b/docs/release-notes.rst
index b28bdc2..474f355 100644
--- a/docs/release-notes.rst
+++ b/docs/release-notes.rst
@@ -5,7 +5,7 @@
 Release Notes
 =============
 
-v2024.03.0 (Unreleased)
+v2024.04.0 (Unreleased)
 -----------------------
 
 Internal Changes
@@ -15,6 +15,11 @@ Internal Changes
 * Add M1 runners to CI by `Katelyn FitzGerald`_ in (:pr:`581`)
 * Reorganize dask compatibility tests by `Anissa Zacharias`_ in (:pr:`568`)
 
+Bug Fixes
+^^^^^^^^^
+* Fix ``delta_pressure`` to handle the case where pressure level(s) are greater than surface pressure by `Katelyn FitzGerald`_ in (:pr:`571`)
+
+
 v2024.02.0 (February 28, 2024)
 ------------------------------
 This release switches the package to use an implicit namespace and adds support
diff --git a/geocat/comp/meteorology.py b/geocat/comp/meteorology.py
index e4a7a52..1bf12dc 100644
--- a/geocat/comp/meteorology.py
+++ b/geocat/comp/meteorology.py
@@ -1401,13 +1401,17 @@ def _delta_pressure1D(pressure_lev, surface_pressure):
         The pressure layer thickness array. Shares dimensions and units of
         `pressure_lev`.
     """
+
+    # Leaving this here to lay groundwork for including it as a separate argument
     pressure_top = min(pressure_lev)
 
     # Safety checks
-    if pressure_top <= 0:
-        warnings.warn("'pressure_lev` values must all be positive.")
+    if pressure_top < 0:
+        raise ValueError(
+            "`pressure_lev` values must all be greater than or equal to 0.")
+
     if pressure_top > surface_pressure:
-        warnings.warn(
+        raise ValueError(
             "`surface_pressure` must be greater than minimum `pressure_lev` value."
         )
 
@@ -1417,15 +1421,25 @@ def _delta_pressure1D(pressure_lev, surface_pressure):
         pressure_lev = np.flip(pressure_lev)
 
     # Calculate delta pressure
-    delta_pressure = np.empty_like(pressure_lev)
+    delta_pressure = np.full_like(pressure_lev, np.nan)
+
+    [indices] = np.nonzero(np.array(pressure_lev) <= surface_pressure)
 
-    delta_pressure[0] = (pressure_lev[0] +
-                         pressure_lev[1]) / 2 - pressure_top  # top level
-    delta_pressure[1:-1] = [
-        (a - b) / 2 for a, b in zip(pressure_lev[2:], pressure_lev[:-1])
+    start_level = min(indices)
+    end_level = max(indices)
+
+    delta_pressure[start_level] = (pressure_lev[start_level] + pressure_lev[
+        start_level + 1]) / 2 - pressure_top  # top level
+
+    delta_pressure[start_level + 1:end_level] = [
+        (a - b) / 2
+        for a, b in zip(pressure_lev[start_level + 2:end_level +
+                                     1], pressure_lev[start_level:end_level])
     ]
-    delta_pressure[-1] = surface_pressure - (
-        pressure_lev[-1] + pressure_lev[-2]) / 2  # bottom level
+
+    delta_pressure[end_level] = surface_pressure - (
+        pressure_lev[end_level] +
+        pressure_lev[end_level - 1]) / 2  # bottom level
 
     # Return delta_pressure to original order
     if is_pressuredecreasing:


Test patch:
diff --git a/test/test_meteorology.py b/test/test_meteorology.py
index d4ed8ab..9315751 100644
--- a/test/test_meteorology.py
+++ b/test/test_meteorology.py
@@ -515,15 +515,29 @@ class Test_Delta_Pressure:
         delta_p = delta_pressure(pressure_lev, self.surface_pressure_scalar)
         assert sum(delta_p) == (self.surface_pressure_scalar - pressure_top)
 
-    def test_negative_pressure_warning(self) -> None:
+    def test_delta_pressure_level_below_surface(self) -> None:
+        pressure_lev = [float(i) for i in self.pressure_lev]
+        surface_pressure_adjusted = 900.0
+        delta_p = delta_pressure(pressure_lev, surface_pressure_adjusted)
+        assert np.nansum(delta_p) == (surface_pressure_adjusted -
+                                      min(pressure_lev))
+
+    def test_delta_pressure_levels_below_surface(self) -> None:
+        pressure_lev = [float(i) for i in self.pressure_lev]
+        surface_pressure_adjusted = 50.0
+        delta_p = delta_pressure(pressure_lev, surface_pressure_adjusted)
+        assert np.nansum(delta_p) == (surface_pressure_adjusted -
+                                      min(pressure_lev))
+
+    def test_negative_pressure_error(self) -> None:
         pressure_lev_negative = self.pressure_lev.copy()
         pressure_lev_negative[0] = -5
-        with pytest.warns(UserWarning):
+        with pytest.raises(ValueError):
             delta_pressure(pressure_lev_negative, self.surface_pressure_scalar)
 
-    def test_relative_pressure_warning(self) -> None:
+    def test_relative_pressure_error(self) -> None:
         surface_pressure_low = 0.5
-        with pytest.warns(UserWarning):
+        with pytest.raises(ValueError):
             delta_pressure(self.pressure_lev, surface_pressure_low)
 
     def test_output_type(self) -> None: