scrum_zsliceuv

  SCRUM_ZSLICE  Returns a matrix containing a horizontal slice of velocity
          at a specified depth at a given time step from a SCRUM NetCDF 
          file.  Regions of grid that are shallower than requested value 
          are returned as NaNs.

       USAGE: [w,x,y]=scrum_zsliceuv(cdf,var,zdepth)

   where zuser is a depth in meters (e.g -10.)

0001 function [w,x,y]=scrum_zsliceuv(cdf,timestep,zdepth)
0002 %  SCRUM_ZSLICE  Returns a matrix containing a horizontal slice of velocity
0003 %          at a specified depth at a given time step from a SCRUM NetCDF
0004 %          file.  Regions of grid that are shallower than requested value
0005 %          are returned as NaNs.
0006 %
0007 %       USAGE: [w,x,y]=scrum_zsliceuv(cdf,var,zdepth)
0008 %
0009 %   where zuser is a depth in meters (e.g -10.)
0010 %
0011 
0012 % John Evans (joevans@usgs.gov)
0013 
0014 if (nargin<2 | nargin>4),
0015   help scrum_zslice; return
0016 end
0017 
0018 ncmex('setopts',0);
0019 ncid=ncmex('open',cdf,'nowrite');
0020 if(ncid==-1),
0021   disp(['file ' cdf ' not found'])
0022   return
0023 end
0024 
0025 %
0026 % Acquire the grid.
0027 %
0028 % If "lon_rho" and "lat_rho" are present, grab them.
0029 % Otherwise, get "x_rho" and "y_rho".
0030 [lon_rho_varid, rcode] = ncmex('VARID', ncid, 'lon_rho');
0031 [lat_rho_varid, rcode] = ncmex('VARID', ncid, 'lat_rho');
0032 if ( (lon_rho_varid >= 0) | (lat_rho_varid >= 0) )
0033     x=ncmex('varget',ncid,'lon_rho',[0 0],[-1 -1]);
0034     y=ncmex('varget',ncid,'lat_rho',[0 0],[-1 -1]);
0035 else
0036     x=ncmex('varget',ncid,'x_rho',[0 0],[-1 -1]);
0037     y=ncmex('varget',ncid,'y_rho',[0 0],[-1 -1]);
0038 end
0039 
0040 x_rho = x';
0041 y_rho = y';
0042 
0043 
0044 [eta_rho_length, xi_rho_length] = size(x_rho);
0045 eta_u_length = eta_rho_length;
0046 eta_v_length = eta_rho_length-1;
0047 xi_u_length = xi_rho_length-1;
0048 xi_v_length = xi_rho_length;
0049 eta_psi_length = eta_rho_length-1;
0050 xi_psi_length = xi_rho_length-1;
0051 
0052 
0053 %
0054 % construct the grids at the psi points.
0055 xtemp = (x_rho(:,1:(xi_rho_length-1)) + x_rho(:,2:xi_rho_length))/2;
0056 x = (xtemp(1:(eta_rho_length-1),:)  + xtemp(2:eta_rho_length,:))/2;
0057 ytemp = (y_rho(:,1:(xi_rho_length-1)) + y_rho(:,2:xi_rho_length))/2;
0058 y = (ytemp(1:(eta_rho_length-1),:)  + ytemp(2:eta_rho_length,:))/2;
0059 
0060 
0061 
0062 
0063 %
0064 % Get all the variables needed to compute z.
0065 % The equation is
0066 %
0067 % z = zeta * (1 + s) + hc*s + (h - hc)*C(s)
0068 
0069 [s_rho_dimid, status] = ncmex ( 'dimid', ncid, 's_rho' );
0070 if ( status == -1 )
0071     fprintf ( 2, 'Could not get s_rho dimid from %s.\n', cdf );
0072     ncmex ( 'close', ncid );
0073     return;
0074 end
0075 [dimname, s_rho_length, status] = ncmex ( 'diminq', ncid, s_rho_dimid );
0076 if ( status == -1 )
0077     fprintf ( 2, 'Could not get s_rho length from %s.\n', cdf );
0078     ncmex ( 'close', ncid );
0079     return;
0080 end
0081 
0082 
0083 [sc_r, status] = ncmex ( 'varget', ncid, 'sc_r', [0], [-1] );
0084 if ( status == -1 )
0085     fprintf ( 2, 'Could not get ''sc_r'' variable from %s.\n', cdf );
0086     ncmex ( 'close', ncid );
0087     return;
0088 end
0089 
0090 [zeta, status] = ncmex ( 'varget', ncid, 'zeta', [timestep 0 0], [1 -1 -1] );
0091 zeta = zeta';
0092 
0093 [hc_varid, status] = ncmex ( 'varid', ncid, 'hc' );
0094 if ( status == -1 )
0095     fprintf ( 2, 'Could not get hc variable from %s.\n', cdf );
0096     ncmex ( 'close', ncid );
0097     return;
0098 end
0099 [hc, status] = ncmex ( 'varget1', ncid, hc_varid, [0] );
0100 if ( status == -1 )
0101     fprintf ( 2, 'Could not get hc variable from %s.\n', cdf );
0102     ncmex ( 'close', ncid );
0103     return;
0104 end
0105 
0106 [h, status] = ncmex ( 'varget', ncid, 'h', [0 0], [-1 -1] );
0107 if ( status == -1 )
0108     fprintf ( 'scrum_zsliceuv:  could not get ''h'' in %s.', cdf );
0109     return;
0110 end
0111 h = h';
0112 
0113 [Cs_r, status] = ncmex ( 'varget', ncid, 'Cs_r', [0], [-1] );
0114 if ( status == -1 )
0115     fprintf ( 'scrum_zsliceuv:  could not get ''Cs_r'' in %s.', cdf );
0116     return;
0117 end
0118 
0119 
0120 
0121 %
0122 % Construct the depth.
0123 n = length(sc_r);
0124 z = [];
0125 
0126 %
0127 % if it looks upside-down, replace with following line
0128 %for i = n:-1:1
0129 for i = 1:n
0130     zi = zeta * (1+sc_r(i)) + hc*sc_r(i) + (h - hc)*Cs_r(i);
0131     z = cat ( 3, z, zi );
0132 end
0133 z = permute ( z, [3 1 2] );
0134     
0135 
0136 %
0137 % Now average z, h, and zeta across the grid cell centers.
0138 z = (z(:,1:eta_rho_length-1,:)+z(:,2:eta_rho_length,:))/2;
0139 z = (z(:,:,1:xi_rho_length-1)+z(:,:,2:xi_rho_length))/2;
0140 h = (h(1:eta_rho_length-1,:)+h(2:eta_rho_length,:))/2;
0141 h = (h(:,1:xi_rho_length-1)+h(:,2:xi_rho_length))/2;
0142 zeta = (zeta(1:eta_rho_length-1,:)+zeta(2:eta_rho_length,:))/2;
0143 zeta = (zeta(:,1:xi_rho_length-1)+zeta(:,2:xi_rho_length))/2;
0144 
0145 
0146 %
0147 % reshape the depth such that each vertical profile is a column.
0148 z = reshape ( z, [s_rho_length eta_psi_length*xi_psi_length] );
0149 
0150 
0151 
0152 %
0153 % Retrieve u at the given timestep.
0154 [data_u, status] = ncmex ( 'varget', ncid, 'u', [timestep 0 0 0], [1 -1 -1 -1] );
0155 if ( status == -1 )
0156     fprintf ( 'scrum_zsliceuv:  could not get ''u'' in %s.', cdf );
0157     return;
0158 end
0159 data_u = permute ( data_u, [3 2 1] );
0160 data_u = (data_u(:,1:eta_u_length-1,:) + data_u(:,2:eta_u_length,:))/2;
0161 data_u = reshape(data_u,[s_rho_length eta_psi_length*xi_psi_length]);
0162 
0163 %
0164 % Retrieve v at the given timestep.
0165 if ( status == -1 )
0166     fprintf ( 'scrum_zsliceuv:  could not get ''u'' in %s.', cdf );
0167     return;
0168 end
0169 [data_v, status] = ncmex ( 'varget', ncid, 'v', [timestep 0 0 0], [1 -1 -1 -1] );
0170 data_v = permute ( data_v, [3 2 1] );
0171 data_v = (data_v(:,:,1:xi_v_length-1) + data_v(:,:,2:xi_v_length))/2;
0172 data_v = reshape(data_v,[s_rho_length eta_psi_length*xi_psi_length]);
0173 
0174 
0175 %
0176 % If there is a mask_rho, we will want to use it later.
0177 [mask_rho_varid, status] = ncmex ( 'varid', ncid, 'mask_rho' );
0178 if ( status ~= -1 )
0179     [mask_rho, status] = ncmex ( 'varget', ncid, 'mask_rho', [0 0], [-1 -1] );
0180     if ( status == -1 )
0181         fprintf ( 'scrum_zsliceuv:  could not get ''mask_rho'' in %s.', cdf );
0182         return;
0183     end
0184 end
0185 
0186 ncmex ( 'close', ncid );
0187 
0188 
0189 
0190 
0191 
0192 velocity = data_u + sqrt(-1)*data_v;
0193 
0194 
0195 %
0196 % define z and velocity at zeta (elevation)
0197 z = [z; reshape(zeta, [1 eta_psi_length*xi_psi_length] ) ];
0198 velocity = [velocity; velocity(s_rho_length,:)];
0199 
0200 %
0201 % define z at values greater than zeta to be Inf
0202 % define velocity to be NaN
0203 z = [z; Inf*ones(1,eta_psi_length*xi_psi_length)];
0204 velocity = [velocity; NaN * ones(1,eta_psi_length*xi_psi_length)];
0205 
0206 %
0207 % define z at bathymetry values to be those bathymetry values
0208 % define velocity at bathymetry values to be same as data at lowest
0209 % existing values.
0210 z = [ -1*reshape(h,[1 eta_psi_length*xi_psi_length]); z];
0211 velocity = [ velocity(1,:); velocity ];
0212 
0213 
0214 %
0215 % define z at depth greater than bathymetry to be -Inf
0216 % define velocity at depth greater than bathymetry to be NaN
0217 z = [ -Inf*ones(1,eta_psi_length*xi_psi_length); z ];
0218 velocity = [ NaN * ones(1,eta_psi_length*xi_psi_length); velocity ];
0219 
0220 
0221 
0222 %
0223 % Find the indices of data values that have just greater depth than
0224 % zdepth.
0225 zgreater = ( z < zdepth );
0226 zg_ind = diff(zgreater);
0227 zg_ind = find(zg_ind~=0);
0228 zg_ind = zg_ind + [0:1:length(zg_ind)-1]';
0229 velocity_greater_z = velocity(zg_ind);
0230 depth_greater_z = z(zg_ind);
0231 
0232 
0233 %
0234 % Find the indices of the data values that have just lesser depth
0235 % than zdepth.
0236 zlesser = ( z > zdepth );
0237 zl_ind = diff(zlesser);
0238 zl_ind = find(zl_ind~=0);
0239 zl_ind = zl_ind + [1:1:length(zg_ind)]';
0240 velocity_lesser_z = velocity(zl_ind);
0241 depth_lesser_z = z(zl_ind);
0242 
0243 
0244 %
0245 % Interpolate between the data values.
0246 alpha = (zdepth - depth_greater_z) ./ ( depth_lesser_z - depth_greater_z );
0247 
0248 velocity_at_depth = (velocity_lesser_z .* alpha) + (velocity_greater_z .* (1-alpha));
0249 
0250 
0251 %
0252 % Now reshape the data into a square grid.
0253 w = reshape ( velocity_at_depth, [eta_psi_length xi_psi_length] );
0254 
0255 
0256 
0257 
0258 w = w.'; x = x.'; y = y.';
0259 
0260 % apply the rotation
0261 %
0262 % Since the "angle" variable is not always present, we
0263 % need to construct it.  In the scrum files, it is apparently
0264 % always in degrees.  Don't bother with degrees here.
0265 angle = zeros(size(x));
0266 
0267 
0268 
0269 [r,c] = size(x);
0270 j = [2:c-1];
0271 for i = 2:r-1
0272     angle(i,j) = atan2(y(i+1,j)-y(i-1,j), x(i+1,j)-x(i-1,j));
0273 end
0274 
0275 % rotate into east and north components
0276  w=w.*exp(sqrt(-1)*angle);
0277 
0278 
0279 %
0280 % Mask out the land.
0281 % Since we average the u and v across cells, the appropriate mask
0282 % dimension is that of "mask_psi".  I don't want to take the chance
0283 % that it may not be present in the file, so I compute it from
0284 % mask_rho.
0285 [umask,vmask,pmask]=uvp_masks ( mask_rho );
0286 mask_inds = find ( pmask == 0 );
0287 w(mask_inds) = NaN * ones(size(mask_inds));
0288 
0289