1 |
function rho = densjmd95(s,t,p); |
2 |
|
3 |
% |
4 |
% DENSJMD95 Density of sea water |
5 |
%========================================================================= |
6 |
% |
7 |
% USAGE: dens = densjmd95(S,Theta,P) |
8 |
% |
9 |
% DESCRIPTION: |
10 |
% Density of Sea Water using Jackett and McDougall 1995 (JAOT 12) |
11 |
% polynomial (modified UNESCO 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 |
% Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388 |
32 |
|
33 |
% created by mlosch on 2002-08-09 |
34 |
% $Header: $ |
35 |
% $Name: $ |
36 |
|
37 |
%---------------------- |
38 |
% CHECK INPUT ARGUMENTS |
39 |
%---------------------- |
40 |
if nargin ~=3 |
41 |
error('densjmd95.m: Must pass 3 parameters') |
42 |
end |
43 |
if ndims(s) > 2 |
44 |
dims = size(s); |
45 |
dimt = size(t); |
46 |
dimp = size(p); |
47 |
if length(dims) ~= length(dimt) | length(dims) ~= length(dimp) ... |
48 |
length(dimt) ~= length(dimp) |
49 |
error(['for more than two dimensions, S, Theta, and P must have the' ... |
50 |
' same number of dimensions']) |
51 |
else |
52 |
for k=length(dims) |
53 |
if dims(k)~=dimt(k) | dims(k)~=dimp(k) | dimt(k)~=dimp(k) |
54 |
error(['for more than two dimensions, S, Theta, and P must have' ... |
55 |
' the same dimensions']) |
56 |
end |
57 |
end |
58 |
end |
59 |
else |
60 |
% CHECK S,T,P dimensions and verify consistent |
61 |
[ms,ns] = size(s); |
62 |
[mt,nt] = size(t); |
63 |
[mp,np] = size(p); |
64 |
|
65 |
% CHECK THAT S & T HAVE SAME SHAPE |
66 |
if (ms~=mt) | (ns~=nt) |
67 |
error('check_stp: S & T must have same dimensions') |
68 |
end %if |
69 |
|
70 |
% CHECK OPTIONAL SHAPES FOR P |
71 |
if mp==1 & np==1 % P is a scalar. Fill to size of S |
72 |
p = p(1)*ones(ms,ns); |
73 |
elseif np==ns & mp==1 % P is row vector with same cols as S |
74 |
p = p( ones(1,ms), : ); % Copy down each column. |
75 |
elseif mp==ms & np==1 % P is column vector |
76 |
p = p( :, ones(1,ns) ); % Copy across each row |
77 |
elseif mp==ms & np==ns % P is a matrix size(S) |
78 |
% shape ok |
79 |
else |
80 |
error('check_stp: P has wrong dimensions') |
81 |
end %if |
82 |
[mp,np] = size(p); |
83 |
% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN. |
84 |
Transpose = 0; |
85 |
if mp == 1 % row vector |
86 |
p = p(:); |
87 |
t = t(:); |
88 |
s = s(:); |
89 |
Transpose = 1; |
90 |
end |
91 |
%***check_stp |
92 |
end |
93 |
|
94 |
% convert pressure to bar |
95 |
p = .1*p; |
96 |
|
97 |
% coefficients nonlinear equation of state in pressure coordinates for |
98 |
% 1. density of fresh water at p = 0 |
99 |
eosJMDCFw(1) = 999.842594; |
100 |
eosJMDCFw(2) = 6.793952e-02; |
101 |
eosJMDCFw(3) = - 9.095290e-03; |
102 |
eosJMDCFw(4) = 1.001685e-04; |
103 |
eosJMDCFw(5) = - 1.120083e-06; |
104 |
eosJMDCFw(6) = 6.536332e-09; |
105 |
% 2. density of sea water at p = 0 |
106 |
eosJMDCSw(1) = 8.244930e-01; |
107 |
eosJMDCSw(2) = - 4.089900e-03; |
108 |
eosJMDCSw(3) = 7.643800e-05 ; |
109 |
eosJMDCSw(4) = - 8.246700e-07; |
110 |
eosJMDCSw(5) = 5.387500e-09; |
111 |
eosJMDCSw(6) = - 5.724660e-03; |
112 |
eosJMDCSw(7) = 1.022700e-04; |
113 |
eosJMDCSw(8) = - 1.654600e-06; |
114 |
eosJMDCSw(9) = 4.831400e-04; |
115 |
|
116 |
t2 = t.*t; |
117 |
t3 = t2.*t; |
118 |
t4 = t3.*t; |
119 |
|
120 |
is = find(s(:) < 0 ); |
121 |
if ~isempty(is) |
122 |
warning('found negative salinity values, reset them to NaN'); |
123 |
s(is) = NaN; |
124 |
end |
125 |
s3o2 = s.*sqrt(s); |
126 |
|
127 |
% density of freshwater at the surface |
128 |
rho = eosJMDCFw(1) ... |
129 |
+ eosJMDCFw(2)*t ... |
130 |
+ eosJMDCFw(3)*t2 ... |
131 |
+ eosJMDCFw(4)*t3 ... |
132 |
+ eosJMDCFw(5)*t4 ... |
133 |
+ eosJMDCFw(6)*t4.*t; |
134 |
% density of sea water at the surface |
135 |
rho = rho ... |
136 |
+ s.*( ... |
137 |
eosJMDCSw(1) ... |
138 |
+ eosJMDCSw(2)*t ... |
139 |
+ eosJMDCSw(3)*t2 ... |
140 |
+ eosJMDCSw(4)*t3 ... |
141 |
+ eosJMDCSw(5)*t4 ... |
142 |
) ... |
143 |
+ s3o2.*( ... |
144 |
eosJMDCSw(6) ... |
145 |
+ eosJMDCSw(7)*t ... |
146 |
+ eosJMDCSw(8)*t2 ... |
147 |
) ... |
148 |
+ eosJMDCSw(9)*s.*s; |
149 |
|
150 |
rho = rho./(1 - p./bulkmodjmd95(s,t,p)); |
151 |
|
152 |
if ndims(s) < 3 & Transpose |
153 |
rho = rho'; |
154 |
end %if |
155 |
|
156 |
return |
157 |
|
158 |
function bulkmod = bulkmodjmd95(s,t,p) |
159 |
%function bulkmod = bulkmodjmd95(s,t,p) |
160 |
|
161 |
dummy = 0; |
162 |
% coefficients in pressure coordinates for |
163 |
% 3. secant bulk modulus K of fresh water at p = 0 |
164 |
eosJMDCKFw(1) = 1.965933e+04; |
165 |
eosJMDCKFw(2) = 1.444304e+02; |
166 |
eosJMDCKFw(3) = - 1.706103e+00; |
167 |
eosJMDCKFw(4) = 9.648704e-03; |
168 |
eosJMDCKFw(5) = - 4.190253e-05; |
169 |
% 4. secant bulk modulus K of sea water at p = 0 |
170 |
eosJMDCKSw(1) = 5.284855e+01; |
171 |
eosJMDCKSw(2) = - 3.101089e-01; |
172 |
eosJMDCKSw(3) = 6.283263e-03; |
173 |
eosJMDCKSw(4) = - 5.084188e-05; |
174 |
eosJMDCKSw(5) = 3.886640e-01; |
175 |
eosJMDCKSw(6) = 9.085835e-03; |
176 |
eosJMDCKSw(7) = - 4.619924e-04; |
177 |
% 5. secant bulk modulus K of sea water at p |
178 |
eosJMDCKP( 1) = 3.186519e+00; |
179 |
eosJMDCKP( 2) = 2.212276e-02; |
180 |
eosJMDCKP( 3) = - 2.984642e-04; |
181 |
eosJMDCKP( 4) = 1.956415e-06; |
182 |
eosJMDCKP( 5) = 6.704388e-03; |
183 |
eosJMDCKP( 6) = - 1.847318e-04; |
184 |
eosJMDCKP( 7) = 2.059331e-07; |
185 |
eosJMDCKP( 8) = 1.480266e-04; |
186 |
eosJMDCKP( 9) = 2.102898e-04; |
187 |
eosJMDCKP(10) = - 1.202016e-05; |
188 |
eosJMDCKP(11) = 1.394680e-07; |
189 |
eosJMDCKP(12) = - 2.040237e-06; |
190 |
eosJMDCKP(13) = 6.128773e-08; |
191 |
eosJMDCKP(14) = 6.207323e-10; |
192 |
|
193 |
t2 = t.*t; |
194 |
t3 = t2.*t; |
195 |
t4 = t3.*t; |
196 |
|
197 |
is = find(s(:) < 0 ); |
198 |
if ~isempty(is) |
199 |
warning('found negative salinity values, reset them to NaN'); |
200 |
s(is) = NaN; |
201 |
end |
202 |
s3o2 = s.*sqrt(s); |
203 |
%p = pressure(i,j,k,bi,bj)*SItoBar |
204 |
p2 = p.*p; |
205 |
% secant bulk modulus of fresh water at the surface |
206 |
bulkmod = eosJMDCKFw(1) ... |
207 |
+ eosJMDCKFw(2)*t ... |
208 |
+ eosJMDCKFw(3)*t2 ... |
209 |
+ eosJMDCKFw(4)*t3 ... |
210 |
+ eosJMDCKFw(5)*t4; |
211 |
% secant bulk modulus of sea water at the surface |
212 |
bulkmod = bulkmod ... |
213 |
+ s.*( eosJMDCKSw(1) ... |
214 |
+ eosJMDCKSw(2)*t ... |
215 |
+ eosJMDCKSw(3)*t2 ... |
216 |
+ eosJMDCKSw(4)*t3 ... |
217 |
) ... |
218 |
+ s3o2.*( eosJMDCKSw(5) ... |
219 |
+ eosJMDCKSw(6)*t ... |
220 |
+ eosJMDCKSw(7)*t2 ... |
221 |
); |
222 |
% secant bulk modulus of sea water at pressure p |
223 |
bulkmod = bulkmod ... |
224 |
+ p.*( eosJMDCKP(1) ... |
225 |
+ eosJMDCKP(2)*t ... |
226 |
+ eosJMDCKP(3)*t2 ... |
227 |
+ eosJMDCKP(4)*t3 ... |
228 |
) ... |
229 |
+ p.*s.*( eosJMDCKP(5) ... |
230 |
+ eosJMDCKP(6)*t ... |
231 |
+ eosJMDCKP(7)*t2 ... |
232 |
) ... |
233 |
+ p.*s3o2*eosJMDCKP(8) ... |
234 |
+ p2.*( eosJMDCKP(9) ... |
235 |
+ eosJMDCKP(10)*t ... |
236 |
+ eosJMDCKP(11)*t2 ... |
237 |
) ... |
238 |
+ p2.*s.*( eosJMDCKP(12) ... |
239 |
+ eosJMDCKP(13)*t ... |
240 |
+ eosJMDCKP(14)*t2 ... |
241 |
); |
242 |
|
243 |
return |
244 |
|