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

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

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


Revision 1.2 - (hide annotations) (download)
Sat Feb 17 23:49:43 2007 UTC (14 years, 5 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 edhill 1.1 function rho = densjmd95(s,t,p);
2    
3 jmc 1.2 %
4 edhill 1.1 % DENSJMD95 Density of sea water
5     %=========================================================================
6 jmc 1.2 %
7 edhill 1.1 % 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 jmc 1.2
30 edhill 1.1
31     % Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388
32    
33     % created by mlosch on 2002-08-09
34 jmc 1.2 % $Header: $
35     % $Name: $
36 edhill 1.1
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