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

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

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


Revision 1.2 - (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.1: +5 -3 lines
add $Header:  $ and $Name:  & (for CVS)

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: $
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

  ViewVC Help
Powered by ViewVC 1.1.22