/[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.1 - (show annotations) (download)
Wed Dec 7 16:59:07 2005 UTC (18 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint58e_post, checkpoint58u_post, checkpoint58r_post, checkpoint57y_post, checkpoint58n_post, checkpoint58t_post, checkpoint58h_post, checkpoint57y_pre, checkpoint58q_post, checkpoint58j_post, checkpoint58, checkpoint58f_post, checkpoint58d_post, checkpoint58c_post, checkpoint58a_post, checkpoint58i_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint58p_post, checkpoint58b_post, checkpoint58m_post
 o add matlab density functions generously provided by Martin Losch in
   this email to the support list:
   http://mitgcm.org/pipermail/mitgcm-support/2005-December/003633.html

1 function rho = densjmd95(s,t,p);
2 %function rho = densjmd95(S,Theta,P);
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
35 %----------------------
36 % CHECK INPUT ARGUMENTS
37 %----------------------
38 if nargin ~=3
39 error('densjmd95.m: Must pass 3 parameters')
40 end
41 if ndims(s) > 2
42 dims = size(s);
43 dimt = size(t);
44 dimp = size(p);
45 if length(dims) ~= length(dimt) | length(dims) ~= length(dimp) ...
46 length(dimt) ~= length(dimp)
47 error(['for more than two dimensions, S, Theta, and P must have the' ...
48 ' same number of dimensions'])
49 else
50 for k=length(dims)
51 if dims(k)~=dimt(k) | dims(k)~=dimp(k) | dimt(k)~=dimp(k)
52 error(['for more than two dimensions, S, Theta, and P must have' ...
53 ' the same dimensions'])
54 end
55 end
56 end
57 else
58 % CHECK S,T,P dimensions and verify consistent
59 [ms,ns] = size(s);
60 [mt,nt] = size(t);
61 [mp,np] = size(p);
62
63 % CHECK THAT S & T HAVE SAME SHAPE
64 if (ms~=mt) | (ns~=nt)
65 error('check_stp: S & T must have same dimensions')
66 end %if
67
68 % CHECK OPTIONAL SHAPES FOR P
69 if mp==1 & np==1 % P is a scalar. Fill to size of S
70 p = p(1)*ones(ms,ns);
71 elseif np==ns & mp==1 % P is row vector with same cols as S
72 p = p( ones(1,ms), : ); % Copy down each column.
73 elseif mp==ms & np==1 % P is column vector
74 p = p( :, ones(1,ns) ); % Copy across each row
75 elseif mp==ms & np==ns % P is a matrix size(S)
76 % shape ok
77 else
78 error('check_stp: P has wrong dimensions')
79 end %if
80 [mp,np] = size(p);
81 % IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
82 Transpose = 0;
83 if mp == 1 % row vector
84 p = p(:);
85 t = t(:);
86 s = s(:);
87 Transpose = 1;
88 end
89 %***check_stp
90 end
91
92 % convert pressure to bar
93 p = .1*p;
94
95 % coefficients nonlinear equation of state in pressure coordinates for
96 % 1. density of fresh water at p = 0
97 eosJMDCFw(1) = 999.842594;
98 eosJMDCFw(2) = 6.793952e-02;
99 eosJMDCFw(3) = - 9.095290e-03;
100 eosJMDCFw(4) = 1.001685e-04;
101 eosJMDCFw(5) = - 1.120083e-06;
102 eosJMDCFw(6) = 6.536332e-09;
103 % 2. density of sea water at p = 0
104 eosJMDCSw(1) = 8.244930e-01;
105 eosJMDCSw(2) = - 4.089900e-03;
106 eosJMDCSw(3) = 7.643800e-05 ;
107 eosJMDCSw(4) = - 8.246700e-07;
108 eosJMDCSw(5) = 5.387500e-09;
109 eosJMDCSw(6) = - 5.724660e-03;
110 eosJMDCSw(7) = 1.022700e-04;
111 eosJMDCSw(8) = - 1.654600e-06;
112 eosJMDCSw(9) = 4.831400e-04;
113
114 t2 = t.*t;
115 t3 = t2.*t;
116 t4 = t3.*t;
117
118 is = find(s(:) < 0 );
119 if ~isempty(is)
120 warning('found negative salinity values, reset them to NaN');
121 s(is) = NaN;
122 end
123 s3o2 = s.*sqrt(s);
124
125 % density of freshwater at the surface
126 rho = eosJMDCFw(1) ...
127 + eosJMDCFw(2)*t ...
128 + eosJMDCFw(3)*t2 ...
129 + eosJMDCFw(4)*t3 ...
130 + eosJMDCFw(5)*t4 ...
131 + eosJMDCFw(6)*t4.*t;
132 % density of sea water at the surface
133 rho = rho ...
134 + s.*( ...
135 eosJMDCSw(1) ...
136 + eosJMDCSw(2)*t ...
137 + eosJMDCSw(3)*t2 ...
138 + eosJMDCSw(4)*t3 ...
139 + eosJMDCSw(5)*t4 ...
140 ) ...
141 + s3o2.*( ...
142 eosJMDCSw(6) ...
143 + eosJMDCSw(7)*t ...
144 + eosJMDCSw(8)*t2 ...
145 ) ...
146 + eosJMDCSw(9)*s.*s;
147
148 rho = rho./(1 - p./bulkmodjmd95(s,t,p));
149
150 if ndims(s) < 3 & Transpose
151 rho = rho';
152 end %if
153
154 return
155
156 function bulkmod = bulkmodjmd95(s,t,p)
157 %function bulkmod = bulkmodjmd95(s,t,p)
158
159 dummy = 0;
160 % coefficients in pressure coordinates for
161 % 3. secant bulk modulus K of fresh water at p = 0
162 eosJMDCKFw(1) = 1.965933e+04;
163 eosJMDCKFw(2) = 1.444304e+02;
164 eosJMDCKFw(3) = - 1.706103e+00;
165 eosJMDCKFw(4) = 9.648704e-03;
166 eosJMDCKFw(5) = - 4.190253e-05;
167 % 4. secant bulk modulus K of sea water at p = 0
168 eosJMDCKSw(1) = 5.284855e+01;
169 eosJMDCKSw(2) = - 3.101089e-01;
170 eosJMDCKSw(3) = 6.283263e-03;
171 eosJMDCKSw(4) = - 5.084188e-05;
172 eosJMDCKSw(5) = 3.886640e-01;
173 eosJMDCKSw(6) = 9.085835e-03;
174 eosJMDCKSw(7) = - 4.619924e-04;
175 % 5. secant bulk modulus K of sea water at p
176 eosJMDCKP( 1) = 3.186519e+00;
177 eosJMDCKP( 2) = 2.212276e-02;
178 eosJMDCKP( 3) = - 2.984642e-04;
179 eosJMDCKP( 4) = 1.956415e-06;
180 eosJMDCKP( 5) = 6.704388e-03;
181 eosJMDCKP( 6) = - 1.847318e-04;
182 eosJMDCKP( 7) = 2.059331e-07;
183 eosJMDCKP( 8) = 1.480266e-04;
184 eosJMDCKP( 9) = 2.102898e-04;
185 eosJMDCKP(10) = - 1.202016e-05;
186 eosJMDCKP(11) = 1.394680e-07;
187 eosJMDCKP(12) = - 2.040237e-06;
188 eosJMDCKP(13) = 6.128773e-08;
189 eosJMDCKP(14) = 6.207323e-10;
190
191 t2 = t.*t;
192 t3 = t2.*t;
193 t4 = t3.*t;
194
195 is = find(s(:) < 0 );
196 if ~isempty(is)
197 warning('found negative salinity values, reset them to NaN');
198 s(is) = NaN;
199 end
200 s3o2 = s.*sqrt(s);
201 %p = pressure(i,j,k,bi,bj)*SItoBar
202 p2 = p.*p;
203 % secant bulk modulus of fresh water at the surface
204 bulkmod = eosJMDCKFw(1) ...
205 + eosJMDCKFw(2)*t ...
206 + eosJMDCKFw(3)*t2 ...
207 + eosJMDCKFw(4)*t3 ...
208 + eosJMDCKFw(5)*t4;
209 % secant bulk modulus of sea water at the surface
210 bulkmod = bulkmod ...
211 + s.*( eosJMDCKSw(1) ...
212 + eosJMDCKSw(2)*t ...
213 + eosJMDCKSw(3)*t2 ...
214 + eosJMDCKSw(4)*t3 ...
215 ) ...
216 + s3o2.*( eosJMDCKSw(5) ...
217 + eosJMDCKSw(6)*t ...
218 + eosJMDCKSw(7)*t2 ...
219 );
220 % secant bulk modulus of sea water at pressure p
221 bulkmod = bulkmod ...
222 + p.*( eosJMDCKP(1) ...
223 + eosJMDCKP(2)*t ...
224 + eosJMDCKP(3)*t2 ...
225 + eosJMDCKP(4)*t3 ...
226 ) ...
227 + p.*s.*( eosJMDCKP(5) ...
228 + eosJMDCKP(6)*t ...
229 + eosJMDCKP(7)*t2 ...
230 ) ...
231 + p.*s3o2*eosJMDCKP(8) ...
232 + p2.*( eosJMDCKP(9) ...
233 + eosJMDCKP(10)*t ...
234 + eosJMDCKP(11)*t2 ...
235 ) ...
236 + p2.*s.*( eosJMDCKP(12) ...
237 + eosJMDCKP(13)*t ...
238 + eosJMDCKP(14)*t2 ...
239 );
240
241 return
242

  ViewVC Help
Powered by ViewVC 1.1.22