function [flux,totin,totout,v,y,z]=iflux(cdf,tstep,yind,x0);
%function [flux,totin,totout,u,y,z]=iflux(cdf,tstep,yind,i0);
%
% flux normal to the line x=i0
%
ny=length(yind);

sigma=mcvgt(cdf,'sigma');
nz=length(sigma)-1;

v=mcvgt(cdf,'u',[tstep-1 0 yind(1)-1 x0-1],[1 nz ny 1] );
dz=-diff(sigma);
sigma=.5*(sigma(1:nz)+sigma(2:nz+1));
d=mcvgt(cdf,'depth',[yind(1)-1 x0-2],[ny 2]); %get 2 rows of d bracketing u 
e=mcvgt(cdf,'elev',[tstep-1 yind(1)-1 x0-2],[1 ny 2]); %get 2 rows of z 
d=mean((d+e)')';  %average depth to get depth at u points
dind=find(d<-11111.);
d(dind)=0*d(dind);
z=d*sigma(1:nz);
%
% vertical flux
out_v=find(v<(-0.0));
in_v=find(v>(-0.0));
vflux=v.*(dz(:)*ones(1,ny)).*(ones(nz,1)*d');
%
%
h2=mcvgt(cdf,'h2',[yind(1)-1 x0-2],[ny 2]); % get two rows of h1 bracketing v
dy=mean(h2')';     %average h2 at elev points to get h1 at v points
%
% multiply by h2 to get total flux
%
flux=vflux.*(ones(nz,1)*dy');    % make sure we don't need to flip dx
%
out=zeros(size(flux));
in=zeros(size(flux));
iout=find(flux>0);   % flux in direction of increasing j index
% out(iout)=flux(iout);
out(out_v)=flux(out_v);
iin=find(flux<0);   % flux in direction of decreasing j index  
%in(iin)=flux(iin);
in(in_v)=flux(in_v);
totin=sum(in(:));
totout=sum(out(:));
din=sum(in);   % southward flux
dout=sum(out);
y=cumsum(dy)/1000;  % convert m to km
plot(y(:),[dout(:) -din(:) -(d(:)*100)]); %increasing flux, 
ylabel('m**3/s')
xlabel(['distance along i=' int2str(x0) '(km)'])
y=y*ones(1,nz);
v=v.';
