#' Seawater Density according to UNESCO formula
#' @description UNESCO (1981) Tenth report of the joint panel on
#' oceanographic tables and standards. UNESCO Technical
#' Papers in Marine Science, Paris, 25p
#'
#' @param S Salinity in psu (which is +/- equal to ppm)
#' @param T Temperature in degrees
#' @param p pressure in Bar
#' @examples
#' rho(S=35,T=0.5,p=10)
#' rho(8,10) #Should be 1005.94659
#' @export
rho = function(S,T,p=0) ifelse(p==0, rho_p0(S,T), rho_p0(S,T) / (1 - p/K(S,T,p)))
#' Standard Mean Ocean Water (SMOW)
#' @param T Temperature in degrees Celsius
#' @export
rho_smow = function(T){
a0 = 999.842594
a1 = 6.793953 * 10^-2
a2 = -9.095290 * 10^-3
a3 = 1.001685 * 10^-4
a4 = -1.120083 * 10^-6
a5 = 6.536332 * 10^-9
a0 + a1*T + a2*T^2 + a3*T^3 + a4*T^4 + a5*T^5
}
#' Density (rho) at normal atmospheric pressure (p = 0)
#' @param S Salinity in psu (~ppm)
#' @param T Temperature in degrees Celsius
#' @export
rho_p0 = function(S,T){
b0 = 8.2449 * 10^-1
b1 = -4.0899 * 10^-3
b2 = 7.6438 * 10^-5
b3 = -8.2467 * 10^-7
b4 = 5.3875 * 10^-9
c0 = -5.7246 * 10^-3
c1 = 1.0227 * 10^-4
c2 = -1.6546 * 10^-6
d0 = 4.8314 * 10^-4
B1 = b0 + b1*T + b2*T^2 + b3*T^3 + b4*T^4
C1 = c0 + c1*T + c2*T^2
rho_smow(T) + B1*S + C1*S^1.5 + d0*S^2
}
#' compressibility at p = 0
#' @param S Salinity in psu (~ppm)
#' @param T Temperature in degrees Celsius
#' @export
K_p0 = function(S, T){
e0 = 19652.210
e1 = 148.4206
e2 = -2.327105
e3 = 1.360477 * 10^-2
e4 = -5.155288 * 10^-5
K_w = e0 + e1*T + e2*T^2 + e3*T^3 + e4*T^4
f0 = 54.6746
f1 = -0.603459
f2 = 1.099870 * 10^-2
f3 = -6.167 * 10^-5
F1 = f0 + f1*T + f2*T^2 + f3*T^3
g0 = 7.944 * 10^-2
g1 = 1.6483 * 10^-2
g2 = -5.3009 * 10^-4
G1 = g0 + g1*T + g2*T^2
K_w + F1*S + G1*S^1.5
}
#' compressibility K
#' @param S Salinity in psu (~ppm)
#' @param T Temperature in degrees Celsius
#' @param p pressure in bar
#' @export
K = function(S,T,p){
h0 = 3.2399
h1 = 1.43713 * 10^-3
h2 = 1.16092 * 10^-4
h3 = -5.77905 * 10^-7
i0 = 2.838 * 10^-3
i1 = -1.0981 * 10^-5
i2 = -1.6078 * 10^-6
j0 = 1.91075 * 10^-4
k0 = 8.50935 * 10^-5
k1 = -6.12293 * 10^-6
k2 = 5.2787 * 10^-8
m0 = -9.9348 * 10^-7
m1 = 2.0816 * 10^-8
m2 = 9.1697 * 10^-10
Bw = k0 + k1*T + k2*T^2
B2 = Bw + (m0 + m1*T + m2*T^2)*S
Aw = h0 + h1*T + h2*T^2 + h3*T^3
A1 = Aw + (i0 + i1*T + i2*T^2) * S + j0*S^1.5
K_p0(S,T) + A1*p + B2*p^2
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.