From 171b379edaafac1524539ab49a1148d0bc9fbe18 Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Fri, 6 Sep 2024 16:44:38 -0700 Subject: [PATCH 01/17] Update GSW-C with SA CT interpolation functions --- src/c_gsw/gsw_oceanographic_toolbox.c | 814 ++++++++++++++++++++++++++ src/c_gsw/gswteos-10.h | 5 + 2 files changed, 819 insertions(+) diff --git a/src/c_gsw/gsw_oceanographic_toolbox.c b/src/c_gsw/gsw_oceanographic_toolbox.c index 371bce2..829306e 100644 --- a/src/c_gsw/gsw_oceanographic_toolbox.c +++ b/src/c_gsw/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 0) { + gsw_util_sort_real(p_obs, prof_len, p_idx); + for (i=0; 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 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 0) { + gsw_util_sort_real(p_obs, prof_len, p_idx); + for (i=0; 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 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 **!========================================================================== diff --git a/src/c_gsw/gswteos-10.h b/src/c_gsw/gswteos-10.h index 0d43caa..fbfa16f 100644 --- a/src/c_gsw/gswteos-10.h +++ b/src/c_gsw/gswteos-10.h @@ -228,6 +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, 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, @@ -297,11 +299,14 @@ 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); DECLSPEC extern double *gsw_util_interp1q_int(int nx, double *x, int *iy, int nxi, double *x_i, double *y_i); +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 02aaa27764a1319579ad00bd50b3a7d79fc51d8a Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Fri, 6 Sep 2024 17:46:42 -0700 Subject: [PATCH 02/17] Update SA-CT interp functions --- src/c_gsw/gsw_oceanographic_toolbox.c | 20 ++++++++++---------- src/c_gsw/gswteos-10.h | 4 ++-- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/c_gsw/gsw_oceanographic_toolbox.c b/src/c_gsw/gsw_oceanographic_toolbox.c index 829306e..45590c4 100644 --- a/src/c_gsw/gsw_oceanographic_toolbox.c +++ b/src/c_gsw/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 Date: Fri, 6 Sep 2024 17:55:27 -0700 Subject: [PATCH 03/17] Add sa_ct_interp function --- src/method_bodies.c | 100 +++++++++++++++++++++++++++++++++++++++ src/method_def_entries.c | 2 + 2 files changed, 102 insertions(+) diff --git a/src/method_bodies.c b/src/method_bodies.c index c9472c4..3fb953e 100644 --- a/src/method_bodies.c +++ b/src/method_bodies.c @@ -219,3 +219,103 @@ util_pchip_interp(PyObject *NPY_UNUSED(self), PyObject *args) } return (PyObject *)yia; } + + +static PyObject * +sa_ct_interp(PyObject *NPY_UNUSED(self), PyObject *args) +{ + PyObject *sa_o, *ct_o, *p_o, *pi_o, *res; + PyArrayObject *sa_a, *ct_a, *p_a, *pi_a, *sai_a, *cti_a; + int np, npi; + int ret = 1; /* error (1) until set to 0 by the C function */ + + if (!PyArg_ParseTuple(args, "OOOO", &sa_o, &ct_o, &p_o, &pi_o)) + return NULL; + + sa_a = (PyArrayObject *)PyArray_ContiguousFromAny(sa_o, NPY_DOUBLE, 1, 1); + if (sa_a == NULL) + return NULL; + + ct_a = (PyArrayObject *)PyArray_ContiguousFromAny(ct_o, NPY_DOUBLE, 1, 1); + if (ct_a == NULL) + { + Py_DECREF(sa_a); + return NULL; + } + p_a = (PyArrayObject *)PyArray_ContiguousFromAny(p_o, NPY_DOUBLE, 1, 1); + if (p_a == NULL) + { + Py_DECREF(sa_a); + Py_DECREF(ct_a); + return NULL; + } + pi_a = (PyArrayObject *)PyArray_ContiguousFromAny(pi_o, NPY_DOUBLE, 1, 1); + if (pi_a == NULL) + { + Py_DECREF(sa_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + return NULL; + } + + np = PyArray_DIM(sa_a, 0); + if (PyArray_DIM(ct_a, 0) != np || PyArray_DIM(p_a, 0) != np) + { + PyErr_SetString(PyExc_ValueError, + "Arguments SA, CT, and p must have the same dimensions."); + Py_DECREF(sa_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + return NULL; + } + + npi = PyArray_DIM(pi_a, 0); + sai_a = (PyArrayObject *)PyArray_NewLikeArray(pi_a, NPY_CORDER, NULL, 0); + if (sai_a == NULL) + { + Py_DECREF(sa_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + return NULL; + } + + cti_a = (PyArrayObject *)PyArray_NewLikeArray(pi_a, NPY_CORDER, NULL, 0); + if (cti_a == NULL) + { + Py_DECREF(sa_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + Py_DECREF(sai_a); + return NULL; + } + + ret = gsw_sa_ct_interp((double *)PyArray_DATA(sa_a), + (double *)PyArray_DATA(ct_a), + (double *)PyArray_DATA(p_a), + np, + (double *)PyArray_DATA(pi_a), + npi, + (double *)PyArray_DATA(sai_a), + (double *)PyArray_DATA(cti_a)); + Py_DECREF(sa_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + + if (ret) + { + PyErr_Format(PyExc_RuntimeError, + "gsw_sa_ct_interp failed with code %d; check input arguments", + ret); + Py_DECREF(sai_a); + Py_DECREF(cti_a); + return NULL; + } + + res = PyTuple_Pack(2, sai_a, cti_a); + + return res; +} diff --git a/src/method_def_entries.c b/src/method_def_entries.c index 0fe48fd..a768898 100644 --- a/src/method_def_entries.c +++ b/src/method_def_entries.c @@ -5,3 +5,5 @@ "geostrophic streamfunction dynamic height"}, {"util_pchip_interp", util_pchip_interp, METH_VARARGS, "PCHIP interpolation"}, +{"sa_ct_interp", sa_ct_interp, METH_VARARGS, + "SA and CT interpolation"}, From cbc21c81cd9ff76de489bdae7a373748ddeb37f1 Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Fri, 6 Sep 2024 18:06:06 -0700 Subject: [PATCH 04/17] Add tracer_ct_interp function --- src/method_bodies.c | 102 +++++++++++++++++++++++++++++++++++++++ src/method_def_entries.c | 2 + 2 files changed, 104 insertions(+) diff --git a/src/method_bodies.c b/src/method_bodies.c index 3fb953e..c4bf92b 100644 --- a/src/method_bodies.c +++ b/src/method_bodies.c @@ -319,3 +319,105 @@ sa_ct_interp(PyObject *NPY_UNUSED(self), PyObject *args) return res; } + + +static PyObject * +tracer_ct_interp(PyObject *NPY_UNUSED(self), PyObject *args) +{ + PyObject *tracer_o, *ct_o, *p_o, *pi_o, *res; + double factor; + PyArrayObject *tracer_a, *ct_a, *p_a, *pi_a, *traceri_a, *cti_a; + int np, npi; + int ret = 1; /* error (1) until set to 0 by the C function */ + + if (!PyArg_ParseTuple(args, "OOOOd", &tracer_o, &ct_o, &p_o, &pi_o, &factor)) + return NULL; + + tracer_a = (PyArrayObject *)PyArray_ContiguousFromAny(tracer_o, NPY_DOUBLE, 1, 1); + if (tracer_a == NULL) + return NULL; + + ct_a = (PyArrayObject *)PyArray_ContiguousFromAny(ct_o, NPY_DOUBLE, 1, 1); + if (ct_a == NULL) + { + Py_DECREF(tracer_a); + return NULL; + } + p_a = (PyArrayObject *)PyArray_ContiguousFromAny(p_o, NPY_DOUBLE, 1, 1); + if (p_a == NULL) + { + Py_DECREF(tracer_a); + Py_DECREF(ct_a); + return NULL; + } + pi_a = (PyArrayObject *)PyArray_ContiguousFromAny(pi_o, NPY_DOUBLE, 1, 1); + if (pi_a == NULL) + { + Py_DECREF(tracer_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + return NULL; + } + + np = PyArray_DIM(tracer_a, 0); + if (PyArray_DIM(ct_a, 0) != np || PyArray_DIM(p_a, 0) != np) + { + PyErr_SetString(PyExc_ValueError, + "Arguments tracer, CT, and p must have the same dimensions."); + Py_DECREF(tracer_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + return NULL; + } + + npi = PyArray_DIM(pi_a, 0); + traceri_a = (PyArrayObject *)PyArray_NewLikeArray(pi_a, NPY_CORDER, NULL, 0); + if (traceri_a == NULL) + { + Py_DECREF(tracer_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + return NULL; + } + + cti_a = (PyArrayObject *)PyArray_NewLikeArray(pi_a, NPY_CORDER, NULL, 0); + if (cti_a == NULL) + { + Py_DECREF(tracer_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + Py_DECREF(traceri_a); + return NULL; + } + + ret = gsw_tracer_ct_interp((double *)PyArray_DATA(tracer_a), + (double *)PyArray_DATA(ct_a), + (double *)PyArray_DATA(p_a), + np, + (double *)PyArray_DATA(pi_a), + npi, + factor, + (double *)PyArray_DATA(traceri_a), + (double *)PyArray_DATA(cti_a)); + Py_DECREF(tracer_a); + Py_DECREF(ct_a); + Py_DECREF(p_a); + Py_DECREF(pi_a); + + if (ret) + { + PyErr_Format(PyExc_RuntimeError, + "gsw_tracer_ct_interp failed with code %d; check input arguments", + ret); + Py_DECREF(traceri_a); + Py_DECREF(cti_a); + return NULL; + } + + res = PyTuple_Pack(2, traceri_a, cti_a); + + return res; +} diff --git a/src/method_def_entries.c b/src/method_def_entries.c index a768898..664693f 100644 --- a/src/method_def_entries.c +++ b/src/method_def_entries.c @@ -7,3 +7,5 @@ "PCHIP interpolation"}, {"sa_ct_interp", sa_ct_interp, METH_VARARGS, "SA and CT interpolation"}, +{"tracer_ct_interp", tracer_ct_interp, METH_VARARGS, + "Tracer and CT interpolation"}, From 70d75c7686c28a43f77c9f8ca4a92e23453c5981 Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Fri, 13 Sep 2024 13:44:49 -0700 Subject: [PATCH 05/17] Update GSW-C source code --- src/c_gsw/gsw_oceanographic_toolbox.c | 35 +++++++++++++++++---------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/src/c_gsw/gsw_oceanographic_toolbox.c b/src/c_gsw/gsw_oceanographic_toolbox.c index 45590c4..c0d355a 100644 --- a/src/c_gsw/gsw_oceanographic_toolbox.c +++ b/src/c_gsw/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,8 +8664,24 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m, sa_i[i] = NAN; ct_i[i] = NAN; } - + // Find NaNs in profile + prof_len = 0; + for (i=0; i 0) { gsw_util_sort_real(p_obs, prof_len, p_idx); for (i=0; i Date: Fri, 13 Sep 2024 14:07:23 -0700 Subject: [PATCH 06/17] Add sa_ct_interp to GSW-Python --- gsw/__init__.py | 1 + gsw/interpolation.py | 92 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 gsw/interpolation.py diff --git a/gsw/__init__.py b/gsw/__init__.py index 937212c..b814eda 100644 --- a/gsw/__init__.py +++ b/gsw/__init__.py @@ -35,6 +35,7 @@ from ._fixed_wrapped_ufuncs import * # noqa from .conversions import t90_from_t68 from .geostrophy import * # noqa +from .interpolation import * # noqa from .stability import * # noqa from .utility import * # noqa diff --git a/gsw/interpolation.py b/gsw/interpolation.py new file mode 100644 index 0000000..6b8687e --- /dev/null +++ b/gsw/interpolation.py @@ -0,0 +1,92 @@ +""" +Functions for vertical interpolation. +""" + +import numpy as np + +from . import _gsw_ufuncs +from ._utilities import indexer, match_args_return + +__all__ = ['sa_ct_interp', + ] + +@match_args_return +def sa_ct_interp(SA, CT, p, p_i, axis=0): + """ + Interpolates vertical casts of values of Absolute Salinity + and Conservative Temperature to the arbitrary pressures p_i. + + Parameters + ---------- + SA : array-like + Absolute Salinity, g/kg + CT : array-like + Conservative Temperature (ITS-90), degrees C + p : array-like + Sea pressure (absolute pressure minus 10.1325 dbar), dbar + p_i : array-like + Sea pressure to interpolate on, dbar + axis : int, optional, default is 0 + The index of the pressure dimension in SA and CT. + + + Returns + ------- + SA_i : array + Values of SA interpolated to p_i along the specified axis. + CT_i : array + Values of CT interpolated to p_i along the specified axis. + + """ + if SA.shape != CT.shape: + raise ValueError(f'Shapes of SA and CT must match; found {SA.shape} and {CT.shape}') + if p.ndim != p_i.ndim: + raise ValueError(f'p and p_i must have the same number of dimensions;\n' + f' found {p.ndim} versus {p_i.ndim}') + if p.ndim == 1 and SA.ndim > 1: + if len(p) != SA.shape[axis]: + raise ValueError( + f'With 1-D p, len(p) must be SA.shape[axis];\n' + f' found {len(p)} versus {SA.shape[axis]} on specified axis, {axis}' + ) + ind = [np.newaxis] * SA.ndim + ind[axis] = slice(None) + p = p[tuple(ind)] + p_i = p_i[tuple(ind)] + elif p.ndim > 1: + if p.shape != SA.shape: + raise ValueError(f'With {p.ndim}-D p, shapes of p and SA must match;\n' + f'found {p.shape} and {SA.shape}') + if any([p.shape[i] != p_i.shape[i] for i in range(p.ndim) if i != axis]): + raise ValueError(f'With {p.ndim}-D p, p and p_i must have the same dimensions outside of axis {axis};\n' + f' found {p.shape} versus {p_i.shape}') + with np.errstate(invalid='ignore'): + # The need for this context seems to be a bug in np.ma.any. + if np.ma.any(np.ma.diff(np.ma.masked_invalid(p_i), axis=axis) <= 0) \ + or np.ma.any(np.ma.diff(np.ma.masked_invalid(p), axis=axis) <= 0): + raise ValueError('p and p_i must be increasing along the specified axis') + p = np.broadcast_to(p, SA.shape) + goodmask = ~(np.isnan(SA) | np.isnan(CT) | np.isnan(p)) + SA_i = np.empty(p_i.shape, dtype=float) + CT_i = np.empty(p_i.shape, dtype=float) + SA_i.fill(np.nan) + CT_i.fill(np.nan) + + try: + order = 'F' if SA.flags.fortran else 'C' + except AttributeError: + order = 'C' # e.g., xarray DataArray doesn't have flags + for ind in indexer(SA.shape, axis, order=order): + # this is needed to support xarray inputs for numpy < 1.23 + igood = np.asarray(goodmask[ind]) + pgood = p[ind][igood] + pi = p_i[ind] + # There must be at least 2 non-NaN values for interpolation + if len(pgood) > 2: + sa = SA[ind][igood] + ct = CT[ind][igood] + sai, cti = _gsw_ufuncs.sa_ct_interp(sa, ct, pgood, pi) + SA_i[ind] = sai + CT_i[ind] = cti + + return (SA_i, CT_i) From 2cf41a9ab850472773fb1217d3d53b50a9b61485 Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Wed, 18 Sep 2024 09:52:50 -0700 Subject: [PATCH 07/17] Add tracer_ct_interp to GSW-Python --- gsw/interpolation.py | 85 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) diff --git a/gsw/interpolation.py b/gsw/interpolation.py index 6b8687e..555074f 100644 --- a/gsw/interpolation.py +++ b/gsw/interpolation.py @@ -8,6 +8,7 @@ from ._utilities import indexer, match_args_return __all__ = ['sa_ct_interp', + 'tracer_ct_interp', ] @match_args_return @@ -90,3 +91,87 @@ def sa_ct_interp(SA, CT, p, p_i, axis=0): CT_i[ind] = cti return (SA_i, CT_i) + +@match_args_return +def tracer_ct_interp(tracer, CT, p, p_i, factor=9., axis=0): + """ + Interpolates vertical casts of values of a tracer + and Conservative Temperature to the arbitrary pressures p_i. + + Parameters + ---------- + tracer : array-like + tracer + CT : array-like + Conservative Temperature (ITS-90), degrees C + p : array-like + Sea pressure (absolute pressure minus 10.1325 dbar), dbar + p_i : array-like + Sea pressure to interpolate on, dbar + factor: float, optional, default is 9. + Ratio between the ranges of Conservative Temperature + and tracer in the world ocean. + axis : int, optional, default is 0 + The index of the pressure dimension in tracer and CT. + + + Returns + ------- + tracer_i : array + Values of tracer interpolated to p_i along the specified axis. + CT_i : array + Values of CT interpolated to p_i along the specified axis. + + """ + if tracer.shape != CT.shape: + raise ValueError(f'Shapes of tracer and CT must match; found {tracer.shape} and {CT.shape}') + if p.ndim != p_i.ndim: + raise ValueError(f'p and p_i must have the same number of dimensions;\n' + f' found {p.ndim} versus {p_i.ndim}') + if p.ndim == 1 and tracer.ndim > 1: + if len(p) != tracer.shape[axis]: + raise ValueError( + f'With 1-D p, len(p) must be tracer.shape[axis];\n' + f' found {len(p)} versus {tracer.shape[axis]} on specified axis, {axis}' + ) + ind = [np.newaxis] * tracer.ndim + ind[axis] = slice(None) + p = p[tuple(ind)] + p_i = p_i[tuple(ind)] + elif p.ndim > 1: + if p.shape != tracer.shape: + raise ValueError(f'With {p.ndim}-D p, shapes of p and tracer must match;\n' + f'found {p.shape} and {tracer.shape}') + if any([p.shape[i] != p_i.shape[i] for i in range(p.ndim) if i != axis]): + raise ValueError(f'With {p.ndim}-D p, p and p_i must have the same dimensions outside of axis {axis};\n' + f' found {p.shape} versus {p_i.shape}') + with np.errstate(invalid='ignore'): + # The need for this context seems to be a bug in np.ma.any. + if np.ma.any(np.ma.diff(np.ma.masked_invalid(p_i), axis=axis) <= 0) \ + or np.ma.any(np.ma.diff(np.ma.masked_invalid(p), axis=axis) <= 0): + raise ValueError('p and p_i must be increasing along the specified axis') + p = np.broadcast_to(p, tracer.shape) + goodmask = ~(np.isnan(tracer) | np.isnan(CT) | np.isnan(p)) + tracer_i = np.empty(p_i.shape, dtype=float) + CT_i = np.empty(p_i.shape, dtype=float) + tracer_i.fill(np.nan) + CT_i.fill(np.nan) + + try: + order = 'F' if tracer.flags.fortran else 'C' + except AttributeError: + order = 'C' # e.g., xarray DataArray doesn't have flags + for ind in indexer(tracer.shape, axis, order=order): + # this is needed to support xarray inputs for numpy < 1.23 + igood = np.asarray(goodmask[ind]) + pgood = p[ind][igood] + pi = p_i[ind] + # There must be at least 2 non-NaN values for interpolation + if len(pgood) > 2: + tr = tracer[ind][igood] + ct = CT[ind][igood] + tri, cti = _gsw_ufuncs.sa_ct_interp(tr, ct, pgood, pi, factor) + tracer_i[ind] = tri + CT_i[ind] = cti + + return (tracer_i, CT_i) From 79ee4f176d13d91fabbc296fb9b21628beda4d9e Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Wed, 18 Sep 2024 10:21:29 -0700 Subject: [PATCH 08/17] Use correct interp function in tracer_ct_interp --- gsw/interpolation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsw/interpolation.py b/gsw/interpolation.py index 555074f..8726478 100644 --- a/gsw/interpolation.py +++ b/gsw/interpolation.py @@ -170,7 +170,7 @@ def tracer_ct_interp(tracer, CT, p, p_i, factor=9., axis=0): if len(pgood) > 2: tr = tracer[ind][igood] ct = CT[ind][igood] - tri, cti = _gsw_ufuncs.sa_ct_interp(tr, ct, pgood, pi, factor) + tri, cti = _gsw_ufuncs.tracer_ct_interp(tr, ct, pgood, pi, factor) tracer_i[ind] = tri CT_i[ind] = cti From 8aa9728184144f418b29921598464ec808f94c1f Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Wed, 18 Sep 2024 10:22:08 -0700 Subject: [PATCH 09/17] Add tests for interpolation functions --- gsw/tests/test_interpolation.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 gsw/tests/test_interpolation.py diff --git a/gsw/tests/test_interpolation.py b/gsw/tests/test_interpolation.py new file mode 100644 index 0000000..fe9ad37 --- /dev/null +++ b/gsw/tests/test_interpolation.py @@ -0,0 +1,29 @@ +import os + +import numpy as np +from numpy.testing import assert_allclose + +import gsw +from gsw._utilities import Bunch + +root_path = os.path.abspath(os.path.dirname(__file__)) + +cv = Bunch(np.load(os.path.join(root_path, 'gsw_cv_v3_0.npz'))) + +def test_sa_ct_interp(): + p = cv.p_chck_cast + CT = cv.CT_chck_cast + SA = cv.SA_chck_cast + p_i = np.repeat(cv.p_i[:, np.newaxis], p.shape[1], axis=1) + SA_i, CT_i = gsw.sa_ct_interp(SA, CT, p, p_i) + assert_allclose(SA_i, cv.SAi_SACTinterp, rtol=0, atol=cv.SAi_SACTinterp_ca) + assert_allclose(CT_i, cv.CTi_SACTinterp, rtol=0, atol=cv.CTi_SACTinterp_ca) + +def test_tracer_ct_interp(): + p = cv.p_chck_cast + CT = cv.CT_chck_cast + tracer = cv.SA_chck_cast + p_i = np.repeat(cv.p_i[:, np.newaxis], p.shape[1], axis=1) + tracer_i, CT_i = gsw.tracer_ct_interp(tracer, CT, p, p_i) + assert_allclose(tracer_i, cv.traceri_tracerCTinterp, rtol=0, atol=cv.traceri_tracerCTinterp_ca) + assert_allclose(CT_i, cv.CTi_SACTinterp, rtol=0, atol=cv.CTi_SACTinterp_ca) From 5edaf1922dd062bec54553a0ffe88df592d3d212 Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Tue, 24 Sep 2024 17:18:39 -0700 Subject: [PATCH 10/17] Update GSW-C source --- src/c_gsw/gsw_oceanographic_toolbox.c | 55 ++++++++++++++++----------- src/c_gsw/gswteos-10.h | 2 +- 2 files changed, 33 insertions(+), 24 deletions(-) diff --git a/src/c_gsw/gsw_oceanographic_toolbox.c b/src/c_gsw/gsw_oceanographic_toolbox.c index c0d355a..1bda1f3 100644 --- a/src/c_gsw/gsw_oceanographic_toolbox.c +++ b/src/c_gsw/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; @@ -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)); @@ -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; @@ -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' @@ -11401,8 +11393,24 @@ 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 0) { gsw_util_sort_real(p_obs, prof_len, p_idx); for (i=0; i Date: Thu, 26 Sep 2024 01:34:56 -0700 Subject: [PATCH 11/17] Update GSW-C source --- src/c_gsw/gsw_oceanographic_toolbox.c | 14 ++++++++++---- src/c_gsw/gswteos-10.h | 1 + 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/c_gsw/gsw_oceanographic_toolbox.c b/src/c_gsw/gsw_oceanographic_toolbox.c index 1bda1f3..de545aa 100644 --- a/src/c_gsw/gsw_oceanographic_toolbox.c +++ b/src/c_gsw/gsw_oceanographic_toolbox.c @@ -4010,8 +4010,9 @@ gsw_geo_strf_dyn_height_1(double *sa, double *ct, double *p, double p_ref, ! ! This function evaluates the pressure integral of specific volume using ! SA and CT interpolated with respect to pressure. The interpolation method -! may be chosen as linear or "PCHIP", piecewise cubic Hermite using a shape- -! preserving algorithm for setting the derivatives. +! may be chosen as linear, "PCHIP" (piecewise cubic Hermite using a shape- +! preserving algorithm for setting the derivatives), or "MRST-PCHIP" +! (Multiply-Rotated Salinity–Temperature PCHIP). ! ! SA = Absolute Salinity [ g/kg ] ! CT = Conservative Temperature (ITS-90) [ deg C ] @@ -4023,7 +4024,7 @@ gsw_geo_strf_dyn_height_1(double *sa, double *ct, double *p, double p_ref, ! geo_strf_dyn_height = dynamic height anomaly [ m^2/s^2 ] ! max_dp_i = maximum pressure difference between points for triggering ! interpolation. -! interp_method = 1 for linear, 2 for PCHIP +! interp_method = 1 for linear, 2 for PCHIP, 3 for MRST-PCHIP ! ! Note. If p_ref falls outside the range of a ! vertical profile, the dynamic height anomaly for each bottle @@ -4269,8 +4270,13 @@ gsw_geo_strf_dyn_height_1(double *sa, double *ct, double *p, double p_ref, err = err || gsw_util_pchip_interp(p, ct, nz, p_i, ct_i, n_i); if (err) err = 6; } + else if (interp_method == INTERP_METHOD_MRST) { + err = gsw_sa_ct_interp(sa, ct, p, nz, + p_i, n_i, sa_i, ct_i); + if (err) err = 7; + } else { - err = 7; + err = 8; } if (err) { free(p_i); diff --git a/src/c_gsw/gswteos-10.h b/src/c_gsw/gswteos-10.h index e273cad..4fb1e76 100644 --- a/src/c_gsw/gswteos-10.h +++ b/src/c_gsw/gswteos-10.h @@ -35,6 +35,7 @@ extern "C" { #define INTERP_METHOD_LINEAR 1 #define INTERP_METHOD_PCHIP 2 +#define INTERP_METHOD_MRST 3 /* ** Prototypes: From 2a650263c55e4519da026ea6687d947b5e289fec Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Mon, 7 Oct 2024 17:02:48 -0700 Subject: [PATCH 12/17] Add MRST as an interpolation method option to geo_strf_dyn_height --- gsw/geostrophy.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gsw/geostrophy.py b/gsw/geostrophy.py index d8c5fda..501f01e 100644 --- a/gsw/geostrophy.py +++ b/gsw/geostrophy.py @@ -36,7 +36,7 @@ def geo_strf_dyn_height(SA, CT, p, p_ref=0, axis=0, max_dp=1.0, If any pressure interval in the input p exceeds max_dp, the dynamic height will be calculated after interpolating to a grid with this spacing. - interp_method : string {'pchip', 'linear'} + interp_method : string {'mrst', 'pchip', 'linear'} Interpolation algorithm. Returns @@ -48,7 +48,7 @@ def geo_strf_dyn_height(SA, CT, p, p_ref=0, axis=0, max_dp=1.0, in an isobaric surface, relative to the reference surface. """ - interp_methods = {'pchip' : 2, 'linear' : 1} + interp_methods = {'mrst' : 3, 'pchip' : 2, 'linear' : 1} if interp_method not in interp_methods: raise ValueError(f'interp_method must be one of {interp_methods.keys()}') if SA.shape != CT.shape: From 489b8ed6a2515c501935db0c02cfd8c276e7d5af Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Mon, 7 Oct 2024 17:43:15 -0700 Subject: [PATCH 13/17] Add test for using geo_strf_dyn_height with the MRST-PCHIP interpolation method --- gsw/tests/test_geostrophy.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/gsw/tests/test_geostrophy.py b/gsw/tests/test_geostrophy.py index a6fefaa..5f427de 100644 --- a/gsw/tests/test_geostrophy.py +++ b/gsw/tests/test_geostrophy.py @@ -1,7 +1,7 @@ import os import numpy as np -from numpy.testing import assert_almost_equal, assert_array_equal +from numpy.testing import assert_almost_equal, assert_array_equal, assert_allclose import gsw from gsw._utilities import Bunch @@ -102,3 +102,14 @@ def test_pz_roundtrip(): zz = gsw.z_from_p(p, 30, 0.5, 0.25) assert_almost_equal(z, zz) +def test_dyn_height_mrst(): + """ + Tests the MRST-PCHIP interpolation method. + """ + p = cv.p_chck_cast + CT = cv.CT_chck_cast + SA = cv.SA_chck_cast + pr = cv.pr + strf = gsw.geo_strf_dyn_height(SA, CT, p, p_ref=pr, interp_method='mrst') + + assert_allclose(strf, cv.geo_strf_dyn_height, rtol=0, atol=cv.geo_strf_dyn_height_ca) From 3bd5d3f4008d7b52e459f9bea2d6e85eaee3b7e0 Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Wed, 6 Nov 2024 08:44:49 -0800 Subject: [PATCH 14/17] Remove unnecessary list comprehensions --- gsw/interpolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gsw/interpolation.py b/gsw/interpolation.py index 8726478..a409066 100644 --- a/gsw/interpolation.py +++ b/gsw/interpolation.py @@ -58,7 +58,7 @@ def sa_ct_interp(SA, CT, p, p_i, axis=0): if p.shape != SA.shape: raise ValueError(f'With {p.ndim}-D p, shapes of p and SA must match;\n' f'found {p.shape} and {SA.shape}') - if any([p.shape[i] != p_i.shape[i] for i in range(p.ndim) if i != axis]): + if any(p.shape[i] != p_i.shape[i] for i in range(p.ndim) if i != axis): raise ValueError(f'With {p.ndim}-D p, p and p_i must have the same dimensions outside of axis {axis};\n' f' found {p.shape} versus {p_i.shape}') with np.errstate(invalid='ignore'): @@ -142,7 +142,7 @@ def tracer_ct_interp(tracer, CT, p, p_i, factor=9., axis=0): if p.shape != tracer.shape: raise ValueError(f'With {p.ndim}-D p, shapes of p and tracer must match;\n' f'found {p.shape} and {tracer.shape}') - if any([p.shape[i] != p_i.shape[i] for i in range(p.ndim) if i != axis]): + if any(p.shape[i] != p_i.shape[i] for i in range(p.ndim) if i != axis): raise ValueError(f'With {p.ndim}-D p, p and p_i must have the same dimensions outside of axis {axis};\n' f' found {p.shape} versus {p_i.shape}') with np.errstate(invalid='ignore'): From cf40756019f705dcc76767772ccfda60c2a661d1 Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Wed, 6 Nov 2024 08:51:17 -0800 Subject: [PATCH 15/17] Apply isort to fix import format --- gsw/tests/test_geostrophy.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gsw/tests/test_geostrophy.py b/gsw/tests/test_geostrophy.py index 5f427de..f72c85b 100644 --- a/gsw/tests/test_geostrophy.py +++ b/gsw/tests/test_geostrophy.py @@ -1,7 +1,8 @@ import os import numpy as np -from numpy.testing import assert_almost_equal, assert_array_equal, assert_allclose +from numpy.testing import (assert_allclose, assert_almost_equal, + assert_array_equal) import gsw from gsw._utilities import Bunch From 12a9564db7e484e70270eb1df36185f24bb644be Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Wed, 6 Nov 2024 12:46:18 -0800 Subject: [PATCH 16/17] Fix import ordering with ruff --- gsw/tests/test_geostrophy.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gsw/tests/test_geostrophy.py b/gsw/tests/test_geostrophy.py index f72c85b..e0d4423 100644 --- a/gsw/tests/test_geostrophy.py +++ b/gsw/tests/test_geostrophy.py @@ -1,8 +1,7 @@ import os import numpy as np -from numpy.testing import (assert_allclose, assert_almost_equal, - assert_array_equal) +from numpy.testing import assert_allclose, assert_almost_equal, assert_array_equal import gsw from gsw._utilities import Bunch From 8e59755ddc2f26d985a4dc60c572611859cf239f Mon Sep 17 00:00:00 2001 From: mauzey1 Date: Mon, 11 Nov 2024 09:43:43 -0800 Subject: [PATCH 17/17] Update GSW-C source code --- gsw/_wrapped_ufuncs.py | 4 ++-- src/c_gsw/gsw_oceanographic_toolbox.c | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/gsw/_wrapped_ufuncs.py b/gsw/_wrapped_ufuncs.py index ebce253..e2bc1f9 100644 --- a/gsw/_wrapped_ufuncs.py +++ b/gsw/_wrapped_ufuncs.py @@ -5135,7 +5135,7 @@ def rho(SA, CT, p): Returns ------- - rho : array-like, kg/m^3^3^3^3 + rho : array-like, kg/m^3 in-situ density @@ -5194,7 +5194,7 @@ def rho_alpha_beta(SA, CT, p): Returns ------- - rho : array-like, kg/m^3^3^3^3 + rho : array-like, kg/m^3 in-situ density alpha : array-like, 1/K thermal expansion coefficient diff --git a/src/c_gsw/gsw_oceanographic_toolbox.c b/src/c_gsw/gsw_oceanographic_toolbox.c index de545aa..36f1918 100644 --- a/src/c_gsw/gsw_oceanographic_toolbox.c +++ b/src/c_gsw/gsw_oceanographic_toolbox.c @@ -8639,7 +8639,7 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m, 1.950903220161283e-1 }; - int i, j, k, prof_len, + int i, 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; @@ -8908,9 +8908,9 @@ gsw_sa_ct_interp(double *sa, double *ct, double *p, int m, break; } - for (i=0; i < i_2_len; ++i) { - if ((i_3[i] - i_frozen) <= 0) { - i_shallower = i; + for (i_shallower=i_2_len-1; i_shallower>=0; --i_shallower) { + if ((i_3[i_shallower] - i_frozen) <= 0) { + break; } } i_above = i_2[i_shallower]; @@ -11369,13 +11369,13 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m, 1.950903220161283e-1 }; - int i, j, k, prof_len, + int i, 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_out_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 d, 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, @@ -11550,7 +11550,7 @@ gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m, } } 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); + 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));