function [w,x,y]=zsliceuv(cdf,time,zlevel)
%  ZSLICEuv 
%           horizontal slice of velocity at a specified depth
%           at a given time step from an ECOMSI.CDF file
%
%       USAGE: [u,x,y]=zsliceuv(cdf,time,zlevel)
%
%   cdf = string containing file name (ex. 
%   where zlevel is a depth in meters (e.g -10.)
%
if (nargin<2 | nargin>3),
  help zsliceuv; return
end
mexcdf('setopts',0);
ncid=mexcdf('open',cdf,'nowrite');

[name, nx]=mexcdf('diminq',ncid,'xpos');
[name, ny]=mexcdf('diminq',ncid,'ypos');
[name, nz]=mexcdf('diminq',ncid,'zpos');
if(nargout==3),
  x=mexcdf('varget',ncid,'x',[0 0],[ny nx]);
  y=mexcdf('varget',ncid,'y',[0 0],[ny nx]);
end
sigma=mexcdf('varget',ncid,'sigma',0,nz);
depth=mexcdf('varget',ncid,'depth',[0 0],[ny nx]);
ang=mexcdf('varget',ncid,'ang',[0 0],[ny nx]);
n=length(sigma);
sigma2=sigma(1:n-1)+0.5*diff(sigma);
sigma2(n)=-1;
uz=mexcdf('varget',ncid,'u',[(time-1) 0 0 0],[1 nz ny nx],1); %profile
vz=mexcdf('varget',ncid,'v',[(time-1) 0 0 0],[1 nz ny nx],1); %profile
mexcdf('close',ncid);

nxy=nx*ny;
nzy=nz*ny;
uz=uz(1:nx-1,:)+.5*diff(uz);
uz(nx,:)=uz(nx-1,:);
vz=vz(:,1:nzy-1)+.5*diff(vz')';
vz(:,nzy)=vz(:,nzy-1);
wz=uz+sqrt(-1)*vz;
wz=reshape(wz,nx*ny,n).';
 
dind=find(depth==-99999);
depth(dind)=depth(dind)*0;
zlev=sigma*(depth(:).');

zind=zlev<zlevel;
dz=diff(zind);
dz(n,:)=zeros(size(dz(n-1,:)));
ik=find(dz==1);

isum=find(sum(dz));

w=zeros(nx*ny,1);
% do interpolation
w(isum)=wz(ik)+(wz(ik+1)-wz(ik)).*(zlevel-zlev(ik))./(zlev(ik+1)-zlev(ik));
w=reshape(w,nx,ny);
% 
dind=find(depth< -zlevel);
w(dind)=w(dind)*NaN;
%
% rotate into true east, north coordinates
%
w=w.*exp(sqrt(-1)*ang);
