function [z_r,z_w,Hz] = get_z(tidx, ncfile)

%function [z_r,z_w,Hz]=get_z(tidx, ncfile)
% z_r is z at rho points
% z_w is a at w points
% tidx is mosdel time step

% Added: cdenham.

if ischar(ncfile)
    nc = netcdf(ncfile);   % Assume it works.
else   % Assume it is a "netcdf" object already.
    nc = ncfile;
end

% End: cdenham.

Tcline = nc{'Tcline'}(:);
h    = nc{'h'}(:);
zz=max(size(h));
if zz==0
  h    = squeeze(nc{'bath'}(tidx,:,:));
end
zz=size(h);
if length(zz)>2
  h    = squeeze(nc{'h'}(tidx,:,:));
end

hmin = min(h(:));
zeta =squeeze(nc{'zeta'}(tidx,:,:));
hc   =nc{'hc'}(:);

%
%  Get z at w points
%
sc_w = [nc{'sc_w'}(:)];
if isempty(sc_w); sc_w=nc{'s_w'}(:); end
Cs_w = [nc{'Cs_w'}(:)];
for  k = 1:length(sc_w)
  z_w(k,:,:) = zeta.*(1+sc_w(k)*hc./h-hc*Cs_w(k)./h+Cs_w(k))+hc*sc_w(k)+(h - hc)*Cs_w(k);
end
%z_w(1,:,:)=-h(:,:);

%
%  Get z at rho points
%
sc_r =nc{'sc_r'}(:);
if isempty(sc_r); sc_r=nc{'s_rho'}(:); end
Cs_r =nc{'Cs_r'}(:);
for  k = 1:length(sc_r)
  z_r(k,:,:) = zeta.*(1+sc_r(k)*hc./h-hc*Cs_r(k)./h+Cs_r(k))+hc*sc_r(k)+(h - hc)*Cs_r(k);
  Hz(k,:,:) = z_w(k+1,:,:)-z_w(k,:,:);
end


% Added: cdenham.

if ischar(ncfile)
    close(nc)   % Assume it works.
end

% End: cdenham.