scrum_ctotal

 CTOTAL computes to total amount of tracer present in a given time
        step of a scrum model run.  

  Usage:  [ctot,cmean]=scrum_ctotal(cdfin,varin,[tind],[bounds])

  where:  cdfin = scrum.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]=scrum_ctotal(cdfin,varin,tind,bounds)
0002 % CTOTAL computes to total amount of tracer present in a given time
0003 %        step of a scrum model run.
0004 %
0005 %  Usage:  [ctot,cmean]=scrum_ctotal(cdfin,varin,[tind],[bounds])
0006 %
0007 %  where:  cdfin = scrum.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 
0018 if(nargin==2),
0019   tind=1;
0020 end
0021 
0022 % suppress netcdf warnings
0023 %mexcdf('setopts', 0);
0024 
0025 %
0026 ncid=mexcdf('open',cdfin,'nowrite');
0027 if(ncid==-1),
0028   disp(['file ' cdf ' not found'])
0029   return
0030 end
0031 
0032 %
0033 % Get the x, y, z dimensions
0034 [nam,nxr]=mexcdf('diminq',ncid,'xi_rho');
0035 [nam,nyr]=mexcdf('diminq',ncid,'eta_rho');
0036 [nam,nz]=mexcdf('diminq',ncid,'s_rho');
0037 
0038 if(nargin< 4),
0039   ix=0;
0040   iy=0;
0041   nx=nxr;
0042   ny=nyr;
0043   ind=1:(nx*ny);
0044 elseif (length(bounds)>4),
0045   ix=0;
0046   iy=0;
0047   nx=nxr;
0048   ny=nyr;
0049   ind=bounds;
0050 else
0051   if(bounds(2)>nxr|bounds(4)>nyr),disp('out of bounds'),return,end
0052   ix=bounds(1)-1;
0053   nx=bounds(2)-bounds(1)+1;
0054   iy=bounds(3)-1;
0055   ny=bounds(4)-bounds(3)+1;
0056   ind=1:(nx*ny);
0057 end
0058 
0059 x_rho = mexcdf ( 'varget', ncid, 'x_rho', [0 0], [-1 -1] );
0060 y_rho = mexcdf ( 'varget', ncid, 'y_rho', [0 0], [-1 -1] );
0061 
0062 [mask,status]=mexcdf('varget',ncid,'mask_rho',[iy ix],[ny nx]);
0063 if ( status == -1 )
0064     fprintf ( 2, 'Could not get ''mask_rho'' from input file.\n' );
0065     return;
0066 end
0067 
0068 [h, status] = mexcdf ( 'varget', ncid, 'h', [iy ix], [ny nx] );
0069 if ( status == -1 )
0070     fprintf ( 2, 'Could not get ''h'' from input file.\n' );
0071     return;
0072 end
0073 [sc_r, status] = mexcdf ( 'varget', ncid, 'sc_r', 0, [nz] );
0074 if ( status == -1 )
0075     fprintf ( 2, 'Could not get ''sc_r'' from input file.\n' );
0076     return;
0077 end
0078 [sc_w, status] = mexcdf ( 'varget', ncid, 'sc_w', 0, [-1] );
0079 if ( status == -1 )
0080     fprintf ( 2, 'Could not get ''sc_w'' from input file.\n' );
0081     return;
0082 end
0083 [zeta, status] = mexcdf ( 'varget', ncid, 'zeta', [tind-1 iy ix], [1 ny nx] );
0084 if ( status == -1 )
0085     fprintf ( 2, 'Could not get ''zeta'' from input file.\n' );
0086     return;
0087 end
0088 [hc, status ] = mexcdf ( 'varget1', ncid, 'hc', [0] );
0089 if ( status == -1 )
0090     fprintf ( 2, 'Could not get ''hc'' from input file.\n' );
0091     return;
0092 end
0093 [Cs_r, status] = mexcdf ( 'varget', ncid, 'Cs_r', [0], [-1] );
0094 if ( status == -1 )
0095     fprintf ( 2, 'Could not get ''Cs_r'' from input file.\n' );
0096     return;
0097 end
0098 Cs_r = Cs_r(:);
0099 [Cs_w, status] = mexcdf ( 'varget', ncid, 'Cs_w', [0], [-1] );
0100 if ( status == -1 )
0101     fprintf ( 2, 'Could not get ''Cs_w'' from input file.\n' );
0102     return;
0103 end
0104 Cs_w = Cs_w(:);
0105 
0106 s = sc_w;
0107 
0108 %
0109 % compute the depths
0110 for  k = 1:length(s) 
0111     slice = zeta*(1+s(k)) + hc*s(k) + (h - hc)*Cs_w(k);
0112     z(k,:,:) = slice;
0113 end
0114 for  k = nz:-1:1
0115     dz(k,:,:) = (z(k+1,:,:) - z(k,:,:));
0116 end
0117 
0118 pm=mexcdf('varget',ncid,'pm',[iy ix],[ny nx]);
0119 h1 = 1./pm;
0120 pn=mexcdf('varget',ncid,'pn',[iy ix],[ny nx]);
0121 h2 = 1./pn;
0122 area=h1.*h2;
0123 
0124 
0125 ctot = zeros(nx,ny);
0126 for k = 1:nz
0127     c=mexcdf('varget',ncid,varin,[tind-1 k-1 iy ix],[1 1 ny nx],1); %get a layer
0128     ctot = ctot + squeeze(dz(k,(ix+1):(ix+nx),(iy+1):(iy+ny))).*c.*mask;
0129 end
0130 ctot = ctot .* area;
0131 
0132 %
0133 % basin total
0134 ctot = sum(ctot(:));
0135 
0136 %
0137 % basin volume
0138 cvol=sum(sum((h+zeta).*area.*mask)); %divide basin total by basin volume
0139 
0140 % basin mean
0141 cmean=ctot/cvol;      %divide basin total by basin volume
0142 
0143 return;