freshtot

 FRESHTOT computes to total amount of fresh water present in a given time
        step of a ecomsi model run.  

  Usage: [ftot,fmean,fthick,x,y]=freshtot(cdfin,tind,bounds,s0)

  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]) if -1 is specified
          bounds may also be a 1 d array containg the indices
            of the matrix 
            (just make sure it has more than 4 elements!) 
           ftot = total sum of scalar quantity c
           fmean = mean value of ctot (ctot/volume)
           fthick = depth-integrated fresh water  (freshwater thickness)

0001 function [ctot,cmean,fthick,x,y]=freshtot(cdfin,tind,bounds,s0)
0002 % FRESHTOT computes to total amount of fresh water present in a given time
0003 %        step of a ecomsi model run.
0004 %
0005 %  Usage: [ftot,fmean,fthick,x,y]=freshtot(cdfin,tind,bounds,s0)
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]) if -1 is specified
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 %           ftot = total sum of scalar quantity c
0015 %           fmean = mean value of ctot (ctot/volume)
0016 %           fthick = depth-integrated fresh water  (freshwater thickness)
0017 %
0018 if(nargin==2),
0019   tind=1;
0020 end
0021 % suppress netcdf warnings
0022 mexcdf('setopts', 0);
0023 %
0024 ncid=mexcdf('open',cdfin,'nowrite');
0025 if(ncid==-1),
0026   disp(['file ' cdf ' not found'])
0027   return
0028 end
0029 [nam,nxr]=mexcdf('diminq',ncid,'xpos');
0030 [nam,nyr]=mexcdf('diminq',ncid,'ypos');
0031 [nam,nz]=mexcdf('diminq',ncid,'zpos');
0032 if(nargin< 4 | bounds==-1),
0033   ix=0;
0034   iy=0;
0035   nx=nxr;
0036   ny=nyr;
0037   ind=1:(nx*ny);
0038 elseif (length(bounds)>4),
0039   ix=0;
0040   iy=0;
0041   nx=nxr;
0042   ny=nyr;
0043   ind=bounds;
0044 else
0045   if(bounds(2)>nxr|bounds(4)>nyr),disp('out of bounds'),return,end
0046   ix=bounds(1)-1;
0047   nx=bounds(2)-bounds(1)+1;
0048   iy=bounds(3)-1;
0049   ny=bounds(4)-bounds(3)+1;
0050   ind=1:(nx*ny);
0051 end
0052 
0053 d=mexcdf('varget',ncid,'depth',[iy ix],[ny nx]);
0054 x=mexcdf('varget',ncid,'x',[iy ix],[ny nx]);
0055 y=mexcdf('varget',ncid,'y',[iy ix],[ny nx]);
0056 iland=find(d==-99999.);
0057 d(iland)=zeros(size(iland));
0058 %
0059 % add surface elevation to depth to get total water depth
0060 e=mexcdf('varget',ncid,'elev',[tind-1 iy ix],[1 ny nx],1);
0061 d=d+e;
0062 
0063 sigma=mexcdf('varget',ncid,'sigma',[0],[nz]);
0064 dsigma=-diff(sigma);
0065 h1=mexcdf('varget',ncid,'h1',[iy ix],[ny nx]);
0066 h2=mexcdf('varget',ncid,'h2',[iy ix],[ny nx]);
0067 area=h1.*h2;
0068 ctot=zeros(nx,ny);
0069 for j=[1:nz-1],
0070   c=mexcdf('varget',ncid,'salt',[tind-1 j-1 iy ix],[1 1 ny nx],1); %get a layer
0071   c=(s0-min(c,s0))/s0;
0072 %  ibad=find(c>s0);   % don't count salinities greater than the reference salinity!
0073 %  c(ibad)=0*c(ibad);
0074   ctot=ctot+dsigma(j)*c;   %multiply layer values by normalized layer thickness
0075 end
0076 mexcdf('close',ncid);
0077 fthick=(ctot.*d);     % integrated fresh water
0078 ctot=fthick.*area;
0079 ctot=sum(ctot(ind));        % find basin total
0080 cmean=ctot/sum(sum(d(ind).*area(ind))); %divide basin total by basin volume
0081 fthick=max(fthick,0.);
0082 fthick(iland)=fthick(iland)*nan;