function [u3d,x,y,jd]=zscrum(cdf,var_name,time,zlevel)
%  ZSLICE  Returns a matrix containing a horizontal slice at a specified depth
%           at a given time step from an ECOMxx.CDF or POM.CDF file.
%           Regions of grid that are shallower than requested value are
%           returned as NaNs.
%
%       USAGE: [u,x,y,jd]=zslice(cdf,var,time,zlevel)
%
%   where zlevel is a depth in meters (e.g -10.)
%
% see also ZSLICEUV, ISLICE, JSLICE, KSLICE, KSLICEUV
% hint: use PSLICE or EXTCONTOUR to plot the results of ZSLICE 

if (nargin<2 | nargin>4),
  help zslice; return
end
ncquiet
nc=netcdf(cdf,'nowrite');
if isempty(nc),
  return
end

% u is an Netcdf object
u=nc{var_name};

if isempty(u),
  close(nc);
  return
end

%bathymetry at RHO points
h=nc{'h'};

% elevation 
zeta=nc{'zeta'};

% S coordinate
Cs_r=nc{'Cs_r'};

% h and zeta are objects (not matrices),
% so (:) gets all the data in right shape

z=(h(:)+zeta(:));

% get all values of var at time "time"
u3d=u(time,:,:,:);

% squeeze out singleton dimensions
u3d=squeeze(u3d);

close(nc);
return
ar
if(nargout>2),
  x=mexcdf('varget',ncid,'x',[0 0],[ny nx]);
  y=mexcdf('varget',ncid,'y',[0 0],[ny nx]);
end
depth=mexcdf('varget',ncid,'depth',[0 0],[ny nx]);
elev=mexcdf('varget',ncid,'elev',[(time-1) 0 0],[1 ny nx],1);
dind=find(depth==-99999);
% use total depth (depth+elev) for sigma model
 depth=depth+elev;
sigma=mexcdf('varget',ncid,'sigma',0,nz);
n=length(sigma);
sigma2=sigma(1:n-1)+0.5*diff(sigma);
sigma2(n)=-1;
uz=mexcdf('varget',ncid,var,[(time-1) 0 0 0],[1 nz ny nx],1); %profile

% get time stuff
if(nargout>3)
[nam,nt]=mexcdf('diminq',ncid,'time');
base_date=[0 0 0 0 0 0];
base_date(1:3)=mexcdf('attget',ncid,'global','base_date');
t=mexcdf('varget',ncid,'time',time-1,1);
jd0=julian(base_date);
jd=jd0+t;
end

%close file
mexcdf('close',ncid);

uz=reshape(uz,nx*ny,n)';

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));

u=zeros(nx*ny,1);
% do interpolation
u(isum)=uz(ik)+(uz(ik+1)-uz(ik)).*(zlevel-zlev(ik))./(zlev(ik+1)-zlev(ik));
u=reshape(u,nx,ny);

dind=find(depth< -zlevel);
u(dind)=u(dind)*NaN;

