# -*- coding: utf-8 -*-
from __future__ import division
import numpy as np
from .constants import sfac, SSO, db2Pascal
from ..utilities import match_args_return, strip_mask, read_data
__all__ = ['deltaSA_atlas',
'enthalpy_SSO_0',
'enthalpy_SSO_0_p',
'entropy_part',
'entropy_part_zerop',
'Fdelta', # Seems to be new in V3.
'gibbs',
'gibbs_pt0_pt0',
'Hill_ratio_at_SP2',
'infunnel',
'interp_ref_cast',
'interp_SA_CT',
'SAAR', # New in V3.
'SA_from_SP_Baltic',
'specvol_SSO_0_p',
'SP_from_SA_Baltic']
h006 = -2.1078768810e-9
h007 = 2.8019291329e-10
class SA_table(object):
"""
TODO: Write docstring.
"""
# Central America barrier
x_ca = np.array([260.0, 272.59, 276.5, 278.65, 280.73, 295.217])
y_ca = np.array([19.55, 13.97, 9.6, 8.1, 9.33, 0.0])
def __init__(self, fname="gsw_data_v3_0.npz", max_p_fudge=10000,
min_frac=0):
self.fname = fname
self.max_p_fudge = max_p_fudge
self.min_frac = min_frac
data = read_data(fname)
self.lon = data.longs_ref.astype(np.float)
self.lat = data.lats_ref.astype(np.float)
self.p = data.p_ref # Depth levels
self.dlon = self.lon[1] - self.lon[0]
self.dlat = self.lat[1] - self.lat[0]
self.i_ca, self.j_ca = self.xy_to_ij(self.x_ca, self.y_ca)
# Make the order x, y, z:
# Start with deltaSA_ref (was delta_SA_ref in V2):
temp = data.deltaSA_ref.transpose((2, 1, 0)).copy()
self.dsa_ref = np.ma.masked_invalid(temp)
self.dsa_ref.data[self.dsa_ref.mask] = 0
# Now SAAR_ref, which did not exist in V2:
temp = data.SAAR_ref.transpose((2, 1, 0)).copy()
self.SAAR_ref = np.ma.masked_invalid(temp)
self.SAAR_ref.data[self.SAAR_ref.mask] = 0
def xy_to_ij(self, x, y):
"""
Convert from lat/lon to grid index coordinates,
without truncation or rounding.
"""
i = (x - self.lon[0]) / self.dlon
j = (y - self.lat[0]) / self.dlat
return i, j
def _central_america(self, di, dj, ii, jj, gm):
"""
Use a line running through Central America to zero
the goodmask for grid points in the Pacific forming
the grid box around input locations in the Atlantic,
and vice-versa.
"""
ix, jy = ii[0] + di, jj[0] + dj # Reconstruction: minor inefficiency.
inear = ((ix >= self.i_ca[0]) & (ix <= self.i_ca[-1])
& (jy >= self.j_ca[-1]) & (jy <= self.j_ca[0]))
if not inear.any():
return gm
inear_ind = inear.nonzero()[0]
ix = ix[inear]
jy = jy[inear]
ii = ii[:, inear]
jj = jj[:, inear]
jy_ca = np.interp(ix, self.i_ca, self.j_ca)
above = jy - jy_ca # > 0 if input point is above dividing line
# Intersections of left and right grid lines with dividing line
jleft_ca = np.interp(ii[0], self.i_ca, self.j_ca)
jright_ca = np.interp(ii[1], self.i_ca, self.j_ca)
jgrid_ca = [jleft_ca, jright_ca, jright_ca, jleft_ca]
# Zero the goodmask for grid points on opposite side of divider
for i in range(4):
opposite = (above * (jj[i] - jgrid_ca[i])) < 0
gm[i, inear_ind[opposite]] = 0
return gm
def xy_interp(self, di, dj, ii, jj, k):
"""
2-D interpolation, bilinear if all 4 surrounding
grid points are present, but treating missing points
as having the average value of the remaining grid
points. This matches the matlab V2 behavior.
"""
# Array of weights, CCW around the grid box
w = np.vstack(((1 - di) * (1 - dj), # lower left
di * (1 - dj), # lower right
di * dj, # upper right
(1 - di) * dj)) # upper left
gm = ~self.dsa.mask[ii, jj, k] # gm is "goodmask"
gm = self._central_america(di, dj, ii, jj, gm)
# Save a measure of real interpolation quality.
frac = (w * gm).sum(axis=0)
# Now loosen the interpolation, allowing a value to
# be calculated on a grid point that is masked.
# This matches the matlab gsw version 2 behavior.
jm_partial = gm.any(axis=0) & (~(gm.all(axis=0)))
# The weights of the unmasked points will be increased
# by the sum of the weights of the masked points divided
# by the number of unmasked points in the grid square.
# This is equivalent to setting the masked data values
# to the average of the unmasked values, and then
# unmasking, which is the matlab v2 implementation.
if jm_partial.any():
w_bad = w * (~gm)
w[:, jm_partial] += (w_bad[:, jm_partial].sum(axis=0) /
gm[:, jm_partial].sum(axis=0))
w *= gm
wsum = w.sum(axis=0)
valid = wsum > 0 # Only need to prevent division by zero here.
w[:, valid] /= wsum[valid]
w[:, ~valid] = 0
vv = self.dsa.data[ii, jj, k]
vv *= w
dsa = vv.sum(axis=0)
return dsa, frac
def _delta_SA(self, p, lon, lat):
"""
Table lookup engine--to be called only from SAAR or SA_ref.
"""
p = np.ma.masked_less(p, 0)
mask_in = np.ma.mask_or(np.ma.getmask(p), np.ma.getmask(lon))
mask_in = np.ma.mask_or(mask_in, np.ma.getmask(lat))
p, lon, lat = [np.ma.filled(a, 0).astype(float) for a in (p, lon, lat)]
p, lon, lat = np.broadcast_arrays(p, lon, lat)
if p.ndim > 1:
shape_in = p.shape
p, lon, lat = list(map(np.ravel, (p, lon, lat)))
reshaped = True
else:
reshaped = False
p_orig = p.copy() # Save for comparison to clipped p.
ix0, iy0 = self.xy_to_ij(lon, lat)
i0raw = np.floor(ix0).astype(int)
i0 = np.clip(i0raw, 0, len(self.lon) - 2)
di = ix0 - i0
j0raw = np.floor(iy0).astype(int)
j0 = np.clip(j0raw, 0, len(self.lat) - 2)
dj = iy0 - j0
# Start at lower left and go CCW; match order in _xy_interp.
ii = np.vstack((i0, i0 + 1, i0 + 1, i0))
jj = np.vstack((j0, j0, j0 + 1, j0 + 1))
k1 = np.searchsorted(self.p, p, side='right')
# Clip p and k1 at max p of grid cell.
kmax = (self.ndepth[ii, jj].max(axis=0) - 1)
mask_out = kmax.mask
kmax = kmax.filled(1)
clip_p = (p >= self.p[kmax])
p[clip_p] = self.p[kmax[clip_p]]
k1[clip_p] = kmax[clip_p]
k0 = k1 - 1
dsa0, frac0 = self.xy_interp(di, dj, ii, jj, k0)
dsa1, frac1 = self.xy_interp(di, dj, ii, jj, k1)
dp = np.diff(self.p)
pfrac = (p - self.p[k0]) / dp[k0]
delta_SA = dsa0 * (1 - pfrac) + dsa1 * pfrac
# Save intermediate results in case we are curious about
# them; the frac values are most likely to be useful.
# We won't bother to reshape them, though, and we may
# delete them later.
self.dsa0 = dsa0
self.frac0 = frac0
self.dsa1 = dsa1
self.frac1 = frac1
self.pfrac = pfrac
self.p_fudge = p_orig - p
# Editing options, in case we don't want to use
# values calculated from the wrong pressure, or from
# an incomplete SA table grid square.
# mask_out |= self.p_fudge > self.max_p_fudge
# mask_out |= self.frac1 < self.min_frac
# delta_SA = np.ma.array(delta_SA, mask=mask_out, copy=False)
# Later on, it is expected to be a masked array.
delta_SA = np.ma.array(delta_SA, copy=False)
if reshaped:
delta_SA.shape = shape_in
self.p_fudge.shape = shape_in
if mask_in is not np.ma.nomask:
delta_SA = np.ma.array(delta_SA, mask=mask_in, copy=False)
return delta_SA
def SAAR(self, p, lon, lat):
"""
Table lookup of salinity anomaly ratio, given pressure, lon, and lat.
"""
self.dsa = self.SAAR_ref
# In V2,
# ndepth from the file disagrees with the unmasked count from
# SAAR_ref in a few places; this should be fixed in the
# file, but for now we will simply calculate ndepth directly from
# SAAR_ref.
# TODO: check to see whether this discrepancy is also found in V3.
# TODO: check: do we even need to calculate ndepth? It doesn't
# appear to be used for anything.
# self.ndepth = np.ma.masked_invalid(data.ndepth_ref.T).astype(np.int8)
ndepth = self.dsa.count(axis=-1)
self.ndepth = np.ma.masked_equal(ndepth, 0)
return self._delta_SA(p, lon, lat)
def delta_SA_ref(self, p, lon, lat):
"""
Table lookup of salinity anomaly reference value, given pressure,
lon, and lat.
"""
self.dsa = self.dsa_ref
# See comment in previous method.
ndepth = self.dsa.count(axis=-1)
self.ndepth = np.ma.masked_equal(ndepth, 0)
return self._delta_SA(p, lon, lat)
[docs]def Fdelta(p, lon, lat):
"""
Fdelta from the Absolute Salinity Anomaly Ratio (SAAR):
Fdelta = (1 + r1)SAAR/(1 - r1*SAAR)
= (SA/Sstar) - 1
with r1 being the constant 0.35 based on the work of Pawlowicz et al.
(2011). Note that since SAAR is everywhere less than 0.001 in the global
ocean, Fdelta is only slightly different to 1.35*SAAR.
Parameters
----------
p : array_like
pressure [dbar]
lon : array_like
decimal degrees east (will be treated modulo 360)
lat : array_like
decimal degrees (+ve N, -ve S) [-90..+90]
Returns
-------
Fdelta : masked array; masked where no nearby ocean is found in data
Ratio of SA to Sstar, minus 1 [unitless]
Notes
-----
The mask is only set when the observation is well and truly on dry
land; often the warning flag is not set until one is several hundred
kilometers inland from the coast.
References
----------
.. [1] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
of seawater - 2010: Calculation and use of thermodynamic properties.
Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
UNESCO (English), 196 pp. Available from http://www.TEOS-10.org
See section 2.5 and appendices A.4 and A.5 of this TEOS-10 Manual.
.. [2] McDougall, T.J., D.R. Jackett, F.J. Millero, R. Pawlowicz and
P.M. Barker, 2012: A global algorithm for estimating Absolute Salinity.
Ocean Science, 8, 1123-1134.
http://www.ocean-sci.net/8/1123/2012/os-8-1123-2012.pdf
.. [3] Pawlowicz, R., D.G. Wright and F.J. Millero, 2011; The effects of
biogeochemical processes on oceanic conductivity/salinty/density
relationships and the characterization of real seawater. Ocean Science,
7, 363-387. http://www.ocean-sci.net/7/363/2011/os-7-363-2011.pdf
"""
r = 0.35
saar, _ = SAAR(p, lon, lat)
Fdelta = ((1 + r) * saar) / (1 - r * saar)
return Fdelta
@match_args_return
[docs]def Hill_ratio_at_SP2(t):
"""
USAGE:
Hill_ratio = Hill_ratio_at_SP2(t)
DESCRIPTION:
Calculates the Hill ratio, which is the adjustment needed to apply for
Practical Salinities smaller than 2. This ratio is defined at a
Practical Salinity = 2 and in-situ temperature, t using PSS-78. The Hill
ratio is the ratio of 2 to the output of the Hill et al. (1986) formula
for Practical Salinity at the conductivity ratio, Rt, at which Practical
Salinity on the PSS-78 scale is exactly 2.
INPUT:
t = in-situ temperature (ITS-90) [ deg C ]
OUTPUT:
Hill_ratio = Hill ratio at SP of 2 [ unitless ]
AUTHOR:
Trevor McDougall and Paul Barker
VERSION NUMBER: 3.0 (26th March, 2011)
"""
SP2 = 2 * np.ones_like(t)
a0 = 0.0080
a1 = -0.1692
a2 = 25.3851
a3 = 14.0941
a4 = -7.0261
a5 = 2.7081
b0 = 0.0005
b1 = -0.0056
b2 = -0.0066
b3 = -0.0375
b4 = 0.0636
b5 = -0.0144
g0 = 2.641463563366498e-1
g1 = 2.007883247811176e-4
g2 = -4.107694432853053e-6
g3 = 8.401670882091225e-8
g4 = -1.711392021989210e-9
g5 = 3.374193893377380e-11
g6 = -5.923731174730784e-13
g7 = 8.057771569962299e-15
g8 = -7.054313817447962e-17
g9 = 2.859992717347235e-19
k = 0.0162
t68 = t * 1.00024
ft68 = (t68 - 15) / (1 + k * (t68 - 15))
# -------------------------------------------------------------------------
# Find the initial estimates of Rtx (Rtx0) and of the derivative dSP_dRtx
# at SP = 2.
# -------------------------------------------------------------------------
Rtx0 = g0 + t68 * (g1 + t68 * (g2 + t68 * (g3 + t68 * (g4 + t68 * (g5
+ t68 * (g6 + t68 * (g7 + t68 * (g8 + t68 * g9))))))))
dSP_dRtx = (a1 + (2 * a2 + (3 * a3 + (4 * a4 + 5 * a5 * Rtx0) * Rtx0) *
Rtx0) * Rtx0 + ft68 * (b1 + (2 * b2 + (3 * b3 + (4 * b4 + 5 *
b5 * Rtx0) * Rtx0) * Rtx0) * Rtx0))
# -------------------------------------------------------------------------
# Begin a single modified Newton-Raphson iteration to find Rt at SP = 2.
# -------------------------------------------------------------------------
SP_est = (a0 + (a1 + (a2 + (a3 + (a4 + a5 * Rtx0) * Rtx0) * Rtx0) * Rtx0) *
Rtx0 + ft68 * (b0 + (b1 + (b2 + (b3 + (b4 + b5 * Rtx0) * Rtx0) *
Rtx0) * Rtx0) * Rtx0))
Rtx = Rtx0 - (SP_est - SP2) / dSP_dRtx
Rtxm = 0.5 * (Rtx + Rtx0)
dSP_dRtx = (a1 + (2 * a2 + (3 * a3 + (4 * a4 + 5 * a5 * Rtxm) * Rtxm) *
Rtxm) * Rtxm + ft68 * (b1 + (2 * b2 + (3 * b3 + (4 * b4 + 5 *
b5 * Rtxm) * Rtxm) * Rtxm) * Rtxm))
Rtx = Rtx0 - (SP_est - SP2) / dSP_dRtx
# This is the end of one full iteration of the modified Newton-Raphson
# iterative equation solver. The error in Rtx at this point is equivalent
# to an error in SP of 9e-16 psu.
x = 400 * Rtx * Rtx
sqrty = 10 * Rtx
part1 = 1 + x * (1.5 + x)
part2 = 1 + sqrty * (1 + sqrty * (1 + sqrty))
SP_Hill_raw_at_SP2 = SP2 - a0 / part1 - b0 * ft68 / part2
return 2. / SP_Hill_raw_at_SP2
@match_args_return
[docs]def SAAR(p, lon, lat):
"""
Absolute Salinity Anomaly Ratio (excluding the Baltic Sea).
Calculates the Absolute Salinity Anomaly Ratio, SAAR, in the open ocean
by spatially interpolating the global reference data set of SAAR to the
location of the seawater sample.
This function uses version 3.0 of the SAAR look up table.
Parameters
----------
p : array_like
pressure [dbar]
lon : array_like
decimal degrees east (will be treated modulo 360)
lat : array_like
decimal degrees (+ve N, -ve S) [-90..+90]
Returns
-------
SAAR : array
Absolute Salinity Anomaly Ratio [unitless]
in_ocean : boolean array
Notes
-----
The Absolute Salinity Anomaly Ratio in the Baltic Sea is evaluated
separately, since it is a function of Practical Salinity, not of space.
The present function returns a SAAR of zero for data in the Baltic Sea.
The correct way of calculating Absolute Salinity in the Baltic Sea is by
calling SA_from_SP.
The in_ocean flag is only set when the observation is well and truly on dry
land; often the warning flag is not set until one is several hundred
kilometers inland from the coast.
The algorithm is taken from the matlab implementation of the references,
but the numpy implementation here differs substantially from the
matlab implementation.
References
----------
.. [1] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
of seawater - 2010: Calculation and use of thermodynamic properties.
Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
UNESCO (English), 196 pp.
.. [2] McDougall, T.J., D.R. Jackett and F.J. Millero, 2010: An algorithm
for estimating Absolute Salinity in the global ocean. Submitted to
Ocean Science. A preliminary version is available at Ocean Sci.
Discuss., 6, 215-242.
http://www.ocean-sci-discuss.net/6/215/2009/osd-6-215-2009-print.pdf
"""
saar = SA_table().SAAR(p, lon, lat)
return saar, ~saar.mask
[docs]def SA_from_SP_Baltic(SP, lon, lat):
"""
Computes absolute salinity from practical in the Baltic Sea.
Parameters
----------
SP : array_like or masked array
Practical salinity (PSS-78)
lon, lat : array_like or masked arrays
geographical position
Returns
-------
SA : masked array, at least 1D
Absolute salinity [g/kg]
masked where inputs are masked or position outside the Baltic
"""
# Handle masked array input
input_mask = False
if np.ma.is_masked(SP):
input_mask = input_mask | SP.mask
if np.ma.is_masked(lon):
input_mask = input_mask | lon.mask
if np.ma.is_masked(lat):
input_mask = input_mask | lat.mask
SP, lon, lat = list(map(np.atleast_1d, (SP, lon, lat)))
SP, lon, lat = np.broadcast_arrays(SP, lon, lat)
inds_baltic = in_Baltic(lon, lat)
# SA_baltic = np.ma.masked_all(SP.shape, dtype=np.float)
all_nans = np.nan + np.zeros_like(SP)
SA_baltic = np.ma.MaskedArray(all_nans, mask=~inds_baltic)
if np.any(inds_baltic):
SA_baltic[inds_baltic] = (((SSO - 0.087) / 35) *
SP[inds_baltic] + 0.087)
SA_baltic.mask = SA_baltic.mask | input_mask | np.isnan(SP)
return SA_baltic
[docs]def SP_from_SA_Baltic(SA, lon, lat):
"""
Calculates Practical Salinity (SP) for the Baltic Sea, from a value
computed analytically from Absolute Salinity.
Parameters
----------
SA : array_like
Absolute salinity [g kg :sup:`-1`]
lon : array_like
decimal degrees east [0..+360]
lat : array_like
decimal degrees (+ve N, -ve S) [-90..+90]
Returns
-------
SP_baltic : array_like
salinity [psu (PSS-78)], unitless
See Also
--------
SP_from_SA, SP_from_Sstar
Notes
-----
This program will only produce Practical Salinity values for the Baltic
Sea.
Examples
--------
TODO
References
----------
.. [1] Feistel, R., S. Weinreben, H. Wolf, S. Seitz, P. Spitzer, B. Adel,
G. Nausch, B. Schneider and D. G. Wright, 2010c: Density and Absolute
Salinity of the Baltic Sea 2006-2009. Ocean Science, 6, 3-24.
http://www.ocean-sci.net/6/3/2010/os-6-3-2010.pdf
.. [2] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
of seawater - 2010: Calculation and use of thermodynamic properties.
Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
UNESCO (English), 196 pp.
.. [3] McDougall, T.J., D.R. Jackett and F.J. Millero, 2010: An algorithm
for estimating Absolute Salinity in the global ocean. Submitted to Ocean
Science. A preliminary version is available at Ocean Sci. Discuss.,
6, 215-242.
http://www.ocean-sci-discuss.net/6/215/2009/osd-6-215-2009-print.pdf
"""
SA, lon, lat = list(map(np.ma.masked_invalid, (SA, lon, lat)))
lon, lat, SA = np.broadcast_arrays(lon, lat, SA)
inds_baltic = in_Baltic(lon, lat)
if not inds_baltic.sum():
return None
SP_baltic = np.ma.masked_all(SA.shape, dtype=np.float)
SP_baltic[inds_baltic] = ((35 / (SSO - 0.087)) *
(SA[inds_baltic] - 0.087))
return SP_baltic
# FIXME: Check if this is still used and remove it.
def SP_from_SA_Baltic_old(SA, lon, lat):
"""
Calculates Practical Salinity (SP) for the Baltic Sea, from a value
computed analytically from Absolute Salinity.
Parameters
----------
SA : array_like
Absolute salinity [g kg :sup:`-1`]
lon : array_like
decimal degrees east [0..+360]
lat : array_like
decimal degrees (+ve N, -ve S) [-90..+90]
Returns
-------
SP_baltic : array_like
salinity [psu (PSS-78)], unitless
See Also
--------
SP_from_SA, SP_from_Sstar
Notes
-----
This program will only produce Practical Salinity values for the Baltic
Sea.
Examples
--------
TODO
References
----------
.. [1] Feistel, R., S. Weinreben, H. Wolf, S. Seitz, P. Spitzer, B. Adel,
G. Nausch, B. Schneider and D. G. Wright, 2010c: Density and Absolute
Salinity of the Baltic Sea 2006-2009. Ocean Science, 6, 3-24.
http://www.ocean-sci.net/6/3/2010/os-6-3-2010.pdf
.. [2] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
of seawater - 2010: Calculation and use of thermodynamic properties.
Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
UNESCO (English), 196 pp.
.. [3] McDougall, T.J., D.R. Jackett and F.J. Millero, 2010: An algorithm
for estimating Absolute Salinity in the global ocean. Submitted to Ocean
Science. A preliminary version is available at Ocean Sci. Discuss.,
6, 215-242.
http://www.ocean-sci-discuss.net/6/215/2009/osd-6-215-2009-print.pdf
"""
SA, lon, lat = list(map(np.ma.masked_invalid, (SA, lon, lat)))
lon, lat, SA = np.broadcast_arrays(lon, lat, SA, subok=True)
xb1, xb2, xb3 = 12.6, 7., 26.
xb1a, xb3a = 45., 26.
yb1, yb2, yb3 = 50., 59., 69.
inds_baltic = (xb2 < lon) & (lon < xb1a) & (yb1 < lat) & (lat < yb3)
if not inds_baltic.sum():
return None
SP_baltic = np.ma.masked_all(SA.shape, dtype=np.float)
xx_left = np.interp(lat[inds_baltic], [yb1, yb2, yb3], [xb1, xb2, xb3])
xx_right = np.interp(lat[inds_baltic], [yb1, yb3], [xb1a, xb3a])
inds_baltic1 = ((xx_left <= lon[inds_baltic]) &
(lon[inds_baltic] <= xx_right))
if not inds_baltic1.sum():
return None
SP_baltic[inds_baltic[inds_baltic1]] = ((35 / (SSO - 0.087)) *
(SA[inds_baltic[inds_baltic1]] -
0.087))
return SP_baltic
@match_args_return
[docs]def deltaSA_atlas(p, lon, lat):
"""
Absolute Salinity anomaly atlas value (excluding the Baltic Sea).
Calculates the Absolute Salinity anomaly atlas value, SA - SR,
in the open ocean by spatially interpolating the global reference
data set of deltaSA_atlas to the location of the seawater sample.
This function uses version 3.0 of the deltaSA_ref look up table.
Parameters
----------
p : array_like
sea pressure (absolute pressure - 10.1325 dbar) [dbar]
lon : array_like
decimal degrees east (will be treated modulo 360)
lat : array_like
decimal degrees (+ve N, -ve S) [-90..+90]
Returns
-------
deltaSA_atlas : masked array; masked where no nearby ocean is found in data
Absolute Salinity anomaly atlas value [g/kg]
Notes
-----
The Absolute Salinity anomaly atlas value in the Baltic Sea is
evaluated separately, since it is a function of Practical Salinity,
not of space. The present function returns a deltaSA_atlas of zero
for data in the Baltic Sea. The correct way of calculating Absolute
Salinity in the Baltic Sea is by calling SA_from_SP.
The mask is only set when the observation is well and truly on dry
land; often the warning flag is not set until one is several hundred
kilometers inland from the coast.
The algorithm is taken from the Matlab implementation of the references,
but the numpy implementation here differs substantially from the
Matlab implementation.
References
----------
.. [1] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
of seawater - 2010: Calculation and use of thermodynamic properties.
Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
UNESCO (English), 196 pp.
.. [2] McDougall, T.J., D.R. Jackett and F.J. Millero, 2010: An algorithm
for estimating Absolute Salinity in the global ocean. Submitted to Ocean
Science. A preliminary version is available at Ocean Sci. Discuss.,
6, 215-242.
http://www.ocean-sci-discuss.net/6/215/2009/osd-6-215-2009-print.pdf
"""
return SA_table().delta_SA_ref(p, lon, lat)
def enthalpy_SSO_0_CT25(p):
"""
Calculates enthalpy at the Standard Ocean Salinity (SSO) and at a
Conservative Temperature of zero degrees C (CT=0), as a function of
pressure (p [dbar]) or enthalpy_CT25(35.16504,0,p).
Parameters
----------
p : array_like
pressure [dbar]
Returns
-------
enthalpy_CT25 : array_like
enthalpy_CT25 at (SSO, CT = 0, p), 25-term equation.
[J kg :sup:`-1`]
Notes
-----
Uses a streamlined version of the 25-term CT version of the Gibbs function,
that is, a streamlined version of the code "enthalpy_CT25(SA,CT,p)"
"""
p = np.asanyarray(p)
mask = np.ma.getmask(p)
p = np.ma.filled(p, 0)
a0 = 1 + SSO * (2.0777716085618458e-3 + np.sqrt(SSO) *
3.4688210757917340e-6)
a1 = 6.8314629554123324e-6
b0 = 9.9984380290708214e2 + SSO * (2.8925731541277653e0 + SSO *
1.9457531751183059e-3)
b1 = 0.5 * (1.1930681818531748e-2 + SSO * 5.9355685925035653e-6)
b2 = -2.5943389807429039e-8
A = b1 - np.sqrt(b1 ** 2 - b0 * b2)
B = b1 + np.sqrt(b1 ** 2 - b0 * b2)
part = (a0 * b2 - a1 * b1) / (b2 * (B - A))
enthalpy_SSO_0_CT25 = (db2Pascal * ((a1 / (2 * b2)) *
np.log(1 + p * (2 * b1 + b2 * p) / b0) + part *
np.log(1 + (b2 * p * (B - A)) /
(A * (B + b2 * p)))))
return np.ma.array(enthalpy_SSO_0_CT25, mask=mask, copy=False)
[docs]def enthalpy_SSO_0(p):
"""
Enthalpy at SSO and CT(T=0) (75-term equation)
This function calculates enthalpy at the Standard Ocean Salinty, SSO,
and at a Conservative Temperature of zero degrees C, as a function of
pressure, p, in dbar, using a streamlined version of the 75-term
computationally-efficient expression for specific volume, that is, a
streamlined version of the code "enthalpy(SA,CT,p)".
Parameters
----------
p : array_like
pressure [dbar]
Returns
-------
enthalpy_SSO_0 : array_like
enthalpy at (SSO, CT=0, p)
[J kg :sup:`-1`]
Examples
--------
>>> import gsw
>>> p = np.array([10, 50, 125, 250, 600, 1000])
>>> gsw.library.enthalpy_SSO_0(p)
array([ 97.26388583, 486.27439853, 1215.47517122, 2430.24907325,
5827.90879421, 9704.32030926])
References
----------
.. [1] Roquet, F., G. Madec, T.J. McDougall, P.M. Barker, 2015: Accurate
polynomial expressions for the density and specifc volume of
seawater using the TEOS-10 standard. Ocean Modelling.
"""
z = p * 1e-4
dynamic_enthalpy_SSO_0_p = z * (9.726613854843870e-4 + z * (
-2.252956605630465e-5 + z * (2.376909655387404e-6 + z * (
-1.664294869986011e-7 + z * (-5.988108894465758e-9 + z * (
h006 + h007 * z))))))
return dynamic_enthalpy_SSO_0_p * db2Pascal * 1e4
# FIXME: Check if this is still used and remove it.
[docs]def enthalpy_SSO_0_p(p):
"""
This function calculates enthalpy at the Standard Ocean Salinty, SSO,
and at a Conservative Temperature of zero degrees C, as a function of
pressure, p, in dbar, using a streamlined version of the 48-term CT
version of the Gibbs function, that is, a streamlined version of the
code "enthalpy(SA,CT,p).
Examples
--------
>>> import gsw
>>> p = np.array([10, 50, 125, 250, 600, 1000])
>>> gsw.library.enthalpy_SSO_0_p(p)
array([ 97.26388276, 486.27439004, 1215.47518168, 2430.24919716,
5827.90973888, 9704.32296903])
References
----------
.. [1] McDougall T.J., P.M. Barker, R. Feistel and D.R. Jackett, 2013: A
computationally efficient 48-term expression for the density of seawater
in terms of Conservative Temperature, and related properties of
seawater. To be submitted to J. Atm. Ocean. Technol., xx, yyy-zzz.
"""
v01 = 9.998420897506056e+2
v05 = -6.698001071123802
v08 = -3.988822378968490e-2
v12 = -2.233269627352527e-2
v15 = -1.806789763745328e-4
v17 = -3.087032500374211e-7
v20 = 1.550932729220080e-10
v21 = 1.0
v26 = -7.521448093615448e-3
v31 = -3.303308871386421e-5
v36 = 5.419326551148740e-6
v37 = -2.742185394906099e-5
v41 = -1.105097577149576e-7
v43 = -1.119011592875110e-10
v47 = -1.200507748551599e-15
a0 = v21 + SSO * (v26 + v36 * SSO + v31 * np.sqrt(SSO))
a1 = v37 + v41 * SSO
a2 = v43
a3 = v47
b0 = v01 + SSO * (v05 + v08 * np.sqrt(SSO))
b1 = 0.5 * (v12 + v15 * SSO)
b2 = v17 + v20 * SSO
b1sq = b1 ** 2
sqrt_disc = np.sqrt(b1sq - b0 * b2)
N = a0 + (2 * a3 * b0 * b1 / b2 - a2 * b0) / b2
M = a1 + (4 * a3 * b1sq / b2 - a3 * b0 - 2 * a2 * b1) / b2
A = b1 - sqrt_disc
B = b1 + sqrt_disc
part = (N * b2 - M * b1) / (b2 * (B - A))
return (db2Pascal * (p * (a2 - 2 * a3 * b1 / b2 + 0.5 * a3 * p) / b2 +
(M / (2 * b2)) * np.log(1 + p * (2 * b1 + b2 * p) / b0) + part *
np.log(1 + (b2 * p * (B - A)) / (A * (B + b2 * p)))))
[docs]def entropy_part(SA, t, p):
"""
Calculates entropy, except that it does not evaluate any terms that are
functions of Absolute Salinity alone.
Parameters
----------
SA : array_like
Absolute salinity [g kg :sup:`-1`]
t : array_like
in situ temperature [:math:`^\circ` C (ITS-90)]
p : array_like
pressure [dbar]
Returns
-------
entropy_part : array_like
entropy minus the terms that due to SA alone
[J kg :sup:`-1` K :sup:`-1`]
Notes
-----
By not calculating these terms, which are a function only of Absolute
Salinity, several unnecessary computations are avoided (including saving
the computation of a natural logarithm). These terms are a necessary part
of entropy, but are not needed when calculating potential temperature from
in situ temperature.
"""
SA, t, p, mask = strip_mask(SA, t, p)
x2 = sfac * SA
x = np.sqrt(x2)
y = t * 0.025
z = p * 1e-4
g03 = (z * (-270.983805184062 + z * (776.153611613101 + z *
(-196.51255088122 + (28.9796526294175 - 2.13290083518327 * z) *
z))) + y * (-24715.571866078 + z * (2910.0729080936 + z *
(-1513.116771538718 + z * (546.959324647056 + z *
(-111.1208127634436 + 8.68841343834394 * z)))) + y *
(2210.2236124548363 + z * (-2017.52334943521 + z *
(1498.081172457456 + z * (-718.6359919632359 +
(146.4037555781616 - 4.9892131862671505 * z) * z))) + y *
(-592.743745734632 + z * (1591.873781627888 + z *
(-1207.261522487504 + (608.785486935364 - 105.4993508931208 * z) *
z)) + y * (290.12956292128547 + z * (-973.091553087975 + z *
(602.603274510125 + z * (-276.361526170076 + 32.40953340386105 *
z))) + y * (-113.90630790850321 + y * (21.35571525415769 -
67.41756835751434 * z) + z * (381.06836198507096 + z *
(-133.7383902842754 + 49.023632509086724 * z))))))))
g08 = (x2 * (z * (729.116529735046 + z * (-343.956902961561 + z *
(124.687671116248 + z * (-31.656964386073 + 7.04658803315449 *
z)))) + x * (x * (y * (-137.1145018408982 + y *
(148.10030845687618 + y * (-68.5590309679152 + 12.4848504784754 *
y))) - 22.6683558512829 * z) + z * (-175.292041186547 +
(83.1923927801819 - 29.483064349429 * z) * z) + y *
(-86.1329351956084 + z * (766.116132004952 + z *
(-108.3834525034224 + 51.2796974779828 * z)) + y *
(-30.0682112585625 - 1380.9597954037708 * z + y *
(3.50240264723578 + 938.26075044542 * z)))) + y *
(1760.062705994408 + y * (-675.802947790203 + y *
(365.7041791005036 + y * (-108.30162043765552 + 12.78101825083098 *
y) + z * (-1190.914967948748 + (298.904564555024 -
145.9491676006352 * z) * z)) + z * (2082.7344423998043 + z *
(-614.668925894709 + (340.685093521782 - 33.3848202979239 * z) *
z))) + z * (-1721.528607567954 + z * (674.819060538734 + z *
(-356.629112415276 + (88.4080716616 - 15.84003094423364 * z) *
z))))))
entropy_part = -(g03 + g08) * 0.025
return np.ma.array(entropy_part, mask=mask, copy=False)
[docs]def entropy_part_zerop(SA, pt0):
"""
Calculates entropy at a sea surface (p = 0 dbar), except that it does
not evaluate any terms that are functions of Absolute Salinity alone.
Parameters
----------
SA : array_like
Absolute salinity [g kg :sup:`-1`]
pt0 : array_like
potential temperature relative to 0 dbar [:math:`^\circ` C (ITS-90)]
Returns
-------
entropy_part_zerop : array_like
[J kg :sup:`-1` K :sup:`-1`]
Notes
-----
By not calculating these terms, which are a function only of Absolute
Salinity, several unnecessary computations are avoided (including saving
the computation of a natural logarithm). These terms are a necessary part
of entropy, but are not needed when calculating potential temperature from
in situ temperature.
"""
SA, pt0, mask = strip_mask(SA, pt0)
x2 = sfac * SA
x = np.sqrt(x2)
y = pt0 * 0.025
g03 = y * (-24715.571866078 + y * (2210.2236124548363 +
y * (-592.743745734632 + y * (290.12956292128547 +
y * (-113.90630790850321 + y * 21.35571525415769)))))
g08 = x2 * (x * (x * (y * (-137.1145018408982 + y *
(148.10030845687618 +
y * (-68.5590309679152 + 12.4848504784754 * y)))) +
y * (-86.1329351956084 + y * (-30.0682112585625 + y *
3.50240264723578))) +
y * (1760.062705994408 + y * (-675.802947790203 +
y * (365.7041791005036 + y * (-108.30162043765552 +
12.78101825083098 * y)))))
entropy_part_zerop = -(g03 + g08) * 0.025
return np.ma.array(entropy_part_zerop, mask=mask, copy=False)
# FIXME: Check if this is still used and remove it.
[docs]def gibbs(ns, nt, npr, SA, t, p):
"""
Calculates specific Gibbs energy and its derivatives up to order 2 for
seawater.
The Gibbs function approach allows the calculation of internal energy,
entropy, enthalpy, potential enthalpy and the chemical potentials of
seawater as well as the freezing temperature, and the latent heats of
freezing and of evaporation. These quantities were not available from
EOS-80 but are essential for the accurate accounting of heat in the ocean
and for the consistent and accurate treatment of air-sea and ice-sea heat
fluxes.
Parameters
----------
ns : int
order of SA derivative [0, 1 or 2 ]
nt : int
order of t derivative [0, 1 or 2 ]
npr : int
order of p derivative [0, 1 or 2 ]
SA : array_like
Absolute salinity [g kg :sup:`-1`]
t : array_like
in situ temperature [:math:`^\circ` C (ITS-90)]
p : array_like
pressure [dbar]
Returns
-------
gibbs : array_like
Specific Gibbs energy or its derivatives.
Gibbs energy (ns=nt=npr=0) has units of:
[J kg :sup:`-1`]
Absolute Salinity derivatives are output in units of:
[(J kg :sup:`-1`) (g kg :sup:`-1`) :sup:`-ns`]
Temperature derivatives are output in units of:
[(J kg :sup:`-1`) K :sup:`-nt`]
Pressure derivatives are output in units of:
[(J kg :sup:`-1`) Pa :sup:`-npr`]
The mixed derivatives are output in units of:
[(J kg :sup:`-1`) (g kg :sup:`-1`) :sup:`-ns` K :sup:`-nt`
Pa :sup:`-npr`]
Notes
-----
The Gibbs function for seawater is that of TEOS-10 (IOC et al., 2010),
being the sum of IAPWS-08 for the saline part and IAPWS-09 for the pure
water part. These IAPWS releases are the officially blessed IAPWS
descriptions of Feistel (2008) and the pure water part of Feistel (2003).
Absolute Salinity, SA, in all of the GSW routines is expressed on the
Reference-Composition Salinity Scale of 2008 (RCSS-08) of Millero et al.
(2008).
The derivatives are taken with respect to pressure in Pa, not withstanding
that the pressure input into this routine is in dbar.
References
----------
.. [1] Feistel, R., 2003: A new extended Gibbs thermodynamic potential of
seawater Progr. Oceanogr., 58, 43-114.
.. [2] Feistel, R., 2008: A Gibbs function for seawater thermodynamics
for -6 to 80 :math:`^\circ` C and salinity up to 120 g kg :sup:`-1`,
Deep-Sea Res. I, 55, 1639-1671.
.. [3] IAPWS, 2008: Release on the IAPWS Formulation 2008 for the
Thermodynamic Properties of Seawater. The International Association for
the Properties of Water and Steam. Berlin, Germany, September 2008,
available from http://www.iapws.org. This Release is referred to as
IAPWS-08.
.. [4] IAPWS, 2009: Supplementary Release on a Computationally Efficient
Thermodynamic Formulation for Liquid Water for Oceanographic Use. The
International Association for the Properties of Water and Steam.
Doorwerth, The Netherlands, September 2009, available from
http://www.iapws.org. This Release is referred to as IAPWS-09.
.. [5] IOC, SCOR and IAPSO, 2010: The international thermodynamic equation
of seawater - 2010: Calculation and use of thermodynamic properties.
Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
UNESCO (English), 196 pp. See section 2.6 and appendices A.6, G and H.
.. [6] Millero, F. J., R. Feistel, D. G. Wright, and T. J. McDougall, 2008:
The composition of Standard Seawater and the definition of the
Reference-Composition Salinity Scale, Deep-Sea Res. I, 55, 50-72.
"""
SA, t, p = np.asanyarray(SA), np.asanyarray(t), np.asanyarray(p)
SA = np.atleast_1d(SA)
nonzero_SA = np.any(SA > 0)
_SA = SA
_t = t
_p = p
SA = np.ma.filled(SA, 0)
t = np.ma.filled(t, 20)
p = np.ma.filled(p, 10)
SA, t, p = np.broadcast_arrays(SA, t, p)
gibbs = np.zeros(SA.shape, dtype=np.float) # Use if all_masked is True
all_masked = False
# Ensure a full mask, so we can set elements if necessary.
mask = np.ma.mask_or(np.ma.getmaskarray(_SA), np.ma.getmask(_t))
mask = np.ma.mask_or(mask, np.ma.getmask(_p))
mask = np.ma.mask_or(mask, SA < 0)
ipos = (SA > 0)
# inpos = ~ipos # FIXME: Assigned but never used.
if np.all(ipos):
ipos = slice(None) # More efficient for usual case.
x2 = sfac * SA
x = np.sqrt(x2)
y = t * 0.025
z = p * 1e-4 # The input pressure (p) is sea pressure in units of dbar.
if (ns == 0) & (nt == 0) & (npr == 0):
g03 = (101.342743139674 + z * (100015.695367145 +
z * (-2544.5765420363 + z * (284.517778446287 +
z * (-33.3146754253611 + (4.20263108803084 - 0.546428511471039 * z)
* z)))) +
y * (5.90578347909402 + z * (-270.983805184062 +
z * (776.153611613101 + z * (-196.51255088122 + (28.9796526294175 -
2.13290083518327 * z) * z))) +
y * (-12357.785933039 + z * (1455.0364540468 +
z * (-756.558385769359 + z * (273.479662323528 + z *
(-55.5604063817218 + 4.34420671917197 * z)))) +
y * (736.741204151612 + z * (-672.50778314507 +
z * (499.360390819152 + z * (-239.545330654412 + (48.8012518593872 -
1.66307106208905 * z) * z))) +
y * (-148.185936433658 + z * (397.968445406972 +
z * (-301.815380621876 + (152.196371733841 - 26.3748377232802 * z) *
z)) +
y * (58.0259125842571 + z * (-194.618310617595 +
z * (120.520654902025 + z * (-55.2723052340152 +
6.48190668077221 * z))) +
y * (-18.9843846514172 + y * (3.05081646487967 -
9.63108119393062 * z) +
z * (63.5113936641785 + z * (-22.2897317140459 +
8.17060541818112 * z)))))))))
if nonzero_SA:
g08 = x2 * (1416.27648484197 + z * (-3310.49154044839 +
z * (384.794152978599 + z * (-96.5324320107458 +
(15.8408172766824 - 2.62480156590992 * z) * z))) +
x * (-2432.14662381794 + x * (2025.80115603697 +
y * (543.835333000098 + y * (-68.5572509204491 +
y * (49.3667694856254 + y * (-17.1397577419788 +
2.49697009569508 * y))) - 22.6683558512829 * z) +
x * (-1091.66841042967 - 196.028306689776 * y +
x * (374.60123787784 - 48.5891069025409 * x +
36.7571622995805 * y) + 36.0284195611086 * z) +
z * (-54.7919133532887 + (-4.08193978912261 -
30.1755111971161 * z) * z)) +
z * (199.459603073901 + z * (-52.2940909281335 +
(68.0444942726459 - 3.41251932441282 * z) * z)) +
y * (-493.407510141682 + z * (-175.292041186547 +
(83.1923927801819 - 29.483064349429 * z) * z) +
y * (-43.0664675978042 + z * (383.058066002476 + z *
(-54.1917262517112 + 25.6398487389914 * z)) +
y * (-10.0227370861875 - 460.319931801257 * z + y *
(0.875600661808945 + 234.565187611355 * z))))) +
y * (168.072408311545 + z * (729.116529735046 +
z * (-343.956902961561 + z * (124.687671116248 + z *
(-31.656964386073 + 7.04658803315449 * z)))) +
y * (880.031352997204 + y * (-225.267649263401 +
y * (91.4260447751259 + y * (-21.6603240875311 +
2.13016970847183 * y) +
z * (-297.728741987187 + (74.726141138756 -
36.4872919001588 * z) * z)) +
z * (694.244814133268 + z * (-204.889641964903 +
(113.561697840594 - 11.1282734326413 * z) * z))) +
z * (-860.764303783977 + z * (337.409530269367 +
z * (-178.314556207638 + (44.2040358308 -
7.92001547211682 * z) * z))))))
g08[ipos] += x2[ipos] * (5812.81456626732 + 851.226734946706 *
y[ipos]) * np.log(x[ipos])
else:
g08 = 0
gibbs = g03 + g08
elif (ns == 1) & (nt == 0) & (npr == 0):
if nonzero_SA:
g08 = (8645.36753595126 + z * (-6620.98308089678 +
z * (769.588305957198 + z * (-193.0648640214916 +
(31.6816345533648 - 5.24960313181984 * z) * z))) +
x * (-7296.43987145382 + x * (8103.20462414788 +
y * (2175.341332000392 + y * (-274.2290036817964 +
y * (197.4670779425016 + y * (-68.5590309679152 +
9.98788038278032 * y))) - 90.6734234051316 * z) +
x * (-5458.34205214835 - 980.14153344888 * y +
x * (2247.60742726704 - 340.1237483177863 * x +
220.542973797483 * y) + 180.142097805543 * z) +
z * (-219.1676534131548 + (-16.32775915649044 -
120.7020447884644 * z) * z)) +
z * (598.378809221703 + z * (-156.8822727844005 +
(204.1334828179377 - 10.23755797323846 * z) * z)) +
y * (-1480.222530425046 + z * (-525.876123559641 +
(249.57717834054571 - 88.449193048287 * z) * z) +
y * (-129.1994027934126 + z * (1149.174198007428 +
z * (-162.5751787551336 + 76.9195462169742 * z)) +
y * (-30.0682112585625 - 1380.9597954037708 * z + y *
(2.626801985426835 + 703.695562834065 * z))))) +
y * (1187.3715515697959 + z * (1458.233059470092 +
z * (-687.913805923122 + z * (249.375342232496 + z *
(-63.313928772146 + 14.09317606630898 * z)))) +
y * (1760.062705994408 + y * (-450.535298526802 +
y * (182.8520895502518 + y * (-43.3206481750622 +
4.26033941694366 * y) +
z * (-595.457483974374 + (149.452282277512 -
72.9745838003176 * z) * z)) +
z * (1388.489628266536 + z * (-409.779283929806 +
(227.123395681188 - 22.2565468652826 * z) * z))) +
z * (-1721.528607567954 + z * (674.819060538734 +
z * (-356.629112415276 + (88.4080716616 -
15.84003094423364 * z) * z))))))
g08[ipos] = g08[ipos] + (11625.62913253464 + 1702.453469893412 *
y[ipos]) * np.log(x[ipos])
gibbs = 0.5 * sfac * g08
else:
all_masked = True
elif (ns == 0) & (nt == 1) & (npr == 0):
g03 = (5.90578347909402 + z * (-270.983805184062 +
z * (776.153611613101 + z * (-196.51255088122 +
(28.9796526294175 - 2.13290083518327 * z) * z))) +
y * (-24715.571866078 + z * (2910.0729080936 +
z * (-1513.116771538718 + z * (546.959324647056 + z *
(-111.1208127634436 + 8.68841343834394 * z)))) +
y * (2210.2236124548363 + z * (-2017.52334943521 +
z * (1498.081172457456 + z * (-718.6359919632359 +
(146.4037555781616 - 4.9892131862671505 * z) * z))) +
y * (-592.743745734632 + z * (1591.873781627888 +
z * (-1207.261522487504 + (608.785486935364 -
105.4993508931208 * z) * z)) +
y * (290.12956292128547 + z * (-973.091553087975 +
z * (602.603274510125 + z * (-276.361526170076 +
32.40953340386105 * z))) +
y * (-113.90630790850321 + y * (21.35571525415769 -
67.41756835751434 * z) +
z * (381.06836198507096 + z * (-133.7383902842754 +
49.023632509086724 * z))))))))
if nonzero_SA:
g08 = x2 * (168.072408311545 + z * (729.116529735046 +
z * (-343.956902961561 + z * (124.687671116248 + z *
(-31.656964386073 + 7.04658803315449 * z)))) +
x * (-493.407510141682 + x * (543.835333000098 + x *
(-196.028306689776 + 36.7571622995805 * x) +
y * (-137.1145018408982 + y * (148.10030845687618 + y *
(-68.5590309679152 + 12.4848504784754 * y))) -
22.6683558512829 * z) + z * (-175.292041186547 +
(83.1923927801819 - 29.483064349429 * z) * z) +
y * (-86.1329351956084 + z * (766.116132004952 + z *
(-108.3834525034224 + 51.2796974779828 * z)) +
y * (-30.0682112585625 - 1380.9597954037708 * z + y *
(3.50240264723578 + 938.26075044542 * z)))) +
y * (1760.062705994408 + y * (-675.802947790203 +
y * (365.7041791005036 + y * (-108.30162043765552 +
12.78101825083098 * y) +
z * (-1190.914967948748 + (298.904564555024 -
145.9491676006352 * z) * z)) +
z * (2082.7344423998043 + z * (-614.668925894709 +
(340.685093521782 - 33.3848202979239 * z) * z))) +
z * (-1721.528607567954 + z * (674.819060538734 +
z * (-356.629112415276 + (88.4080716616 -
15.84003094423364 * z) * z)))))
g08[ipos] += 851.226734946706 * x2[ipos] * np.log(x[ipos])
gibbs = (g03 + g08) * 0.025
else:
gibbs = g03
elif (ns == 0) & (nt == 0) & (npr == 1):
g03 = (100015.695367145 + z * (-5089.1530840726 +
z * (853.5533353388611 + z * (-133.2587017014444 +
(21.0131554401542 - 3.278571068826234 * z) * z))) +
y * (-270.983805184062 + z * (1552.307223226202 +
z * (-589.53765264366 + (115.91861051767 -
10.664504175916349 * z) * z)) +
y * (1455.0364540468 + z * (-1513.116771538718 +
z * (820.438986970584 + z * (-222.2416255268872 +
21.72103359585985 * z))) +
y * (-672.50778314507 + z * (998.720781638304 +
z * (-718.6359919632359 + (195.2050074375488 -
8.31535531044525 * z) * z)) +
y * (397.968445406972 + z * (-603.630761243752 +
(456.589115201523 - 105.4993508931208 * z) * z) +
y * (-194.618310617595 + y * (63.5113936641785 -
9.63108119393062 * y +
z * (-44.5794634280918 + 24.511816254543362 * z)) +
z * (241.04130980405 + z * (-165.8169157020456 +
25.92762672308884 * z))))))))
if nonzero_SA:
g08 = x2 * (-3310.49154044839 + z * (769.588305957198 +
z * (-289.5972960322374 + (63.3632691067296 -
13.1240078295496 * z) * z)) +
x * (199.459603073901 + x * (-54.7919133532887 +
36.0284195611086 * x - 22.6683558512829 * y +
(-8.16387957824522 - 90.52653359134831 * z) * z) +
z * (-104.588181856267 + (204.1334828179377 -
13.65007729765128 * z) * z) +
y * (-175.292041186547 + (166.3847855603638 -
88.449193048287 * z) * z +
y * (383.058066002476 + y * (-460.319931801257 +
234.565187611355 * y) +
z * (-108.3834525034224 + 76.9195462169742 * z)))) +
y * (729.116529735046 + z * (-687.913805923122 +
z * (374.063013348744 + z * (-126.627857544292 +
35.23294016577245 * z))) +
y * (-860.764303783977 + y * (694.244814133268 +
y * (-297.728741987187 + (149.452282277512 -
109.46187570047641 * z) * z) +
z * (-409.779283929806 + (340.685093521782 -
44.5130937305652 * z) * z)) +
z * (674.819060538734 + z * (-534.943668622914 +
(176.8161433232 - 39.600077360584095 * z) * z)))))
else:
g08 = 0
# Pressure derivative of the Gibbs function
# in units of (J kg :sup:`-1`) (Pa :sup:`-1`) = m :sup:`3` kg :sup:`-1`
gibbs = (g03 + g08) * 1e-8
elif (ns == 1) & (nt == 1) & (npr == 0):
if nonzero_SA:
g08 = (1187.3715515697959 + z * (1458.233059470092 +
z * (-687.913805923122 + z * (249.375342232496 + z *
(-63.313928772146 + 14.09317606630898 * z)))) +
x * (-1480.222530425046 + x * (2175.341332000392 + x *
(-980.14153344888 + 220.542973797483 * x) +
y * (-548.4580073635929 + y * (592.4012338275047 + y *
(-274.2361238716608 + 49.9394019139016 * y))) -
90.6734234051316 * z) + z * (-525.876123559641 +
(249.57717834054571 - 88.449193048287 * z) * z) +
y * (-258.3988055868252 + z * (2298.348396014856 +
z * (-325.1503575102672 + 153.8390924339484 * z)) +
y * (-90.2046337756875 - 4142.8793862113125 * z + y *
(10.50720794170734 + 2814.78225133626 * z)))) +
y * (3520.125411988816 + y * (-1351.605895580406 +
y * (731.4083582010072 + y * (-216.60324087531103 +
25.56203650166196 * y) +
z * (-2381.829935897496 + (597.809129110048 -
291.8983352012704 * z) * z)) +
z * (4165.4688847996085 + z * (-1229.337851789418 +
(681.370187043564 - 66.7696405958478 * z) * z))) +
z * (-3443.057215135908 + z * (1349.638121077468 +
z * (-713.258224830552 + (176.8161433232 -
31.68006188846728 * z) * z)))))
g08[ipos] = g08[ipos] + 1702.453469893412 * np.log(x[ipos])
gibbs = 0.5 * sfac * 0.025 * g08
# FIXME: commented by FF, g110 without nan did not pass
# mask[inpos] = True
else:
all_masked = True
elif (ns == 1) & (nt == 0) & (npr == 1):
g08 = (-6620.98308089678 + z * (1539.176611914396 +
z * (-579.1945920644748 + (126.7265382134592 -
26.2480156590992 * z) * z)) +
x * (598.378809221703 + x * (-219.1676534131548 +
180.142097805543 * x - 90.6734234051316 * y +
(-32.65551831298088 - 362.10613436539325 * z) * z) +
z * (-313.764545568801 + (612.4004484538132 -
40.95023189295384 * z) * z) +
y * (-525.876123559641 + (499.15435668109143 -
265.347579144861 * z) * z +
y * (1149.174198007428 + y * (-1380.9597954037708 +
703.695562834065 * y) +
z * (-325.1503575102672 + 230.7586386509226 * z)))) +
y * (1458.233059470092 + z * (-1375.827611846244 +
z * (748.126026697488 + z * (-253.255715088584 +
70.4658803315449 * z))) +
y * (-1721.528607567954 + y * (1388.489628266536 +
y * (-595.457483974374 + (298.904564555024 -
218.92375140095282 * z) * z) +
z * (-819.558567859612 + (681.370187043564 -
89.0261874611304 * z) * z)) +
z * (1349.638121077468 + z * (-1069.887337245828 +
(353.6322866464 - 79.20015472116819 * z) * z)))))
# Derivative of the Gibbs function is in units of
# (m :sup:`3` kg :sup:`-1`) / (g kg :sup:`-1`) = m :sup:`3` g :sup:`-1`
# that is, it is the derivative of specific volume with respect to
# Absolute Salinity measured in g kg :sup:`-1`
gibbs = g08 * sfac * 0.5e-8
elif (ns == 0) & (nt == 1) & (npr == 1):
g03 = (-270.983805184062 + z * (1552.307223226202 + z *
(-589.53765264366 +
(115.91861051767 - 10.664504175916349 * z) * z)) +
y * (2910.0729080936 + z * (-3026.233543077436 +
z * (1640.877973941168 + z * (-444.4832510537744 +
43.4420671917197 * z))) +
y * (-2017.52334943521 + z * (2996.162344914912 +
z * (-2155.907975889708 + (585.6150223126464 -
24.946065931335752 * z) * z)) +
y * (1591.873781627888 + z * (-2414.523044975008 +
(1826.356460806092 - 421.9974035724832 * z) * z) +
y * (-973.091553087975 + z * (1205.20654902025 + z *
(-829.084578510228 + 129.6381336154442 * z)) +
y * (381.06836198507096 - 67.41756835751434 * y + z *
(-267.4767805685508 + 147.07089752726017 * z)))))))
if nonzero_SA:
g08 = x2 * (729.116529735046 + z * (-687.913805923122 +
z * (374.063013348744 + z * (-126.627857544292 +
35.23294016577245 * z))) +
x * (-175.292041186547 - 22.6683558512829 * x +
(166.3847855603638 - 88.449193048287 * z) * z +
y * (766.116132004952 + y * (-1380.9597954037708 +
938.26075044542 * y) +
z * (-216.7669050068448 + 153.8390924339484 * z))) +
y * (-1721.528607567954 + y * (2082.7344423998043 +
y * (-1190.914967948748 + (597.809129110048 -
437.84750280190565 * z) * z) +
z * (-1229.337851789418 + (1022.055280565346 -
133.5392811916956 * z) * z)) +
z * (1349.638121077468 + z * (-1069.887337245828 +
(353.6322866464 - 79.20015472116819 * z) * z))))
else:
g08 = 0
# Derivative of the Gibbs function is in units of (m :sup:`3` (K kg))
# that is, the pressure of the derivative in Pa.
gibbs = (g03 + g08) * 2.5e-10
elif (ns == 2) & (nt == 0) & (npr == 0):
g08 = 2.0 * (8103.20462414788 +
y * (2175.341332000392 + y * (-274.2290036817964 +
y * (197.4670779425016 + y * (-68.5590309679152 +
9.98788038278032 * y))) - 90.6734234051316 * z) +
1.5 * x * (-5458.34205214835 - 980.14153344888 * y +
(4.0 / 3.0) * x * (2247.60742726704 - 340.1237483177863 * 1.25 *
x + 220.542973797483 * y) +
180.142097805543 * z) +
z * (-219.1676534131548 + (-16.32775915649044 -
120.7020447884644 * z) * z))
if nonzero_SA:
tmp = ((-7296.43987145382 + z * (598.378809221703 +
z * (-156.8822727844005 + (204.1334828179377 -
10.23755797323846 * z) * z)) +
y * (-1480.222530425046 + z * (-525.876123559641 +
(249.57717834054571 - 88.449193048287 * z) * z) +
y * (-129.1994027934126 + z * (1149.174198007428 +
z * (-162.5751787551336 + 76.9195462169742 * z)) +
y * (-30.0682112585625 - 1380.9597954037708 * z +
y * (2.626801985426835 + 703.695562834065 * z))))) / x +
(11625.62913253464 + 1702.453469893412 * y) / x2)
g08[ipos] += tmp[ipos]
gibbs = 0.25 * sfac ** 2 * g08
elif (ns == 0) & (nt == 2) & (npr == 0):
g03 = (-24715.571866078 + z * (2910.0729080936 + z *
(-1513.116771538718 + z * (546.959324647056 + z *
(-111.1208127634436 + 8.68841343834394 * z)))) +
y * (4420.4472249096725 + z * (-4035.04669887042 +
z * (2996.162344914912 + z * (-1437.2719839264719 +
(292.8075111563232 - 9.978426372534301 * z) * z))) +
y * (-1778.231237203896 + z * (4775.621344883664 +
z * (-3621.784567462512 + (1826.356460806092 -
316.49805267936244 * z) * z)) +
y * (1160.5182516851419 + z * (-3892.3662123519 +
z * (2410.4130980405 + z * (-1105.446104680304 +
129.6381336154442 * z))) +
y * (-569.531539542516 + y * (128.13429152494615 -
404.50541014508605 * z) +
z * (1905.341809925355 + z * (-668.691951421377 +
245.11816254543362 * z)))))))
if nonzero_SA:
g08 = x2 * (1760.062705994408 + x * (-86.1329351956084 +
x * (-137.1145018408982 + y * (296.20061691375236 +
y * (-205.67709290374563 + 49.9394019139016 * y))) +
z * (766.116132004952 + z * (-108.3834525034224 +
51.2796974779828 * z)) +
y * (-60.136422517125 - 2761.9195908075417 * z +
y * (10.50720794170734 + 2814.78225133626 * z))) +
y * (-1351.605895580406 + y * (1097.1125373015109 +
y * (-433.20648175062206 + 63.905091254154904 * y) +
z * (-3572.7449038462437 + (896.713693665072 -
437.84750280190565 * z) * z)) +
z * (4165.4688847996085 + z * (-1229.337851789418 +
(681.370187043564 - 66.7696405958478 * z) * z))) +
z * (-1721.528607567954 + z * (674.819060538734 +
z * (-356.629112415276 + (88.4080716616 -
15.84003094423364 * z) * z))))
else:
g08 = 0
gibbs = (g03 + g08) * 0.000625
elif (ns == 0) & (nt == 0) & (npr == 2):
g03 = (-5089.1530840726 + z * (1707.1066706777221 +
z * (-399.7761051043332 + (84.0526217606168 -
16.39285534413117 * z) * z)) +
y * (1552.307223226202 + z * (-1179.07530528732 +
(347.75583155301 - 42.658016703665396 * z) * z) +
y * (-1513.116771538718 + z * (1640.877973941168 +
z * (-666.7248765806615 + 86.8841343834394 * z)) +
y * (998.720781638304 + z * (-1437.2719839264719 +
(585.6150223126464 - 33.261421241781 * z) * z) +
y * (-603.630761243752 + (913.178230403046 -
316.49805267936244 * z) * z +
y * (241.04130980405 + y * (-44.5794634280918 +
49.023632509086724 * z) +
z * (-331.6338314040912 + 77.78288016926652 * z)))))))
if nonzero_SA:
g08 = x2 * (769.588305957198 + z * (-579.1945920644748 +
(190.08980732018878 - 52.4960313181984 * z) * z) +
x * (-104.588181856267 + x * (-8.16387957824522 -
181.05306718269662 * z) +
(408.2669656358754 - 40.95023189295384 * z) * z +
y * (166.3847855603638 - 176.898386096574 * z + y *
(-108.3834525034224 + 153.8390924339484 * z))) +
y * (-687.913805923122 + z * (748.126026697488 +
z * (-379.883572632876 + 140.9317606630898 * z)) +
y * (674.819060538734 + z * (-1069.887337245828 +
(530.4484299696 - 158.40030944233638 * z) * z) +
y * (-409.779283929806 + y * (149.452282277512 -
218.92375140095282 * z) +
(681.370187043564 - 133.5392811916956 * z) * z))))
else:
g08 = 0
# Second derivative of the Gibbs function with respect to pressure,
# measured in Pa; units of (J kg :sup:`-1`) (Pa :sup:`-2`).
gibbs = (g03 + g08) * 1e-16
else:
raise ValueError('Illegal derivative of the Gibbs function')
gibbs = np.ma.array(gibbs, mask=mask, copy=False)
# BÅ: Code below is not needed? EF: it is needed.
if all_masked:
gibbs[:] = np.ma.masked
# Do not allow zero salinity with salinity derivatives
if ns > 0:
gibbs = np.ma.masked_where(SA == 0, gibbs)
return gibbs
[docs]def gibbs_pt0_pt0(SA, pt0):
"""
Calculates the second derivative of the specific Gibbs function with
respect to temperature at zero sea pressure or _gibbs(0,2,0,SA,t,0).
Parameters
----------
SA : array_like
Absolute salinity [g kg :sup:`-1`]
pt0 : array_like
potential temperature relative to 0 dbar [:math:`^\circ` C (ITS-90)]
Returns
-------
gibbs_pt0_pt0 : array_like
TODO: write the eq. for the second derivative of the
specific Gibbs function. FIXME: [units]
Notes
-----
This library function is called by both "pt_from_CT(SA,CT)"
and "pt0_from_t(SA,t,p)".
"""
SA, pt0, mask = strip_mask(SA, pt0)
x2 = sfac * SA
x = np.sqrt(x2)
y = pt0 * 0.025
g03 = (-24715.571866078 +
y * (4420.4472249096725 +
y * (-1778.231237203896 +
y * (1160.5182516851419 +
y * (-569.531539542516 + y * 128.13429152494615)))))
g08 = (x2 * (1760.062705994408 + x * (-86.1329351956084 +
x * (-137.1145018408982 + y * (296.20061691375236 +
y * (-205.67709290374563 + 49.9394019139016 * y))) +
y * (-60.136422517125 + y * 10.50720794170734)) +
y * (-1351.605895580406 + y * (1097.1125373015109 +
y * (-433.20648175062206 + 63.905091254154904 * y)))))
gibbs_pt0_pt0 = (g03 + g08) * 0.000625
return np.ma.array(gibbs_pt0_pt0, mask=mask, copy=False)
def in_Baltic(lon, lat):
"""
Check if positions are in the Baltic Sea
Parameters
----------
lon, lat : array_like or masked arrays
Returns
-------
in_Baltic : boolean array (at least 1D)
True for points in the Baltic Sea
False for points outside, masked or NaN
"""
lon, lat = np.atleast_1d(lon, lat)
# Polygon bounding the Baltic, (xb, yb)
# Effective boundary is the intersection of this polygon
# with rectangle defined by xmin, xmax, ymin, ymax
# start with southwestern point and go round cyclonically
xb = np.array([12.6, 45.0, 26.0, 7.0, 12.6])
yb = np.array([50.0, 50.0, 69.0, 59.0, 50.0])
# Enclosing rectangle
# xmin, xmax = xb.min(), xb.max()
# ymin, ymax = yb.min(), yb.max()
xmin, xmax = 7.0, 32.0
ymin, ymax = 52.0, 67.0
# First check if outside the rectangle
in_rectangle = ((xmin < lon) & (lon < xmax) &
(ymin < lat) & (lat < ymax))
# Masked values are also considered outside the rectangle
if np.ma.is_masked(in_rectangle):
in_rectangle = in_rectangle.data & ~in_rectangle.mask
# Closer check for points in the rectangle
if np.any(in_rectangle):
lon, lat = np.broadcast_arrays(lon, lat, subok=True)
in_baltic = np.zeros(lon.shape, dtype='bool')
lon1 = lon[in_rectangle]
lat1 = lat[in_rectangle]
# There are general ways of testing for point in polygon
# This works for this special configuration of points
xx_right = np.interp(lat1, yb[1:3], xb[1:3])
xx_left = np.interp(lat1, yb[-1:1:-1], xb[-1:1:-1])
in_baltic[in_rectangle] = (xx_left <= lon1) & (lon1 <= xx_right)
return in_baltic
else: # Nothing inside the rectangle, return the False array.
return in_rectangle
[docs]def infunnel(SA, CT, p):
"""
Oceanographic funnel check for the 25-term equation
Parameters
----------
SA : array_like
Absolute Salinity [g/kg]
CT : array_like
Conservative Temperature [°C]
p : array_like
sea pressure [dbar]
(ie. absolute pressure - 10.1325 dbar)
Returns
-------
in_funnel : boolean ndarray or scalar
True, if SA, CT and p are inside the "funnel"
False, if SA, CT and p are outside the "funnel",
or one of the values are NaN or masked
Note. The term "funnel" describes the range of SA, CT and p over which
the error in the fit of the computationally-efficient 25-term
expression for density in terms of SA, CT and p was calculated
(McDougall et al., 2010).
"""
# Check variables and re-size if necessary.
scalar = np.isscalar(SA) and np.isscalar(CT) and np.isscalar(p)
SA, CT, p = np.broadcast_arrays(SA, CT, p, subok=True)
input_nan = np.isnan(SA) | np.isnan(CT) | np.isnan(p)
infunnel = ((p <= 8000) &
(SA >= 0) &
(SA <= 42.2) &
(CT >= (-0.3595467 - 0.0553734 * SA)) &
((p >= 5500) | (SA >= 0.006028 * (p - 500))) &
((p >= 5500) | (CT <= (33.0 - 0.003818181818182 * p))) &
((p <= 5500) | (SA >= 30.14)) &
((p <= 5500) | (CT <= 12.0)))
infunnel = infunnel & np.logical_not(input_nan)
if scalar:
infunnel = bool(infunnel)
return infunnel
[docs]def interp_SA_CT(SA, CT, p, p_i):
"""
function [SA_i, CT_i] = interp_SA_CT(SA,CT,p,p_i)
interp_SA_CT linear interpolation to p_i on a cast
==========================================================================
This function interpolates the cast with respect to the interpolating
variable p. This function finds the values of SA, CT at p_i on this cast.
"""
return interp_S_T(SA, CT, p, p_i)
def interp_S_T(S, T, z, znew, P=None):
"""
Linear interpolation of ndarrays *S* and *T* from *z* to *znew*.
Optionally interpolate a third ndarray, *P*.
*z* must be strictly increasing or strictly decreasing. It must
be a 1-D array, and its length must match the last dimension
of *S* and *T*.
*znew* may be a scalar or a sequence.
It is assumed, but not checked, that *S*, *T*, and *z* are
all plain ndarrays, not masked arrays or other sequences.
Out-of-range values of *znew*, and *nan* in *S* and *T*,
yield corresponding *nan* in the output.
The basic algorithm is from scipy.interpolate.
"""
isscalar = False
if not np.iterable(znew):
isscalar = True
znew = [znew]
znew = np.asarray(znew)
inverted = False
if z[1] - z[0] < 0:
inverted = True
z = z[::-1]
S = S[..., ::-1]
T = T[..., ::-1]
if P is not None:
P = P[..., ::-1]
if (np.diff(z) <= 0).any():
raise ValueError("z must be strictly increasing or decreasing")
hi = np.searchsorted(z, znew)
hi = hi.clip(1, len(z) - 1).astype(int)
lo = hi - 1
z_lo = z[lo]
z_hi = z[hi]
S_lo = S[lo]
S_hi = S[hi]
T_lo = T[lo]
T_hi = T[hi]
zratio = (znew - z_lo) / (z_hi - z_lo)
Si = S_lo + (S_hi - S_lo) * zratio
Ti = T_lo + (T_hi - T_lo) * zratio
if P is not None:
Pi = P[lo] + (P[hi] - P[lo]) * zratio
if inverted:
Si = Si[..., ::-1]
Ti = Ti[..., ::-1]
if P is not None:
Pi = Pi[..., ::-1]
outside = (znew < z.min()) | (znew > z.max())
if np.any(outside):
Si[..., outside] = np.nan
Ti[..., outside] = np.nan
if P is not None:
Pi[..., outside] = np.nan
if isscalar:
Si = Si[0]
Ti = Ti[0]
if P is not None:
Pi = Pi[0]
if P is None:
return Si, Ti
return Si, Ti, Pi
[docs]def interp_ref_cast(spycnl, A="gn"):
"""
Translation of:
[SA_iref_cast, CT_iref_cast, p_iref_cast] = interp_ref_cast(spycnl, A)
interp_ref_cast linear interpolation of the reference cast
==========================================================================
This function interpolates the reference cast with respect to the
interpolating variable "spycnl". This reference cast is at the location
188E,4N from the reference data set which underlies the Jackett &
McDougall (1997) Neutral Density computer code. This function finds the
values of SA, CT and p on this reference cast which correspond to the
value of isopycnal which is passed to this function from the function
"geo_strf_isopycnal_CT". The isopycnal could be either gamma_n or
sigma_2. If A is set to any of the following 's2','S2','sigma2','sigma_2'
the interpolation will take place in sigma 2 space, any other input
will result in the programme working in gamma_n space.
VERSION NUMBER: 3.0 (14th April, 2011)
REFERENCE:
Jackett, D. R. and T. J. McDougall, 1997: A neutral density variable
for the world<92>s oceans. Journal of Physical Oceanography, 27, 237-263.
FIXME? Do we need argument checking here to handle masked arrays,
etc.? I suspect not, since I don't think this is intended to be
user-callable, but is instead used internally by user-callable
functions.
Note: The v3.03 matlab code is incorrectly using approximate numbers
for the gamma_n case, even when the sigma_2 case is in effect.
That bug is fixed here.
"""
if A.lower() in ["s2", "sigma2", "sigma_2"]:
A = "s2"
gsw_data = read_data("gsw_data_v3_0.npz")
SA_ref = gsw_data.SA_ref_cast
CT_ref = gsw_data.CT_ref_cast
p_ref = gsw_data.p_ref_cast
if A == "s2":
zvar_ref = gsw_data.sigma_2_ref_cast
else:
zvar_ref = gsw_data.gamma_n_ref_cast
zvar_new = spycnl
Si, Ci, Pi = interp_S_T(SA_ref, CT_ref, zvar_ref, zvar_new, P=p_ref)
shallower = spycnl <= zvar_ref[0]
deeper = spycnl >= zvar_ref[-1]
if shallower.any():
Si[shallower] = SA_ref[0]
Ci[shallower] = CT_ref[0]
Pi[shallower] = p_ref[0]
if deeper.any():
Si[deeper] = SA_ref[-1]
Ci[deeper] = CT_ref[-1]
Pi[deeper] = p_ref[-1]
return Si, Ci, Pi
def specvol_SSO_0_CT25(p):
"""
Calculates specific volume at the Standard Ocean Salinity (SSO) and
Conservative Temperature of zero degrees C (CT=0), as a function of
pressure (p [dbar]) or spec_vol_CT25(35.16504,0,p).
Parameters
----------
p : array_like
pressure [dbar]
Returns
-------
specvol_SSO_0_CT25 : array_like
Specific volume at (SSO, CT=0, p), 25-term equation.
[m :sup:`3` kg :sup:`-1`]
Notes
-----
It uses a streamlined version of the 25-term CT version of specific volume
that is, a streamlined version of the code "rho_alpha_beta_CT25(SA,CT,p)"
Modifications
"""
p = np.asanyarray(p)
# No need to strip mask and replace it here; the calculation is simple.
specvol_SSO_0_CT25 = ((1. + SSO * (2.0777716085618458e-003 + np.sqrt(SSO) *
3.4688210757917340e-006) + p *
6.8314629554123324e-006) /
(9.9984380290708214e+002 + SSO *
(2.8925731541277653e+000 + SSO *
1.9457531751183059e-003) + p *
(1.1930681818531748e-002 + SSO *
5.9355685925035653e-006 + p *
-2.5943389807429039e-008)))
return specvol_SSO_0_CT25
# Salinity lib functions.
[docs]def specvol_SSO_0_p(p):
"""
This function calculates specific volume at the Standard Ocean
Salinity, SSO, and at a Conservative Temperature of zero degrees C, as a
function of pressure, p, in dbar, using a streamlined version of the
48-term CT version of specific volume, that is, a streamlined version of
the code "specvol(SA, CT, p)".
"""
v01 = 9.998420897506056e+2
v05 = -6.698001071123802
v08 = -3.988822378968490e-2
v12 = -2.233269627352527e-2
v15 = -1.806789763745328e-4
v17 = -3.087032500374211e-7
v20 = 1.550932729220080e-10
v21 = 1.0
v26 = -7.521448093615448e-3
v31 = -3.303308871386421e-5
v36 = 5.419326551148740e-6
v37 = -2.742185394906099e-5
v41 = -1.105097577149576e-7
v43 = -1.119011592875110e-10
v47 = -1.200507748551599e-15
return ((v21 + SSO * (v26 + v36 * SSO + v31 * np.sqrt(SSO)) + p *
(v37 + v41 * SSO + p * (v43 + v47 * p))) / (v01 + SSO *
(v05 + v08 * np.sqrt(SSO)) + p * (v12 + v15 * SSO + p *
(v17 + v20 * SSO))))
if __name__ == '__main__':
import doctest
doctest.testmod()