function [u,x,y,jd]=kslice(cdf,var,time,layer)
%  KSLICE  returns a matrix containing a horizontal slice at a 
%          specified sigma level (K level) at a given time step
%          from an ECOMxx.CDF or POM.CDF file.  
%
%          This function can also be used
%          to read in 2D and 3D fields such as depth and heat_flux.
%          The coordinates of u are returned as x and y.
%
%       USAGE: [u,x,y,jd]=kslice(cdf,var,[time],[layer])
%
%          where 
%                cdf= file name for netCDf file (e.g. 'pom.cdf')
%                layer = sigma layer (e.g 1 for top layer)
%                var =  the variable to select (eg. 'salt' for salinity)
%                time = time step 
%
%    
%       Examples: 
%
%          [s,x,y]=kslice('ecomsi.cdf','salt',2,3);
%              returns the salinity field from the 3rd sigma level
%              at the 2nd time step.
%
%          [elev,x,y]=kslice('ecomsi.cdf','depth',4);
%              returns the elevation field from the 4th time step
%
%          [depth,x,y]=kslice('ecomsi.cdf','depth');
%              returns the depth field
%
if (nargin<2 | nargin>4),
  help kslice; return
end
% turn off warnings from NetCDf
mexcdf('setopts',0);
%
% open existing file
ncid=mexcdf('open',cdf,'nowrite');
if(ncid==-1),
  disp(['file ' cdf ' not found'])
  return
end
[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
%
% use depth=-99999  to mask points on land
%
depth=mexcdf('varget',ncid,'depth',[0 0],[ny nx]);
land=find(depth==-99999);
%
% allow for using kslice on 2D, 3D and 4D variables
%
if(nargin==4),
 [u,ierr]=mexcdf('varget',ncid,var,[(time-1) layer-1 0 0],[1 1 ny nx],1); 
elseif(nargin==3),
 [u,ierr]=mexcdf('varget',ncid,var,[(time-1) 0 0],[1 ny nx],1); 
else
 [u,ierr]=mexcdf('varget',ncid,var,[0 0],[ny nx],1); 
end
if(nargin>2),
 [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

mexcdf('close',ncid);
u(land)=u(land)*NaN;

