Warning, /utils/matlab/densjmd95.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
4a6ce70b57 Ed H*0001 function rho = densjmd95(s,t,p);
0002
7f0effbc3a Jean*0003 %
4a6ce70b57 Ed H*0004 % DENSJMD95 Density of sea water
0005 %=========================================================================
7f0effbc3a Jean*0006 %
4a6ce70b57 Ed H*0007 % USAGE: dens = densjmd95(S,Theta,P)
0008 %
0009 % DESCRIPTION:
0010 % Density of Sea Water using Jackett and McDougall 1995 (JAOT 12)
0011 % polynomial (modified UNESCO polynomial).
0012 %
0013 % INPUT: (all must have same dimensions)
0014 % S = salinity [psu (PSS-78)]
0015 % Theta = potential temperature [degree C (IPTS-68)]
0016 % P = pressure [dbar]
0017 % (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
0018 %
0019 % OUTPUT:
0020 % dens = density [kg/m^3]
0021 %
0022 % AUTHOR: Martin Losch 2002-08-09 (mlosch@mit.edu)
0023 %
0024 % check value
0025 % S = 35.5 PSU
0026 % Theta = 3 degC
0027 % P = 3000 dbar
0028 % rho = 1041.83267 kg/m^3
7f0effbc3a Jean*0029
4a6ce70b57 Ed H*0030
0031 % Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388
0032
0033 % created by mlosch on 2002-08-09
0034 %----------------------
0035 % CHECK INPUT ARGUMENTS
0036 %----------------------
0037 if nargin ~=3
0038 error('densjmd95.m: Must pass 3 parameters')
0039 end
0040 if ndims(s) > 2
0041 dims = size(s);
0042 dimt = size(t);
0043 dimp = size(p);
0044 if length(dims) ~= length(dimt) | length(dims) ~= length(dimp) ...
0045 length(dimt) ~= length(dimp)
0046 error(['for more than two dimensions, S, Theta, and P must have the' ...
0047 ' same number of dimensions'])
0048 else
0049 for k=length(dims)
0050 if dims(k)~=dimt(k) | dims(k)~=dimp(k) | dimt(k)~=dimp(k)
0051 error(['for more than two dimensions, S, Theta, and P must have' ...
0052 ' the same dimensions'])
0053 end
0054 end
0055 end
0056 else
0057 % CHECK S,T,P dimensions and verify consistent
0058 [ms,ns] = size(s);
0059 [mt,nt] = size(t);
0060 [mp,np] = size(p);
0061
0062 % CHECK THAT S & T HAVE SAME SHAPE
0063 if (ms~=mt) | (ns~=nt)
0064 error('check_stp: S & T must have same dimensions')
0065 end %if
0066
0067 % CHECK OPTIONAL SHAPES FOR P
0068 if mp==1 & np==1 % P is a scalar. Fill to size of S
0069 p = p(1)*ones(ms,ns);
0070 elseif np==ns & mp==1 % P is row vector with same cols as S
0071 p = p( ones(1,ms), : ); % Copy down each column.
0072 elseif mp==ms & np==1 % P is column vector
0073 p = p( :, ones(1,ns) ); % Copy across each row
0074 elseif mp==ms & np==ns % P is a matrix size(S)
0075 % shape ok
0076 else
0077 error('check_stp: P has wrong dimensions')
0078 end %if
0079 [mp,np] = size(p);
0080 % IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
0081 Transpose = 0;
0082 if mp == 1 % row vector
0083 p = p(:);
0084 t = t(:);
0085 s = s(:);
0086 Transpose = 1;
0087 end
0088 %***check_stp
0089 end
0090
0091 % convert pressure to bar
0092 p = .1*p;
0093
0094 % coefficients nonlinear equation of state in pressure coordinates for
0095 % 1. density of fresh water at p = 0
0096 eosJMDCFw(1) = 999.842594;
0097 eosJMDCFw(2) = 6.793952e-02;
0098 eosJMDCFw(3) = - 9.095290e-03;
0099 eosJMDCFw(4) = 1.001685e-04;
0100 eosJMDCFw(5) = - 1.120083e-06;
0101 eosJMDCFw(6) = 6.536332e-09;
0102 % 2. density of sea water at p = 0
0103 eosJMDCSw(1) = 8.244930e-01;
0104 eosJMDCSw(2) = - 4.089900e-03;
0105 eosJMDCSw(3) = 7.643800e-05 ;
0106 eosJMDCSw(4) = - 8.246700e-07;
0107 eosJMDCSw(5) = 5.387500e-09;
0108 eosJMDCSw(6) = - 5.724660e-03;
0109 eosJMDCSw(7) = 1.022700e-04;
0110 eosJMDCSw(8) = - 1.654600e-06;
0111 eosJMDCSw(9) = 4.831400e-04;
0112
0113 t2 = t.*t;
0114 t3 = t2.*t;
0115 t4 = t3.*t;
0116
0117 is = find(s(:) < 0 );
0118 if ~isempty(is)
0119 warning('found negative salinity values, reset them to NaN');
0120 s(is) = NaN;
0121 end
0122 s3o2 = s.*sqrt(s);
0123
0124 % density of freshwater at the surface
0125 rho = eosJMDCFw(1) ...
0126 + eosJMDCFw(2)*t ...
0127 + eosJMDCFw(3)*t2 ...
0128 + eosJMDCFw(4)*t3 ...
0129 + eosJMDCFw(5)*t4 ...
0130 + eosJMDCFw(6)*t4.*t;
0131 % density of sea water at the surface
0132 rho = rho ...
0133 + s.*( ...
0134 eosJMDCSw(1) ...
0135 + eosJMDCSw(2)*t ...
0136 + eosJMDCSw(3)*t2 ...
0137 + eosJMDCSw(4)*t3 ...
0138 + eosJMDCSw(5)*t4 ...
0139 ) ...
0140 + s3o2.*( ...
0141 eosJMDCSw(6) ...
0142 + eosJMDCSw(7)*t ...
0143 + eosJMDCSw(8)*t2 ...
0144 ) ...
0145 + eosJMDCSw(9)*s.*s;
0146
0147 rho = rho./(1 - p./bulkmodjmd95(s,t,p));
0148
0149 if ndims(s) < 3 & Transpose
0150 rho = rho';
0151 end %if
0152
0153 return
0154
0155 function bulkmod = bulkmodjmd95(s,t,p)
0156 %function bulkmod = bulkmodjmd95(s,t,p)
0157
0158 dummy = 0;
0159 % coefficients in pressure coordinates for
0160 % 3. secant bulk modulus K of fresh water at p = 0
0161 eosJMDCKFw(1) = 1.965933e+04;
0162 eosJMDCKFw(2) = 1.444304e+02;
0163 eosJMDCKFw(3) = - 1.706103e+00;
0164 eosJMDCKFw(4) = 9.648704e-03;
0165 eosJMDCKFw(5) = - 4.190253e-05;
0166 % 4. secant bulk modulus K of sea water at p = 0
0167 eosJMDCKSw(1) = 5.284855e+01;
0168 eosJMDCKSw(2) = - 3.101089e-01;
0169 eosJMDCKSw(3) = 6.283263e-03;
0170 eosJMDCKSw(4) = - 5.084188e-05;
0171 eosJMDCKSw(5) = 3.886640e-01;
0172 eosJMDCKSw(6) = 9.085835e-03;
0173 eosJMDCKSw(7) = - 4.619924e-04;
0174 % 5. secant bulk modulus K of sea water at p
0175 eosJMDCKP( 1) = 3.186519e+00;
0176 eosJMDCKP( 2) = 2.212276e-02;
0177 eosJMDCKP( 3) = - 2.984642e-04;
0178 eosJMDCKP( 4) = 1.956415e-06;
0179 eosJMDCKP( 5) = 6.704388e-03;
0180 eosJMDCKP( 6) = - 1.847318e-04;
0181 eosJMDCKP( 7) = 2.059331e-07;
0182 eosJMDCKP( 8) = 1.480266e-04;
0183 eosJMDCKP( 9) = 2.102898e-04;
0184 eosJMDCKP(10) = - 1.202016e-05;
0185 eosJMDCKP(11) = 1.394680e-07;
0186 eosJMDCKP(12) = - 2.040237e-06;
0187 eosJMDCKP(13) = 6.128773e-08;
0188 eosJMDCKP(14) = 6.207323e-10;
0189
0190 t2 = t.*t;
0191 t3 = t2.*t;
0192 t4 = t3.*t;
0193
0194 is = find(s(:) < 0 );
0195 if ~isempty(is)
0196 warning('found negative salinity values, reset them to NaN');
0197 s(is) = NaN;
0198 end
0199 s3o2 = s.*sqrt(s);
0200 %p = pressure(i,j,k,bi,bj)*SItoBar
0201 p2 = p.*p;
0202 % secant bulk modulus of fresh water at the surface
0203 bulkmod = eosJMDCKFw(1) ...
0204 + eosJMDCKFw(2)*t ...
0205 + eosJMDCKFw(3)*t2 ...
0206 + eosJMDCKFw(4)*t3 ...
0207 + eosJMDCKFw(5)*t4;
0208 % secant bulk modulus of sea water at the surface
0209 bulkmod = bulkmod ...
0210 + s.*( eosJMDCKSw(1) ...
0211 + eosJMDCKSw(2)*t ...
0212 + eosJMDCKSw(3)*t2 ...
0213 + eosJMDCKSw(4)*t3 ...
0214 ) ...
0215 + s3o2.*( eosJMDCKSw(5) ...
0216 + eosJMDCKSw(6)*t ...
0217 + eosJMDCKSw(7)*t2 ...
0218 );
0219 % secant bulk modulus of sea water at pressure p
0220 bulkmod = bulkmod ...
0221 + p.*( eosJMDCKP(1) ...
0222 + eosJMDCKP(2)*t ...
0223 + eosJMDCKP(3)*t2 ...
0224 + eosJMDCKP(4)*t3 ...
0225 ) ...
0226 + p.*s.*( eosJMDCKP(5) ...
0227 + eosJMDCKP(6)*t ...
0228 + eosJMDCKP(7)*t2 ...
0229 ) ...
0230 + p.*s3o2*eosJMDCKP(8) ...
0231 + p2.*( eosJMDCKP(9) ...
0232 + eosJMDCKP(10)*t ...
0233 + eosJMDCKP(11)*t2 ...
0234 ) ...
0235 + p2.*s.*( eosJMDCKP(12) ...
0236 + eosJMDCKP(13)*t ...
0237 + eosJMDCKP(14)*t2 ...
0238 );
0239
0240 return
0241
