function div = roms_divergence(u,v,g)
% Simple and approximate calculation of the divergence of a ROMS
% vector field
%
% div = roms_divergence(u,v,grd)
%
% Inputs
%       u is the xi-direction component on the c-grid u points
%       v is the eta-direction component on the c-grid v points
%
% Optional inputs
%       grd is a ROMS grid structure from function roms_get_grid
%           (pm,pn metrics are used in calculation)
%
%           If this argument is omitted, it is assumed dx=dy=1;
%
% John Wilkin Nov 17, 2005

if nargin == 3
  dudx = (u(:,2:end)-u(:,1:end-1)).*g.pm(:,2:end-1);
  dvdy = (v(2:end,:)-v(1:end-1,:)).*g.pn(2:end-1,:);
else
  dudx = (u(:,2:end)-u(:,1:end-1));
  dvdy = (v(2:end,:)-v(1:end-1,:));
  % dudx = zeros(size(dudx));
end

% There is some noisiness to this calculation and it needs more checking
if 1
  dudx = dudx(:,[1 1:end end]);
  dvdy = dvdy([1 1:end end],:);
  div = dudx + dvdy;
end

div(:,[1:2 end-1 end]) = 0;
div([1 2 end-1 end],:) = 0;

if 0
  % disable the original way of doing this
  div = zeros([size(u,1) size(v,2)]);
  div(2:end-1,2:end-1) = dudx(2:end-1,:)+dvdy(:,2:end-1);
  div(:,[1 end]) = div(:,[2 end-1]);
  div([1 end],:) = div([2 end-1],:);
end
