function [flux,s,v,area,x,z]=jsflux(cdf,var,tstep,j0,irange,s0);
%
% JSFLUX finds flux of a scalar normal to J= constant in
%   the Blumberg/Mellor model output.  
%
%  Usage:  
%       [flux,s,v,area,x,z]=jsflux(cdf,var,tstep,j0,iind,[s0]);
%     iind=indices that indicate extent of I along section j=j0, eg. [5:45]; 
%     
%    averages two lines of "v" so that flux is centered on "tracer" points
%
nc=netcdf(cdf);
sigma=nc{'sigma'}(:);
d=nc{'depth'}(j0,:);
if( exist('irange') & isfinite(irange) ),
  d=d(irange);
else
  irange=1:length(d);
end
iland=find(d==-99999);
d(iland)=0;
v=nc{'v',1}(tstep,:,[j0:j0+1],irange);
v=squeeze(mean(v,2));  % take mean along j dimension
s=squeeze(nc{var,1}(tstep,:,j0,irange));
s(:,iland)=nan;
v(:,iland)=nan;
%                                                         bracketing v 
e=nc{'elev',1}(tstep,j0,irange);
                                                               %rows of z 
h1=nc{'h1'}(j0,irange);
%                                                         h1 bracketing v
dz=-diff(sigma);
sigma=.5*(sigma(1:end-1)+sigma(2:end));
td=d+e;  %average elevation and depth to get total depth at "var" points
z=sigma*td;
%
% area of each grid cell face in section
area=dz*(d.*h1);
v(end,:)=[];
s(end,:)=[];

%
% if optional argument, compute "freshwater" flux
if(exist('s0')),
  s=max(0,(s0-s)/s0);
end
flux=s.*v.*area;
%
x=cumsum(h1)/1000;  % convert m to km 
x=ones(size(sigma))*x;
