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

  ViewVC Help
Powered by ViewVC 1.1.22