function nsq = bvf(S, T, P, lat)
% function nsq = bvf(S, T, P, lat)
% Brunt-Vaisala frequency squared; calculated using centered differences.

if nargin < 3, help('bvf'); error('Not enough inputs.'); end
if nargin == 3, lat = 30; end

[mS, nS] = size(S);
[mT, nT] = size(T);
[mP, nP] = size(P);

% Put in all the shape checking later; maybe even a subroutine.

m = mS; n = nS;
ii0 = 1:(m-2);
ii1 = 3:m;
iim = 2:(m-1);

% At this point, S, T, P must be column vectors or matrices with
% pressure varying down the column.
nsq = NaN * zeros(m,n);
[Z, g] =  depth(P, lat); %%%New version of depth.m to output gravity.

[sv, sig0] = svan(S(ii0,:), theta(S(ii0,:), T(ii0,:), P(ii0,:),P(iim,:)), P(iim,:) );
[sv, sig1] = svan(S(ii1,:), theta(S(ii1,:), T(ii1,:), P(ii1,:),P(iim,:)), P(iim,:) );
sig  = 0.5 * (sig0 + sig1);
rho = 1000 + sig;

nsq(iim,:) =  (g(iim,:)./rho) .*(sig1 - sig0) ./ (Z(ii1,:) - Z(ii0,:));

