function [u,v] = uv_from_h_s(lon,lat,h)
% [u,v] = uv_from_h_s(lon,lat,h)
% Calculate velocities from surface height h by geostrophy
% [u,v] will be at staggered grid locations with respect to the h grid
% For h in metres u,v will be m/s
%
% John Wilkin 2 Oct 96

% Indexing convention is that the first index of a matrix is
% latitude (y), and second index is longitude (x),
% i.e. h(lat,lon), lon(lat,lon), lat(lat,lon) etc
% Take care in differencing operations to get this right

% Staggered grid arrangement:
%
%    h------v------h
%    |      |
%    u(1,1)--------h
%    |      |
%    h(1,1)-v(1,1)-h...

global DX_V DY_U GFoDY GFoDX

if isempty(DX_V) | isempty(DY_U)
  % dx at v points and dy at u points
  [DX_V,DY_U] = dxdy_uv(lon,lat);
end

% We use a rectilinear (NOT curvilinear) lat,lon grid, so lat
% is the same for all longitudes, and lon is the same for all
% latitudes.  This let's us cut some corners in the indexing
% of variables

[nlat nlon] = size(h);

if isempty (GFoDY) | isempty(GFoDX)
  % constants
  g = 9.81;
  twoOmega = 2*7.292e-5;

  % coriolis f = 2*omega*sin(lat)
  lat_u = 0.5*(lat(2:nlat,:)+lat(1:nlat-1,:));
  f = twoOmega*sin(pi/180*lat_u);
  f_u = f(:,ones([1 nlon]));

  % factor to multiply diff(h) to get u component
  GFoDY = g./f_u./DY_U;
  
  lat_v = lat(:,1);
  f = twoOmega*sin(pi/180*lat_v);
  f_v = f(:,ones([1 nlon-1]));
  
  % factor to multiply diff(h) to get v component
  GFoDX = g./f_v./DX_V;
end
  
% u component
% coriolis f = 2*omega*sin(lat)
u = -GFoDY.*(h(2:nlat,:)-h(1:nlat-1,:));

% v component
v = GFoDX.*(h(:,2:nlon)-h(:,1:nlon-1));
