function [u,x,z,jd] =jslice(cdf,var,time,jindex)
%  JSLICE  Extracts a vertical slice of a scalar variable 
%          along row J from a 4D ECOM NetCDF file. 
%
%  USAGE: [u,x,z,jd]=jslice(cdf,var,time,j_index)
%    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)
%    j_index = J index along which to take vertical slice
%
%    u = 2D array of the variable "var"
%    x = distance along slice in km (assuming units in netCDF file are m)
%    z = depth in m
%   jd = Julian time of slice (type "help julian" for details)
%
%
%   Example:  to extract salinity along J=10
%             from the 3rd time step in the file 'ecom3d.cdf', do:
%
%          [salt,x,z]=jslice('ecom3d.cdf','salt',3,10);
%
%   Tip: Visualize these arrays using EXTCONTOUR,
%       [e.g. EXTCONTOUR(x,z,temp,'label','fill');]

% Rich Signell rsignell@usgs.gov

if nargin~=4 ,
  help jslice; return
end
mexcdf('setopts',0);
ncid=mexcdf('open',cdf,'nowrite');
[name, nx]=mexcdf('diminq',ncid,'xpos');
[name, nz]=mexcdf('diminq',ncid,'zpos');
  u=mexcdf('varget',ncid,var,[(time-1) 0 (jindex-1) 2],[1 nz-1 1 nx-4],1);
if nargout>1,
  d=mexcdf('varget',ncid,'depth',[(jindex-1) 2],[1 nx-4],1);
  ind=find(d==-99999.);
  if(~isempty(ind)),
   d(ind)=d(ind)*nan;
   u(ind,:)=u(ind,:)*NaN;
  end
  sigma=mexcdf('varget',ncid,'sigma',0,nz,1);
% assume we are slicing a variable defined at 
% at the center of the sigma level (t,s,u,v,km)
  sigma=0.5*(sigma(1:nz-1)+sigma(2:nz));
  h1=mexcdf('varget',ncid,'h1',[(jindex-1) 0],[1 nx],1);
  x=cumsum(h1)/1000.;
  x=x(3:nx-2);
  x=x*ones(1,nz-1);
  z=d*(sigma');
  ind=find(isnan(z));
  z(ind)=zeros(size(ind));
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);
