scrum_zslice

SCRUM_ZSLICE:  Returns horizontal slice from SCRUM Netcdf file.

 Returns a matrix containing a horizontal slice 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: [u,x,y]=scrum_zslice(cdf,var,time,depth)
    cdf:     name of SCRUM/ROMS NetCDF file (or object).
    var:     name of SCRUM variable.
    time:    time index, must be zero-based
    depth:   depth in meters.

0001 function [u,x,y]=scrum_zslice(cdf,var,timestep,depth)
0002 %SCRUM_ZSLICE:  Returns horizontal slice from SCRUM Netcdf file.
0003 %
0004 % Returns a matrix containing a horizontal slice at a specified depth
0005 % at a given time step from a SCRUM NetCDF file.  Regions of grid that
0006 % are shallower than requested value are returned as NaNs.
0007 %
0008 % USAGE: [u,x,y]=scrum_zslice(cdf,var,time,depth)
0009 %    cdf:     name of SCRUM/ROMS NetCDF file (or object).
0010 %    var:     name of SCRUM variable.
0011 %    time:    time index, must be zero-based
0012 %    depth:   depth in meters.
0013 %
0014 
0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0016 % Load in NetCDF file
0017 if (nargin~=4)
0018   help scrum_zslice; return
0019 end
0020 if isa(cdf, 'netcdf')
0021   nc = cdf;
0022 elseif isa(cdf, 'char')
0023   nc = netcdf(cdf,'nowrite');
0024 end
0025 if(isempty(nc))
0026   disp(['ERROR:  File ',cdfin,' not found.']);
0027   return;
0028 end
0029 ncmex('setopts',0);
0030 nc=quick(nc,1);
0031 
0032 ncid=ncmex('open',cdf,'nowrite');
0033 if(ncid==-1),
0034   disp(['file ' cdf ' not found'])
0035   return
0036 end
0037 depth = -abs(depth);             % account for depth or z
0038 
0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0040 % Acquire the grid. If "lon_rho" and "lat_rho" are present, grab them,
0041 % otherwise get "x_rho" and "y_rho".
0042 [lon_rho_varid, rcode] = ncmex('VARID', ncid, 'lon_rho');
0043 [lat_rho_varid, rcode] = ncmex('VARID', ncid, 'lat_rho');
0044 if ( (lon_rho_varid >= 0) | (lat_rho_varid >= 0) )
0045     x=ncmex('varget',ncid,'lon_rho',[0 0],[-1 -1]);
0046     y=ncmex('varget',ncid,'lat_rho',[0 0],[-1 -1]);
0047 else
0048     x=ncmex('varget',ncid,'x_rho',[0 0],[-1 -1]);
0049     y=ncmex('varget',ncid,'y_rho',[0 0],[-1 -1]);
0050 end
0051 x = x';
0052 y = y';
0053 [eta_rho_length, xi_rho_length] = size(x);
0054 
0055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0056 % The grid upon which the requested variable is based may not be
0057 % [eta_rho xi_rho].  If the requested variable is, for example, 'u', then
0058 % the grid will be [eta_u xi_u].  If this is the case, then x and y
0059 % need to be altered to reflect this.
0060 [varid, rcode] = ncmex('VARID', ncid, var);
0061 [dud, dud, ndims, vardims, dud, status] = ncmex('varinq', ncid, varid);
0062 
0063 y_dimid = vardims(ndims-1);
0064 x_dimid = vardims(ndims);
0065 
0066 [y_dim_name, y_length, status] = ncmex('diminq', ncid, y_dimid);
0067 [x_dim_name, x_length, status] = ncmex('diminq', ncid, x_dimid);
0068 
0069 [y_length, x_length] = size(x);
0070 if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
0071     ;    
0072 elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
0073     x_length = x_length-1;
0074     eta_u_length = eta_rho_length;
0075     xi_u_length = xi_rho_length-1;
0076     x = (x(:,1:xi_rho_length-1) + x(:,2:xi_rho_length))/2;
0077     y = (y(:,1:xi_rho_length-1) + y(:,2:xi_rho_length))/2;
0078 elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
0079     y_length = y_length-1;
0080     eta_v_length = eta_rho_length-1;
0081     xi_v_length = xi_rho_length;
0082     x = (x(1:eta_rho_length-1,:) + x(2:eta_rho_length,:))/2;
0083     y = (y(1:eta_rho_length-1,:) + y(2:eta_rho_length,:))/2;
0084 else
0085     disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
0086     help scrum_zslice;
0087     return;
0088 end
0089 
0090 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0091 % Get variables for z eqn [z = zeta * (1 + s) + hc*s + (h - hc)*C(s)]
0092 
0093 [s_rho_dimid, status] = ncmex ( 'dimid', ncid, 's_rho' );
0094 if ( status == -1 )
0095     fprintf ( 2, 'Could not get s_rho dimid from %s.\n', cdf );
0096     ncmex ( 'close', ncid );
0097     return;
0098 end
0099 [dimname, s_rho_length, status] = ncmex ( 'diminq', ncid, s_rho_dimid );
0100 if ( status == -1 )
0101     fprintf ( 2, 'Could not get s_rho length from %s.\n', cdf );
0102     ncmex ( 'close', ncid );
0103     return;
0104 end
0105 
0106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0107 % w is defined at different locations than the others
0108 if ( strcmp(var,'w') )
0109     [sc, status] = ncmex ( 'varget', ncid, 'sc_w', [0], [-1] );
0110 else
0111     [sc, status] = ncmex ( 'varget', ncid, 'sc_r', [0], [-1] );
0112 end
0113 
0114 [zeta, status] = ncmex ( 'varget', ncid, 'zeta', [timestep 0 0], [1 -1 -1] );
0115 zeta = zeta';
0116 
0117 [hc_varid, status] = ncmex ( 'varid', ncid, 'hc' );
0118 if ( status == -1 )
0119     fprintf ( 2, 'Could not get hc varid from %s.\n', cdf );
0120     ncmex ( 'close', ncid );
0121     return;
0122 end
0123 [hc, status] = ncmex ( 'varget1', ncid, hc_varid, [0] );
0124 if ( status == -1 )
0125     fprintf ( 2, 'Could not get hc from %s.\n', cdf );
0126     ncmex ( 'close', ncid );
0127     return;
0128 end
0129 
0130 [h, status] = ncmex ( 'varget', ncid, 'h', [0 0], [-1 -1] );
0131 if ( status == -1 )
0132     fprintf ( 'scrum_zslice:  could not get ''h'' in %s.', cdf );
0133     return;
0134 end
0135 h = h';
0136 
0137 [Cs_r, status] = ncmex ( 'varget', ncid, 'Cs_r', [0], [-1] );
0138 if ( status == -1 )
0139     fprintf ( 'scrum_zslice:  could not get ''Cs_r'' in %s.', cdf );
0140     return;
0141 end
0142 
0143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0144 % Construct the depth.
0145 n = length(sc);
0146 z = zeta * (1+sc(1)) + hc*sc(1) + (h - hc)*Cs_r(1);
0147 for i = 2:n
0148   zi = zeta * (1+sc(i)) + hc*sc(i) + (h - hc)*Cs_r(i);
0149   z = cat ( 3, z, zi );
0150 end
0151 z = permute ( z, [3 1 2] );
0152 
0153 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0154 % reshape the depth such that each vertical profile is a column.
0155 if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
0156     ;
0157 elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
0158     z = (z(:,:,1:xi_rho_length-1) + z(:,:,2:xi_rho_length))/2;
0159     zeta = (zeta(:,1:xi_rho_length-1) + zeta(:,2:xi_rho_length))/2;
0160     h = (h(:,1:xi_rho_length-1) + h(:,2:xi_rho_length))/2;
0161 elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
0162     z = (z(:,1:eta_rho_length-1,:) + z(:,2:eta_rho_length,:))/2;
0163     zeta = (zeta(1:eta_rho_length-1,:) + zeta(2:eta_rho_length,:))/2;
0164     h = (h(1:eta_rho_length-1,:) + h(2:eta_rho_length,:))/2;
0165 else
0166     disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
0167     help scrum_zslice;
0168     return;
0169 end
0170 z = reshape ( z, [s_rho_length y_length*x_length] );
0171 
0172 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0173 % Retrieve the requested variable at the given timestep.
0174 [data,status]=ncmex('varget',ncid,var,[(timestep-1) 0 0 0],[1 -1 -1 -1] );
0175 if ( status == -1 )
0176     fprintf ( 'scrum_zslice:  could not get ''%s'' in %s.', var, cdf );
0177     return;
0178 end
0179 data = permute ( data, [3 2 1] );
0180 data = reshape(data,[s_rho_length y_length*x_length]);
0181 
0182 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0183 % define z and data at zeta (elevation)
0184 z = [z; reshape(zeta, [1 y_length*x_length] ) ];
0185 data = [data; data(s_rho_length,:)];
0186 
0187 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0188 % define z at values greater than zeta to be Inf
0189 % define data to be NaN
0190 z = [z; Inf*ones(1,y_length*x_length)];
0191 data = [data; NaN * ones(1,y_length*x_length)];
0192 
0193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0194 % define data at bathymetry values to be same as data at lowest
0195 % existing values.
0196 z = [ -1*reshape(h,[1 y_length*x_length]); z];
0197 data = [ data(1,:); data ];
0198 
0199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0200 % define z at depth greater than bathymetry to be -Inf
0201 % define data at depth greater than bathymetry to be NaN
0202 z = [ -Inf*ones(1,y_length*x_length); z ];
0203 data = [ NaN * ones(1,y_length*x_length); data ];
0204 
0205 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0206 % Find the indices of data values that have just greater depth than
0207 % depth.
0208 zgreater = ( z < depth );
0209 zg_ind = diff(zgreater);
0210 zg_ind = find(zg_ind~=0);
0211 zg_ind = zg_ind + [0:1:length(zg_ind)-1]';
0212 data_greater_z = data(zg_ind);
0213 depth_greater_z = z(zg_ind);
0214 
0215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0216 % Find the indices of the data values that have just lesser depth
0217 % than depth.
0218 zlesser = ( z > depth );
0219 zl_ind = diff(zlesser);
0220 zl_ind = find(zl_ind~=0);
0221 zl_ind = zl_ind + [1:1:length(zg_ind)]';
0222 data_lesser_z = data(zl_ind);
0223 depth_lesser_z = z(zl_ind);
0224 
0225 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0226 % Interpolate between the data values.
0227 alpha = (depth - depth_greater_z) ./ ( depth_lesser_z - depth_greater_z );
0228 
0229 data_at_depth = (data_lesser_z .* alpha) + (data_greater_z .* (1-alpha));
0230 
0231 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0232 % Now reshape the data into a square grid.
0233 u = reshape ( data_at_depth, [y_length x_length] );
0234 
0235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0236 % If there is a suitable mask variable, we will want to mask out certain cells.
0237 % For now, mask out the land cells.
0238 if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
0239   %Get the mask_rho variable
0240   [mask_rho_varid, status] = ncmex ( 'varid', ncid, 'mask_rho' );
0241   if ( status ~= -1 )
0242     [mask_rho, status] = ncmex ( 'varget', ncid, 'mask_rho', [0 0], [-1 -1] );
0243     mask_rho = mask_rho';
0244     mask_inds = find(mask_rho==0);
0245     u(mask_inds) = NaN * ones(size(mask_inds));
0246   end
0247 elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
0248   %
0249   %Get the mask_u variable
0250   [mask_u_varid, status] = ncmex ( 'varid', ncid, 'mask_u' );
0251   if ( status ~= -1 )
0252     [mask_u, status] = ncmex ( 'varget', ncid, 'mask_u', [0 0], [-1 -1] );
0253     mask_u = mask_u';
0254     mask_inds = find(mask_u==0);
0255     u(mask_inds) = NaN * ones(size(mask_inds));
0256   end
0257 elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
0258   %Get the mask_u variable
0259   [mask_v_varid, status] = ncmex ( 'varid', ncid, 'mask_u' );
0260   if ( status ~= -1 )
0261     [mask_v, status] = ncmex ( 'varget', ncid, 'mask_v', [0 0], [-1 -1] );
0262     mask_v = mask_v';
0263     mask_inds = find(mask_v==0);
0264     u(mask_inds) = NaN * ones(size(mask_inds));
0265   end
0266 else
0267   disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
0268   help scrum_zslice;
0269   return;
0270 end
0271 
0272 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0273 % all done with netcdf file operations
0274 ncmex ( 'close', ncid );
0275 
0276 u = u'; x = x'; y = y';