ctotal

 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])

  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])
          bounds may also be a 1 d array containg the indices
            of the matrix 
            (just make sure it has more than 4 elements!) 
           ctot = total sum of scalar quantity c
           cmean = mean value of ctot (ctot/volume)

0001 function [ctot,cmean]=ctotal(cdfin,varin,tind,bounds)
0002 % CTOTAL computes to total amount of tracer present in a given time
0003 %        step of a ecomsi model run.
0004 %
0005 %  Usage:  [ctot,cmean]=ctotal(cdfin,varin,[tind],[bounds])
0006 %
0007 %  where:  cdfin = ecomsi.cdf run
0008 %          tind = single index of time (defaults to [1])
0009 %          bounds = [imin imax jmin jmax] limits
0010 %                   (defaults to [1 nx 1 ny])
0011 %          bounds may also be a 1 d array containg the indices
0012 %            of the matrix
0013 %            (just make sure it has more than 4 elements!)
0014 %           ctot = total sum of scalar quantity c
0015 %           cmean = mean value of ctot (ctot/volume)
0016 %
0017 if(nargin==2),
0018   tind=1;
0019 end
0020 % suppress netcdf warnings
0021 mexcdf('setopts', 0);
0022 %
0023 ncid=mexcdf('open',cdfin,'nowrite');
0024 if(ncid==-1),
0025   disp(['file ' cdf ' not found'])
0026   return
0027 end
0028 [nam,nxr]=mexcdf('diminq',ncid,'xpos');
0029 [nam,nyr]=mexcdf('diminq',ncid,'ypos');
0030 [nam,nz]=mexcdf('diminq',ncid,'zpos');
0031 if(nargin< 4),
0032   ix=0;
0033   iy=0;
0034   nx=nxr;
0035   ny=nyr;
0036   ind=1:(nx*ny);
0037 elseif (length(bounds)>4),
0038   ix=0;
0039   iy=0;
0040   nx=nxr;
0041   ny=nyr;
0042   ind=bounds;
0043 else
0044   if(bounds(2)>nxr|bounds(4)>nyr),disp('out of bounds'),return,end
0045   ix=bounds(1)-1;
0046   nx=bounds(2)-bounds(1)+1;
0047   iy=bounds(3)-1;
0048   ny=bounds(4)-bounds(3)+1;
0049   ind=1:(nx*ny);
0050 end
0051 
0052 d=mexcdf('varget',ncid,'depth',[iy ix],[ny nx]);
0053 d=replace(d,-99999,0.);
0054 %
0055 % add surface elevation to depth to get total water depth
0056 e=mexcdf('varget',ncid,'elev',[tind-1 iy ix],[1 ny nx],1);
0057 d=d+e;
0058 
0059 sigma=mexcdf('varget',ncid,'sigma',[0],[nz]);
0060 dsigma=-diff(sigma);
0061 h1=mexcdf('varget',ncid,'h1',[iy ix],[ny nx]);
0062 h2=mexcdf('varget',ncid,'h2',[iy ix],[ny nx]);
0063 area=h1.*h2;
0064 ctot=zeros(nx,ny);
0065 for j=[1:nz-1],
0066   c=mexcdf('varget',ncid,varin,[tind-1 j-1 iy ix],[1 1 ny nx],1); %get a layer
0067   ctot=ctot+dsigma(j)*c;   %multiply layer values by normalized layer thickness
0068 end
0069 ctot=(ctot.*d).*area;     % multiply depth-averaged values by depth and area
0070 ctot=sum(ctot(ind));        % find basin total
0071 cmean=ctot/sum(sum(d(ind).*area(ind))); %divide basin total by basin volume
0072 mexcdf('close',ncid);