| 1 |
function rho = densmdjwf(s,t,p); |
| 2 |
%function rho = densmdjwf(S,Theta,P); |
| 3 |
|
| 4 |
% DENSMDJWF Density of sea water |
| 5 |
%========================================================================= |
| 6 |
|
| 7 |
% USAGE: dens = densmdjwf(S,Theta,P) |
| 8 |
% |
| 9 |
% DESCRIPTION: |
| 10 |
% Density of Sea Water using the McDougall et al. 2003 (JAOT 20) |
| 11 |
% polynomial. |
| 12 |
% |
| 13 |
% INPUT: (all must have same dimensions) |
| 14 |
% S = salinity [psu (PSS-78)] |
| 15 |
% Theta = potential temperature [degree C (IPTS-68)] |
| 16 |
% P = pressure [dbar] |
| 17 |
% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) ) |
| 18 |
% |
| 19 |
% OUTPUT: |
| 20 |
% dens = density [kg/m^3] |
| 21 |
% |
| 22 |
% AUTHOR: Martin Losch 2002-08-09 (mlosch@mit.edu) |
| 23 |
% |
| 24 |
% check value |
| 25 |
% S = 35.5 PSU |
| 26 |
% Theta = 3 degC |
| 27 |
% P = 3000 dbar |
| 28 |
% rho = 1041.83267 kg/m^3 |
| 29 |
% |
| 30 |
|
| 31 |
% McDougall et al., 2003, JAOT 20(5), pp. 730-741 |
| 32 |
|
| 33 |
% created by mlosch on 2002-08-09 |
| 34 |
|
| 35 |
%---------------------- |
| 36 |
% CHECK INPUT ARGUMENTS |
| 37 |
%---------------------- |
| 38 |
if nargin ~=3 |
| 39 |
error('densjmd95.m: Must pass 3 parameters') |
| 40 |
end |
| 41 |
if ndims(s) > 2 |
| 42 |
dims = size(s); |
| 43 |
dimt = size(t); |
| 44 |
dimp = size(p); |
| 45 |
if length(dims) ~= length(dimt) | length(dims) ~= length(dimp) ... |
| 46 |
length(dimt) ~= length(dimp) |
| 47 |
error(['for more than two dimensions, S, Theta, and P must have the' ... |
| 48 |
' same number of dimensions']) |
| 49 |
else |
| 50 |
for k=length(dims) |
| 51 |
if dims(k)~=dimt(k) | dims(k)~=dimp(k) | dimt(k)~=dimp(k) |
| 52 |
error(['for more than two dimensions, S, Theta, and P must have' ... |
| 53 |
' the same dimensions']) |
| 54 |
end |
| 55 |
end |
| 56 |
end |
| 57 |
else |
| 58 |
% CHECK S,T,P dimensions and verify consistent |
| 59 |
[ms,ns] = size(s); |
| 60 |
[mt,nt] = size(t); |
| 61 |
[mp,np] = size(p); |
| 62 |
|
| 63 |
% CHECK THAT S & T HAVE SAME SHAPE |
| 64 |
if (ms~=mt) | (ns~=nt) |
| 65 |
error('check_stp: S & T must have same dimensions') |
| 66 |
end %if |
| 67 |
|
| 68 |
% CHECK OPTIONAL SHAPES FOR P |
| 69 |
if mp==1 & np==1 % P is a scalar. Fill to size of S |
| 70 |
p = p(1)*ones(ms,ns); |
| 71 |
elseif np==ns & mp==1 % P is row vector with same cols as S |
| 72 |
p = p( ones(1,ms), : ); % Copy down each column. |
| 73 |
elseif mp==ms & np==1 % P is column vector |
| 74 |
p = p( :, ones(1,ns) ); % Copy across each row |
| 75 |
elseif mp==ms & np==ns % P is a matrix size(S) |
| 76 |
% shape ok |
| 77 |
else |
| 78 |
error('check_stp: P has wrong dimensions') |
| 79 |
end %if |
| 80 |
[mp,np] = size(p); |
| 81 |
% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN. |
| 82 |
Transpose = 0; |
| 83 |
if mp == 1 % row vector |
| 84 |
p = p(:); |
| 85 |
t = t(:); |
| 86 |
s = s(:); |
| 87 |
Transpose = 1; |
| 88 |
end |
| 89 |
%***check_stp |
| 90 |
end |
| 91 |
|
| 92 |
% coefficients nonlinear equation of state in pressure coordinates for |
| 93 |
eosMDJWFnum(12) = 9.99843699e+02; |
| 94 |
eosMDJWFnum( 1) = 7.35212840e+00; |
| 95 |
eosMDJWFnum( 2) = -5.45928211e-02; |
| 96 |
eosMDJWFnum( 3) = 3.98476704e-04; |
| 97 |
eosMDJWFnum( 4) = 2.96938239e+00; |
| 98 |
eosMDJWFnum( 5) = -7.23268813e-03; |
| 99 |
eosMDJWFnum( 6) = 2.12382341e-03; |
| 100 |
eosMDJWFnum( 7) = 1.04004591e-02; |
| 101 |
eosMDJWFnum( 8) = 1.03970529e-07; |
| 102 |
eosMDJWFnum( 9) = 5.18761880e-06; |
| 103 |
eosMDJWFnum(10) = -3.24041825e-08; |
| 104 |
eosMDJWFnum(11) = -1.23869360e-11; |
| 105 |
|
| 106 |
eosMDJWFden(13) = 1.00000000e+00; |
| 107 |
eosMDJWFden( 1) = 7.28606739e-03; |
| 108 |
eosMDJWFden( 2) = -4.60835542e-05; |
| 109 |
eosMDJWFden( 3) = 3.68390573e-07; |
| 110 |
eosMDJWFden( 4) = 1.80809186e-10; |
| 111 |
eosMDJWFden( 5) = 2.14691708e-03; |
| 112 |
eosMDJWFden( 6) = -9.27062484e-06; |
| 113 |
eosMDJWFden( 7) = -1.78343643e-10; |
| 114 |
eosMDJWFden( 8) = 4.76534122e-06; |
| 115 |
eosMDJWFden( 9) = 1.63410736e-09; |
| 116 |
eosMDJWFden(10) = 5.30848875e-06; |
| 117 |
eosMDJWFden(11) = -3.03175128e-16; |
| 118 |
eosMDJWFden(12) = -1.27934137e-17; |
| 119 |
|
| 120 |
p1 = p; |
| 121 |
|
| 122 |
t1 = t; |
| 123 |
t2 = t.*t; |
| 124 |
|
| 125 |
s1=s; |
| 126 |
is = find(s1(:) < 0 ); |
| 127 |
if ~isempty(is) |
| 128 |
warning('found negative salinity values, reset them to NaN'); |
| 129 |
s1(is) = NaN; |
| 130 |
end |
| 131 |
sp5 = sqrt(s1); |
| 132 |
p1t1=p1.*t1; |
| 133 |
% |
| 134 |
num = eosMDJWFnum(12) ... |
| 135 |
+ t1.*(eosMDJWFnum(1) ... |
| 136 |
+ t1.*(eosMDJWFnum(2) + eosMDJWFnum(3).*t1) ) ... |
| 137 |
+ s1.*(eosMDJWFnum(4) ... |
| 138 |
+ eosMDJWFnum(5).*t1 + eosMDJWFnum(6).*s1) ... |
| 139 |
+ p1.*(eosMDJWFnum(7) + eosMDJWFnum(8).*t2 ... |
| 140 |
+ eosMDJWFnum(9).*s1 ... |
| 141 |
+ p1.*(eosMDJWFnum(10) + eosMDJWFnum(11).*t2) ); |
| 142 |
den = eosMDJWFden(13) ... |
| 143 |
+ t1.*(eosMDJWFden(1) ... |
| 144 |
+ t1.*(eosMDJWFden(2) ... |
| 145 |
+ t1.*(eosMDJWFden(3) + t1.*eosMDJWFden(4) ) ) ) ... |
| 146 |
+ s1.*(eosMDJWFden(5) ... |
| 147 |
+ t1.*(eosMDJWFden(6) ... |
| 148 |
+ eosMDJWFden(7).*t2) ... |
| 149 |
+ sp5.*(eosMDJWFden(8) + eosMDJWFden(9).*t2) ) ... |
| 150 |
+ p1.*(eosMDJWFden(10) ... |
| 151 |
+ p1t1.*(eosMDJWFden(11).*t2 + eosMDJWFden(12).*p1) ); |
| 152 |
|
| 153 |
epsln = 0; |
| 154 |
denom = 1.0./(epsln+den) ; |
| 155 |
|
| 156 |
|
| 157 |
rho = num.*denom; |
| 158 |
|
| 159 |
return |
| 160 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph