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

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

v=mcvgt(cdf,'v',[tstep-1 0 yind(1)-1 xind(1)-1],[1 nz 1 nx]);
sigma=-diff(sigma);
sigma(nz)=0;
sigma=(sigma'*ones(1,length(xind)));
d=mcvgt(cdf,'depth',[yind(1)-1 xind-1],[2 nx]); % 
d=mean(d);     %average depth at elev points to get depth at u points
dind=find(d<-11111.);
d(dind)=0*d(dind);
d2=ones(nz,1)*d;
vflux=v.*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
%vflux=vflux-v2flux;
%
h1=mcvgt(cdf,'h1',[yind(1)-1 xind-1],[2 nx]); % get h1
h1=mean(h1);     %average h1 at elev points to get h1 at v points
% multiply by h1 to get total flux
flux=vflux.*(ones(nz,1)*h1);    % make sure we don't need to flip h1
out=zeros(size(flux));
in=zeros(size(flux));
iout=find(flux>0);   % flux in direction of increasing j index
out(iout)=flux(iout);
iin=find(flux<0);   % flux in direction of decreasing j index  
in(iin)=flux(iin);
totin=sum(in(:));
totout=sum(out(:));
din=sum(in);   % southward flux
dout=sum(out);
x=cumsum(h1)/1000;
plot(x,[dout' -din' -(d*100)']);   % increasing j flux, - decreasing j, depth
hold on
ylabel('m**3/s')
xlabel(['distance along j=' int2str(yind) '(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
