function [s1,s2,jd]=baymean(cdfin,var,tind)
%function [s1,s2,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))
vol=h(ind).*h1(ind).*h2(ind);
volume=sum(sum(vol))
sigma(nz)=-1;
for i=tind,
  stot1=0;
  stot2=0;
  for j=[1:nz-1],
    z=.5*(sigma(j)+sigma(j+1))*h;
    u=mexcdf('varget',cdfid1,var,[i-1 j-1 0 0],[1 1 ny nx]);
    iup=find(z(ind)>=-15);
    ilw=find(z(ind)<-15);
if(length(iup)>0),
    s1tot=s1tot+u(ind(iup)).*(sigma(j)-sigma(j+1)).*(h(ind(iup))).*...
       h1(ind(iup)).*h2(ind(iup));
    end
    if(length(ilw)>0),
     s2tot=s2tot+u(ind(ilw)).*(sigma(j)-sigma(j+1)).*(h(ind(ilw))).*...
       h1(ind(ilw)).*h2(ind(ilw));
    end
  end
  s1(i)=s1tot;
  s2(i)=s2tot;
end
mexcdf('close',cdfid1);
