function nsq = bvfsq(zvec);
%
% returns Brunt-Vaisala frequency squared as a function of z
% It is assumed that zvec is negative.
% JVM fiddled to have salinity as well as temperature.

[m n] = size(zvec);

if n == 1
  z = zvec;
else
  z = zvec';
end

rho0 = 1000.;
g    = 9.81;

% Specify the temperature and salinity profiles.

T = 17.4 + z.*( 0.01727 + z.*( 5.707e-06 + z*5.921e-10 ));
S = 35.08 + z.*( 0.000756 + z.*( 3.185e-07 + z*3.702e-11 ));

% Calculate the density.  Note that Phil Morgan's routine density
% gives the same result as Nathan's sw_dens.

%rhobar  = density(-z, T, S);
rhobar  = sw_dens(S, T, -z);

% Calculate the Brunt Vaisalla frequency.

nsq = -g/rho0*diff(rhobar)./diff(z);
nsq = [ nsq(1); nsq];


