matlab-grid-generator/bin/n2_poly3.m

function [n2] = dens_poly3(poly3,t,s,dzc)
% D=DENS_POLY3(P,T,S,dzc)
%
% Calculates Brunt-Vasala frequecy 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);
% >> N2=n2_poly3(P,T,S);
%
% or to work within a single model level
% >> N2=n2_poly3(P,T(:,:,3),S(:,:,3),3);

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]);

n2=0*t;
for k=2:nz,
 d1=dens_poly3(poly3(k-1),t(:,k-1),s(:,k-1));
 d2=dens_poly3(poly3(k-1),t(:,k  ),s(:,k  ));
 dRdz=(d1-d2)/dzc(k);
 dRdz( find( t(:,k)==0 ) )=0;
 n2(:,k)=-9.81*dRdz;
end

n2=reshape(n2,n);
1	cnh	1.1	function [n2] = dens_poly3(poly3,t,s,dzc)
2			% D=DENS_POLY3(P,T,S,dzc)
3			%
4			% Calculates Brunt-Vasala frequecy 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			% >> N2=n2_poly3(P,T,S);
15			%
16			% or to work within a single model level
17			% >> N2=n2_poly3(P,T(:,:,3),S(:,:,3),3);
18
19			if size(t) ~= size(s)
20			error('T and S must be the same shape and size')
21			end
22			%if size(t,ndims(t)) ~= size(poly3,2)
23			% error('Last dimension of T and S must be the number of levels in P')
24			%end
25
26			n=size(t);
27			nz=size(poly3,2);
28
29			t=reshape(t,[prod(size(t))/nz nz]);
30			s=reshape(s,[prod(size(t))/nz nz]);
31
32			n2=0*t;
33			for k=2:nz,
34			d1=dens_poly3(poly3(k-1),t(:,k-1),s(:,k-1));
35			d2=dens_poly3(poly3(k-1),t(:,k ),s(:,k ));
36			dRdz=(d1-d2)/dzc(k);
37			dRdz( find( t(:,k)==0 ) )=0;
38			n2(:,k)=-9.81*dRdz;
39			end
40
41			n2=reshape(n2,n);