| 1 |
function rho = eh3densjmd95(s,t,p); |
| 2 |
|
| 3 |
% DENSJMD95 Density of sea water |
| 4 |
%========================================================================= |
| 5 |
|
| 6 |
% USAGE: dens = densjmd95(S,T,P) |
| 7 |
% |
| 8 |
% DESCRIPTION: |
| 9 |
% Density of Sea Water using Jackett and McDougall 1995 (JAOT 12) |
| 10 |
% polynomial (modified UNESCO polynomial). |
| 11 |
% Modified from |
| 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 |
|
| 35 |
%---------------------- |
| 36 |
% CHECK INPUT ARGUMENTS |
| 37 |
%---------------------- |
| 38 |
if nargin ~=3 |
| 39 |
error('densjmd95.m: Must pass 3 parameters') |
| 40 |
end |
| 41 |
if length(p) ~= 1 |
| 42 |
error('pressure must be a scalar -- otherwise use densjmd95()') |
| 43 |
end |
| 44 |
|
| 45 |
% convert pressure to bar |
| 46 |
p = 0.1 * p; |
| 47 |
|
| 48 |
% coefficients nonlinear equation of state in pressure coordinates for |
| 49 |
% 1. density of fresh water at p = 0 |
| 50 |
eosJMDCFw(1) = 999.842594; |
| 51 |
eosJMDCFw(2) = 6.793952e-02; |
| 52 |
eosJMDCFw(3) = - 9.095290e-03; |
| 53 |
eosJMDCFw(4) = 1.001685e-04; |
| 54 |
eosJMDCFw(5) = - 1.120083e-06; |
| 55 |
eosJMDCFw(6) = 6.536332e-09; |
| 56 |
% 2. density of sea water at p = 0 |
| 57 |
eosJMDCSw(1) = 8.244930e-01; |
| 58 |
eosJMDCSw(2) = - 4.089900e-03; |
| 59 |
eosJMDCSw(3) = 7.643800e-05; |
| 60 |
eosJMDCSw(4) = - 8.246700e-07; |
| 61 |
eosJMDCSw(5) = 5.387500e-09; |
| 62 |
eosJMDCSw(6) = - 5.724660e-03; |
| 63 |
eosJMDCSw(7) = 1.022700e-04; |
| 64 |
eosJMDCSw(8) = - 1.654600e-06; |
| 65 |
eosJMDCSw(9) = 4.831400e-04; |
| 66 |
|
| 67 |
t2 = t.*t; |
| 68 |
t3 = t2.*t; |
| 69 |
t4 = t3.*t; |
| 70 |
|
| 71 |
is = find(s(:) < 0 ); |
| 72 |
if ~isempty(is) |
| 73 |
warning('found negative salinity values, reset them to NaN'); |
| 74 |
s(is) = NaN; |
| 75 |
end |
| 76 |
s3o2 = s.*sqrt(s); |
| 77 |
|
| 78 |
% density of freshwater at the surface |
| 79 |
rho = eosJMDCFw(1) + eosJMDCFw(2)*t + eosJMDCFw(3)*t2 ... |
| 80 |
+ eosJMDCFw(4)*t3 + eosJMDCFw(5)*t4 + eosJMDCFw(6)*t4.*t; |
| 81 |
% density of sea water at the surface |
| 82 |
rho = rho ... |
| 83 |
+ s.*( eosJMDCSw(1) + eosJMDCSw(2)*t + eosJMDCSw(3)*t2 ... |
| 84 |
+ eosJMDCSw(4)*t3 + eosJMDCSw(5)*t4 ) ... |
| 85 |
+ s3o2.*( eosJMDCSw(6) + eosJMDCSw(7)*t + eosJMDCSw(8)*t2 ) ... |
| 86 |
+ eosJMDCSw(9)*s.*s; |
| 87 |
|
| 88 |
rho = rho./(1 - p/bulkmodjmd95(s,t,p)); |
| 89 |
|
| 90 |
return |
| 91 |
|
| 92 |
%======================================================= |
| 93 |
|
| 94 |
function bulkmod = bulkmodjmd95(s,t,p) |
| 95 |
%function bulkmod = bulkmodjmd95(s,t,p) |
| 96 |
|
| 97 |
dummy = 0; |
| 98 |
% coefficients in pressure coordinates for |
| 99 |
% 3. secant bulk modulus K of fresh water at p = 0 |
| 100 |
eosJMDCKFw(1) = 1.965933e+04; |
| 101 |
eosJMDCKFw(2) = 1.444304e+02; |
| 102 |
eosJMDCKFw(3) = - 1.706103e+00; |
| 103 |
eosJMDCKFw(4) = 9.648704e-03; |
| 104 |
eosJMDCKFw(5) = - 4.190253e-05; |
| 105 |
% 4. secant bulk modulus K of sea water at p = 0 |
| 106 |
eosJMDCKSw(1) = 5.284855e+01; |
| 107 |
eosJMDCKSw(2) = - 3.101089e-01; |
| 108 |
eosJMDCKSw(3) = 6.283263e-03; |
| 109 |
eosJMDCKSw(4) = - 5.084188e-05; |
| 110 |
eosJMDCKSw(5) = 3.886640e-01; |
| 111 |
eosJMDCKSw(6) = 9.085835e-03; |
| 112 |
eosJMDCKSw(7) = - 4.619924e-04; |
| 113 |
% 5. secant bulk modulus K of sea water at p |
| 114 |
eosJMDCKP( 1) = 3.186519e+00; |
| 115 |
eosJMDCKP( 2) = 2.212276e-02; |
| 116 |
eosJMDCKP( 3) = - 2.984642e-04; |
| 117 |
eosJMDCKP( 4) = 1.956415e-06; |
| 118 |
eosJMDCKP( 5) = 6.704388e-03; |
| 119 |
eosJMDCKP( 6) = - 1.847318e-04; |
| 120 |
eosJMDCKP( 7) = 2.059331e-07; |
| 121 |
eosJMDCKP( 8) = 1.480266e-04; |
| 122 |
eosJMDCKP( 9) = 2.102898e-04; |
| 123 |
eosJMDCKP(10) = - 1.202016e-05; |
| 124 |
eosJMDCKP(11) = 1.394680e-07; |
| 125 |
eosJMDCKP(12) = - 2.040237e-06; |
| 126 |
eosJMDCKP(13) = 6.128773e-08; |
| 127 |
eosJMDCKP(14) = 6.207323e-10; |
| 128 |
|
| 129 |
t2 = t.*t; |
| 130 |
t3 = t2.*t; |
| 131 |
t4 = t3.*t; |
| 132 |
|
| 133 |
is = find(s(:) < 0 ); |
| 134 |
if ~isempty(is) |
| 135 |
warning('found negative salinity values, reset them to NaN'); |
| 136 |
s(is) = NaN; |
| 137 |
end |
| 138 |
s3o2 = s.*sqrt(s); |
| 139 |
%p = pressure(i,j,k,bi,bj)*SItoBar |
| 140 |
p2 = p*p; |
| 141 |
% secant bulk modulus of fresh water at the surface |
| 142 |
bulkmod = eosJMDCKFw(1) + eosJMDCKFw(2)*t + eosJMDCKFw(3)*t2 ... |
| 143 |
+ eosJMDCKFw(4)*t3 + eosJMDCKFw(5)*t4; |
| 144 |
% secant bulk modulus of sea water at the surface |
| 145 |
bulkmod = bulkmod ... |
| 146 |
+ s.*( eosJMDCKSw(1) + eosJMDCKSw(2)*t ... |
| 147 |
+ eosJMDCKSw(3)*t2 + eosJMDCKSw(4)*t3 ) ... |
| 148 |
+ s3o2.*( eosJMDCKSw(5) + eosJMDCKSw(6)*t ... |
| 149 |
+ eosJMDCKSw(7)*t2 ); |
| 150 |
% secant bulk modulus of sea water at pressure p |
| 151 |
bulkmod = bulkmod ... |
| 152 |
+ p *( eosJMDCKP(1) + eosJMDCKP(2)*t ... |
| 153 |
+ eosJMDCKP(3)*t2 + eosJMDCKP(4)*t3 ) ... |
| 154 |
+ p *s.*( eosJMDCKP(5) + eosJMDCKP(6)*t + eosJMDCKP(7)*t2 ) ... |
| 155 |
+ p *s3o2*eosJMDCKP(8) ... |
| 156 |
+ p2 *( eosJMDCKP(9) + eosJMDCKP(10)*t + eosJMDCKP(11)*t2 ) ... |
| 157 |
+ p2 *s.*( eosJMDCKP(12) + eosJMDCKP(13)*t + eosJMDCKP(14)*t2 ); |
| 158 |
|
| 159 |
return |
| 160 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph