#
# converted from matlab version to python in Jan 2018
import sys
import numpy as np
__doc__ = """
Density of Sea Water using McDougall et al. 2003 (JAOT 20) polynomial
"""
# coefficients nonlinear equation of state in pressure coordinates for
eosMDJWFnum = [ 7.35212840e+00,
-5.45928211e-02,
3.98476704e-04,
2.96938239e+00,
-7.23268813e-03,
2.12382341e-03,
1.04004591e-02,
1.03970529e-07,
5.18761880e-06,
-3.24041825e-08,
-1.23869360e-11,
9.99843699e+02 ]
eosMDJWFden = [ 7.28606739e-03,
-4.60835542e-05,
3.68390573e-07,
1.80809186e-10,
2.14691708e-03,
-9.27062484e-06,
-1.78343643e-10,
4.76534122e-06,
1.63410736e-09,
5.30848875e-06,
-3.03175128e-16,
-1.27934137e-17,
1.00000000e+00 ]
[docs]def densmdjwf(s,theta,p):
"""
Computes in-situ density of sea water
Density of Sea Water using McDougall et al. 2003 (JAOT 20)
polynomial (Gibbs Potential).
Parameters
----------
s : array_like
salinity [psu (PSS-78)]
theta : array_like
potential temperature [degree C (IPTS-68)];
same shape as s
p : array_like
pressure [dbar]; broadcastable to shape of s
Returns
-------
dens : array
density [kg/m^3]
Example
-------
>>> densmdjwf(35., 25., 2000.)
1031.654229
Notes
-----
AUTHOR: Martin Losch 2002-08-09 (Martin.Losch@awi.de)
McDougall et al., 2003, JAOT 20(5), pp. 730-741
"""
# make sure arguments are floating point
s = np.asfarray(s)
t = np.asfarray(theta)
p = np.asfarray(p)
p1 = np.copy(p);
t1 = np.copy(t);
t2 = t*t;
s1 = np.copy(s);
if np.any(s1<0):
sys.stderr.write('negative salinity values! setting to nan\n')
# # the sqrt will take care of this
# if s.ndim > 0:
# s[s<0] = np.nan
# else:
# s = np.nan
sp5 = np.sqrt(s1)
p1t1=p1*t1
num = ( eosMDJWFnum[11]
+ t1*(eosMDJWFnum[0]
+ t1*(eosMDJWFnum[1] + eosMDJWFnum[2]*t1) )
+ s1*(eosMDJWFnum[3]
+ eosMDJWFnum[4]*t1 + eosMDJWFnum[5]*s1)
+ p1*(eosMDJWFnum[6] + eosMDJWFnum[7]*t2
+ eosMDJWFnum[8]*s1
+ p1*(eosMDJWFnum[9] + eosMDJWFnum[10]*t2) )
)
den = ( eosMDJWFden[12]
+ t1*(eosMDJWFden[0]
+ t1*(eosMDJWFden[1]
+ t1*(eosMDJWFden[2] + t1*eosMDJWFden[3] ) ) )
+ s1*(eosMDJWFden[4]
+ t1*(eosMDJWFden[5]
+ eosMDJWFden[6]*t2)
+ sp5*(eosMDJWFden[7] + eosMDJWFden[8]*t2) )
+ p1*(eosMDJWFden[9]
+ p1t1*(eosMDJWFden[10]*t2 + eosMDJWFden[11]*p1) )
)
epsln = 0
denom = 1.0/(epsln+den)
rho = num*denom;
return rho
# aliases
dens = densmdjwf