% calculate the flushing time of Mass Bay based on the mass
% flux normal to the line j=33;
%
tstep=3;  % which time step to do flux
yind=33;
xind=18:50;  % xindex 
nz=11;    %number of vertical layers

ny=length(yind);
nx=length(xind);
vflux=xzslice('ecomsi.cdf','vflux',tstep,yind);  % units m^2/s
vflux=vflux(1:nz,xind);
%
% if flux is due to tidal rectification, might want to subtract tidal flux
h1=mcvgt('ecomsi.cdf','h1',[yind-1 xind(1)-1],[2 nx]); % get h2 at j=33,34
d=mcvgt('ecomsi.cdf','depth',[yind-1 xind(1)-1],[2 nx]); % get h2 at j=33,34
h1=mean(h1);     %average h2 at elev points to get h2 at u points
d=mean(d);     %average h2 at elev points to get h2 at u points
% multiply by h1 to get total flux
flux=vflux.*(ones(nz,1)*h1);    % make sure we don't need to flip h2
out=zeros(flux);
in=zeros(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(:));
x=cumsum(h1)/1000;
din=sum(in);
dout=sum(out);
plot(x,[dout' din' (dout+din)' -(d*100)'],x,dout','x',...
   x,din','+',x,(dout+din)','o',x,-(d*100)','*')
hold on
ax=axis;
ax=axis;
plot([ax(1) ax(2)],[0 0])
grid
ylabel('m**3/s')
xlabel('distance along j=33 (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
