function [u,x,y,jd]=kslice(cdf,var,time,layer)
%  KSLICE  Extracts a horizontal slice of a scalar variable 
%          at a given sigma level from a 4D ECOM NetCDF file. 
%          This function can also be used
%          to read in 2D and 3D fields such as depth and heat_flux.
%
%  USAGE: [u,x,y,jd]=kslice(cdf,var,[time],[layer])
%    cdf = NetCDF file name (string) (e.g. 'ecomsi.cdf')
%    var = variable name (string) (e.g. 'salt' for salinity)
%    time = time step (integer) (e.g. 1 for 1st time step)
%    layer = sigma layer to slice (e.g 1 for top layer)
%
%    u = 2D array of the variable "var"
%    x = 2D array of x grid locations
%    y = 2D array of y grid locations
%   jd = Julian time of slice (type "help julian" for details)
%
%   Example:  to extract temperature from the 1st sigma layer 
%             from the 3rd time step in the file 'ecom3d.cdf', do:
%
%          [temp,x,y]=kslice('ecom3d.cdf','temp',3,1);
%
%             or to extract the depth array, do:
%
%          [depth,x,y]=kslice('ecom3d.cdf','depth');
%
%   Tip: Visualize these arrays using PSLICE or EXTCONTOUR,
%       [e.g. PSLICE(x,y,temp); or EXTCONTOUR(x,y,temp,'label','fill');]

% Rich Signell rsignell@usgs.gov

if (nargin<2 | nargin>4),
  help zslice; 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>1),
  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]);
if(nargin==4),
 u=mexcdf('varget',ncid,var,[(time-1) layer-1 0 0],[1 1 ny nx],1); 
elseif(nargin==3),
 u=mexcdf('varget',ncid,var,[(time-1) 0 0],[1 ny nx],1); 
else
 u=mexcdf('varget',ncid,var,[0 0],[ny nx],1); 
end
% 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

mexcdf('close',ncid);
dind=find(depth==-99999 | depth==0);
u(dind)=u(dind)*NaN;

