function [flux,s,u,area,y,z]=isflux(cdf,var,tstep,i0,jrange,s0);
%
% isflux finds flux of a scalar normal to I= constant in
%   the Blumberg/Mellor model output.  
%
%  Usage:  
%       [flux,s,u,area,y,z]=isflux(cdf,var,tstep,i0,jind,[s0]);
%     iind=indices that indicate extent of J along section i=i0, eg. [5:45]; 
%     
%    averages two lines of "u" so that flux is centered on "tracer" points
%
nc=netcdf(cdf);
sigma=nc{'sigma'}(:);
d=nc{'depth'}(:,i0);
if( exist('jrange') & isfinite(jrange) ),
  d=d(jrange);
else
  jrange=1:length(d);
end
iland=find(d==-99999);
d(iland)=0;
u=nc{'u',1}(tstep,:,jrange,[i0:i0+1]);
u=squeeze(mean(u,3));  % take mean along i dimension
s=squeeze(nc{var,1}(tstep,:,jrange,i0));
s(:,iland)=nan;
u(:,iland)=nan;
%                                                         bracketing v 
e=nc{'elev',1}(tstep,jrange,i0);
                                             %rows of z 
h2=nc{'h2'}(jrange,i0);
%                                                         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.*h2).');
u(end,:)=[];
s(end,:)=[];

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