function [n2, e] = bm_bvfrq( S, T, P, LAT, DP)
%*************************************************************
%                ***** brunt-vaisala freq *****
%
%  uses 1980 equation of state
%
% in :
%       P   pressure         decibars
%       T   temperature      deg celsius (ipts-68)
%       S   salinity         psu (pss-78)
%       LAT latitude
%       DP  intervalle de pression (db) sur lequel est calcule n2
%           
% out :
%       n2  bouyancy freq    cph
%       e   stability        1/m
% 
% r. millard   feb. 1991
%
%*************************************************************
%------
% BEGIN
%------

% change call list xlat to glat & pass gravity m/s^2
% after formulation of n. p. fofonoff & breck owen's
%

% Calcul de la gravit�
% note that sw_g expects height as argument
Z  = sw_dpth(P,LAT);
gp = sw_g(LAT,-Z);
g2 = gp .* gp;
g2 = g2 * 1.e-4;

% compute least squares estimate of specific volume anamoly gradient

[m,n] = size(P);
n2    = zeros(size(P));

% Calcul du nombre d'observation sur lequel est calcule n2

dp = P(2)-P(1);
nobs = DP / dp;

% Debut du profil

for i = 1:(nobs/2)
  j = 1:i+nobs/2;
  [n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end

% Fin du profil

for i = m-(nobs/2):m
  j = i-nobs/2:m;
  [n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end

% Profil

for i = (nobs/2+1):m-(nobs/2+1)
  j = i-nobs/2:i+nobs/2;
  [n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end

% Corps du calcul

function [n2, e] = cal_bvfrq( S, T, P, g2, gp, nobs, j)
  
%  default gravity & rad/sec to cph conversion
%
radsec = 572.9578;

  pav  = mean(P(j));
  data = sw_svan( S(j), sw_ptmp( S(j), T(j), P(j), pav), pav);
  cxy  = sum( data .* (P(j)-pav) );
  cy   = sum( data );
  cxx  = sum( (P(j)-pav).*(P(j)-pav) );
  
  sigma = sw_pden(S(j), sw_ptmp( S(j), T(j), P(j), pav), P(j), pav);

  a0    = cxy/cxx;

  v350p = 1./sigma - data;
  vbar  = v350p(end) + cy/(nobs+1);
  dvdp  = a0;

  e     = -g2*dvdp/(vbar)^2;
  n2    = radsec*sqrt(abs(e))*sign(e);

  %  define stability parameter units (1/m)
  e = e/gp;

  return