| 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: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/AUTO/input/densjmd95.m,v 1.2 2015/10/13 15:54:17 dgoldberg Exp $ |
| 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 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph