Repository: NCAR/geocat-comp (python)

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)


Pull request description:
Fix delta_pressure to handle the case when pressure levels are greater than the surface pressure
## PR Summary
Adds some code to handle the case when pressure levels are greater than the surface pressure.

I'm not necessarily convinced this is a the best / most performant way to handle this, but should do the job and we'll likely want to revisit for adding the pressure top (which will require even more logic to handle well).  
 
Closes #538

## PR Checklist
<!-- Please mark any checkboxes that do not apply to this PR as [N/A]. If an entire section doesn't
apply to this PR, comment it out or delete it. -->
**General**
- [x] Make an issue if one doesn't already exist
- [x] Link the issue this PR resolves by adding `closes #XXX` to the PR description where XXX is the number of the issue.
- [x] Add a brief summary of changes to `docs/release-notes.rst` in a relevant section for the next unreleased release. Possible sections include: Documentation, New Features, Bug Fixes, Internal Changes, Breaking Changes/Deprecated
- [x] Add appropriate labels to this PR
- [x] Make your changes in a forked repository rather than directly in this repo
- [x] Open this PR as a draft if it is not ready for review
- [x] Convert this PR from a draft to a full PR before requesting reviewers
- [x] Passes `precommit`. To set up on your local, run `pre-commit install` from the top level of the repository. To manually run pre-commits, use `pre-commit run --all-files` and re-add any changed files before committing again and pushing.
- [ ] If needed, squash and merge PR commits into a single commit to clean up commit history

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: