diff --git a/gsw/geostrophy.py b/gsw/geostrophy.py
index 501f01e..ac0c2dc 100644
--- a/gsw/geostrophy.py
+++ b/gsw/geostrophy.py
@@ -9,6 +9,7 @@
 from .conversions import z_from_p
 
 __all__ = ['geo_strf_dyn_height',
+           'geo_strf_steric_height',
            'distance',
            'f',
            'geostrophic_velocity',
@@ -105,6 +106,52 @@ def geo_strf_dyn_height(SA, CT, p, p_ref=0, axis=0, max_dp=1.0,
     return dh
 
 
+@match_args_return
+def geo_strf_steric_height(SA, CT, p, p_ref=0, axis=0, max_dp=1.0, interp_method="pchip"):
+    """
+    Steric height anomaly as a function of pressure.
+
+    Parameters
+    ----------
+    SA : array-like
+        Absolute Salinity, g/kg
+    CT : array-like
+        Conservative Temperature (ITS-90), degrees C
+    p : array-like
+        Sea pressure (absolute pressure minus 10.1325 dbar), dbar
+    p_ref : float or array-like, optional
+        Reference pressure, dbar
+    axis : int, optional, default is 0
+        The index of the pressure dimension in SA and CT.
+    max_dp : float
+        If any pressure interval in the input p exceeds max_dp, the dynamic
+        height will be calculated after interpolating to a grid with this
+        spacing.
+    interp_method : string {'mrst', 'pchip', 'linear'}
+        Interpolation algorithm.
+
+    Returns
+    -------
+    steric_height : array
+        This is the integral of specific volume anomaly with respect
+        to pressure, from each pressure in p to the specified
+        reference pressure, divided by the global mean surface value of
+        gravitational acceleration, 9.7963 m s^-2. (see page 46 of Griffies, 2004) 
+        It is not exactly the height of an isobaric surface above a
+        geopotential surface, but is an exact geostrophic streamfunction.
+
+    References
+    ----------
+    Griffies, S. M., 2004: Fundamentals of Ocean Climate Models. Princeton, 
+    NJ: Princeton University Press, 518 pp + xxxiv.
+
+    """
+    return (
+        geo_strf_dyn_height(SA, CT, p, p_ref, axis=axis, max_dp=max_dp, interp_method=interp_method)
+        / 9.7963
+    )
+
+
 def unwrap(lon, centered=True, copy=True):
     """
     Unwrap a sequence of longitudes or headings in degrees.
diff --git a/gsw/tests/geo_strf_steric_height.npy b/gsw/tests/geo_strf_steric_height.npy
new file mode 100644
index 0000000..be24cf0
Binary files /dev/null and b/gsw/tests/geo_strf_steric_height.npy differ
diff --git a/gsw/tests/test_check_functions.py b/gsw/tests/test_check_functions.py
index f9166c7..2cf0c6c 100644
--- a/gsw/tests/test_check_functions.py
+++ b/gsw/tests/test_check_functions.py
@@ -33,6 +33,7 @@
 
 # Substitute new check values for the pchip interpolation version.
 cv.geo_strf_dyn_height = np.load(os.path.join(root_path,'geo_strf_dyn_height.npy'))
+cv.geo_strf_steric_height = np.load(os.path.join(root_path,'geo_strf_steric_height.npy'))
 cv.geo_strf_velocity = np.load(os.path.join(root_path,'geo_strf_velocity.npy'))
 
 cf = Bunch()
diff --git a/gsw/tests/test_geostrophy.py b/gsw/tests/test_geostrophy.py
index e0d4423..8d69bb6 100644
--- a/gsw/tests/test_geostrophy.py
+++ b/gsw/tests/test_geostrophy.py
@@ -113,3 +113,15 @@ def test_dyn_height_mrst():
     strf = gsw.geo_strf_dyn_height(SA, CT, p, p_ref=pr, interp_method='mrst')
 
     assert_allclose(strf, cv.geo_strf_dyn_height, rtol=0, atol=cv.geo_strf_dyn_height_ca)
+
+def test_steric_height_mrst():
+    """
+    Tests the MRST-PCHIP interpolation method.
+    """
+    p = cv.p_chck_cast
+    CT = cv.CT_chck_cast
+    SA = cv.SA_chck_cast
+    pr = cv.pr
+    strf = gsw.geo_strf_steric_height(SA, CT, p, p_ref=pr, interp_method='mrst')
+
+    assert_allclose(strf, cv.geo_strf_steric_height, rtol=0, atol=cv.geo_strf_steric_height_ca)
\ No newline at end of file
diff --git a/gsw/tests/test_xarray.py b/gsw/tests/test_xarray.py
index 431fd75..2a76e3a 100644
--- a/gsw/tests/test_xarray.py
+++ b/gsw/tests/test_xarray.py
@@ -42,6 +42,7 @@
 
 # Substitute new check values for the pchip interpolation version.
 cv.geo_strf_dyn_height = np.load(os.path.join(root_path,'geo_strf_dyn_height.npy'))
+cv.geo_strf_steric_height = np.load(os.path.join(root_path,'geo_strf_steric_height.npy'))
 cv.geo_strf_velocity = np.load(os.path.join(root_path,'geo_strf_velocity.npy'))
 
 for name in ['SA_chck_cast', 't_chck_cast', 'p_chck_cast']:
diff --git a/gsw/tests/write_geo_npyfiles.py b/gsw/tests/write_geo_npyfiles.py
index d91e72d..d816b37 100644
--- a/gsw/tests/write_geo_npyfiles.py
+++ b/gsw/tests/write_geo_npyfiles.py
@@ -25,6 +25,12 @@
                                      cv.pr)
 np.save('geo_strf_dyn_height.npy', dyn_height)
 
+steric_height = gsw.geo_strf_steric_height(cv.SA_chck_cast,
+                                     cv.CT_chck_cast,
+                                     cv.p_chck_cast,
+                                     cv.pr)
+np.save('geo_strf_steric_height.npy', steric_height)
+
 lon = cv.long_chck_cast
 lat = cv.lat_chck_cast
 p = cv.p_chck_cast