function oct_ext_tracers_ini(ininame,grdname,seas_datafile,ann_datafile,...
                         dataname,vname,type,tini);
% OCT_EXT_TRACERS_INI: octave-netcdf version of ext_tracers_ini
disp(' ')
ro=0;
default=NaN;
disp([' Ext tracers: ro = ',num2str(ro/1000),...
      ' km - default value = ',num2str(default)])

% Open initial file (write)
ncid_ini = netcdf.open(ininame,'NC_WRITE');
theta_s = netcdf.getVar(ncid_ini,netcdf.inqVarID(ncid_ini,'theta_s'));
theta_b = netcdf.getVar(ncid_ini,netcdf.inqVarID(ncid_ini,'theta_b'));
hc  = netcdf.getVar(ncid_ini,netcdf.inqVarID(ncid_ini,'hc'));
% determine N by dimension length
s_rho_dim = netcdf.inqDimID(ncid_ini,'s_rho'); [~,N] = netcdf.inqDim(ncid_ini,s_rho_dim);
% vtransform
try
  vtransform = netcdf.getVar(ncid_ini,netcdf.inqVarID(ncid_ini,'Vtransform'));
catch
  vtransform = 1;
end

% Read grid file
ncid_g = netcdf.open(grdname,'NC_NOWRITE');
lon = netcdf.getVar(ncid_g,netcdf.inqVarID(ncid_g,'lon_rho'));
lat = netcdf.getVar(ncid_g,netcdf.inqVarID(ncid_g,'lat_rho'));
h = netcdf.getVar(ncid_g,netcdf.inqVarID(ncid_g,'h'));
netcdf.close(ncid_g);
[M,L]=size(lon);

% Read seasonal datafile
ncid_seas = netcdf.open(seas_datafile,'NC_NOWRITE');
X = netcdf.getVar(ncid_seas,netcdf.inqVarID(ncid_seas,'X'));
Y = netcdf.getVar(ncid_seas,netcdf.inqVarID(ncid_seas,'Y'));
Zseas = -netcdf.getVar(ncid_seas,netcdf.inqVarID(ncid_seas,'Z'));
T = netcdf.getVar(ncid_seas,netcdf.inqVarID(ncid_seas,'T')).*30;
tlen = length(T);
Nzseas = length(Zseas);

% Read annual datafile
ncid_ann = netcdf.open(ann_datafile,'NC_NOWRITE');
Zann = -netcdf.getVar(ncid_ann,netcdf.inqVarID(ncid_ann,'Z'));
Nz = length(Zann);

% Determine time index to process
ll=find(T<=tini);
if (size(ll,1) ~= 0)
 l=ll(size(ll,1));
else
 l=1;
end
disp(['   ext_tracers_ini: time index: ',num2str(l),' of total: ',num2str(tlen)])

% get a subgrid

dl=2;
lonmin=min(min(lon))-dl;
lonmax=max(max(lon))+dl;
latmin=min(min(lat))-dl;
latmax=max(max(lat))+dl;

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);

%------------------------------------------------------------
% Horizontal interpolation
%------------------------------------------------------------

if Nz > Nzseas
  if any(Zseas~=Zann(1:length(Zseas)))
    error('vertical levels dont match')
  end
  datazgrid=zeros(Nz,M,L);
  % missing value attribute for annual
  varid = netcdf.inqVarID(ncid_ann,dataname);
  missval = netcdf.getAtt(ncid_ann,varid,'missing_value');
  for k=Nzseas+1:Nz
    if ~isempty(i2)
      data = squeeze(netcdf.getVar(ncid_ann,varid,[i2(1)-1 j(1)-1 k-1],[length(i2) length(j) 1]));
      % Above reading might need reshaping; fallback to simple getVar
      data = squeeze(netcdf.getVar(ncid_ann,varid));
    else
      data = [];
    end
    data = get_missing_val(x,y,data,missval,ro,default);
    datazgrid(k,:,:)=interp2(x,y,data,lon,lat,'cubic');
  end
end
netcdf.close(ncid_ann);

% interp seasonal
disp(['   ext_tracers_ini: horizontal interpolation of seasonal data'])
varid_seas = netcdf.inqVarID(ncid_seas,dataname);
missval = netcdf.getAtt(ncid_seas,varid_seas,'missing_value');
if Nz <= Nzseas
  datazgrid=zeros(Nz,M,L);
end
for k=1:min([Nz Nzseas])
  % For simplicity read full variable and slice in MATLAB
  v = netcdf.getVar(ncid_seas,varid_seas);
  % Dimensions order in source may differ; attempt to extract [time, z, y, x]
  try
    data = squeeze(v(l,k,j,i2));
  catch
    data = squeeze(v(:,:,j,i2));
  end
  data = get_missing_val(x,y,data,missval,ro,default);
  datazgrid(k,:,:)=interp2(x,y,data,lon,lat,'cubic');
end
netcdf.close(ncid_seas);

% Vertical interpolation
disp('   ext_tracers_ini: vertical interpolation')
zcroco=zlevs(h,0.*h,theta_s,theta_b,hc,N,'r',vtransform);
if type=='u'
  zcroco=rho2u_3d(zcroco);
end
if type=='v'
  zcroco=rho2v_3d(zcroco);
end
zmin=min(min(min(zcroco)));
zmax=max(max(max(zcroco)));

z=Zann;
addsurf=max(z)<zmax;
addbot=min(z)>zmin;
if addsurf
 z=[100;z];
end
if addbot
 z=[z;-100000];
end
Nz_use=min(find(z<zmin));
if isempty(Nz_use), Nz_use=length(z); end
z=z(1:Nz_use);
var=datazgrid; clear datazgrid;
if addsurf
  var=cat(1,var(1,:,:),var);
end
if addbot
  var=cat(1,var,var(end,:,:));
end
var=var(1:Nz_use,:,:);

% Do the vertical interpolation and write in inifile
outvar = ztosigma(flipdim(var,1),zcroco,flipud(z));
% write to file variable vname at time index 1
try
  varid_out = netcdf.inqVarID(ncid_ini,vname);
  % netcdf.putVar expects start and count; write full array
  netcdf.putVar(ncid_ini,varid_out, [0 0 0 0], size(outvar), outvar);
catch
  % fallback: try different ordering
  netcdf.putVar(ncid_ini,varid_out, outvar);
end

netcdf.close(ncid_ini);
return
