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,y,z]=isflux(cdf,var,tstep,i0,jrange);
%     jrange=extents of flux calculation along i=i0, eg. [5 45]; 
%     

ncid=mexcdf('open',cdf);
[dim_name,y_size,rcode]=mexcdf('diminq',ncid,'ypos');
[dim_name,k_size,rcode]=mexcdf('diminq',ncid,'zpos');
kbm=k_size-1;
if(~exist('jrange')),
  jmin=0;
  jm=y_size;
else
  jmin=jrange(1)-1;
  jm=min(y_size,jrange(2)-jrange(1)+1);
end
sigma=mexcdf('varget',ncid,'sigma',0,kbm+1);
sigma=sigma';
u=mexcdf('varget',ncid,'u',[tstep-1 0 jmin i0-1],[1 kbm jm 1],1);
s=mexcdf('varget',ncid,var,[tstep-1 0 jmin i0-2],[1 kbm jm 2],1);
%
s=mean(s);  %average scalar at j and j-1 to get value at u(j) points.
%
s=reshape(s,jm,kbm);
d=mexcdf('varget',ncid,'depth',[jmin i0-2],[jm 2]); %get 2 rows of depth 
% after averaging
iland=find(mean(d)<-11111);  
s(iland,:)=s(iland,:)*nan;
%                                                         bracketing u 
dind=find(d==-99999);
d(dind)=zeros(size(dind));
%
e=mexcdf('varget',ncid,'elev',[tstep-1 jmin i0-2],[1 jm 2],1); %get 2 
%                                                                  rows of z 
h2=mexcdf('varget',ncid,'h2',[jmin i0-2],[jm 2]); % get 2 rows of
dy=mean(h2);     %average h2 at elev points to get h1 at v points
dy=dy';
%                                                         h1 bracketing v
mexcdf('close',ncid);
dz=-diff(sigma);
sigma=.5*(sigma(1:kbm)+sigma(2:kbm+1));
d=mean(d+e);  %average elevation and depth to get total depth at u points
d=d';
z=d*sigma;
%
% area of each grid cell face in section
area=(d.*dy)*dz;
%
% flux
if(exist('s0')),
  s=(s0-s)/s0;
end
flux=s.*u.*area;
%
y=cumsum(dy)/1000;  % convert m to km 
y=y*ones(1,kbm);
%
