function [U,V]=SSHvel(SSH,lat,lon);
%SSHvel.m-----Computes velocity from sea-surface height data using geostrophic relationship
%
%[U,V]=SSHvel(SSH,lat,lon);
%
%SSH should be an m-by-n array of SSH observations on a grid of points.
%The east-west spacing of the grid is defined by lon (length=n), while the
%north-south spacing is defined by lat (length=m).  

[m,n]=size(SSH);
p=length(lat);
if(p~=m);error('Length of lat must equal the number of rows in SSH');end
p=length(lon);
if(p~=n);error('Length of lon must equal the number of columns in SSH');end

I=1:m-1;J=1:n-1; %index vectors for differences

%Compute dx and dy
%
%
[x,y]=ll2xy(lat(:),ones(m,1)*lon(1),[lat(1),lon(1)]);%ll2xy converts lat-lon to meters from ref.
dy=diff(y);
DY=dy(:)*ones(1,n);
DX=zeros(m,n-1);

for j=1:m;
	%Since distance between longitude lines changes with lat, must
	%compute dx for every lat.
	[x,y]=ll2xy(lat(j)*ones(1,n),lon,[lat(j),lon(1)]);
	DX(j,:)=diff(x(:))';
end

dZdx=diff(SSH,1,2)./DX;
dZdy=diff(SSH,1,1)./DY;
f=2*7.29e-5*sin(pi/180*lat(:));
g=9.8;
U=(-g./f(I)*ones(1,length(J))).*dZdy(I,J);
V=(g./f(I)*ones(1,length(J))).*dZdx(I,J);