% calculate the flushing time of Mass Bay based on the mass
% flux normal to the line i=48;
%
tstep=5;  % which time step to do flux
xind=48;  % xindex 
yind=18:49;  
nz=12;    %number of vertical layers

cdf='ecomsi.cdf';
ny=length(yind);
u=yzslice(cdf,'u',tstep,xind);  % units m^2/s
u=u(yind,1:12)';
sigma=mcvgt(cdf,'sigma');
sigma=-diff(sigma);
sigma(nz)=0;
sigma=(sigma'*ones(1,32));
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
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(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(:));
din=sum(in);
dout=sum(out);
y=cumsum(h2)/1000;
plot(y,[dout' din' (dout+din)' -(d*100)'],y,dout','x',...
   y,din','+',y,(dout+din)','o',y,-(d*100)','*')
hold on
ax=axis;
ax=axis;
plot([ax(1) ax(2)],[0 0])
ylabel('m**3/s')
xlabel('distance along i=48 (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
