function [w,x,y]=roms_uv2d(gridfile,cdf,var1,var2,timestep)
%  ROMS_UV2D  Read a ROMS 2D vector field at a particular time step 
%             
%   USAGE: [w,x,y]=roms_uv2d(gfile,ffile,var1,var2,timestep);
%
%   w = complex vector field (u+sqrt(-1)*v)
%   x = lon or x
%   y = lat or y
%   var1 = u component of 2d field (e.g. 'ubar')
%   var2 = v component of 2d field (e.g. 'vbar');
%   gfile = grid file name  (e.g. 'ocean_grd.nc')
%   ffile = field file name (e.g. 'ocean_his.nc', 'ocean_avg.nc', 'ocean_rst.nc')
%   
%   Example: [w,lon,lat]=roms_uv2d('ocean_grd.nc','ocean_his.nc','ubar','vbar',10);
%     read the depth-averaged velocity at time step 10 from file 'ocean_his.nc');

% John Evans & Rich Signell (rsignell@usgs.gov);

if (nargin<4 | nargin>5),
  help roms_uv2d; return
end

ncmex('setopts',0);
ncid=ncmex('open',cdf,'nowrite');
ncid2=ncmex('open',gridfile,'nowrite');
if(ncid==-1),
  disp(['file ' cdf ' not found'])
  return
end

%
% Acquire the grid.
%
% If "lon_rho" and "lat_rho" are present, grab them.
% Otherwise, get "x_rho" and "y_rho".
[lon_rho_varid, rcode] = ncmex('VARID', ncid, 'lon_rho');
[lat_rho_varid, rcode] = ncmex('VARID', ncid, 'lat_rho');
if ( (lon_rho_varid >= 0) | (lat_rho_varid >= 0) )
    x=ncmex('varget',ncid,'lon_rho',[0 0],[-1 -1]);
    y=ncmex('varget',ncid,'lat_rho',[0 0],[-1 -1]);
else
    x=ncmex('varget',ncid,'x_rho',[0 0],[-1 -1]);
    y=ncmex('varget',ncid,'y_rho',[0 0],[-1 -1]);
end

% Retrieve u at the given timestep.
data_u=nan*x;
[utmp, status] = ncmex ( 'varget', ncid, var1, [timestep-1 0 0], [1  -1 -1] );
if ( status == -1 )
	fprintf ( 'roms_depave_uv:  could not get ''var1'' in %s.', cdf );
	return;
end

data_u(2:end-1,2:end-1) = (utmp(1:end-1,2:end-1) + utmp(2:end,2:end-1))/2;

%
% Retrieve v at the given timestep.
data_v = nan*x;
[vtmp, status] = ncmex ( 'varget', ncid, var2, [timestep-1 0 0], [1 -1 -1] );
if ( status == -1 )
	fprintf ( 'roms_depave_uv:  could not get ''var2'' in %s.', cdf );
	return;
end
data_v(2:end-1,2:end-1) = (vtmp(2:end-1,1:end-1) + vtmp(2:end-1,2:end))/2;


%
% If there is a mask_rho, we will want to use it later.
[mask_rho_varid, status] = ncmex ( 'varid', ncid, 'mask_rho' );
if ( status ~= -1 )
	[mask_rho, status] = ncmex ( 'varget', ncid, 'mask_rho', [0 0], [-1 -1] );
	if ( status == -1 )
		fprintf ( 'roms_depave_uv:  could not get ''mask_rho'' in %s.', cdf );
		return;
	end
end

angle=ncmex('varget',ncid2,'angle',[0 0],[-1 -1]);

ncmex ( 'close', ncid );
ncmex ( 'close', ncid2 );

w = data_u + sqrt(-1)*data_v;

%w = w.'; x = x.'; y = y.';

% apply the rotation
%
% Since the "angle" variable is not always present, we
% need to construct it.  In the roms files, it is apparently
% always in degrees.  Don't bother with degrees here.
%angle = zeros(size(x));
%
%
%
%[r,c] = size(x);
%j = [2:c-1];
%for i = 2:r-1
%    angle(i,j) = atan2(y(i+1,j)-y(i-1,j), x(i+1,j)-x(i-1,j));
%end

% rotate into east and north components
 w=w.*exp(sqrt(-1)*angle);

% Mask out the land.
w(~mask_rho) = NaN ;
w=w.';
x=x.';
y=y.';