SCRUM_KSLICE: returns horizontal slice at particular layer.
Horizontal slice is at specified s_rho layer at given time-step for a
SCRUM file.
This function can also be used to read in 2D and 3D fields such as
bathymetry and heat_flux. The coordinates of u are returned as x and y.
USAGE: [u,x,y]=scrum_kslice(cdf,var,[time],[layer])
where
cdf: file name for netCDf file (e.g. 'scrum_c.nc')
var: the variable to select (eg. 'salt' for salinity)
time: time step
layer: in the NetCDF file, s_rho is dimensioned as [0 ... s_rho_length-1].
The layer variables which are a function of s_rho run from -1 to 0,
which is reversed from ECOM. In order to have conformity, we will
have the parameter 'layer' be 1-based, and have layer 1 correspond
to the top layer instead of the bottom.
Examples:
[s,x,y]=scrum_kslice('scrum_c.nc','salt',2,3);
returns the salinity field from the 3rd sigma level
at the 2nd time step.
[elev,x,y]=scrum_kslice('scrum_c.nc','zeta',4);
returns the elevation field from the 4th time step
[depth,x,y]=scrum_kslice('scrum_c.nc','h');
returns the depth field0001 function [u,x,y]=scrum_kslice(cdf,var,time,layer) 0002 %SCRUM_KSLICE: returns horizontal slice at particular layer. 0003 % 0004 % Horizontal slice is at specified s_rho layer at given time-step for a 0005 % SCRUM file. 0006 % 0007 % This function can also be used to read in 2D and 3D fields such as 0008 % bathymetry and heat_flux. The coordinates of u are returned as x and y. 0009 % 0010 % USAGE: [u,x,y]=scrum_kslice(cdf,var,[time],[layer]) 0011 % 0012 % where 0013 % cdf: file name for netCDf file (e.g. 'scrum_c.nc') 0014 % var: the variable to select (eg. 'salt' for salinity) 0015 % time: time step 0016 % layer: in the NetCDF file, s_rho is dimensioned as [0 ... s_rho_length-1]. 0017 % The layer variables which are a function of s_rho run from -1 to 0, 0018 % which is reversed from ECOM. In order to have conformity, we will 0019 % have the parameter 'layer' be 1-based, and have layer 1 correspond 0020 % to the top layer instead of the bottom. 0021 % 0022 % 0023 % Examples: 0024 % 0025 % [s,x,y]=scrum_kslice('scrum_c.nc','salt',2,3); 0026 % returns the salinity field from the 3rd sigma level 0027 % at the 2nd time step. 0028 % 0029 % [elev,x,y]=scrum_kslice('scrum_c.nc','zeta',4); 0030 % returns the elevation field from the 4th time step 0031 % 0032 % [depth,x,y]=scrum_kslice('scrum_c.nc','h'); 0033 % returns the depth field 0034 % 0035 if (nargin<2 | nargin>4), 0036 help scrum_kslice; return 0037 end 0038 0039 % turn off warnings from NetCDf 0040 mexcdf('setopts',0); 0041 0042 % 0043 % open existing file 0044 ncid=mexcdf('open',cdf,'nowrite'); 0045 if(ncid==-1), 0046 disp(['file ' cdf ' not found']) 0047 return 0048 end 0049 0050 0051 % 0052 % Acquire the grid. 0053 % 0054 % If "lon_rho" and "lat_rho" are present, grab them. 0055 % Otherwise, get "x_rho" and "y_rho". 0056 [lon_rho_varid, rcode] = ncmex('VARID', ncid, 'lon_rho'); 0057 [lat_rho_varid, rcode] = ncmex('VARID', ncid, 'lat_rho'); 0058 if ( (lon_rho_varid >= 0) | (lat_rho_varid >= 0) ) 0059 x=ncmex('varget',ncid,'lon_rho',[0 0],[-1 -1]); 0060 y=ncmex('varget',ncid,'lat_rho',[0 0],[-1 -1]); 0061 else 0062 x=ncmex('varget',ncid,'x_rho',[0 0],[-1 -1]); 0063 y=ncmex('varget',ncid,'y_rho',[0 0],[-1 -1]); 0064 end 0065 0066 0067 [xi_rho_length, eta_rho_length] = size(x); 0068 0069 0070 % 0071 % Determine which variables make up the grid. 0072 % This assumes that the variable being asked for is 0073 % defined something like this. 0074 % 0075 % salt ( time, depth, ydimension, xdimension ); 0076 % 0077 % The x dimension must be last, and the y dimension must 0078 % be second to last. 0079 [varid, rcode] = ncmex('VARID', ncid, var); 0080 [dud, dud, ndims, vardims, dud, status] = ncmex('varinq', ncid, varid); 0081 0082 y_dimid = vardims(ndims-1); 0083 x_dimid = vardims(ndims); 0084 0085 [dud, y_length, status] = ncmex('diminq', ncid, y_dimid); 0086 [dud, x_length, status] = ncmex('diminq', ncid, x_dimid); 0087 0088 % 0089 % if the lengths of the x and y dimensions do not match up with 0090 % the eta_rho and xi_rho lengths, then we must have something 0091 % that depends upon eta_u, xi_v, or similar. Adjust the grid 0092 % accordingly. 0093 if ( x_length == (xi_rho_length-1) ) 0094 x = (x(1:xi_rho_length-1,:) + x(2:xi_rho_length,:))/2; 0095 y = (y(1:xi_rho_length-1,:) + y(2:xi_rho_length,:))/2; 0096 end 0097 if ( y_length == (eta_rho_length-1) ) 0098 x = (x(:,1:eta_rho_length-1) + x(:,2:eta_rho_length))/2; 0099 y = (y(:,1:eta_rho_length-1) + y(:,2:eta_rho_length))/2; 0100 end 0101 0102 0103 0104 0105 % 0106 % allow for using kslice on 2D, 3D and 4D variables 0107 switch ( nargin ) 0108 0109 % 0110 % 2D variable, such as bathymetry (h) 0111 case 2 0112 [u,ierr]=ncmex('varget',ncid,var,[0 0],[-1 -1]); 0113 0114 % 0115 % 3D 0116 case 3 0117 [u,ierr]=ncmex('varget',ncid,var,[(time-1) 0 0],[1 -1 -1]); 0118 0119 % 0120 % 4D 0121 case 4 0122 [s_rho_dimid, status] = ncmex ( 'dimid', ncid, 's_rho' ); 0123 [dud, s_rho_length, status] = ncmex ( 'diminq', ncid, s_rho_dimid ); 0124 [u,ierr]=ncmex('varget',ncid,var,[(time-1) (s_rho_length-layer) 0 0],[1 1 -1 -1]); 0125 0126 end 0127 0128 0129 0130 % 0131 % If the appropriate mask is present, use it to mask out the land. 0132 % Find the correct mask by comparing the dimension ids of the 0133 % variable in question and that of each mask variable. 0134 mask_vars = ['mask_rho'; 'mask_u '; 'mask_v '; 'mask_psi' ]; 0135 for ind = 1:4 0136 mask_var = deblank(mask_vars(ind,:)); 0137 0138 [mask_varid, status] = ncmex ( 'varid', ncid, mask_var ); 0139 if ( mask_varid ~= -1 ) 0140 [dud, dud, ndims, mask_dimids, dud, status] = ncmex('varinq', ncid, mask_varid); 0141 0142 if ( ~isempty(find(vardims==mask_dimids(1))) & ~isempty(find(vardims==mask_dimids(2))) ) 0143 [mask, status] = ncmex ( 'varget', ncid, mask_var, [0 0], [-1 -1] ); 0144 land = find(mask==0); 0145 u(land) = u(land)*NaN; 0146 end 0147 0148 end 0149 0150 end 0151 0152 ncmex('close',ncid); 0153