From 2db27fb70872a3dec0a2d0556692f9b22d864263 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 9 Aug 2024 16:39:49 -0700
Subject: [PATCH 01/27] Adding gsw_sa_ct_interp

---
 Makefile           |   6 +-
 gsw_sa_ct_interp.c | 533 +++++++++++++++++++++++++++++++++++++++++++++
 gswteos-10.h       |   4 +
 3 files changed, 541 insertions(+), 2 deletions(-)
 create mode 100644 gsw_sa_ct_interp.c

diff --git a/Makefile b/Makefile
index aa0ff06..f3b43f0 100644
--- a/Makefile
+++ b/Makefile
@@ -3,9 +3,11 @@ Program         = gsw_check
 
 PROGRAM_SOURCES = gsw_check_functions.c\
                   gsw_oceanographic_toolbox.c\
-                  gsw_saar.c
+                  gsw_saar.c\
+				  gsw_sa_ct_interp.c
 LIBRARY_SRCS    = gsw_oceanographic_toolbox.c \
-                  gsw_saar.c
+                  gsw_saar.c \
+                  gsw_sa_ct_interp.c
 
 # \
 !ifdef 0 #
diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
new file mode 100644
index 0000000..95114fd
--- /dev/null
+++ b/gsw_sa_ct_interp.c
@@ -0,0 +1,533 @@
+#include "gswteos-10.h"
+#include "gsw_internal_const.h"
+
+/*
+!==========================================================================
+subroutine gsw_util_intersect (x, nx, y, ny, ix, iy)
+!==========================================================================
+!
+! Find unique values common to both x and y arrays.
+!
+!  x      :  x array of values
+!  nx     :  length of x array
+!  y      :  y array of values
+!  ny     :  length of y array
+!  xi     :  indexes of x of values common to y
+!  yi     :  indexes of y of values common to x
+!
+!  ni     :  length of ix and iy
+!
+*/
+int
+gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy)
+{
+        int *sort_ix, *sort_iy, *unique_ix, *unique_iy;
+        int i, uix, uiy, nxs, nys, itx, ity, ni;
+
+        double xx, yy;
+
+        if (nx <= 0 || ny <= 0) {
+            return 0;
+        }
+
+        sort_ix = malloc(nx * sizeof(int));
+        sort_iy = malloc(ny * sizeof(int));
+        unique_ix = malloc(nx * sizeof(int));
+        unique_iy = malloc(ny * sizeof(int));
+
+        gsw_util_sort_real(x, nx, sort_ix);
+        gsw_util_sort_real(y, ny, sort_iy);
+
+        nxs = 0;
+        uix = sort_ix[0];
+        for (i=1; i<nx; ++i) {
+            if (x[uix] != x[sort_ix[i]]) {
+                unique_ix[nxs] = uix;
+                ++nxs;
+                uix = sort_ix[i];
+            } else if (sort_ix[i] < uix) {
+                uix = sort_ix[i];
+            }
+        }
+        unique_ix[nxs] = uix;
+        ++nxs;
+
+        nys = 0;
+        uiy = sort_iy[0];
+        for (i=1; i<ny; ++i) {
+            if (y[uiy] != y[sort_iy[i]]) {
+                unique_iy[nys] = uiy;
+                ++nys;
+                uiy = sort_iy[i];
+            } else if (sort_iy[i] < uiy) {
+                uiy = sort_iy[i];
+            }
+        }
+        unique_iy[nys] = uiy;
+        ++nys;
+
+        itx = ity = ni = 0;
+        while(itx < nxs && ity < nys) {
+            xx = x[unique_ix[itx]];
+            yy = y[unique_iy[ity]];
+            if (xx < yy) {
+                ++itx;
+            } else if (xx > yy) {
+                ++ity;
+            } else {
+                ix[ni] = unique_ix[itx];
+                iy[ni] = unique_iy[ity];
+                ++itx;
+                ++ity;
+                ++ni;
+            }
+        }
+
+        free(sort_ix);
+        free(sort_iy);
+        free(unique_ix);
+        free(unique_iy);
+
+        return ni;
+}
+
+/*
+!==========================================================================
+subroutine gsw_SA_CT_interp (sa,ct,p,p_i)
+!==========================================================================
+!
+! SA and CT interpolation to p_i on a cast
+!
+!  SA   =  Absolute Salinity                                       [ g/kg ]
+!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
+!  p    =  sea pressure                                            [ dbar ]
+!           ( i.e. absolute pressure - 10.1325 dbar )
+!  p_i  =  specific query points at which the interpolated SA_i and CT_i
+!            are required                                          [ dbar ]
+!
+!  SA & CT need to have the same dimensions.
+!  p may have dimensions Mx1 or 1xN or MxN, where SA & CT are MxN.
+!  p_i needs to be either a vector or a matrix and have dimensions M_ix1
+!  or M_ixN.
+!
+!  SA_i  =  interpolated SA values at pressures p_i                [ g/kg ]
+!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
+!
+!--------------------------------------------------------------------------
+*/
+void
+gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p, int mp,
+        int np, double *p_i, int mp_i, int np_i, double *sa_i, double *ct_i)
+{
+        double  factor      = 9.,
+                rec_factor  = 1./factor,
+
+                sin_kpi_on_16[] = {
+                    1.950903220161283e-1,
+                    3.826834323650898e-1,
+                    5.555702330196022e-1,
+                    7.071067811865475e-1,
+                    8.314696123025452e-1,
+                    9.238795325112867e-1,
+                    9.807852804032304e-1
+                },
+                cos_kpi_on_16[] = {
+                    9.807852804032304e-1,
+                    9.238795325112867e-1,
+                    8.314696123025452e-1,
+                    7.071067811865476e-1,
+                    5.555702330196023e-1,
+                    3.826834323650898e-1,
+                    1.950903220161283e-1
+                };
+
+        int i, j, k, i_prof, prof_len, interp_profile_length, 
+            num_interp_profiles, not_monotonic, unique_count, new_len, p_all_len,
+            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
+            i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
+        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
+            *i_out, *i_1, *i_2, *i_3;
+        double d, ct_f, unique_p, ct_sum, sa_sum, min_p_obs, max_p_obs;
+        double *p_mat, *p_i_mat, *sa_obs, *ct_obs, *p_obs, *p_i_tmp,
+               *p_sort, *sa_sort, *ct_sort, *p_all, *p_all_sort,
+               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
+               *independent_variable, *independent_variable_obs_plus_interp,
+               *scaled_sa_obs, *v_tmp, *q_tmp, *v_i, *q_i,
+               *sa_i_obs_plus_interp, *ct_i_obs_plus_interp,
+               *sa_i_limiting_obs_plus_interp, *ct_i_limiting_obs_plus_interp,
+               *ctf_i_tointerp, *sa_i_tooutput, *ct_i_tooutput;
+
+        p_mat = (double *) malloc(m*n*sizeof (double));
+        p_i_mat = (double *) malloc(mp_i*np_i*sizeof (double));
+        sa_obs = (double *) malloc(m*sizeof (double));
+        ct_obs = (double *) malloc(m*sizeof (double));
+        p_obs = (double *) malloc(m*sizeof (double));
+        p_i_tmp = (double *) malloc(m*sizeof (double));
+        p_idx = (int *) malloc(m*sizeof (int));
+        p_sort = (double *) malloc(m*sizeof (double));
+        sa_sort = (double *) malloc(m*sizeof (double));
+        ct_sort = (double *) malloc(m*sizeof (double));
+        p_all = (double *) malloc((m + mp_i)*sizeof (double));
+        p_all_sort = (double *) malloc((m + mp_i)*sizeof (double));
+        p_all_idx = (int *) malloc((m + mp_i)*sizeof (int));
+        i_obs_plus_interp = (int *) malloc((m + mp_i)*sizeof (int));
+        i_surf_and_obs_plus_interp = (int *) malloc((m + mp_i)*sizeof (int));
+        p_obs_plus_interp = (double *) malloc((m + mp_i)*sizeof (double));
+        p_surf_and_obs_plus_interp = (double *) malloc((m + mp_i)*sizeof (double));
+        i_out = (int *) malloc((m + mp_i)*sizeof (int));
+        i_1 = (int *) malloc((m + mp_i)*sizeof (int));
+        i_2 = (int *) malloc((m + mp_i)*sizeof (int));
+        i_3 = (int *) malloc((m + mp_i)*sizeof (int));
+
+        if (m < 4 && np == 1)
+            return; // There must be at least 4 bottles'
+        else if ((n == np && mp == 1) || (n == mp && np == 1)){
+            for (j=0; j<n; j++) {
+                for (i=0; i<m; i++) {
+                    p_mat[m*j + i] = p[j];
+                }
+            }
+        }
+        else if ((m == mp && np == 1) || (m == np && mp == 1)){
+            for (j=0; j<n; j++) {
+                for (i=0; i<m; i++) {
+                    p_mat[m*j + i] = p[i];
+                }
+            }
+        }
+        else if (m == np && n == mp) {
+            for (j=0; j<n; j++) {
+                for (i=0; i<m; i++) {
+                    p_mat[m*j + i] = p[n*i + j];
+                }
+            }
+        }
+        else if (m == mp && n == np) {
+            memcpy(p_mat, p, m*n*sizeof (double));
+        }
+        else
+            return;
+
+        // Get the dimensions of the interpolated arrays
+        if ((mp_i == 1 && np_i > 1) || (mp_i == n && np_i != n)) {
+            interp_profile_length = np_i;
+            num_interp_profiles = mp_i;
+            for (j=0; j<np_i; j++) {
+                for (i=0; i<mp_i; i++) {
+                    p_i_mat[np_i*i + j] = p_i[mp_i*j + i];
+                }
+            }
+        } else {
+            interp_profile_length = mp_i;
+            num_interp_profiles = np_i;
+            memcpy(p_i_mat, p_i, mp_i*np_i*sizeof (double));
+        }
+
+        // Check if interpolating pressure is monotonic
+        for (j=0; j<num_interp_profiles; ++j) {
+            for (i=0; i<interp_profile_length-1; ++i) {
+                d = p_i_mat[j*interp_profile_length + i + 1]
+                    - p_i_mat[j*interp_profile_length + i];
+                if (d < 0.) {
+                    return;
+                }
+            }
+        }
+        
+        for (i_prof=0; i_prof<n; ++i_prof) {
+            // Find NaNs in profile
+            prof_len = 0;
+            for (i=0; i<m; ++i) {
+                d = sa[m*i_prof + i] + ct[m*i_prof + i] + p_mat[m*i_prof + i];
+                if (!isnan(d)) {
+                    sa_obs[prof_len] = sa[m*i_prof + i];
+                    ct_obs[prof_len] = ct[m*i_prof + i];
+
+                    p_obs[prof_len] = 1e-3 * round(1e3 * p_mat[m*i_prof + i]);
+
+                    ct_f = gsw_ct_freezing_poly(sa_obs[prof_len],
+                                                p_obs[prof_len], 
+                                                0.);
+                    if (ct_obs[prof_len] < (ct_f - 0.1)) {
+                        ct_obs[prof_len] = ct_f;
+                    }
+
+                    ++prof_len;
+                }
+            }
+
+            if (prof_len < 2) {
+                for (i=0; i<mp_i; i++) {
+                    sa_i[mp_i*i_prof + i] = NAN;
+                    ct_i[mp_i*i_prof + i] = NAN;
+                }
+                continue;
+            }
+
+            for (i=0; i<interp_profile_length; ++i) {
+                p_i_tmp[i] = 1e-3 * round(1e3 * p_i_mat[interp_profile_length*i_prof + i]);
+            }
+
+            // Check if profile pressure values are monotonic
+            // If they are not, then sort pressure values along with SA and CT
+            not_monotonic = 0;
+            for (i=0; i<prof_len-1; ++i) {
+                d = p_obs[i+1] - p_obs[i];
+                if (d < 0.) {
+                    ++not_monotonic;
+                }
+            }
+
+            if (not_monotonic > 0) {
+                gsw_util_sort_real(p_obs, prof_len, p_idx);
+                for (i=0; i<prof_len; ++i) {
+                    p_sort[i] = p_obs[p_idx[i]];
+                    sa_sort[i] = sa_obs[p_idx[i]];
+                    ct_sort[i] = ct_obs[p_idx[i]];
+                }
+
+                // Once sorted, only save unique pressure values.
+                // SA and CT observations with the same pressure
+                // will be averaged.
+                new_len = 0;
+                unique_count = 1;
+                unique_p = p_sort[0];
+                sa_sum = sa_sort[0];
+                ct_sum = ct_sort[0];
+                for (i=1; i<prof_len; ++i) {
+                    if (unique_p == p_sort[i]) {
+                        sa_sum += sa_sort[i];
+                        ct_sum += ct_sort[i];
+                        ++unique_count;
+                    } else {
+                        p_obs[new_len] = unique_p;
+                        sa_obs[new_len] = sa_sum/(double)unique_count;
+                        ct_obs[new_len] = ct_sum/(double)unique_count;
+                        ++new_len;
+                        unique_p = p_sort[i];
+                        sa_sum = sa_sort[i];
+                        ct_sum = ct_sort[i];
+                        unique_count = 1;
+                    }
+                }
+                p_obs[new_len] = unique_p;
+                sa_obs[new_len] = sa_sum/(double)unique_count;
+                ct_obs[new_len] = ct_sum/(double)unique_count;
+                ++new_len;
+                prof_len = new_len;
+            }
+
+            // Combine pressure values of observed and interpolated SA and CT
+            p_all_len = prof_len + interp_profile_length;
+            memcpy(p_all, p_obs, prof_len*sizeof (double));
+            memcpy(&p_all[prof_len], p_i_tmp, interp_profile_length*sizeof (double));
+            gsw_util_sort_real(p_all, p_all_len, p_all_idx);
+            for (i=0; i<p_all_len; ++i) {
+                p_all_sort[i] = p_all[p_all_idx[i]];
+            }
+            new_len = 0;
+            unique_p = p_all_sort[0];
+            for (i=1; i<p_all_len; ++i) {
+                if (unique_p != p_all_sort[i]) {
+                    p_all[new_len] = unique_p;
+                    ++new_len;
+                    unique_p = p_all_sort[i];
+                }
+            }
+            p_all[new_len] = unique_p;
+            ++new_len;
+            p_all_len = new_len;
+
+            i_min_p_obs = 0;
+            min_p_obs = p_obs[0];
+            max_p_obs = p_obs[0];
+            for (i=1; i<prof_len; ++i) {
+                if (min_p_obs < p_obs[i]) {
+                    i_min_p_obs = i;
+                    min_p_obs = p_obs[i];
+                }
+                max_p_obs = max(max_p_obs, p_obs[i]);
+            }
+
+            i_obs_plus_interp_len = 0;
+            i_surf_and_obs_plus_interp_len = 0;
+            for (i=0; i<p_all_len; ++i) {
+                if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
+                    i_obs_plus_interp[i_obs_plus_interp_len] = i;
+                    p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
+                    i_obs_plus_interp_len++;
+                }
+                if (p_all[i] <= max_p_obs) {
+                    i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
+                    p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
+                    i_surf_and_obs_plus_interp_len++;
+                }
+            }
+            i_out_len = gsw_util_intersect(p_i_tmp, interp_profile_length, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
+            i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
+
+            independent_variable = (double *) malloc(prof_len*sizeof (double));
+            independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+
+            for(i=0; i<prof_len; ++i) {
+                independent_variable[i] = (float)i;
+            }
+            gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
+                                  p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+            
+            scaled_sa_obs = (double *) malloc(prof_len*sizeof (double));
+            v_tmp = (double *) malloc(prof_len*sizeof (double));
+            q_tmp = (double *) malloc(prof_len*sizeof (double));
+            v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+            q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+            sa_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+            ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+            sa_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+            ct_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+            for(i=0; i<prof_len; ++i) {
+                scaled_sa_obs[i] = factor * sa_obs[i];
+            }
+
+            gsw_util_pchip_interp(ct_obs, independent_variable, prof_len,
+                                  ct_i_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(scaled_sa_obs, independent_variable, prof_len,
+                                  q_i, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+
+            for(i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_obs_plus_interp[i] = rec_factor * q_i[i];
+            }
+
+            for(k=0; k<7; ++k) {
+                for(i=0; i<prof_len; ++i) {
+                    v_tmp[i] = scaled_sa_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
+                    q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] + ct_obs[i] * sin_kpi_on_16[k];
+                }
+                gsw_util_pchip_interp(v_tmp, independent_variable, prof_len,
+                                      v_i, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+                gsw_util_pchip_interp(q_tmp, independent_variable, prof_len,
+                                      q_i, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+                for(i=0; i<i_obs_plus_interp_len; ++i) {
+                    ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
+                    sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[i]);
+                }
+            }
+
+            for(i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_obs_plus_interp[i] *= 0.125;
+                ct_i_obs_plus_interp[i] *= 0.125;
+            }
+
+            gsw_linear_interp_sa_ct(sa_obs, ct_obs, independent_variable, prof_len,
+                                    independent_variable_obs_plus_interp, i_obs_plus_interp_len,
+                                    sa_i_limiting_obs_plus_interp, ct_i_limiting_obs_plus_interp);
+
+            ctf_i_tointerp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+            while (1) {
+                for (i=0; i<i_obs_plus_interp_len; ++i) {
+                    d = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_obs_plus_interp[i], 0.);
+                    if (ct_i_limiting_obs_plus_interp[i] < (d - 0.1)) {
+                        ctf_i_tointerp[i] = ct_i_limiting_obs_plus_interp[i];
+                    } else {
+                        ctf_i_tointerp[i] = d;
+                    }
+                }
+
+                for (i_frozen=0; i_frozen<i_obs_plus_interp_len; ++i_frozen) {
+                    if (ct_i_obs_plus_interp[i_frozen] < (ctf_i_tointerp[i_frozen] - 0.1)) {
+                        break;
+                    }
+                }
+
+                if (i_frozen == i_obs_plus_interp_len) {
+                    break;
+                }
+
+                for (i=0; i < i_obs_plus_interp_len && (i_3[i] - i_frozen) <= 0; ++i) {
+                    i_shallower = i;
+                }
+                i_above = i_2[i_shallower];
+                i_above_i = i_3[i_shallower];
+                if ((i_above + 1) > i_3[i_shallower]) {
+                    i_below_i = i_3[i_shallower];
+                } else {
+                    i_below_i = i_3[i_above + 1];
+                }
+
+                for (i=i_above_i; i<=i_below_i; ++i) {
+                    sa_i_obs_plus_interp[i] = sa_i_limiting_obs_plus_interp[i];
+                    ct_i_obs_plus_interp[i] = ct_i_limiting_obs_plus_interp[i];
+                    ctf_i_tointerp[i] = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_all[i], 0.);
+                }
+            }
+
+            sa_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+            ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+
+            if (min_p_obs != 0.) {
+                for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
+                    sa_i_tooutput[i] = NAN;
+                    ct_i_tooutput[i] = NAN;
+                    if (p_i_tmp[i] < min_p_obs) {
+                        sa_i_tooutput[i] = sa_i_obs_plus_interp[i_3[i_min_p_obs]];
+                        ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
+                    }
+                }
+                for (i=0; i<i_obs_plus_interp_len; ++i) {
+                    sa_i_tooutput[i_obs_plus_interp[i]] = sa_i_obs_plus_interp[i];
+                    ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
+                }
+            } else {
+                for (i=0; i<i_obs_plus_interp_len; ++i) {
+                    sa_i_tooutput[i] = sa_i_obs_plus_interp[i];
+                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
+                }
+
+            }
+
+            for(i=0; i<i_out_len; ++i) {
+                sa_i[i_prof*interp_profile_length + i_out[i]] = sa_i_tooutput[i_1[i]];
+                ct_i[i_prof*interp_profile_length + i_out[i]] = ct_i_tooutput[i_1[i]];
+            }
+
+            free(independent_variable);
+            free(independent_variable_obs_plus_interp);
+            free(scaled_sa_obs);
+            free(v_tmp);
+            free(q_tmp);
+            free(v_i);
+            free(q_i);
+            free(sa_i_obs_plus_interp);
+            free(ct_i_obs_plus_interp);
+            free(sa_i_limiting_obs_plus_interp);
+            free(ct_i_limiting_obs_plus_interp);
+            free(ctf_i_tointerp);
+            free(sa_i_tooutput);
+            free(ct_i_tooutput);
+        }
+
+        free(p_mat);
+        free(p_i_mat);
+        free(sa_obs);
+        free(ct_obs);
+        free(p_obs);
+        free(p_i_tmp);
+        free(p_idx);
+        free(p_sort);
+        free(sa_sort);
+        free(ct_sort);
+        free(p_all);
+        free(p_all_sort);
+        free(p_all_idx);
+        free(i_obs_plus_interp);
+        free(i_surf_and_obs_plus_interp);
+        free(p_obs_plus_interp);
+        free(p_surf_and_obs_plus_interp);
+        free(i_out);
+        free(i_1);
+        free(i_2);
+        free(i_3);
+
+        return;
+}
\ No newline at end of file
diff --git a/gswteos-10.h b/gswteos-10.h
index 26463ba..8f122db 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -221,6 +221,9 @@ extern double gsw_rho_t_exact(double sa, double t, double p);
 extern void   gsw_rr68_interp_sa_ct(double *sa, double *ct, double *p, int mp,
                 double *p_i, int mp_i, double *sa_i, double *ct_i);
 extern double gsw_saar(double p, double lon, double lat);
+extern void   gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p,
+                int mp, int np, double *p_i, int mp_i, int np_i, double *sa_i,
+                double *ct_i);
 extern double gsw_sa_freezing_estimate(double p, double saturation_fraction,
                 double *ct, double *t);
 extern double gsw_sa_freezing_from_ct(double ct, double p,
@@ -295,6 +298,7 @@ extern void   gsw_turner_rsubrho(double *sa, double *ct, double *p, int nz,
 extern int    gsw_util_indx(double *x, int n, double z);
 extern double *gsw_util_interp1q_int(int nx, double *x, int *iy, int nxi,
                 double *x_i, double *y_i);
+extern int    gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy);
 extern double *gsw_util_linear_interp(int nx, double *x, int ny, double *y,
                 int nxi, double *x_i, double *y_i);
 extern void   gsw_util_sort_real(double *rarray, int nx, int *iarray);

From 7727e41c5c8fa74d1338a90605dfa44ab00aac4e Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 9 Aug 2024 18:36:27 -0700
Subject: [PATCH 02/27] Fix several mistakes in gsw_sa_ct_interp code

---
 gsw_sa_ct_interp.c | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
index 95114fd..4d45314 100644
--- a/gsw_sa_ct_interp.c
+++ b/gsw_sa_ct_interp.c
@@ -342,7 +342,7 @@ gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p, int mp,
             min_p_obs = p_obs[0];
             max_p_obs = p_obs[0];
             for (i=1; i<prof_len; ++i) {
-                if (min_p_obs < p_obs[i]) {
+                if (p_obs[i] < min_p_obs) {
                     i_min_p_obs = i;
                     min_p_obs = p_obs[i];
                 }
@@ -389,10 +389,10 @@ gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p, int mp,
                 scaled_sa_obs[i] = factor * sa_obs[i];
             }
 
-            gsw_util_pchip_interp(ct_obs, independent_variable, prof_len,
-                                  ct_i_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
-            gsw_util_pchip_interp(scaled_sa_obs, independent_variable, prof_len,
-                                  q_i, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
+                                  independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, scaled_sa_obs, prof_len,
+                                  independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
 
             for(i=0; i<i_obs_plus_interp_len; ++i) {
                 sa_i_obs_plus_interp[i] = rec_factor * q_i[i];
@@ -401,15 +401,15 @@ gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p, int mp,
             for(k=0; k<7; ++k) {
                 for(i=0; i<prof_len; ++i) {
                     v_tmp[i] = scaled_sa_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
-                    q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] + ct_obs[i] * sin_kpi_on_16[k];
+                    q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
                 }
-                gsw_util_pchip_interp(v_tmp, independent_variable, prof_len,
-                                      v_i, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
-                gsw_util_pchip_interp(q_tmp, independent_variable, prof_len,
-                                      q_i, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+                gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
+                                      independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
+                gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
+                                      independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
                 for(i=0; i<i_obs_plus_interp_len; ++i) {
                     ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
-                    sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[i]);
+                    sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
                 }
             }
 

From 26599efedf09a9bd5187c3312588b6455a5d94ef Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 9 Aug 2024 19:10:04 -0700
Subject: [PATCH 03/27] Add DECLSPEC to gsw_sa_ct_interp function

---
 gswteos-10.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/gswteos-10.h b/gswteos-10.h
index e712b30..5f732ac 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -228,7 +228,7 @@ DECLSPEC extern double gsw_rho_t_exact(double sa, double t, double p);
 DECLSPEC extern void   gsw_rr68_interp_sa_ct(double *sa, double *ct, double *p, int mp,
                 double *p_i, int mp_i, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_saar(double p, double lon, double lat);
-extern void   gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p,
+DECLSPEC extern void   gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p,
                 int mp, int np, double *p_i, int mp_i, int np_i, double *sa_i,
                 double *ct_i);
 DECLSPEC extern double gsw_sa_freezing_estimate(double p, double saturation_fraction,
@@ -305,7 +305,7 @@ DECLSPEC extern void   gsw_turner_rsubrho(double *sa, double *ct, double *p, int
 DECLSPEC extern int    gsw_util_indx(double *x, int n, double z);
 DECLSPEC extern double *gsw_util_interp1q_int(int nx, double *x, int *iy, int nxi,
                 double *x_i, double *y_i);
-extern int    gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy);
+DECLSPEC extern int    gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy);
 DECLSPEC extern double *gsw_util_linear_interp(int nx, double *x, int ny, double *y,
                 int nxi, double *x_i, double *y_i);
 DECLSPEC extern void   gsw_util_sort_real(double *rarray, int nx, int *iarray);

From 99d141951f9474cff1714054605a79ba249be663 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 11:09:33 -0700
Subject: [PATCH 04/27] Rework gsw_sa_ct_interp to process one profile at a
 time

---
 gsw_sa_ct_interp.c | 552 ++++++++++++++++++++-------------------------
 gswteos-10.h       |   5 +-
 2 files changed, 252 insertions(+), 305 deletions(-)

diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
index 4d45314..bbbfe16 100644
--- a/gsw_sa_ct_interp.c
+++ b/gsw_sa_ct_interp.c
@@ -104,11 +104,7 @@ subroutine gsw_SA_CT_interp (sa,ct,p,p_i)
 !           ( i.e. absolute pressure - 10.1325 dbar )
 !  p_i  =  specific query points at which the interpolated SA_i and CT_i
 !            are required                                          [ dbar ]
-!
-!  SA & CT need to have the same dimensions.
-!  p may have dimensions Mx1 or 1xN or MxN, where SA & CT are MxN.
-!  p_i needs to be either a vector or a matrix and have dimensions M_ix1
-!  or M_ixN.
+
 !
 !  SA_i  =  interpolated SA values at pressures p_i                [ g/kg ]
 !  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
@@ -116,8 +112,8 @@ subroutine gsw_SA_CT_interp (sa,ct,p,p_i)
 !--------------------------------------------------------------------------
 */
 void
-gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p, int mp,
-        int np, double *p_i, int mp_i, int np_i, double *sa_i, double *ct_i)
+gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
+        double *p_i, int m_i, double *sa_i, double *ct_i)
 {
         double  factor      = 9.,
                 rec_factor  = 1./factor,
@@ -141,14 +137,14 @@ gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p, int mp,
                     1.950903220161283e-1
                 };
 
-        int i, j, k, i_prof, prof_len, interp_profile_length, 
-            num_interp_profiles, not_monotonic, unique_count, new_len, p_all_len,
+        int i, j, k, prof_len, 
+            not_monotonic, unique_count, new_len, p_all_len,
             i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
             i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
         int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
             *i_out, *i_1, *i_2, *i_3;
         double d, ct_f, unique_p, ct_sum, sa_sum, min_p_obs, max_p_obs;
-        double *p_mat, *p_i_mat, *sa_obs, *ct_obs, *p_obs, *p_i_tmp,
+        double *sa_obs, *ct_obs, *p_obs, *p_i_tmp,
                *p_sort, *sa_sort, *ct_sort, *p_all, *p_all_sort,
                *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
                *independent_variable, *independent_variable_obs_plus_interp,
@@ -157,358 +153,310 @@ gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p, int mp,
                *sa_i_limiting_obs_plus_interp, *ct_i_limiting_obs_plus_interp,
                *ctf_i_tointerp, *sa_i_tooutput, *ct_i_tooutput;
 
-        p_mat = (double *) malloc(m*n*sizeof (double));
-        p_i_mat = (double *) malloc(mp_i*np_i*sizeof (double));
         sa_obs = (double *) malloc(m*sizeof (double));
         ct_obs = (double *) malloc(m*sizeof (double));
         p_obs = (double *) malloc(m*sizeof (double));
-        p_i_tmp = (double *) malloc(m*sizeof (double));
         p_idx = (int *) malloc(m*sizeof (int));
         p_sort = (double *) malloc(m*sizeof (double));
         sa_sort = (double *) malloc(m*sizeof (double));
         ct_sort = (double *) malloc(m*sizeof (double));
-        p_all = (double *) malloc((m + mp_i)*sizeof (double));
-        p_all_sort = (double *) malloc((m + mp_i)*sizeof (double));
-        p_all_idx = (int *) malloc((m + mp_i)*sizeof (int));
-        i_obs_plus_interp = (int *) malloc((m + mp_i)*sizeof (int));
-        i_surf_and_obs_plus_interp = (int *) malloc((m + mp_i)*sizeof (int));
-        p_obs_plus_interp = (double *) malloc((m + mp_i)*sizeof (double));
-        p_surf_and_obs_plus_interp = (double *) malloc((m + mp_i)*sizeof (double));
-        i_out = (int *) malloc((m + mp_i)*sizeof (int));
-        i_1 = (int *) malloc((m + mp_i)*sizeof (int));
-        i_2 = (int *) malloc((m + mp_i)*sizeof (int));
-        i_3 = (int *) malloc((m + mp_i)*sizeof (int));
-
-        if (m < 4 && np == 1)
+
+        if (m < 4)
             return; // There must be at least 4 bottles'
-        else if ((n == np && mp == 1) || (n == mp && np == 1)){
-            for (j=0; j<n; j++) {
-                for (i=0; i<m; i++) {
-                    p_mat[m*j + i] = p[j];
-                }
-            }
-        }
-        else if ((m == mp && np == 1) || (m == np && mp == 1)){
-            for (j=0; j<n; j++) {
-                for (i=0; i<m; i++) {
-                    p_mat[m*j + i] = p[i];
-                }
+
+        // Check if interpolating pressure is monotonic
+        for (i=0; i<m_i-1; ++i) {
+            d = p_i[i + 1] - p_i[i];
+            if (d < 0.) {
+                return;
             }
         }
-        else if (m == np && n == mp) {
-            for (j=0; j<n; j++) {
-                for (i=0; i<m; i++) {
-                    p_mat[m*j + i] = p[n*i + j];
+        
+        // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = sa[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                sa_obs[prof_len] = sa[i];
+                ct_obs[prof_len] = ct[i];
+
+                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
+
+                ct_f = gsw_ct_freezing_poly(sa_obs[prof_len],
+                                            p_obs[prof_len], 
+                                            0.);
+                if (ct_obs[prof_len] < (ct_f - 0.1)) {
+                    ct_obs[prof_len] = ct_f;
                 }
+
+                ++prof_len;
             }
         }
-        else if (m == mp && n == np) {
-            memcpy(p_mat, p, m*n*sizeof (double));
-        }
-        else
-            return;
 
-        // Get the dimensions of the interpolated arrays
-        if ((mp_i == 1 && np_i > 1) || (mp_i == n && np_i != n)) {
-            interp_profile_length = np_i;
-            num_interp_profiles = mp_i;
-            for (j=0; j<np_i; j++) {
-                for (i=0; i<mp_i; i++) {
-                    p_i_mat[np_i*i + j] = p_i[mp_i*j + i];
-                }
+        if (prof_len < 2) {
+            for (i=0; i<m_i; i++) {
+                sa_i[i] = NAN;
+                ct_i[i] = NAN;
             }
-        } else {
-            interp_profile_length = mp_i;
-            num_interp_profiles = np_i;
-            memcpy(p_i_mat, p_i, mp_i*np_i*sizeof (double));
+            return;
         }
 
-        // Check if interpolating pressure is monotonic
-        for (j=0; j<num_interp_profiles; ++j) {
-            for (i=0; i<interp_profile_length-1; ++i) {
-                d = p_i_mat[j*interp_profile_length + i + 1]
-                    - p_i_mat[j*interp_profile_length + i];
-                if (d < 0.) {
-                    return;
-                }
-            }
+        p_i_tmp = (double *) malloc(m_i*sizeof (double));
+        for (i=0; i<m_i; ++i) {
+            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
         }
-        
-        for (i_prof=0; i_prof<n; ++i_prof) {
-            // Find NaNs in profile
-            prof_len = 0;
-            for (i=0; i<m; ++i) {
-                d = sa[m*i_prof + i] + ct[m*i_prof + i] + p_mat[m*i_prof + i];
-                if (!isnan(d)) {
-                    sa_obs[prof_len] = sa[m*i_prof + i];
-                    ct_obs[prof_len] = ct[m*i_prof + i];
-
-                    p_obs[prof_len] = 1e-3 * round(1e3 * p_mat[m*i_prof + i]);
-
-                    ct_f = gsw_ct_freezing_poly(sa_obs[prof_len],
-                                                p_obs[prof_len], 
-                                                0.);
-                    if (ct_obs[prof_len] < (ct_f - 0.1)) {
-                        ct_obs[prof_len] = ct_f;
-                    }
-
-                    ++prof_len;
-                }
-            }
 
-            if (prof_len < 2) {
-                for (i=0; i<mp_i; i++) {
-                    sa_i[mp_i*i_prof + i] = NAN;
-                    ct_i[mp_i*i_prof + i] = NAN;
-                }
-                continue;
+        // Check if profile pressure values are monotonic
+        // If they are not, then sort pressure values along with SA and CT
+        not_monotonic = 0;
+        for (i=0; i<prof_len-1; ++i) {
+            d = p_obs[i+1] - p_obs[i];
+            if (d < 0.) {
+                ++not_monotonic;
             }
+        }
 
-            for (i=0; i<interp_profile_length; ++i) {
-                p_i_tmp[i] = 1e-3 * round(1e3 * p_i_mat[interp_profile_length*i_prof + i]);
+        if (not_monotonic > 0) {
+            gsw_util_sort_real(p_obs, prof_len, p_idx);
+            for (i=0; i<prof_len; ++i) {
+                p_sort[i] = p_obs[p_idx[i]];
+                sa_sort[i] = sa_obs[p_idx[i]];
+                ct_sort[i] = ct_obs[p_idx[i]];
             }
 
-            // Check if profile pressure values are monotonic
-            // If they are not, then sort pressure values along with SA and CT
-            not_monotonic = 0;
-            for (i=0; i<prof_len-1; ++i) {
-                d = p_obs[i+1] - p_obs[i];
-                if (d < 0.) {
-                    ++not_monotonic;
+            // Once sorted, only save unique pressure values.
+            // SA and CT observations with the same pressure
+            // will be averaged.
+            new_len = 0;
+            unique_count = 1;
+            unique_p = p_sort[0];
+            sa_sum = sa_sort[0];
+            ct_sum = ct_sort[0];
+            for (i=1; i<prof_len; ++i) {
+                if (unique_p == p_sort[i]) {
+                    sa_sum += sa_sort[i];
+                    ct_sum += ct_sort[i];
+                    ++unique_count;
+                } else {
+                    p_obs[new_len] = unique_p;
+                    sa_obs[new_len] = sa_sum/(double)unique_count;
+                    ct_obs[new_len] = ct_sum/(double)unique_count;
+                    ++new_len;
+                    unique_p = p_sort[i];
+                    sa_sum = sa_sort[i];
+                    ct_sum = ct_sort[i];
+                    unique_count = 1;
                 }
             }
+            p_obs[new_len] = unique_p;
+            sa_obs[new_len] = sa_sum/(double)unique_count;
+            ct_obs[new_len] = ct_sum/(double)unique_count;
+            ++new_len;
+            prof_len = new_len;
+        }
 
-            if (not_monotonic > 0) {
-                gsw_util_sort_real(p_obs, prof_len, p_idx);
-                for (i=0; i<prof_len; ++i) {
-                    p_sort[i] = p_obs[p_idx[i]];
-                    sa_sort[i] = sa_obs[p_idx[i]];
-                    ct_sort[i] = ct_obs[p_idx[i]];
-                }
+        // Combine pressure values of observed and interpolated SA and CT
+        p_all_len = prof_len + m_i;
 
-                // Once sorted, only save unique pressure values.
-                // SA and CT observations with the same pressure
-                // will be averaged.
-                new_len = 0;
-                unique_count = 1;
-                unique_p = p_sort[0];
-                sa_sum = sa_sort[0];
-                ct_sum = ct_sort[0];
-                for (i=1; i<prof_len; ++i) {
-                    if (unique_p == p_sort[i]) {
-                        sa_sum += sa_sort[i];
-                        ct_sum += ct_sort[i];
-                        ++unique_count;
-                    } else {
-                        p_obs[new_len] = unique_p;
-                        sa_obs[new_len] = sa_sum/(double)unique_count;
-                        ct_obs[new_len] = ct_sum/(double)unique_count;
-                        ++new_len;
-                        unique_p = p_sort[i];
-                        sa_sum = sa_sort[i];
-                        ct_sum = ct_sort[i];
-                        unique_count = 1;
-                    }
-                }
-                p_obs[new_len] = unique_p;
-                sa_obs[new_len] = sa_sum/(double)unique_count;
-                ct_obs[new_len] = ct_sum/(double)unique_count;
-                ++new_len;
-                prof_len = new_len;
-            }
+        p_all = (double *) malloc(p_all_len*sizeof (double));
+        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
+        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
 
-            // Combine pressure values of observed and interpolated SA and CT
-            p_all_len = prof_len + interp_profile_length;
-            memcpy(p_all, p_obs, prof_len*sizeof (double));
-            memcpy(&p_all[prof_len], p_i_tmp, interp_profile_length*sizeof (double));
-            gsw_util_sort_real(p_all, p_all_len, p_all_idx);
-            for (i=0; i<p_all_len; ++i) {
-                p_all_sort[i] = p_all[p_all_idx[i]];
+        memcpy(p_all, p_obs, prof_len*sizeof (double));
+        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
+        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
+        for (i=0; i<p_all_len; ++i) {
+            p_all_sort[i] = p_all[p_all_idx[i]];
+        }
+        new_len = 0;
+        unique_p = p_all_sort[0];
+        for (i=1; i<p_all_len; ++i) {
+            if (unique_p != p_all_sort[i]) {
+                p_all[new_len] = unique_p;
+                ++new_len;
+                unique_p = p_all_sort[i];
             }
-            new_len = 0;
-            unique_p = p_all_sort[0];
-            for (i=1; i<p_all_len; ++i) {
-                if (unique_p != p_all_sort[i]) {
-                    p_all[new_len] = unique_p;
-                    ++new_len;
-                    unique_p = p_all_sort[i];
-                }
+        }
+        p_all[new_len] = unique_p;
+        ++new_len;
+        p_all_len = new_len;
+
+        i_min_p_obs = 0;
+        min_p_obs = p_obs[0];
+        max_p_obs = p_obs[0];
+        for (i=1; i<prof_len; ++i) {
+            if (p_obs[i] < min_p_obs) {
+                i_min_p_obs = i;
+                min_p_obs = p_obs[i];
             }
-            p_all[new_len] = unique_p;
-            ++new_len;
-            p_all_len = new_len;
+            max_p_obs = max(max_p_obs, p_obs[i]);
+        }
 
-            i_min_p_obs = 0;
-            min_p_obs = p_obs[0];
-            max_p_obs = p_obs[0];
-            for (i=1; i<prof_len; ++i) {
-                if (p_obs[i] < min_p_obs) {
-                    i_min_p_obs = i;
-                    min_p_obs = p_obs[i];
-                }
-                max_p_obs = max(max_p_obs, p_obs[i]);
+        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        i_out = (int *) malloc(p_all_len*sizeof (int));
+        i_1 = (int *) malloc(p_all_len*sizeof (int));
+        i_2 = (int *) malloc(p_all_len*sizeof (int));
+        i_3 = (int *) malloc(p_all_len*sizeof (int));
+
+        i_obs_plus_interp_len = 0;
+        i_surf_and_obs_plus_interp_len = 0;
+        for (i=0; i<p_all_len; ++i) {
+            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
+                i_obs_plus_interp[i_obs_plus_interp_len] = i;
+                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
+                i_obs_plus_interp_len++;
             }
-
-            i_obs_plus_interp_len = 0;
-            i_surf_and_obs_plus_interp_len = 0;
-            for (i=0; i<p_all_len; ++i) {
-                if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
-                    i_obs_plus_interp[i_obs_plus_interp_len] = i;
-                    p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
-                    i_obs_plus_interp_len++;
-                }
-                if (p_all[i] <= max_p_obs) {
-                    i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
-                    p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
-                    i_surf_and_obs_plus_interp_len++;
-                }
+            if (p_all[i] <= max_p_obs) {
+                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
+                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
+                i_surf_and_obs_plus_interp_len++;
             }
-            i_out_len = gsw_util_intersect(p_i_tmp, interp_profile_length, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
-            i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
+        }
+        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
+        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
 
-            independent_variable = (double *) malloc(prof_len*sizeof (double));
-            independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        independent_variable = (double *) malloc(prof_len*sizeof (double));
+        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
 
-            for(i=0; i<prof_len; ++i) {
-                independent_variable[i] = (float)i;
-            }
-            gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
-                                  p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
-            
-            scaled_sa_obs = (double *) malloc(prof_len*sizeof (double));
-            v_tmp = (double *) malloc(prof_len*sizeof (double));
-            q_tmp = (double *) malloc(prof_len*sizeof (double));
-            v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-            q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-            sa_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-            ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-            sa_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-            ct_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-
-            for(i=0; i<prof_len; ++i) {
-                scaled_sa_obs[i] = factor * sa_obs[i];
-            }
+        for(i=0; i<prof_len; ++i) {
+            independent_variable[i] = (float)i;
+        }
+        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
+                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+        
+        scaled_sa_obs = (double *) malloc(prof_len*sizeof (double));
+        v_tmp = (double *) malloc(prof_len*sizeof (double));
+        q_tmp = (double *) malloc(prof_len*sizeof (double));
+        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        sa_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        sa_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            scaled_sa_obs[i] = factor * sa_obs[i];
+        }
 
-            gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
-                                  independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
-            gsw_util_pchip_interp(independent_variable, scaled_sa_obs, prof_len,
-                                  independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
+                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
+        gsw_util_pchip_interp(independent_variable, scaled_sa_obs, prof_len,
+                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
 
-            for(i=0; i<i_obs_plus_interp_len; ++i) {
-                sa_i_obs_plus_interp[i] = rec_factor * q_i[i];
-            }
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            sa_i_obs_plus_interp[i] = rec_factor * q_i[i];
+        }
 
-            for(k=0; k<7; ++k) {
-                for(i=0; i<prof_len; ++i) {
-                    v_tmp[i] = scaled_sa_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
-                    q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
-                }
-                gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
-                                      independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
-                gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
-                                      independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
-                for(i=0; i<i_obs_plus_interp_len; ++i) {
-                    ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
-                    sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
-                }
+        for(k=0; k<7; ++k) {
+            for(i=0; i<prof_len; ++i) {
+                v_tmp[i] = scaled_sa_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
+                q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
             }
-
+            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
             for(i=0; i<i_obs_plus_interp_len; ++i) {
-                sa_i_obs_plus_interp[i] *= 0.125;
-                ct_i_obs_plus_interp[i] *= 0.125;
+                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
+                sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
             }
+        }
 
-            gsw_linear_interp_sa_ct(sa_obs, ct_obs, independent_variable, prof_len,
-                                    independent_variable_obs_plus_interp, i_obs_plus_interp_len,
-                                    sa_i_limiting_obs_plus_interp, ct_i_limiting_obs_plus_interp);
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            sa_i_obs_plus_interp[i] *= 0.125;
+            ct_i_obs_plus_interp[i] *= 0.125;
+        }
 
-            ctf_i_tointerp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        gsw_linear_interp_sa_ct(sa_obs, ct_obs, independent_variable, prof_len,
+                                independent_variable_obs_plus_interp, i_obs_plus_interp_len,
+                                sa_i_limiting_obs_plus_interp, ct_i_limiting_obs_plus_interp);
 
-            while (1) {
-                for (i=0; i<i_obs_plus_interp_len; ++i) {
-                    d = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_obs_plus_interp[i], 0.);
-                    if (ct_i_limiting_obs_plus_interp[i] < (d - 0.1)) {
-                        ctf_i_tointerp[i] = ct_i_limiting_obs_plus_interp[i];
-                    } else {
-                        ctf_i_tointerp[i] = d;
-                    }
-                }
+        ctf_i_tointerp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
 
-                for (i_frozen=0; i_frozen<i_obs_plus_interp_len; ++i_frozen) {
-                    if (ct_i_obs_plus_interp[i_frozen] < (ctf_i_tointerp[i_frozen] - 0.1)) {
-                        break;
-                    }
+        while (1) {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                d = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_obs_plus_interp[i], 0.);
+                if (ct_i_limiting_obs_plus_interp[i] < (d - 0.1)) {
+                    ctf_i_tointerp[i] = ct_i_limiting_obs_plus_interp[i];
+                } else {
+                    ctf_i_tointerp[i] = d;
                 }
+            }
 
-                if (i_frozen == i_obs_plus_interp_len) {
+            for (i_frozen=0; i_frozen<i_obs_plus_interp_len; ++i_frozen) {
+                if (ct_i_obs_plus_interp[i_frozen] < (ctf_i_tointerp[i_frozen] - 0.1)) {
                     break;
                 }
+            }
 
-                for (i=0; i < i_obs_plus_interp_len && (i_3[i] - i_frozen) <= 0; ++i) {
-                    i_shallower = i;
-                }
-                i_above = i_2[i_shallower];
-                i_above_i = i_3[i_shallower];
-                if ((i_above + 1) > i_3[i_shallower]) {
-                    i_below_i = i_3[i_shallower];
-                } else {
-                    i_below_i = i_3[i_above + 1];
-                }
-
-                for (i=i_above_i; i<=i_below_i; ++i) {
-                    sa_i_obs_plus_interp[i] = sa_i_limiting_obs_plus_interp[i];
-                    ct_i_obs_plus_interp[i] = ct_i_limiting_obs_plus_interp[i];
-                    ctf_i_tointerp[i] = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_all[i], 0.);
-                }
+            if (i_frozen == i_obs_plus_interp_len) {
+                break;
             }
 
-            sa_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
-            ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
-
-            if (min_p_obs != 0.) {
-                for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
-                    sa_i_tooutput[i] = NAN;
-                    ct_i_tooutput[i] = NAN;
-                    if (p_i_tmp[i] < min_p_obs) {
-                        sa_i_tooutput[i] = sa_i_obs_plus_interp[i_3[i_min_p_obs]];
-                        ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
-                    }
-                }
-                for (i=0; i<i_obs_plus_interp_len; ++i) {
-                    sa_i_tooutput[i_obs_plus_interp[i]] = sa_i_obs_plus_interp[i];
-                    ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
-                }
+            for (i=0; i < i_obs_plus_interp_len && (i_3[i] - i_frozen) <= 0; ++i) {
+                i_shallower = i;
+            }
+            i_above = i_2[i_shallower];
+            i_above_i = i_3[i_shallower];
+            if ((i_above + 1) > i_3[i_shallower]) {
+                i_below_i = i_3[i_shallower];
             } else {
-                for (i=0; i<i_obs_plus_interp_len; ++i) {
-                    sa_i_tooutput[i] = sa_i_obs_plus_interp[i];
-                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
-                }
+                i_below_i = i_3[i_above + 1];
+            }
 
+            for (i=i_above_i; i<=i_below_i; ++i) {
+                sa_i_obs_plus_interp[i] = sa_i_limiting_obs_plus_interp[i];
+                ct_i_obs_plus_interp[i] = ct_i_limiting_obs_plus_interp[i];
+                ctf_i_tointerp[i] = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_all[i], 0.);
             }
+        }
 
-            for(i=0; i<i_out_len; ++i) {
-                sa_i[i_prof*interp_profile_length + i_out[i]] = sa_i_tooutput[i_1[i]];
-                ct_i[i_prof*interp_profile_length + i_out[i]] = ct_i_tooutput[i_1[i]];
+        sa_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+
+        if (min_p_obs != 0.) {
+            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i] = NAN;
+                ct_i_tooutput[i] = NAN;
+                if (p_i_tmp[i] < min_p_obs) {
+                    sa_i_tooutput[i] = sa_i_obs_plus_interp[i_3[i_min_p_obs]];
+                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
+                }
             }
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i_obs_plus_interp[i]] = sa_i_obs_plus_interp[i];
+                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
+            }
+        } else {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i] = sa_i_obs_plus_interp[i];
+                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
+            }
+
+        }
 
-            free(independent_variable);
-            free(independent_variable_obs_plus_interp);
-            free(scaled_sa_obs);
-            free(v_tmp);
-            free(q_tmp);
-            free(v_i);
-            free(q_i);
-            free(sa_i_obs_plus_interp);
-            free(ct_i_obs_plus_interp);
-            free(sa_i_limiting_obs_plus_interp);
-            free(ct_i_limiting_obs_plus_interp);
-            free(ctf_i_tointerp);
-            free(sa_i_tooutput);
-            free(ct_i_tooutput);
+        for(i=0; i<i_out_len; ++i) {
+            sa_i[i_out[i]] = sa_i_tooutput[i_1[i]];
+            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
         }
 
-        free(p_mat);
-        free(p_i_mat);
+        free(independent_variable);
+        free(independent_variable_obs_plus_interp);
+        free(scaled_sa_obs);
+        free(v_tmp);
+        free(q_tmp);
+        free(v_i);
+        free(q_i);
+        free(sa_i_obs_plus_interp);
+        free(ct_i_obs_plus_interp);
+        free(sa_i_limiting_obs_plus_interp);
+        free(ct_i_limiting_obs_plus_interp);
+        free(ctf_i_tointerp);
+        free(sa_i_tooutput);
+        free(ct_i_tooutput);
         free(sa_obs);
         free(ct_obs);
         free(p_obs);
diff --git a/gswteos-10.h b/gswteos-10.h
index 5f732ac..9ae0db4 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -228,9 +228,8 @@ DECLSPEC extern double gsw_rho_t_exact(double sa, double t, double p);
 DECLSPEC extern void   gsw_rr68_interp_sa_ct(double *sa, double *ct, double *p, int mp,
                 double *p_i, int mp_i, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_saar(double p, double lon, double lat);
-DECLSPEC extern void   gsw_sa_ct_interp(double *sa, double *ct, int m, int n, double *p,
-                int mp, int np, double *p_i, int mp_i, int np_i, double *sa_i,
-                double *ct_i);
+DECLSPEC extern void   gsw_sa_ct_interp(double *sa, double *ct, double *p,
+                int m, double *p_i, int m_i, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_sa_freezing_estimate(double p, double saturation_fraction,
                 double *ct, double *t);
 DECLSPEC extern double gsw_sa_freezing_from_ct(double ct, double p,

From 83411d06bee462b1fdff4584ba8e969e1c8d4ad7 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 11:43:08 -0700
Subject: [PATCH 05/27] Initialize interpolated SA and CT arrays with NaNs

---
 gsw_sa_ct_interp.c | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
index bbbfe16..dbfc7eb 100644
--- a/gsw_sa_ct_interp.c
+++ b/gsw_sa_ct_interp.c
@@ -171,6 +171,11 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
                 return;
             }
         }
+
+        for (i=0; i<m_i; i++) {
+            sa_i[i] = NAN;
+            ct_i[i] = NAN;
+        }
         
         // Find NaNs in profile
         prof_len = 0;
@@ -194,10 +199,6 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         }
 
         if (prof_len < 2) {
-            for (i=0; i<m_i; i++) {
-                sa_i[i] = NAN;
-                ct_i[i] = NAN;
-            }
             return;
         }
 

From 5ae5ea2083e2730e16cf19f773d06362a93cdf1f Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 13:29:46 -0700
Subject: [PATCH 06/27] Move gsw_util_intersect to gsw_oceanographic_toolbox.c

---
 gsw_oceanographic_toolbox.c | 90 +++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index 371bce2..e783077 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -11489,6 +11489,96 @@ gsw_z_from_p(double p, double lat, double geo_strf_dyn_height,
         return (-2.0*c/(b + sqrt(b*b - 4.0*a*c)));
 }
 
+/*
+!==========================================================================
+subroutine gsw_util_intersect (x, nx, y, ny, ix, iy)
+!==========================================================================
+!
+! Find unique values common to both x and y arrays.
+!
+!  x      :  x array of values
+!  nx     :  length of x array
+!  y      :  y array of values
+!  ny     :  length of y array
+!  xi     :  indexes of x of values common to y
+!  yi     :  indexes of y of values common to x
+!
+!  ni     :  length of ix and iy
+!
+*/
+int
+gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy)
+{
+        int *sort_ix, *sort_iy, *unique_ix, *unique_iy;
+        int i, uix, uiy, nxs, nys, itx, ity, ni;
+
+        double xx, yy;
+
+        if (nx <= 0 || ny <= 0) {
+            return 0;
+        }
+
+        sort_ix = malloc(nx * sizeof(int));
+        sort_iy = malloc(ny * sizeof(int));
+        unique_ix = malloc(nx * sizeof(int));
+        unique_iy = malloc(ny * sizeof(int));
+
+        gsw_util_sort_real(x, nx, sort_ix);
+        gsw_util_sort_real(y, ny, sort_iy);
+
+        nxs = 0;
+        uix = sort_ix[0];
+        for (i=1; i<nx; ++i) {
+            if (x[uix] != x[sort_ix[i]]) {
+                unique_ix[nxs] = uix;
+                ++nxs;
+                uix = sort_ix[i];
+            } else if (sort_ix[i] < uix) {
+                uix = sort_ix[i];
+            }
+        }
+        unique_ix[nxs] = uix;
+        ++nxs;
+
+        nys = 0;
+        uiy = sort_iy[0];
+        for (i=1; i<ny; ++i) {
+            if (y[uiy] != y[sort_iy[i]]) {
+                unique_iy[nys] = uiy;
+                ++nys;
+                uiy = sort_iy[i];
+            } else if (sort_iy[i] < uiy) {
+                uiy = sort_iy[i];
+            }
+        }
+        unique_iy[nys] = uiy;
+        ++nys;
+
+        itx = ity = ni = 0;
+        while(itx < nxs && ity < nys) {
+            xx = x[unique_ix[itx]];
+            yy = y[unique_iy[ity]];
+            if (xx < yy) {
+                ++itx;
+            } else if (xx > yy) {
+                ++ity;
+            } else {
+                ix[ni] = unique_ix[itx];
+                iy[ni] = unique_iy[ity];
+                ++itx;
+                ++ity;
+                ++ni;
+            }
+        }
+
+        free(sort_ix);
+        free(sort_iy);
+        free(unique_ix);
+        free(unique_iy);
+
+        return ni;
+}
+
 /*
 **  The End
 **!==========================================================================

From b16962367d16db2bca7c7c9e7ee5070d06dddb65 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 13:40:41 -0700
Subject: [PATCH 07/27] Remove gsw_util_intersect from gsw_sa_ct_interp.c

---
 gsw_sa_ct_interp.c | 90 ----------------------------------------------
 1 file changed, 90 deletions(-)

diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
index dbfc7eb..7972ed3 100644
--- a/gsw_sa_ct_interp.c
+++ b/gsw_sa_ct_interp.c
@@ -1,96 +1,6 @@
 #include "gswteos-10.h"
 #include "gsw_internal_const.h"
 
-/*
-!==========================================================================
-subroutine gsw_util_intersect (x, nx, y, ny, ix, iy)
-!==========================================================================
-!
-! Find unique values common to both x and y arrays.
-!
-!  x      :  x array of values
-!  nx     :  length of x array
-!  y      :  y array of values
-!  ny     :  length of y array
-!  xi     :  indexes of x of values common to y
-!  yi     :  indexes of y of values common to x
-!
-!  ni     :  length of ix and iy
-!
-*/
-int
-gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy)
-{
-        int *sort_ix, *sort_iy, *unique_ix, *unique_iy;
-        int i, uix, uiy, nxs, nys, itx, ity, ni;
-
-        double xx, yy;
-
-        if (nx <= 0 || ny <= 0) {
-            return 0;
-        }
-
-        sort_ix = malloc(nx * sizeof(int));
-        sort_iy = malloc(ny * sizeof(int));
-        unique_ix = malloc(nx * sizeof(int));
-        unique_iy = malloc(ny * sizeof(int));
-
-        gsw_util_sort_real(x, nx, sort_ix);
-        gsw_util_sort_real(y, ny, sort_iy);
-
-        nxs = 0;
-        uix = sort_ix[0];
-        for (i=1; i<nx; ++i) {
-            if (x[uix] != x[sort_ix[i]]) {
-                unique_ix[nxs] = uix;
-                ++nxs;
-                uix = sort_ix[i];
-            } else if (sort_ix[i] < uix) {
-                uix = sort_ix[i];
-            }
-        }
-        unique_ix[nxs] = uix;
-        ++nxs;
-
-        nys = 0;
-        uiy = sort_iy[0];
-        for (i=1; i<ny; ++i) {
-            if (y[uiy] != y[sort_iy[i]]) {
-                unique_iy[nys] = uiy;
-                ++nys;
-                uiy = sort_iy[i];
-            } else if (sort_iy[i] < uiy) {
-                uiy = sort_iy[i];
-            }
-        }
-        unique_iy[nys] = uiy;
-        ++nys;
-
-        itx = ity = ni = 0;
-        while(itx < nxs && ity < nys) {
-            xx = x[unique_ix[itx]];
-            yy = y[unique_iy[ity]];
-            if (xx < yy) {
-                ++itx;
-            } else if (xx > yy) {
-                ++ity;
-            } else {
-                ix[ni] = unique_ix[itx];
-                iy[ni] = unique_iy[ity];
-                ++itx;
-                ++ity;
-                ++ni;
-            }
-        }
-
-        free(sort_ix);
-        free(sort_iy);
-        free(unique_ix);
-        free(unique_iy);
-
-        return ni;
-}
-
 /*
 !==========================================================================
 subroutine gsw_SA_CT_interp (sa,ct,p,p_i)

From fe0cc28d5d9df44b49294208291aaad21ba2a4f0 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 14:03:30 -0700
Subject: [PATCH 08/27] Add gsw_tracer_ct_interp

---
 gsw_tracer_ct_interp.c | 334 +++++++++++++++++++++++++++++++++++++++++
 gswteos-10.h           |   2 +
 2 files changed, 336 insertions(+)
 create mode 100644 gsw_tracer_ct_interp.c

diff --git a/gsw_tracer_ct_interp.c b/gsw_tracer_ct_interp.c
new file mode 100644
index 0000000..3cf4350
--- /dev/null
+++ b/gsw_tracer_ct_interp.c
@@ -0,0 +1,334 @@
+#include "gswteos-10.h"
+#include "gsw_internal_const.h"
+
+/*
+!==========================================================================
+subroutine gsw_tracer_CT_interp (sa,ct,p,p_i)
+!==========================================================================
+!
+! SA and CT interpolation to p_i on a cast
+!
+!  tracer   =  tracer                                                 [ ? ]
+!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
+!  p    =  sea pressure                                            [ dbar ]
+!           ( i.e. absolute pressure - 10.1325 dbar )
+!  p_i  =  specific query points at which the interpolated tracer_i and CT_i
+!            are required                                          [ dbar ]
+
+!
+!  tracer_i  =  interpolated SA values at pressures p_i                [ g/kg ]
+!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
+!
+!--------------------------------------------------------------------------
+*/
+void
+gsw_tracer_ct_interp(double *sa, double *ct, double *p, int m,
+        double *p_i, int m_i, double factor, double *tracer_i, double *ct_i)
+{
+        double  rec_factor  = 1./factor,
+
+                sin_kpi_on_16[] = {
+                    1.950903220161283e-1,
+                    3.826834323650898e-1,
+                    5.555702330196022e-1,
+                    7.071067811865475e-1,
+                    8.314696123025452e-1,
+                    9.238795325112867e-1,
+                    9.807852804032304e-1
+                },
+                cos_kpi_on_16[] = {
+                    9.807852804032304e-1,
+                    9.238795325112867e-1,
+                    8.314696123025452e-1,
+                    7.071067811865476e-1,
+                    5.555702330196023e-1,
+                    3.826834323650898e-1,
+                    1.950903220161283e-1
+                };
+
+        int i, j, k, prof_len, 
+            not_monotonic, unique_count, new_len, p_all_len,
+            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
+            i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
+        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
+            *i_out, *i_1, *i_2, *i_3;
+        double d, ct_f, unique_p, ct_sum, tracer_sum, min_p_obs, max_p_obs;
+        double *tracer_obs, *ct_obs, *p_obs, *p_i_tmp,
+               *p_sort, *tracer_sort, *ct_sort, *p_all, *p_all_sort,
+               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
+               *independent_variable, *independent_variable_obs_plus_interp,
+               *scaled_tracer_obs, *v_tmp, *q_tmp, *v_i, *q_i,
+               *tracer_i_obs_plus_interp, *ct_i_obs_plus_interp,
+               *tracer_i_tooutput, *ct_i_tooutput;
+
+        tracer_obs = (double *) malloc(m*sizeof (double));
+        ct_obs = (double *) malloc(m*sizeof (double));
+        p_obs = (double *) malloc(m*sizeof (double));
+        p_idx = (int *) malloc(m*sizeof (int));
+        p_sort = (double *) malloc(m*sizeof (double));
+        tracer_sort = (double *) malloc(m*sizeof (double));
+        ct_sort = (double *) malloc(m*sizeof (double));
+
+        if (m < 4)
+            return; // There must be at least 4 bottles'
+
+        // Check if interpolating pressure is monotonic
+        for (i=0; i<m_i-1; ++i) {
+            d = p_i[i + 1] - p_i[i];
+            if (d < 0.) {
+                return;
+            }
+        }
+
+        for (i=0; i<m_i; i++) {
+            tracer_i[i] = NAN;
+            ct_i[i] = NAN;
+        }
+        
+        // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = sa[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                tracer_obs[prof_len] = sa[i];
+                ct_obs[prof_len] = ct[i];
+
+                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
+
+                ++prof_len;
+            }
+        }
+
+        if (prof_len < 2) {
+            return;
+        }
+
+        p_i_tmp = (double *) malloc(m_i*sizeof (double));
+        for (i=0; i<m_i; ++i) {
+            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
+        }
+
+        // Check if profile pressure values are monotonic
+        // If they are not, then sort pressure values along with SA and CT
+        not_monotonic = 0;
+        for (i=0; i<prof_len-1; ++i) {
+            d = p_obs[i+1] - p_obs[i];
+            if (d < 0.) {
+                ++not_monotonic;
+            }
+        }
+
+        if (not_monotonic > 0) {
+            gsw_util_sort_real(p_obs, prof_len, p_idx);
+            for (i=0; i<prof_len; ++i) {
+                p_sort[i] = p_obs[p_idx[i]];
+                tracer_sort[i] = tracer_obs[p_idx[i]];
+                ct_sort[i] = ct_obs[p_idx[i]];
+            }
+
+            // Once sorted, only save unique pressure values.
+            // SA and CT observations with the same pressure
+            // will be averaged.
+            new_len = 0;
+            unique_count = 1;
+            unique_p = p_sort[0];
+            tracer_sum = tracer_sort[0];
+            ct_sum = ct_sort[0];
+            for (i=1; i<prof_len; ++i) {
+                if (unique_p == p_sort[i]) {
+                    tracer_sum += tracer_sort[i];
+                    ct_sum += ct_sort[i];
+                    ++unique_count;
+                } else {
+                    p_obs[new_len] = unique_p;
+                    tracer_obs[new_len] = tracer_sum/(double)unique_count;
+                    ct_obs[new_len] = ct_sum/(double)unique_count;
+                    ++new_len;
+                    unique_p = p_sort[i];
+                    tracer_sum = tracer_sort[i];
+                    ct_sum = ct_sort[i];
+                    unique_count = 1;
+                }
+            }
+            p_obs[new_len] = unique_p;
+            tracer_obs[new_len] = tracer_sum/(double)unique_count;
+            ct_obs[new_len] = ct_sum/(double)unique_count;
+            ++new_len;
+            prof_len = new_len;
+        }
+
+        // Combine pressure values of observed and interpolated SA and CT
+        p_all_len = prof_len + m_i;
+
+        p_all = (double *) malloc(p_all_len*sizeof (double));
+        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
+        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
+
+        memcpy(p_all, p_obs, prof_len*sizeof (double));
+        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
+        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
+        for (i=0; i<p_all_len; ++i) {
+            p_all_sort[i] = p_all[p_all_idx[i]];
+        }
+        new_len = 0;
+        unique_p = p_all_sort[0];
+        for (i=1; i<p_all_len; ++i) {
+            if (unique_p != p_all_sort[i]) {
+                p_all[new_len] = unique_p;
+                ++new_len;
+                unique_p = p_all_sort[i];
+            }
+        }
+        p_all[new_len] = unique_p;
+        ++new_len;
+        p_all_len = new_len;
+
+        i_min_p_obs = 0;
+        min_p_obs = p_obs[0];
+        max_p_obs = p_obs[0];
+        for (i=1; i<prof_len; ++i) {
+            if (p_obs[i] < min_p_obs) {
+                i_min_p_obs = i;
+                min_p_obs = p_obs[i];
+            }
+            max_p_obs = max(max_p_obs, p_obs[i]);
+        }
+
+        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        i_out = (int *) malloc(p_all_len*sizeof (int));
+        i_1 = (int *) malloc(p_all_len*sizeof (int));
+        i_2 = (int *) malloc(p_all_len*sizeof (int));
+        i_3 = (int *) malloc(p_all_len*sizeof (int));
+
+        i_obs_plus_interp_len = 0;
+        i_surf_and_obs_plus_interp_len = 0;
+        for (i=0; i<p_all_len; ++i) {
+            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
+                i_obs_plus_interp[i_obs_plus_interp_len] = i;
+                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
+                i_obs_plus_interp_len++;
+            }
+            if (p_all[i] <= max_p_obs) {
+                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
+                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
+                i_surf_and_obs_plus_interp_len++;
+            }
+        }
+        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
+        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
+
+        independent_variable = (double *) malloc(prof_len*sizeof (double));
+        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            independent_variable[i] = (float)i;
+        }
+        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
+                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+        
+        scaled_tracer_obs = (double *) malloc(prof_len*sizeof (double));
+        v_tmp = (double *) malloc(prof_len*sizeof (double));
+        q_tmp = (double *) malloc(prof_len*sizeof (double));
+        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        tracer_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            scaled_tracer_obs[i] = factor * tracer_obs[i];
+        }
+
+        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
+                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
+        gsw_util_pchip_interp(independent_variable, scaled_tracer_obs, prof_len,
+                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            tracer_i_obs_plus_interp[i] = rec_factor * q_i[i];
+        }
+
+        for(k=0; k<7; ++k) {
+            for(i=0; i<prof_len; ++i) {
+                v_tmp[i] = scaled_tracer_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
+                q_tmp[i] = scaled_tracer_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
+            }
+            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+            for(i=0; i<i_obs_plus_interp_len; ++i) {
+                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
+                tracer_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
+            }
+        }
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            tracer_i_obs_plus_interp[i] *= 0.125;
+            ct_i_obs_plus_interp[i] *= 0.125;
+        }
+
+        tracer_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+
+        if (min_p_obs != 0.) {
+            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i] = NAN;
+                ct_i_tooutput[i] = NAN;
+                if (p_i_tmp[i] < min_p_obs) {
+                    tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i_3[i_min_p_obs]];
+                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
+                }
+            }
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i_obs_plus_interp[i]] = tracer_i_obs_plus_interp[i];
+                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
+            }
+        } else {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i];
+                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
+            }
+
+        }
+
+        for(i=0; i<i_out_len; ++i) {
+            tracer_i[i_out[i]] = tracer_i_tooutput[i_1[i]];
+            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
+        }
+
+        free(independent_variable);
+        free(independent_variable_obs_plus_interp);
+        free(scaled_tracer_obs);
+        free(v_tmp);
+        free(q_tmp);
+        free(v_i);
+        free(q_i);
+        free(tracer_i_obs_plus_interp);
+        free(ct_i_obs_plus_interp);
+        free(tracer_i_tooutput);
+        free(ct_i_tooutput);
+        free(tracer_obs);
+        free(ct_obs);
+        free(p_obs);
+        free(p_i_tmp);
+        free(p_idx);
+        free(p_sort);
+        free(tracer_sort);
+        free(ct_sort);
+        free(p_all);
+        free(p_all_sort);
+        free(p_all_idx);
+        free(i_obs_plus_interp);
+        free(i_surf_and_obs_plus_interp);
+        free(p_obs_plus_interp);
+        free(p_surf_and_obs_plus_interp);
+        free(i_out);
+        free(i_1);
+        free(i_2);
+        free(i_3);
+
+        return;
+}
\ No newline at end of file
diff --git a/gswteos-10.h b/gswteos-10.h
index 9ae0db4..d2c8155 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -230,6 +230,8 @@ DECLSPEC extern void   gsw_rr68_interp_sa_ct(double *sa, double *ct, double *p,
 DECLSPEC extern double gsw_saar(double p, double lon, double lat);
 DECLSPEC extern void   gsw_sa_ct_interp(double *sa, double *ct, double *p,
                 int m, double *p_i, int m_i, double *sa_i, double *ct_i);
+DECLSPEC extern void   gsw_tracer_ct_interp(double *sa, double *ct, double *p,
+                int m, double *p_i, int m_i, double factor, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_sa_freezing_estimate(double p, double saturation_fraction,
                 double *ct, double *t);
 DECLSPEC extern double gsw_sa_freezing_from_ct(double ct, double p,

From 766c5db2f1b1e369fe763818aa3696f7399af6a2 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 14:06:53 -0700
Subject: [PATCH 09/27] gsw_tracer_ct_interp changes

---
 Makefile               |  6 ++++--
 gsw_tracer_ct_interp.c | 22 +++++++++++-----------
 2 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/Makefile b/Makefile
index f3b43f0..2a2a2fb 100644
--- a/Makefile
+++ b/Makefile
@@ -4,10 +4,12 @@ Program         = gsw_check
 PROGRAM_SOURCES = gsw_check_functions.c\
                   gsw_oceanographic_toolbox.c\
                   gsw_saar.c\
-				  gsw_sa_ct_interp.c
+				  gsw_sa_ct_interp.c\
+				  gsw_tracer_ct_interp.c
 LIBRARY_SRCS    = gsw_oceanographic_toolbox.c \
                   gsw_saar.c \
-                  gsw_sa_ct_interp.c
+                  gsw_sa_ct_interp.c\
+				  gsw_tracer_ct_interp.c
 
 # \
 !ifdef 0 #
diff --git a/gsw_tracer_ct_interp.c b/gsw_tracer_ct_interp.c
index 3cf4350..d230a88 100644
--- a/gsw_tracer_ct_interp.c
+++ b/gsw_tracer_ct_interp.c
@@ -3,26 +3,26 @@
 
 /*
 !==========================================================================
-subroutine gsw_tracer_CT_interp (sa,ct,p,p_i)
+subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i)
 !==========================================================================
 !
-! SA and CT interpolation to p_i on a cast
+!  Tracer and CT interpolation to p_i on a cast
 !
 !  tracer   =  tracer                                                 [ ? ]
 !  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
 !  p    =  sea pressure                                            [ dbar ]
 !           ( i.e. absolute pressure - 10.1325 dbar )
-!  p_i  =  specific query points at which the interpolated tracer_i and CT_i
-!            are required                                          [ dbar ]
+!  p_i  =  specific query points at which the interpolated tracer_i
+!          and CT_i are required                                   [ dbar ]
 
 !
-!  tracer_i  =  interpolated SA values at pressures p_i                [ g/kg ]
+!  tracer_i  =  interpolated tracer values at pressures p_i        [ g/kg ]
 !  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
 !
 !--------------------------------------------------------------------------
 */
 void
-gsw_tracer_ct_interp(double *sa, double *ct, double *p, int m,
+gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         double *p_i, int m_i, double factor, double *tracer_i, double *ct_i)
 {
         double  rec_factor  = 1./factor,
@@ -88,9 +88,9 @@ gsw_tracer_ct_interp(double *sa, double *ct, double *p, int m,
         // Find NaNs in profile
         prof_len = 0;
         for (i=0; i<m; ++i) {
-            d = sa[i] + ct[i] + p[i];
+            d = tracer[i] + ct[i] + p[i];
             if (!isnan(d)) {
-                tracer_obs[prof_len] = sa[i];
+                tracer_obs[prof_len] = tracer[i];
                 ct_obs[prof_len] = ct[i];
 
                 p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
@@ -109,7 +109,7 @@ gsw_tracer_ct_interp(double *sa, double *ct, double *p, int m,
         }
 
         // Check if profile pressure values are monotonic
-        // If they are not, then sort pressure values along with SA and CT
+        // If they are not, then sort pressure values along with tracer and CT
         not_monotonic = 0;
         for (i=0; i<prof_len-1; ++i) {
             d = p_obs[i+1] - p_obs[i];
@@ -127,7 +127,7 @@ gsw_tracer_ct_interp(double *sa, double *ct, double *p, int m,
             }
 
             // Once sorted, only save unique pressure values.
-            // SA and CT observations with the same pressure
+            // tracer and CT observations with the same pressure
             // will be averaged.
             new_len = 0;
             unique_count = 1;
@@ -157,7 +157,7 @@ gsw_tracer_ct_interp(double *sa, double *ct, double *p, int m,
             prof_len = new_len;
         }
 
-        // Combine pressure values of observed and interpolated SA and CT
+        // Combine pressure values of observed and interpolated tracer and CT
         p_all_len = prof_len + m_i;
 
         p_all = (double *) malloc(p_all_len*sizeof (double));

From 6be4d4c0e34dda34e9dbc70219e6b28f042c005f Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 14:11:56 -0700
Subject: [PATCH 10/27] Remove unused variables

---
 gsw_tracer_ct_interp.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gsw_tracer_ct_interp.c b/gsw_tracer_ct_interp.c
index d230a88..6cdb9b5 100644
--- a/gsw_tracer_ct_interp.c
+++ b/gsw_tracer_ct_interp.c
@@ -49,7 +49,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         int i, j, k, prof_len, 
             not_monotonic, unique_count, new_len, p_all_len,
             i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
-            i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
+            i_out_len, i_2_len;
         int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
             *i_out, *i_1, *i_2, *i_3;
         double d, ct_f, unique_p, ct_sum, tracer_sum, min_p_obs, max_p_obs;

From 93dc45828fbdbfe3e30cc20ee00c716bde313a32 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 30 Aug 2024 18:04:57 -0700
Subject: [PATCH 11/27] Update comment for gsw_tracer_ct_interp with note on
 factor parameter

---
 gsw_tracer_ct_interp.c | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/gsw_tracer_ct_interp.c b/gsw_tracer_ct_interp.c
index 6cdb9b5..2d9cd03 100644
--- a/gsw_tracer_ct_interp.c
+++ b/gsw_tracer_ct_interp.c
@@ -3,7 +3,7 @@
 
 /*
 !==========================================================================
-subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i)
+subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i,factor)
 !==========================================================================
 !
 !  Tracer and CT interpolation to p_i on a cast
@@ -14,6 +14,8 @@ subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i)
 !           ( i.e. absolute pressure - 10.1325 dbar )
 !  p_i  =  specific query points at which the interpolated tracer_i
 !          and CT_i are required                                   [ dbar ]
+!  factor  =  ratio between the ranges of Conservative Temperature 
+!             and tracer in the world ocean.                          [ ? ]
 
 !
 !  tracer_i  =  interpolated tracer values at pressures p_i        [ g/kg ]

From 33e18540392852d983b1860c5b4dfbf114155f3d Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 6 Sep 2024 12:24:47 -0700
Subject: [PATCH 12/27] Fix finding of index

---
 gsw_sa_ct_interp.c | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
index 7972ed3..14713d9 100644
--- a/gsw_sa_ct_interp.c
+++ b/gsw_sa_ct_interp.c
@@ -307,13 +307,15 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
                 break;
             }
 
-            for (i=0; i < i_obs_plus_interp_len && (i_3[i] - i_frozen) <= 0; ++i) {
-                i_shallower = i;
+            for (i=0; i < i_2_len; ++i) {
+                if ((i_3[i] - i_frozen) <= 0) {
+                    i_shallower = i;
+                }
             }
             i_above = i_2[i_shallower];
             i_above_i = i_3[i_shallower];
-            if ((i_above + 1) > i_3[i_shallower]) {
-                i_below_i = i_3[i_shallower];
+            if ((i_above + 1) > i_3[i_2_len - 1]) {
+                i_below_i = i_3[i_2_len - 1];
             } else {
                 i_below_i = i_3[i_above + 1];
             }

From cd589c2e755fbc1d9f8f539d23280cbaba4000d0 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 6 Sep 2024 12:36:46 -0700
Subject: [PATCH 13/27] Clean up indentations and add newlines

---
 Makefile               | 6 +++---
 gsw_sa_ct_interp.c     | 2 +-
 gsw_tracer_ct_interp.c | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/Makefile b/Makefile
index 2a2a2fb..bd76dd0 100644
--- a/Makefile
+++ b/Makefile
@@ -4,12 +4,12 @@ Program         = gsw_check
 PROGRAM_SOURCES = gsw_check_functions.c\
                   gsw_oceanographic_toolbox.c\
                   gsw_saar.c\
-				  gsw_sa_ct_interp.c\
-				  gsw_tracer_ct_interp.c
+                  gsw_sa_ct_interp.c\
+                  gsw_tracer_ct_interp.c
 LIBRARY_SRCS    = gsw_oceanographic_toolbox.c \
                   gsw_saar.c \
                   gsw_sa_ct_interp.c\
-				  gsw_tracer_ct_interp.c
+                  gsw_tracer_ct_interp.c
 
 # \
 !ifdef 0 #
diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
index 14713d9..dc94dd6 100644
--- a/gsw_sa_ct_interp.c
+++ b/gsw_sa_ct_interp.c
@@ -391,4 +391,4 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         free(i_3);
 
         return;
-}
\ No newline at end of file
+}
diff --git a/gsw_tracer_ct_interp.c b/gsw_tracer_ct_interp.c
index 2d9cd03..9e79b54 100644
--- a/gsw_tracer_ct_interp.c
+++ b/gsw_tracer_ct_interp.c
@@ -333,4 +333,4 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         free(i_3);
 
         return;
-}
\ No newline at end of file
+}

From 5fdd25e4b227c144a0f53cec2284705cc15ad1e7 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 6 Sep 2024 14:29:36 -0700
Subject: [PATCH 14/27] Place gsw_tracer_ct_interp in alphabetical order in
 header file

---
 gswteos-10.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/gswteos-10.h b/gswteos-10.h
index d2c8155..fbfa16f 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -230,8 +230,6 @@ DECLSPEC extern void   gsw_rr68_interp_sa_ct(double *sa, double *ct, double *p,
 DECLSPEC extern double gsw_saar(double p, double lon, double lat);
 DECLSPEC extern void   gsw_sa_ct_interp(double *sa, double *ct, double *p,
                 int m, double *p_i, int m_i, double *sa_i, double *ct_i);
-DECLSPEC extern void   gsw_tracer_ct_interp(double *sa, double *ct, double *p,
-                int m, double *p_i, int m_i, double factor, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_sa_freezing_estimate(double p, double saturation_fraction,
                 double *ct, double *t);
 DECLSPEC extern double gsw_sa_freezing_from_ct(double ct, double p,
@@ -301,6 +299,8 @@ DECLSPEC extern double gsw_t_freezing_poly(double sa, double p,
 DECLSPEC extern double gsw_t_from_ct(double sa, double ct, double p);
 DECLSPEC extern double gsw_t_from_pt0_ice(double pt0_ice, double p);
 DECLSPEC extern double gsw_thermobaric(double sa, double ct, double p);
+DECLSPEC extern void   gsw_tracer_ct_interp(double *sa, double *ct, double *p,
+                int m, double *p_i, int m_i, double factor, double *sa_i, double *ct_i);
 DECLSPEC extern void   gsw_turner_rsubrho(double *sa, double *ct, double *p, int nz,
                 double *tu, double *rsubrho, double *p_mid);
 DECLSPEC extern int    gsw_util_indx(double *x, int n, double z);

From 9c083e4482a59098de106c556847b37f782bc231 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 6 Sep 2024 14:35:13 -0700
Subject: [PATCH 15/27] Move gsw_tracer_ct_interp into main toolbox code

---
 gsw_oceanographic_toolbox.c | 333 +++++++++++++++++++++++++++++++++++
 gsw_tracer_ct_interp.c      | 336 ------------------------------------
 2 files changed, 333 insertions(+), 336 deletions(-)
 delete mode 100644 gsw_tracer_ct_interp.c

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index e783077..d479aec 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -10918,6 +10918,339 @@ gsw_thermobaric(double sa, double ct, double p)
 }
 /*
 !==========================================================================
+subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i,factor)
+!==========================================================================
+!
+!  Tracer and CT interpolation to p_i on a cast
+!
+!  tracer   =  tracer                                                 [ ? ]
+!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
+!  p    =  sea pressure                                            [ dbar ]
+!           ( i.e. absolute pressure - 10.1325 dbar )
+!  p_i  =  specific query points at which the interpolated tracer_i
+!          and CT_i are required                                   [ dbar ]
+!  factor  =  ratio between the ranges of Conservative Temperature 
+!             and tracer in the world ocean.                          [ ? ]
+
+!
+!  tracer_i  =  interpolated tracer values at pressures p_i        [ g/kg ]
+!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
+!
+!--------------------------------------------------------------------------
+*/
+void
+gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
+        double *p_i, int m_i, double factor, double *tracer_i, double *ct_i)
+{
+        double  rec_factor  = 1./factor,
+
+                sin_kpi_on_16[] = {
+                    1.950903220161283e-1,
+                    3.826834323650898e-1,
+                    5.555702330196022e-1,
+                    7.071067811865475e-1,
+                    8.314696123025452e-1,
+                    9.238795325112867e-1,
+                    9.807852804032304e-1
+                },
+                cos_kpi_on_16[] = {
+                    9.807852804032304e-1,
+                    9.238795325112867e-1,
+                    8.314696123025452e-1,
+                    7.071067811865476e-1,
+                    5.555702330196023e-1,
+                    3.826834323650898e-1,
+                    1.950903220161283e-1
+                };
+
+        int i, j, k, prof_len, 
+            not_monotonic, unique_count, new_len, p_all_len,
+            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
+            i_out_len, i_2_len;
+        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
+            *i_out, *i_1, *i_2, *i_3;
+        double d, ct_f, unique_p, ct_sum, tracer_sum, min_p_obs, max_p_obs;
+        double *tracer_obs, *ct_obs, *p_obs, *p_i_tmp,
+               *p_sort, *tracer_sort, *ct_sort, *p_all, *p_all_sort,
+               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
+               *independent_variable, *independent_variable_obs_plus_interp,
+               *scaled_tracer_obs, *v_tmp, *q_tmp, *v_i, *q_i,
+               *tracer_i_obs_plus_interp, *ct_i_obs_plus_interp,
+               *tracer_i_tooutput, *ct_i_tooutput;
+
+        tracer_obs = (double *) malloc(m*sizeof (double));
+        ct_obs = (double *) malloc(m*sizeof (double));
+        p_obs = (double *) malloc(m*sizeof (double));
+        p_idx = (int *) malloc(m*sizeof (int));
+        p_sort = (double *) malloc(m*sizeof (double));
+        tracer_sort = (double *) malloc(m*sizeof (double));
+        ct_sort = (double *) malloc(m*sizeof (double));
+
+        if (m < 4)
+            return; // There must be at least 4 bottles'
+
+        // Check if interpolating pressure is monotonic
+        for (i=0; i<m_i-1; ++i) {
+            d = p_i[i + 1] - p_i[i];
+            if (d < 0.) {
+                return;
+            }
+        }
+
+        for (i=0; i<m_i; i++) {
+            tracer_i[i] = NAN;
+            ct_i[i] = NAN;
+        }
+        
+        // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = tracer[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                tracer_obs[prof_len] = tracer[i];
+                ct_obs[prof_len] = ct[i];
+
+                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
+
+                ++prof_len;
+            }
+        }
+
+        if (prof_len < 2) {
+            return;
+        }
+
+        p_i_tmp = (double *) malloc(m_i*sizeof (double));
+        for (i=0; i<m_i; ++i) {
+            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
+        }
+
+        // Check if profile pressure values are monotonic
+        // If they are not, then sort pressure values along with tracer and CT
+        not_monotonic = 0;
+        for (i=0; i<prof_len-1; ++i) {
+            d = p_obs[i+1] - p_obs[i];
+            if (d < 0.) {
+                ++not_monotonic;
+            }
+        }
+
+        if (not_monotonic > 0) {
+            gsw_util_sort_real(p_obs, prof_len, p_idx);
+            for (i=0; i<prof_len; ++i) {
+                p_sort[i] = p_obs[p_idx[i]];
+                tracer_sort[i] = tracer_obs[p_idx[i]];
+                ct_sort[i] = ct_obs[p_idx[i]];
+            }
+
+            // Once sorted, only save unique pressure values.
+            // tracer and CT observations with the same pressure
+            // will be averaged.
+            new_len = 0;
+            unique_count = 1;
+            unique_p = p_sort[0];
+            tracer_sum = tracer_sort[0];
+            ct_sum = ct_sort[0];
+            for (i=1; i<prof_len; ++i) {
+                if (unique_p == p_sort[i]) {
+                    tracer_sum += tracer_sort[i];
+                    ct_sum += ct_sort[i];
+                    ++unique_count;
+                } else {
+                    p_obs[new_len] = unique_p;
+                    tracer_obs[new_len] = tracer_sum/(double)unique_count;
+                    ct_obs[new_len] = ct_sum/(double)unique_count;
+                    ++new_len;
+                    unique_p = p_sort[i];
+                    tracer_sum = tracer_sort[i];
+                    ct_sum = ct_sort[i];
+                    unique_count = 1;
+                }
+            }
+            p_obs[new_len] = unique_p;
+            tracer_obs[new_len] = tracer_sum/(double)unique_count;
+            ct_obs[new_len] = ct_sum/(double)unique_count;
+            ++new_len;
+            prof_len = new_len;
+        }
+
+        // Combine pressure values of observed and interpolated tracer and CT
+        p_all_len = prof_len + m_i;
+
+        p_all = (double *) malloc(p_all_len*sizeof (double));
+        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
+        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
+
+        memcpy(p_all, p_obs, prof_len*sizeof (double));
+        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
+        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
+        for (i=0; i<p_all_len; ++i) {
+            p_all_sort[i] = p_all[p_all_idx[i]];
+        }
+        new_len = 0;
+        unique_p = p_all_sort[0];
+        for (i=1; i<p_all_len; ++i) {
+            if (unique_p != p_all_sort[i]) {
+                p_all[new_len] = unique_p;
+                ++new_len;
+                unique_p = p_all_sort[i];
+            }
+        }
+        p_all[new_len] = unique_p;
+        ++new_len;
+        p_all_len = new_len;
+
+        i_min_p_obs = 0;
+        min_p_obs = p_obs[0];
+        max_p_obs = p_obs[0];
+        for (i=1; i<prof_len; ++i) {
+            if (p_obs[i] < min_p_obs) {
+                i_min_p_obs = i;
+                min_p_obs = p_obs[i];
+            }
+            max_p_obs = max(max_p_obs, p_obs[i]);
+        }
+
+        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        i_out = (int *) malloc(p_all_len*sizeof (int));
+        i_1 = (int *) malloc(p_all_len*sizeof (int));
+        i_2 = (int *) malloc(p_all_len*sizeof (int));
+        i_3 = (int *) malloc(p_all_len*sizeof (int));
+
+        i_obs_plus_interp_len = 0;
+        i_surf_and_obs_plus_interp_len = 0;
+        for (i=0; i<p_all_len; ++i) {
+            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
+                i_obs_plus_interp[i_obs_plus_interp_len] = i;
+                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
+                i_obs_plus_interp_len++;
+            }
+            if (p_all[i] <= max_p_obs) {
+                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
+                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
+                i_surf_and_obs_plus_interp_len++;
+            }
+        }
+        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
+        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
+
+        independent_variable = (double *) malloc(prof_len*sizeof (double));
+        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            independent_variable[i] = (float)i;
+        }
+        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
+                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+        
+        scaled_tracer_obs = (double *) malloc(prof_len*sizeof (double));
+        v_tmp = (double *) malloc(prof_len*sizeof (double));
+        q_tmp = (double *) malloc(prof_len*sizeof (double));
+        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        tracer_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            scaled_tracer_obs[i] = factor * tracer_obs[i];
+        }
+
+        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
+                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
+        gsw_util_pchip_interp(independent_variable, scaled_tracer_obs, prof_len,
+                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            tracer_i_obs_plus_interp[i] = rec_factor * q_i[i];
+        }
+
+        for(k=0; k<7; ++k) {
+            for(i=0; i<prof_len; ++i) {
+                v_tmp[i] = scaled_tracer_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
+                q_tmp[i] = scaled_tracer_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
+            }
+            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+            for(i=0; i<i_obs_plus_interp_len; ++i) {
+                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
+                tracer_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
+            }
+        }
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            tracer_i_obs_plus_interp[i] *= 0.125;
+            ct_i_obs_plus_interp[i] *= 0.125;
+        }
+
+        tracer_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+
+        if (min_p_obs != 0.) {
+            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i] = NAN;
+                ct_i_tooutput[i] = NAN;
+                if (p_i_tmp[i] < min_p_obs) {
+                    tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i_3[i_min_p_obs]];
+                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
+                }
+            }
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i_obs_plus_interp[i]] = tracer_i_obs_plus_interp[i];
+                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
+            }
+        } else {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i];
+                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
+            }
+
+        }
+
+        for(i=0; i<i_out_len; ++i) {
+            tracer_i[i_out[i]] = tracer_i_tooutput[i_1[i]];
+            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
+        }
+
+        free(independent_variable);
+        free(independent_variable_obs_plus_interp);
+        free(scaled_tracer_obs);
+        free(v_tmp);
+        free(q_tmp);
+        free(v_i);
+        free(q_i);
+        free(tracer_i_obs_plus_interp);
+        free(ct_i_obs_plus_interp);
+        free(tracer_i_tooutput);
+        free(ct_i_tooutput);
+        free(tracer_obs);
+        free(ct_obs);
+        free(p_obs);
+        free(p_i_tmp);
+        free(p_idx);
+        free(p_sort);
+        free(tracer_sort);
+        free(ct_sort);
+        free(p_all);
+        free(p_all_sort);
+        free(p_all_idx);
+        free(i_obs_plus_interp);
+        free(i_surf_and_obs_plus_interp);
+        free(p_obs_plus_interp);
+        free(p_surf_and_obs_plus_interp);
+        free(i_out);
+        free(i_1);
+        free(i_2);
+        free(i_3);
+
+        return;
+}
+/*
+!==========================================================================
 subroutine gsw_turner_rsubrho(sa,ct,p,nz,tu,rsubrho,p_mid)
 !==========================================================================
 
diff --git a/gsw_tracer_ct_interp.c b/gsw_tracer_ct_interp.c
deleted file mode 100644
index 9e79b54..0000000
--- a/gsw_tracer_ct_interp.c
+++ /dev/null
@@ -1,336 +0,0 @@
-#include "gswteos-10.h"
-#include "gsw_internal_const.h"
-
-/*
-!==========================================================================
-subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i,factor)
-!==========================================================================
-!
-!  Tracer and CT interpolation to p_i on a cast
-!
-!  tracer   =  tracer                                                 [ ? ]
-!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
-!  p    =  sea pressure                                            [ dbar ]
-!           ( i.e. absolute pressure - 10.1325 dbar )
-!  p_i  =  specific query points at which the interpolated tracer_i
-!          and CT_i are required                                   [ dbar ]
-!  factor  =  ratio between the ranges of Conservative Temperature 
-!             and tracer in the world ocean.                          [ ? ]
-
-!
-!  tracer_i  =  interpolated tracer values at pressures p_i        [ g/kg ]
-!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
-!
-!--------------------------------------------------------------------------
-*/
-void
-gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
-        double *p_i, int m_i, double factor, double *tracer_i, double *ct_i)
-{
-        double  rec_factor  = 1./factor,
-
-                sin_kpi_on_16[] = {
-                    1.950903220161283e-1,
-                    3.826834323650898e-1,
-                    5.555702330196022e-1,
-                    7.071067811865475e-1,
-                    8.314696123025452e-1,
-                    9.238795325112867e-1,
-                    9.807852804032304e-1
-                },
-                cos_kpi_on_16[] = {
-                    9.807852804032304e-1,
-                    9.238795325112867e-1,
-                    8.314696123025452e-1,
-                    7.071067811865476e-1,
-                    5.555702330196023e-1,
-                    3.826834323650898e-1,
-                    1.950903220161283e-1
-                };
-
-        int i, j, k, prof_len, 
-            not_monotonic, unique_count, new_len, p_all_len,
-            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
-            i_out_len, i_2_len;
-        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
-            *i_out, *i_1, *i_2, *i_3;
-        double d, ct_f, unique_p, ct_sum, tracer_sum, min_p_obs, max_p_obs;
-        double *tracer_obs, *ct_obs, *p_obs, *p_i_tmp,
-               *p_sort, *tracer_sort, *ct_sort, *p_all, *p_all_sort,
-               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
-               *independent_variable, *independent_variable_obs_plus_interp,
-               *scaled_tracer_obs, *v_tmp, *q_tmp, *v_i, *q_i,
-               *tracer_i_obs_plus_interp, *ct_i_obs_plus_interp,
-               *tracer_i_tooutput, *ct_i_tooutput;
-
-        tracer_obs = (double *) malloc(m*sizeof (double));
-        ct_obs = (double *) malloc(m*sizeof (double));
-        p_obs = (double *) malloc(m*sizeof (double));
-        p_idx = (int *) malloc(m*sizeof (int));
-        p_sort = (double *) malloc(m*sizeof (double));
-        tracer_sort = (double *) malloc(m*sizeof (double));
-        ct_sort = (double *) malloc(m*sizeof (double));
-
-        if (m < 4)
-            return; // There must be at least 4 bottles'
-
-        // Check if interpolating pressure is monotonic
-        for (i=0; i<m_i-1; ++i) {
-            d = p_i[i + 1] - p_i[i];
-            if (d < 0.) {
-                return;
-            }
-        }
-
-        for (i=0; i<m_i; i++) {
-            tracer_i[i] = NAN;
-            ct_i[i] = NAN;
-        }
-        
-        // Find NaNs in profile
-        prof_len = 0;
-        for (i=0; i<m; ++i) {
-            d = tracer[i] + ct[i] + p[i];
-            if (!isnan(d)) {
-                tracer_obs[prof_len] = tracer[i];
-                ct_obs[prof_len] = ct[i];
-
-                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
-
-                ++prof_len;
-            }
-        }
-
-        if (prof_len < 2) {
-            return;
-        }
-
-        p_i_tmp = (double *) malloc(m_i*sizeof (double));
-        for (i=0; i<m_i; ++i) {
-            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
-        }
-
-        // Check if profile pressure values are monotonic
-        // If they are not, then sort pressure values along with tracer and CT
-        not_monotonic = 0;
-        for (i=0; i<prof_len-1; ++i) {
-            d = p_obs[i+1] - p_obs[i];
-            if (d < 0.) {
-                ++not_monotonic;
-            }
-        }
-
-        if (not_monotonic > 0) {
-            gsw_util_sort_real(p_obs, prof_len, p_idx);
-            for (i=0; i<prof_len; ++i) {
-                p_sort[i] = p_obs[p_idx[i]];
-                tracer_sort[i] = tracer_obs[p_idx[i]];
-                ct_sort[i] = ct_obs[p_idx[i]];
-            }
-
-            // Once sorted, only save unique pressure values.
-            // tracer and CT observations with the same pressure
-            // will be averaged.
-            new_len = 0;
-            unique_count = 1;
-            unique_p = p_sort[0];
-            tracer_sum = tracer_sort[0];
-            ct_sum = ct_sort[0];
-            for (i=1; i<prof_len; ++i) {
-                if (unique_p == p_sort[i]) {
-                    tracer_sum += tracer_sort[i];
-                    ct_sum += ct_sort[i];
-                    ++unique_count;
-                } else {
-                    p_obs[new_len] = unique_p;
-                    tracer_obs[new_len] = tracer_sum/(double)unique_count;
-                    ct_obs[new_len] = ct_sum/(double)unique_count;
-                    ++new_len;
-                    unique_p = p_sort[i];
-                    tracer_sum = tracer_sort[i];
-                    ct_sum = ct_sort[i];
-                    unique_count = 1;
-                }
-            }
-            p_obs[new_len] = unique_p;
-            tracer_obs[new_len] = tracer_sum/(double)unique_count;
-            ct_obs[new_len] = ct_sum/(double)unique_count;
-            ++new_len;
-            prof_len = new_len;
-        }
-
-        // Combine pressure values of observed and interpolated tracer and CT
-        p_all_len = prof_len + m_i;
-
-        p_all = (double *) malloc(p_all_len*sizeof (double));
-        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
-        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
-
-        memcpy(p_all, p_obs, prof_len*sizeof (double));
-        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
-        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
-        for (i=0; i<p_all_len; ++i) {
-            p_all_sort[i] = p_all[p_all_idx[i]];
-        }
-        new_len = 0;
-        unique_p = p_all_sort[0];
-        for (i=1; i<p_all_len; ++i) {
-            if (unique_p != p_all_sort[i]) {
-                p_all[new_len] = unique_p;
-                ++new_len;
-                unique_p = p_all_sort[i];
-            }
-        }
-        p_all[new_len] = unique_p;
-        ++new_len;
-        p_all_len = new_len;
-
-        i_min_p_obs = 0;
-        min_p_obs = p_obs[0];
-        max_p_obs = p_obs[0];
-        for (i=1; i<prof_len; ++i) {
-            if (p_obs[i] < min_p_obs) {
-                i_min_p_obs = i;
-                min_p_obs = p_obs[i];
-            }
-            max_p_obs = max(max_p_obs, p_obs[i]);
-        }
-
-        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
-        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
-        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
-        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
-        i_out = (int *) malloc(p_all_len*sizeof (int));
-        i_1 = (int *) malloc(p_all_len*sizeof (int));
-        i_2 = (int *) malloc(p_all_len*sizeof (int));
-        i_3 = (int *) malloc(p_all_len*sizeof (int));
-
-        i_obs_plus_interp_len = 0;
-        i_surf_and_obs_plus_interp_len = 0;
-        for (i=0; i<p_all_len; ++i) {
-            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
-                i_obs_plus_interp[i_obs_plus_interp_len] = i;
-                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
-                i_obs_plus_interp_len++;
-            }
-            if (p_all[i] <= max_p_obs) {
-                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
-                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
-                i_surf_and_obs_plus_interp_len++;
-            }
-        }
-        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
-        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
-
-        independent_variable = (double *) malloc(prof_len*sizeof (double));
-        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-
-        for(i=0; i<prof_len; ++i) {
-            independent_variable[i] = (float)i;
-        }
-        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
-                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
-        
-        scaled_tracer_obs = (double *) malloc(prof_len*sizeof (double));
-        v_tmp = (double *) malloc(prof_len*sizeof (double));
-        q_tmp = (double *) malloc(prof_len*sizeof (double));
-        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-        tracer_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-
-        for(i=0; i<prof_len; ++i) {
-            scaled_tracer_obs[i] = factor * tracer_obs[i];
-        }
-
-        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
-                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
-        gsw_util_pchip_interp(independent_variable, scaled_tracer_obs, prof_len,
-                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
-
-        for(i=0; i<i_obs_plus_interp_len; ++i) {
-            tracer_i_obs_plus_interp[i] = rec_factor * q_i[i];
-        }
-
-        for(k=0; k<7; ++k) {
-            for(i=0; i<prof_len; ++i) {
-                v_tmp[i] = scaled_tracer_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
-                q_tmp[i] = scaled_tracer_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
-            }
-            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
-                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
-            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
-                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
-            for(i=0; i<i_obs_plus_interp_len; ++i) {
-                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
-                tracer_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
-            }
-        }
-
-        for(i=0; i<i_obs_plus_interp_len; ++i) {
-            tracer_i_obs_plus_interp[i] *= 0.125;
-            ct_i_obs_plus_interp[i] *= 0.125;
-        }
-
-        tracer_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
-        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
-
-        if (min_p_obs != 0.) {
-            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
-                tracer_i_tooutput[i] = NAN;
-                ct_i_tooutput[i] = NAN;
-                if (p_i_tmp[i] < min_p_obs) {
-                    tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i_3[i_min_p_obs]];
-                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
-                }
-            }
-            for (i=0; i<i_obs_plus_interp_len; ++i) {
-                tracer_i_tooutput[i_obs_plus_interp[i]] = tracer_i_obs_plus_interp[i];
-                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
-            }
-        } else {
-            for (i=0; i<i_obs_plus_interp_len; ++i) {
-                tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i];
-                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
-            }
-
-        }
-
-        for(i=0; i<i_out_len; ++i) {
-            tracer_i[i_out[i]] = tracer_i_tooutput[i_1[i]];
-            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
-        }
-
-        free(independent_variable);
-        free(independent_variable_obs_plus_interp);
-        free(scaled_tracer_obs);
-        free(v_tmp);
-        free(q_tmp);
-        free(v_i);
-        free(q_i);
-        free(tracer_i_obs_plus_interp);
-        free(ct_i_obs_plus_interp);
-        free(tracer_i_tooutput);
-        free(ct_i_tooutput);
-        free(tracer_obs);
-        free(ct_obs);
-        free(p_obs);
-        free(p_i_tmp);
-        free(p_idx);
-        free(p_sort);
-        free(tracer_sort);
-        free(ct_sort);
-        free(p_all);
-        free(p_all_sort);
-        free(p_all_idx);
-        free(i_obs_plus_interp);
-        free(i_surf_and_obs_plus_interp);
-        free(p_obs_plus_interp);
-        free(p_surf_and_obs_plus_interp);
-        free(i_out);
-        free(i_1);
-        free(i_2);
-        free(i_3);
-
-        return;
-}

From 72abf4b0e5b6540dcfc97eb69d5bb6d6fd611aa7 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 6 Sep 2024 14:38:38 -0700
Subject: [PATCH 16/27] Move gsw_sa_ct_interp into main toolbox code

---
 Makefile                    |   8 +-
 gsw_oceanographic_toolbox.c | 391 +++++++++++++++++++++++++++++++++++
 gsw_sa_ct_interp.c          | 394 ------------------------------------
 3 files changed, 393 insertions(+), 400 deletions(-)
 delete mode 100644 gsw_sa_ct_interp.c

diff --git a/Makefile b/Makefile
index bd76dd0..d6740a6 100644
--- a/Makefile
+++ b/Makefile
@@ -3,13 +3,9 @@ Program         = gsw_check
 
 PROGRAM_SOURCES = gsw_check_functions.c\
                   gsw_oceanographic_toolbox.c\
-                  gsw_saar.c\
-                  gsw_sa_ct_interp.c\
-                  gsw_tracer_ct_interp.c
+                  gsw_saar.c
 LIBRARY_SRCS    = gsw_oceanographic_toolbox.c \
-                  gsw_saar.c \
-                  gsw_sa_ct_interp.c\
-                  gsw_tracer_ct_interp.c
+                  gsw_saar.c 
 
 # \
 !ifdef 0 #
diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index d479aec..829306e 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -8589,6 +8589,397 @@ rr68_interp_section(int sectnum, double *sa, double *ct, double *p, int mp,
 }
 /*
 !==========================================================================
+subroutine gsw_SA_CT_interp (sa,ct,p,p_i)
+!==========================================================================
+!
+! SA and CT interpolation to p_i on a cast
+!
+!  SA   =  Absolute Salinity                                       [ g/kg ]
+!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
+!  p    =  sea pressure                                            [ dbar ]
+!           ( i.e. absolute pressure - 10.1325 dbar )
+!  p_i  =  specific query points at which the interpolated SA_i and CT_i
+!            are required                                          [ dbar ]
+
+!
+!  SA_i  =  interpolated SA values at pressures p_i                [ g/kg ]
+!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
+!
+!--------------------------------------------------------------------------
+*/
+void
+gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
+        double *p_i, int m_i, double *sa_i, double *ct_i)
+{
+        double  factor      = 9.,
+                rec_factor  = 1./factor,
+
+                sin_kpi_on_16[] = {
+                    1.950903220161283e-1,
+                    3.826834323650898e-1,
+                    5.555702330196022e-1,
+                    7.071067811865475e-1,
+                    8.314696123025452e-1,
+                    9.238795325112867e-1,
+                    9.807852804032304e-1
+                },
+                cos_kpi_on_16[] = {
+                    9.807852804032304e-1,
+                    9.238795325112867e-1,
+                    8.314696123025452e-1,
+                    7.071067811865476e-1,
+                    5.555702330196023e-1,
+                    3.826834323650898e-1,
+                    1.950903220161283e-1
+                };
+
+        int i, j, k, prof_len, 
+            not_monotonic, unique_count, new_len, p_all_len,
+            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
+            i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
+        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
+            *i_out, *i_1, *i_2, *i_3;
+        double d, ct_f, unique_p, ct_sum, sa_sum, min_p_obs, max_p_obs;
+        double *sa_obs, *ct_obs, *p_obs, *p_i_tmp,
+               *p_sort, *sa_sort, *ct_sort, *p_all, *p_all_sort,
+               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
+               *independent_variable, *independent_variable_obs_plus_interp,
+               *scaled_sa_obs, *v_tmp, *q_tmp, *v_i, *q_i,
+               *sa_i_obs_plus_interp, *ct_i_obs_plus_interp,
+               *sa_i_limiting_obs_plus_interp, *ct_i_limiting_obs_plus_interp,
+               *ctf_i_tointerp, *sa_i_tooutput, *ct_i_tooutput;
+
+        sa_obs = (double *) malloc(m*sizeof (double));
+        ct_obs = (double *) malloc(m*sizeof (double));
+        p_obs = (double *) malloc(m*sizeof (double));
+        p_idx = (int *) malloc(m*sizeof (int));
+        p_sort = (double *) malloc(m*sizeof (double));
+        sa_sort = (double *) malloc(m*sizeof (double));
+        ct_sort = (double *) malloc(m*sizeof (double));
+
+        if (m < 4)
+            return; // There must be at least 4 bottles'
+
+        // Check if interpolating pressure is monotonic
+        for (i=0; i<m_i-1; ++i) {
+            d = p_i[i + 1] - p_i[i];
+            if (d < 0.) {
+                return;
+            }
+        }
+
+        for (i=0; i<m_i; i++) {
+            sa_i[i] = NAN;
+            ct_i[i] = NAN;
+        }
+        
+        // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = sa[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                sa_obs[prof_len] = sa[i];
+                ct_obs[prof_len] = ct[i];
+
+                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
+
+                ct_f = gsw_ct_freezing_poly(sa_obs[prof_len],
+                                            p_obs[prof_len], 
+                                            0.);
+                if (ct_obs[prof_len] < (ct_f - 0.1)) {
+                    ct_obs[prof_len] = ct_f;
+                }
+
+                ++prof_len;
+            }
+        }
+
+        if (prof_len < 2) {
+            return;
+        }
+
+        p_i_tmp = (double *) malloc(m_i*sizeof (double));
+        for (i=0; i<m_i; ++i) {
+            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
+        }
+
+        // Check if profile pressure values are monotonic
+        // If they are not, then sort pressure values along with SA and CT
+        not_monotonic = 0;
+        for (i=0; i<prof_len-1; ++i) {
+            d = p_obs[i+1] - p_obs[i];
+            if (d < 0.) {
+                ++not_monotonic;
+            }
+        }
+
+        if (not_monotonic > 0) {
+            gsw_util_sort_real(p_obs, prof_len, p_idx);
+            for (i=0; i<prof_len; ++i) {
+                p_sort[i] = p_obs[p_idx[i]];
+                sa_sort[i] = sa_obs[p_idx[i]];
+                ct_sort[i] = ct_obs[p_idx[i]];
+            }
+
+            // Once sorted, only save unique pressure values.
+            // SA and CT observations with the same pressure
+            // will be averaged.
+            new_len = 0;
+            unique_count = 1;
+            unique_p = p_sort[0];
+            sa_sum = sa_sort[0];
+            ct_sum = ct_sort[0];
+            for (i=1; i<prof_len; ++i) {
+                if (unique_p == p_sort[i]) {
+                    sa_sum += sa_sort[i];
+                    ct_sum += ct_sort[i];
+                    ++unique_count;
+                } else {
+                    p_obs[new_len] = unique_p;
+                    sa_obs[new_len] = sa_sum/(double)unique_count;
+                    ct_obs[new_len] = ct_sum/(double)unique_count;
+                    ++new_len;
+                    unique_p = p_sort[i];
+                    sa_sum = sa_sort[i];
+                    ct_sum = ct_sort[i];
+                    unique_count = 1;
+                }
+            }
+            p_obs[new_len] = unique_p;
+            sa_obs[new_len] = sa_sum/(double)unique_count;
+            ct_obs[new_len] = ct_sum/(double)unique_count;
+            ++new_len;
+            prof_len = new_len;
+        }
+
+        // Combine pressure values of observed and interpolated SA and CT
+        p_all_len = prof_len + m_i;
+
+        p_all = (double *) malloc(p_all_len*sizeof (double));
+        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
+        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
+
+        memcpy(p_all, p_obs, prof_len*sizeof (double));
+        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
+        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
+        for (i=0; i<p_all_len; ++i) {
+            p_all_sort[i] = p_all[p_all_idx[i]];
+        }
+        new_len = 0;
+        unique_p = p_all_sort[0];
+        for (i=1; i<p_all_len; ++i) {
+            if (unique_p != p_all_sort[i]) {
+                p_all[new_len] = unique_p;
+                ++new_len;
+                unique_p = p_all_sort[i];
+            }
+        }
+        p_all[new_len] = unique_p;
+        ++new_len;
+        p_all_len = new_len;
+
+        i_min_p_obs = 0;
+        min_p_obs = p_obs[0];
+        max_p_obs = p_obs[0];
+        for (i=1; i<prof_len; ++i) {
+            if (p_obs[i] < min_p_obs) {
+                i_min_p_obs = i;
+                min_p_obs = p_obs[i];
+            }
+            max_p_obs = max(max_p_obs, p_obs[i]);
+        }
+
+        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        i_out = (int *) malloc(p_all_len*sizeof (int));
+        i_1 = (int *) malloc(p_all_len*sizeof (int));
+        i_2 = (int *) malloc(p_all_len*sizeof (int));
+        i_3 = (int *) malloc(p_all_len*sizeof (int));
+
+        i_obs_plus_interp_len = 0;
+        i_surf_and_obs_plus_interp_len = 0;
+        for (i=0; i<p_all_len; ++i) {
+            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
+                i_obs_plus_interp[i_obs_plus_interp_len] = i;
+                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
+                i_obs_plus_interp_len++;
+            }
+            if (p_all[i] <= max_p_obs) {
+                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
+                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
+                i_surf_and_obs_plus_interp_len++;
+            }
+        }
+        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
+        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
+
+        independent_variable = (double *) malloc(prof_len*sizeof (double));
+        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            independent_variable[i] = (float)i;
+        }
+        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
+                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+        
+        scaled_sa_obs = (double *) malloc(prof_len*sizeof (double));
+        v_tmp = (double *) malloc(prof_len*sizeof (double));
+        q_tmp = (double *) malloc(prof_len*sizeof (double));
+        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        sa_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        sa_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            scaled_sa_obs[i] = factor * sa_obs[i];
+        }
+
+        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
+                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
+        gsw_util_pchip_interp(independent_variable, scaled_sa_obs, prof_len,
+                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            sa_i_obs_plus_interp[i] = rec_factor * q_i[i];
+        }
+
+        for(k=0; k<7; ++k) {
+            for(i=0; i<prof_len; ++i) {
+                v_tmp[i] = scaled_sa_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
+                q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
+            }
+            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+            for(i=0; i<i_obs_plus_interp_len; ++i) {
+                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
+                sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
+            }
+        }
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            sa_i_obs_plus_interp[i] *= 0.125;
+            ct_i_obs_plus_interp[i] *= 0.125;
+        }
+
+        gsw_linear_interp_sa_ct(sa_obs, ct_obs, independent_variable, prof_len,
+                                independent_variable_obs_plus_interp, i_obs_plus_interp_len,
+                                sa_i_limiting_obs_plus_interp, ct_i_limiting_obs_plus_interp);
+
+        ctf_i_tointerp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        while (1) {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                d = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_obs_plus_interp[i], 0.);
+                if (ct_i_limiting_obs_plus_interp[i] < (d - 0.1)) {
+                    ctf_i_tointerp[i] = ct_i_limiting_obs_plus_interp[i];
+                } else {
+                    ctf_i_tointerp[i] = d;
+                }
+            }
+
+            for (i_frozen=0; i_frozen<i_obs_plus_interp_len; ++i_frozen) {
+                if (ct_i_obs_plus_interp[i_frozen] < (ctf_i_tointerp[i_frozen] - 0.1)) {
+                    break;
+                }
+            }
+
+            if (i_frozen == i_obs_plus_interp_len) {
+                break;
+            }
+
+            for (i=0; i < i_2_len; ++i) {
+                if ((i_3[i] - i_frozen) <= 0) {
+                    i_shallower = i;
+                }
+            }
+            i_above = i_2[i_shallower];
+            i_above_i = i_3[i_shallower];
+            if ((i_above + 1) > i_3[i_2_len - 1]) {
+                i_below_i = i_3[i_2_len - 1];
+            } else {
+                i_below_i = i_3[i_above + 1];
+            }
+
+            for (i=i_above_i; i<=i_below_i; ++i) {
+                sa_i_obs_plus_interp[i] = sa_i_limiting_obs_plus_interp[i];
+                ct_i_obs_plus_interp[i] = ct_i_limiting_obs_plus_interp[i];
+                ctf_i_tointerp[i] = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_all[i], 0.);
+            }
+        }
+
+        sa_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+
+        if (min_p_obs != 0.) {
+            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i] = NAN;
+                ct_i_tooutput[i] = NAN;
+                if (p_i_tmp[i] < min_p_obs) {
+                    sa_i_tooutput[i] = sa_i_obs_plus_interp[i_3[i_min_p_obs]];
+                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
+                }
+            }
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i_obs_plus_interp[i]] = sa_i_obs_plus_interp[i];
+                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
+            }
+        } else {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i] = sa_i_obs_plus_interp[i];
+                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
+            }
+
+        }
+
+        for(i=0; i<i_out_len; ++i) {
+            sa_i[i_out[i]] = sa_i_tooutput[i_1[i]];
+            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
+        }
+
+        free(independent_variable);
+        free(independent_variable_obs_plus_interp);
+        free(scaled_sa_obs);
+        free(v_tmp);
+        free(q_tmp);
+        free(v_i);
+        free(q_i);
+        free(sa_i_obs_plus_interp);
+        free(ct_i_obs_plus_interp);
+        free(sa_i_limiting_obs_plus_interp);
+        free(ct_i_limiting_obs_plus_interp);
+        free(ctf_i_tointerp);
+        free(sa_i_tooutput);
+        free(ct_i_tooutput);
+        free(sa_obs);
+        free(ct_obs);
+        free(p_obs);
+        free(p_i_tmp);
+        free(p_idx);
+        free(p_sort);
+        free(sa_sort);
+        free(ct_sort);
+        free(p_all);
+        free(p_all_sort);
+        free(p_all_idx);
+        free(i_obs_plus_interp);
+        free(i_surf_and_obs_plus_interp);
+        free(p_obs_plus_interp);
+        free(p_surf_and_obs_plus_interp);
+        free(i_out);
+        free(i_1);
+        free(i_2);
+        free(i_3);
+
+        return;
+}
+/*
+!==========================================================================
 elemental function gsw_sa_freezing_estimate (p, saturation_fraction, ct, t)
 !==========================================================================
 !
diff --git a/gsw_sa_ct_interp.c b/gsw_sa_ct_interp.c
deleted file mode 100644
index dc94dd6..0000000
--- a/gsw_sa_ct_interp.c
+++ /dev/null
@@ -1,394 +0,0 @@
-#include "gswteos-10.h"
-#include "gsw_internal_const.h"
-
-/*
-!==========================================================================
-subroutine gsw_SA_CT_interp (sa,ct,p,p_i)
-!==========================================================================
-!
-! SA and CT interpolation to p_i on a cast
-!
-!  SA   =  Absolute Salinity                                       [ g/kg ]
-!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
-!  p    =  sea pressure                                            [ dbar ]
-!           ( i.e. absolute pressure - 10.1325 dbar )
-!  p_i  =  specific query points at which the interpolated SA_i and CT_i
-!            are required                                          [ dbar ]
-
-!
-!  SA_i  =  interpolated SA values at pressures p_i                [ g/kg ]
-!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
-!
-!--------------------------------------------------------------------------
-*/
-void
-gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
-        double *p_i, int m_i, double *sa_i, double *ct_i)
-{
-        double  factor      = 9.,
-                rec_factor  = 1./factor,
-
-                sin_kpi_on_16[] = {
-                    1.950903220161283e-1,
-                    3.826834323650898e-1,
-                    5.555702330196022e-1,
-                    7.071067811865475e-1,
-                    8.314696123025452e-1,
-                    9.238795325112867e-1,
-                    9.807852804032304e-1
-                },
-                cos_kpi_on_16[] = {
-                    9.807852804032304e-1,
-                    9.238795325112867e-1,
-                    8.314696123025452e-1,
-                    7.071067811865476e-1,
-                    5.555702330196023e-1,
-                    3.826834323650898e-1,
-                    1.950903220161283e-1
-                };
-
-        int i, j, k, prof_len, 
-            not_monotonic, unique_count, new_len, p_all_len,
-            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
-            i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
-        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
-            *i_out, *i_1, *i_2, *i_3;
-        double d, ct_f, unique_p, ct_sum, sa_sum, min_p_obs, max_p_obs;
-        double *sa_obs, *ct_obs, *p_obs, *p_i_tmp,
-               *p_sort, *sa_sort, *ct_sort, *p_all, *p_all_sort,
-               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
-               *independent_variable, *independent_variable_obs_plus_interp,
-               *scaled_sa_obs, *v_tmp, *q_tmp, *v_i, *q_i,
-               *sa_i_obs_plus_interp, *ct_i_obs_plus_interp,
-               *sa_i_limiting_obs_plus_interp, *ct_i_limiting_obs_plus_interp,
-               *ctf_i_tointerp, *sa_i_tooutput, *ct_i_tooutput;
-
-        sa_obs = (double *) malloc(m*sizeof (double));
-        ct_obs = (double *) malloc(m*sizeof (double));
-        p_obs = (double *) malloc(m*sizeof (double));
-        p_idx = (int *) malloc(m*sizeof (int));
-        p_sort = (double *) malloc(m*sizeof (double));
-        sa_sort = (double *) malloc(m*sizeof (double));
-        ct_sort = (double *) malloc(m*sizeof (double));
-
-        if (m < 4)
-            return; // There must be at least 4 bottles'
-
-        // Check if interpolating pressure is monotonic
-        for (i=0; i<m_i-1; ++i) {
-            d = p_i[i + 1] - p_i[i];
-            if (d < 0.) {
-                return;
-            }
-        }
-
-        for (i=0; i<m_i; i++) {
-            sa_i[i] = NAN;
-            ct_i[i] = NAN;
-        }
-        
-        // Find NaNs in profile
-        prof_len = 0;
-        for (i=0; i<m; ++i) {
-            d = sa[i] + ct[i] + p[i];
-            if (!isnan(d)) {
-                sa_obs[prof_len] = sa[i];
-                ct_obs[prof_len] = ct[i];
-
-                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
-
-                ct_f = gsw_ct_freezing_poly(sa_obs[prof_len],
-                                            p_obs[prof_len], 
-                                            0.);
-                if (ct_obs[prof_len] < (ct_f - 0.1)) {
-                    ct_obs[prof_len] = ct_f;
-                }
-
-                ++prof_len;
-            }
-        }
-
-        if (prof_len < 2) {
-            return;
-        }
-
-        p_i_tmp = (double *) malloc(m_i*sizeof (double));
-        for (i=0; i<m_i; ++i) {
-            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
-        }
-
-        // Check if profile pressure values are monotonic
-        // If they are not, then sort pressure values along with SA and CT
-        not_monotonic = 0;
-        for (i=0; i<prof_len-1; ++i) {
-            d = p_obs[i+1] - p_obs[i];
-            if (d < 0.) {
-                ++not_monotonic;
-            }
-        }
-
-        if (not_monotonic > 0) {
-            gsw_util_sort_real(p_obs, prof_len, p_idx);
-            for (i=0; i<prof_len; ++i) {
-                p_sort[i] = p_obs[p_idx[i]];
-                sa_sort[i] = sa_obs[p_idx[i]];
-                ct_sort[i] = ct_obs[p_idx[i]];
-            }
-
-            // Once sorted, only save unique pressure values.
-            // SA and CT observations with the same pressure
-            // will be averaged.
-            new_len = 0;
-            unique_count = 1;
-            unique_p = p_sort[0];
-            sa_sum = sa_sort[0];
-            ct_sum = ct_sort[0];
-            for (i=1; i<prof_len; ++i) {
-                if (unique_p == p_sort[i]) {
-                    sa_sum += sa_sort[i];
-                    ct_sum += ct_sort[i];
-                    ++unique_count;
-                } else {
-                    p_obs[new_len] = unique_p;
-                    sa_obs[new_len] = sa_sum/(double)unique_count;
-                    ct_obs[new_len] = ct_sum/(double)unique_count;
-                    ++new_len;
-                    unique_p = p_sort[i];
-                    sa_sum = sa_sort[i];
-                    ct_sum = ct_sort[i];
-                    unique_count = 1;
-                }
-            }
-            p_obs[new_len] = unique_p;
-            sa_obs[new_len] = sa_sum/(double)unique_count;
-            ct_obs[new_len] = ct_sum/(double)unique_count;
-            ++new_len;
-            prof_len = new_len;
-        }
-
-        // Combine pressure values of observed and interpolated SA and CT
-        p_all_len = prof_len + m_i;
-
-        p_all = (double *) malloc(p_all_len*sizeof (double));
-        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
-        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
-
-        memcpy(p_all, p_obs, prof_len*sizeof (double));
-        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
-        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
-        for (i=0; i<p_all_len; ++i) {
-            p_all_sort[i] = p_all[p_all_idx[i]];
-        }
-        new_len = 0;
-        unique_p = p_all_sort[0];
-        for (i=1; i<p_all_len; ++i) {
-            if (unique_p != p_all_sort[i]) {
-                p_all[new_len] = unique_p;
-                ++new_len;
-                unique_p = p_all_sort[i];
-            }
-        }
-        p_all[new_len] = unique_p;
-        ++new_len;
-        p_all_len = new_len;
-
-        i_min_p_obs = 0;
-        min_p_obs = p_obs[0];
-        max_p_obs = p_obs[0];
-        for (i=1; i<prof_len; ++i) {
-            if (p_obs[i] < min_p_obs) {
-                i_min_p_obs = i;
-                min_p_obs = p_obs[i];
-            }
-            max_p_obs = max(max_p_obs, p_obs[i]);
-        }
-
-        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
-        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
-        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
-        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
-        i_out = (int *) malloc(p_all_len*sizeof (int));
-        i_1 = (int *) malloc(p_all_len*sizeof (int));
-        i_2 = (int *) malloc(p_all_len*sizeof (int));
-        i_3 = (int *) malloc(p_all_len*sizeof (int));
-
-        i_obs_plus_interp_len = 0;
-        i_surf_and_obs_plus_interp_len = 0;
-        for (i=0; i<p_all_len; ++i) {
-            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
-                i_obs_plus_interp[i_obs_plus_interp_len] = i;
-                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
-                i_obs_plus_interp_len++;
-            }
-            if (p_all[i] <= max_p_obs) {
-                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
-                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
-                i_surf_and_obs_plus_interp_len++;
-            }
-        }
-        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
-        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
-
-        independent_variable = (double *) malloc(prof_len*sizeof (double));
-        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-
-        for(i=0; i<prof_len; ++i) {
-            independent_variable[i] = (float)i;
-        }
-        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
-                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
-        
-        scaled_sa_obs = (double *) malloc(prof_len*sizeof (double));
-        v_tmp = (double *) malloc(prof_len*sizeof (double));
-        q_tmp = (double *) malloc(prof_len*sizeof (double));
-        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
-        sa_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-        sa_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-        ct_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-
-        for(i=0; i<prof_len; ++i) {
-            scaled_sa_obs[i] = factor * sa_obs[i];
-        }
-
-        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
-                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
-        gsw_util_pchip_interp(independent_variable, scaled_sa_obs, prof_len,
-                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
-
-        for(i=0; i<i_obs_plus_interp_len; ++i) {
-            sa_i_obs_plus_interp[i] = rec_factor * q_i[i];
-        }
-
-        for(k=0; k<7; ++k) {
-            for(i=0; i<prof_len; ++i) {
-                v_tmp[i] = scaled_sa_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
-                q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
-            }
-            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
-                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
-            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
-                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
-            for(i=0; i<i_obs_plus_interp_len; ++i) {
-                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
-                sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
-            }
-        }
-
-        for(i=0; i<i_obs_plus_interp_len; ++i) {
-            sa_i_obs_plus_interp[i] *= 0.125;
-            ct_i_obs_plus_interp[i] *= 0.125;
-        }
-
-        gsw_linear_interp_sa_ct(sa_obs, ct_obs, independent_variable, prof_len,
-                                independent_variable_obs_plus_interp, i_obs_plus_interp_len,
-                                sa_i_limiting_obs_plus_interp, ct_i_limiting_obs_plus_interp);
-
-        ctf_i_tointerp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
-
-        while (1) {
-            for (i=0; i<i_obs_plus_interp_len; ++i) {
-                d = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_obs_plus_interp[i], 0.);
-                if (ct_i_limiting_obs_plus_interp[i] < (d - 0.1)) {
-                    ctf_i_tointerp[i] = ct_i_limiting_obs_plus_interp[i];
-                } else {
-                    ctf_i_tointerp[i] = d;
-                }
-            }
-
-            for (i_frozen=0; i_frozen<i_obs_plus_interp_len; ++i_frozen) {
-                if (ct_i_obs_plus_interp[i_frozen] < (ctf_i_tointerp[i_frozen] - 0.1)) {
-                    break;
-                }
-            }
-
-            if (i_frozen == i_obs_plus_interp_len) {
-                break;
-            }
-
-            for (i=0; i < i_2_len; ++i) {
-                if ((i_3[i] - i_frozen) <= 0) {
-                    i_shallower = i;
-                }
-            }
-            i_above = i_2[i_shallower];
-            i_above_i = i_3[i_shallower];
-            if ((i_above + 1) > i_3[i_2_len - 1]) {
-                i_below_i = i_3[i_2_len - 1];
-            } else {
-                i_below_i = i_3[i_above + 1];
-            }
-
-            for (i=i_above_i; i<=i_below_i; ++i) {
-                sa_i_obs_plus_interp[i] = sa_i_limiting_obs_plus_interp[i];
-                ct_i_obs_plus_interp[i] = ct_i_limiting_obs_plus_interp[i];
-                ctf_i_tointerp[i] = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_all[i], 0.);
-            }
-        }
-
-        sa_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
-        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
-
-        if (min_p_obs != 0.) {
-            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
-                sa_i_tooutput[i] = NAN;
-                ct_i_tooutput[i] = NAN;
-                if (p_i_tmp[i] < min_p_obs) {
-                    sa_i_tooutput[i] = sa_i_obs_plus_interp[i_3[i_min_p_obs]];
-                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
-                }
-            }
-            for (i=0; i<i_obs_plus_interp_len; ++i) {
-                sa_i_tooutput[i_obs_plus_interp[i]] = sa_i_obs_plus_interp[i];
-                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
-            }
-        } else {
-            for (i=0; i<i_obs_plus_interp_len; ++i) {
-                sa_i_tooutput[i] = sa_i_obs_plus_interp[i];
-                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
-            }
-
-        }
-
-        for(i=0; i<i_out_len; ++i) {
-            sa_i[i_out[i]] = sa_i_tooutput[i_1[i]];
-            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
-        }
-
-        free(independent_variable);
-        free(independent_variable_obs_plus_interp);
-        free(scaled_sa_obs);
-        free(v_tmp);
-        free(q_tmp);
-        free(v_i);
-        free(q_i);
-        free(sa_i_obs_plus_interp);
-        free(ct_i_obs_plus_interp);
-        free(sa_i_limiting_obs_plus_interp);
-        free(ct_i_limiting_obs_plus_interp);
-        free(ctf_i_tointerp);
-        free(sa_i_tooutput);
-        free(ct_i_tooutput);
-        free(sa_obs);
-        free(ct_obs);
-        free(p_obs);
-        free(p_i_tmp);
-        free(p_idx);
-        free(p_sort);
-        free(sa_sort);
-        free(ct_sort);
-        free(p_all);
-        free(p_all_sort);
-        free(p_all_idx);
-        free(i_obs_plus_interp);
-        free(i_surf_and_obs_plus_interp);
-        free(p_obs_plus_interp);
-        free(p_surf_and_obs_plus_interp);
-        free(i_out);
-        free(i_1);
-        free(i_2);
-        free(i_3);
-
-        return;
-}

From a858ed0d5ad15338503c5fd94798e8e1472e287b Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 6 Sep 2024 17:45:22 -0700
Subject: [PATCH 17/27] Make SA-CT interp functions return error code

---
 gsw_oceanographic_toolbox.c | 20 ++++++++++----------
 gswteos-10.h                |  4 ++--
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index 829306e..45590c4 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -8607,7 +8607,7 @@ subroutine gsw_SA_CT_interp (sa,ct,p,p_i)
 !
 !--------------------------------------------------------------------------
 */
-void
+int
 gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         double *p_i, int m_i, double *sa_i, double *ct_i)
 {
@@ -8658,13 +8658,13 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         ct_sort = (double *) malloc(m*sizeof (double));
 
         if (m < 4)
-            return; // There must be at least 4 bottles'
+            return 1; // There must be at least 4 bottles'
 
         // Check if interpolating pressure is monotonic
         for (i=0; i<m_i-1; ++i) {
             d = p_i[i + 1] - p_i[i];
             if (d < 0.) {
-                return;
+                return 1;
             }
         }
 
@@ -8695,7 +8695,7 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         }
 
         if (prof_len < 2) {
-            return;
+            return 1;
         }
 
         p_i_tmp = (double *) malloc(m_i*sizeof (double));
@@ -8976,7 +8976,7 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         free(i_2);
         free(i_3);
 
-        return;
+        return 0;
 }
 /*
 !==========================================================================
@@ -11329,7 +11329,7 @@ subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i,factor)
 !
 !--------------------------------------------------------------------------
 */
-void
+int
 gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         double *p_i, int m_i, double factor, double *tracer_i, double *ct_i)
 {
@@ -11378,13 +11378,13 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         ct_sort = (double *) malloc(m*sizeof (double));
 
         if (m < 4)
-            return; // There must be at least 4 bottles'
+            return 1; // There must be at least 4 bottles'
 
         // Check if interpolating pressure is monotonic
         for (i=0; i<m_i-1; ++i) {
             d = p_i[i + 1] - p_i[i];
             if (d < 0.) {
-                return;
+                return 1;
             }
         }
 
@@ -11408,7 +11408,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         }
 
         if (prof_len < 2) {
-            return;
+            return 1;
         }
 
         p_i_tmp = (double *) malloc(m_i*sizeof (double));
@@ -11638,7 +11638,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         free(i_2);
         free(i_3);
 
-        return;
+        return 0;
 }
 /*
 !==========================================================================
diff --git a/gswteos-10.h b/gswteos-10.h
index fbfa16f..74be239 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -228,7 +228,7 @@ DECLSPEC extern double gsw_rho_t_exact(double sa, double t, double p);
 DECLSPEC extern void   gsw_rr68_interp_sa_ct(double *sa, double *ct, double *p, int mp,
                 double *p_i, int mp_i, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_saar(double p, double lon, double lat);
-DECLSPEC extern void   gsw_sa_ct_interp(double *sa, double *ct, double *p,
+DECLSPEC extern int    gsw_sa_ct_interp(double *sa, double *ct, double *p,
                 int m, double *p_i, int m_i, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_sa_freezing_estimate(double p, double saturation_fraction,
                 double *ct, double *t);
@@ -299,7 +299,7 @@ DECLSPEC extern double gsw_t_freezing_poly(double sa, double p,
 DECLSPEC extern double gsw_t_from_ct(double sa, double ct, double p);
 DECLSPEC extern double gsw_t_from_pt0_ice(double pt0_ice, double p);
 DECLSPEC extern double gsw_thermobaric(double sa, double ct, double p);
-DECLSPEC extern void   gsw_tracer_ct_interp(double *sa, double *ct, double *p,
+DECLSPEC extern int    gsw_tracer_ct_interp(double *sa, double *ct, double *p,
                 int m, double *p_i, int m_i, double factor, double *sa_i, double *ct_i);
 DECLSPEC extern void   gsw_turner_rsubrho(double *sa, double *ct, double *p, int nz,
                 double *tu, double *rsubrho, double *p_mid);

From 994dba7d865ac5fd662efc94528723afae03d8bc Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Thu, 12 Sep 2024 12:45:03 -0700
Subject: [PATCH 18/27] Update check data with interpolation test values

---
 gsw_check_data.h      | 205 +++++++++++++++++++++++++++++++++++++++++-
 gsw_check_data.nc     | Bin 268800 -> 274652 bytes
 gsw_data_ncdump.txt   |  11 +++
 make_data_from_mat.py |   8 ++
 4 files changed, 223 insertions(+), 1 deletion(-)

diff --git a/gsw_check_data.h b/gsw_check_data.h
index 8a288c7..07b6b81 100644
--- a/gsw_check_data.h
+++ b/gsw_check_data.h
@@ -5,7 +5,7 @@
 **  version_date: 15th_May_2011
 **  version_number: 3.06.16
 **  mat_zip_version: 3_06_16
-**  This file created on: 2024-02-12 11:26:21.045021
+**  This file created on: 2024-09-12 12:38:26.783150
 */
 
 /*
@@ -24,6 +24,8 @@
 #define	cast_ice_n	3
 #define	cast_mpres_m	44
 #define	cast_mpres_n	3
+#define	interp_m	51
+#define	interp_n	3
 
 static UNUSED double	ct[135] = {
 27.996436412058213, 27.993857176639086, 27.944017615967979, 27.948372399994319,
@@ -296,6 +298,13 @@ static UNUSED double	p[135] = {
 8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90
 };
 
+static UNUSED double	p_i[51] = {
+0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400,
+1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700,
+2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000,
+4100, 4200, 4300, 4400, 4500, 4600, 4700, 4800, 4900, 5000
+};
+
 static UNUSED double	delta_p[135] = {
 0, 10, 10, 10, 10, 10, 26, 25, 25, 25, 25, 26, 50, 51, 50, 51, 101, 101, 101,
 101, 101, 101, 101, 102, 101, 102, 101, 254, 254, 254, 255, 255, 256, 255, 256,
@@ -7933,6 +7942,200 @@ static UNUSED double	rsubrho[132] = {
 
 #define	rsubrho_ca	1.690921180852456e-08
 
+static UNUSED double	sai_sactinterp[153] = {
+34.468236430490606, 34.981684011594275, 34.845532135029757, 34.594880552563986,
+34.646910055942165, 34.658726028373025, 34.667610684451823, 34.681241443337626,
+34.693152815925814, 34.704242814066355, 34.715310018759446, 34.727191420710305,
+34.737970181514079, 34.749908756622247, 34.75865652413632, 34.767977853907084,
+34.776667058040957, 34.784334257348284, 34.79178015056219, 34.799429517990411,
+34.806552614204548, 34.812691689382405, 34.818169026804789, 34.822997593767191,
+34.827267561289105, 34.83107566345695, 34.834475961281619, 34.837525533487018,
+34.840341326091348, 34.843005908156698, 34.845452404272699, 34.847569222979239,
+34.849391445147802, 34.851249802860892, 34.853492505612103, 34.855706803475442,
+34.857247772994597, 34.858497276933569, 34.859092057229148, 34.858381605762077,
+34.856823030966915, 34.855679183654935, 34.854739459543588, 34.854308742400669,
+34.855299304684273, 34.857513577115284, 34.85900682590129, 34.859483671997552,
+34.859866852752297, 34.860716166459341, 34.862075375111715, 34.556977983037378,
+35.143552028278961, 34.743955515099316, 34.823482547808247, 34.804942319288521,
+34.760842876649718, 34.73197711758587, 34.714860568930348, 34.713553992644862,
+34.718238136890584, 34.723558346430032, 34.732932331364822, 34.741835167338756,
+34.751656516251295, 34.760987253055077, 34.770634831831927, 34.779831462444207,
+34.787896487807899, 34.795973697286236, 34.804846003583471, 34.812671598162424,
+34.817995350850836, 34.822145849128248, 34.826452359442619, 34.831415297815248,
+34.836065183487634, 34.839513990464418, 34.842264529794129, 34.845044087445629,
+34.848173560330608, 34.851227079351254, 34.853788960522365, 34.856002592079761,
+34.858008149143096, 34.859816041891456, 34.861511533268335, 34.863314423422317,
+34.865265218099864, 34.86708175309419, 34.868590932972623, 34.869910820929661,
+34.871133587707021, 34.872546396459661, 34.873351690456062, 34.872775691608311,
+34.871547239164059, 34.870721688384769, 34.870383374148055, 34.870265396867914,
+34.870338679141682, 34.870546744141642, 6.6699043409245728, 10.346695553446363,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	sai_sactinterp_ca	1.3001510978938313e-10
+
+static UNUSED double	cti_sactinterp[153] = {
+27.996436412058227, 25.533144569030263, 16.019861520340768, 10.341891502382861,
+8.4456966374532545, 7.2721994214211767, 6.5074004434211012, 5.8840829086292556,
+5.3600631828344829, 4.8832999981632534, 4.4343297157729946, 4.0194186111174233,
+3.6464959786663744, 3.3302175416713258, 3.0360389453491519, 2.7908185312784872,
+2.5904631148025259, 2.423166010494068, 2.2715684781542738, 2.1281996714820086,
+2.0038859800808941, 1.9057066559935123, 1.8235890129026773, 1.7482666532235847,
+1.6762447258413102, 1.6115589460702915, 1.5572356514838051, 1.5101427761616257,
+1.4678690688374059, 1.4291052039518444, 1.3940324397786752, 1.3631384058258678,
+1.3360645039300703, 1.3092326083968704, 1.2796554533485569, 1.2504454968981085,
+1.2262396689330433, 1.2054680503482365, 1.187164480536496, 1.1684829755830797,
+1.1477518433658003, 1.1286972420975527, 1.1093760918020004, 1.0943939612593574,
+1.0873483431052355, 1.0817477224915897, 1.0745478755635225, 1.0634152009122975,
+1.0515515736003886, 1.0426620383447947, 1.0351918856149576, 27.322890736362599,
+22.011254836061603, 11.560227514375722, 9.8757443413275237, 9.0022892545844027,
+8.05506023421108, 7.1412969334874559, 6.2985045169753038, 5.6689611002002005,
+5.1619651842686372, 4.6499362388415042, 4.3058993935480157, 3.9759763533444907,
+3.6420020781766036, 3.3196655649536524, 2.9994634344611657, 2.7820843247989426,
+2.6147405296234192, 2.4435788959532885, 2.2471193199195505, 2.0822572605932823,
+1.9899103600091417, 1.9259139687763396, 1.8589909978970756, 1.7804728271451022,
+1.705937917028272, 1.6475130925660368, 1.5987301233149647, 1.5517208830746407,
+1.5039888179174263, 1.4579677965960154, 1.4151089811281323, 1.3749691776165345,
+1.3359607264588433, 1.2978547591325253, 1.2603848521798859, 1.2218814621548661,
+1.1820088935554813, 1.1431502947536156, 1.1066031829996259, 1.0714954834787538,
+1.0377190422510292, 1.0030329830608418, 0.97360861725787962,
+0.9534590677768886, 0.93592175165578984, 0.92023524972450943,
+0.90565719701788883, 0.8921688460133359, 0.8795567615391916,
+0.86702210345731012, 10.502767735302228, 4.5888811896604285,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	cti_sactinterp_ca	6.2611249518340628e-10
+
+static UNUSED double	traceri_tracerctinterp[153] = {
+34.468236430490606, 34.981684011594275, 34.845532135029757, 34.594880552563986,
+34.646910055942165, 34.658726028373025, 34.667610684451823, 34.681241443337626,
+34.693152815925814, 34.704242814066355, 34.715310018759446, 34.727191420710305,
+34.737970181514079, 34.749908756622247, 34.75865652413632, 34.767977853907084,
+34.776667058040957, 34.784334257348284, 34.79178015056219, 34.799429517990411,
+34.806552614204548, 34.812691689382405, 34.818169026804789, 34.822997593767191,
+34.827267561289105, 34.83107566345695, 34.834475961281619, 34.837525533487018,
+34.840341326091348, 34.843005908156698, 34.845452404272699, 34.847569222979239,
+34.849391445147802, 34.851249802860892, 34.853492505612103, 34.855706803475442,
+34.857247772994597, 34.858497276933569, 34.859092057229148, 34.858381605762077,
+34.856823030966915, 34.855679183654935, 34.854739459543588, 34.854308742400669,
+34.855299304684273, 34.857513577115284, 34.85900682590129, 34.859483671997552,
+34.859866852752297, 34.860716166459341, 34.862075375111715, 34.556977983037378,
+35.143552028278961, 34.743955515099316, 34.823482547808247, 34.804942319288521,
+34.760842876649718, 34.73197711758587, 34.714860568930348, 34.713553992644862,
+34.718238136890584, 34.723558346430032, 34.732932331364822, 34.741835167338756,
+34.751656516251295, 34.760987253055077, 34.770634831831927, 34.779831462444207,
+34.787896487807899, 34.795973697286236, 34.804846003583471, 34.812671598162424,
+34.817995350850836, 34.822145849128248, 34.826452359442619, 34.831415297815248,
+34.836065183487634, 34.839513990464418, 34.842264529794129, 34.845044087445629,
+34.848173560330608, 34.851227079351254, 34.853788960522365, 34.856002592079761,
+34.858008149143096, 34.859816041891456, 34.861511533268335, 34.863314423422317,
+34.865265218099864, 34.86708175309419, 34.868590932972623, 34.869910820929661,
+34.871133587707021, 34.872546396459661, 34.873351690456062, 34.872775691608311,
+34.871547239164059, 34.870721688384769, 34.870383374148055, 34.870265396867914,
+34.870338679141682, 34.870546744141642, 6.6699043409245728, 10.346695553446363,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	traceri_tracerctinterp_ca	1.3001510978938313e-10
+
+static UNUSED double	cti_tracerctinterp[153] = {
+27.996436412058227, 25.533144569030263, 16.019861520340768, 10.341891502382861,
+8.4456966374532545, 7.2721994214211767, 6.5074004434211012, 5.8840829086292556,
+5.3600631828344829, 4.8832999981632534, 4.4343297157729946, 4.0194186111174233,
+3.6464959786663744, 3.3302175416713258, 3.0360389453491519, 2.7908185312784872,
+2.5904631148025259, 2.423166010494068, 2.2715684781542738, 2.1281996714820086,
+2.0038859800808941, 1.9057066559935123, 1.8235890129026773, 1.7482666532235847,
+1.6762447258413102, 1.6115589460702915, 1.5572356514838051, 1.5101427761616257,
+1.4678690688374059, 1.4291052039518444, 1.3940324397786752, 1.3631384058258678,
+1.3360645039300703, 1.3092326083968704, 1.2796554533485569, 1.2504454968981085,
+1.2262396689330433, 1.2054680503482365, 1.187164480536496, 1.1684829755830797,
+1.1477518433658003, 1.1286972420975527, 1.1093760918020004, 1.0943939612593574,
+1.0873483431052355, 1.0817477224915897, 1.0745478755635225, 1.0634152009122975,
+1.0515515736003886, 1.0426620383447947, 1.0351918856149576, 27.322890736362599,
+22.011254836061603, 11.560227514375722, 9.8757443413275237, 9.0022892545844027,
+8.05506023421108, 7.1412969334874559, 6.2985045169753038, 5.6689611002002005,
+5.1619651842686372, 4.6499362388415042, 4.3058993935480157, 3.9759763533444907,
+3.6420020781766036, 3.3196655649536524, 2.9994634344611657, 2.7820843247989426,
+2.6147405296234192, 2.4435788959532885, 2.2471193199195505, 2.0822572605932823,
+1.9899103600091417, 1.9259139687763396, 1.8589909978970756, 1.7804728271451022,
+1.705937917028272, 1.6475130925660368, 1.5987301233149647, 1.5517208830746407,
+1.5039888179174263, 1.4579677965960154, 1.4151089811281323, 1.3749691776165345,
+1.3359607264588433, 1.2978547591325253, 1.2603848521798859, 1.2218814621548661,
+1.1820088935554813, 1.1431502947536156, 1.1066031829996259, 1.0714954834787538,
+1.0377190422510292, 1.0030329830608418, 0.97360861725787962,
+0.9534590677768886, 0.93592175165578984, 0.92023524972450943,
+0.90565719701788883, 0.8921688460133359, 0.8795567615391916,
+0.86702210345731012, 10.502767735302228, 4.5888811896604285,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	cti_tracerctinterp_ca	6.2611249518340628e-10
+
 static UNUSED double	p_mid_tursr[132] = {
 5, 15, 25, 35, 45, 63, 88.5, 113.5, 138.5, 163.5, 189, 227, 277.5, 328, 378.5,
 454.5, 555.5, 656.5, 757.5, 858.5, 959.5, 1060.5, 1162, 1263.5, 1365, 1466.5,
diff --git a/gsw_check_data.nc b/gsw_check_data.nc
index 75b54eb47fc5992d7a669402296c3f88605e51f8..e019372aaad1bd94d7e9f8537e21c71b9a996e6a 100644
GIT binary patch
delta 8478
zcmeHLZCq64wwE=-yz%}j?+&mTU^8xEzzGrAdlV5A5djsH4xmGVj!HP1c1Tds)I5qT
z&dG|Baw@CQFsanY=;SG1$SBE^MTt43*_fzkl&E*@8Bh0~d+&#PK6O7Fen0r{|61!=
z&-1MHJZtSuLx*?Yd2jDYnOQC)oIff8xBcFsE@OSomCMa#)l1ElwKZj?%B2g<RSuJ6
z`uk71r7M=qudH?GBpBzS@gFxjNvBTIr6s4P>Ga$KsmV!W6LG%Fvo5{Oz`=udOeC1k
zU^&4<rPxHUs0F(SmUiO+!E!08C7kAwbaE=GDXSKJiwjhW5BGm$B-%M*#2?L=FVVwe
z?zo9yZ6<a}bg-cg2PGVA?m$M$DK>Z05apyI+#-#a+r;NW1J0Gk!Jjg*L&`#v1^WqF
zr0DL-!IMTzcV*>IA0<P@BWQ4qaeI!TXIF~Ff=aBAM#A%r*yYNC#fEl*J((CJ<K(Yd
z$VgU&k+K+gU5aK|6ue=?oiYyIFyS$o2HvvZ0I|2MsFrKsunn{292~aeDme@9F=&x<
ziubhS!-H1$u`C4J+*I&^v00&W`@lqE@51Hgidxb{!bcV?S8(u=6`K?qIBCO!3LTuZ
zqn(6342CHQw3w@8q1TA3i1nIq53w&S*sJ6eUsxTp<l_>RbCy0ECaQGMXGf!oh4TzH
z5j(HNqr|>5;($u0_|Cb6d_19I;0FulyK(S?6&u`GxMD*qu`71$bJM_e2Bq$t$92a}
zCKZ)ERJ&{mu8j%A334UeGGe(q2e(Yv<j%sc7VISUs}=3;8qd38W)zP)29(vHH$=c)
zI|?2#5*G=BRUYw#T5R#qfZT|k9y*YlaKM8DHw&sgbqY6eQ@%B3JmMJ#UN&s>)Pa{B
z4|;Oo&)_9u{#p$4;vmR~xn4R5GGPO;FblSOWk8q}2fWfiZ9}a$3u-$Sdut$u!3JV6
zTC{pcLA()ryfvzLF%a_kmF5a;^9X^_77X*z!DuVy`>>E?!#W?<Ey*r!Wo6mQiuqXK
z6AXF=yL~v&YjM~|<C$g@wW0@yL61iSq?s_^S0m35kC1dFy2Pj@E*VyA^No^E6vINY
z=i^8=k#g)C1ydOG_R}b)h%uI}EW`8i06ga@hiOJE^<!1j#3ia$meD!778Qu!cq(89
zZPSm18CLA`)2IqXpYoMubR~KcA|)<`bR+#WP|RSlKdUMh-6|c!`~4%J)QDDpR$eMP
zRj$OVAwl@HrxMC7%>f*gTQM<!1)~j112l4@oor;Z{;+VUW^hk{4yv`-M+Qrbs0d_L
zOPo%`mAQmDSN5<4O9NTu!&VwfA*&*Vzhg5i=vW1&$ySIygRng?z!Azv>3EOSC?3^{
zGu2jM(j4d3*BG&CB<rz8oIq7uwXCd`Vn4<!6xPup8_B^sI%Fd?o*TuTJ7X^zh_T;D
zhb%|~PtX|)VwF#bj+NxV5jt_B+v&8CV!P<ssqhVSL;x-59n1n+F+Dg6p0r_kFjWsb
zHU(>-mBE9-oU&E)<}P|W?-%6QRhtpLLpW$Np%9`0iv_DfSPzSp%-qCP7aGKDIs>rB
zhMgf)UhH^@O!rdGLpj*1MIn@h{YEr}vdaD9q7<V4e?;v~+HR;0-lX$P)`x7U4rAev
z9dpApZbulhmW$I=VMRbB9MNK1m<Il8#A9Ks`(GU!Twc9!Nkt-V^ACd$>1c<u@Szp6
z!!>ZihBd@a*s(QS<MD~;O6D+59)X{DjDk<JI6!8njTjWcx}W|}XDJS6f4871f`h+X
zu{DB)uWWcMf>nKG7ni@dqNWDReL|q0!Jttb^lLGDlm=`@tRZGIVcRGcE?Kaf*d;3t
z6SLbePR)Ux?x~uEYYaB3b#P6K2Z`M@qK(*169z?c@QVerBU#xmR+`BT+oI^S-?3q9
zqz>-bu_rQK0>I#KWD=nk6QekgQHhGufy{(0QJh?5p>a|=lrfQ@q9YU)1)etaj@HOM
z#kI?++{cH5AA_aQtjbSIB30Et?%i{Y0!A9KJ(`7)ChUvWs6s_=@+#-Oi-J(PJd6fL
z*|3;l6{8%%peuPTDgduC3Y8dQSE`^?(<+!q1tWS_kbZ4U0PbSs5NE_8hINZGk=%{!
z9M?S@v=+>bVL@xfRWTY!v|$^uL_2oJXq00eX%lPLIc1at#%fWm(ZDz(W@}jGIMI<(
zR9>SK{w$+}2^MVCa4^A&-5QP3VEc#uy<m+B4B~p-MRW0NiwnwPFgKQiEG^c>Y7|qA
zq<5vo|AAp(stJ2yqulZ=PMH`1+7skaFN}*)6o{S`WO|oO3+z}Er-3;P?u=s<bHrV{
zk?Oo(uw#;WMzqJVFptt3uYocP3h^A2S+ORbg$0z^cpWUD%#yH*X=dp!rIlsDM5U3{
zfr&Z~32Q9a#b!W_wK+jAGmEoFU_(L#&P`B**^YGy8L*nc&ICQI*5Y7-4jwaNoHkAN
zm>8l^)M|sUK^p>VDLz_y>9At2mR>q+s2H6F8|;`lS`S+otQ)O^Em~|J&B8V#_KoIX
zn+e@H7Md-X%ITq*R?X?)FE+Fg`wR6?j)P|ybRR>Y#q=?(e3uwRH#0^i#mjezlY7!(
zNQ#l~76a*H#%Gg~U^m?f8u2pqu0$QYY{&9M4y@F~6Enc7#h%18Xg8vJQU<h}FgGa;
zIxMuq2>dBYj(d|Lp~H$jq&!IHS~pR4P)uty<|l<=qE3pdbRqCAgDpBeysO1-oethN
zHYe-heG{f8>)>y6t%?23ik4(fbzF?3bYX49O7vzT;J6)!lcV4igX$DLoYG=(3J0GV
zaZ`#8J~Lrw3adIJ?!&`zBbriJIAg`Iu`GN^??z*F@FitpEC=Tpv=Td~#eQNJsDSEq
zstaN;{2T)#$3~!d*DhMntml*$t)%f3Yod30s6+9i4SV!Dj~~UI`PM8mm;G1e0D~zH
z=wU#MrUy6}G-4~UK@;{qpo1X`DpCok!l!aDY(rD39)|5`NlkK>xG>_<#o|?gh6my#
z5*IDHk7L2rh^gar;A%qCI1ZF_L`0w~L^=-}4iWRPV_X`m@)bR(=;r1(r)gBaTG5}#
zFflC(0;zJOWk4Xk<fNrRhy}Idvt%LS-eNIvJT-&yArL`te&aJB!j41Z(;%95lCFbj
zEmozop0P$Um1P+#P3DRbBd{ke0%A?rm9BvV3tkeziry2VJjZmCouZ~}Nwu+TsSXFy
zBVbG~mQT<qk_SlxMJw>9bTuRoWBY_CkJK=d$qSaRm@mfgYWgThWtuZI+3A@i_bfD2
zQx%$2<osqJ35};;@sr}ey}3(|jgX}0lHBK?$6#JoShm5R6rSYfUvOp+7v}TNV?<}g
zlWn93G}l%vsI09nbK?JXg_G?xHu!%ZD>~1g+Cc`s;;jFg<EfoAF5rKT6P?10MvTkg
zV8&J~&S2dN#gxl(@@GGkhgZ@gp|AtnGpIcF;umBzrw<i8Rl{LS;q_1|#Wp?&8~G3@
zPp8K;D9^+`UN1A&k&`g&=SSigAzEf^Bs$7LmkFbwsuP<87OIY7r=S57gEoOGr51w>
z94y<4*&?*yDgz6(gJ?0NLG3X18}zWU6oV$}U?sg!Pvl_rQLGc88?6&H@R$~RC(>tG
z3c63?prHlRC+VS~70r`$uznc#OyXpXQd)%<4ouOwKW?N)x5}mF+GRDEH;IAAO_-9&
z!j^6<7oiuoW@=>HQpk!b@Av@xI5QHqrQ<-R9(J^1+++g!CYr2+KTEM?avJ>E9eYT)
z(}apF9qcq?dKT5xUNjMVrVrb)SlE@0-NbffHc#Q;1q&ulNrM-xXquvzz2xqgs6RUh
zyQW0IOa9nT(pStR4af5Aa2%H%39r;)X?B{-+UJn=-<Mi#*q2QOm%a&e^w6G*`8hf`
zV8zCqG&s<KT_k+lABS=>;O#I>oth>)ET$$JyQW6tGgC*w;YRG7IuYKrHBXxg@6vnQ
zv@Cc(7uC5D=uzx}?@bGZ_lvM^+C(_oLDDQZ+DR+|GjbKUEjJ8~g`qV!5024G4Jl7H
z5@Ya_5_dFDkC2_*N^?;|%2ng2xuG~_oE%P3*O{IrJ2gm}SnQupEj%v-J}bhiJSx(q
z*qX<|=bhM1?DM0|Ggvswpmv51&T6roSl?FMIztbAEqIiK=LT_z*tua$%-1O|l#)+h
zBenW_GscK@yER`2-yg-E{4Ck`-H!RnXNKWGz5<1rA#j;mRlXYQW@==Y#YD$q`Aqt5
zni(Ry(n3^2Gz=95tn5lF(Rc?fD$u~dFxD01!8Iv%7i5VAysjVy#}&rQ22EtA#fE}N
z+*}w5Khwc1%!8kMaky|Q+)lylq71m5j+=_o;BG4(C3e?>(plpqQjnrBYXYG=Hi(>R
zd8f!-%{V-3oLs8tCC{U<VfHAjpA#XGD*CW)HtF5dv2%8gT<V_bSYg*}4Stp5kI$7T
zz>7Ya=S-JKz38iMPL4$C=Z^L{Wa{Tnrg1nhCkB5h2$Dzx%~)4V@<8!e#9>2m45if<
z&&^Y~^r3r69t7L4v}7uTXA+C0%*M*XbIEiJdY42~R-<KVE74eyHs{H)c5WD`JFsbP
z7BK$!#avd&gppn`L$@3&y!~-k3B5=bVca|p;!3f29=$lWVB<U;Xj}0h2@`t{cppub
z9#EH;1lO0)Ncs^R|BoT8?Q-(f4j3WS3r@bS%ZVHMow(TnBZRGLC*LkOaYv#Po4cIY
zlJCSn_dD^)8Ye#GfD!e=(@jo6tJ;arwmWf`;KVlPf2f4#6P^6|At%1j<;2~aoM_2+
z;)_A|vEPX=?Q!DEHBQ`<>cm$ZP(MO=<(O0OYLgRRD|Vt)eIGA5ac{d5_pNf`>w**8
zrT4MwK8p8%BDp{DKDOM)At$~e-XYSz(M5kf=-~h(a>7mP7q5PN@oJa17q5n__v#Jf
zuO@umUoIFPd-=1u51kbZrnJiEKR+TEa%0>+U~USAf{(9%)Ay-hD4H3y^kI`=m}Sl^
z8W}AZ=IDN0d?{Zr6dxLIIO-`F<}z&;mz)+1^R`bZ-aJPzJm~+-(Se@@!yi_h+OfVw
zFqGYaowJ`73>Be!5*03jVgB^sgsCsp3x);Nr|Sdq1cR~CyDjxS!BAP^k;m*43=b7{
zoNfI<Ff5!9{+6LpFjUF1U3zK=4}FmFpsQf0o-?ewNsg;8T+CkZf?!xw{`$?0IfOl-
zhZ;z}I7PYOIxTc@)AXSYzJwQ_3%Tk~8+_o{$FtWGJ{EK9yZSPEIQZ4Mlup8%U6))=
z2!<tmhS!>rgpYS$z0MFGp0@Es0pW!SY9F5P*5S3X!-VH*Per5>TFW<$PbS>>UB>4V
z2&+OeUqix@JDQim2_LQf%it4)AGio7WrR1bZerdN3`-^3pQ)d{k{-OnZyt#yOt`*5
zwTQ5|Zb4<SV3;z^<JLFTf}y$Why4$JNFf||b<BDj?Ra5-PVWVZjH0B)FpomIAlEY8
zOd;RDB;-&JrJ&(Yf#sCKd9ym-NmxgCHI19|x?m`McKfVZ!DQESSo%gTrDo7vzqrgx
zFqG@3e7*e!IjWNXb#go9#uC4NAca!%Y2(8u>nJs6Pki)k5vAtaM+)PWG_h;at}i!J
zDqMH1;%f;ze2Yg?ZmZM%_$-Q8_3^}W?~~o4#h*&wrv)ySjtNYDm2gGUhbJNlKPmXd
zmClBVtzY<L+huxay0L!s88HQrV=+*OvIj3VH4*MulN9Yvcv>%uIzW-}Tz;6pN0<d$
z44)IO4V+SOobY{k$72=Y@TQ;mcLc+-(Zgp;LkLUK&#^4w+C1(w`C0b*Q~MeZ5gs4@
zb!-&j<hR2$!6)cpTOZT(DB<!`KYnnLu<GE|)pUrK8NZsGM4?-@aH7jz3hA<?eidIw
z3n6=NFU;#dA~5^c`IkK`d%NwPWRpw%g!A7tI4;@vdx7T5vla{EHr2)NlszOQx(2)`
zoiQvhue&cykZ%;!HeJd3r<6i?;5!fB_~oh)BFEb${oTUI&C})kyif?J=~$THn=JTk
zsy|t4^$~nJyRMw>&<oyv->&-UomGN|Ru-^mpjS{?a)J^|dId$!spAhC>Z=7=;m%hd
zpY}&VdMNWy=Z*t{WWJ?u&8bxa(3fXaxmj?L{b}n9!v<l5{ASsr^LP2Xk9{~PWlR@;
zXY^ugdSWksd+2*}=7t{r*H6E@Q!>rO|8nfQVtc|Ce(3VDpHf=Q{80Q8tByRH$p2ht
z%=-Hy1^lhApIsa@=|}$7^ue1i)j#zVfAf<?4{Q(4<8PL4(9U}<fgim7tS#f+WBlOi
zNA59y%;yJPPDjTd>f~=cx3}&M7mmM?v-9iY-?j4B@3wLK6pj4#H$qd&@6P70KbWMf
z?QiC<dogFenNi7K>vD<l_qxnqTemwyaZJKrtNf^Z#`WX;wY)P2;&0F9uO*d))eo&*
z$zPMax2NpmIDX)C>=Cx}VSZrWkp*u*v6~-gig~cFhY`jX?D@IvtIq{JoZs1Vs#b`N
zIKwE_c|vqo#r7LBrwLImTF>3gZ9yHD_~o~2@`SLac4q&gWFa`|wFJc~S0Qk(G<wah
z<AVQ-UnXx1x+nNPwdvbe(yt0W<0jO%cR#mJ@H*Aqet2@7;1Tud@btckf@+)f)3ZJ!
z1o_?XUpGE-QIHwlTr&Dpx*)xLCqU-|f@JIRkxy?}F1Y*}y#Lkk7-7Vom229*Y~=5q
zv1VK-v+#F&o?uOVgZ!N%ONvt$U*d=BKJ^p)hWOjlHf(z}`y2k(Qhj^o7w7q3TIwHC
zTC?x+L#NMXeXifZ4|&gSiSvJz|GDstDqs2s{+1;tfA<zY{^q@N`$k_qLMaNGdn#DX
z4<41w4ZIw}53-x{9zT}P-`IK@R~<0(*MC?SH+w0iV(Oe*FZu4|uXVh*MWrt0uQ6ZL
zzWL&t{J_)o*?Fs{@dNIr{KN}3{%T(A#}lVytl)oI*ZJgR%K(4nhOf3vILu$E-I4Sm
zSI*mSGanahRPgq#FYc38`ttUeqSQm|J^sfdFGg2QIw{1sEc)AVeUIQfPx9Jp*`M9N
zR2i=89lcWceL7YBZ>{S0=~VT<wW{BzQ`P_0s(znNRsUP7`h7Z8{co-6_vuu5zqP7=
acdCDPs(*K?|G%9|Dt4;3t*Mft^Zx-Usm(P2

delta 4299
zcmYM1e^6BC)yL&5-o5(+xFUuH7FZUzELXf)Zg7Ju5(L3j5w8IgjQ9dZlXYIyRSD6I
zKQ>L8?AzGMm9Nv-w`NKv>6CUfFGgE!LoG?ucH21h&LnAVYu@0DZEaumC7EeU8{53+
zS@KTzpR=Fyoafwg&U2pUKKI&nC3#*^+E#UB&71qVS-8$BC0TehtFpjYQf~z6>!U^_
zYSioDhNWSn9{-v3a@M<OBwgNeldZiE4x^bu@<2d9_sR@x@6Euv48*Mrm;z`hjgUik
z?q$_HPXbFEwe&<IPCFuWFo7E^4v)j1ZP3%xFtc@f&WEHBVN3`S!!(Qk7>6=nqc6*7
zVqwJ*pRW_rI5}T)z7%4|@{@RZzM7T743EEj6MjjfS7L~<cx@cRlI|Rr*-VP#u%u?i
zVN1ch*VC+Vwe>%?Z_lp&ZWC`gDx9y|tn;?vbyqR=xSVuOMyFG!a|#YQRk{$vYfKkR
zyzkWL&D5x@)0-B;vPM5lqfa*Ihc+yxcVtZR*gFbt$ttB|C~z5+HW72_bUB3sLRc7g
zX=JDIK9g-D+pW<DGQw`1K2R{g^y?UoxefZYi7BQJQ@F+Sp@o1)qw8sOdJMX5Bgyo8
z8K*ru{a(Q}ra#3Xdv)ngrpN`?Jtdh`?xZl_)#;7}%d62RX-qPGV&fLmUuF1mbosA}
zn2{L9151iBb3TnBk)zS4CQjv4=gs>p#p7~fu3}FioLxqrS(wRD$&rRH*C2<D&Rm^j
z8Hbo;1rtnOUPC6YiGn<xloVFw8KhV^z*LaNI8%X*YfMElyh<|_Dd<$9R2oA{(Y&Q1
z4fp=Ujz>rK;CR;}995iDk-{}aqY4Z1f*@6=QMW*KSBu4SZrd?D2y<aEX)+EhFi2By
zdV%KDV>~#21oipVPF=W^M=+RQN=6Eye4Pvnz4<Cd(n#iO6tQtSUv)1RJ(E4s&pmGE
zl+khpck_evKn#I}s`CNi#5sbNf_pL{wWTn`>Nbnj(g?D`en-~oG%haGoU27(?h)Ku
z<HKRaMV*{yfkvGQ8Vgi1V;EpEO&lvw-94hlQs$Y;q9FBHm}Z;aG~`8^yH_}Jjr6zm
zu>X5-M3HH$jJ`!WZB=lL4F*M)j*<TOyiCM{CZ^fRjuf)*G2A;uXiCQQ@kJ#f=-i#@
z|9=`!nVYkdTlgN;vs>no`Ru7bj_V65Ji8TUd5#w>D58DbWQ9@MXJS>MA&sPX)B!6H
zq!9~ag~6PsglFfz_@4fKgWGVlFpy<qrckA$+-g3Zjw*=xRC<OR%%{;a++aR6_j%z(
z5_7s&ejd&9++S?>0{2&u=6WHMKO1n>3p43kGKP4-%0wpyvcd(K<v0<;WRXrMOx!9`
z^G*uSt~t-5dYqB#q$r11%vn+3%gi~I#!#_Nr)-=mR&!sGd0hT{j+v?FD+*?cRr*E@
zz7oy*4U-K#!XUFAbDjAuZi5n?zGdN5iRL*YsxV)e%vJG>NFd`g4pGH7xEg+g-r$b+
z>-60i4*50uu8DEKO5abxX8OK`yM8tIlJF$$-}(5i!SLs`xn#pERq03E=cStWM+zIb
zG6wayljo<O#&Ee*qo10XWvh2n2$X5`u7%DrHTTNDyfu98?TU@DGL_z!kuKBeeFd{+
zn&)av#Q*PuiNv6&8C^ADmIw2$rNpxy9o)8U$I#B*kBgzaYZmLg|I&%{rqj4wuG6%Q
zyX6}FPDUu8)9)1Y1~mDP!nKDhrw2mvA5Df_hJOoG(d`tb0|wo;P*9=K$7!@w81%7?
z;R-$XFESg-J9ZDp2Zs^|GmG*U1$Ko>e~)37?f!0}E?Dcx`ezEg!FrY!js<lhKJS7C
z5jP=|Q%0atbvi}KoN;6=^gGBC!w~B{A`;e#ZT(-BPRdK+Vx{KH6WO`rnJs*#(w&ut
zw@RahHX5r`^2r!r^6?(5Qa%1m&Dwgzv2nbig8aM*tAg%;sM&m0Ctx)MDZuZI>R@hF
znzeGqSlBDo<-AXl)hgA>7_ZhmwIW2Wk{GO4OKC|AcdIp8Vj`rfRG)&WYE*CGkgB?t
z3SW{q{`TBgNlR_mY!a0*tLhY05UNqBIR>*vr)CpLrq&e3Yjjtu$jvL}9p{bk^%^-V
zJsQ%e!$w^wNb6*nOzZd+$h0AbF{TYBE{9Z4x5!CO42$UYR#s3qznp4S&%>FdSj`K{
zYWZPXB*E{M;{`>TqV~x+RjX5<g3GlU4e-j<8Z^Kwr|Gmkg+3v8-87Bje1d8Q#rZwO
z%3U)2i<@beg1*HOX>W|(mSJFVG2U2QCG9mCD>C@0=ELiiGCh%kuP#DQSTO4hI+(_@
zbx}HKV~UkWxG+K}2rn_{3o&e85~VMg7+Vsd=Tf-J%I7S2b&bB1Mx(CNmu$qDUX*b{
zkI;)8yKWF-qxVL*Lj33g&i1Wvg-Fo(k|>OK7+LkT&X+`a^M<jqzS{XUQ5nxL4zvDi
zqAIz=*sYglf*s?g3`d+}3abm)+^BWFCNflp@wBmoUNbQn)@hut0pVsEw@_Cfl}-zP
zB}f=0c)q?|IxSq5G5T>(&lg^$vofnH5H>2{GAih-f;yv_zQfIKw9<D>OfY?q?@dO8
zzGoqHpFuz1tIT~7`hkrTOmFjf!1T5<+Mv@<VhA^=?w^?K)A{(m!4ZrFi|OYn9BPQr
z&n=v1n_r}Hqrq_gLWH8kF}<`*RK)$Bh$AO1TH;oN=6+9Psbp}9cf+svX0_CyUzxbN
zRHxsjAV&=Pt%atDP9O29%=8hTPfRyt*g`0y_iJ=BhPwL=x@ltIew}WmVBH^~TNctx
zGil5+&DaPx8uVuw+l5eYvQekM#gJ}{(BI6_W%arAkCX^4bIA1%<9{1!=pPp1%QWI$
zeqxzHyd%@gbdtDXgy5YN)yN$~PgHliMJaN4fA*k0T8ul(WXerpBC1m^zq#2ipHGVA
z&6LkKqvcU5l#yEAL4|x|EpMd~z8E$&Q;CV5rYMzj_)P|tTX?5Qqsla9gs|aXp*lmN
zdbAtITPi4|U}!~<>cX&AXjIpT$rY*`7WZ%Jc>-%!1SmWJZ*!0uk_b0z)W936S*6Bt
z)V1i;IDx(vjhbwvS~O47Ri5TrV?Oj)7xu0wrIuM-W22Sx;9aRvdjO3qBh(&3Vx>WA
z23YtoxiUy=cwKp5{TO_$&9r_Ty{%E&IE@pn8g0CZsaBOX#c;FLpiL(HZ94T%pi_v8
zIMk+5pAX}0QK>J$4==&#Hs9!~V(L%AzsjI~KGIg{wCyI6Ld;-dl}e9sU#yDIW4%z?
zb=r9mE$tE7IfVo52JH&qG%I(7aIIaV1h-CyP6-Q*9h$s<hMyokn!ta@s%ZZ#tPWLr
zQdG(ZW%(k!)m9-rIgklLV9x4NdTI)Ft0Pk8_EL!O>H-{ET`V05&COBf<{Syn&0)OF
zbDrj7BE}yuC(#nq>DgH%Vn)sv*<_x=XQ24JKzVyPr3PTLSt^d~HM;adI^#01)(3M<
z8NFa5u_hwDC_MWx#@fmDGU>%e#uCON+-@(T(Ub76ZKcsM46KdPSLWdq(^q6nGkrA<
z|GJ3u)dbHdL||P3hSn8Jr)+_vorQRHUAc5htakBerw`NXifOD7ch@!3SPPmuqx9++
z4s>?VtK*pJY?Zz#`YnV*os~GUSw?wh34K$6xxPdCws=4giLPS&bG?j9>jU)d1DIal
zDxDGD#Yhlqy8?7(nzd!P{YdfXayOlA!9-UVosErd=%DlC3{@D~pwan>(G69Y+z_OT
z3Z^%7NEdzVr3M@Ns$iPM^yUE!nFhU?#F(kmrD@nqmu@0^qegG@X|&OBy>0RgWjh}%
zo{g3C;|W+BBlP2oNN+UgXFl9z`dI)C4;pkiiTHyN@dY-<bmb-nx&vtHmT~vNGP*JY
ze|Iar$7e}*2ff$J6hM*(?e1dw<we};UQ556Lb4|SZ;u;Nm!E#a20dN$n-J`t4)KAO
z?5W0zO|{}9&E~&f-c&);+`*f=XnJ<^p|y0Q7o86^(~SWfdnihOm_nNA4{68`FLOBl
zCxlfGuV5L*vx4tfm=b*F<Y+JB$FpeaUFmRqG7rgKR{wbbQ@v|g#^K%E?sUwGhFT0|
za}ADe4l;|S*E4)&bD6_&H;kFhY=5^AO^>jj&rTxw2=ixSxco>fIp#sx(nSuw%4}In
zQk<!VSICQxwv>>RU~LUAQwesvoX%{Q@Lhw6EtS~QS1Gx~%4-?yGz*1^BRi)N*ZMSf
zP74po8;t$p+if9#R7#HQ1qV>KRi_0>^ljB>(M6=T8nkE%ldLS7J$GyC@zQgLk9<Ar
jT<76<)-(R(^ZsPk1%Hz2Cog<(ucIdW!uf7T=)?a57GWw3

diff --git a/gsw_data_ncdump.txt b/gsw_data_ncdump.txt
index f840ab1..20397be 100644
--- a/gsw_data_ncdump.txt
+++ b/gsw_data_ncdump.txt
@@ -15,6 +15,8 @@ dimensions:
 	test_cast_midpressure_number = 3 ;
 	test_cast_midlocation_length = 45 ;
 	test_cast_midlocation_number = 2 ;
+    test_interp_length = 51 ;
+    test_interp_number = 3 ;
 variables:
 	double p_ref(nz) ;
 		p_ref:standard_name = "pressure_reference_atlas" ;
@@ -85,6 +87,7 @@ variables:
 	double SP_chck_cast(test_cast_number, test_cast_length) ;
 	double t_chck_cast(test_cast_number, test_cast_length) ;
 	double p_chck_cast(test_cast_number, test_cast_length) ;
+	double p_i(test_interp_length) ;
 	double lat_chck_cast(test_cast_number) ;
 	double long_chck_cast(test_cast_number) ;
 	double pr(value_test_cast) ;
@@ -320,6 +323,14 @@ variables:
 		Tu:computation_accuracy = 2.19097096021414e-08 ;
 	double Rsubrho(test_cast_number, test_cast_midpressure_length) ;
 		Rsubrho:computation_accuracy = 1.70965819279445e-08 ;
+	double SAi_SACTinterp(test_interp_number, test_interp_length) ;
+		SAi_SACTinterp:computation_accuracy = 1.300151097893831e-10 ;
+	double CTi_SACTinterp(test_interp_number, test_interp_length) ;
+		CTi_SACTinterp:computation_accuracy = 6.261124951834063e-10 ;
+	double traceri_tracerCTinterp(test_interp_number, test_interp_length) ;
+		traceri_tracerCTinterp:computation_accuracy = 1.300151097893831e-10 ;
+	double CTi_tracerCTinterp(test_interp_number, test_interp_length) ;
+		CTi_tracerCTinterp:computation_accuracy = 6.261124951834063e-10 ;
 	double p_mid_TuRsr(test_cast_number, test_cast_midpressure_length) ;
 		p_mid_TuRsr:computation_accuracy = 2.30002115131356e-08 ;
 	double n2(test_cast_number, test_cast_midpressure_length) ;
diff --git a/make_data_from_mat.py b/make_data_from_mat.py
index dede788..6dc8b7a 100644
--- a/make_data_from_mat.py
+++ b/make_data_from_mat.py
@@ -97,6 +97,7 @@ def get_arraydims():
     ["sp", "SP_chck_cast"],
     ["t", "t_chck_cast"],
     ["p", "p_chck_cast"],
+    ["p_i", "p_i"],
     ["delta_p", "delta_p_chck_cast"],
     ["p_shallow", "p_chck_cast_shallow"],
     ["p_deep", "p_chck_cast_deep"],
@@ -271,6 +272,10 @@ def get_arraydims():
     ['p_mid_n2', ""],
     ['Tu', ""],
     ['Rsubrho', ""],
+    ['SAi_SACTinterp', ""],
+    ['CTi_SACTinterp', ""],
+    ['traceri_tracerCTinterp', ""],
+    ['CTi_tracerCTinterp', ""],
     ['p_mid_TuRsr', ""],
     ['IPVfN2', ""],
     ['p_mid_IPVfN2', ""],
@@ -380,6 +385,7 @@ def get_arraydims():
 work_dims["cast_m"], work_dims["cast_n"] = cv['p_chck_cast'].shape
 work_dims["cast_ice_m"], work_dims["cast_ice_n"] = cv["p_Arctic"].shape
 work_dims["cast_mpres_m"], work_dims["cast_mpres_n"] = cv["p_mid_n2"].shape
+work_dims["interp_m"], work_dims["interp_n"] = cv['SAi_SACTinterp'].shape
 
 for key, val in work_dims.items():
     out.write("#define\t%s\t%d\n" % (key, val))
@@ -410,6 +416,8 @@ def get_arraydims():
     test_cast_midpressure_number = 3,
     test_cast_midlocation_length = 45,
     test_cast_midlocation_number = 2,
+    test_interp_length = 51,
+    test_interp_number = 3
 )
 arraydims = get_arraydims()
 

From 832030a037b885af1b43033a5b7588f54d1965be Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 13 Sep 2024 12:28:56 -0700
Subject: [PATCH 19/27] Make sure memory is freed after returning from an error

---
 gsw_oceanographic_toolbox.c | 22 +++++++++++++---------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index 45590c4..ce81871 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -8649,14 +8649,6 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
                *sa_i_limiting_obs_plus_interp, *ct_i_limiting_obs_plus_interp,
                *ctf_i_tointerp, *sa_i_tooutput, *ct_i_tooutput;
 
-        sa_obs = (double *) malloc(m*sizeof (double));
-        ct_obs = (double *) malloc(m*sizeof (double));
-        p_obs = (double *) malloc(m*sizeof (double));
-        p_idx = (int *) malloc(m*sizeof (int));
-        p_sort = (double *) malloc(m*sizeof (double));
-        sa_sort = (double *) malloc(m*sizeof (double));
-        ct_sort = (double *) malloc(m*sizeof (double));
-
         if (m < 4)
             return 1; // There must be at least 4 bottles'
 
@@ -8672,7 +8664,11 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
             sa_i[i] = NAN;
             ct_i[i] = NAN;
         }
-        
+
+        sa_obs = (double *) malloc(m*sizeof (double));
+        ct_obs = (double *) malloc(m*sizeof (double));
+        p_obs = (double *) malloc(m*sizeof (double));
+
         // Find NaNs in profile
         prof_len = 0;
         for (i=0; i<m; ++i) {
@@ -8695,6 +8691,9 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         }
 
         if (prof_len < 2) {
+            free(sa_obs);
+            free(ct_obs);
+            free(p_obs);
             return 1;
         }
 
@@ -8713,6 +8712,11 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
             }
         }
 
+        p_idx = (int *) malloc(prof_len*sizeof (int));
+        p_sort = (double *) malloc(prof_len*sizeof (double));
+        sa_sort = (double *) malloc(prof_len*sizeof (double));
+        ct_sort = (double *) malloc(prof_len*sizeof (double));
+
         if (not_monotonic > 0) {
             gsw_util_sort_real(p_obs, prof_len, p_idx);
             for (i=0; i<prof_len; ++i) {

From f2a25960c20d69c12d83f21d3586a0b074758eb1 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Fri, 13 Sep 2024 12:40:53 -0700
Subject: [PATCH 20/27] Allocate arrays after NaN check

---
 gsw_oceanographic_toolbox.c | 27 ++++++++++++++++-----------
 1 file changed, 16 insertions(+), 11 deletions(-)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index ce81871..c0d355a 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -8665,11 +8665,23 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
             ct_i[i] = NAN;
         }
 
-        sa_obs = (double *) malloc(m*sizeof (double));
-        ct_obs = (double *) malloc(m*sizeof (double));
-        p_obs = (double *) malloc(m*sizeof (double));
-
         // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = sa[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                ++prof_len;
+            }
+        }
+
+        if (prof_len < 2) {
+            return 1;
+        }
+
+        sa_obs = (double *) malloc(prof_len*sizeof (double));
+        ct_obs = (double *) malloc(prof_len*sizeof (double));
+        p_obs = (double *) malloc(prof_len*sizeof (double));
+
         prof_len = 0;
         for (i=0; i<m; ++i) {
             d = sa[i] + ct[i] + p[i];
@@ -8690,13 +8702,6 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
             }
         }
 
-        if (prof_len < 2) {
-            free(sa_obs);
-            free(ct_obs);
-            free(p_obs);
-            return 1;
-        }
-
         p_i_tmp = (double *) malloc(m_i*sizeof (double));
         for (i=0; i<m_i; ++i) {
             p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);

From 532c4bcb5a9d2ad1899afc7354921ed45cf983c1 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Mon, 23 Sep 2024 11:27:48 -0700
Subject: [PATCH 21/27] Apply the same memory allocation change from
 gsw_sa_ct_interp to gsw_tracer_ct_interp

---
 gsw_oceanographic_toolbox.c | 33 +++++++++++++++++++++------------
 1 file changed, 21 insertions(+), 12 deletions(-)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index c0d355a..3f51c75 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -11378,14 +11378,6 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
                *tracer_i_obs_plus_interp, *ct_i_obs_plus_interp,
                *tracer_i_tooutput, *ct_i_tooutput;
 
-        tracer_obs = (double *) malloc(m*sizeof (double));
-        ct_obs = (double *) malloc(m*sizeof (double));
-        p_obs = (double *) malloc(m*sizeof (double));
-        p_idx = (int *) malloc(m*sizeof (int));
-        p_sort = (double *) malloc(m*sizeof (double));
-        tracer_sort = (double *) malloc(m*sizeof (double));
-        ct_sort = (double *) malloc(m*sizeof (double));
-
         if (m < 4)
             return 1; // There must be at least 4 bottles'
 
@@ -11403,6 +11395,22 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         }
         
         // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = tracer[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                ++prof_len;
+            }
+        }
+
+        if (prof_len < 2) {
+            return 1;
+        }
+
+        tracer_obs = (double *) malloc(prof_len*sizeof (double));
+        ct_obs = (double *) malloc(prof_len*sizeof (double));
+        p_obs = (double *) malloc(prof_len*sizeof (double));
+        
         prof_len = 0;
         for (i=0; i<m; ++i) {
             d = tracer[i] + ct[i] + p[i];
@@ -11416,10 +11424,6 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
             }
         }
 
-        if (prof_len < 2) {
-            return 1;
-        }
-
         p_i_tmp = (double *) malloc(m_i*sizeof (double));
         for (i=0; i<m_i; ++i) {
             p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
@@ -11435,6 +11439,11 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
             }
         }
 
+        p_idx = (int *) malloc(prof_len*sizeof (int));
+        p_sort = (double *) malloc(prof_len*sizeof (double));
+        tracer_sort = (double *) malloc(prof_len*sizeof (double));
+        ct_sort = (double *) malloc(prof_len*sizeof (double));
+
         if (not_monotonic > 0) {
             gsw_util_sort_real(p_obs, prof_len, p_idx);
             for (i=0; i<prof_len; ++i) {

From 73200b4830073d0d172e7f61e7641f28bb698104 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Mon, 23 Sep 2024 13:31:31 -0700
Subject: [PATCH 22/27] Add tests for vertical interpolation functions

---
 gsw_check_functions.c | 35 ++++++++++++++++++++++++++++++++++-
 1 file changed, 34 insertions(+), 1 deletion(-)

diff --git a/gsw_check_functions.c b/gsw_check_functions.c
index a396a4e..06fc48c 100644
--- a/gsw_check_functions.c
+++ b/gsw_check_functions.c
@@ -130,7 +130,8 @@ main(int argc, char **argv)
         int     count = cast_m*cast_n, i, j, k, l, n;
         double  saturation_fraction, value[cast_m*cast_n], lat[cast_m*cast_n],
                 lon[cast_m*cast_n], val1[cast_m*cast_n], val2[cast_m*cast_n], val3[cast_m*cast_n],
-                val4[cast_m*cast_n], val5[cast_m*cast_n];
+                val4[cast_m*cast_n], val5[cast_m*cast_n], val6[interp_m*interp_n],
+                val7[interp_m*interp_n];
 
         if (argc==2 && !strcmp(argv[1],"-debug"))
             debug       = 1;
@@ -583,6 +584,38 @@ main(int argc, char **argv)
         test_func(o2sol, (sa[i],ct[i],p[i],lon[i],lat[i]), value, o2sol);
         test_func(o2sol_sp_pt, (sp[i],pt[i]), value, o2sol_sp_pt);
 
+        section_title("Vertical Interpolation");
+
+        for (j = 0; j<cast_n; j++) {
+            k = j*cast_m;
+            l = j*interp_m;
+            for (n=0; n<cast_m; n++)
+                if (isnan(sa[k+n]) || fabs(sa[k+n]) >= GSW_ERROR_LIMIT)
+                    break;
+            if (gsw_sa_ct_interp(&sa[k],&ct[k],&p[k],n,
+                p_i,interp_m,&val6[l],&val7[l]) == 1)
+                printf("gsw_sa_ct_interp returned error.\n");
+        }
+        check_accuracy("gsw_sa_ct_interp",sai_sactinterp_ca,
+                "sai_sactinterp",interp_n*interp_m, val6, sai_sactinterp);
+        check_accuracy("gsw_sa_ct_interp",cti_sactinterp_ca,
+                "cti_sactinterp",interp_n*interp_m, val7, cti_sactinterp);
+
+        for (j = 0; j<cast_n; j++) {
+            k = j*cast_m;
+            l = j*interp_m;
+            for (n=0; n<cast_m; n++)
+                if (isnan(sa[k+n]) || fabs(sa[k+n]) >= GSW_ERROR_LIMIT)
+                    break;
+            if (gsw_tracer_ct_interp(&sa[k],&ct[k],&p[k],n,
+                p_i,interp_m,9.,&val6[l],&val7[l]) == 1)
+                printf("gsw_sa_ct_interp returned error.\n");
+        }
+        check_accuracy("gsw_tracer_ct_interp",traceri_tracerctinterp_ca,
+                "traceri_tracerctinterp",interp_n*interp_m, val6, traceri_tracerctinterp);
+        check_accuracy("gsw_tracer_ct_interp",cti_tracerctinterp_ca,
+                "cti_tracerctinterp",interp_n*interp_m, val7, cti_tracerctinterp);
+
         test_infunnel();
 
         if (gsw_error_flag)

From d5cad06ecec1ff1d43098089f596b7c8cecec04d Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Mon, 23 Sep 2024 13:54:16 -0700
Subject: [PATCH 23/27] clean up

---
 Makefile            | 2 +-
 gsw_data_ncdump.txt | 4 ++--
 gswteos-10.h        | 2 +-
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Makefile b/Makefile
index 75cb1cb..1fe0b19 100644
--- a/Makefile
+++ b/Makefile
@@ -5,7 +5,7 @@ PROGRAM_SOURCES = gsw_check_functions.c\
                   gsw_oceanographic_toolbox.c\
                   gsw_saar.c
 LIBRARY_SRCS    = gsw_oceanographic_toolbox.c \
-                  gsw_saar.c 
+                  gsw_saar.c
 
 # This includes TOOLS.gcc if make (unix) is used
 # The #\ logic causes the include to be ignored by nmake (windows)
diff --git a/gsw_data_ncdump.txt b/gsw_data_ncdump.txt
index 20397be..4b6b4b2 100644
--- a/gsw_data_ncdump.txt
+++ b/gsw_data_ncdump.txt
@@ -15,8 +15,8 @@ dimensions:
 	test_cast_midpressure_number = 3 ;
 	test_cast_midlocation_length = 45 ;
 	test_cast_midlocation_number = 2 ;
-    test_interp_length = 51 ;
-    test_interp_number = 3 ;
+	test_interp_length = 51 ;
+	test_interp_number = 3 ;
 variables:
 	double p_ref(nz) ;
 		p_ref:standard_name = "pressure_reference_atlas" ;
diff --git a/gswteos-10.h b/gswteos-10.h
index 74be239..e273cad 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -299,7 +299,7 @@ DECLSPEC extern double gsw_t_freezing_poly(double sa, double p,
 DECLSPEC extern double gsw_t_from_ct(double sa, double ct, double p);
 DECLSPEC extern double gsw_t_from_pt0_ice(double pt0_ice, double p);
 DECLSPEC extern double gsw_thermobaric(double sa, double ct, double p);
-DECLSPEC extern int    gsw_tracer_ct_interp(double *sa, double *ct, double *p,
+DECLSPEC extern int    gsw_tracer_ct_interp(double *tracer, double *ct, double *p,
                 int m, double *p_i, int m_i, double factor, double *sa_i, double *ct_i);
 DECLSPEC extern void   gsw_turner_rsubrho(double *sa, double *ct, double *p, int nz,
                 double *tu, double *rsubrho, double *p_mid);

From e18e4c3ac70470ed8ee70976a800e8689ba78a9e Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Mon, 23 Sep 2024 13:56:42 -0700
Subject: [PATCH 24/27] Use correct function name in error message in test

---
 gsw_check_functions.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gsw_check_functions.c b/gsw_check_functions.c
index 06fc48c..b8bab4b 100644
--- a/gsw_check_functions.c
+++ b/gsw_check_functions.c
@@ -609,7 +609,7 @@ main(int argc, char **argv)
                     break;
             if (gsw_tracer_ct_interp(&sa[k],&ct[k],&p[k],n,
                 p_i,interp_m,9.,&val6[l],&val7[l]) == 1)
-                printf("gsw_sa_ct_interp returned error.\n");
+                printf("gsw_tracer_ct_interp returned error.\n");
         }
         check_accuracy("gsw_tracer_ct_interp",traceri_tracerctinterp_ca,
                 "traceri_tracerctinterp",interp_n*interp_m, val6, traceri_tracerctinterp);

From 3999a39c068fe6b05d668c157eda71a368a90938 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Mon, 23 Sep 2024 14:21:19 -0700
Subject: [PATCH 25/27] Fix pointer type issue with malloc in
 gsw_util_intersect

---
 gsw_oceanographic_toolbox.c | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index 3f51c75..e8bb04c 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -12260,10 +12260,10 @@ gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy)
             return 0;
         }
 
-        sort_ix = malloc(nx * sizeof(int));
-        sort_iy = malloc(ny * sizeof(int));
-        unique_ix = malloc(nx * sizeof(int));
-        unique_iy = malloc(ny * sizeof(int));
+        sort_ix = (int *) malloc(nx * sizeof(int));
+        sort_iy = (int *) malloc(ny * sizeof(int));
+        unique_ix = (int *) malloc(nx * sizeof(int));
+        unique_iy = (int *) malloc(ny * sizeof(int));
 
         gsw_util_sort_real(x, nx, sort_ix);
         gsw_util_sort_real(y, ny, sort_iy);

From 5e6015b24cadc51b41905cc5ba672a7972395263 Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Mon, 23 Sep 2024 14:37:19 -0700
Subject: [PATCH 26/27] Trim trailing whitespace

---
 gsw_oceanographic_toolbox.c | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index e8bb04c..806d998 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -8633,7 +8633,7 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
                     1.950903220161283e-1
                 };
 
-        int i, j, k, prof_len, 
+        int i, j, k, prof_len,
             not_monotonic, unique_count, new_len, p_all_len,
             i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
             i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
@@ -8692,7 +8692,7 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
                 p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
 
                 ct_f = gsw_ct_freezing_poly(sa_obs[prof_len],
-                                            p_obs[prof_len], 
+                                            p_obs[prof_len],
                                             0.);
                 if (ct_obs[prof_len] < (ct_f - 0.1)) {
                     ct_obs[prof_len] = ct_f;
@@ -11329,7 +11329,7 @@ subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i,factor)
 !           ( i.e. absolute pressure - 10.1325 dbar )
 !  p_i  =  specific query points at which the interpolated tracer_i
 !          and CT_i are required                                   [ dbar ]
-!  factor  =  ratio between the ranges of Conservative Temperature 
+!  factor  =  ratio between the ranges of Conservative Temperature
 !             and tracer in the world ocean.                          [ ? ]
 
 !
@@ -11363,7 +11363,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
                     1.950903220161283e-1
                 };
 
-        int i, j, k, prof_len, 
+        int i, j, k, prof_len,
             not_monotonic, unique_count, new_len, p_all_len,
             i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
             i_out_len, i_2_len;
@@ -11393,7 +11393,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
             tracer_i[i] = NAN;
             ct_i[i] = NAN;
         }
-        
+
         // Find NaNs in profile
         prof_len = 0;
         for (i=0; i<m; ++i) {
@@ -11410,7 +11410,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         tracer_obs = (double *) malloc(prof_len*sizeof (double));
         ct_obs = (double *) malloc(prof_len*sizeof (double));
         p_obs = (double *) malloc(prof_len*sizeof (double));
-        
+
         prof_len = 0;
         for (i=0; i<m; ++i) {
             d = tracer[i] + ct[i] + p[i];
@@ -11554,7 +11554,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
         }
         gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
                                 p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
-        
+
         scaled_tracer_obs = (double *) malloc(prof_len*sizeof (double));
         v_tmp = (double *) malloc(prof_len*sizeof (double));
         q_tmp = (double *) malloc(prof_len*sizeof (double));

From 78655753e38aa8babda13b047e7457c489259ecd Mon Sep 17 00:00:00 2001
From: mauzey1 <mauzey1@llnl.gov>
Date: Mon, 23 Sep 2024 16:46:51 -0700
Subject: [PATCH 27/27] Remove more trailing whitespace

---
 gsw_oceanographic_toolbox.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index 806d998..1bda1f3 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -8832,7 +8832,7 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
         }
         gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
                                 p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
-        
+
         scaled_sa_obs = (double *) malloc(prof_len*sizeof (double));
         v_tmp = (double *) malloc(prof_len*sizeof (double));
         q_tmp = (double *) malloc(prof_len*sizeof (double));