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

ncid = netcdf.open(fname,'NC_NOWRITE');
% read zeta and hmorph
try
  zeta = squeeze(netcdf.getVar(ncid,netcdf.inqVarID(ncid,'zeta')));
  if ndims(zeta)==3, zeta = squeeze(zeta(tindex,:,:)); end
catch
  zeta = [];
end
try
  hmorph = squeeze(netcdf.getVar(ncid,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 = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_s')); end
if isempty(theta_s)
  try, theta_s = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_s')); end
end
try, theta_b = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'theta_b')); end
try, Tcline = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'Tcline')); end
try, hc = netcdf.getVar(ncid,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 = netcdf.inqDimID(ncid,'s_rho'); [~,N] = netcdf.inqDim(ncid,did);
catch
  % fallback: look for s_rho variable
  try, s_rho = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'s_rho')); N = length(s_rho); catch, N = 0; end
end

% Vertical transform detection
s_coord = 1;
try
  VertCoordType = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'VertCoordType');
catch
  VertCoordType = '';
end
if isempty(VertCoordType)
  try
    vtrans = netcdf.getVar(ncid,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

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
