function [w,x,y]=scrum_zsliceuv(cdf,timestep)
%  SCRUM_ZSLICE  Returns a matrix containing a horizontal slice of velocity
%          at a specified depth at a given time step from a SCRUM NetCDF 
%          file.  Regions of grid that are shallower than requested value 
%          are returned as NaNs.
%
%       USAGE: [w,x,y]=scrum_zsliceuv(cdf,timestep)
%

% John Evans (joevans@usgs.gov)

if (nargin<2 | nargin>4),
  help scrum_zslice; return
end

ncmex('setopts',0);
ncid=ncmex('open',cdf,'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

x_rho = x';
y_rho = y';


[eta_rho_length, xi_rho_length] = size(x_rho);
eta_u_length = eta_rho_length;
eta_v_length = eta_rho_length-1;
xi_u_length = xi_rho_length-1;
xi_v_length = xi_rho_length;
eta_psi_length = eta_rho_length-1;
xi_psi_length = xi_rho_length-1;


%
% construct the grids at the psi points.
xtemp = (x_rho(:,1:(xi_rho_length-1)) + x_rho(:,2:xi_rho_length))/2;
x = (xtemp(1:(eta_rho_length-1),:)  + xtemp(2:eta_rho_length,:))/2;
ytemp = (y_rho(:,1:(xi_rho_length-1)) + y_rho(:,2:xi_rho_length))/2;
y = (ytemp(1:(eta_rho_length-1),:)  + ytemp(2:eta_rho_length,:))/2;


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

%
% Retrieve v at the given timestep.
if ( status == -1 )
	fprintf ( 'scrum_zsliceuv:  could not get ''u'' in %s.', cdf );
	return;
end
[data_v, status] = ncmex ( 'varget', ncid, 'svstr', [timestep 0 0], [1 -1 -1] );
data_u = permute ( data_u, [2 1] );
data_u = (data_u(1:eta_u_length-1,:) + data_u(2:eta_u_length,:))/2;
data_v = permute ( data_v, [2 1] );
data_v = (data_v(:,1:xi_v_length-1) + data_v(:,2:xi_v_length))/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 ( 'scrum_zsliceuv:  could not get ''mask_rho'' in %s.', cdf );
		return;
	end
end

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

ncmex ( 'close', ncid );



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

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

[r,c] = size(x);

% calculate the angle if not present in the netcdf file
%
%angle = zeros(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(1:r,1:c));

%
% Mask out the land.
% Since we average the u and v across cells, the appropriate mask
% dimension is that of "mask_psi".  I don't want to take the chance
% that it may not be present in the file, so I compute it from 
% mask_rho.
[umask,vmask,pmask]=uvp_masks ( mask_rho );
mask_inds = find ( pmask == 0 );
w(mask_inds) = NaN * ones(size(mask_inds));