New version: gsw_geo_strf_dyn_height_1.c

[efiring]

New function: gsw_geo_strf_dyn_height_1.c

This is a complete re-implementation of the basic geostrophic
streamfunction calculation.  It should replace gsw_geo_strf_dyn_height.c,
or the latter should be replaced with a shim that calls the new
version.

This function differs from gsw_geo_strf_dyn_height in several ways:

1) Most importantly, to limited the extent that I have tested it,
it does not leak memory, use uninitialized memory, or otherwise
behave badly, unlike the original.

2) Interpolation can be selected as linear or PCHIP, and the
dp threshold for interpolation is an input argument.

Note that the original (and ancient) RR68 algorithm, for which
I think the implementation is buggy, is no longer used in
GSW-Matlab.  I don't think it is worth the effort it would take
to debug it.

3) It requires the reference pressure to be within the range of the
profile; the original allowed it to be less than the minimum pressure,
which has no use case that I can see, but added complexity to the code.

4) It returns an integer error code instead of either the input
dynamic height pointer or NULL.  This makes it easy for the calling
routine to safely free memory in the event of an error.

[efiring]

This is intended as a step toward solving https://github.com/TEOS-10/GSW-Python/issues/21 and #6. Obviously this can be squashed down to get rid of my earlier attempts to diagnose and fix the original version. Whether or when to do that may depend on how we decide to handle the original version.

If this re-implementation didn't change the API--the function's return type and value--then it could simply replace the original, although I suspect that would lead to somewhat different numerical results and consequent failures of the Matlab-derived test code. I think that is using a rather coarse sample grid, so the interpolation method will matter.

Or we could leave the old, broken code in place, and just choose to use the new version in GSW-Python, and, I would hope, in GSW-R (since the old code really is a minefield).

I consider the API change to be essential for decent library behavior; returning a pointer to the original array, or NULL upon error, seems to me to be a poor strategy.

[PaulMBarker]

I am finalising a variation to the interpolation method that is currently implemented in GSW-matlab. It does not require a stable cast

[efiring]

Also: I am aware that I chose a different code style for the new function; if others are adamant that something closer to the original style should be maintained, I could do it, but I strongly prefer the style I chose because, at some cost in loss of compactness, it makes it easier to see matching braces and corresponding structure.

[dankelley]

Re interpolation: in the oce package, we provide several ways to interpolate. One is Reiniger and Ross, but another is the related method promoted subsequently by UNESCO. Below are the citations that are provided by help(oceApprox, package="oce"). According to that documentation, the method is defined in references 2 and 3, but I've not gone back to check on the details (and I don't remember much, since git blame tells me that the underlying code for that goes back to 2009).

1. R.F. Reiniger and C.K. Ross, 1968. A method of interpolation with application to oceanographic data. Deep Sea Research, 15, 185-193.

2. Daphne R. Johnson, Tim P. Boyer, Hernan E. Garcia, Ricardo A. Locarnini, Olga K. Baranova, and Melissa M. Zweng, 2011. World Ocean Database 2009 Documentation. NODC Internal report 20. Ocean Climate Laboratory, National Oceanographic Data Center. Silver Spring, Maryland.

3. UNESCO, 1991. Processing of oceanographic station data, 138 pp., Imprimerie des Presses Universitaires de France, United Nations Educational, Scientific and Cultural Organization, France.

[dankelley]

I copied the new gsw_oceanographic_toolbox.c code (from branch efiring-strf_new, commit b4ed88d11df00d6124e5cdc09d5e4f28d74765a4) into GSW-R/src (not pushed to git though), built it up, and ran the test code from the R documentation (below), which uses check values as listed on the TEOS-10. The old check values still work. I'm not quite sure how that can be, with so many digits being checked and with what I am inferring is a different interpolation scheme. I wonder if this means I've done something wrong

> library(gsw)
> SA <- c(34.7118, 34.8915, 35.0256, 34.8472, 34.7366, 34.7324)
> CT <- c(28.8099, 28.4392, 22.7862, 10.2262,  6.8272,  4.3236)
> p <- c(      10,      50,     125,     250,     600,    1000)
> p_ref <- 1000
> dh <- gsw_geo_strf_dyn_height(SA, CT, p, p_ref)
> expect_equal(dh, c(17.039204557769487, 14.665853784722286, 10.912861136923812,
+                    7.567928838774945, 3.393524055565328, 0))
> SA <- c(34.7118, 34.8915, 35.0256, 34.8472, 34.7366, 34.7324)
> CT <- c(28.8099, 28.4392, 22.7862, 10.2262,  6.8272,  4.3236)
> delta_p <- c(10,      40,      75,     125,     350,     400)
> r <- gsw_geo_strf_dyn_height_pc(SA, CT, delta_p)
> expect_equal(r$dyn_height, c(-0.300346215853487, -1.755165998114308, -4.423531083131365,
+                              -6.816659136254657, -9.453175257818430, -12.721009624991439))
> expect_equal(r$p_mid/1e2, c(0.050000000000000, 0.300000000000000, 0.875000000000000,
+                             1.875000000000000, 4.250000000000000, 8.000000000000000))
> 
> 
> SA <- c(34.7118, 34.8915, 35.0256, 34.8472, 34.7366, 34.7324)
> CT <- c(28.8099, 28.4392, 22.7862, 10.2262,  6.8272,  4.3236)
> p <- c(      10,      50,     125,     250,     600,    1000)
> p_ref <- 1000
> dh <- gsw_geo_strf_dyn_height(SA, CT, p, p_ref)
> expect_equal(dh, c(17.039204557769487, 14.665853784722286, 10.912861136923812,
+                    7.567928838774945, 3.393524055565328, 0))

[efiring]

@dankelley Thanks for the prompt response and attempt to incorporate the new code. The reason you didn't see any change in the check values is that the changes I made are half-baked; instead of modifying gsw_geo_strf_dyn_height, as I originally started to do, I left it in place, and made a new function, gsw_geo_strf_dyn_height_1. So I suspect your R is using the old function, not the new one. Similarly, I made no attempt to run the C check code on the new function. I have been testing it only with the sanity checker (and memory checker, "valgrind") that I added for that purpose, and that is only built when using debug_makefile. More testing is certainly needed, but I wanted to get the "work in progress" exposed for review before going too far.

[efiring]

I have made some more changes, then squashed and force-pushed those involving the new function. I think this is ready-to-use, although of course much more testing, including corner cases, would be nice. Testing at the C level is laborious; it will be easier to incorporate this in the Python (or R, for the R people) and test there.

One more step that could be taken in this PR or in a follow-up would be to replace the current gsw_geo_strf_dyn_height.c with a simple shim that runs the new version. This would be for the sake of compatibility for people using the original, but I would recommend deprecating it. The new version's arguments and return value make more sense than the original.

I expect that this shimmed version will get closer to the original test values with PCHIP than with linear interpolation, but it will not reproduce them closely enough for the test to pass; hence it will be necessary to use the obtained values as the new test values.

Discussion of such strategies, and help with testing, is invited.

[efiring]

I'm closing this; better to go straight to #23.