1 
function [dens] = dens_poly3(poly3,t,s) 
2 
% D=DENS_POLY3(P,T,S) 
3 
% 
4 
% Calculates insitu density as approximated by the POLY3 method 
5 
% used in the MITgcm. 
6 
% P  coefficients read from file 'POLY3.COEFFS' using INI_POLY3 
7 
% T  potential temperature 
8 
% S  salinity 
9 
% 
10 
% eg. 
11 
% >> P=ini_poly3; 
12 
% >> T=rdmds('T',100); 
13 
% >> S=rdmds('S',100); 
14 
% >> D=dens_poly3(P,T,S); 
15 
% 
16 
% or to work within a single model level 
17 
% >> D=dens_poly3(P(3),T(:,:,3),S(:,:,3)); 
18 

19 
% $Header: $ 
20 
% $Name: $ 
21 

22 
if size(t) ~= size(s) 
23 
error('T and S must be the same shape and size') 
24 
end 
25 
%if size(t,ndims(t)) ~= size(poly3,2) 
26 
% error('Last dimension of T and S must be the number of levels in P') 
27 
%end 
28 

29 
n=size(t); 
30 
nz=size(poly3,2); 
31 

32 
t=reshape(t,[prod(size(t))/nz nz]); 
33 
s=reshape(s,[prod(size(t))/nz nz]); 
34 

35 
for k=1:nz, 
36 
tRef=poly3(k).t; 
37 
sRef=poly3(k).s; 
38 
dRef=poly3(k).dens; 
39 
tp=t(:,k)tRef; 
40 
sp=s(:,k)sRef; 
41 

42 
deltaSig= ... 
43 
poly3(k).coeffs(1) .*tp ... 
44 
+poly3(k).coeffs(2) .*sp ... 
45 
+poly3(k).coeffs(3) .*tp.*tp ... 
46 
+poly3(k).coeffs(4) .*tp .*sp ... 
47 
+poly3(k).coeffs(5) .*sp.*sp ... 
48 
+poly3(k).coeffs(6) .*tp.*tp.*tp ... 
49 
+poly3(k).coeffs(7) .*tp.*tp .*sp ... 
50 
+poly3(k).coeffs(8) .*tp .*sp.*sp ... 
51 
+poly3(k).coeffs(9) .*sp.*sp.*sp ... 
52 
; 
53 
dens(:,k)=deltaSig+dRef; 
54 
end 
55 

56 
dens=reshape(dens,n); 
57 
dens( find(t==0 & s==0) )=0; 