utils/matlab/dens_poly3.m

function [dens] = dens_poly3(poly3,t,s)
% D=DENS_POLY3(P,T,S)
%
% Calculates in-situ density as approximated by the POLY3 method
% used in the MITgcm.
%  P - coefficients read from file 'POLY3.COEFFS' using INI_POLY3
%  T - potential temperature
%  S - salinity
%
% eg.
% >> P=ini_poly3;
% >> T=rdmds('T',100);
% >> S=rdmds('S',100);
% >> D=dens_poly3(P,T,S);
%
% or to work within a single model level
% >> D=dens_poly3(P(3),T(:,:,3),S(:,:,3));

% $Header:  $
% $Name:  $

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

n=size(t);
nz=size(poly3,2);

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

for k=1:nz,
 tRef=poly3(k).t;
 sRef=poly3(k).s;
 dRef=poly3(k).dens;
 tp=t(:,k)-tRef;
 sp=s(:,k)-sRef;

 deltaSig=                                    ...
  poly3(k).coeffs(1) .*tp                     ...
 +poly3(k).coeffs(2)             .*sp         ...
 +poly3(k).coeffs(3) .*tp.*tp                 ...
 +poly3(k).coeffs(4) .*tp        .*sp         ...
 +poly3(k).coeffs(5)             .*sp.*sp     ...
 +poly3(k).coeffs(6) .*tp.*tp.*tp             ...
 +poly3(k).coeffs(7) .*tp.*tp    .*sp         ...
 +poly3(k).coeffs(8) .*tp        .*sp.*sp     ...
 +poly3(k).coeffs(9)             .*sp.*sp.*sp ...
  ;
 dens(:,k)=deltaSig+dRef;
end

dens=reshape(dens,n);
dens( find(t==0 & s==0) )=0;
1	adcroft	1.1	function [dens] = dens_poly3(poly3,t,s)
2			% D=DENS_POLY3(P,T,S)
3			%
4			% Calculates in-situ 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	jmc	1.3	% $Header: $
20			% $Name: $
21
22	adcroft	1.1	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	adcroft	1.2	dens( find(t==0 & s==0) )=0;