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,x,z]=jsflux(cdf,var,tstep,j0,irange);
%     irange=extent of I along section j=j0, eg. [5 45]; 
%     

ncid=mexcdf('open',cdf);
[dim_name,x_size,rcode]=mexcdf('diminq',ncid,'xpos');
[dim_name,k_size,rcode]=mexcdf('diminq',ncid,'zpos');
kbm=k_size-1;
if(~exist('irange')|(isnan(irange))),
  imin=0;
  im=x_size;
else
  imin=irange(1)-1;
  im=min(x_size,irange(2)-irange(1)+1);
end
sigma=mexcdf('varget',ncid,'sigma',0,kbm+1);
sigma=sigma';
d=mexcdf('varget',ncid,'depth',[j0-2 imin],[2 im]); %get 2 rows of d 
iland=find(d==-99999);
d1=d(:,1);
d2=d(:,2);
iland1=find(d1==-99999);
iland2=find(d2==-99999);
d(iland)==zeros(size(iland));
v=mexcdf('varget',ncid,'v',[tstep-1 0 j0-1 imin],[1 kbm 1 im],1);
s=mexcdf('varget',ncid,var,[tstep-1 0 j0-1 imin],[1 kbm 2 im],1);
s1=s(:,1:2:2*kbm);
s1(iland1,:)=s1(iland1,:)*nan;
s2=s(:,2:2:2*kbm);
s2(iland2,:)=s2(iland2,:)*nan;
s=0.5*(s1+s2);
%                                                         bracketing v 
e=mexcdf('varget',ncid,'elev',[tstep-1 j0-2 imin],[1 2 im],1); %get 2 
                                                               %rows of z 
h1=mexcdf('varget',ncid,'h1',[j0-2 imin],[2 im]); % get 2 rows of
dx=mean(h1')';     %average h1 at elev points to get h1 at v points
%                                                         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 v points
z=d*sigma;
%
% area of each grid cell face in section
area=(d.*dx)*dz;
%
% flux
if(exist('s0')),
  s=min(0,(s0-s)/s0);
end
flux=s.*v.*area;
%
x=cumsum(dx)/1000;  % convert m to km 
x=x*ones(1,kbm);
%
