function [w,x,y]=roms_uv100(gname,fname,timestep)
%  ROMS_UV100  Returns a horizontal slice of velocity 1m above bottom
%        at a given time step from a ROMS NetCDF file
%
%       USAGE: [w,x,y]=roms_uv100(gname,fname,timestep);
%           where var1 and var2 are names of u,v 3d fields  ('u','v')
%

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

if (nargin~=3)
  help roms_uv100; return
end
var1='u';
var2='v';
klevel=1;
ncmex('setopts',0);
ncid=ncmex('open',fname,'nowrite');
ncid2=ncmex('open',gname,'nowrite');
if(ncid==-1),
  disp(['file ' fname ' 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


[data_u, status] = ncmex ( 'varget', ncid, var1, [timestep-1 klevel-1 0 0], [1 1  -1 -1],1);
if ( status == -1 )
	fprintf ( 'roms_depave_uv:  could not get ''var1'' in %s.', fname );
	return;
end
data_u = 0.5*(data_u(1:(end-1),:) + data_u(2:end,:));

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


%
% 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_uv100:  could not get ''mask_rho'' in %s.', fname );
		return;
	end
end
[vid, status] = ncmex ( 'varid', ncid, 'angle' );
if ( status ~= -1 )
	[ang, status] = ncmex ( 'varget', ncid, vid, [0 0], [-1 -1] );
	if ( status == -1 )
		fprintf ( 'roms_uv100:  could not get ''angle'' in %s.', fname );
		return;
	end
end
% try to get Apparent Bottom roughness
[vid, status] = ncmex ( 'varid', ncid, 'Zoa' );
if ( status ~= -1 )
	[Zob, status] = ncmex ( 'varget', ncid, vid, [0 0], [-1 -1] );
	if ( status == -1 )
		fprintf ( 'roms_uv100:  could not get ''Zoa'' in %s.', fname );
		return;
	end
else
% but if doesn't exist, use Constant bottom roughness
  [vid, status] = ncmex ( 'varid', ncid, 'Zob' );
  Zob=ncmex('varget1',ncid,'Zob',[0],[1]);
		fprintf ( 'roms_uv100:  using Zob for Zoa %s.', fname );
end

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

w=nan*ones(size(x));

% center w on rho points
w(2:(end-1),2:(end-1))=data_u(:,2:(end-1)) + sqrt(-1)*data_v(2:(end-1),:);
% rotate vectors
w=w.*exp(sqrt(-1)*ang);

% compute angle from lon/lat
%
% [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


% adjust bottom level current to 1 m above bottom using law-of-wall
z1=depths(fname,fname,1,0,timestep);
z2=depths(fname,fname,5,0,timestep);
zr=squeeze(z1(:,:,1)-z2(:,:,1));  % thickness of bottom cell
cd=z0tocd(Zob,zr);
ustar=sqrt(cd).*abs(w);
ur = (ustar/0.4).*log(ones(size(zr))./Zob);
w=w.*(ur)./(abs(w)+eps);

% mask w where there is land
w(~mask_rho)=NaN;