Fix constant integration offset in geo_strf_dyn_height

[romeroqe]

This pull request ensures that geo_strf_dyn_height returns a dynamic height anomaly that is exactly zero at the reference pressure (p_ref).

In the MATLAB GSW toolbox, the dynamic height is explicitly anchored by subtracting its value at p_ref, ensuring:

dh(p_ref) = 0

In the current Python implementation, a constant integration offset may remain due to the implicit integration constant of the numerical integral. As a result, the dynamic height at p_ref can differ from zero even when p_ref is included in the integration grid.

This modification removes that constant by subtracting the dynamic height evaluated at p_ref within each processed profile.

This prevents discrepancies between MATLAB and Python implementations when comparing dynamic height or derived quantities.

[DocOtak]

Can you provide an example where this doesn't return 0 where p==p_ref? The python method just calls the underlying c-function which appears to be handling this situation already.

[romeroqe]

I am currently performing a global computation using DMQC Argo profiles and detected that this behavior occurs when p_ref is equal to the deepest pressure level in the profile. In these cases, even though p_ref is explicitly part of the pressure grid, geo_strf_dyn_height returns a non-zero value at p == p_ref, indicating a constant integration offset.

So far, I have not identified additional scenarios where this occurs. However, in the MATLAB GSW toolbox the dynamic height is explicitly anchored by subtracting its value at p_ref immediately before returning the result. This anchoring is applied both in the “perfect resolution” branch:

geo_strf_dyn_height = geo_strf_dyn_height0 - geo_strf_dyn_height_p_ref;

and in the interpolated branch:

geo_strf_dyn_height(Ibdata,Iprofile) = D0_i(Iidata) - D0_i(Ibpr);

In both cases, the value at p_ref is explicitly removed, guaranteeing dh(p_ref) = 0.

Applying the same anchoring step in the Python implementation resolves the deepest-level case and ensures consistency with MATLAB. Importantly, this modification does not affect profiles where dh(p_ref) is already zero — in those cases the subtraction has no effect.

Below I share a figure showing three DMQC Argo profiles where p_ref corresponds to the deepest pressure level. The dynamic height was computed using:

files = [
	"aoml/1900047/profiles/D1900047_062.nc",
	"aoml/1900047/profiles/D1900047_095.nc",
	"coriolis/6904086/profiles/D6904086_032.nc",
]

fig, axs = plt.subplots(1, 3, figsize=(17, 6))

for k in range(3):
  profile = getProfileDataFromArgoNc(files[k])
  profile = (profile.sort_values("pres").drop_duplicates(subset="pres", keep="first").reset_index(drop=True))
  
  file = files[k].split("/")[-1]
  lat = np.round(profile.latitude.values[0], 2)
  lon = np.round(profile.longitude.values[0], 2)
  p = profile.pres.values
  SA = profile.asal.values
  CT = profile.ctemp.values
  p_ref = profile.pres.values[-1]
  
  iref = np.where(np.isclose(p, p_ref))[0][0]
  gsdh = gsw.geo_strf_dyn_height(SA, CT, p, p_ref, interp_method="mrst")

  axs[k].plot(gsdh, p, label=f"Original (dh(p_ref) = {gsdh[iref]:.3f} m² s⁻²)")
  axs[k].plot(gsdh-gsdh[iref], p, label=f"Anchored (dh(p_ref) = {(gsdh-gsdh[iref])[iref]:.3f} m² s⁻²)")
  axs[k].axhline(y=p[iref], ls="--", lw=1, c="red", label=f"p_ref = {p_ref:.3f}")
  axs[k].invert_yaxis()
  axs[k].legend(loc="upper left")
  axs[k].set_title(f"{file} [lat: {lat}; lon: {lon}]")
		
fig.savefig('gsdh.png', bbox_inches="tight", dpi=300)

[efiring]

I think this does need to be solved in the C, not in the Python. The C is using exact equality tests for floating point numbers. This could cause a problem like this if one side of a test is promoted from single precision, for example.

We need a simple easy-to-use test case. @romeroqe, please make a complete self-contained example we can use to investigate this, and later in an actual test. Maybe save all the inputs to gsw.geo_strf_dyn_height from the first profile in a single netcdf file, preserving the types of all variables, so that a single read of that file followed by a single call to gsw.geo_strf_dyn_height will reproduce the error.