/[MITgcm]/MITgcm/utils/matlab/densmdjwf.m
ViewVC logotype

Contents of /MITgcm/utils/matlab/densmdjwf.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (show annotations) (download)
Sat Feb 17 23:49:43 2007 UTC (14 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint58w_post, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint60, checkpoint61, checkpoint62, checkpoint63, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint58x_post, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint59, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint58y_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y, HEAD
Changes since 1.2: +5 -3 lines
add $Header:  $ and $Name:  & (for CVS)

1 function rho = densmdjwf(s,t,p);
2
3 %
4 % DENSMDJWF Density of sea water
5 %=========================================================================
6 %
7 % USAGE: dens = densmdjwf(S,Theta,P)
8 %
9 % DESCRIPTION:
10 % Density of Sea Water using the McDougall et al. 2003 (JAOT 20)
11 % 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 % McDougall et al., 2003, JAOT 20(5), pp. 730-741
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 % coefficients nonlinear equation of state in pressure coordinates for
95 eosMDJWFnum(12) = 9.99843699e+02;
96 eosMDJWFnum( 1) = 7.35212840e+00;
97 eosMDJWFnum( 2) = -5.45928211e-02;
98 eosMDJWFnum( 3) = 3.98476704e-04;
99 eosMDJWFnum( 4) = 2.96938239e+00;
100 eosMDJWFnum( 5) = -7.23268813e-03;
101 eosMDJWFnum( 6) = 2.12382341e-03;
102 eosMDJWFnum( 7) = 1.04004591e-02;
103 eosMDJWFnum( 8) = 1.03970529e-07;
104 eosMDJWFnum( 9) = 5.18761880e-06;
105 eosMDJWFnum(10) = -3.24041825e-08;
106 eosMDJWFnum(11) = -1.23869360e-11;
107
108 eosMDJWFden(13) = 1.00000000e+00;
109 eosMDJWFden( 1) = 7.28606739e-03;
110 eosMDJWFden( 2) = -4.60835542e-05;
111 eosMDJWFden( 3) = 3.68390573e-07;
112 eosMDJWFden( 4) = 1.80809186e-10;
113 eosMDJWFden( 5) = 2.14691708e-03;
114 eosMDJWFden( 6) = -9.27062484e-06;
115 eosMDJWFden( 7) = -1.78343643e-10;
116 eosMDJWFden( 8) = 4.76534122e-06;
117 eosMDJWFden( 9) = 1.63410736e-09;
118 eosMDJWFden(10) = 5.30848875e-06;
119 eosMDJWFden(11) = -3.03175128e-16;
120 eosMDJWFden(12) = -1.27934137e-17;
121
122 p1 = p;
123
124 t1 = t;
125 t2 = t.*t;
126
127 s1=s;
128 is = find(s1(:) < 0 );
129 if ~isempty(is)
130 warning('found negative salinity values, reset them to NaN');
131 s1(is) = NaN;
132 end
133 sp5 = sqrt(s1);
134 p1t1=p1.*t1;
135 %
136 num = eosMDJWFnum(12) ...
137 + t1.*(eosMDJWFnum(1) ...
138 + t1.*(eosMDJWFnum(2) + eosMDJWFnum(3).*t1) ) ...
139 + s1.*(eosMDJWFnum(4) ...
140 + eosMDJWFnum(5).*t1 + eosMDJWFnum(6).*s1) ...
141 + p1.*(eosMDJWFnum(7) + eosMDJWFnum(8).*t2 ...
142 + eosMDJWFnum(9).*s1 ...
143 + p1.*(eosMDJWFnum(10) + eosMDJWFnum(11).*t2) );
144 den = eosMDJWFden(13) ...
145 + t1.*(eosMDJWFden(1) ...
146 + t1.*(eosMDJWFden(2) ...
147 + t1.*(eosMDJWFden(3) + t1.*eosMDJWFden(4) ) ) ) ...
148 + s1.*(eosMDJWFden(5) ...
149 + t1.*(eosMDJWFden(6) ...
150 + eosMDJWFden(7).*t2) ...
151 + sp5.*(eosMDJWFden(8) + eosMDJWFden(9).*t2) ) ...
152 + p1.*(eosMDJWFden(10) ...
153 + p1t1.*(eosMDJWFden(11).*t2 + eosMDJWFden(12).*p1) );
154
155 epsln = 0;
156 denom = 1.0./(epsln+den) ;
157
158
159 rho = num.*denom;
160
161 return
162

  ViewVC Help
Powered by ViewVC 1.1.22