function [inflush,outflush]=xzflush(cdf,tstep,xind,yind);
%function [inflush,outflush]=xzflush(cdf,tstep,xind,yind);
% calculate the flushing time of Mass Bay based on the mass
% flux normal to the line i=xind
%
if(nargin==3),
  yind=1:68;
end
ny=length(yind);

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

u=mcvgt(cdf,'u',[tstep-1 0 yind(1)-1 xind(1)-1],[1 nz ny 1]);
sigma=-diff(sigma);
sigma(nz)=0;
sigma=(sigma'*ones(1,length(yind)));
d=mcvgt(cdf,'depth',[yind(1)-1 xind-1],[ny 2]); % get h2 at i=48,49
d=mean(d');     %average h2 at elev points to get h2 at u points
dind=find(d<-11111.);
d(dind)=0*d(dind);
d2=ones(nz,1)*d;
uflux=u.*sigma.*d2;
%
% if flux is due to tidal rectification, might want to subtract tidal flux
%load m2flux   % flux from no wind run with xind=48, yind=18:49, nz=11
%uflux=uflux-m2flux;

h2=mcvgt(cdf,'h2',[yind(1)-1 xind-1],[ny 2]); % get h2 at i=48,49
h2=mean(h2');     %average h2 at elev points to get h2 at u points
%
% multiply by h2 to get total flux
%
flux=uflux.*(ones(nz,1)*h2);    % make sure we don't need to flip h2
out=zeros(size(flux));
in=zeros(size(flux));
iout=find(flux>0);   % flux out of bay
out(iout)=flux(iout);
iin=find(flux<0);    % flux in to bay
in(iin)=flux(iin);
totin=sum(in(:));
totout=sum(out(:));
din=sum(in);
dout=sum(out);
y=cumsum(h2)/1000;
plot(y,[dout' -din' -(d*100)']);
hold on
ylabel('m**3/s')
xlabel(['distance along i=' int2str(xind) '(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
hold off
