function [z]=oct_get_depths(fname,gname,tindex,type)
% OCT_GET_DEPTHS Procedural variant of get_depths using octave-oct_netcdf
% Read bathymetry from grid file
ngid = oct_netcdf.open(gname,'NC_NOWRITE');
try
  h = oct_netcdf.getVar(ngid,oct_netcdf.inqVarID(ngid,'h'));
catch
  h = [];
end
oct_netcdf.close(ngid);

ncid = oct_netcdf.open(fname,'NC_NOWRITE');
% read zeta and hmorph
try
  zeta = squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'zeta')));
  if ndims(zeta)==3, zeta = squeeze(zeta(tindex,:,:)); end
catch
  zeta = [];
end
try
  hmorph = squeeze(oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'hmorph')));
  if ndims(hmorph)==3, hmorph = squeeze(hmorph(tindex,:,:)); end
catch
  hmorph = [];
end
if ~isempty(hmorph), h=hmorph; end

% read vertical coordinate parameters
theta_s = []; theta_b = []; Tcline = [];
try, theta_s = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'theta_s')); end
if isempty(theta_s)
  try, theta_s = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'theta_s')); end
end
try, theta_b = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'theta_b')); end
try, Tcline = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'Tcline')); end
try, hc = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'hc')); catch, hc = [] ; end

if isempty(Tcline)
  if ~isempty(h)
    hmin = min(min(h));
    if isempty(Tcline), Tcline = hc; end
    if isempty(hc), hc = Tcline; end
    if isempty(hc), hc = min(hmin,0); end
  end
end

% determine N from dimension s_rho
try
  did = oct_netcdf.inqDimID(ncid,'s_rho'); [~,N] = oct_netcdf.inqDim(ncid,did);
catch
  % fallback: look for s_rho variable
  try, s_rho = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'s_rho')); N = length(s_rho); catch, N = 0; end
end

% Vertical transform detection
s_coord = 1;
try
  VertCoordType = oct_netcdf.getAtt(ncid,oct_netcdf.getConstant('NC_GLOBAL'),'VertCoordType');
catch
  VertCoordType = '';
end
if isempty(VertCoordType)
  try
    vtrans = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'Vtransform'));
    if ~isempty(vtrans), s_coord = vtrans; end
  catch
  end
elseif strcmp(VertCoordType,'NEW')
  s_coord = 2;
end
if s_coord==2, hc = Tcline; end

oct_netcdf.close(ncid);

if isempty(zeta)
  zeta = 0 .* h;
end

vtype = type;
if (type=='u') || (type=='v'), vtype='r'; end
z = zlevs(h,zeta,theta_s,theta_b,hc,N,vtype,s_coord);
if type=='u', z = rho2u_3d(z); end
if type=='v', z = rho2v_3d(z); end
return