function [s,jd]=baymean(cdfin,var,tind)
%function [s,jd]=baymean(cdfin,var,tind)
%
% baymean computes the average salinity in mass bay over time steps tind 
%    of a ecomsi model run.  cdfin is the ecomsi.cdf run
%
%  tind is the indices to do the tranport calculation (1 is first step)
if(nargin==2),
  tind=1:1;
end
cdfid1=mexcdf('open',cdfin,'nowrite');

[nam,nx]=mexcdf('diminq',cdfid1,'xpos');
[nam,ny]=mexcdf('diminq',cdfid1,'ypos');
[nam,nz]=mexcdf('diminq',cdfid1,'zpos');

base_date=mexcdf('attget',cdfid1,'global','base_date');
time=mexcdf('varget',cdfid1,'time',tind(1)-1,length(tind));
jd1=julian([base_date 0 0 0]);
jd=jd1+time;

h=mexcdf('varget',cdfid1,'depth',[0 0],[ny nx]);
h1=mexcdf('varget',cdfid1,'h1',[0 0],[ny nx]);
h2=mexcdf('varget',cdfid1,'h2',[0 0],[ny nx]);
sigma=mexcdf('varget',cdfid1,'sigma',[0],[nz]);
h=replace(h,-99999.,0);
%%%II are the indices of elements in Mass Bay, and are
%%% site and grid specific:  
load ii
ind=ii;
%ind=find(h>0);
hmean=mean(h(ind))
s=zeros(1,length(tind));
vol=h(ind).*h1(ind).*h2(ind);
volume=sum(sum(vol))
for i=tind,
  stot=zeros(nx,ny);
  for j=[1:nz-1],
    u=mexcdf('varget',cdfid1,var,[i-1 j-1 0 0],[1 1 ny nx]);
    stot=stot+u.*(sigma(j)-sigma(j+1));
  end
  s(i)=sum(stot(ind).* vol)/volume;
end
mexcdf('close',cdfid1);
