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..6b23c01
--- /dev/null
+++ b/gsw/tests/test_fixups.py
@@ -0,0 +1,35 @@
+"""
+Tests of function modifications by _fixed_wrapped_ufuncs.py
+"""
+
+import numpy as np
+import pytest
+
+import gsw
+
+zvals_ok = [
+    -5500,
+    np.linspace(-100, 0, 11),
+    np.nan,
+    [np.nan, -100],
+    np.ma.masked_invalid([np.nan, -100]),
+]
+
+zvals_bad = [
+    5500,
+    np.linspace(0, 100, 11),
+    [np.nan, 100],
+    np.ma.masked_invalid([np.nan, 100]),
+]
+
+
+@pytest.mark.parametrize("z", zvals_ok)
+def test_p_from_z_ok(z):
+    # smoke test: doesn't raise an exception
+    gsw.p_from_z(z, 30)
+
+
+@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