Linearly interpolate CTD S, t, and p (optional) from `z` to `z_new`.

    Args
    ----
    S: ndarray
        CTD salinities
    t: ndarray
        CTD temperatures
    z: ndarray
        CTD Depths, must be a strictly increasing or decreasing 1-D array, and
        its length must match the last dimension of `S` and `t`.
    z_new: ndarray
        Depth to interpolate `S`, `t`, and `p` to. May be a scalar or a sequence.
    p: ndarray (optional)
        CTD pressures

    Returns
    -------
    S_i: ndarray
        Interpolated salinities
    t_i: ndarray
        Interpolated temperatures
    p_i: ndarray
        Interpolated pressures. `None` returned if `p` not passed.

    Note
    ----
    It is assumed, but not checked, that `S`, `t`, and `z` are all plain
    ndarrays, not masked arrays or other sequences.

    Out-of-range values of `z_new`, and `nan` in `S` and `t`, yield
    corresponding `numpy.nan` in the output.

    This method is adapted from the the legacy `python-gsw` package, where
    their basic algorithm is from scipy.interpolate.

def interp_S_t(S, t, z, z_new, p=None):
    ''' Linearly interpolate CTD S, t, and p (optional) from `z` to `z_new`.

    Args
    ----
    S: ndarray
        CTD salinities
    t: ndarray
        CTD temperatures
    z: ndarray
        CTD Depths, must be a strictly increasing or decreasing 1-D array, and
        its length must match the last dimension of `S` and `t`.
    z_new: ndarray
        Depth to interpolate `S`, `t`, and `p` to. May be a scalar or a sequence.
    p: ndarray (optional)
        CTD pressures

    Returns
    -------
    S_i: ndarray
        Interpolated salinities
    t_i: ndarray
        Interpolated temperatures
    p_i: ndarray
        Interpolated pressures. `None` returned if `p` not passed.

    Note
    ----
    It is assumed, but not checked, that `S`, `t`, and `z` are all plain
    ndarrays, not masked arrays or other sequences.

    Out-of-range values of `z_new`, and `nan` in `S` and `t`, yield
    corresponding `numpy.nan` in the output.

    This method is adapted from the the legacy `python-gsw` package, where
    their basic algorithm is from scipy.interpolate.
    '''
    import numpy

    # Create array-like query depth if single value passed
    isscalar = False
    if not numpy.iterable(z_new):
        isscalar = True
        z_new = [z_new]

    # Type cast to numpy array
    z_new = numpy.asarray(z_new)

    # Determine if depth direction is inverted
    inverted = False
    if z[1] - z[0] < 0:
        inverted = True
        z = z[::-1]
        S = S[..., ::-1]
        t = t[..., ::-1]
        if p is not None:
            p = p[..., ::-1]

    # Ensure query depths are ordered
    if (numpy.diff(z) <= 0).any():
        raise ValueError("z must be strictly increasing or decreasing")

    # Find z indices where z_new elements can insert with maintained order
    hi = numpy.searchsorted(z, z_new)

    # Replaces indices below/above with 1 or len(z)-1
    hi = hi.clip(1, len(z) - 1).astype(int)
    lo = hi - 1

    # Get CTD depths above and below query depths
    z_lo = z[lo]
    z_hi = z[hi]
    S_lo = S[lo]
    S_hi = S[hi]
    t_lo = t[lo]
    t_hi = t[hi]

    # Calculate distance ratio between CTD depths
    z_ratio = (z_new - z_lo) / (z_hi - z_lo)

    # Interpolate salinity and temperature with delta and ratio
    S_i = S_lo + (S_hi - S_lo) * z_ratio
    t_i = t_lo + (t_hi - t_lo) * z_ratio
    if p is None:
        p_i = None
    else:
        p_i = p[lo] + (p[hi] - p[lo]) * z_ratio

    # Invert interp values if passed depths inverted
    if inverted:
        S_i = S_i[..., ::-1]
        t_i = t_i[..., ::-1]
        if p is not None:
            p_i = p_i[..., ::-1]

    # Replace values not within CTD sample range with `nan`s
    outside = (z_new < z.min()) | (z_new > z.max())
    if numpy.any(outside):
        S_i[..., outside] = numpy.nan
        t_i[..., outside] = numpy.nan
        if p is not None:
            p_i[..., outside] = numpy.nan

    # Return single interp values if single query depth passed
    if isscalar:
        S_i = S_i[0]
        t_i = t_i[0]
        if p is not None:
            p_i = p_i[0]

    return S_i, t_i, p_i