From f935d4a75fe01191a7117524c311e96520e69ecd Mon Sep 17 00:00:00 2001
From: Eric Firing <efiring@hawaii.edu>
Date: Mon, 2 Mar 2026 10:29:42 -1000
Subject: [PATCH] p_from_z: raise ValueError if any z > 5

The criterion for raising the exception is taken
from the Matlab version.

Closes #223.
---
 gsw/_fixed_wrapped_ufuncs.py | 24 ++++++++++++---
 gsw/_wrapped_ufuncs.py       |  2 +-
 gsw/tests/test_fixups.py     | 58 ++++++++++++++++++++++++++++++++++++
 tools/docstring_parts.py     |  2 +-
 4 files changed, 80 insertions(+), 6 deletions(-)
 create mode 100644 gsw/tests/test_fixups.py

diff --git a/gsw/_fixed_wrapped_ufuncs.py b/gsw/_fixed_wrapped_ufuncs.py
index 4294f87..5b662ab 100644
--- a/gsw/_fixed_wrapped_ufuncs.py
+++ b/gsw/_fixed_wrapped_ufuncs.py
@@ -8,34 +8,50 @@
 from ._wrapped_ufuncs import *
 
 _p_from_z = p_from_z
+
+
 def p_from_z(z, lat, geo_strf_dyn_height=0, sea_surface_geopotential=0):
+    if numpy.ma.any(numpy.ma.asarray(z) > 5):
+        raise ValueError(f"z should be generally negative; found max(z) = {numpy.ma.max(z)}")
     return _p_from_z(z, lat, geo_strf_dyn_height, sea_surface_geopotential)
+
+
 p_from_z.__doc__ = _p_from_z.__doc__
 
 _z_from_p = z_from_p
+
+
 def z_from_p(p, lat, geo_strf_dyn_height=0, sea_surface_geopotential=0):
     return _z_from_p(p, lat, geo_strf_dyn_height, sea_surface_geopotential)
+
+
 z_from_p.__doc__ = _z_from_p.__doc__
 
 _gibbs = gibbs
+
+
 def gibbs(ns, nt, np, SA, t, p):
     params = {"ns": ns, "nt": nt, "np": np}
     for k, v in params.items():
         u = numpy.unique(v)
         if u.min() < 0 or u.max() > 2 or u.dtype.kind != "i":
-            raise ValueError("ns, nt, np must contain integers 0, 1, or 2;"
-                             f" found {k}={v}")
+            raise ValueError(f"ns, nt, np must contain integers 0, 1, or 2; found {k}={v}")
     return _gibbs(ns, nt, np, SA, t, p)
+
+
 gibbs.__doc__ = _gibbs.__doc__
 
 
 _gibbs_ice = gibbs_ice
+
+
 def gibbs_ice(nt, np, t, p):
     params = {"nt": nt, "np": np}
     for k, v in params.items():
         u = numpy.unique(v)
         if u.min() < 0 or u.max() > 2 or u.dtype.kind != "i":
-            raise ValueError("nt, np must contain integers 0, 1, or 2;"
-                             f" found {k}={v}")
+            raise ValueError(f"nt, np must contain integers 0, 1, or 2; found {k}={v}")
     return _gibbs_ice(nt, np, t, p)
+
+
 gibbs_ice.__doc__ = _gibbs_ice.__doc__
diff --git a/gsw/_wrapped_ufuncs.py b/gsw/_wrapped_ufuncs.py
index e2bc1f9..8ec7a14 100644
--- a/gsw/_wrapped_ufuncs.py
+++ b/gsw/_wrapped_ufuncs.py
@@ -4302,7 +4302,7 @@ def p_from_z(z, lat, geo_strf_dyn_height, sea_surface_geopotential):
     Parameters
     ----------
     z : array-like
-        Depth, positive up, m
+        Height, positive up (so z = -depth), m
     lat : array-like
         Latitude, -90 to 90 degrees
     geo_strf_dyn_height : array-like
diff --git a/gsw/tests/test_fixups.py b/gsw/tests/test_fixups.py
new file mode 100644
index 0000000..88c8413
--- /dev/null
+++ b/gsw/tests/test_fixups.py
@@ -0,0 +1,58 @@
+"""
+Tests of function modifications by _fixed_wrapped_ufuncs.py.
+"""
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+import gsw
+
+# In the following pairs of z, p, the expected values for p are simply the
+# values returned by z_from_p on an M4 Mac with gsw version
+# 3.6.22.dev3+g2bad68c16 (prior to the commit in which this test is added).
+
+zpvals_ok = [
+    (-5500, 5609.875428946266),
+    (
+        np.linspace(-100, 0, 11),
+        np.array(
+            [
+                100.70968879,
+                90.6365065,
+                80.56381587,
+                70.49161697,
+                60.41990988,
+                50.34869469,
+                40.27797148,
+                30.20774032,
+                20.1380013,
+                10.0687545,
+                -0.0,
+            ]
+        ),
+    ),
+    (np.nan, np.nan),
+    ([np.nan, -100], np.array([np.nan, 100.70968879])),
+    (np.ma.masked_invalid([np.nan, -100]), np.ma.masked_invalid([np.nan, 100.70968879])),
+]
+
+zvals_bad = [
+    5500,
+    np.linspace(0, 100, 11),
+    [np.nan, 100],
+    np.ma.masked_invalid([np.nan, 100]),
+]
+
+
+@pytest.mark.parametrize("zp", zpvals_ok)
+def test_p_from_z_ok(zp):
+    z, expected = zp
+    p = gsw.p_from_z(z, 30)
+    assert_allclose(p, expected)
+
+
+@pytest.mark.parametrize("z", zvals_bad)
+def test_p_from_z_bad(z):
+    with pytest.raises(ValueError):
+        gsw.p_from_z(z, 30)
diff --git a/tools/docstring_parts.py b/tools/docstring_parts.py
index 7afea1f..deafe9e 100644
--- a/tools/docstring_parts.py
+++ b/tools/docstring_parts.py
@@ -36,7 +36,7 @@
 pt = "Potential temperature referenced to a sea pressure, degrees C",
 rho = "Seawater density (not anomaly) in-situ, e.g., 1026 kg/m^3.",
 t_Ih = "In-situ temperature of ice (ITS-90), degrees C",
-z = "Depth, positive up, m",
+z = "Height, positive up (so z = -depth), m",
 SA_bulk = "bulk Absolute Salinity of the seawater and ice mixture, g/kg",
 w_Ih =
 """mass fraction of ice: the mass of ice divided by the