Warning, /utils/matlab/densmdjwf.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 = densmdjwf(s,t,p);
0002
7f0effbc3a Jean*0003 %
4a6ce70b57 Ed H*0004 % DENSMDJWF Density of sea water
0005 %=========================================================================
7f0effbc3a Jean*0006 %
4a6ce70b57 Ed H*0007 % USAGE: dens = densmdjwf(S,Theta,P)
0008 %
0009 % DESCRIPTION:
a88d221e3f Mart*0010 % Density of Sea Water using the McDougall et al. 2003 (JAOT 20)
0011 % polynomial.
4a6ce70b57 Ed H*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
a88d221e3f Mart*0031 % McDougall et al., 2003, JAOT 20(5), pp. 730-741
4a6ce70b57 Ed H*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 % coefficients nonlinear equation of state in pressure coordinates for
0092 eosMDJWFnum(12) = 9.99843699e+02;
0093 eosMDJWFnum( 1) = 7.35212840e+00;
0094 eosMDJWFnum( 2) = -5.45928211e-02;
0095 eosMDJWFnum( 3) = 3.98476704e-04;
0096 eosMDJWFnum( 4) = 2.96938239e+00;
0097 eosMDJWFnum( 5) = -7.23268813e-03;
0098 eosMDJWFnum( 6) = 2.12382341e-03;
0099 eosMDJWFnum( 7) = 1.04004591e-02;
0100 eosMDJWFnum( 8) = 1.03970529e-07;
0101 eosMDJWFnum( 9) = 5.18761880e-06;
0102 eosMDJWFnum(10) = -3.24041825e-08;
0103 eosMDJWFnum(11) = -1.23869360e-11;
0104
0105 eosMDJWFden(13) = 1.00000000e+00;
0106 eosMDJWFden( 1) = 7.28606739e-03;
0107 eosMDJWFden( 2) = -4.60835542e-05;
0108 eosMDJWFden( 3) = 3.68390573e-07;
0109 eosMDJWFden( 4) = 1.80809186e-10;
0110 eosMDJWFden( 5) = 2.14691708e-03;
0111 eosMDJWFden( 6) = -9.27062484e-06;
0112 eosMDJWFden( 7) = -1.78343643e-10;
0113 eosMDJWFden( 8) = 4.76534122e-06;
0114 eosMDJWFden( 9) = 1.63410736e-09;
0115 eosMDJWFden(10) = 5.30848875e-06;
0116 eosMDJWFden(11) = -3.03175128e-16;
0117 eosMDJWFden(12) = -1.27934137e-17;
0118
0119 p1 = p;
0120
0121 t1 = t;
0122 t2 = t.*t;
0123
0124 s1=s;
0125 is = find(s1(:) < 0 );
0126 if ~isempty(is)
0127 warning('found negative salinity values, reset them to NaN');
0128 s1(is) = NaN;
0129 end
0130 sp5 = sqrt(s1);
0131 p1t1=p1.*t1;
0132 %
0133 num = eosMDJWFnum(12) ...
0134 + t1.*(eosMDJWFnum(1) ...
0135 + t1.*(eosMDJWFnum(2) + eosMDJWFnum(3).*t1) ) ...
0136 + s1.*(eosMDJWFnum(4) ...
0137 + eosMDJWFnum(5).*t1 + eosMDJWFnum(6).*s1) ...
0138 + p1.*(eosMDJWFnum(7) + eosMDJWFnum(8).*t2 ...
0139 + eosMDJWFnum(9).*s1 ...
0140 + p1.*(eosMDJWFnum(10) + eosMDJWFnum(11).*t2) );
0141 den = eosMDJWFden(13) ...
0142 + t1.*(eosMDJWFden(1) ...
0143 + t1.*(eosMDJWFden(2) ...
0144 + t1.*(eosMDJWFden(3) + t1.*eosMDJWFden(4) ) ) ) ...
0145 + s1.*(eosMDJWFden(5) ...
0146 + t1.*(eosMDJWFden(6) ...
0147 + eosMDJWFden(7).*t2) ...
0148 + sp5.*(eosMDJWFden(8) + eosMDJWFden(9).*t2) ) ...
0149 + p1.*(eosMDJWFden(10) ...
0150 + p1t1.*(eosMDJWFden(11).*t2 + eosMDJWFden(12).*p1) );
0151
0152 epsln = 0;
0153 denom = 1.0./(epsln+den) ;
0154
0155
0156 rho = num.*denom;
0157
0158 return
0159
