# 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, "\") START\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(...) START\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, "END 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::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 \dQuote{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(...) START\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, "END 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() START\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 START\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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.