function [totin,totout,hflux,hflux_tot,v,x,z]=xzflush(cdf,tstep,yind,x0);
%function [inflush,outflush,hflux,hflux_tot,v,x,z]=xzflush(cdf,tstep,yind,x0);
% calculate the flushing time of Mass Bay based on the mass
% flux normal to the line i=i0
%
if(nargin==3),
  x0=yind;
  yind=1:68;
end
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] );
t1=mcvgt(cdf,'temp',[tstep-1 0 yind(1)-1 x0-1],[1 nz ny 1]);
t2=mcvgt(cdf,'temp',[tstep-1 0 yind(1)-1 x0-2],[1 nz ny 1]);
t=0.5*(t1+t2);
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 v 
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
%
% heat flux
%
hflux=flux.*t;
hflux_tot=sum(hflux(:));
%
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)'])
volume=1.2773e11;
disp(' flushing time in days')
inflush=volume/totin/3600/24  % flushing time in days
outflush=volume/totout/3600/24  % flushing time in days
y=y*ones(1,nz);
