function oct_nested_restart(child_grd,parent_rst,child_rst,...
                        vertical_correc,extrapmask)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Compute the restart file of the embedded grid
%
%  Further Information:  
%  http://www.croco-ocean.org
%  
%  This file is part of CROCOTOOLS
%
%  CROCOTOOLS is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published
%  by the Free Software Foundation; either version 2 of the License,
%  or (at your option) any later version.
%
%  CROCOTOOLS is distributed in the hope that it will be useful, but
%  WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
%  MA  02111-1307  USA
%
%  Copyright (c) 2004-2006 by Pierrick Penven 
%  e-mail:Pierrick.Penven@ird.fr  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Check the number of variables in the clim parents file
a_id = netcdf.open(parent_rst, 'NC_NOWRITE');
Vars = var(a_id);
Varnames = [ncnames(Vars)];
nvar=length(Varnames);
%
namebiol={''};
unitbiol={''};
namepisces={''};
unitpisces={''};
biol=0;
pisces=0;
%biol_NPZD=0;
%biol_NPZD_wO2=0;
%biol_N2PZD2=0;

if nvar <= 37
    disp(['No biology'])
elseif nvar == 42
    %biol_NPZD=1;
    namebiol={'NO3';'CHLA';'PHYTO';'ZOO';'DET'};
    unitbiol={'mMol N m-3' ; 'mg C l-1' ; 'mMol N m-3';...
        'mMol N m-3' ; 'mMol N m-3'};
    biol=1;
    disp('Biol NPZD is on')
    disp('==========')
elseif nvar == 43
    %biol_NPZD_wO2=1;
    namebiol={'NO3';'O2';'CHLA';'PHYTO';'ZOO';'DET'};
    unitbiol={'mMol N m-3' ;'mMol O m-3' ; 'mg C l-1' ; 'mMol N m-3';...
        'mMol N m-3' ; 'mMol N m-3'};
    biol=1;
    disp('Biol NPZD_wO2 is on')
    disp('==========')
elseif nvar == 44
    %biol_N2PZD2=1;
    namebiol={'NO3' ; 'NH4' ; 'CHLA' ; 'PHYTO' ; 'ZOO' ; 'SDET' ; 'LDET'};
    unitbiol={'mMol N m-3' ; 'mMol N m-3' ; 'mg C l-1' ; ...
        'mMol N m-3' ; 'mMol N m-3' ; 'mMol N m-3' ; ...
        'mMol N m-3'};
    biol=1;
    disp('Biol N2PZD2 is on')
    disp('==========')
elseif nvar == 49
    %biol_bioebus=1;
    namebiol={'NO3';'NO2';'NH4';'SPHYTO';'LPHYTO';'SZOO';'LZOO';'SDET';'LDET';'DON';'O2';'N2O'};
    unitbiol={'mMol N m-3' ; 'mMol N m-3'  ; 'mMol N m-3' ; 'mMol N m-3' ; 'mMol N m-3' ; 'mMol N m-3'; ...
              'mMol N m-3' ; 'mMol N m-3'  ; 'mMol N m-3' ; 'mMol N m-3' ; 'mMol N m-3' ; 'mMol N m-3'};
    biol=1;
    disp('BioEBUS is on')
    disp('==========')   
else
    pisces=1;
    namepisces={'DIC' ; 'TALK' ; 'O2' ; 'CACO3'  ;  'PO4' ;  'POC' ;
        'Si'  ; 'NANO' ; 'ZOO'  ;  'DOC' ;  'DIA' ; 'MESO' ; 'BSI' ;
        'FER' ; 'BFE' ; 'GOC'; 'SFE' ; 'DFE' ; 'DSI' ; 'NFE' ; ...
        'NCHL' ;  'DCHL';  'NO3' ; 'NH4'};
    unitpisces={'umol C L-1' ; 'umol C L-1' ; 'umol O L-1' ;  'umol C L-1' ;  'umol P L-1' ;  'umol C L-1' ; ...
        'umol Si L-1' ;  'umol C L-1' ; 'umol C L-1'  ; 'umol C L-1' ; 'umol C L-1' ; 'umol C L-1' ; 'umol Si L-1' ; ...
        'umol Fe L-1' ; 'umol Fe L-1' ;  'umol C L-1' ;  'umol Fe L-1' ; 'umol Fe L-1' ; 'umol Si L-1';'umol Fe L-1' ;...
        'mg Chl m-3' ; 'mg Chl m-3' ; 'umol N L-1' ; 'umol N L-1'};
    disp('Pisces is on')
    disp('==========')
end

% if  (biol_NPZD | biol_NPZD_wO2 | biol_N2PZD2 )
%   biol=1;
% end

if extrapmask==1
  disp('Extrapolation under mask is on')
  disp('====================')
end

if vertical_correc==1
  disp('Vertical correction is on')
  disp('====================')
end

if pisces & (biol_NPZD | biol_NPZD_wO2 | biol_N2PZD2)
  error(['Both Biol NPZD type  and Pisces are ON, no possible yet !'])
end

%
% Title
%
title=['restart file for the embedded grid :',child_rst,...
       ' using oct_parent restart file: ',parent_rst];
disp(' ')
disp(title)
%
% Read in the embedded grid
%
disp(' ')
disp(' Read in the embedded grid...')
ncid = netcdf.open(child_grd, 'NC_NOWRITE');
parent_grd=ncid.parent_grid(:);
imin=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'grd_pos'));
imax=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'grd_pos'));
jmin=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'grd_pos'));
jmax=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'grd_pos'));
refinecoeff=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'refine_coef'));
result=close(ncid);
ncid = netcdf.open(parent_grd, 'NC_NOWRITE');
Lp=netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 'xi_rho'));
Mp=netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 'eta_rho'));
if extrapmask==1
  mask=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'mask_rho'));
else
  mask=[];
end
result=close(ncid);
%
% Read in the oct_parent restart file
%
disp(' ')
disp(' Read in the oct_parent restart file...')
ncid = netcdf.open(parent_rst, 'NC_NOWRITE');
theta_s=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'theta_s'));
vtransform=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'Vtransform'));
if isempty(theta_s)
  theta_s=ncid.theta_s(:);
  theta_b=ncid.theta_b(:);
  hc=ncid.hc(:);
else
  theta_b=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'theta_b'));
  hc=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'hc'));
end
if  ~exist('vtransform') | isempty(vtransform)
    vtransform=1; %Old Vtransform
    disp([' NO VTRANSFORM parameter found'])
    disp([' USE TRANSFORM default value vtransform = 1'])
end
N=netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 's_rho'));
thetime = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'scrum_time'));
result=close(ncid);
%
% Create the restart file
%
disp(' ')
disp(' Create the restart file...')
ncrst=oct_create_nestedrestart(child_rst,child_grd,parent_rst,title,...
			   'clobber',biol,pisces,namebiol,namepisces,unitbiol,unitpisces,hc,vtransform);
           
%
% Get the oct_parent indices
%
[igrd_r,jgrd_r]=meshgrid((1:1:Lp),(1:1:Mp));
[igrd_p,jgrd_p]=meshgrid((1:1:Lp-1),(1:1:Mp-1));
[igrd_u,jgrd_u]=meshgrid((1:1:Lp-1),(1:1:Mp));
[igrd_v,jgrd_v]=meshgrid((1:1:Lp),(1:1:Mp-1));
%
% Get the children indices
%
ipchild=(imin:1/refinecoeff:imax);
jpchild=(jmin:1/refinecoeff:jmax);
irchild=(imin+0.5-0.5/refinecoeff:1/refinecoeff:imax+0.5+0.5/refinecoeff);
jrchild=(jmin+0.5-0.5/refinecoeff:1/refinecoeff:jmax+0.5+0.5/refinecoeff);
[ichildgrd_p,jchildgrd_p]=meshgrid(ipchild,jpchild);
[ichildgrd_r,jchildgrd_r]=meshgrid(irchild,jrchild);
[ichildgrd_u,jchildgrd_u]=meshgrid(ipchild,jrchild);
[ichildgrd_v,jchildgrd_v]=meshgrid(irchild,jpchild);
%
% interpolations
% 
disp(' ')
disp(' Do the interpolations...')                     
np_id = netcdf.open(parent_rst, 'NC_NOWRITE');
for tindex=1:length(thetime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(thetime))])                     
  ncrst{'scrum_time'}(tindex)=thetime(tindex);
  ncrst{'time_step'}(tindex,:)= netcdf.getVar(np_id, netcdf.inqVarID(np_id, 'time_step'));
  %
  disp('zeta...')
  oct_interpvar3d(np_id,ncrst,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,'zeta',mask,tindex)
  disp('ubar...')
  oct_interpvar3d(np_id,ncrst,igrd_u,jgrd_u,ichildgrd_u,jchildgrd_u,'ubar',mask,tindex)
  disp('vbar...')
  oct_interpvar3d(np_id,ncrst,igrd_v,jgrd_v,ichildgrd_v,jchildgrd_v,'vbar',mask,tindex)
  disp('u...')
  oct_interpvar4d(np_id,ncrst,igrd_u,jgrd_u,ichildgrd_u,jchildgrd_u,'u',mask,tindex,N)
  disp('v...')
  oct_interpvar4d(np_id,ncrst,igrd_v,jgrd_v,ichildgrd_v,jchildgrd_v,'v',mask,tindex,N)
  disp('temp...')
  oct_interpvar4d(np_id,ncrst,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,'temp',mask,tindex,N)
  disp('salt...')
  oct_interpvar4d(np_id,ncrst,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,'salt',mask,tindex,N)
  %
  if (biol==1)
    for k=1:length(namebiol)
      disp(char(namebiol(k)))
      oct_interpvar4d(np_id,ncrst,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,char(namebiol(k)),mask,tindex,N)
    end
  end
  %
  if (pisces==1)
    for k=1:length(namepisces)
      disp(['K=',num2str(k)])
      disp(char(namepisces(k)))
      oct_interpvar4d(np_id,ncrst,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,char(namepisces(k)),mask,tindex,N)
    end
  end
end
result=close(np_id);
result=close(ncrst);
%
%  Vertical corrections
%
if (vertical_correc==1)
  for tindex=1:length(thetime)
    disp([' Time index : ',num2str(tindex),' of ',num2str(length(thetime))])                     
    oct_vert_correc(child_rst,tindex,biol,pisces,namebiol,namepisces)
  end
end

%
% Make a plot
%
if ~isempty(thetime)
  disp(' ')
  disp(' Make a plot...')
  figure(1)
  oct_plot_nestclim(child_rst,child_grd,'temp',1)
else
  disp(' ')
  disp(' Warning : no restart variable to plot...')
end

return