function data=oct_ext_data(datafile,dataname,tindex,lon,lat,time,Roa,savefile)
% oct_ext_data: procedural replacement for ext_data using octave-netcdf
if nargin < 8 
  savefile=2;
end

disp(['Getting ',dataname,' for time index ',num2str(tindex)])
default=NaN;
dl=1;
lonmin=min(min(lon))-dl;
lonmax=max(max(lon))+dl;
latmin=min(min(lat))-dl;
latmax=max(max(lat))+dl;

ncid = netcdf.open(datafile,'NC_NOWRITE');
% get var id for dataname
try
  varid = netcdf.inqVarID(ncid,dataname);
catch
  netcdf.close(ncid);
  error(['Variable ',dataname,' not found in ',datafile]);
end
% get attributes safely
missval = NaN;
try
  missval = netcdf.getAtt(ncid,varid,'missing_value');
catch
end
add_offset = [];
scale_factor = [];
try
  add_offset = netcdf.getAtt(ncid,varid,'add_offset');
  scale_factor = netcdf.getAtt(ncid,varid,'scale_factor');
catch
end

% find lon var
X = [];
for name = {'X','lon','longitude'}
  try
    vid = netcdf.inqVarID(ncid,name{1});
    X = netcdf.getVar(ncid,vid);
    break
  catch
  end
end
if isempty(X)
  netcdf.close(ncid);
  error(['Empty longitude in ',datafile])
end
% find lat var
Y = [];
for name = {'Y','lat','latitude'}
  try
    vid = netcdf.inqVarID(ncid,name{1});
    Y = netcdf.getVar(ncid,vid);
    break
  catch
  end
end
if isempty(Y)
  netcdf.close(ncid);
  error(['Empty latitude in ',datafile])
end

% get a subgrid
j=find(Y>=latmin & Y<=latmax);
i1=find(X-360>=lonmin & X-360<=lonmax);
i2=find(X>=lonmin & X<=lonmax);
i3=find(X+360>=lonmin & X+360<=lonmax);
x=cat(1,X(i1)-360,X(i2),X(i3)+360);
y=Y(j);

% determine ndims by inquiring var
[~,xtype,dimids,~] = netcdf.inqVar(ncid,varid);
ndims = length(dimids);

data = [];
if ~isempty(i2)
  if ndims==3
    % (tindex, j, i2) -> zero-based start/count
    start = [tindex-1, j(1)-1, i2(1)-1];
    count = [1, length(j), length(i2)];
    data = squeeze(netcdf.getVar(ncid,varid,start,count));
  elseif ndims==4
    start = [tindex-1, 0, j(1)-1, i2(1)-1];
    count = [1, 1, length(j), length(i2)];
    data = squeeze(netcdf.getVar(ncid,varid,start,count));
  else
    netcdf.close(ncid);
    error(['Bad dimension number ',num2str(ndims)])
  end
end
if ~isempty(i1)
  if ndims==3
    start = [tindex-1, j(1)-1, i1(1)-1];
    count = [1, length(j), length(i1)];
    d1 = squeeze(netcdf.getVar(ncid,varid,start,count));
  elseif ndims==4
    start = [tindex-1, 0, j(1)-1, i1(1)-1];
    count = [1,1,length(j), length(i1)];
    d1 = squeeze(netcdf.getVar(ncid,varid,start,count));
  else
    netcdf.close(ncid);
    error(['Bad dimension number ',num2str(ndims)])
  end
  data = cat(2,d1,data);
end
if ~isempty(i3)
  if ndims==3
    start = [tindex-1, j(1)-1, i3(1)-1];
    count = [1, length(j), length(i3)];
    d3 = squeeze(netcdf.getVar(ncid,varid,start,count));
  elseif ndims==4
    start = [tindex-1, 0, j(1)-1, i3(1)-1];
    count = [1,1,length(j), length(i3)];
    d3 = squeeze(netcdf.getVar(ncid,varid,start,count));
  else
    netcdf.close(ncid);
    error(['Bad dimension number ',num2str(ndims)])
  end
  data = cat(2,data,d3);
end

netcdf.close(ncid);

% Perform the extrapolation
if savefile==2
  data=get_missing_val(x,y,data,missval,Roa,default,2);
else
  if tindex==1
    data=get_missing_val(x,y,data,missval,Roa,default,1);
  else
    data=get_missing_val(x,y,data,missval,Roa,default,0);
  end
end
% Interpolation on the CROCO grid
data=interp2(x,y,data,lon,lat,'cubic');
% Apply offset
if ~isempty(add_offset)
  data=add_offset+data*scale_factor;
end
% Salinity test
if strcmp(dataname,'salinity')
  data(data<2)=2;
end
return
