function [u,x,y]=scrum_zslice(cdf,var,timestep,depth)
%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.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load in NetCDF file
if (nargin~=4)
  help scrum_zslice; return
end
if isa(cdf, 'netcdf')
  nc = cdf;
elseif isa(cdf, 'char')
  nc = netcdf(cdf,'nowrite');
end
if(isempty(nc))
  disp(['ERROR:  File ',cdfin,' not found.']);
  return;
end
ncmex('setopts',0);
nc=quick(nc,1);

ncid=ncmex('open',cdf,'nowrite');
if(ncid==-1),
  disp(['file ' cdf ' not found'])
  return
end
depth = -abs(depth); 			% account for depth or z

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Acquire the grid. If "lon_rho" and "lat_rho" are present, grab them,
% otherwise get "x_rho" and "y_rho".
[lon_rho_varid, rcode] = ncmex('VARID', ncid, 'lon_rho');
[lat_rho_varid, rcode] = ncmex('VARID', ncid, 'lat_rho');
if ( (lon_rho_varid >= 0) | (lat_rho_varid >= 0) )
    x=ncmex('varget',ncid,'lon_rho',[0 0],[-1 -1]);
    y=ncmex('varget',ncid,'lat_rho',[0 0],[-1 -1]);
else
    x=ncmex('varget',ncid,'x_rho',[0 0],[-1 -1]);
    y=ncmex('varget',ncid,'y_rho',[0 0],[-1 -1]);
end
x = x';
y = y';
[eta_rho_length, xi_rho_length] = size(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The grid upon which the requested variable is based may not be
% [eta_rho xi_rho].  If the requested variable is, for example, 'u', then
% the grid will be [eta_u xi_u].  If this is the case, then x and y
% need to be altered to reflect this.
[varid, rcode] = ncmex('VARID', ncid, var);
[dud, dud, ndims, vardims, dud, status] = ncmex('varinq', ncid, varid);

y_dimid = vardims(ndims-1);
x_dimid = vardims(ndims);

[y_dim_name, y_length, status] = ncmex('diminq', ncid, y_dimid);
[x_dim_name, x_length, status] = ncmex('diminq', ncid, x_dimid);

[y_length, x_length] = size(x);
if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
	;	
elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
	x_length = x_length-1;
	eta_u_length = eta_rho_length;
	xi_u_length = xi_rho_length-1;
	x = (x(:,1:xi_rho_length-1) + x(:,2:xi_rho_length))/2;
	y = (y(:,1:xi_rho_length-1) + y(:,2:xi_rho_length))/2;
elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
	y_length = y_length-1;
	eta_v_length = eta_rho_length-1;
	xi_v_length = xi_rho_length;
	x = (x(1:eta_rho_length-1,:) + x(2:eta_rho_length,:))/2;
	y = (y(1:eta_rho_length-1,:) + y(2:eta_rho_length,:))/2;
else
	disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
	help scrum_zslice;
	return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get variables for z eqn [z = zeta * (1 + s) + hc*s + (h - hc)*C(s)]

[s_rho_dimid, status] = ncmex ( 'dimid', ncid, 's_rho' );
if ( status == -1 )
	fprintf ( 2, 'Could not get s_rho dimid from %s.\n', cdf );
	ncmex ( 'close', ncid );
	return;
end
[dimname, s_rho_length, status] = ncmex ( 'diminq', ncid, s_rho_dimid );
if ( status == -1 )
	fprintf ( 2, 'Could not get s_rho length from %s.\n', cdf );
	ncmex ( 'close', ncid );
	return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% w is defined at different locations than the others
if ( strcmp(var,'w') )
	[sc, status] = ncmex ( 'varget', ncid, 'sc_w', [0], [-1] );
else
	[sc, status] = ncmex ( 'varget', ncid, 'sc_r', [0], [-1] );
end

[zeta, status] = ncmex ( 'varget', ncid, 'zeta', [timestep 0 0], [1 -1 -1] );
zeta = zeta';

[hc_varid, status] = ncmex ( 'varid', ncid, 'hc' );
if ( status == -1 )
	fprintf ( 2, 'Could not get hc varid from %s.\n', cdf );
	ncmex ( 'close', ncid );
	return;
end
[hc, status] = ncmex ( 'varget1', ncid, hc_varid, [0] );
if ( status == -1 )
	fprintf ( 2, 'Could not get hc from %s.\n', cdf );
	ncmex ( 'close', ncid );
	return;
end

[h, status] = ncmex ( 'varget', ncid, 'h', [0 0], [-1 -1] );
if ( status == -1 )
	fprintf ( 'scrum_zslice:  could not get ''h'' in %s.', cdf );
	return;
end
h = h';

[Cs_r, status] = ncmex ( 'varget', ncid, 'Cs_r', [0], [-1] );
if ( status == -1 )
	fprintf ( 'scrum_zslice:  could not get ''Cs_r'' in %s.', cdf );
	return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the depth.
n = length(sc);
z = zeta * (1+sc(1)) + hc*sc(1) + (h - hc)*Cs_r(1);
for i = 2:n
  zi = zeta * (1+sc(i)) + hc*sc(i) + (h - hc)*Cs_r(i);
  z = cat ( 3, z, zi );
end
z = permute ( z, [3 1 2] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reshape the depth such that each vertical profile is a column.
if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
	;
elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
	z = (z(:,:,1:xi_rho_length-1) + z(:,:,2:xi_rho_length))/2;
	zeta = (zeta(:,1:xi_rho_length-1) + zeta(:,2:xi_rho_length))/2;
	h = (h(:,1:xi_rho_length-1) + h(:,2:xi_rho_length))/2;
elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
	z = (z(:,1:eta_rho_length-1,:) + z(:,2:eta_rho_length,:))/2;
	zeta = (zeta(1:eta_rho_length-1,:) + zeta(2:eta_rho_length,:))/2;
	h = (h(1:eta_rho_length-1,:) + h(2:eta_rho_length,:))/2;
else
	disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
	help scrum_zslice;
	return;
end
z = reshape ( z, [s_rho_length y_length*x_length] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Retrieve the requested variable at the given timestep.
[data,status]=ncmex('varget',ncid,var,[(timestep-1) 0 0 0],[1 -1 -1 -1] );
if ( status == -1 )
	fprintf ( 'scrum_zslice:  could not get ''%s'' in %s.', var, cdf );
	return;
end
data = permute ( data, [3 2 1] );
data = reshape(data,[s_rho_length y_length*x_length]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define z and data at zeta (elevation)
z = [z; reshape(zeta, [1 y_length*x_length] ) ];
data = [data; data(s_rho_length,:)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define z at values greater than zeta to be Inf
% define data to be NaN
z = [z; Inf*ones(1,y_length*x_length)];
data = [data; NaN * ones(1,y_length*x_length)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define data at bathymetry values to be same as data at lowest 
% existing values.
z = [ -1*reshape(h,[1 y_length*x_length]); z];
data = [ data(1,:); data ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define z at depth greater than bathymetry to be -Inf
% define data at depth greater than bathymetry to be NaN
z = [ -Inf*ones(1,y_length*x_length); z ];
data = [ NaN * ones(1,y_length*x_length); data ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the indices of data values that have just greater depth than
% depth.
zgreater = ( z < depth );
zg_ind = diff(zgreater);
zg_ind = find(zg_ind~=0);
zg_ind = zg_ind + [0:1:length(zg_ind)-1]';
data_greater_z = data(zg_ind);
depth_greater_z = z(zg_ind);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the indices of the data values that have just lesser depth
% than depth.
zlesser = ( z > depth );
zl_ind = diff(zlesser);
zl_ind = find(zl_ind~=0);
zl_ind = zl_ind + [1:1:length(zg_ind)]';
data_lesser_z = data(zl_ind);
depth_lesser_z = z(zl_ind);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interpolate between the data values.
alpha = (depth - depth_greater_z) ./ ( depth_lesser_z - depth_greater_z );

data_at_depth = (data_lesser_z .* alpha) + (data_greater_z .* (1-alpha));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now reshape the data into a square grid.
u = reshape ( data_at_depth, [y_length x_length] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If there is a suitable mask variable, we will want to mask out certain cells.
% For now, mask out the land cells.
if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
  %Get the mask_rho variable
  [mask_rho_varid, status] = ncmex ( 'varid', ncid, 'mask_rho' );
  if ( status ~= -1 )
    [mask_rho, status] = ncmex ( 'varget', ncid, 'mask_rho', [0 0], [-1 -1] );
    mask_rho = mask_rho';
    mask_inds = find(mask_rho==0);
    u(mask_inds) = NaN * ones(size(mask_inds));
  end
elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
  %
  %Get the mask_u variable
  [mask_u_varid, status] = ncmex ( 'varid', ncid, 'mask_u' );
  if ( status ~= -1 )
    [mask_u, status] = ncmex ( 'varget', ncid, 'mask_u', [0 0], [-1 -1] );
    mask_u = mask_u';
    mask_inds = find(mask_u==0);
    u(mask_inds) = NaN * ones(size(mask_inds));
  end
elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
  %Get the mask_u variable
  [mask_v_varid, status] = ncmex ( 'varid', ncid, 'mask_u' );
  if ( status ~= -1 )
    [mask_v, status] = ncmex ( 'varget', ncid, 'mask_v', [0 0], [-1 -1] );
    mask_v = mask_v';
    mask_inds = find(mask_v==0);
    u(mask_inds) = NaN * ones(size(mask_inds));
  end
else
  disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
  help scrum_zslice;
  return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% all done with netcdf file operations
ncmex ( 'close', ncid );

u = u'; x = x'; y = y';