function [cmean,x,y]=depthave(cdfin,varin,tind,bounds)
% DEPTHAVE computes to depth-averaged variable in a given time
%        step of a ecomsi model run.  
%
%  Usage:  [cmean]=depthave(cdf,var,tstep)
%
%  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)
%  
if(nargin==2),
  tind=1;
end
% suppress netcdf warnings
mexcdf('setopts', 0);
%
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
depth=mexcdf('varget',cdfid1,'depth',[ix iy],[ny nx]);
x=mexcdf('varget',cdfid1,'x',[ix iy],[ny nx]);
y=mexcdf('varget',cdfid1,'y',[ix iy],[ny nx]);
sigma=mexcdf('varget',cdfid1,'sigma',[0],[nz]);
dsigma=-diff(sigma);
h1=mexcdf('varget',cdfid1,'h1',[ix iy],[ny nx]);
h2=mexcdf('varget',cdfid1,'h2',[ix iy],[ny nx]);
area=h1.*h2;
cmean=zeros(nx,ny);
for j=[1:nz-1],
  c=mexcdf('varget',cdfid1,varin,[tind-1 j-1 iy ix],[1 1 ny nx],1); %get a layer
  cmean=cmean+dsigma(j)*c;   %multiply layer values by normalized layer thickness
end
mexcdf('close',cdfid1);
dind=find(depth==-99999 | depth==0);
cmean(dind)=cmean(dind)*NaN;

