function [ctot,cmean]=ctotal2(cdfin,varin,tind,bounds,depths)
% CTOTAL computes to total amount of tracer present in a given time
%        step of a ecomsi model run.  
%
%  Usage:  [ctot,cmean]=ctotal(cdfin,varin,[tind],[bounds],[depths])
%
%  where:  cdfin = ecomsi.cdf run
%          tind = single index of time (defaults to [1])
%          bounds = [imin imax jmin jmax] limits  
%                   (defaults to [1 nx 1 ny])
%          
%           ctot = total sum of scalar quantity c
%           cmean = mean value of ctot (ctot/volume)
%  
%  tind is the indices of time steps to average over (1 is first step)
if(nargin==2),
  tind=1;
end
depths=abs(depths);
cdfid1=mexcdf('open',cdfin,'nowrite');
[nam,nxr]=mexcdf('diminq',cdfid1,'xpos');
[nam,nyr]=mexcdf('diminq',cdfid1,'ypos');
[nam,nz]=mexcdf('diminq',cdfid1,'zpos');
if(nargin< 4),
  ix=0;
  iy=0;
  nx=nxr;
  ny=nyr;
else
  if(bounds(2)>nxr|bounds(4)>nyr),disp('out of bounds'),return,end
  ix=bounds(1)-1;
  nx=bounds(2)-bounds(1)+1;
  iy=bounds(3)-1;
  ny=bounds(4)-bounds(3)+1;
end
d=mexcdf('varget',cdfid1,'depth',[iy ix],[ny nx]);
sigma=mexcdf('varget',cdfid1,'sigma',[0],[nz]);
sigma2=.5*(sigma(1:nz-1)+sigma(2:nz));
d=replace(d,-99999,0.);
h1=mexcdf('varget',cdfid1,'h1',[iy ix],[ny nx]);
h2=mexcdf('varget',cdfid1,'h2',[iy ix],[ny nx]);
area=h1.*h2;
ctot=0.;
for j=bounds(3):bounds(4);
for i=bounds(1):bounds(2);
  if(d(i,j)~=0.0),
  c=mexcdf('varget',cdfid1,varin,[tind-1 0 j-1 i-1],[1 nz-1 1 1]);%get a profile
  z=sigma*d(i,j);
  dmax=max(depths(2),d(i,j)))/50;
  zi=-[abs(depths(1)):dzi:abs(depths(2))];
  ci=spline(z,c,zi);
  ctot=ctot+sum(ci)*dzi*area(i,j);
  end
end
cmean=ctot/sum(sum(d.*area)); %divide basin total by basin volume
mexcdf('close',cdfid1);
