function svol=fresh(cdfin,tind,s0)
%function svol=fresh(cdfin,tind,s0)
%
% fresh computes the volume of fresh water for the whole 
% domain time steps tind of a ecomsi model run.  
% 
%    cdfin = netCDF model output file name
%    tind = is the indices to do the tranport calculation (1 is first step)
%    s0 = the reference salinity 

if(nargin==1),
  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');

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);
load ii
ind=ii;
%ind=find(h>0);
hmean=mean(h(ind))
svol=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,'salt',[i-1 j-1 0 0],[1 1 ny nx]);
    stot=stot+u.*(sigma(j)-sigma(j+1));
  end
  svol(i)=sum(sum( ((s0-stot(ind))/s0) .* vol));
end
mexcdf('close',cdfid1);
