function oct_ext_tracers_ini(ininame,grdname,seas_datafile,ann_datafile,...
                         dataname,vname,type,tini);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% P. Marchesiello - 2005. Adapted from P. Penven's oct_ext_tracers.m 
%
%  Ext tracers in a CROCO initial file
%  take seasonal data for the upper levels and annual data for the
%  lower levels
%
%  input:
%    ininame       : CROCO initial file name
%    grdname       : CROCO grid file name    
%    seas_datafile : regular longitude - latitude - z seasonal data 
%                    file used for the upper levels  (oct_netcdf)
%    ann_datafile  : regular longitude - latitude - z annual data 
%                    file used for the lower levels  (oct_netcdf)
%    dataname      : variable name in data file
%    vname         : variable name in CROCO file
%    type          : position on C-grid ('r', 'u', 'v', 'p')
%    tini          : initialisation time [days]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
%
% set the value of ro (oa decorrelation scale [m]) 
% and default (value if no data)
%
ro=0;
default=NaN;
disp([' Ext tracers: ro = ',num2str(ro/1000),...
      ' km - default value = ',num2str(default)])

% Open initial file
%
ncid = netcdf.open(ininame, 'NC_WRITE');
theta_s = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'theta_s'));
theta_b =  netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'theta_b'));
hc  =  netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'hc'));
N =  netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 's_rho'));
vtransform = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'Vtransform'));
if  ~exist('vtransform')
    vtransform=1; %Old Vtransform
    disp([' NO VTRANSFORM parameter found'])
    disp([' USE TRANSFORM default value vtransform = 1'])
end
%
% Open and Read grid file  
% 
ng_id = netcdf.open(grdname, 'NC_NOWRITE');
lon=netcdf.getVar(ng_id, netcdf.inqVarID(ng_id, 'lon_rho'));
lat=netcdf.getVar(ng_id, netcdf.inqVarID(ng_id, 'lat_rho'));
h=netcdf.getVar(ng_id, netcdf.inqVarID(ng_id, 'h'));
netcdf.close(ng_id);
[M,L]=size(lon);
%
% Read seasonal datafile 
%
ncidseas = netcdf.open(seas_datafile, 'NC_NOWRITE');
X=netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'X'));
Y=netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'Y'));
Zseas=-netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'Z'));
T=netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'T')).*30;
tlen=length(T);
Nzseas=length(Zseas);
%
% Read annual datafile
%
ncidann = netcdf.open(ann_datafile, 'NC_NOWRITE');
Zann=-netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, '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(['   oct_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
%------------------------------------------------------------
%
%
% Interpole annual dataset on horizontal CROCO grid
%
if Nz > Nzseas
  if Zseas~=Zann(1:length(Zseas)) 
    error('vertical levels dont match')
  end
  datazgrid=zeros(Nz,M,L);
  missval=netcdf.getAtt(ncidann, netcdf.inqVarID(ncidann, dataname), 'missing_value');
  for k=Nzseas+1:Nz
    if ~isempty(i2)
      data=squeeze(netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, dataname)));
    else
      data=[];
    end
    if ~isempty(i1)
      data=cat(2,squeeze(netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, dataname))),data);
    end
    if ~isempty(i3)
      data=cat(2,data,squeeze(netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, dataname))));
    end
    data=oct_get_missing_val(x,y,data,missval,ro,default);
    datazgrid(k,:,:)=interp2(x,y,data,lon,lat,'cubic');
  end
end
netcdf.close(ncidann);
%
% interpole seasonal dataset on horizontal croco grid
%
disp(['   oct_ext_tracers_ini: horizontal interpolation of seasonal data'])
missval=netcdf.getAtt(ncidseas, netcdf.inqVarID(ncidseas, dataname), 'missing_value');
if Nz <= Nzseas
  datazgrid=zeros(Nz,M,L);
end
for k=1:min([Nz Nzseas])
  if ~isempty(i2)
    data=squeeze(netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, dataname)));
  else
    data=[];
  end
  if ~isempty(i1)
    data=cat(2,squeeze(netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, dataname))),data);
  end
  if ~isempty(i3)
    data=cat(2,data,squeeze(netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, dataname))));
  end
  data=oct_get_missing_val(x,y,data,missval,ro,default);
  datazgrid(k,:,:)=interp2(x,y,data,lon,lat,'cubic');
end
netcdf.close(ncidseas);
%
%----------------------------------------------------
%  Vertical interpolation
%-----------------------------------------------------
%
disp('   oct_ext_tracers_ini: vertical interpolation')
%
% Get the sigma depths
%
zcroco=oct_zlevs(h,0.*h,theta_s,theta_b,hc,N,'r',vtransform);
if type=='u'
  zcroco=oct_rho2u_3d(zcroco);
end
if type=='v'
  zcroco=oct_rho2v_3d(zcroco);
end
zmin=min(min(min(zcroco)));
zmax=max(max(max(zcroco)));
%
% Check if the min z level is below the min sigma level
%    (if not add a deep layer)
%
z=Zann;
addsurf=max(z)<zmax;
addbot=min(z)>zmin;
if addsurf
 z=[100;z];
end
if addbot
 z=[z;-100000];
end
Nz=min(find(z<zmin));
z=z(1:Nz);
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,:,:);
%
% Do the vertical interpolation and write in inifile
%
netcdf.putVar(ncid, netcdf.inqVarID(ncid, vname), 1,:,:,:-1, 1, oct_ztosigma(flipdim(var,1),zcroco,flipud(z)));  % [conv] 0-based
netcdf.close(ncid);

return