function [totin,totout,hflux,hflux_tot,v,x,z]=yzflush(cdf,tstep,xind,y0);
%function [inflush,outflush,hflux,hflux_tot,v,x,z]=yzflush(cdf,tstep,xind,j0);
% calculate the flushing time of Mass Bay based on the mass
% flux normal to the line j=j0
%
if(nargin==3),
  y0=xind;
  xind=1:68;
end
nx=length(xind);

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

v=mcvgt(cdf,'v',[tstep-1 0 y0-1 xind(1)-1],[1 nz 1 nx]);
t1=mcvgt(cdf,'temp',[tstep-1 0 y0-1 xind(1)-1],[1 nz 1 nx]);
t2=mcvgt(cdf,'temp',[tstep-1 0 y0-2 xind(1)-1],[1 nz 1 nx]);
t=0.5*(t1+t2);
dz=-diff(sigma);
sigma=.5*(sigma(1:nz)+sigma(2:nz+1));
d=mcvgt(cdf,'depth',[y0-2 xind(1)-1],[2 nx]); %get 2 rows of d bracketing v 
e=mcvgt(cdf,'elev',[tstep-1 y0-2 xind(1)-1],[1 2 nx]); %get 2 rows of z 
d=mean(d+e);  %average depth (defined at elev points) to get depth at v points
dind=find(d<-11111.);
d(dind)=0*d(dind);
z=(sigma(1:nz)')*d;
%
% vertical flux
out_v=find(v<(-0.0));
in_v=find(v>(-0.0));
vflux=v.*(dz(:)*ones(1,nx)).*(ones(nz,1)*d);
%
%
h1=mcvgt(cdf,'h1',[y0-2 xind(1)-1],[2 nx]); % get two rows of h1 bracketing v
dx=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)*dx);    % 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);
x=cumsum(dx)/1000;  % convert m to km
plot(x,[dout(:) -din(:) -(d(:)*100)]); %increasing u flux, - decreasing j, depth
ylabel('m**3/s')
xlabel(['distance along j=' int2str(y0) '(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
x=ones(nz,1)*x;
