Description of scrum

0001 function [wmean,x,y]=scrum_depaveuv(cdfin,tind,bounds)
0002 % ECOM_DEPAVEUV computes the depth-averaged value of velocity at
0003 %        a given time step of a SCRUM model run.
0004 %
0005 %  Usage:  [wmean,x,y]=scrum_depaveuv(cdf,[tstep],[bounds])
0006 %
0007 %  where:  cdf = ecomsi.cdf run
0008 %          tstep = time step (default = 1)
0009 %          bounds = [imin imax jmin jmax] limits
0010 %                   (default = [1 nx 1 ny])
0011 %
0012 %          wmean = depth-averaged velocity
0013 %          x = x locations of the returned array wmean
0014 %          y = y locations of the returned array wmean
0015 %
0016 %  Example 1:  [wmean,x,y]=scrum_depaveuv('scrum.cdf');
0017 %
0018 %       computes the depth-averaged velocity at the 1st time step
0019 %       over the entire domain.
0020 %
0021 %  Example 2:  [wmean,x,y]=scrum_depaveuv('scrum.cdf',10);
0022 %
0023 %       computes the depth-averaged velocity at the 10th time step
0024 %       over the entire domain.
0025 %
0026 %  Example 3:  [wmean,x,y]=scrum_depaveuv('scrum.cdf',10,[10 30 30 50]);
0027 %
0028 %       computes the depth-averaged velocity at the 10th time step
0029 %       in the subdomain defined by i=10:30 and j=30:50.
0030 %
0031 %
0032 
0033 if(nargin==1),
0034   tind=1;
0035 end
0036 
0037 ncid=mexcdf('open',cdfin,'nowrite');
0038 if(ncid==-1),
0039   disp(['file ' cdf ' not found'])
0040   return
0041 end
0042 
0043 
0044 %
0045 % Acquire the grid.
0046 %
0047 % If "lon_rho" and "lat_rho" are present, grab them.
0048 % Otherwise, get "x_rho" and "y_rho".
0049 [lon_rho_varid, rcode] = ncmex('VARID', ncid, 'lon_rho');
0050 [lat_rho_varid, rcode] = ncmex('VARID', ncid, 'lat_rho');
0051 if ( (lon_rho_varid >= 0) | (lat_rho_varid >= 0) )
0052     x=ncmex('varget',ncid,'lon_rho',[0 0],[-1 -1]);
0053     y=ncmex('varget',ncid,'lat_rho',[0 0],[-1 -1]);
0054 else
0055     x=ncmex('varget',ncid,'x_rho',[0 0],[-1 -1]);
0056     y=ncmex('varget',ncid,'y_rho',[0 0],[-1 -1]);
0057 end
0058 
0059 x_rho = x;
0060 y_rho = y;
0061 
0062 [eta_rho_length, xi_rho_length] = size(x_rho);
0063 eta_u_length = eta_rho_length;
0064 eta_v_length = eta_rho_length-1;
0065 xi_u_length = xi_rho_length-1;
0066 xi_v_length = xi_rho_length;
0067 eta_psi_length = eta_rho_length-1;
0068 xi_psi_length = xi_rho_length-1;
0069 
0070 
0071 %
0072 % construct the grids at the psi points.
0073 xtemp = (x_rho(:,1:(xi_rho_length-1)) + x_rho(:,2:xi_rho_length))/2;
0074 x = (xtemp(1:(eta_rho_length-1),:)  + xtemp(2:eta_rho_length,:))/2;
0075 ytemp = (y_rho(:,1:(xi_rho_length-1)) + y_rho(:,2:xi_rho_length))/2;
0076 y = (ytemp(1:(eta_rho_length-1),:)  + ytemp(2:eta_rho_length,:))/2;
0077 
0078 
0079 
0080 %
0081 % Get ubar and vbar.
0082 [u_bar, status] = ncmex ( 'varget', ncid, 'ubar', [tind-1 0 0], [1 -1 -1] );
0083 if ( status == -1 )
0084     fprintf ( 'scrum_depaveuv:  could not get ubar in %s.', cdf );
0085     return;
0086 end
0087 [v_bar, status] = ncmex ( 'varget', ncid, 'vbar', [tind-1 0 0], [1 -1 -1] );
0088 if ( status == -1 )
0089     fprintf ( 'scrum_depaveuv:  could not get vbar in %s.', cdf );
0090     return;
0091 end
0092 
0093 %
0094 % If there is a mask_rho, we will want to use it later.
0095 [mask_rho_varid, status] = ncmex ( 'varid', ncid, 'mask_rho' );
0096 if ( status ~= -1 )
0097     [mask_rho, status] = ncmex ( 'varget', ncid, 'mask_rho', [0 0], [-1 -1] );
0098     if ( status == -1 )
0099         fprintf ( 'scrum_zsliceuv:  could not get ''mask_rho'' in %s.', cdf );
0100         return;
0101     end
0102 end
0103 
0104 
0105 ncmex ( 'close', ncid );
0106 
0107 %
0108 % Average the velocities onto the psi points.
0109 [r,c] = size(u_bar);
0110 u_bar = (u_bar(:,1:c-1) + u_bar(:,2:c))/2;
0111 [r,c] = size(v_bar);
0112 v_bar = (v_bar(1:r-1,:) + v_bar(2:r,:))/2;
0113 
0114 
0115 %
0116 % Use complex numbers.
0117 wmean = u_bar + sqrt(-1) * v_bar;
0118 
0119 
0120 %
0121 % Cut out the requested part.
0122 [nx,ny] = size(wmean);
0123 if ( nargin < 3 )
0124     x1 = 1;
0125     y1 = 1;
0126     x2 = nx;
0127     y2 = ny;
0128 else
0129     if (min(bounds(:))<=0)
0130     disp('out of bounds');
0131     return;
0132     end
0133     if (length(bounds)~=4)
0134     disp('out of bounds');
0135     return;
0136     end
0137     if (bounds(2) >= nx )
0138     disp('out of bounds');
0139     return;
0140     end
0141     if (bounds(4) >= ny )
0142     disp('out of bounds');
0143     return;
0144     end
0145 
0146     x1 = bounds(1);
0147     x2 = bounds(2);
0148     y1 = bounds(3);
0149     y2 = bounds(4);
0150 
0151 end
0152 
0153 %
0154 % apply the rotation
0155 %
0156 % Since the "angle" variable is not always present, we
0157 % need to construct it.  In the scrum files, it is apparently
0158 % always in degrees.  Don't bother with degrees here.
0159 angle = zeros(size(x));
0160 
0161 
0162 
0163 [r,c] = size(x);
0164 j = [2:c-1];
0165 for i = 2:r-1
0166     angle(i,j) = atan2(y(i+1,j)-y(i-1,j), x(i+1,j)-x(i-1,j));
0167 end
0168 
0169 % rotate into east and north components
0170   wmean=wmean.*exp(sqrt(-1)*angle);
0171 
0172 
0173 wmean = wmean(x1:x2,y1:y2);
0174 x = x(x1:x2,y1:y2);
0175 y = y(x1:x2,y1:y2);
0176 
0177 %
0178 % Mask out the land.
0179 % Since we average the u and v across cells, the appropriate mask
0180 % dimension is that of "mask_psi".  I don't want to take the chance
0181 % that it may not be present in the file, so I compute it from
0182 % mask_rho.
0183 [umask,vmask,pmask]=uvp_masks ( mask_rho );
0184 mask_inds = find ( pmask == 0 );
0185 wmean(mask_inds) = NaN * ones(size(mask_inds));
0186 
0187 
0188 return;
scrum_depaveuv

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
