Fix reference pressure in geo_strf_dyn_height

gsw_check_functions.c

@@ -132,6 +132,35 @@ void test_infunnel() {
     assert_equal(result, 0, "gsw_infunnel -> invalid salinity");
 }
 
+void test_refine_grid_for_dh() {
+    int nz = 4;
+    double p[]  = {4.9000001 ,   5.69999981,  15.80000019, 995.90002441};
+    double p_ref = 995.90002441;
+    double max_dp_i = 1.0;  // grid delta-p; the problem appeared with *exactly* 1.0...
+    double *p_i;
+    int *p_indices;
+    int ni_max;
+    int ipref;
+    int n_i;
+    int result;
+
+    // Chunk taken from gsw_geo_strf_dyn_height_1.
+    ni_max = nz + (int) ceil((p[nz-1] - p[0]) / max_dp_i) + 2;
+    /* Maximum possible size of new grid: Original grid size plus
+       the number of dp intervals plus 1 for the p_ref,
+       plus 1 so that we can know we exited the loop before we hit it.
+    */
+    p_i = (double *) malloc(ni_max * sizeof(double));
+    p_indices = (int *) malloc(nz * sizeof(int));
+
+    n_i = refine_grid_for_dh(p, p_ref, nz, max_dp_i,
+                             p_i, ni_max,
+                             p_indices, &ipref);
+    // End of chunk from gsw_geo_strf_dyn_height_1
+
+    assert_equal(ipref, n_i - 1, "refine_grid_for_dh -> reference pressure index");
+}
+
 int
 main(int argc, char **argv)
 {
@@ -625,6 +654,7 @@ main(int argc, char **argv)
                 "cti_tracerctinterp",interp_n*interp_m, val7, cti_tracerctinterp);
 
         test_infunnel();
+        test_refine_grid_for_dh();
 
         if (gsw_error_flag)
         {

gsw_oceanographic_toolbox.c

@@ -4057,7 +4057,7 @@ gsw_geo_strf_dyn_height_1(double *sa, double *ct, double *p, double p_ref,
 
     Its return argument is the number of elements used in p_i, or -1 on error.
 */
-static int refine_grid_for_dh(double *p, double p_ref, int nz,
+int refine_grid_for_dh(double *p, double p_ref, int nz,
     double dp,
     double *p_i, int ni_max,  /* size of p_i array; larger than needed */
     int *p_indices, int *p_ref_ind_ptr)
@@ -4094,10 +4094,16 @@ static int refine_grid_for_dh(double *p, double p_ref, int nz,
         /* We did not insert p_ref, so insert either p_next or p[iorig]. */
         if (p_next < p[iorig] - p_tol) {
             p_i[i] = p_next;
+            if (p_ref == p_next) {
+                *p_ref_ind_ptr = i;
+            }
             iuniform++;
         }
         else {
             p_i[i] = p[iorig];
+            if (p_ref == p[iorig]) {
+                *p_ref_ind_ptr = i;
+            }
             p_indices[iorig] = i;
             /* Skip this p_next if it is close to the point we just added. */
             if (p_next < p[iorig] + p_tol) {

gswteos-10.h

@@ -318,6 +318,8 @@ DECLSPEC extern double gsw_z_from_p(double p, double lat, double geo_strf_dyn_he
                 double sea_surface_geopotential);
 DECLSPEC extern double gsw_p_from_z(double z, double lat, double geo_strf_dyn_height,
                 double sea_surface_geopotential);
+DECLSPEC extern int refine_grid_for_dh(double *p, double p_ref, int nz,
+                double dp, double *p_i, int ni_max, int *p_indices, int *p_ref_ind_ptr);
 
 #ifdef __cplusplus
 }