function oct_nested_clim(child_grd,parent_clim,child_clim,...
  vertical_correc,extrapmask,biol,bioebus,pisces)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  compute the climatology 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_clim, 'NC_NOWRITE');
Vars = var(a_id);
Varnames = [ncnames(Vars)];
nvar =length(Varnames);
%
namebiol={''};
unitbiol={''};
namepisces={''};
unitpisces={''};
timebiol=[]; cyclebiol=[];
timepisces=[];cyclepisces=[];
tbiol=[];cbiol=[];
tpisces=[];cpisces=[];


% Check type of clim file
if biol
  if (bioebus~=1) %% for biol + npzd
    timebiol={'no3_time';'o2_time';'chla_time';'phyto_time';'zoo_time'};
    cyclebiol={'no3_cyle';'o2_cycle';'chla_cycle';'phyto_cycle';'zoo_cycle'};
    namebiol={'NO3';'O2';'CHLA';'PHYTO';'ZOO'};
    unitbiol={'mMol N m-3';'mMol O m-3';'mg C l-1';'mMol N m-3';'mMol N m-3'};
    %timebiol={'no3_time';'chla_time';'phyto_time';'zoo_time'};
    disp(['Compute Biogeochemical variables type NPZD : '])
    disp(['NChlPZD or N2ChlPZD2                     '])
    disp('==========================')
  else % for biol + bioebus)
    timebiol={'no3_time';'o2_time';'chla_time';'sphyto_time';'lphyto_time';'szoo_time';'lzoo_time';'n2o_time'};
    cyclebiol={'no3_cyle';'o2_cycle';'chla_cycle';'sphyto_cycle';'lphyto_cycle';'szoo_cycle';'lzoo_cycle';'n2o_cycle'};
    namebiol={'NO3';'O2';'CHLA';'SPHYTO';'LPHYTO';'SZOO';'LZOO';'N2O'};
    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';'mMol N m-3';'mMol N m-3'};
    disp(['Compute Biological variables for BIOEBUS : '])
    disp('==========================')
  end
elseif (pisces)
  timepisces={'no3_time';'po4_time';'si_time';'o2_time';'dic_time';'talk_time';'doc_time';'fer_time'};
  cyclepisces={'no3_cycle';'po4_cycle';'si_cycle';'o2_cycle';'dic_cycle';'talk_cycle';'doc_cycle';'fer_cycle'};
  namepisces={'NO3';'PO4';'Si';'O2';'DIC';'TALK';'DOC';'FER'};
  unitpisces={'mMol N m-3';'mMol P m-3';'mMol Si m-3';'mMol O m-3';'mMol C m-3';'mMol C m-3';'mMol C m-3';'uMol Fe m-3'};
  disp('Compute Biogeochemical variables for PISCES')
  disp('=========================')
end

% Title
%
title=['Climatology file for the embedded grid :',child_clim,...
  ' using oct_parent forcing file: ',parent_clim];
disp(' ')
disp(title)
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
  error(['Both Biol NPZD and Pisces are ON, no possible yet...!'])
end
%
% 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 climatology file
%
disp(' ')
disp(' Read in the oct_parent climatology file...')
ncid = netcdf.open(parent_clim, 'NC_NOWRITE');
theta_s = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'theta_s'));
theta_b = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'theta_b'));
Tcline = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'Tcline'));
N = netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 's_rho'));
vtransform=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'Vtransform'));
hc = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'hc'));
disp([' Use oct_parent VTRANSFORM = ',num2str(vtransform)])
if ~exist('vtransform') | isempty(vtransform)
  disp([' No VTRANSFORM parameter found'])
  disp([' Use the default one VTRANSFORM = 1'])
  vtransform=1;
end
ttime = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'tclm_time'));
tcycle = netcdf.getAtt(ncid, netcdf.inqVarID(ncid, 'tclm_time'), 'cycle_length');
stime = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'sclm_time'));
scycle = netcdf.getAtt(ncid, netcdf.inqVarID(ncid, 'sclm_time'), 'cycle_length');
utime = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ssh_time'));
ucycle = netcdf.getAtt(ncid, netcdf.inqVarID(ncid, 'ssh_time'), 'cycle_length');
vtime = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ssh_time'));
vcycle = netcdf.getAtt(ncid, netcdf.inqVarID(ncid, 'ssh_time'), 'cycle_length');
sshtime = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ssh_time'));
sshcycle = netcdf.getAtt(ncid, netcdf.inqVarID(ncid, 'ssh_time'), 'cycle_length');
%
if biol==1
  for k=1:length(namebiol)
    if strcmp(char(namebiol{k}), 'NO3') || strcmp(char(namebiol{k}), 'O2') || strcmp(char(namebiol{k}), 'N2O')
      skip = 3;  % for no3_time + o2_time + n2o_time % a_id bit sloppy
                % wa tale value every 3 month to have kinf of seasonal ...
    else
      skip = 1;
    end

    eval([char(timebiol(k)),'=netcdf.getVar(ncid, netcdf.inqVarID(ncid, ''',char(timebiol(k)),'''));']);
    eval([char(cyclebiol(k)),'=netcdf.getAtt(ncid, netcdf.inqVarID(ncid, ''',char(timebiol(k)),'''), 'cycle_length');']);

    eval(['tbiol(k,:)=',char(timebiol(k)),';']);
    eval(['cbiol(k,:)=',char(cyclebiol(k)),';']);
  end
end
%
if pisces==1
  for k=1:length(namepisces)
    eval([char(timepisces(k)),'=netcdf.getVar(ncid, netcdf.inqVarID(ncid, ''',char(timepisces(k)),'''));']);
    eval([char(cyclepisces(k)),'=netcdf.getAtt(ncid, netcdf.inqVarID(ncid, ''',char(timepisces(k)),'''), 'cycle_length');']);
    %
    eval(['tpisces(k,:)=',char(timepisces(k)),';']);
    eval(['cpisces(k,:)=',char(cyclepisces(k)),';']);
  end
end
%
result=close(ncid);

climtime=ttime;
if stime~=climtime | utime~=climtime | vtime~=climtime | sshtime~=climtime
  error('Nested_clim: different times... ')
end
%
% Create the climatology file
%
disp(' ')
disp(' Create the climatology file...')
%
ncclim=oct_create_nestedclim(child_clim,child_grd,parent_clim,title,...
  theta_s,theta_b,Tcline,N,...
  ttime,stime,utime,vtime,sshtime,...
  tcycle,scycle,ucycle,vcycle,sshcycle,...
  tbiol,cbiol,tpisces, cpisces, ...
  'clobber',...
  biol,pisces,...
  timebiol, cyclebiol,timepisces,cyclepisces,...
  namebiol,namepisces,unitbiol,unitpisces,hc,vtransform);
%
% 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));
%
% 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_clim, 'NC_NOWRITE');
disp('u...')
for tindex=1:length(climtime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
  oct_interpvar4d(np_id,ncclim,igrd_u,jgrd_u,ichildgrd_u,jchildgrd_u,'u',mask,tindex,N)
end
disp('v...')
for tindex=1:length(climtime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
  oct_interpvar4d(np_id,ncclim,igrd_v,jgrd_v,ichildgrd_v,jchildgrd_v,'v',mask,tindex,N)
end
disp('zeta...')
for tindex=1:length(climtime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
  oct_interpvar3d(np_id,ncclim,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,'SSH',mask,tindex)
end
disp('ubar...')
for tindex=1:length(climtime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
  oct_interpvar3d(np_id,ncclim,igrd_u,jgrd_u,ichildgrd_u,jchildgrd_u,'ubar',mask,tindex)
end
disp('vbar...')
for tindex=1:length(climtime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
  oct_interpvar3d(np_id,ncclim,igrd_v,jgrd_v,ichildgrd_v,jchildgrd_v,'vbar',mask,tindex)
end
disp('temp...')
for tindex=1:length(climtime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
  oct_interpvar4d(np_id,ncclim,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,'temp',mask,tindex,N)
end
disp('salt...')
for tindex=1:length(climtime)
  disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
  oct_interpvar4d(np_id,ncclim,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,'salt',mask,tindex,N)
end
%
%%
%
if biol
  for k=1:length(namebiol)
    disp(char(namebiol(k)))
    for tindex=1:length(tbiol(k,:))
      disp([' Time index : ',num2str(tindex),' of ',num2str(length(tbiol(k,:)))])
      oct_interpvar4d(np_id,ncclim,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,char(namebiol(k)),mask,tindex,N)
    end
  end
end
%
if pisces
  for k=1:length(namepisces)
    disp(char(namepisces(k)))
    for tindex=1:length(tpisces(k,:))
      disp([' Time index : ',num2str(tindex),' of ',num2str(length(tpisces(k,:)))])
      oct_interpvar4d(np_id,ncclim,igrd_r,jgrd_r,ichildgrd_r,jchildgrd_r,char(namepisces(k)),mask,tindex,N)
    end
  end
end
netcdf.close(np_id);
close(ncclim);
%
%  Vertical corrections
%
if (vertical_correc==1)
  disp('Process variable physical variables')
  for tindex=1:length(climtime)
    disp([' Time index : ',num2str(tindex),' of ',num2str(length(climtime))])
    oct_vert_correc(child_clim,tindex,0,0,namebiol,namepisces)
  end
  %%
  if biol
    disp('Process variable NPZD')
    for k=1:length(namebiol)
      disp(char(namebiol(k)))
      for tindex=1:length(tbiol(k,:))
        disp([' Time index : ',num2str(tindex),' of ',num2str(length(tbiol(k,:)))])
        oct_vert_correc_onefield(child_clim,tindex,char(namebiol(k)))
      end
    end
  end
  %%
  if pisces
    disp('Process variable PISCES')
    for k=1:length(namepisces)
      disp(char(namepisces(k)))
      for tindex=1:length(tpisces(k,:))
        disp([' Time index : ',num2str(tindex),' of ',num2str(length(tpisces(k,:)))])
        oct_vert_correc_onefield(child_clim,tindex,char(namepisces(k)))
      end
    end
  end
end
%
% Make a plot
%
disp(' ')
disp(' Make a plot...')
figure(1)
oct_plot_nestclim(child_clim,child_grd,'temp',4)
return