R/sw.R

In dankelley/oce: Analysis of Oceanographic Data

# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4

#' Alter Longitude and Latitude for gsw Computations
#'
#' This function repeats location information as required by some seawater
#' functions, e.g. [swAbsoluteSalinity()], that use the `gsw` package to compute
#' seawater properties in the Gibbs Seawater formulation. It seems unlikely that
#' users will need to call this function directly in routine work.
#'
#' Several `gsw` functions require location information to be matched up with
#' hydrographic information.  The scheme depends on the dimensionality of the
#' hydrographic variables and the location variables. For example, the
#' [ctd-class] stores `salinity` etc in vectors, an stores just one
#' longitude-latitude pair for each vector.  By contrast, the [argo-class]
#' stores `salinity` etc as matrices, and stores e.g. `longitude` as a vector of
#' length matching the first dimension of `salinity`.
#'
#' @param x an [oce-class] object.
#'
#' @return `locationForGsw` returns a list containing `longitude` and
#' `latitude`, with dimensionality matching `pressure` in the `data` slot of
#' `x`.  If `x` lacks location information (in either its `metadata` or `data`
#' slot) or lacks `pressure` in its data slot, then the returned list will hold
#' NULL values for both `longitude` and `latitude`.
#'
#' @author Dan Kelley
#'
#' @family functions that calculate seawater properties
locationForGsw <- function(x) {
    if (!inherits(x, "oce")) {
        stop("x must be an oce-class object")
    }
    pressure <- x@data$pressure
    if (is.null(pressure)) {
        return(list(longitude = NULL, latitude = NULL))
    }
    longitude <- x[["longitude"]]
    latitude <- x[["latitude"]]
    if (is.null(longitude) || is.null(latitude)) {
        return(list(longitude = NULL, latitude = NULL))
    }
    if (is.array(pressure)) {
        dim <- dim(pressure)
        longitude <- rep(longitude, each = dim[1])
        dim(longitude) <- dim
        latitude <- rep(latitude, each = dim[1])
        dim(latitude) <- dim
    } else {
        np <- length(pressure)
        longitude <- rep(longitude[1], length.out = np)
        latitude <- rep(latitude[1], length.out = np)
    }
    dim <- dim(pressure)
    list(longitude = longitude, latitude = latitude)
}

#' Determine Available Derived Water Properties
#'
#' This determines what things can be derived from the supplied
#' variables.  For example, if `salinity`, `temperature`,
#' and `pressure` are supplied, then potential temperature, sound
#' speed, and several other things can be derived.  If, in addition,
#' `longitude` and `latitude` are supplied, then Absolute Salinity,
#' Conservative Temperature, and some other things can be derived.
#' Similarly, `nitrate` can be computed from `NO2+NO3` together
#' with `nitrate`, and `nitrite` can be computed from `NO2+NO3`
#' together with `nitrate`.
#' See the \dQuote{Examples} for a full listing.
#'
#' @param x a specification of the names of known variables.  This
#' may be (a) an [oce-class] object, in which case the names are
#' determined by calling [names()] on the `data` slot of `x`, or
#' (b) a vector of character values indicating the names.
#'
#' @return [computableWaterProperties()] returns a sorted
#' character vector holding the names of computable
#' water properties, or NULL, if there are no computable values.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' # Example 1
#' data(ctd)
#' computableWaterProperties(ctd)
#' # Example 2: nothing an be computed from just salinity
#' computableWaterProperties("salinity")
#' # Example 3: quite a lot can be computed from this trio of values
#' computableWaterProperties(c("salinity", "temperature", "pressure"))
#' # Example 4: now we can get TEOS-10 values as well
#' computableWaterProperties(c(
#'     "salinity", "temperature", "pressure",
#'     "longitude", "latitude"
#' ))
#'
#' @family functions that calculate seawater properties
computableWaterProperties <- function(x) {
    res <- NULL
    names <- if (inherits(x, "oce")) unique(c(names(x@metadata), names(x@data))) else x
    haveSTP <- all(c("salinity", "temperature", "pressure") %in% names)
    haveLocation <- all(c("longitude", "latitude") %in% names)
    if (haveSTP) {
        res <- c(res, c(
            "theta", paste("potential", "temperature"), "z",
            "depth", "spice", "Rrho", "RrhoSF", "sigmaTheta", "SP",
            "density", "N2", paste("sound", "speed")
        ))
        if (haveLocation) {
            res <- c(
                res, "SR", "Sstar",
                paste0("sigma", 0:4),
                "SA", paste("Absolute", "Salinity"),
                "CT", paste("Conservative", "Temperature"),
                paste0("spiciness", 0:2)
            )
        }
    }
    # It is possible to compute nitrate from NO2+NO3 and nitrite, if
    # it's not stored in the object already.  This occurs in data(section).
    # I also added a similar scheme for nitrite, although I don't know
    # whether that condition ever happens, in practice.
    if (!("nitrate" %in% names) && ("NO2+NO3" %in% names) && ("nitrite" %in% names)) {
        res <- c(res, "nitrate")
    }
    if (!("nitrite" %in% names) && ("NO2+NO3" %in% names) && ("nitrate" %in% names)) {
        res <- c(res, "nitrite")
    }
    # Remove things that were in the original data.
    res <- res[!(res %in% names)]
    # remove duplicates, and sort
    sort(unique(res))
}

#' Convert From ITS-90 to IPTS-68 Temperature
#'
#' @template temperatureConversionTemplate
#' @param temperature numeric vector of temperatures]in \eqn{^\circ}{deg}C on the ITS-90 scale.
#' @return Corresponding temperatures in \eqn{^\circ}{deg}C on the IPTS-68 scale.
T68fromT90 <- function(temperature) temperature * 1.00024

#' Convert From IPTS-68 to ITS-90 Temperature
#'
#' @template temperatureConversionTemplate
#' @param temperature numeric vector of temperatures in \eqn{^\circ}{deg}C on the IPTS-68 scale.
#' @return Corresponding temperatures in \eqn{^\circ}{deg}C on the ITS-90 scale.
T90fromT68 <- function(temperature) temperature / 1.00024

#' Convert From ITS-48 to ITS-90 Temperature
#'
#' @template temperatureConversionTemplate
#' @param temperature Vector of temperatures in \eqn{^\circ}{deg}C on the IPTS-48 scale.
#' @return Corresponding temperatures in \eqn{^\circ}{deg}C on the ITS-90 scale.
T90fromT48 <- function(temperature) (temperature - 4.4e-6 * temperature * (100 - temperature)) / 1.00024

#' Look Within the First Element of a List for Replacement Values
#'
#' This is a helper function used by some seawater functions
#' (with names starting with `sw`) to
#' facilitate the specification of water properties either with
#' distinct arguments, or with data stored within an `oce`
#' object that is the first argument.
#'
#' If `list[1]` is not an `oce` object, then the
#' return value of `lookWithin` is the same as the input
#' value, except that (a) `eos` is completed to either
#' `"gsw"` or `"unesco"` and (b) if `longitude`
#' and `latitude` are within `list[1]`, then they
#' are possibly lengthened, to have the same length as the first
#' item in the `data` slot of `list[1]`.
#'
#' The examples may clarify this somewhat.
#'
#' @param list A list of elements, typically arguments that will be used in sw functions.
#'
#' @return A list with elements of the same names but possibly filled in from the first element.
#'
#' @examples
#' # 1. If first item is not a CTD object, just return the input
#' lookWithin(list(a = 1, b = 2)) # returns a list
#' # 2. Extract salinity from a CTD object
#' data(ctd)
#' str(lookWithin(list(salinity = ctd)))
#' # 3. Extract salinity and temperature. Note that the
#' # value specified for temperature is ignored; all that matters
#' # is that temperature is named.
#' str(lookWithin(list(salinity = ctd, temperature = NULL)))
#' # 4. How it is used by swRho()
#' rho1 <- swRho(ctd, eos = "unesco")
#' rho2 <- swRho(ctd[["salinity"]], ctd[["temperature"]], ctd[["pressure"]], eos = "unesco")
#' stopifnot(all.equal(rho1, rho2))
lookWithin <- function(list) {
    n <- length(list)
    names <- names(list)
    list1 <- list[[1]]
    if (inherits(list[[1]], "oce")) {
        for (i in 1:n) {
            if ("eos" != names[i]) {
                # Note: the accessor [[]] will return temperature
                # in ITS-90 regardless of how it is stored, and similarly pressure
                # in dbar and salinity in FIXME: what unit to use??
                try(
                    {
                        list[[i]] <- list1[[names[i], "nowarn"]]
                    },
                    silent = TRUE
                )
            }
        }
        if (inherits(list1, "ctd")) {
            nrows <- length(list[[names[1]]])
            if (length(list[["longitude"]])) {
                list[["longitude"]] <- rep(mean(list[["longitude"]], na.rm = TRUE), nrows)
            }
            if (length(list[["latitude"]])) {
                list[["latitude"]] <- rep(mean(list[["latitude"]], na.rm = TRUE), nrows)
            }
        }
        # FIXME: should special-case some other object types
    }
    if ("eos" %in% names) {
        list[["eos"]] <- match.arg(list[["eos"]], c("unesco", "gsw"))
    }
    list
}


#' Density Ratio
#'
#' Compute density ratio for a `ctd` object.  An error (perhaps with some hints)
#' is issued for any other type of object.
#'
#' If `eos="unesco"`, the work is done by calculating salinity and
#' potential-temperature derivatives from smoothing splines whose properties
#' are governed by `smoothingLength` or `df`.  If
#' `sense="diffusive"` the definition is
#' \eqn{(beta*dS/dz)/(alpha*d(theta)/dz)}{(beta*dS/dz)/(alpha*d(theta)/dz)} and
#' the reciprocal for `"finger"`.
#'
#' If `eos="gsw"`, the work is done by extracting absolute salinity and
#' conservative temperature, smoothing with a smoothing spline as in the
#' `"unesco"` case, and then calling [gsw::gsw_Turner_Rsubrho()]
#' on these smoothed fields. Since the gsw function works on mid-point
#' pressures, [approx()] is used to interpolate back to the original
#' pressures.
#'
#' If the default arguments are acceptable, `ctd[["Rrho"]]` may be used
#' instead of `swRrho(ctd)`.
#'
#' @param ctd an [oce-class] object that holds `salinity`, `temperature`, and
#' `pressure`.  If `eos` is `"gsw"`, then it must also hold `longitude` and
#' `latitude`.
#'
#' @param sense an indication of the sense of double diffusion under study and
#' therefore of the definition of Rrho; see \dQuote{Details}
#'
#' @param smoothingLength ignored if `df` supplied, but otherwise the
#' latter is calculated as the number of data points, divided by the number
#' within a depth interval of `smoothingLength` metres.
#'
#' @param df if given, this is provided to [smooth.spline()].
#'
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#'
#' @template debugTemplate
#'
#' @return Density ratio defined in either the `"diffusive"` or
#' `"finger"` sense.
#'
#' @author Dan Kelley and Chantelle Layton
#'
#' @examples
#' library(oce)
#' data(ctd)
#' u <- swRrho(ctd, eos = "unesco")
#' g <- swRrho(ctd, eos = "gsw")
#' p <- ctd[["p"]]
#' plot(u, p, ylim = rev(range(p)), type = "l", xlab = expression(R[rho]))
#' lines(g, p, lty = 2, col = "red")
#' legend("topright", lty = 1:2, legend = c("unesco", "gsw"), col = c("black", "red"))
#'
#' @family functions that calculate seawater properties
swRrho <- function(
    ctd,
    sense = c("diffusive", "finger"),
    smoothingLength = 10, df,
    eos = getOption("oceEOS", default = "gsw"),
    debug = getOption("oceDebug")) {
    if (!inherits(ctd, "oce")) {
        stop("method is only for objects of class 'oce'")
    }
    oceDebug(debug, "swRho(..., sense=\"", sense, "\", eos=\"", eos, "\")\n")
    temperature <- ctd[["temperature"]]
    pressure <- ctd[["pressure"]]
    salinity <- ctd[["salinity"]] # overwrites
    location <- locationForGsw(ctd)
    longitude <- location$longitude
    latitude <- location$latitude
    sense <- match.arg(sense)
    eos <- match.arg(eos, c("unesco", "gsw"))
    if (debug > 0) {
        cat(vectorShow(salinity))
        cat(vectorShow(temperature))
        cat(vectorShow(pressure))
        cat(vectorShow(longitude))
        cat(vectorShow(latitude))
    }
    # If the data are in matrix form, split into columns and call this again
    # to compute each column.  We need this because it is derivative-based, so
    # casting the matrices into vectors will not work properly.
    if (is.matrix(salinity)) {
        rval <- matrix(NA, nrow = nrow(salinity), ncol = ncol(salinity))
        for (col in seq_len(ncol(salinity))) {
            # message("col=", col)
            CTD <- as.ctd(
                salinity = salinity[, col],
                temperature = temperature[, col],
                pressure = pressure[, col],
                longitude = longitude[, col],
                latitude = latitude[, col]
            )
            rval[, col] <- swRrho(CTD, sense = sense, debug = debug)
        }
        # message("all ok")
        return(rval)
    }
    ok <- !is.na(pressure) & !is.na(salinity) & !is.na(temperature)
    nok <- sum(ok)
    np <- length(pressure)
    # message("nok = ",nok, ", np = ",np)
    # smooth.spline issues warnings if under 4 good data, and we
    # don't want that noise.
    if (nok < 4L) {
        return(rep(NA, length.out = np))
    }
    A <- smoothingLength / mean(diff(pressure), na.rm = TRUE)
    if (missing(df)) {
        df <- nok / A
    }
    if (df > nok) {
        df <- nok / 2
    }
    # message("df=", df, ", eos=", eos)
    if (eos == "unesco") {
        theta <- ctd[["theta"]]
        ok <- !is.na(pressure) & !is.na(salinity) & !is.na(temperature)
        # infer d(theta)/dp and d(salinity)/dp from smoothing splines
        temperatureSpline <- smooth.spline(pressure[ok], temperature[ok], df = df)
        salinitySpline <- smooth.spline(pressure[ok], salinity[ok], df = df)
        # Smooth temperature and salinity to get smoothed alpha and beta
        CTD <- as.ctd(predict(salinitySpline, pressure)$y, predict(temperatureSpline, pressure)$y, pressure)
        alpha <- swAlpha(CTD, eos = "unesco")
        beta <- swBeta(CTD, eos = "unesco")
        # Using alpha ... is that right, since we have theta?
        thetaSpline <- smooth.spline(pressure[ok], theta[ok], df = df)
        dthetadp <- predict(thetaSpline, pressure, deriv = 1)$y
        dsalinitydp <- predict(salinitySpline, pressure, deriv = 1)$y
        Rrho <- if (sense == "diffusive") {
            (beta * dsalinitydp) / (alpha * dthetadp)
        } else {
            (alpha * dthetadp) / (beta * dsalinitydp)
        }
    } else if (eos == "gsw") {
        SA <- ctd[["SA"]]
        CT <- ctd[["CT"]]
        ok <- is.finite(pressure) & is.finite(SA) & is.finite(CT)
        # 'ok' is used since the gsw routine and approx have trouble with NA
        SA <- SA[ok]
        CT <- CT[ok]
        pressure <- pressure[ok]
        oceDebug(debug, "about to smooth.spline() for SA and CT, with df =", df, "\n")
        SA <- predict(smooth.spline(pressure, SA, df = df), pressure)$y
        CT <- predict(smooth.spline(pressure, CT, df = df), pressure)$y
        # message("did SA and CT splines with df=", df)
        # message("DAN 1 gsw part (sw.R:366)")
        a <- gsw::gsw_Turner_Rsubrho(SA, CT, pressure)
        # message("DAN 2")
        Rrho <- a$Rsubrho
        # message("DAN 4")
        Rrho <- approx(a$p_mid, a$Rsubrho, pressure, rule = 2)$y
        # message("DAN 5")
        if (sense == "diffusive") {
            Rrho <- 1 / Rrho
        }
        res <- rep(NA, np)
        res[ok] <- Rrho
        res[res == 9e15] <- NA
        Rrho <- res
        # message("DAN 6 (all done with gsw part)")
    }
    Rrho
}


#' Squared Buoyancy Frequency for Seawater
#'
#' Compute \eqn{N^2}{N^2}, the square of the buoyancy frequency for a seawater
#' profile.
#'
#' Smoothing is often useful prior to computing buoyancy frequency, and so this
#' may optionally be done with [smooth.spline()], unless
#' `df=NA`, in which case raw data are used.  If `df` is not
#' provided, a possibly reasonable value computed from an analysis of the
#' profile, based on the number of pressure levels.
#'
#' The core of the method involves computing potential density referenced to median
#' pressure, using the UNESCO-style [swSigmaTheta] function, and then differentiating
#' this with respect to pressure. The `derivs` argument is used
#' to control how this is done, as follows.
#'
#' * If `derivs` is not supplied, the action is as though it were
#' given as the string `"smoothing"`
#'
#' * If `derivs` equals `"simple"`, then the derivative of
#' density with respect to pressure is calculated as the ratio of first-order
#' derivatives of density and pressure, each calculated using
#' [diff()].  (A zero is appended at the top level.)
#'
#' * If `derivs` equals `"smoothing"`, then the processing
#' depends on the number of data in the profile, and on whether `df` is
#' given as an optional argument.  When the number of points exceeds 4, and
#' when `df` exceeds 1, [smooth.spline()] is used to calculate
#' smoothing spline representation the variation of density as a function of
#' pressure, and derivatives are extracted from the spline using
#' `predict`.  Otherwise, density is smoothed using [smooth()],
#' and derivatives are calculated as with the `"simple"` method.
#'
#' * If `derivs` is a function taking two arguments (first pressure,
#' then density) then that function is called directly to calculate the
#' derivative, and no smoothing is done before or after that call.
#'
#' For precise work, it makes sense to skip `swN2` entirely, choosing
#' whether, what, and how to smooth based on an understanding of fundamental
#' principles as well as data practicalities.
#'
#' @param pressure either pressure (dbar) (in which case `sigmaTheta` must
#' be provided) *or* an object of class `ctd` object (in which case
#' `sigmaTheta` is inferred from the object.
#'
#' @param sigmaTheta Surface-referenced potential density minus 1000
#' (kg/m\eqn{^3}{^3}).
#'
#' @param derivs optional argument to control how the derivative
#' \eqn{d\sigma_\theta/dp}{d(sigmaTheta)/d(pressure)} is calculated.  This may
#' be a character string or a function of two arguments.  See \dQuote{Details}.
#'
#' @param df argument passed to [smooth.spline()] if this function is
#' used for smoothing; set to `NA` to prevent smoothing.
#'
#' @param \dots additional argument, passed to [smooth.spline()], in
#' the case that `derivs="smoothing"`.  See \dQuote{Details}.
#'
#' @template debugTemplate
#'
#' @seealso The [gsw::gsw_Nsquared()] function of the \CRANpkg{gsw}
#' provides an alternative to this, as formulated in the GSW system. It
#' has a more sophisticated treatment of potential density, but it is based
#' on simple first-difference derivatives, so its results may require
#' smoothing, depending on the dataset and purpose of the analysis.
#'
#' @return Square of buoyancy frequency (\eqn{radian^2/s^2}{radian^2/s^2}).
#'
#' @author Dan Kelley
#'
#' @examples
#'
#' library(oce)
#' data(ctd)
#' # Left panel: density
#' p <- ctd[["pressure"]]
#' ylim <- rev(range(p))
#' par(mfrow = c(1, 2), mar = c(3, 3, 1, 1), mgp = c(2, 0.7, 0))
#' plot(ctd[["sigmaTheta"]], p, ylim = ylim, type = "l", xlab = expression(sigma[theta]))
#' # Right panel: N2, with default settings (black) and with df=2 (red)
#' N2 <- swN2(ctd)
#' plot(N2, p, ylim = ylim, xlab = "N2 [1/s^2]", ylab = "p", type = "l")
#' lines(swN2(ctd, df = 3), p, col = 2)
#'
#' @section Deprecation Notice:
#' Until 2019 April 11, `swN2` had an argument named `eos`. However,
#' this did not work as stated, unless the first argument was a `ctd`
#' object. Besides, the argument name was inherently deceptive, because the UNESCO
#' scheme does not specify how N2 is to be calculated.
#' Nothing is really lost by making this change, because the new default is the
#' same as was previously available with the `eos="unesco"`
#' setup, and the gsw-formulated estimate of N2 is provided
#' by [gsw::gsw_Nsquared()] in the \CRANpkg{gsw} package.
#'
#' @family functions that calculate seawater properties
swN2 <- function(pressure, sigmaTheta = NULL, derivs, df, debug = getOption("oceDebug"), ...) {
    oceDebug(debug, "swN2(...) ...\n", sep = "", unindent = 1)
    # cat("swN2(..., df=", df, ")\n",sep="")
    # useSmoothing <- !missing(df) && is.finite(df)
    if (inherits(pressure, "oce")) {
        oceDebug(debug, "swN2() with an oce object as first argument ...\n")
        p <- pressure[["pressure"]]
        pref <- 0 # median(p, na.rm = TRUE)
        sigmaTheta <- swSigmaTheta(pressure, referencePressure = pref, eos = "unesco")
        pressure <- p
    } else {
        oceDebug(debug, "swN2() with vectors for pressure, sigmaTheta\n")
    }
    #>if (inherits(pressure, "ctd")) {
    #>    if (!is.null(sigmaTheta))
    #>        warning("sigmaTheta is ignored, since first argument was a ctd object\n")
    #>    oceDebug(debug, "first parameter is a 'ctd', so get data from it\n")
    #>    pref <- median(pressure[["pressure"]], na.rm=TRUE)
    #>    oceDebug(debug, "setting referencePressure to ", pref, "\n", sep="")
    #>    sigmaTheta <- swSigmaTheta(pressure, referencePressure=pref, eos="unesco") # NOTE: UNESCO used
    #>    pressure <- pressure[["pressure"]] # over-writes pressure
    #>}
    if (missing(derivs)) {
        derivs <- "smoothing"
    }
    if (is.matrix(pressure)) {
        oceDebug(debug, "pressure is a matrix\n")
        res <- matrix(NA, nrow = nrow(pressure), ncol = ncol(pressure))
        for (icol in seq_len(ncol(p))) {
            oceDebug(debug, "   handling column ", icol)
            res[, icol] <- swN2(pressure[, icol], sigmaTheta[, icol], derivs, df, debug = debug, ...)
        }
        oceDebug(debug, "returning a matrix of N2\n")
        return(res)
    }
    ok <- !is.na(pressure) & !is.na(sigmaTheta)
    if (is.character(derivs)) {
        if (derivs == "simple") {
            sigmaThetaDeriv <- c(0, diff(sigmaTheta) / diff(pressure))
        } else if (derivs == "smoothing") {
            depthsAll <- sum(ok)
            depths <- length(unique(pressure[ok]))
            if (missing(df)) {
                df <- if (depths > 100) {
                    floor(depths / 10)
                } # at least 10
                else if (depths > 20) {
                    floor(depths / 3)
                } # at least 7
                else if (depths > 10) {
                    floor(depths / 2)
                } # at least 5
                else {
                    depths
                }
                oceDebug(debug, "df not supplied; set to ", df, ", given ", depthsAll, " depths, ", depths, " of which are distinct\n", sep = "")
            }
            if (depths > 4 && df > 5) {
                oceDebug(debug, "using smooth.spline with df=", df, "\n", sep = "")
                sigmaThetaSmooth <- smooth.spline(pressure[ok], sigmaTheta[ok], df = df)
                sigmaThetaDeriv <- rep(NA, length(pressure))
                sigmaThetaDeriv[ok] <- predict(sigmaThetaSmooth, pressure[ok], deriv = 1)$y
            } else {
                oceDebug(debug, "using smooth.spline with df not specified\n")
                sigmaThetaSmooth <- as.numeric(smooth(sigmaTheta[ok]))
                sigmaThetaDeriv <- rep(NA, length(pressure))
                sigmaThetaDeriv[ok] <- c(0, diff(sigmaThetaSmooth) / diff(pressure[ok]))
            }
        } else {
            stop("derivs must be 'simple', 'smoothing', or a function")
        }
    } else {
        if (!is.function(derivs)) {
            stop("derivs must be 'smoothing', 'simple', or a function")
        }
        sigmaThetaDeriv <- derivs(pressure, sigmaTheta)
    }
    # FIXME (DK 2016-05-04) I am not sure I like the following since it
    # uses a standardized rho_0. But it's from some official source I think.
    # Must check this. (UNESCO book?)
    res <- ifelse(ok, 9.8 * 9.8 * 1e-4 * sigmaThetaDeriv, NA)
    oceDebug(debug, "... # swN2()\n", sep = "", unindent = 1)
    res
}


#' Water Pressure
#'
#' Compute seawater pressure from depth by inverting [swDepth()]
#' using [uniroot()].
#'
#' If `eos="unesco"` this is done by numerical inversion of
#' [swDepth()] is done using [uniroot()]. If
#' `eos="gsw"`, it is done using [gsw::gsw_p_from_z()] in the
#' \CRANpkg{gsw} package.
#'
#' @param depth distance below the surface in metres.
#'
#' @param latitude Latitude in \eqn{^\circ}{deg}N.
#'
#' @param eos indication of formulation to be used, either `"unesco"` or
#' `"gsw"`.
#'
#' @return Pressure in dbar.
#'
#' @author Dan Kelley
#'
#' @references Unesco 1983. Algorithms for computation of fundamental
#' properties of seawater, 1983. *Unesco Tech. Pap. in Mar. Sci.*, No. 44,
#' 53 pp.
#'
#' @examples
#' swPressure(9712.653, 30, eos = "unesco") # 10000
#' swPressure(9712.653, 30, eos = "gsw") #  9998.863
#'
#' @family functions that calculate seawater properties
swPressure <- function(depth, latitude = 45, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(depth)) {
        stop("must supply depth")
    }
    # FIXME-gsw add gsw version
    ndepth <- length(depth)
    if (length(latitude) < ndepth) {
        latitude <- rep(latitude, ndepth)
    }
    res <- vector("numeric", ndepth)
    eos <- match.arg(eos, c("unesco", "gsw"))
    # Takes 3.55s for 15225 points
    if (eos == "unesco") {
        for (i in 1:ndepth) {
            # FIXME: this loop is slow and should be done in C, like swCStp()
            res[i] <- if (depth[i] == 0) {
                0
            } else {
                uniroot(function(p) depth[i] - swDepth(p, latitude[i], eos), interval = depth[i] * c(0.9, 1.1))$root
            }
        }
    } else if (eos == "gsw") {
        res <- gsw::gsw_p_from_z(-depth, latitude)
    } else {
        stop("eos must be 'unesco' or 'gsw'")
    }
    res
}


#' Electrical Conductivity Ratio From Salinity, Temperature and Pressure
#'
#' Compute electrical conductivity ratio based on salinity, temperature, and
#' pressure (relative to the conductivity of seawater with salinity=35,
#' temperature68=15, and pressure=0).
#'
#' If `eos="unesco"`, the calculation is done by a bisection root search
#' on the UNESCO formula relating salinity to conductivity, temperature, and
#' pressure (see [swSCTp()]).  If it is `"gsw"` then the
#' Gibbs-SeaWater formulation is used, via [gsw_C_from_SP()].
#'
#' @param salinity practical salinity, or a CTD object (in which case its
#' temperature and pressure are used, and the next two arguments are ignored)
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see the examples, as well as the
#' \dQuote{Temperature units} section in the documentation for [swRho()].
#'
#' @param pressure pressure (dbar)
#'
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#'
#' @return Conductivity ratio (unitless), i.e. the ratio of conductivity to the
#' conductivity at salinity=35, temperature=15 (IPTS-68 scale) and pressure=0,
#' which has numerical value 42.9140 mS/cm = 4.29140 S/m (see
#' Culkin and Smith, 1980, in the regression result cited at the bottom of
#' the left-hand column on page 23).
#'
#' @author Dan Kelley
#'
#' @seealso For thermal (as opposed to electrical) conductivity, see
#' [swThermalConductivity()]. For computation of salinity from
#' electrical conductivity, see [swSCTp()].
#'
#' @references
#' 1. Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for
#' computation of fundamental properties of seawater. *Unesco Technical
#' Papers in Marine Science*, *44*, 53 pp.
#'
#' 2. Culkin, F., and Norman D. Smith, 1980. Determination of the concentration of
#' potassium chloride solution having the same electrical conductivity, at 15 C
#' and infinite frequency, as standard seawater of salinity 35.0000 ppt
#' (Chlorinity 19.37394 ppt). *IEEE Journal of Oceanic Engineering*,
#' *5*, pp 22-23.
#'
#' @examples
#' stopifnot(abs(1.0 - swCSTp(35, T90fromT68(15), 0, eos = "unesco")) < 1e-7)
#' stopifnot(abs(1.0 - swCSTp(34.25045, T90fromT68(15), 2000, eos = "unesco")) < 1e-7)
#' stopifnot(abs(1.0 - swCSTp(34.25045, T90fromT68(15), 2000, eos = "gsw")) < 1e-7)
#'
#' @family functions that calculate seawater properties
swCSTp <- function(salinity, temperature = 15, pressure = 0, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        ctd <- salinity
        salinity <- ctd[["salinity"]]
        temperature <- ctd[["temperature"]]
        pressure <- ctd[["pressure"]]
    }
    dim <- dim(salinity)
    salinity <- as.vector(salinity)
    temperature <- as.vector(temperature)
    pressure <- as.vector(pressure)
    n <- length(salinity)
    if (length(temperature) != n) {
        temperature <- rep(temperature, length.out = n)
    }
    if (length(pressure) != n) {
        pressure <- rep(pressure, length.out = n)
    }
    eos <- match.arg(eos, c("unesco", "gsw"))
    if (eos == "unesco") {
        # cat("S= ", paste(salinity, collapse=" "), "\n")
        # cat("T= ", paste(temperature, collapse=" "), "\n")
        # cat("p= ", paste(pressure, collapse=" "), "\n")
        res <- .C("sw_CSTp",
            as.integer(n), as.double(salinity), as.double(T68fromT90(temperature)), as.double(pressure),
            C = double(n), NAOK = TRUE, PACKAGE = "oce"
        )$C
    } else {
        # for the use of a constant, as opposed to a function call with (35,15,0), see
        # https://github.com/dankelley/oce/issues/746
        res <- gsw::gsw_C_from_SP(SP = salinity, t = temperature, p = pressure) / 42.9140
    }
    dim(res) <- dim
    res
}


#' Practical Salinity From Electrical Conductivity, Temperature and Pressure
#'
#' Calculate salinity from what is actually measured by a CTD, *i.e.*
#' conductivity, *in-situ* temperature and pressure.  Often this is done
#' by the CTD processing software, but sometimes it is helpful to do this
#' directly, *e.g.* when there is a concern about mismatches in sensor
#' response times.
#'
#' Two variants are provided. First, if `eos` is
#' `"unesco"`, then salinity is calculated using
#' the UNESCO algorithm described by Fofonoff and Millard (1983) as in
#' reference 1. Second, if `eos` is `"gsw"`, then the
#' Gibbs-SeaWater formulation is used, via [gsw::gsw_SP_from_C()]
#' in the \CRANpkg{gsw} package. The latter starts with the same formula
#' as the former, but if this yields a Practical Salinity less than 2,
#' then the result is instead calculated using
#' formulae provided by Hill et al. (1986; reference 2), modified to match the
#' `"unesco"` value at Practical salinity equal to 2 (reference 3).
#'
#' @param conductivity a measure of conductivity (see also `conductivityUnit`)
#' or an `oce` object holding hydrographic information. In the second case,
#' all the other arguments to `swSCTp` are ignored.
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure pressure (dbar).
#'
#' @param conductivityUnit string indicating the unit used for conductivity.
#' This may be `"ratio"` or `""` (meaning conductivity ratio),
#' `"mS/cm"` or `"S/m"`.  Note that the ratio mode assumes that
#' measured conductivity has been divided by the standard conductivity
#' of 4.2914 S/m.  In dealing with unfamiliar data for which the measurement
#' unit has not been recorded, it can be sensible to try all three possibilities
#' for `conductivityUnit`, to see which yields the most sensible salinities.
#'
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#'
#' @return Practical Salinity.
#'
#' @author Dan Kelley
#'
#' @seealso For thermal (as opposed to electrical) conductivity, see
#' [swThermalConductivity()].  For computation of electrical
#' conductivity from salinity, see [swCSTp()].
#'
#' @references
#' 1. Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for
#' computation of fundamental properties of seawater. *Unesco Technical
#' Papers in Marine Science*, *44*, 53 pp.
#'
#' 2.  K. Hill, T. Dauphinee, and D. Woods. \dQuote{The Extension of the Practical
#' Salinity Scale 1978 to Low Salinities.} IEEE Journal of Oceanic Engineering 11,
#' no. 1 (January 1986): 109-12.
#' \doi{10.1109/JOE.1986.1145154}
#'
#' 3. `gsw_from_SP` online documentation, available at
#' `http://www.teos-10.org/pubs/gsw/html/gsw_C_from_SP.html`
#'
#' @examples
#' # 1. Demonstrate agreement with test value in UNESCO documents
#' swSCTp(1, T90fromT68(15), 0, eos = "unesco") # expect 35
#' # 2. Demonstrate agreement of gsw and unesco, S>2 case
#' swSCTp(1, T90fromT68(15), 0, eos = "gsw") # again, expect 35
#' # 3. Demonstrate close values even in very brackish water
#' swSCTp(0.02, 10, 100, eos = "gsw") # 0.6013981
#' swSCTp(0.02, 10, 100, eos = "unesco") # 0.6011721
#'
#' @family functions that calculate seawater properties
swSCTp <- function(
    conductivity, temperature = NULL, pressure = NULL,
    conductivityUnit, eos = getOption("oceEOS", default = "gsw")) {
    C0 <- 42.9140 # Culkin and Smith (1980)
    # FIXME-gsw add gsw version
    if (missing(conductivity)) {
        stop("must supply conductivity (which may be S or a CTD object)")
    }
    if (missing(conductivityUnit)) {
        conductivityUnit <- ""
    } else {
        if (is.list(conductivityUnit) && "unit" %in% names(conductivityUnit)) {
            conductivityUnit <- conductivityUnit$unit
        }
        if (is.expression(conductivityUnit)) {
            conductivityUnit <- as.character(conductivityUnit)
        }
        if (conductivityUnit == "ratio") {
            conductivityUnit <- ""
        }
    }
    if (conductivityUnit != "" && conductivityUnit != "mS/cm" && conductivityUnit != "S/m") {
        stop("conductivity unit must be \"\", \"mS/cm\", or \"S/m\"")
    }
    if (inherits(conductivity, "oce")) {
        if (inherits(conductivity, "rsk")) {
            ctd <- as.ctd(conductivity)
        } else {
            ctd <- conductivity
        }
        # cat("< ", paste(names(ctd@data), collapse=" "), " >\n", sep="")
        conductivity <- ctd[["conductivity"]]
        if (is.null(conductivity)) {
            stop("this CTD object has no conductivity")
        }
        # Use unit from within the object, but may be overridden after this block.
        tmp <- ctd[["conductivityUnit"]]
        if (is.list(tmp) && "unit" %in% names(tmp)) {
            conductivityUnit <- as.character(tmp$unit)
        }
        temperature <- ctd[["temperature"]]
        pressure <- ctd[["pressure"]]
    }
    if (is.list(conductivityUnit)) {
        conductivityUnit <- as.character(conductivityUnit$unit)
    }
    if (!length(conductivityUnit)) {
        conductivityUnit <- ""
    }
    if (conductivityUnit == "mS/cm") {
        conductivity <- conductivity / C0
    } else if (conductivityUnit == "S/m") {
        conductivity <- conductivity / (C0 / 10)
    } else {
        conductivity <- conductivity
    }
    # Now, "conductivity" is in ratio form
    dim <- dim(conductivity)
    nC <- length(conductivity)
    nT <- length(temperature)
    if (nC != nT) {
        stop("lengths of conductivity and temperature must agree, but they are ", nC, " and ", nT)
    }
    if (is.null(pressure)) {
        pressure <- rep(0, nC)
    }
    np <- length(pressure)
    if (nC != np) {
        stop("lengths of conductivity and pressure must agree, but they are ", nC, " and ", np)
    }
    if (eos == "unesco") {
        #> message("swSCTp() unesco; conductivity[1]=", conductivity[1], ", temperature[1]=", temperature[1], ", pressure[1]=", pressure[1])
        res <- .C("sw_salinity",
            as.integer(nC),
            as.double(conductivity),
            as.double(T68fromT90(temperature)), # original formula is in IPTS-68 but we now use ITS-90
            as.double(pressure),
            value = double(nC),
            NAOK = TRUE, PACKAGE = "oce"
        )$value
    } else if (eos == "gsw") {
        # we don't need to convert to IPTS-68 for the gsw formulation, because it is already formulated
        # to work with ITS-90
        #> message("swSCTp() gsw; conductivity[1]=", conductivity[1], ", temperature[1]=", temperature[1], ", pressure[1]=", pressure[1])
        res <- gsw::gsw_SP_from_C(C0 * conductivity, temperature, pressure)
    }
    dim(res) <- dim
    res
}


#' Seawater Salinity From Temperature and Density
#'
#' Compute Practical or Absolute Salinity, given in-situ or Conservative
#' Temperature, density, and pressure. This is mainly used to draw isopycnal
#' lines on TS diagrams, hence the dual meanings for salinity and temperature,
#' depending on the value of `eos`.
#'
#' For `eos="unesco"`, finds the practical salinity that yields the given
#' density, with the given in-situ temperature and pressure.  The method is a
#' bisection search with a salinity tolerance of 0.001.  For `eos="gsw"`,
#' the function [gsw::gsw_SA_from_rho()] in the \CRANpkg{gsw}
#' package is used
#' to infer Absolute Salinity from Conservative Temperature.
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param density *in-situ* density or sigma value (\eqn{kg/m^3}{kg/m^3})
#'
#' @param pressure *in-situ* pressure (dbar)
#'
#' @param eos equation of state, either `"unesco"` (see references 1 and 2)
#' or `"gsw"` (see references 3 and 4).
#'
#' @return Practical Salinity, if `eos="unesco"`, or Absolute Salinity, if
#' `eos="gsw"`.
#'
#' @author Dan Kelley
#'
#' @seealso [swTSrho()]
#'
#' @references
#' 1. Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for computation of
#' fundamental properties of seawater. *Unesco Technical Papers in Marine
#' Science*, *44*, 53 pp
#'
#' 2. Gill, A.E., 1982. *Atmosphere-ocean Dynamics*, Academic Press, New
#' York, 662 pp.
#'
#' 3. IOC, SCOR, and IAPSO (2010). The international thermodynamic equation of
#' seawater-2010: Calculation and use of thermodynamic properties.  Technical
#' Report 56, Intergovernmental Oceanographic Commission, Manuals and Guide.
#'
#' 4. McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127,
#' ISBN 978-0-646-55621-5.
#'
#' @examples
#' swSTrho(10, 22, 0, eos = "gsw") # 28.76285
#' swSTrho(10, 22, 0, eos = "unesco") # 28.651625
#'
#' @family functions that calculate seawater properties
swSTrho <- function(temperature, density, pressure, eos = getOption("oceEOS", default = "gsw")) {
    # cat(file=stderr(), "swSTrho(", temperature[1], ", ", density[1], ", ", pressure[1], ", ", eos, "\n", sep="")
    # cat(file=stderr(), vectorShow(length(temperature)))
    # FIXME-gsw add gsw version
    eos <- match.arg(eos, c("unesco", "gsw"))
    teos <- eos == "gsw" # FIXME still the best way?
    dim <- dim(temperature)
    nt <- length(temperature)
    nrho <- length(density)
    np <- length(pressure)
    if (nrho == 1) density <- rep(density, nt)
    if (np == 1) pressure <- rep(pressure, nt)
    if (nt == 1) temperature <- rep(temperature, nt)
    sigma <- ifelse(density > 500, density - 1000, density)
    if (eos == "unesco") {
        # cat(file=stderr(), "  about to call C function sw_strho()\n")
        res <- .C("sw_strho",
            as.integer(nt),
            as.double(T68fromT90(temperature)),
            as.double(sigma),
            as.double(pressure),
            as.integer(teos),
            S = double(nt),
            NAOK = TRUE, PACKAGE = "oce"
        )$S
        # NAOK=TRUE)$S # permits dyn.load() on changing .so
    } else if (eos == "gsw") {
        density <- ifelse(density < 900, density + 1000, density)
        res <- gsw::gsw_SA_from_rho(density, temperature, pressure) # assumes temperature=CT
    }
    dim(res) <- dim
    res
}


#' Seawater Temperature from Salinity and Density
#'
#' Compute *in-situ* temperature, given salinity, density, and pressure.
#'
#' Finds the temperature that yields the given density, with the given salinity
#' and pressure.  The method is a bisection search with temperature tolerance
#' 0.001 \eqn{^\circ}{deg}C.
#'
#' @param salinity *in-situ* salinity (PSU)
#'
#' @param density *in-situ* density or sigma value (kg/m\eqn{^3}{^3})
#'
#' @param pressure *in-situ* pressure (dbar)
#'
#' @param eos equation of state to be used, either `"unesco"` or
#' `"gsw"` (ignored at present).
#'
#' @return *In-situ* temperature in \eqn{^\circ}{deg}C on the ITS-90
#' scale.
#'
#' @author Dan Kelley
#'
#' @seealso [swSTrho()]
#'
#' @references Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for
#' computation of fundamental properties of seawater. *Unesco Technical
#' Papers in Marine Science*, *44*, 53 pp
#'
#' Gill, A.E., 1982. *Atmosphere-ocean Dynamics*, Academic Press, New
#' York, 662 pp.
#'
#' @examples
#' swTSrho(35, 23, 0, eos = "unesco") # 26.11301
#'
#' @family functions that calculate seawater properties
swTSrho <- function(salinity, density, pressure = NULL, eos = getOption("oceEOS", default = "gsw")) {
    # cat(file=stderr(), "swTSrho(", salinity[1],", ", density[1], ", ", pressure[1], ", ", eos, "\n")
    # cat(file=stderr(), vectorShow(length(salinity)))
    # FIXME-gsw add gsw version
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    eos <- match.arg(eos, c("unesco", "gsw"))
    teos <- eos == "gsw"
    dim <- dim(salinity)
    nS <- length(salinity)
    nrho <- length(density)
    if (is.null(pressure)) {
        pressure <- rep(0, nS)
    }
    if (length(pressure) == 1) {
        pressure <- rep(pressure[1], length.out = nS)
    }
    np <- length(pressure)
    if (nS != nrho) {
        stop("lengths of salinity and rho must agree, but they are ", nS, " and ", nrho, ", respectively")
    }
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    for (i in 1:nS) {
        # FIXME: avoid loops
        sig <- density[i]
        if (sig > 500) {
            sig <- sig - 1000
        }
        # FIXME: is this right for all equations of state? I doubt it
        this.T <- .C("sw_tsrho",
            as.double(salinity[i]),
            as.double(sig),
            as.double(pressure[i]),
            as.integer(teos),
            temperature = double(1),
            NAOK = TRUE, PACKAGE = "oce"
        )$t
        this.T <- T90fromT68(this.T)
        if (i == 1) res <- this.T else res <- c(res, this.T)
    }
    dim(res) <- dim
    res
}


#' Seawater Freezing Temperature
#'
#' Compute in-situ freezing temperature of seawater, using either the UNESCO
#' formulation (computed as in Section 5 of Fofonoff and Millard, 1983) or the
#' GSW formulation (computed by using [gsw::gsw_SA_from_SP()] to get Absolute
#' Salinity, and then [gsw::gsw_t_freezing()] to get the freezing temperature).
#'
#' If the first argument is an `oce` object, and if the `pressure` argument is
#' `NULL`, then the pressure is sought within the first argument. In the case of
#' `eos="gsw"`, then a similar procedure also applies to the `longitude` and
#' `latitude` arguments.
#'
#' @param salinity Either practical salinity (PSU) or a `ctd` object from which
#' practical salinity and pressure (plus in the `eos="gsw"` case,
#' longitude and latitude) are inferred.
#'
#' @param pressure Seawater pressure (dbar).
#'
#' @param longitude Longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#'
#' @param latitude Latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#'
#' @param saturation_fraction The saturation fraction of dissolved air in seawater,
#' ignored if `eos="unesco"`).
#'
#' @param eos The equation of state, either `"unesco"` (Fofonoff and Millard,
#' 1983; Gill 1982) or `"gsw"` (IOC, SCOR and IAPSO 2010; McDougall and Barker
#' 2011).
#'
#' @return Temperature (degC), defined on the ITS-90 scale.
#'
#' @author Dan Kelley
#'
#' @references
#'
#' Fofonoff, N. P., and R. C. Millard. Algorithms for Computation of
#' Fundamental Properties of Seawater. UNESCO Technical Papers in Marine
#' Research. SCOR working group on Evaluation of CTD data; UNESCO/ICES/SCOR/IAPSO
#' Joint Panel on Oceanographic Tables and Standards, 1983.
#'
#' Gill, A E. Atmosphere-Ocean Dynamics. New York, NY, USA: Academic Press,
#' 1982.
#'
#' IOC, SCOR, and IAPSO (2010). The international thermodynamic equation of
#' seawater-2010: Calculation and use of thermodynamic properties.  Technical
#' Report 56, Intergovernmental Oceanographic Commission, Manuals and Guide, 2010.
#'
#' McDougall, Trevor J., and Paul M. Barker. Getting Started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox. SCOR/IAPSO WG127, 2011.
#'
#' @examples
#' # 1. Test for a check-value given in reference 1. This value, -2.588567 degC,
#' # is in the 1968 temperature scale (IPTS-68), but swTFreeze reports
#' # in the newer ITS-90 scale, so we must convert before checking.
#' Tcheck <- -2.588567 # IPTS-68
#' T <- swTFreeze(salinity = 40, pressure = 500, eos = "unesco")
#' stopifnot(abs(Tcheck - T68fromT90(T)) < 1e-6)
#'
#' # 2. Compare unesco and gsw formulations.
#' data(ctd)
#' p <- ctd[["pressure"]]
#' par(mfrow = c(1, 2), mar = c(3, 3, 1, 2), mgp = c(2, 0.7, 0))
#' plot(swTFreeze(ctd, eos = "unesco"),
#'     p,
#'     xlab = "unesco", ylim = rev(range(p))
#' )
#' plot(swTFreeze(ctd, eos = "unesco") - swTFreeze(ctd, eos = "gsw"),
#'     p,
#'     xlab = "unesco-gsw", ylim = rev(range(p))
#' )
#'
#' @family functions that calculate seawater properties
swTFreeze <- function(
    salinity, pressure = NULL, longitude = NULL, latitude = NULL,
    saturation_fraction = 1, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must supply salinity (which may be S or a CTD object)")
    }
    if (inherits(salinity, "oce")) {
        if (is.null(pressure)) {
            pressure <- salinity[["pressure"]]
        }
    }
    if (is.null(pressure)) {
        stop("must supply pressure")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        # Note: the pressure in the next line is for computing SA; see below.
        l <- lookWithin(list(salinity = salinity, latitude = latitude, longitude = longitude, pressure = pressure))
    } else {
        l <- lookWithin(list(salinity = salinity, pressure = pressure))
    }
    dim <- dim(l$salinity)
    if (eos == "gsw") {
        # Note that l$pressure is used for computing SA, but not for gsw_t_freezing().
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        res <- gsw::gsw_t_freezing(SA = SA, p = l$pressure, saturation_fraction = saturation_fraction)
    } else if (eos == "unesco") {
        res <- (-.0575 + 1.710523e-3 * sqrt(abs(l$salinity)) - 2.154996e-4 * l$salinity) * l$salinity - 7.53e-4 * l$pressure
        res <- T90fromT68(res)
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Thermal Expansion Coefficient
#'
#' Compute \eqn{\alpha}{alpha}, the thermal expansion coefficient for seawater.
#'
#' @param salinity either practical salinity (in which case `temperature`
#' and `pressure` must be provided) *or* an `oce` object (in
#' which case `salinity`, etc. are inferred from the object).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure pressure (dbar)
#'
#' @param longitude longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#'
#' @param latitude latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#'
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#'
#' @return Value in 1/degC.
#'
#' @author Dan Kelley
#'
#' @references The `eos="unesco"` formulae are based on the UNESCO
#' equation of state, but are formulated empirically by Trevor J. McDougall,
#' 1987, Neutral Surfaces, Journal of Physical Oceanography, volume 17, pages
#' 1950-1964. The `eos="gsw"` formulae come from GSW; see references in
#' the [swRho()] documentation.
#'
#' @family functions that calculate seawater properties
swAlpha <- function(
    salinity, temperature = NULL, pressure = 0,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
    }
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    if (length(l$pressure) == 1) l$pressure <- rep(l$pressure, length.out = nS)
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    if (l$eos == "unesco") {
        res <- swAlphaOverBeta(l$salinity, l$temperature, l$pressure, eos = "unesco") * swBeta(l$salinity, l$temperature, l$pressure, eos = "unesco")
    } else if (l$eos == "gsw") {
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        CT <- gsw::gsw_CT_from_t(SA = SA, t = l$temperature, p = l$pressure)
        res <- gsw::gsw_alpha(SA = SA, CT = CT, p = l$pressure)
    }
    res
}


#' Ratio of Seawater Thermal Expansion Coefficient to Haline Contraction Coefficient
#'
#' Compute \eqn{\alpha/\beta}{alpha/beta} using McDougall's (1987) algorithm.
#'
#' @param salinity either practical salinity (in which case `temperature`
#' and `pressure` must be provided) *or* an `oce` object (in
#' which case `salinity`, etc. are inferred from the object).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C)
#'
#' @param pressure pressure (dbar)
#'
#' @param longitude longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#'
#' @param latitude latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#'
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#'
#' @return Value in psu/\eqn{^\circ}{deg}C.
#'
#' @author Dan Kelley
#'
#' @references The `eos="unesco"` formulae are based on the UNESCO
#' equation of state, but are formulated empirically by Trevor J. McDougall,
#' 1987, Neutral Surfaces, Journal of Physical Oceanography, volume 17, pages
#' 1950-1964. The `eos="gsw"` formulae come from GSW; see references in
#' the [swRho()] documentation.
#'
#' @examples
#' swAlphaOverBeta(40, 10, 4000, eos = "unesco") # 0.3476
#'
#' @family functions that calculate seawater properties
swAlphaOverBeta <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            temperature <- salinity[["temperature"]]
            pressure <- salinity[["pressure"]]
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
    }
    dim <- dim(l$salinity)
    if (is.null(l$temperature)) {
        stop("must provide temperature")
    }
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    if (is.null(l$pressure)) pressure <- 0
    if (length(l$pressure) != nS) l$pressure <- rep(l$pressure, length.out = nS)
    if (l$eos == "gsw") {
        # not likely to be called since gsw has a direct function for alpha, but put this here anyway
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        CT <- gsw::gsw_CT_from_t(SA = SA, t = l$temperature, p = l$pressure)
        res <- gsw::gsw_alpha_on_beta(SA = SA, CT = CT, p = l$pressure)
    } else if (l$eos == "unesco") {
        theta <- swTheta(l$salinity, l$temperature, l$pressure, eos = "unesco")
        res <- .C("sw_alpha_over_beta", as.integer(nS),
            as.double(l$salinity), as.double(theta), as.double(l$pressure),
            value = double(nS), NAOK = TRUE, PACKAGE = "oce"
        )$value
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Haline Contraction Coefficient
#'
#' Compute \eqn{\beta}{beta}, the haline contraction coefficient for seawater.
#'
#' @param salinity either practical salinity (in which case `temperature`
#' and `pressure` must be provided) *or* an `oce` object (in
#' which case `salinity`, etc. are inferred from the object).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure seawater pressure (dbar)
#'
#' @param longitude longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#'
#' @param latitude latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#'
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#'
#' @return Value in 1/psu.
#'
#' @author Dan Kelley
#'
#' @references The `eos="unesco"` formulae are based on the UNESCO
#' equation of state, but are formulated empirically by Trevor J. McDougall,
#' 1987, Neutral Surfaces, Journal of Physical Oceanography, volume 17, pages
#' 1950-1964. The `eos="gsw"` formulae come from GSW; see references in
#' the [swRho()] documentation.
#'
#' @family functions that calculate seawater properties
swBeta <- function(
    salinity, temperature = NULL, pressure = 0,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        temperature <- salinity[["temperature"]]
        pressure <- salinity[["pressure"]]
        if (is.null(longitude)) {
            longitude <- salinity[["longitude"]]
        }
        if (is.null(latitude)) {
            latitude <- salinity[["latitude"]]
        }
    }
    if (eos == "gsw") {
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
    }
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    if (length(l$pressure) == 1) l$pressure <- rep(l$pressure, length.out = nS)
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    if (eos == "gsw") {
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        CT <- gsw::gsw_CT_from_t(SA = SA, t = l$temperature, p = l$pressure)
        res <- gsw::gsw_beta(SA = SA, CT = CT, p = l$pressure)
    } else if (eos == "unesco") {
        theta <- swTheta(l$salinity, l$temperature, l$pressure, eos = "unesco") # the formula is i.t.o. theta
        res <- .C("sw_beta", as.integer(nS),
            as.double(l$salinity), as.double(theta), as.double(l$pressure),
            value = double(nS), NAOK = TRUE, PACKAGE = "oce"
        )$value
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}

# thermal (not electrical) conductivity, using Caldwell (1974) as of 2015-jan-9
# NOTE: no gsw equivalent


#' Seawater Thermal Conductivity
#'
#' Compute seawater thermal conductivity, in \eqn{W
#' m^{-1\circ}C^{-1}}{W/(m*degC)}
#'
#' Caldwell's (1974) detailed formulation is used.  To be specific, his
#' equation 6 to calculate K, and his two sentences above that equation are
#' used to infer this to be K(0,T,S) in his notation of equation 7. Then,
#' application of his equations 7 and 8 is straightforward. He states an
#' accuracy for this method of 0.3 percent.  (See the check against his Table 1
#' in the \dQuote{Examples}.)
#'
#' @param salinity salinity (PSU), or a `ctd` object, in which case
#' `temperature` and `pressure` will be ignored.
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure pressure (dbar)
#'
#' @return Conductivity of seawater in \eqn{W m^{-1} {^\circ} C^{-1}}{W/(m*degC)}.
#' To calculate thermal diffusivity in \eqn{m^2/s^2}, divide by the
#' product of density and specific heat, as in the example.
#'
#' @author Dan Kelley
#'
#' @references Caldwell, Douglas R., 1974. Thermal conductivity of seawater,
#' *Deep-sea Research*, *21*, 131-137.
#'
#' @examples
#' library(oce)
#' # Values in m^2/s, a unit that is often used instead of W/(m*degC).
#' swThermalConductivity(35, 10, 100) / (swRho(35, 10, 100) * swSpecificHeat(35, 10, 100)) # ocean
#' swThermalConductivity(0, 20, 0) / (swRho(0, 20, 0) * swSpecificHeat(0, 20, 0)) # lab
#' # Caldwell Table 1 gives 1478e-6 cal/(cm*sec*degC) at 31.5 o/oo, 10degC, 1kbar
#' joulePerCalorie <- 4.18400
#' cmPerM <- 100
#' swThermalConductivity(31.5, 10, 1000) / joulePerCalorie / cmPerM
#'
#' @family functions that calculate seawater properties
swThermalConductivity <- function(salinity, temperature = NULL, pressure = NULL) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure))
    # below is formula used prior to 2015-jan-9
    # return(0.57057 * (1 + l$temperature * (0.003 - 1.025e-05 * l$temperature) + 0.000653 * l$pressure - 0.00029 * l$salinity))
    S <- l$salinity
    # nolint start T_and_F_symbol_linter
    T <- T68fromT90(l$temperature)
    p <- l$pressure / 1e3 # Caldwell formula is for kbar, not dbar
    if (TRUE) {
        K <- 1.3507e-3 + 4.061e-6 * T + (-1.412e-8) * T^2 # Caldwell eq 6
        f <- 0.0690 + (-8e-5) * T + (-0.0020) * p + (-0.00010) * S # Caldwell eq 8
        cond <- K * (1 + f) # Caldwell eq 7
    } else {
        cond <- 0.001365 * (1 + 0.003 * T - 1.025e-5 * T^2 + 0.0653 * p - 0.00029 * S)
    }
    # nolint end T_and_F_symbol_linter
    418.400 * cond # convert from cal/(cm*sec*degC) to J/(m*sec*degC)
}


#' Water Depth
#'
#' Retrieve or compute depth below the surface, i.e. a positive number within
#' the water column.  If the first parameter is an oce object that has an
#' element named `"depth"` in its `data` slot, then return that value.
#' Otherwise, compute depth from a formula that includes pressure and latitude,
#' as explained in \sQuote{Details}.
#'
#' For the calculated case, the method depends on the value of `eos` parameter.
#' If this is `"unesco"`, then depth is calculated from pressure using Saunders
#' and Fofonoff's method, with the formula refitted for 1980 UNESCO equation of
#' state (reference 1).  On the other hand, if it is `eos="gsw"`, then
#' [gsw::gsw_z_from_p()] from the \CRANpkg{gsw} package (references 2 and 3) is
#' used.
#'
#' @param pressure either pressure (dbar), in which case `latitude` must also
#' be given, or a `ctd` object, in which case `latitude` will be inferred
#' from the object.
#'
#' @param latitude numeric value for latitude in degrees North.
#'
#' @param eos character value indicating the formulation to be used, either
#' `"unesco"` or `"gsw"`.
#'
#' @template debugTemplate
#'
#' @return [swDepth] returns depth below the ocean surface, in metres.
#'
#' @author Dan Kelley
#'
#' @references
#' 1. Unesco 1983. Algorithms for computation of fundamental
#' properties of seawater, 1983. *Unesco Tech. Pap. in Mar. Sci.*, No. 44,
#' 53 pp.
#'
#' 2. IOC, SCOR, and IAPSO (2010). The international thermodynamic equation of
#' seawater-2010: Calculation and use of thermodynamic properties.  Technical
#' Report 56, Intergovernmental Oceanographic Commission, Manuals and Guide.
#'
#' 3. McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127,
#' ISBN 978-0-646-55621-5.
#'
#' @examples
#' d <- swDepth(10, 45)
#'
#' @family functions that calculate seawater properties
swDepth <- function(pressure, latitude = 45, eos = getOption("oceEOS", default = "gsw"),
                    debug = getOption("oceDebug")) {
    # FIXME-gsw need a gsw version but it is not in the C library as of Dec 2014
    oceDebug(debug, "swDepth(...) ...\n", sep = "", unindent = 1)
    if (missing(pressure)) {
        stop("must provide pressure")
    }
    if (inherits(pressure, "oce")) {
        x <- pressure # store this for clarity
        if ("depth" %in% names(x@data)) {
            return(x@data$depth)
        }
        oceDebug(debug, "swDepth() with an oce object as first argument ...\n")
        location <- locationForGsw(x)
        if (!is.null(location$latitude)) {
            latitude <- location$latitude
        }
        oceDebug(debug, vectorShow(latitude))
    }
    dim <- dim(pressure)
    l <- lookWithin(list(pressure = pressure, eos = eos))
    l$latitude <- latitude
    if (any(is.na(l$latitude))) {
        l$latitude <- 45
    } # default to mid latitudes
    if (l$eos == "unesco") {
        l$latitude <- l$latitude * atan2(1, 1) / 45
        x <- sin(l$latitude)^2
        gr <- 9.780318 * (1.0 + (5.2788e-3 + 2.36e-5 * x) * x) + 1.092e-6 * l$pressure
        res <- (((-1.82e-15 * l$pressure + 2.279e-10) * l$pressure - 2.2512e-5) * l$pressure + 9.72659) * l$pressure / gr
    } else if (l$eos == "gsw") {
        res <- -gsw::gsw_z_from_p(
            p = as.vector(l$pressure),
            latitude = as.vector(l$latitude)
        )
        dim(res) <- dim(l$pressure)
    }
    oceDebug(debug, "... # swDepth()\n", sep = "", unindent = 1)
    res
}

#' Vertical Coordinate
#'
#' Compute height above the surface. This is the negative of depth,
#' and so is defined simply in terms of [swDepth()].
#'
#' @inheritParams swDepth
#'
#' @family functions that calculate seawater properties
swZ <- function(pressure, latitude = 45, eos = getOption("oceEOS", default = "gsw"),
                debug = getOption("oceDebug")) {
    # FIXME-gsw need a gsw version but it is not in the C library as of Dec 2014
    if (missing(pressure)) {
        stop("must provide pressure")
    }
    -swDepth(pressure = pressure, latitude = latitude, eos = eos)
}


#' Dynamic Height of a Seawater Profile
#'
#' Compute the dynamic height of a column of seawater.
#'
#' If the first argument is a `section`, then dynamic height is calculated
#' for each station within a section, and returns a list containing distance
#' along the section along with dynamic height.
#'
#' If the first argument is a `ctd`, then this returns just a single
#' value, the dynamic height.
#'
#' If `eos="unesco"`, processing is as follows.  First, a piecewise-linear
#' model of the density variation with pressure is developed using
#' [stats::approxfun()]. (The option `rule=2` is used to
#' extrapolate the uppermost density up to the surface, preventing a possible a
#' bias for bottle data, in which the first depth may be a few metres below the
#' surface.) A second function is constructed as the density of water with
#' salinity 35PSU, temperature of 0\eqn{^\circ}{deg}C, and pressure as in the
#' `ctd`. The difference of the reciprocals of these densities, is then
#' integrated with [stats::integrate()] with pressure limits `0`
#' to `referencePressure`.  (For improved numerical results, the variables
#' are scaled before the integration, making both independent and dependent
#' variables be of order one.)
#'
#' If `eos="gsw"`, [gsw::gsw_geo_strf_dyn_height()] is used
#' to calculate a result in m^2/s^2, and this is divided by
#' 9.7963\eqn{m/s^2}{m/s^2}.
#' If pressures are out of order, the data are sorted. If any pressure
#' is repeated, only the first level is used.
#' If there are under 4 remaining distinct
#' pressures, `NA` is returned, with a warning.
#'
#' @param x a [section-class] object.
#'
#' @param referencePressure reference pressure (dbar). If this exceeds the
#' highest pressure supplied to [swDynamicHeight()], then that highest
#' pressure is used, instead of the supplied value of
#' `referencePressure`.
#'
#' @param subdivisions number of subdivisions for call to
#' [integrate()].  (The default value is considerably larger than the
#' default for [integrate()], because otherwise some test profiles
#' failed to integrate.
#'
#' @param rel.tol absolute tolerance for call to [integrate()].  Note
#' that this call is made in scaled coordinates, i.e. pressure is divided by
#' its maximum value, and dz/dp is also divided by its maximum.
#'
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#'
#' @return In the first form, a list containing `distance`, the distance
#' (km( from the first station in the section and `height`, the dynamic
#' height (m). In the second form, a single value, containing the
#' dynamic height (m).
#'
#' @author Dan Kelley
#'
#' @references Gill, A.E., 1982. *Atmosphere-ocean Dynamics*, Academic
#' Press, New York, 662 pp.
#'
#' @section Sample of Usage:
#' \preformatted{
#' library(oce)
#' data(section)
#'
#' # Dynamic height and geostrophy
#' par(mfcol=c(2, 2))
#' par(mar=c(4.5, 4.5, 2, 1))
#'
#' # Left-hand column: whole section
#' # (The smoothing lowers Gulf Stream speed greatly)
#' westToEast <- subset(section, 1<=stationId&stationId<=123)
#' dh <- swDynamicHeight(westToEast)
#' plot(dh$distance, dh$height, type="p", xlab="", ylab="dyn. height [m]")
#' ok <- !is.na(dh$height)
#' smu <- supsmu(dh$distance, dh$height)
#' lines(smu, col="blue")
#' f <- coriolis(section[["station", 1]][["latitude"]])
#' g <- gravity(section[["station", 1]][["latitude"]])
#' v <- diff(smu$y)/diff(smu$x) * g / f / 1e3 # 1e3 converts to m
#' plot(smu$x[-1], v, type="l", col="blue", xlab="distance [km]", ylab="velocity (m/s)")
#'
#' # right-hand column: gulf stream region, unsmoothed
#' gs <- subset(section, 102<=stationId&stationId<=124)
#' dh.gs <- swDynamicHeight(gs)
#' plot(dh.gs$distance, dh.gs$height, type="b", xlab="", ylab="dyn. height [m]")
#' v <- diff(dh.gs$height)/diff(dh.gs$distance) * g / f / 1e3
#' plot(dh.gs$distance[-1], v, type="l", col="blue",
#'   xlab="distance [km]", ylab="velocity (m/s)")
#' }
#'
#' @family functions that calculate seawater properties
swDynamicHeight <- function(
    x, referencePressure = 2000,
    subdivisions = 500, rel.tol = .Machine$double.eps^0.25,
    eos = getOption("oceEOS", default = "gsw")) {
    eos <- match.arg(eos, c("unesco", "gsw"))
    height <- function(ctd, referencePressure, subdivisions, rel.tol, eos = getOption("oceEOS", default = "gsw")) {
        if (sum(!is.na(ctd@data$pressure)) < 2) {
            return(NA)
        }
        g <- if (is.na(ctd@metadata$latitude)) 9.8 else gravity(ctd@metadata$latitude)
        p <- ctd[["pressure"]]
        np <- length(p)
        p_ref <- min(max(p, na.rm = TRUE), referencePressure)
        if (eos == "unesco") {
            rho <- swRho(ctd, eos = eos)
            if (sum(!is.na(rho)) < 2) {
                return(NA)
            }
            # 1e4 converts decibar to Pa
            dzdp <- ((1 / rho - 1 / swRho(rep(35, np), rep(0, np), p, eos = eos)) / g) * 1e4
            # Scale both pressure and dz/dp to make integration work better (issue 499)
            max <- max(dzdp, na.rm = TRUE)
            integrand <- approxfun(p / p_ref, dzdp / max, rule = 2)
            # plot(dzdp/max, ctd@data$pressure/referencePressure, type="l")
            res <- integrate(integrand, 0, 1,
                subdivisions = subdivisions, rel.tol = rel.tol
            )$value * p_ref * max
        } else { # "gsw"
            if (np > 3) {
                o <- order(p)
                p <- p[o]
                SA <- ctd[["SA"]][o]
                CT <- ctd[["CT"]][o]
                # handle repeated values
                dp0 <- diff(p) == 0
                if (any(dp0)) {
                    SA <- SA[!dp0]
                    CT <- CT[!dp0]
                    p <- p[!dp0]
                    g <- 9.7963 # NOTE: do not use local gravity; the gsw is defined to use this value
                    if (length(p) < 4) {
                        warning("In swDynamicHeight() : returning NA since < 4 levels", call. = FALSE)
                        res <- NA
                    } else {
                        res <- gsw::gsw_geo_strf_dyn_height(SA = SA, CT = CT, p = p, p_ref = p_ref)[1] / g
                    }
                } else {
                    res <- gsw::gsw_geo_strf_dyn_height(SA = SA, CT = CT, p = p, p_ref = p_ref)[1] / g
                }
                res[is.nan(res)] <- NA
            } else {
                warning("In swDynamicHeight() : returning NA since < 4 levels", call. = FALSE)
                res <- NA
            }
        }
        res
    }
    if (inherits(x, "section")) {
        lon0 <- x@data$station[[1]]@metadata$longitude
        lat0 <- x@data$station[[1]]@metadata$latitude
        ns <- length(x@data$station)
        d <- vector("numeric", ns)
        h <- vector("numeric", ns)
        for (i in 1:ns) {
            d[i] <- geodDist(x@data$station[[i]]@metadata$longitude, x@data$station[[i]]@metadata$latitude, lon0, lat0)
            h[i] <- height(x@data$station[[i]], referencePressure, subdivisions = subdivisions, rel.tol = rel.tol, eos = eos)
        }
        return(list(distance = d, height = h))
    } else if (inherits(x, "ctd")) {
        return(height(x, referencePressure, subdivisions = subdivisions, rel.tol = rel.tol, eos = eos))
    } else {
        stop("method is only for objects of class '", "section", " or '", "ctd", "'")
    }
}


#' Seawater Lapse Rate
#'
#' Compute adiabatic lapse rate
#'
#' If `eos="unesco"`, the density is calculated using the UNESCO equation
#' of state for seawater (references 1 and 2), and if `eos="gsw"`, the GSW formulation
#' (references 3 and 4) is used.
#'
#' @param salinity either salinity (PSU) (in which case `temperature` and
#' `pressure` must be provided) *or* a `ctd` object (in which
#' case `salinity`, `temperature` and `pressure` are determined
#' from the object, and must not be provided in the argument list).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure pressure (dbar)
#'
#' @param longitude longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#'
#' @param latitude latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#'
#' @param eos equation of state, either `"unesco"` (references 1 and 2)
#' or `"gsw"` (references 3 and 4).
#'
#' @return Lapse rate (\eqn{deg}{deg}C/m).
#'
#' @author Dan Kelley
#'
#' @references Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for
#' computation of fundamental properties of seawater.
#' *Unesco Technical Papers in Marine Science*,
#' *44*, 53 pp.  (Section 7, pages 38-40)
#'
#' @examples
#' lr <- swLapseRate(40, 40, 10000) # 3.255976e-4
#'
#' @family functions that calculate seawater properties
swLapseRate <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
    }
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    if (is.null(l$temperature)) {
        stop("must provide temperature")
    }
    nt <- length(l$temperature)
    if (is.null(l$pressure)) {
        l$pressure <- rep(0, length.out = nS)
    }
    if (length(l$pressure) == 1) {
        l$pressure <- rep(l$pressure[1], length.out = nS)
    }
    np <- length(l$pressure)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    if (eos == "unesco") {
        res <- .C("sw_lapserate", as.integer(nS), as.double(l$salinity), as.double(T68fromT90(l$temperature)), as.double(l$pressure),
            value = double(nS), NAOK = TRUE, PACKAGE = "oce"
        )$value
    } else if (eos == "gsw") {
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        CT <- gsw::gsw_CT_from_t(SA = SA, t = l$temperature, p = l$pressure)
        # the 1e4 is to convert from 1/Pa to 1/dbar
        res <- 1e4 * gsw::gsw_adiabatic_lapse_rate_from_CT(SA = SA, CT = CT, p = l$pressure)
    } else {
        stop("eos must be either \"unesco\" or \"eos\"")
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
} # swLapseRate


#' Seawater Density
#'
#' Compute \eqn{\rho}{rho}, the *in-situ* density of seawater.
#'
#' If `eos="unesco"`, the density is calculated using the UNESCO equation
#' of state for seawater (references 1 and 2), and if `eos="gsw"`, the GSW formulation
#' (references 3 and 4) is used.
#'
#' @param salinity either practical salinity (in which case `temperature`
#' and `pressure` must be provided) *or* an `oce` object, in
#' which case `salinity`, `temperature` (in the ITS-90 scale; see
#' next item), etc. are inferred from the object, ignoring the
#' other parameters, if they are supplied.
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale.  This scale is used by GSW-style calculation (as
#' requested by setting `eos="gsw"`), and is the value contained within
#' `ctd` objects (and probably most other objects created with data
#' acquired in the past decade or two). Since the UNESCO-style calculation is
#' based on IPTS-68, the temperature is converted within the present function,
#' using [T68fromT90()].
#'
#' @param pressure pressure (dbar)
#'
#' @param longitude longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#'
#' @param latitude latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#'
#' @param eos equation of state, either `"unesco"` (references 1 and 2)
#' or `"gsw"` (references 3 and 4).
#'
#' @template debugTemplate
#'
#' @return *In-situ* density (kg/m\eqn{^3}{^3}).
#'
#' @section Temperature units: The UNESCO formulae are defined in terms of
#' temperature measured on the IPTS-68 scale, whereas the replacement GSW
#' formulae are based on the ITS-90 scale. Prior to the addition of GSW
#' capabilities, the various `sw*` functions took temperature to be in
#' IPTS-68 units. As GSW capabilities were added in early 2015, the assumed
#' unit of `temperature` was taken to be ITS-90.  This change means that
#' old code has to be modified, by replacing e.g. `swRho(S, T, p)` with
#' `swRho(S, T90fromT68(T), p)`. At typical oceanic values, the difference
#' between the two scales is a few millidegrees.
#'
#' @author Dan Kelley
#'
#' @seealso Related density routines include [swSigma0()] (and
#' equivalents at other pressure horizons), [swSigmaT()], and
#' [swSigmaTheta()].
#'
#' @references
#' 1. Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for computation of
#' fundamental properties of seawater.
#' *Unesco Technical Papers in Marine Science*,
#' *44*, 53 pp.
#'
#' 2. Gill, A.E., 1982. *Atmosphere-ocean Dynamics*, Academic Press, New
#' York, 662 pp.
#'
#' 3. IOC, SCOR, and IAPSO (2010). The international thermodynamic equation of
#' seawater-2010: Calculation and use of thermodynamic properties.  Technical
#' Report 56, Intergovernmental Oceanographic Commission, Manuals and Guide.
#'
#' 4. McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127,
#' ISBN 978-0-646-55621-5.
#'
#' @examples
#' library(oce)
#' # The numbers in the comments are the check values listed in reference 1;
#' # note that temperature in that reference was on the T68 scale, but that
#' # the present function works with the ITS-90 scale, so a conversion
#' # is required.
#' swRho(35, T90fromT68(5), 0, eos = "unesco") # 1027.67547
#' swRho(35, T90fromT68(5), 10000, eos = "unesco") # 1069.48914
#' swRho(35, T90fromT68(25), 0, eos = "unesco") # 1023.34306
#' swRho(35, T90fromT68(25), 10000, eos = "unesco") # 1062.53817
#'
#' @family functions that calculate seawater properties
swRho <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw"),
    debug = getOption("oceDebug")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        if (eos == "gsw") { # do not need these for UNESCO calculations
            longitude <- salinity[["longitude"]]
            latitude <- salinity[["latitude"]]
        }
        temperature <- salinity[["temperature"]]
        pressure <- salinity[["pressure"]]
        salinity <- salinity[["salinity"]]
    }
    if (eos == "gsw") {
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
    } else { # must be "unesco"
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
    }
    if (is.null(l$temperature)) {
        stop("must provide temperature")
    }
    if (is.null(l$pressure)) {
        stop("must provide pressure")
    }
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    if (is.null(l$pressure)) l$pressure <- 0
    if (length(l$pressure) == 1) l$pressure <- rep(l$pressure, length.out = nS)
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    if (eos == "unesco") {
        res <- .C("sw_rho", as.integer(nS), as.double(l$salinity),
            as.double(T68fromT90(l$temperature)),
            as.double(l$pressure),
            value = double(nS), NAOK = TRUE, PACKAGE = "oce"
        )$value
    } else if (eos == "gsw") {
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        CT <- gsw::gsw_CT_from_t(SA = SA, t = l$temperature, p = l$pressure)
        res <- gsw::gsw_rho(SA, CT, p = l$pressure)
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Density Anomaly
#'
#' Compute \eqn{\sigma_\theta}{sigma}, the density of seawater, minus 1000
#' kg/m\eqn{^3}{^3}.
#'
#' @inheritParams swRho
#'
#' @return Density anomaly (kg/m\eqn{^3}{^3}), as computed with [swRho()], minus-
#' 1000 kg/m\eqn{^3}{^3}.
#'
#' @author Dan Kelley
#'
#' @references See citations provided in the [swRho()] documentation.
#'
#' @examples
#' library(oce)
#' swSigma(35, 13, 1000, longitude = 300, latitude = 30, eos = "gsw") # 30.82374
#' swSigma(35, T90fromT68(13), 1000, eos = "unesco") # 30.8183
#'
#' @family functions that calculate seawater properties
swSigma <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    swRho(salinity, temperature, pressure,
        longitude = longitude, latitude = latitude, eos = eos
    ) - 1000
}


#' Seawater Quasi-Potential Density Anomaly
#'
#' Compute \eqn{\sigma_t}{sigma-t}, a rough estimate of potential density of
#' seawater, minus 1000 kg/m\eqn{^3}{^3}.
#'
#' If the first argument is an [oce-class] object, then salinity, etc., are
#' extracted from it, and used for the calculation.
#'
#' @inheritParams swRho
#'
#' @return Quasi-potential density anomaly (kg/m\eqn{^3}{^3}), defined as the
#' density calculated with pressure set to zero.
#'
#' @author Dan Kelley
#'
#' @references See citations provided in the [swRho()] documentation.
#'
#' @examples
#' swSigmaT(35, 13, 1000, longitude = 300, latitude = 30, eos = "gsw") # 26.39623
#' swSigmaT(35, T90fromT68(13), 1000, eos = "unesco") # 26.39354
#'
#' @family functions that calculate seawater properties
swSigmaT <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
        swRho(l$salinity, l$temperature,
            pressure = rep(0, length(l$salinity)),
            longitude = l$longitude, latitude = l$latitude, eos = l$eos
        ) - 1000
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
        swRho(l$salinity, l$temperature,
            pressure = rep(0, length(l$salinity)),
            eos = l$eos
        ) - 1000
    }
}


#' Seawater Potential Density Anomaly
#'
#' Compute the potential density (minus 1000 kg/m^3) that seawater would have if raised
#' adiabatically to the surface.  In the UNESCO system, this quantity is
#' is denoted \eqn{\sigma_\theta}{sigma-theta} (hence the function name), but in
#' the GSW system, a somewhat related quantity is denoted `sigma0`. (In a
#' deep-water CTD cast, the RMS deviation between sigma-theta and sigma0 is
#' typically of order 0.0003 kg/m^3, corresponding to a temperature shift of
#' about 0.002C, so the distinction between the quantities is not large.)
#'
#' If the first argument is an `oce` object, then salinity, etc., are
#' extracted from it, and used for the calculation instead of any values
#' provided in the other arguments.
#'
#' @inheritParams swRho
#'
#' @param referencePressure The reference pressure, in dbar.
#'
#' @return Potential density anomaly (kg/m\eqn{^3}{^3}), defined as
#' \eqn{\sigma_\theta=\rho(S,\theta(S,t,p),0}{sigma_theta=rho(S,theta(S,t,p),0)}
#' - 1000 kg/m\eqn{^3}{^3}.
#'
#' @author Dan Kelley
#'
#' @references See citations provided in the [swRho()] documentation.
#'
#' @examples
#' stopifnot(abs(26.4212790994 - swSigmaTheta(35, 13, 1000, eos = "unesco")) < 1e-7)
#'
#' @family functions that calculate seawater properties
swSigmaTheta <- function(
    salinity, temperature = NULL, pressure = NULL, referencePressure = 0,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw"),
    debug = getOption("oceDebug")) {
    # message("DEBUG: in swSigmaTheta with eos=", eos, ", referencePressure=", referencePressure)
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    oceDebug(debug, "swSigmaTheta() ...\n")
    if (inherits(salinity, "oce")) {
        x <- salinity # store this for clarity
        oceDebug(debug, "  first parameter is an oce object, so extracting data from it\n")
        temperature <- x[["temperature"]]
        pressure <- x[["pressure"]]
        location <- locationForGsw(x)
        # cat(vectorShow(location))
        if (is.null(longitude)) {
            longitude <- location$longitude
        }
        if (is.null(latitude)) {
            latitude <- location$latitude
        }
        oceDebug(debug, vectorShow(longitude))
        oceDebug(debug, vectorShow(latitude))
        # we're done extracting other things, so now it's safe to rewrite 'salinity'
        salinity <- x[["salinity"]]
    }
    if (eos == "gsw") {
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
    }
    if (is.null(temperature)) {
        stop("must provide temperature")
    }
    if (is.null(pressure)) {
        stop("must provide pressure")
    }
    dim <- dim(salinity)
    nS <- length(salinity)
    nt <- length(temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    np <- length(pressure)
    if (np == 1) {
        pressure <- rep(pressure, length.out = nS)
    }
    np <- length(pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    referencePressure <- rep(referencePressure[1], length.out = nS)
    if (eos == "unesco") {
        theta <- swTheta(
            salinity = salinity, temperature = temperature, pressure = pressure,
            referencePressure = referencePressure,
            longitude = longitude, latitude = latitude, eos = eos
        )
        res <- swRho(
            salinity = salinity, temperature = theta, pressure = referencePressure,
            longitude = longitude, latitude = latitude, eos = eos
        ) - 1000.0
    } else if (eos == "gsw") {
        oceDebug(debug, "swSigmaTheta() handing following to gsw_SA_from_SP()\n")
        oceDebug(debug, vectorShow(as.vector(salinity)))
        oceDebug(debug, vectorShow(as.vector(pressure)))
        oceDebug(debug, vectorShow(as.vector(longitude)))
        oceDebug(debug, vectorShow(as.vector(latitude)))
        SA <- gsw::gsw_SA_from_SP(
            SP = as.vector(salinity),
            p = as.vector(pressure),
            longitude = as.vector(longitude),
            latitude = as.vector(latitude)
        )
        # CT <- gsw::gsw_CT_from_t(SA=SA, p=pressure, t=temperature)
        oceDebug(debug, "in swSigmaTheta about to call gsw_rho_t_exact()\n")
        oceDebug(debug, vectorShow(as.vector(SA)))
        oceDebug(debug, vectorShow(as.vector(temperature)))
        oceDebug(debug, vectorShow(as.vector(referencePressure)))
        res <- gsw::gsw_pot_rho_t_exact(
            SA = as.vector(SA),
            t = as.vector(temperature),
            p = as.vector(pressure),
            p_ref = as.vector(referencePressure)
        ) - 1000.0
        # res <- gsw::gsw_sigma0(SA=SA, CT=CT)
    } else {
        stop("eos must be either \"gsw\" or \"unesco\"; \"", eos, "\" is not acceptable.")
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}

#' Seawater Potential Density Anomaly Referenced to Surface Pressure
#'
#' Compute the potential density of seawater (minus
#' 1000 kg/m\eqn{^3}{^3}), referenced to surface pressure. This is done using
#' [gsw::gsw_sigma0()] if `eos="gsw"`, or using [swSigmaTheta()] if it is
#' `"unesco"`. (The difference between the formulations is typically
#' under 0.01 kg/m^3, corresponding to a few millidegrees of temperature.)
#'
#' @inheritParams swRho
#'
#' @return Potential density anomaly (kg/m\eqn{^3}{^3}).
#'
#' @author Dan Kelley
#'
#' @references See citations provided in the [swRho()] documentation.
#'
#' @family functions that calculate seawater properties
swSigma0 <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            temperature <- salinity[["temperature"]]
            pressure <- salinity[["pressure"]]
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
            salinity <- salinity[["salinity"]]
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude
        ))
        SA <- gsw::gsw_SA_from_SP(l$salinity, l$pressure, l$longitude, l$latitude)
        CT <- gsw::gsw_CT_from_t(SA, l$temperature, l$pressure)
        gsw::gsw_sigma0(SA = SA, CT = CT)
    } else if (eos == "unesco") {
        swSigmaTheta(
            salinity = salinity, temperature = temperature, pressure = pressure, referencePressure = 0,
            longitude = longitude, latitude = latitude, eos = eos
        )
    } else {
        stop("eos must be either \"gsw\" or \"unesco\", but it is \"", eos, "\"")
    }
}

#' Seawater Potential Density Anomaly Referenced to 1000db Pressure
#'
#' This is analogous to [swSigma0()], but referenced to 1000db pressure.
#'
#' @inheritParams swRho
#' @return Potential density anomaly (kg/m\eqn{^3}{^3}).
#' @author Dan Kelley
#' @references See citations provided in the [swRho()] documentation.
#' @family functions that calculate seawater properties
swSigma1 <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            temperature <- salinity[["temperature"]]
            pressure <- salinity[["pressure"]]
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
            salinity <- salinity[["salinity"]]
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude
        ))
        SA <- gsw::gsw_SA_from_SP(l$salinity, l$pressure, l$longitude, l$latitude)
        CT <- gsw::gsw_CT_from_t(SA, l$temperature, l$pressure)
        gsw::gsw_sigma1(SA = SA, CT = CT)
    } else if (eos == "unesco") {
        swSigmaTheta(
            salinity = salinity, temperature = temperature, pressure = pressure, referencePressure = 1000,
            longitude = longitude, latitude = latitude, eos = eos
        )
    } else {
        stop("eos must be either \"gsw\" or \"unesco\", but it is \"", eos, "\"")
    }
}

#' Seawater Potential Density Anomaly Referenced to 2000db Pressure
#'
#' This is analogous to [swSigma0()], but referenced to 2000db pressure.
#'
#' @inheritParams swRho
#' @return Potential density anomaly (kg/m\eqn{^3}{^3}).
#' @author Dan Kelley
#' @references See citations provided in the [swRho()] documentation.
#' @family functions that calculate seawater properties
swSigma2 <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            temperature <- salinity[["temperature"]]
            pressure <- salinity[["pressure"]]
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
            salinity <- salinity[["salinity"]]
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude
        ))
        SA <- gsw::gsw_SA_from_SP(l$salinity, l$pressure, l$longitude, l$latitude)
        CT <- gsw::gsw_CT_from_t(SA, l$temperature, l$pressure)
        gsw::gsw_sigma2(SA = SA, CT = CT)
    } else if (eos == "unesco") {
        swSigmaTheta(
            salinity = salinity, temperature = temperature, pressure = pressure, referencePressure = 2000,
            longitude = longitude, latitude = latitude, eos = eos
        )
    } else {
        stop("eos must be either \"gsw\" or \"unesco\", but it is \"", eos, "\"")
    }
}

#' Seawater Potential Density Anomaly Referenced to 3000db Pressure
#'
#' This is analogous to [swSigma0()], but referenced to 3000db pressure.
#'
#' @inheritParams swRho
#' @return Potential density anomaly (kg/m\eqn{^3}{^3}).
#' @author Dan Kelley
#' @references See citations provided in the [swRho()] documentation.
#' @family functions that calculate seawater properties
swSigma3 <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            temperature <- salinity[["temperature"]]
            pressure <- salinity[["pressure"]]
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
            salinity <- salinity[["salinity"]]
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude
        ))
        SA <- gsw::gsw_SA_from_SP(l$salinity, l$pressure, l$longitude, l$latitude)
        CT <- gsw::gsw_CT_from_t(SA, l$temperature, l$pressure)
        gsw::gsw_sigma3(SA = SA, CT = CT)
    } else if (eos == "unesco") {
        swSigmaTheta(
            salinity = salinity, temperature = temperature, pressure = pressure, referencePressure = 3000,
            longitude = longitude, latitude = latitude, eos = eos
        )
    } else {
        stop("eos must be either \"gsw\" or \"unesco\", but it is \"", eos, "\"")
    }
}

#' Seawater Potential Density Anomaly Referenced to 4000db Pressure
#'
#' This is analogous to [swSigma0()], but referenced to 4000db pressure.
#'
#' @inheritParams swRho
#' @return Potential density anomaly (kg/m\eqn{^3}{^3}).
#' @author Dan Kelley
#' @references See citations provided in the [swRho()] documentation.
#' @family functions that calculate seawater properties
swSigma4 <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            temperature <- salinity[["temperature"]]
            pressure <- salinity[["pressure"]]
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
            salinity <- salinity[["salinity"]]
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude
        ))
        SA <- gsw::gsw_SA_from_SP(l$salinity, l$pressure, l$longitude, l$latitude)
        CT <- gsw::gsw_CT_from_t(SA, l$temperature, l$pressure)
        gsw::gsw_sigma4(SA = SA, CT = CT)
    } else if (eos == "unesco") {
        swSigmaTheta(
            salinity = salinity, temperature = temperature, pressure = pressure, referencePressure = 4000,
            longitude = longitude, latitude = latitude, eos = eos
        )
    } else {
        stop("eos must be either \"gsw\" or \"unesco\", but it is \"", eos, "\"")
    }
}


#' Seawater Sound Absorption
#'
#' Compute the sound absorption of seawater, in dB/m
#'
#' Salinity and pH are ignored in this formulation.  Several formulae may be
#' found in the literature, and they give results differing by 10 percent, as
#' shown in reference 3 for example.  For this reason, it is likely that more
#' formulations will be added to this function, and entirely possible that the
#' default may change.
#'
#' @inheritParams swRho
#'
#' @param frequency The frequency of sound, in Hz.
#'
#' @param formulation character string indicating the formulation to use, either
#' of `"fischer-simmons"` or `"francois-garrison"`; see \dQuote{References}.
#'
#' @param pH seawater pH
#'
#' @return Sound absorption in dB/m.
#'
#' @author Dan Kelley
#'
#' @references
#' 1. F. H. Fisher and V. P. Simmons, 1977.  Sound absorption in
#' sea water.  Journal of the Acoustical Society of America, 62(3), 558-564.
#'
#' 2. R. E. Francois and G. R. Garrison, 1982.  Sound absorption based on
#' ocean measurements.  Part II: Boric acid contribution and equation for total
#' absorption.  Journal of the Acoustical Society of America, 72(6):1879-1890.
#'
#' 3. `http://resource.npl.co.uk/acoustics/techguides/seaabsorption/`
#'
#' @examples
#' # Fisher & Simmons (1977 table IV) gives 0.52 dB/km for 35 PSU, 5 degC, 500 atm
#' # (4990 dbar of water)a and 10 kHz
#' alpha <- swSoundAbsorption(35, 4, 4990, 10e3)
#'
#' # reproduce part of Fig 8 of Francois and Garrison (1982 Fig 8)
#' f <- 1e3 * 10^(seq(-1, 3, 0.1)) # in KHz
#' plot(f / 1000, 1e3 * swSoundAbsorption(f, 35, 10, 0, formulation = "fr"),
#'     xlab = " Freq [kHz]", ylab = " dB/km", type = "l", log = "xy"
#' )
#' lines(f / 1000, 1e3 * swSoundAbsorption(f, 0, 10, 0, formulation = "fr"), lty = "dashed")
#' legend("topleft", lty = c("solid", "dashed"), legend = c("S=35", "S=0"))
#'
#' @family functions that calculate seawater properties
swSoundAbsorption <- function(
    frequency, salinity, temperature, pressure, pH = 8,
    formulation = c("fisher-simmons", "francois-garrison")) {
    formulation <- match.arg(formulation)
    # nolint start T_and_F_symbol_linter
    if (formulation == "fisher-simmons") {
        # Equation numbers are from Fisher & Simmons (1977); see help page for ref
        p <- 1 + pressure / 10 # add atmophere, then convert water part from dbar
        S <- salinity
        T <- T68fromT90(temperature)
        f <- frequency
        A1 <- 1.03e-8 + 2.36e-10 * T - 5.22e-12 * T^2 # (5)
        f1 <- 1.32e3 * (T + 273.1) * exp(-1700 / (T + 273.1)) # (6)
        A2 <- 5.62e-8 + 7.52e-10 * T # (7)
        f2 <- 1.55e7 * (T + 273.1) * exp(-3052 / (T + 273.1)) # (8)
        P2 <- 1 - 10.3e-4 * p + 3.7e-7 * p^2 # (9)
        A3 <- (55.9 - 2.37 * T + 4.77e-2 * T^2 - 3.48e-4 * T^3) * 1e-15 # (10)
        P3 <- 1 - 3.84e-4 * p + 7.57e-8 * p^2 # (11)
        alpha <- (A1 * f1 * f^2) / (f1^2 + f^2) + (A2 * P2 * f2 * f^2) / (f2^2 + f^2) + A3 * P3 * f^2 # (3a)
        alpha <- alpha * 8686 / 1000 # dB/m
    } else if (formulation == "francois-garrison") {
        S <- salinity
        T <- T68fromT90(temperature)
        D <- pressure # FIXME: approximation
        c <- 1412 + 3.21 * T + 1.19 * S + 0.0167 * D # sound speed m/s
        f <- frequency / 1000 # convert to kHz
        theta <- 273 + T
        # f1 in kHz
        f1 <- 2.8 * sqrt(S / 35) * 10^(4 - 1245 / theta) # nolint (space before left parenthesis)
        # subscript 1 for boric acid contribution
        # A1 in dB / km / kHz
        A1 <- 8.86 / c * 10^(0.78 * pH - 5) # nolint (space before left parenthesis)
        P1 <- 1
        # MgSO4 contribution
        A2 <- 21.44 * (S / c) * (1 + 0.025 * T) # dB / km / kHz
        P2 <- 1 - 1.37e-4 * D + 6.2e-9 * D^2
        # f2 in kHz
        f2 <- (8.17 * 10^(8 - 1990 / theta)) / (1 + 0.0018 * (S - 35)) # nolint (space before left parenthesis)
        # pure water contribution
        A3 <- 3.964e-4 - 1.146e-5 * T + 1.45e-7 * T^2 - 6.5e-10 * T^3 # dB / km / kHz^2
        P3 <- 1 - 3.83e-5 * D + 4.9e-10 * D^2
        alpha <- (A1 * P1 * f1 * f^2) / (f^2 + f1^2) + (A2 * P2 * f2 * f^2) / (f^2 + f2^2) + A3 * P3 * f^2
        alpha <- alpha / 1000
    }
    # nolint end T_and_F_symbol_linter
    alpha
}


#' Seawater Sound Speed
#'
#' Compute the seawater speed of sound.
#'
#' If `eos="unesco"`, the sound speed is calculated using the formulation
#' in section 9 of Fofonoff and Millard (1983). If `eos="gsw"`, then the
#' [gsw::gsw_sound_speed()] function from the
#' \CRANpkg{gsw} package is used.
#'
#' @inheritParams swRho
#'
#' @return Sound speed (m/s).
#'
#' @author Dan Kelley
#'
#' @references Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for
#' computation of fundamental properties of seawater.
#' *Unesco Technical Papers in Marine Science*,
#' *44*, 53 pp.  (See section 9.)
#'
#' @examples
#' swSoundSpeed(40, T90fromT68(40), 10000) # 1731.995 (p48 of Fofonoff + Millard 1983)
#'
#' @family functions that calculate seawater properties
swSoundSpeed <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
    }
    dim <- dim(l$salinity)
    if (is.null(l$temperature)) {
        stop("must provide temperature")
    }
    if (is.null(l$pressure)) {
        stop("must provide pressure")
    }
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    # np <- length(l$pressure)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    l$pressure <- rep(l$pressure, length.out = nS)
    if (eos == "unesco") {
        res <- .C("sw_svel", as.integer(nS), as.double(l$salinity), as.double(T68fromT90(l$temperature)), as.double(l$pressure),
            value = double(nS), NAOK = TRUE, PACKAGE = "oce"
        )$value
    } else if (eos == "gsw") {
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        CT <- gsw::gsw_CT_from_t(SA = SA, t = l$temperature, p = l$pressure)
        res <- gsw::gsw_sound_speed(SA = SA, CT = CT, p = l$pressure)
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Specific Heat
# Source= http://sam.ucsd.edu/sio210/propseawater/ppsw_fortran/ppsw.f
# check value: cpsw = 3849.500 j/(kg deg. c) for s = 40 (ipss-78),
#'
#' Compute specific heat of seawater.
#'
#' If the first argument is a `ctd` object, then salinity, etc, are
#' extracted from it, and used for the calculation.
#'
#' @param salinity either practical salinity (in which case `temperature`
#' and `pressure` must be provided) *or* an `oce` object (in
#' which case `salinity`, etc. are inferred from the object).
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale.
#' @param pressure seawater pressure (dbar)
#' @param longitude longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#' @param latitude latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#' @param eos equation of state, either `"unesco"` or `"gsw"`.
#' @return Specific heat (J/kg/degC).
#' @author Dan Kelley
#' @references Millero et. al., J. Geophys. Res. 78 (1973), 4499-4507
#'
#' Millero et. al., UNESCO report 38 (1981), 99-188.
#' @examples
#' swSpecificHeat(40, T90fromT68(40), 10000, eos = "unesco") # 3949.499
#'
#' @family functions that calculate seawater properties
swSpecificHeat <- function(
    salinity, temperature = NULL, pressure = 0,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (eos == "gsw") {
        if (inherits(salinity, "oce")) {
            if (is.null(longitude)) {
                longitude <- salinity[["longitude"]]
            }
            if (is.null(latitude)) {
                latitude <- salinity[["latitude"]]
            }
        }
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure,
            longitude = longitude, latitude = latitude, eos = eos
        ))
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure, eos = eos))
    }
    if (is.null(l$temperature)) {
        stop("must provide temperature")
    }
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    if (length(l$pressure) == 1) l$pressure <- rep(l$pressure, length.out = nS)
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    if (eos == "unesco") {
        res <- .Fortran("cp_driver", as.double(l$salinity), as.double(T68fromT90(l$temperature)), as.double(l$pressure),
            as.integer(nS),
            CP = double(nS)
        )$CP
    } else {
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        res <- gsw::gsw_cp_t_exact(SA = SA, t = l$temperature, p = l$pressure)
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Spiciness
#'
#' Compute seawater "spice", a variable that is in some sense orthogonal to
#' density in TS space.  Larger spice values correspond to relative warm and
#' salty water, compared with smaller spice values. Two distinct variants exist.
#' If `eos="unesco"` then Flament's (2002) formulation is used. If `eos="gsw"`
#' then [gsw::gsw_spiciness0()] is used to compute a newer variant that is part
#' of the Gibbs SeaWater formulation (McDougall and Krzysik, 2015). See the
#' \dQuote{Examples} section for a graphical illustration of the difference in a
#' typical coastal case.
#'
#' If the first argument is a `ctd` object, then salinity, temperature and
#' pressure values are extracted from it, and used for the calculation. For
#' the `eos="gsw"` case, longitude and latitude are also extracted, because
#' these are required by [gsw::gsw_spiciness0()].
#'
#' Roughly speaking, seawater with a high spiciness is relatively warm and
#' salty compared with less spicy water. Another interpretation is that spice
#' is a variable measuring distance orthogonal to isopycnal lines on TS
#' diagrams (if the diagrams are scaled to make the isopycnals run at 45
#' degrees). Note that pressure, longitude and latitude are all
#' ignored in the Flament definition.
#'
#' @param salinity either salinity (PSU) (in which case `temperature` and
#' `pressure` must be provided) *or* a `ctd` object (in which
#' case `salinity`, `temperature` and `pressure` are determined
#' from the object, and must not be provided in the argument list).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C) on the
#' ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure Seawater pressure (dbar) (only used if `eos` is
#' `"gsw"`); see \dQuote{Details}..
#'
#' @param longitude longitude of observation (only used if `eos` is
#' `"gsw"`; see \dQuote{Details}).
#'
#' @param latitude latitude of observation (only used if `eos` is
#' `"gsw"`; see \dQuote{Details}).
#'
#' @param eos Character value specifying the equation of state, either
#' `"unesco"` (for the Flament formulation, although this is not actually part
#' of UNESCO) or `"gsw"` for the Gibbs SeaWater formulation.
#'
#' @template debugTemplate
#'
#' @return Flament-formulated spice \eqn{kg/m^3} if `eos` is `"unesco"`
#' or surface-referenced GSW spiciness0 \eqn{kg/m^3} if `eos` is `"gsw"`,
#' the latter provided by [gsw::gsw_spiciness0()], and hence aimed
#' at application within the top half-kilometre of the ocean.
#'
#' @author Dan Kelley coded this, merely an interface to the code described
#' by references 1 and 2.
#'
#' @examples
#' # Compare unesco and gsw formulations
#' library(oce)
#' data(ctd)
#' p <- ctd[["pressure"]]
#' U <- swSpice(ctd, eos = "unesco")
#' G <- swSpice(ctd, eos = "gsw")
#' xlim <- range(c(U, G), na.rm = TRUE)
#' ylim <- rev(range(p))
#' plot(U, p,
#'     xlim = xlim, ylim = ylim,
#'     xlab = "Measure of Spiciness", ylab = "Pressure (dbar)"
#' )
#' points(G, p, col = 2)
#' legend("topleft", col = 1:2, pch = 1, legend = c("unesco", "gsw"))
#'
#' @references
#' 1. Flament, P. \dQuote{A State Variable for Characterizing Water Masses and Their
#' Diffusive Stability: Spiciness.} Progress in Oceanography, Observations of the
#' 1997-98 El Nino along the West Coast of North America, 54, no. 1
#' (July 1, 2002):493-501.
#' \doi{10.1016/S0079-6611(02)00065-4}
#'
#' 2. McDougall, Trevor J., and Oliver A. Krzysik. \dQuote{Spiciness.}
#' Journal of Marine Research 73, no. 5 (September 1, 2015): 141-52.
#'
#' @family functions that calculate seawater spiciness
#' @family functions that calculate seawater properties
swSpice <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw"),
    debug = getOption("oceDebug")) {
    oceDebug(debug, vectorShow(eos))
    if (eos != "gsw" && eos != "unesco") {
        stop("eos must be \"gsw\" or \"unesco\", but \"", eos, "\" was given")
    }
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        temperature <- salinity[["temperature"]]
        pressure <- salinity[["pressure"]]
        if (is.null(longitude)) {
            longitude <- salinity[["longitude"]]
        }
        if (is.null(latitude)) {
            latitude <- salinity[["latitude"]]
        }
        salinity <- salinity[["salinity"]]
    }
    if (eos == "gsw") {
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
    }
    if (eos == "gsw" && is.null(pressure)) {
        stop("must provide pressure")
    }
    dim <- dim(salinity)
    nS <- length(salinity)
    nt <- length(temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    if (is.null(pressure)) pressure <- rep(0, nS)
    if (length(pressure) == 1) pressure <- rep(pressure, length.out = nS)
    np <- length(pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    if (eos == "unesco") {
        res <- .C("sw_spice", as.integer(nS), as.double(salinity),
            as.double(T68fromT90(temperature)), as.double(pressure),
            value = double(nS), NAOK = TRUE, PACKAGE = "oce"
        )$value
    } else if (eos == "gsw") {
        SA <- gsw::gsw_SA_from_SP(SP = salinity, p = pressure, longitude = longitude, latitude = latitude)
        CT <- gsw::gsw_CT_from_t(SA = SA, t = temperature, p = pressure)
        res <- gsw::gsw_spiciness0(SA, CT)
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}

# Internal function used by swSpiciness0(), swSpiciness1() and swSpiciness2()
swSpiciness <- function(
    salinity = NULL, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL, referencePressure = NULL) {
    if (is.null(salinity)) stop("must provide salinity")
    if (inherits(salinity, "oce")) {
        temperature <- salinity[["temperature"]]
        pressure <- salinity[["pressure"]]
        longitude <- salinity[["longitude"]]
        latitude <- salinity[["latitude"]]
        salinity <- salinity[["salinity"]] # overwrites
    }
    if (is.null(temperature)) stop("must supply temperature")
    if (is.null(longitude)) stop("must supply longitude")
    if (is.null(latitude)) stop("must supply latitude")
    if (is.null(pressure)) stop("must provide pressure")
    if (is.null(referencePressure)) stop("must provide referencePressure")
    if (length(referencePressure) != 1) stop("referencePressure must be of length 1")
    if (referencePressure != 0 && referencePressure != 1000 && referencePressure != 2000) {
        stop("referencePressure must be 0, 1000 or 2000, but it is ", referencePressure)
    }
    dim <- dim(salinity)
    nS <- length(salinity)
    nT <- length(temperature)
    np <- length(pressure)
    if (nS != nT) stop("lengths of salinity (", nS, ") and temperature (", nT, ") disagree")
    if (nS != np) stop("lengths of salinity (", nS, ") and pressure (", np, ") disagree")
    # FIXME: what if e.g. length(latitude) == length(salinity) etc?
    SA <- gsw::gsw_SA_from_SP(SP = salinity, p = pressure, longitude = longitude, latitude = latitude)
    CT <- gsw::gsw_CT_from_t(SA = SA, t = temperature, p = pressure)
    res <- if (referencePressure[1] == 0) {
        gsw::gsw_spiciness0(SA, CT)
    } else if (referencePressure[1] == 1000) {
        gsw::gsw_spiciness1(SA, CT)
    } else if (referencePressure[1] == 2000) {
        gsw::gsw_spiciness2(SA, CT)
    }
    if (!is.null(dim)) dim(res) <- dim
    res
}

#' Spiciness in gsw System, Referenced to Surface Pressure
#'
#' Computes seawater spiciness using [gsw::gsw_spiciness0()] for surface
#' referenced computation.
#'
#' @param salinity either salinity, or an oce object that contains salinity,
#' temperature, pressure, longitude and latitude.
#'
#' @param temperature in-situ temperature (ignored if `salinity` is an oce object)
#'
#' @param pressure seawater pressure in dbar (ignored if `salinity` is an oce object)
#'
#' @param longitude,latitude observation location (ignored if `salinity` is an oce object).
#'
#' @return seawater spiciness with respect to a reference pressure of 0 dbar
#' (that is, the sea surface), as defined in the `gsw` (TEOS-10)
#' system (McDougall et al, 2011).
#'
#' @references
#' McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127,
#' ISBN 978-0-646-55621-5.
#'
#' @family functions that calculate seawater spiciness
#' @family functions that calculate seawater properties
#'
#' @author Dan Kelley
swSpiciness0 <- function(salinity, temperature, pressure, longitude, latitude) {
    swSpiciness(
        salinity = salinity, temperature = temperature, pressure = pressure,
        longitude = longitude, latitude = latitude, referencePressure = 0
    )
}


#' Spiciness in gsw System, Referenced to 1000 dbar Pressure
#'
#' Computes seawater spiciness using [gsw::gsw_spiciness1()] for
#' a reference pressure of 1000 dbar.
#'
#' @param salinity either salinity, or an oce object that contains salinity,
#' temperature, pressure, longitude and latitude.
#'
#' @param temperature in-situ temperature (ignored if `salinity` is an oce object)
#'
#' @param pressure seawater pressure in dbar (ignored if `salinity` is an oce object)
#'
#' @param longitude,latitude observation location (ignored if `salinity` is an oce object).
#'
#' @return seawater spiciness with respect to a reference pressure of 1000 dbar,
#' as defined in the `gsw` (TEOS-10)
#' system (McDougall et al, 2011) and computed with [gsw::gsw_spiciness1()].
#'
#' @references
#' McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127,
#' ISBN 978-0-646-55621-5.
#'
#' @family functions that calculate seawater spiciness
#' @family functions that calculate seawater properties
#'
#' @author Dan Kelley
swSpiciness1 <- function(salinity, temperature, pressure, longitude, latitude) {
    swSpiciness(
        salinity = salinity, temperature = temperature, pressure = pressure,
        longitude = longitude, latitude = latitude, referencePressure = 1000
    )
}

#' Spiciness in gsw System, Referenced to 2000 dbar Pressure
#'
#' Computes seawater spiciness using [gsw::gsw_spiciness2()] for
#' a reference pressure of 2000 dbar.
#'
#' @param salinity either salinity, or an oce object that contains salinity,
#' temperature, pressure, longitude and latitude.
#'
#' @param temperature in-situ temperature (ignored if `salinity` is an oce object)
#'
#' @param pressure seawater pressure in dbar (ignored if `salinity` is an oce object)
#'
#' @param longitude,latitude observation location (ignored if `salinity` is an oce object).
#'
#' @return seawater spiciness with respect to a reference pressure of 2000 dbar,
#' as defined in the `gsw` (TEOS-10)
#' system (McDougall et al, 2011) and computed with [gsw::gsw_spiciness2()].
#'
#' @references
#' McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127,
#' ISBN 978-0-646-55621-5.
#'
#' @family functions that calculate seawater spiciness
#' @family functions that calculate seawater properties
#'
#' @author Dan Kelley
swSpiciness2 <- function(salinity, temperature, pressure, longitude, latitude) {
    swSpiciness(
        salinity = salinity, temperature = temperature, pressure = pressure,
        longitude = longitude, latitude = latitude, referencePressure = 2000
    )
}


#' Seawater Potential Temperature (UNESCO Version)
#'
#' Compute the potential temperature of seawater, denoted \eqn{\theta}{theta}
#' in the UNESCO system, and `pt` in the GSW system.
#'
#' Different formulae are used depending on the equation of state. If `eos`
#' is `"unesco"`, the method of Fofonoff *et al.* (1983) is used
#' (see references 1 and 2).
#' Otherwise, `swTheta` uses [gsw::gsw_pt_from_t()] from
#' the \CRANpkg{gsw} package.
#'
#' If the first argument is a `ctd` or `section` object, then values
#' for salinity, etc., are extracted from it, and used for the calculation, and
#' the corresponding arguments to the present function are ignored.
#'
#' @param salinity either salinity (PSU) (in which case `temperature` and
#' `pressure` must be provided) *or* an `oce` object (in which
#' case `salinity`, etc. are inferred from the object).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()], and the examples below.
#'
#' @param pressure pressure (dbar)
#'
#' @param referencePressure reference pressure (dbar)
#'
#' @param longitude longitude of observation (only used if `eos="gsw"`;
#' see \dQuote{Details}).
#'
#' @param latitude latitude of observation (only used if `eos="gsw"`; see
#' \dQuote{Details}).
#'
#' @param eos equation of state, either `"unesco"` (references 1 and 2) or `"gsw"`
#' (references 3 and 4).
#'
#' @template debugTemplate
#'
#' @return Potential temperature (\eqn{^\circ}{deg}C) of seawater, referenced
#' to pressure `referencePressure`.
#'
#' @author Dan Kelley
#'
#' @references
#'
#' 1. Fofonoff, P. and R. C. Millard Jr, 1983. Algorithms for computation of
#' fundamental properties of seawater.
#' *Unesco Technical Papers in Marine Science*, *44*, 53 pp
#'
#' 2. Gill, A.E., 1982. *Atmosphere-ocean Dynamics*, Academic Press, New
#' York, 662 pp.
#'
#' 3. IOC, SCOR, and IAPSO (2010). The international thermodynamic equation of
#' seawater-2010: Calculation and use of thermodynamic properties.  Technical
#' Report 56, Intergovernmental Oceanographic Commission, Manuals and Guide.
#'
#' 4. McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and
#' the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127,
#' ISBN 978-0-646-55621-5.
#'
#' @examples
#' library(oce)
#' # Example 1: test value from Fofonoff et al., 1983
#' stopifnot(abs(36.8818748026 - swTheta(40, T90fromT68(40), 10000, 0, eos = "unesco")) < 0.0000000001)
#'
#' # Example 2: a deep-water station. Note that theta and CT are
#' # visually identical on this scale.
#' data(section)
#' stn <- section[["station", 70]]
#' plotProfile(stn, "temperature", ylim = c(6000, 1000))
#' lines(stn[["theta"]], stn[["pressure"]], col = 2)
#' lines(stn[["CT"]], stn[["pressure"]], col = 4, lty = 2)
#' legend("bottomright",
#'     lwd = 1, col = c(1, 2, 4), lty = c(1, 1, 2),
#'     legend = c("in-situ", "theta", "CT"),
#'     title = sprintf("MAD(theta-CT)=%.4f", mean(abs(stn[["theta"]] - stn[["CT"]])))
#' )
#'
#' @family functions that calculate seawater properties
swTheta <- function(
    salinity, temperature = NULL, pressure = NULL, referencePressure = 0,
    longitude = NULL, latitude = NULL, eos = getOption("oceEOS", default = "gsw"),
    debug = getOption("oceDebug")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    oceDebug(debug, "swTheta(..., referencePressure=", referencePressure, ", eos=\"", eos, "\")\n")
    if (inherits(salinity, "oce")) {
        oceDebug(debug, "  first parameter is an oce object, so extracting data from it\n")
        x <- salinity # store this for clarity
        temperature <- x[["temperature"]]
        pressure <- x[["pressure"]]
        location <- locationForGsw(x)
        longitude <- location$longitude
        latitude <- location$latitude
        oceDebug(debug, vectorShow(longitude))
        oceDebug(debug, vectorShow(latitude))
        salinity <- x[["salinity"]]
    }
    if (eos == "gsw") {
        if (is.null(longitude)) {
            stop("must supply longitude")
        }
        if (is.null(latitude)) {
            stop("must supply latitude")
        }
        l <- lookWithin(list(
            salinity = salinity, temperature = temperature, pressure = pressure
        ))
        l$latitude <- latitude
        l$longitude <- longitude
    } else {
        l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure))
    }
    # message("sw.R:3147")
    # cat(str(l))
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    np <- length(l$pressure)
    if (np == 1L) {
        l$pressure <- rep(l$pressure, length.out = nS)
    }
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    referencePressure <- rep(referencePressure[1], length.out = nS)
    if (eos == "gsw") {
        SA <- gsw::gsw_SA_from_SP(SP = l$salinity, p = l$pressure, longitude = l$longitude, latitude = l$latitude)
        res <- gsw::gsw_pt_from_t(SA = SA, t = l$temperature, p = l$pressure, p_ref = referencePressure)
    } else if (eos == "unesco") {
        # message("DAN it is unesco")
        # Note the conversion to the T68 scale, because that's the scale
        # used by the UNESCO formula.
        res <- .C("theta_UNESCO_1983",
            as.integer(nS),
            as.double(l$salinity),
            as.double(T68fromT90(l$temperature)),
            as.double(l$pressure),
            as.double(referencePressure),
            value = double(nS), NAOK = TRUE, PACKAGE = "oce"
        )$value
        res <- T90fromT68(res)
    }
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Viscosity
#'
#' Compute viscosity of seawater, in \eqn{Pa\cdot s}{Pa*s}
#'
#' If the first argument is a `ctd` object, then salinity, temperature and
#' pressure values are extracted from it, and used for the calculation.
#'
#' The result is determined from a regression of the data provided in Table 87
#' of Dorsey (1940). The fit matches the table to within 0.2 percent at worst,
#' and with average absolute error of 0.07 percent. The maximum deviation from
#' the table is one unit in the last decimal place.
#'
#' No pressure dependence was reported by Dorsey (1940).
#'
#' @param salinity either salinity (PSU) (in which case `temperature` and
#' `pressure` must be provided) *or* a `ctd` object (in which
#' case `salinity`, `temperature` and `pressure` are determined
#' from the object, and must not be provided in the argument list).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()], and the examples below.
#'
#' @return Viscosity of seawater in \eqn{Pa\cdot s}{Pa*s}.  Divide by density
#' to get kinematic viscosity in \eqn{m^2/s}{m^2/s}.
#'
#' @author Dan Kelley
#'
#' @references N. Ernest Dorsey (1940),
#' *Properties of ordinary Water-substance*,
#' American Chemical Society Monograph Series.  Reinhold
#' Publishing.
#'
#' @examples
#' swViscosity(30, 10) # 0.001383779
#'
#' @family functions that calculate seawater properties
swViscosity <- function(salinity, temperature) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (missing(temperature)) {
        stop("must provide temperature")
    }
    l <- lookWithin(list(salinity = salinity, temperature = temperature))
    0.001798525 + l$salinity * (2.634749e-06 - 7.088328e-10 *
        l$temperature^2 + l$salinity * (-4.702342e-09 + l$salinity *
        (5.32178e-11))) + l$temperature * (-6.293088e-05 +
        l$temperature * (1.716685e-06 + l$temperature * (-3.479273e-08
        + l$temperature * (+3.566255e-10))))
}


#' Seawater Conservative Temperature (GSW Formulation)
#'
#' Compute seawater Conservative Temperature, according to the GSW/TEOS-10
#' formulation.
#'
#' If the first argument is an `oce` object, then values for salinity,
#' etc., are extracted from it, and used for the calculation, and the
#' corresponding arguments to the present function are ignored.
#'
#' The conservative temperature is calculated using the TEOS-10 function
#' [gsw::gsw_CT_from_t] from the \CRANpkg{gsw} package.
#'
#' @param salinity either practical salinity (in which case `temperature`
#' and `pressure` must be provided) *or* an `oce` object (in
#' which case `salinity`, etc. are inferred from the object).
#'
#' @param temperature *in-situ* temperature (\eqn{^\circ}{deg}C), defined
#' on the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure pressure (dbar)
#'
#' @param longitude longitude of observation.
#'
#' @param latitude latitude of observation.
#'
#' @template debugTemplate
#'
#' @return Conservative temperature in degrees Celcius.
#'
#' @author Dan Kelley
#'
#' @seealso The related TEOS-10 quantity ``absolute salinity'' may be computed
#' with [swAbsoluteSalinity()].  For a ctd object, conservative
#' temperature may also be recovered by indexing as e.g.
#' \code{ctd[["conservativeTemperature"]]} or \code{ctd[["CT"]]}.
# NOTE: the markdown-Rd translator balks on the above if backticks are used
#'
#' @references McDougall, T.J. and P.M. Barker, 2011: Getting started with
#' TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp.,
#' SCOR/IAPSO WG127, ISBN 978-0-646-55621-5.
#'
#' @examples
#' swConservativeTemperature(35, 10, 1000, 188, 4) # 9.86883
#'
#' @family functions that calculate seawater properties
swConservativeTemperature <- function(
    salinity, temperature = NULL, pressure = NULL,
    longitude = NULL, latitude = NULL,
    debug = getOption("oceDebug")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        x <- salinity # store this for clarity
        oceDebug(debug, "swConservativeTemperature() with an oce object as first argument ...\n")
        location <- locationForGsw(x)
        if (is.null(longitude)) {
            longitude <- location$longitude
        }
        if (is.null(latitude)) {
            latitude <- location$latitude
        }
        oceDebug(debug, vectorShow(longitude))
        oceDebug(debug, vectorShow(latitude))
    }
    if (is.null(longitude) || is.null(latitude)) {
        stop("need longitude and latitude to compute CT")
    }
    # we already have properly dimensioned longitude and latitude, so don't look
    # within the object for them.
    l <- lookWithin(list(salinity = salinity, temperature = temperature, pressure = pressure))
    l$longitude <- longitude
    l$latitude <- latitude
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    nt <- length(l$temperature)
    if (nS != nt) {
        stop("lengths of salinity and temperature must agree, but they are ", nS, " and ", nt, ", respectively")
    }
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    # make longitude and latitude dimensions match salinity
    if (length(l$longitude) < nS) {
        l$longitude <- rep(l$longitude, length.out = nS)
        dim(l$longitude) <- dim
    }
    if (length(l$latitude) < nS) {
        l$latitude <- rep(l$latitude, length.out = nS)
        dim(l$latitude) <- dim
    }
    bad <- is.na(l$salinity) | is.na(l$temperature) | is.na(l$pressure)
    SA <- gsw::gsw_SA_from_SP(
        SP = l$salinity[!bad], p = l$pressure[!bad],
        longitude = l$longitude[!bad], latitude = l$latitude[!bad]
    )
    good <- gsw::gsw_CT_from_t(SA = SA, t = l$temperature[!bad], p = l$pressure[!bad])
    res <- rep(NA, nS)
    res[!bad] <- good
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Absolute Salinity (GSW Formulation)
#'
#' Compute the seawater Absolute Salinity, according to the GSW/TEOS-10
#' formulation with [gsw::gsw_SA_from_SP()] in the \CRANpkg{gsw} package.
#' Typically, this is a fraction of a unit
#' higher than practical salinity as defined in the UNESCO formulae.
#'
#' @param salinity either practical salinity (in which case `temperature`
#' and `pressure` must be provided) *or* an `oce` object (in
#' which case `salinity`, etc. are inferred from the object).
#'
#' @param pressure pressure in dbar.
#'
#' @param longitude longitude of observation.
#'
#' @param latitude latitude of observation.
#'
#' @template debugTemplate
#'
#' @return Absolute Salinity in \eqn{g/kg}{g/kg}.
#'
#' @author Dan Kelley
#'
#' @seealso The related TEOS-10 quantity ``conservative temperature'' may be
#' computed with [swConservativeTemperature()].  For a ctd object,
#' absolute salinity may also be recovered by indexing as e.g.
#' \code{ctd[["absoluteSalinity"]]} or \code{ctd[["SA"]]}.
# NOTE: the markdown-Rd translator balks on the above if backticks are used
#'
#' @references McDougall, T.J. and P.M. Barker, 2011: Getting started with
#' TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp.,
#' SCOR/IAPSO WG127, ISBN 978-0-646-55621-5.
#'
#' @examples
#' swAbsoluteSalinity(35.5, 300, 260, 16) # 35.67136
#'
#' @family functions that calculate seawater properties
swAbsoluteSalinity <- function(salinity, pressure = NULL, longitude = NULL, latitude = NULL,
                               debug = getOption("oceDebug")) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        x <- salinity # store this for clarity
        if (!"salinity" %in% names(x@data)) {
            stop("this oce object lacks salinity, so SA cannot be computed")
        }
        oceDebug(debug, "swAbsoluteSalinity() with an oce object as first argument ...\n")
        location <- locationForGsw(x)
        # cat(vectorShow(location))
        if (is.null(longitude)) {
            longitude <- location$longitude
        }
        if (is.null(latitude)) {
            latitude <- location$latitude
        }
        oceDebug(debug, vectorShow(longitude))
        oceDebug(debug, vectorShow(latitude))
        if (is.null(pressure)) {
            pressure <- x[["pressure"]]
        }
    }
    if (is.null(longitude) || is.null(latitude)) {
        stop("need longitude and latitude to compute SA")
    }
    if (is.null(pressure)) {
        stop("need pressure to compute SA")
    }
    if (length(longitude) != length(pressure)) {
        stop("lengths of longitude (", length(longitude), ") and pressure (", length(pressure), ") do not match")
    }
    if (length(latitude) != length(pressure)) {
        stop("lengths of latitude (", length(latitude), ") and pressure (", length(pressure), ") do not match")
    }
    # we already have properly dimensioned longitude and latitude, so don't look
    # within the object for them.
    l <- lookWithin(list(salinity = salinity, pressure = pressure))
    l$longitude <- longitude
    l$latitude <- latitude
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity and pressure must agree, but they are ", nS, " and ", np, ", respectively")
    }
    # make longitude and latitude dimensions match salinity
    if (length(l$longitude) < nS) {
        l$longitude <- rep(l$longitude, length.out = nS)
        dim(l$longitude) <- dim
    }
    if (length(l$latitude) < nS) {
        l$latitude <- rep(l$latitude, length.out = nS)
        dim(l$latitude) <- dim
    }
    bad <- is.na(l$salinity) | is.na(l$pressure) | is.na(l$longitude) | is.na(l$latitude)
    good <- gsw::gsw_SA_from_SP(l$salinity[!bad], l$pressure[!bad], l$longitude[!bad], l$latitude[!bad])
    res <- rep(NA, nS)
    res[!bad] <- good
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Preformed Salinity (GSW Formulation)
#'
#' Compute seawater Preformed Salinity (S*), according to the GSW/TEOS-10
#' formulation with [gsw::gsw_Sstar_from_SA()] in the \CRANpkg{gsw} package.
#'
#' @param salinity either practical salinity (in which case
#' `pressure` must be provided) *or* an `oce` object with
#' `salinity` and `pressure` in its data slot, and with
#' `longitude` and `latitude` either there, or in the metadata
#' slot.
#'
#' @param pressure pressure in dbar.
#'
#' @param longitude longitude of observation.
#'
#' @param latitude latitude of observation.
#'
#' @return Preformed Salinity, S*, in \eqn{g/kg}{g/kg}.
#'
#' @author Dan Kelley
#'
#' @seealso For some objects, S-star may also be recovered by
#' indexing as e.g. \code{ctd[["Sstar"]]}.
#'
#' @references McDougall, T.J. and P.M. Barker, 2011: Getting started with
#' TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp.,
#' SCOR/IAPSO WG127, ISBN 978-0-646-55621-5.
#'
#' @examples
#' swSstar(35.5, 300, 260, 16) # 35.66601
#'
#' @family functions that calculate seawater properties
swSstar <- function(salinity, pressure = NULL, longitude = NULL, latitude = NULL) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        x <- salinity # store this for clarity
        for (item in c("salinity", "pressure")) {
            if (!item %in% names(x@data)) {
                stop("this oce object lacks ", item, ", so SA cannot be computed")
            }
        }
        # https://github.com/dankelley/oce/issues/1911
        # Create new longitude and latitude, to match the length or
        # dimensionality of pressure.  The argo-class is a special case,
        # because it has vectors of longitude and latitude, of length
        # that matches dim(pressure)[2].  However, it seems prudent
        # to base our method on the shapes of the data, not on
        # the object class, as a way to address not only argo but
        # also perhaps any future-defined object that has this
        # characteristic.
        pressure <- x@data$pressure
        # temperature <- x@data$temperature
        salinity <- x@data$salinity
        if (is.null(longitude)) {
            longitude <- x[["longitude"]]
        }
        if (is.null(latitude)) {
            latitude <- x[["latitude"]]
        }
        if (is.null(longitude) || is.null(latitude)) {
            stop("need longitude and latitude to compute Sstar")
        }
        if (is.array(pressure)) {
            nlevels <- dim(pressure)[1]
            longitude <- rep(longitude, each = nlevels)
            latitude <- rep(latitude, each = nlevels)
        } else {
            np <- length(pressure)
            longitude <- rep(longitude[1], length.out = np)
            latitude <- rep(latitude[1], length.out = np)
        }
    }
    if (is.null(longitude) || is.null(latitude)) {
        stop("need longitude and latitude to compute Sstar")
    }
    if (length(longitude) != length(pressure)) {
        stop("lengths of longitude and pressure must match")
    }
    if (length(latitude) != length(pressure)) {
        stop("lengths of latitude and pressure must match")
    }
    l <- lookWithin(list(salinity = salinity, pressure = pressure, longitude = longitude, latitude = latitude))
    dim <- dim(l$salinity)
    nS <- length(l$salinity)
    np <- length(l$pressure)
    if (nS != np) {
        stop("lengths of salinity (", nS, ") and pressure (", np, ") disagree")
    }
    bad <- is.na(l$salinity) | is.na(l$pressure) | is.na(l$longitude) | is.na(l$latitude)
    resGood <- gsw::gsw_Sstar_from_SP(SP = l$salinity[!bad], l$pressure[!bad], longitude = l$longitude[!bad], latitude = l$latitude[!bad])
    res <- rep(NA, nS)
    res[!bad] <- resGood
    if (!is.null(dim)) {
        dim(res) <- dim
    }
    res
}


#' Seawater Reference Salinity (GSW Formulation)
#'
#' Compute seawater Reference Salinity (SR), according to the GSW/TEOS-10
#' formulation with [gsw::gsw_SR_from_SP()] in the \CRANpkg{gsw} package.
#'
#' @param salinity either practical salinity or an `oce` object
#' that holds salinity in its data slot.
#'
#' @return Reference Salinity, SR, in \eqn{g/kg}{g/kg}.
#'
#' @author Dan Kelley
#'
#' @seealso For some objects, SR may also be recovered by indexing as e.g.
#' \code{ctd[["SR"]]}.
#'
#' @references McDougall, T.J. and P.M. Barker, 2011: Getting started with
#' TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp.,
#' SCOR/IAPSO WG127, ISBN 978-0-646-55621-5.
#'
#' @examples
#' SR <- swSR(35.0) # 35.16504
#'
#' @family functions that calculate seawater properties
swSR <- function(salinity) {
    if (missing(salinity)) {
        stop("must provide salinity")
    }
    if (inherits(salinity, "oce")) {
        x <- salinity # store this for clarity
        if (!"salinity" %in% names(x@data)) {
            stop("this oce object lacks salinity, so SR cannot be computed")
        }
        salinity <- x@data$salinity
    }
    gsw::gsw_SR_from_SP(salinity)
}