#if defined key_agrif
  SUBROUTINE Agrif_InitWorkspace
!
!     Modules used:
!
   USE parameters
   USE comvarp3d
   USE comvarp2d
   USE comvars2d
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,     ONLY : MUSTANG_alloc,cv_sed,l_repsed,l_morphocoupl
   USE MUSTANG_MARS,   ONLY : sed_init_fromfile,sed_init
   USE sedimorpho,     ONLY : sed_morphoinit,sed_alloc_morpho
#endif
!
!     Declarations:
!
   IMPLICIT NONE
!
!     Variables
!
!     Begin
!
!
   IF ( .NOT. Agrif_Root() ) THEN
     imax=imaxm2+2
     jmax=jmaxm2+2
#if ! defined key_MPI_2D
     ljmin=jmin
     ljmax=jmax
     ljminp1=jmin+1
     ljminp2=jmin+2
     ljmaxm1=jmax-1
     ljminm2=jmin
     ljmaxp2=jmax
     ljminm1=jmin
     ljmaxp1=jmax
     ljmaxm1p1=jmax-1
     ljminp1m1=jmin+1
     ljminp1m2=jmin+2
     ljmaxp3=jmax
!
! pour les i
     limin=imin
     limax=imax
     liminp1=imin+1
     liminp2=imin+2
     limaxm1=imax-1
     liminm2=imin
     limaxp2=imax
     liminm1=imin
     limaxp1=imax
     limaxm1p1=imax-1
     liminp1m1=imin+1
     liminp1m2=imin+2
     limaxp3=imax
#endif

! no shift of the obc
     ishift_obc_west=0
     ishift_obc_east=0
     jshift_obc_south=0
     jshift_obc_north=0

   END IF

   RETURN

  END SUBROUTINE Agrif_InitWorkspace

   !!======================================================================

  SUBROUTINE Agrif_InitValues
!     ------------------------------------------------------------------
!     You should declare the variable which has to be interpolated here
!     -----------------------------------------------------------------
!
#include "toolcpp.h"

!     Modules used:
!
   USE comvarp2d
   USE comvarp3d
   USE comvars2d
   USE comvars3d
   USE comdem,       ONLY : nsdem
   USE comcolig2d,   ONLY : igdv_old
   USE comcolig3d,   ONLY : intz_exp_vstar
   USE comfluxdo
   USE comturb,      ONLY : eps
   USE agrifcommonid
   USE comdiag,      ONLY : l_diag
   USE output, ONLY : output_init,output_mng,l_out_demer,nb_out
   USE casinit,     ONLY : cas_init,cas_init_misc
#ifdef key_nudging
   USE nudging,     ONLY : nud_init,l_nud_ssh,l_nud_sal,l_nud_temp,l_nud_boz,l_nud_u3d,l_nud_ubt,l_nud_zonalmean
#endif
#ifdef key_ibm
   USE ibm,         ONLY : ibm_init
#endif
#ifdef key_diagdyn
   USE diagdyn,     ONLY : diagdyn_init
#endif
#ifdef key_substance
   USE subreaddat,   ONLY : sub_read_var,sub_read_vardiag
   USE subreadflxatm, ONLY : sub_sflx_read
   USE comsubstance, ONLY : nv_dri,cv_wat,l_cvrain_readfile,l_subflxatm_readfile,l_bilan
   USE usersubstance,   ONLY : user_init_vdriv
   USE submassbalance,  ONLY : sub_budget_readdomain
#ifdef key_biolo
   USE bioloinit,  ONLY : biolo_userinit
#endif
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,     ONLY : MUSTANG_alloc,cv_sed,l_repsed,l_morphocoupl
   USE MUSTANG_MARS,   ONLY : sed_init_fromfile,sed_init
   USE sedimorpho,     ONLY : sed_morphoinit
#endif
#endif
   USE sflxatm, ONLY : drag,exch,imeteo_dragtype
   USE sflxturb, ONLY : sflx_coef_init
   USE trackpoints, ONLY : l_points,track_points
#ifdef key_trajec3d
   USE traject3d, ONLY : traj_init3d
#endif
#ifdef key_trajec2d
   USE traject2d, ONLY : traj_init2d
#endif
   USE_MPI toolmpi, ONLY : MPI_START1,MPI_START2,MPI_START3,MPI_END1
   USE_MPI toolmpi, ONLY : ADD_TAG,MASTER
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh
!
!
!     Declarations:
!
   Implicit none
!
!     Variables
!
   INTEGER :: i,j,nb,narg
!
!
!     Begin
!
!     Specific fine grid Initializations
!
!
!     We should declare the type of variable which have to be interpolated
!
   integer :: k,l
   External :: interpssh,interp3dtraceur

!
   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter Agrif_InitValues routine'
!
! Variables initialization
!
   nracin='zzzz'
   root='./'
   nbouc=0
   nbouc_macho=0
!
! Call routines for initializations
!
#ifdef key_MPI_2D
   CALL MPI_START1
#endif

   IF (l_testcase) THEN

   ! Initialize test case configuration
     CALL cas_init
     CALL cas_init_misc

   ! Initialize coefficients of turbulence
     IF (turb_nbeq==1 .OR. turb_nbeq==2) CALL init_turbcoef(ect,psi)

     t=Agrif_Parent(t)
     dt=Agrif_Parent(dt)/(Agrif_rhot())
     halfdt=dt*0.5_rlg

   ! Initialize the prognostic variables
   ! - from a restart file ( restart file from MARS or climatological file )
   ! - from a spatial interpolation of mother fields
     CALL Agrif_Init_mainVars

   ELSE   ! general cas (realistic configuration)

!cval utile ? car reinitialise vi paraspec1.txt  ????
     dt=Agrif_Parent(dt)/(Agrif_rhot())
     IF (Agrif_Grandmother()) THEN
       dt = Agrif_Parent(dt) * Agrif_irhot() ! For the grandmother grid, Agrif_Parent is Agrif_Root
       ! A multiplication is required (instead of a division)
     END IF
     PRINT_DBG*,'dt = ',Agrif_Parent(dt),agrif_rhot(),dt
     halfdt=dt*0.5_rlg
!cval si dt reinitialise par paraspec1.txt , haldt n est plu s correct car pas reinitialise !!!!

     narg=0
     CALL init(narg)
     tsauv=Agrif_parent(tsauv)  ! cval pourquoi ? car zoom peut demarrer a une date ulterieure au demarrage du rang precedent ?
     l_obc2drank_save=.false.   ! cval mouais pour pallier erreur user

   ! Read parameters of XXXXX
#ifdef key_waveforcing
     CALL lec_para_swan
#endif

   ! Initialize flux exchange coefficients
     IF (imeteo_exchtype==1) CALL sflx_coef_init(drag,exch,imeteo_dragtype)

   ! Read the bathymetry
     CALL bathy

   ! Initialize ig,id,jb,jh
     CALL grid_mask

   ! Divide the domain for MPI run
#ifdef key_MPI_2D
     CALL MPI_START2(1,1)
     CALL MPI_START3
#endif

   ! Connection of the bathymetries
   ! ------------------------------
     CALL Agrif_Init_Connect_Bathy

   ! COMMENT
   ! ------------------------------
     CALL Agrif_Declare_Variables

   ! Allocation of variables
   ! -----------------------
   ! Dynamical variables
     CALL ALLOC
   ! Intermediate variables
     CALL alloc_inter
   ! Stationnary variable (grid metrics)
     CALL alloc_stat

   ! Initialize metrics (dx,dy,sig,lon,lat...)
     CALL grid_metric

#ifdef key_substance
   ! Read file variables.dat describing substances
   ! Allocate substance variables
   ! Make variables in order
     CALL sub_read_var

#ifdef key_sedim_MUSTANG
!toseelater utile cela ou dramatique ? (cval)
     CALL MUSTANG_alloc
     IF (l_repsed) THEN
       CALL sed_init_fromfile
     ELSE
       CALL sed_init
     END IF
     !!CALL sed_settlveloc_ini
     IF (l_morphocoupl)CALL sed_morphoinit
#endif

#ifdef key_biolo
   ! Read file vardiag.dat describing diagnoses
   ! Allocate diagnoses variables
     CALL sub_read_vardiag
#endif
#ifdef key_eco3m
   ! Initialization of ECO3M
     IF (nbouc == 0) THEN
       nbcallbio = -1
       CALL main_bio
       CALL alloc_var_eco3m
       !In order to compute the modulo of the call to ECO3M
       dt_phy_ini= dt
     END IF
#endif
#endif


     dt=Agrif_Parent(dt)/(Agrif_rhot())
     t=Agrif_Parent(t)
     tdebagrif_messat=Agrif_Parent(tdebagrif_messat)
!cval ne devrait on pas remettre halfdt= ici ?
!
!********************************

   ! Initialize viscosity coefficients
     CALL init_viscoef
   ! Initialize coefficients of turbulence
     CALL init_turbcoef(ect,psi)
   ! Initialize sal, temp, cv_wat
     CALL init_tssub(sal,temp   &
#ifdef key_substance
                        ,cv_wat &
#endif
                             )


   ! Initialize the prognostic variables
   ! - from a restart file ( restart file from MARS or climatological file )
   ! - from a spatial interpolation of mother fields
     CALL Agrif_Init_mainVars


   ! Initialize wave fields
#ifdef key_waveforcing
     CALL init_wave_forcing_fric
     CALL lect_wave
!cval ?     CALL bathy
     CALL compute_wave_forcing
#endif

   ! Initialize ecosystem variables
#ifdef key_substance
#ifdef key_biolo
   ! User special initialization of cv_wat.
     CALL biolo_userinit(cv_wat)
#endif
!defined key_contaminant
#if defined key_contaminant
   !!!!! deja fait avec sedim??     IF (.NOT.l_initfromfile) CALL
!sed_settlveloc_ini
   CALL conta_read_react
   !!!!!     CALL conta_userinit
#endif
#endif

#ifdef key_substance
   ! initialization of driving variables
     IF (nv_dri > 0) &
       CALL user_init_vdriv(cv_wat                       &
#ifdef key_sedim_MUSTANG
                           ,cv_sed                       &
#endif
                           )
#endif

   ! Read output.dat and allocate variables with the correct number of output files
     CALL output_init

   END IF   ! ends test on l_testcase

#if defined key_chainop
   ! Read output.dat and allocate variables with the correct number of output files
!mettre ou pas CALL output_init
#endif

#ifdef key_nudging
   CALL nud_init(l_nud_ssh,l_nud_sal,l_nud_temp,l_nud_boz,l_nud_u3d,l_nud_ubt,l_nud_zonalmean)
#endif

#ifdef key_substance
   !!! introduction of submassbalance
   IF (l_bilan) CALL sub_budget_readdomain(ssh)
   PRINT_DBG*, 'call read_limits'
   !!! end submassbalance
#endif

   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2)) ! for dyn2d_dt val
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2)) ! for dyn2d_dt val

   ! Save summary of simulation
   IF_MPI (MASTER) THEN
     CALL summary
   ENDIF_MPI

!cval on ne fait rien pour les listing files ?

   ! Initialize meteorological variables if l_meteo_stat or l_meteo_hom
   IF (.NOT. l_testcase) call meteo

   ! Read file river.dat, initialize river discharges inputs with time interpolation
   CALL river(ssh)

#ifdef key_substance
   ! Read file outflow.dat, file describing discharges
   ! initialize discharges inputs with time interpolation
   CALL outflow

   ! Read files sflx_cvrain and/or sflx_subflxatm
   ! if temporal variations of cvrain or of atm. fluxes for at least one
   ! substance
   !  (l_cvrain_readfile=.true. or l_subflxatm_readfile=.true.)
   IF (l_cvrain_readfile .or. l_subflxatm_readfile) CALL sub_sflx_read
#endif

   ! Estimate buoyancy from sal and temp
   IF (l_modele2d) THEN
     IF (icon == 2) CALL ts_stateeq2d(sal)
   ELSE
     IF (icon >= 1 .AND. l_stateeq_lin) CALL ts_stateeq_linear (sal,temp,boz,pn2)
     IF (icon >= 1 .AND. .NOT.l_stateeq_lin) CALL ts_stateeq_blumberg(sal,temp,boz,pn2)
   END IF


!cval j ai deplace pour faire comme dans main mais je crois que j avais ete obligee
! de mettre ces appels apres demerliac sans que j y comprenne rien

   ! Intialize variables for traking
   IF (l_points) CALL track_points(u,v,uz,vz,ssh,sal,temp,          &
#ifdef key_substance
                                   cv_wat,        &
#endif
                                   ect,boz)

   ! Initialize diagnoses variables
   IF (l_diag) CALL diagnostic(u,v,uz,vz,ssh,sal,temp,            &
#ifdef key_substance
                               cv_wat,          &
#endif
                               ect,boz)

! Must be called on because initialization for nbouc=0
! Moved because variables may be saved at initial step

   CALL output_mng

   ! Initialize demerliac filter
   DO nb=1,nb_out
!cval bizarre : meme variables
     IF (l_out_demer(nb)) CALL out_demerliac(ssh(limin:limax,ljmin:ljmax), &
                                             u(limin:limax,ljmin:ljmax),uz(:,limin:limax,ljmin:ljmax), &
                                             v(limin:limax,ljmin:ljmax),vz(:,limin:limax,ljmin:ljmax), &
                                             sal(:,limin:limax,ljmin:ljmax),temp(:,limin:limax,ljmin:ljmax), &
                                             ssh(limin:limax,ljmin:ljmax), &
                                             u(limin:limax,ljmin:ljmax),uz(:,limin:limax,ljmin:ljmax), &
                                             v(limin:limax,ljmin:ljmax),vz(:,limin:limax,ljmin:ljmax), &
                                             sal(:,limin:limax,ljmin:ljmax),temp(:,limin:limax,ljmin:ljmax), &
#ifdef key_substance
                                             cv_wat(:,:,limin:limax,ljmin:ljmax), cv_wat(:,:,limin:limax,ljmin:ljmax),    &
#endif
                                             nsdem)
   END DO


   ! Initialize diagnoses variables relative to terms of the momentum equation
   ! after output_mng to be allowed to complete the main output file ( fileresu(1) )
#ifdef key_diagdyn
   CALL diagdyn_init
#endif

!
#ifdef key_trajec3d
   IF (.not.l_modele2d) CALL traj_init3d(ssh)
#endif
#ifdef key_trajec2d
   IF (l_modele2d) CALL traj_init2d
#endif
!
#ifdef key_ibm
   IF (.NOT.l_modele2d) CALL ibm_init(ssh,sal,temp, &
#ifdef key_substance
                                      cv_wat,       &
#endif
                                      boz)
#endif

   PRINT_DBG*, 'exit Agrif_InitValues routine'
!
   RETURN

  END SUBROUTINE Agrif_InitValues

   !!======================================================================

  SUBROUTINE Agrif_Init_Connect_Bathy
   USE parameters
   USE agrifcommonid
   USE ionc4, ONLY : ionc4_openr,ionc4_read_xy,ionc4_close
   USE_MPI toolmpi, ONLY : NPROCS,MPI_COMM_MARS
!
!     Declarations:
!
   IMPLICIT NONE
!
!     Variables
!
   REAL(kind=riolg),dimension(imin:imax,jmin:jmax) :: tmp2drlg
   integer :: glob_imin, glob_jmin, glob_imax, glob_jmax
   INTEGER :: ncid, varid, ierr, i, nb
#ifdef key_MPI_2D
   integer,dimension(:,:),allocatable :: i_connect_loc
#endif
#ifdef key_MPI_2D
   integer :: ierr_mpi
   include 'mpif.h'
#endif

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'Enter Agrif_Init_Connect_Bathy routine'

! Read the connection areas
! -------------------------

   IF (nb_connection_areas .ne. 0) THEN
!cval Franck a faire vraiment modif pgarreau


     ! Read from the bathymetry file les indices de zones de recouvrement dans
     ! la grille mere
     ! Pas generique ne marche que pour le premier d enfant. Donc il ne peut y
     ! avoir recouvrement de plusieurs grilles uniquement pour le premier niveau
     IF (.not.allocated(imin_connect)) THEN
       ierr = nf90_open(file_bathy,NF90_NOWRITE,ncid)
       ierr = nf90_inq_varid(ncid,'H0_weights',varid)
       ierr = nf90_get_att(ncid,varid,'nb_connection_areas',nb_connection_areas)
       allocate(imin_connect(nb_connection_areas))
       allocate(imax_connect(nb_connection_areas))
       allocate(jmin_connect(nb_connection_areas))
       allocate(jmax_connect(nb_connection_areas))
       ierr = nf90_get_att(ncid,varid,'imin_connect',imin_connect)
       ierr = nf90_get_att(ncid,varid,'imax_connect',imax_connect)
       ierr = nf90_get_att(ncid,varid,'jmin_connect',jmin_connect)
       ierr = nf90_get_att(ncid,varid,'jmax_connect',jmax_connect)
       ierr = nf90_close(ncid)
       allocate(connected_values(nb_connection_areas))
       DO i=1,nb_connection_areas
         allocate(connected_values(i)%temp(kmax+2,imin_connect(i):imax_connect(i),jmin_connect(i):jmax_connect(i)))
         allocate(connected_values(i)%sal(kmax+2,imin_connect(i):imax_connect(i),jmin_connect(i):jmax_connect(i)))
         allocate(connected_values(i)%uz(kmax+2,imin_connect(i):imax_connect(i),jmin_connect(i):jmax_connect(i)))
         allocate(connected_values(i)%vz(kmax+2,imin_connect(i):imax_connect(i),jmin_connect(i):jmax_connect(i)))
       END DO
     END IF

   END IF  ! ends test on nb_connection_areas

!SPETTON TODO : verifie que ne fait le tout que si nb_connection_areas=1


   ! Read the weights from the bathymetric file (after connect_topo)

   allocate(update_at_Tpoints(liminm2:limaxp2,ljminm2:ljmaxp2))
   allocate(update_at_Upoints(liminm2:limaxp2,ljminm2:ljmaxp2))
   allocate(update_at_Vpoints(liminm2:limaxp2,ljminm2:ljmaxp2))

   CALL ionc4_openr(file_bathy)
   tmp2drlg(:,:) = 0.
   CALL ionc4_read_xy(file_bathy,'H0_weights',tmp2drlg,imin,imax,jmin,jmax)
   update_at_Tpoints = tmp2drlg(liminm2:limaxp2,ljminm2:ljmaxp2)
   tmp2drlg = 0.
   CALL ionc4_read_xy(file_bathy,'HX_weights',tmp2drlg,imin,imax,jmin,jmax)
   update_at_Upoints = tmp2drlg(liminm2:limaxp2,ljminm2:ljmaxp2)
   tmp2drlg = 0.
   CALL ionc4_read_xy(file_bathy,'HY_weights',tmp2drlg,imin,imax,jmin,jmax)
   update_at_Vpoints = tmp2drlg(liminm2:limaxp2,ljminm2:ljmaxp2)
   CALL ionc4_close(file_bathy)

   ! TODO comment 

   IF (nb_connection_areas .ne. 0) THEN

#ifdef key_MPI_2D

   ALLOCATE(i_connect(nb_connection_areas,4,0:NPROCS-1))
   ALLOCATE(i_connect_loc(nb_connection_areas,4))
   glob_imin = (Agrif_ix()-1)*agrif_irhox()+limin
   glob_jmin = (Agrif_iy()-1)*agrif_irhoy()+ljmin
   glob_imax = (Agrif_ix()-1)*agrif_irhox()+limax
   glob_jmax = (Agrif_iy()-1)*agrif_irhoy()+ljmax

   DO nb=1,nb_connection_areas
     i_connect_loc(nb,1)=max(glob_imin,imin_connect(nb))
     i_connect_loc(nb,2)=min(glob_imax,imax_connect(nb))
     i_connect_loc(nb,3)=max(glob_jmin,jmin_connect(nb))
     i_connect_loc(nb,4)=min(glob_jmax,jmax_connect(nb))
   END DO

   CALL MPI_ALLGATHER(i_connect_loc,4*nb_connection_areas,MPI_INTEGER,i_connect, &
                      4*nb_connection_areas,MPI_INTEGER,MPI_COMM_MARS,ierr_mpi)
   DEALLOCATE(i_connect_loc)
   ALLOCATE(has_connection_areas(nb_connection_areas,0:NPROCS-1))
   ALLOCATE(thearea_inthegrid(nb_connection_areas,0:NPROCS-1))
   thearea_inthegrid = .FALSE.

   DO i=1,nb_connection_areas
     DO nb=0,NPROCS-1
       has_connection_areas(i,nb) = (i_connect(i,1,nb)<=i_connect(i,2,nb)).AND.(i_connect(i,3,nb)<=i_connect(i,4,nb))
       IF (has_connection_areas(i,nb)) THEN
         thearea_inthegrid(i,:) = .TRUE.
       END IF
     END DO
   END DO

#else

   ALLOCATE(i_connect(nb_connection_areas,4,0:0))
   ALLOCATE(thearea_inthegrid(nb_connection_areas,0:0))
   ALLOCATE(has_connection_areas(nb_connection_areas,0:0))

   glob_imin = (Agrif_ix()-1)*agrif_irhox()+limin
   glob_jmin = (Agrif_iy()-1)*agrif_irhoy()+ljmin
   glob_imax = (Agrif_ix()-1)*agrif_irhox()+limax
   glob_jmax = (Agrif_iy()-1)*agrif_irhoy()+ljmax

   thearea_inthegrid = .FALSE.

   DO nb=1,nb_connection_areas
     i_connect(nb,1,0)=max(glob_imin,imin_connect(nb))
     i_connect(nb,2,0)=min(glob_imax,imax_connect(nb))
     i_connect(nb,3,0)=max(glob_jmin,jmin_connect(nb))
     i_connect(nb,4,0)=min(glob_jmax,jmax_connect(nb))
     has_connection_areas(nb,0)=(i_connect(nb,1,0)<=i_connect(nb,2,0)).AND.(i_connect(nb,3,0)<=i_connect(nb,4,0))
     IF (has_connection_areas(nb,0)) THEN
       thearea_inthegrid(nb,:) = .TRUE.
     END IF
   END DO

#endif

   END IF  ! ends test on nb_connection_areas

   PRINT_DBG*, 'exit Agrif_Init_Connect_Bathy routine'

   RETURN

  END SUBROUTINE Agrif_Init_Connect_Bathy

   !!======================================================================

  SUBROUTINE Agrif_Init_mainVars

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE step3dxy  ***
   !&E
   !&E ** Purpose : Intialize child grids (and update mother grid)
   !&E
   !&E ** Description :
   !&E
   !&E ** Reference :
   !&E                See MARS3D documentation (AGRIF part)
   !&E
   !&E ** History :
   !&E       !           (L. Debreu)  Original code
   !&E       !  2016-01  (S. Petton) Reorganisation of the initialization
   !&E
   !&E---------------------------------------------------------------------


   USE parameters
   USE comvars2d
   USE comvarp2d, ONLY : ssh,temp,sal,u,v
   USE comvarp3d, ONLY : uz,vz
   USE comvars3d, ONLY : nz,kz,wz
   USE agrifcommonid
   USE comfluxdo, ONLY : xeadv
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh
   USE tooldebug

   External :: updatetemp,updatesal,updatessh,updateturb
   External :: updatekz, updatenz,updateuz,updatevz,updatecv_wat, updatesurf, updatefwet
   External :: interpu_stepu,interpv_stepu,interpssh,interpxeu,interpxev,interpuz_stepu
   External :: interpvz_stepu,interp3dturb,interp3dkz, interp3dnz
   External :: interp3dtraceur,updatewz,updatehzex,updatehzey
#ifdef key_substance
   External :: interp3dsub
#endif
#if defined key_agrif_obc_sponge_ts || defined key_agrif_obc_sponge_uv
   Integer, parameter :: MaxSearchAfterInit=5
#else
   Integer, parameter :: MaxSearchAfterInit=2
#endif

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter Agrif_Init_mainVars routine'

   ! Use of gridmother cells within a 5 cells radius to interpolate (15 cells for 1st interpolation)
   CALL Agrif_Set_MaskMaxSearch(MaxSearchAfterInit)

   ! Allocation of temporary arrays used in AGRIF procname for temporary variables
   ! for the treatment (interpolation) of intermediate variables like xeuv,uv,uvz.
   ! (direct use of procname for ssh,u,v,uz,vz variables)
   ALLOCATE(nesting_tabtemp_2d(liminm2:limaxp2,ljminm2:ljmaxp2,2))
   ALLOCATE(nesting_tabtemp_3d(kmax,liminm2:limaxp2,ljminm2:ljmaxp2))

   ! Restart from a file ( restart file from MARS or climatological file )
   IF (l_initfromfile) THEN
     CALL init_fromfile(u,uz,v,vz,ssh,sal,temp,        &
                        intz_exp_vstar,igdv_old,  &
#ifdef key_substance
                        cv_wat,      &
#endif
                        ect,psi,eps)
   END IF

   ! Now there are two possibilities :
   ! 
   ! 1) Variable has not already been read above
   !  -> CALL Agrif_Init_Variable with c_agrif_interp='  _init' 
   !     + Initialisation over the whole domain : Agrif_Interp_variable
   !     + Initialisation of Boundary Condition arrays (child % oldvalues2d in AGRIF) at this time step
   !          (There is CALL Agrif_Bc_Variable(id, procname, 1.) in the subroutine Agrif_Init_Variable ; 
   !           Here, callweights are set to 1. for the first time interpolation because c_agrif_interp will change result in procname) 
   ! 
   ! 2) Variable is already set up
   !  -> CALL Agrif_Bc_Variable with c_agrif_interp='do_nothing'
   !     + Initialisation of Boundary Condition array (child % oldvalues2d in AGRIF)
   !          (calledweigths are not important because it doesn t change anything to child grid due to c_agrif_interp)


   ! We use only mothergrid values unequal to Agrif_SpecialValue to interpolate spatially
   Agrif_UseSpecialValue=.TRUE.
   Agrif_SpecialValue=valmanq

   ! Le Max Search n est pris en compte :
   ! - qu a la premiere utilisation du profil_id par AGRIF
   ! - que pour les directions ou la variable est centree
   ! - que pour les methodes  interpolation : PPM, ENO, PPM_LIM et WENO
   ! Donc ne concerne que SSH, T et S et les substances.
   ! On met Max Search a 15 pour les Agrif_Init_Variables et a 5 par defaut pour le traitement des obc


   ! SSH
   IF (.NOT. ((l_initfromfile .AND. l_init_rssh) .OR. (l_initfromfile .AND. l_init_restart)) ) THEN
     ! None restart file, interpolation of data over the full grid using Max Search set to 15
     CALL Agrif_Set_MaskMaxSearch(15)
     c_agrif_interp='ssh_init'
     nesting_tabtemp_2d(:,:,1) = ssh   ! necessaire pour garder les memes resultats sur PRONY. Non compris. 
     CALL Agrif_Init_Variable(ssh_interp_init_id,procname=interpssh)
     ssh = nesting_tabtemp_2d(:,:,1)
     CALL Agrif_Set_MaskMaxSearch(MaxSearchAfterInit)
   END IF
   c_agrif_interp='do_nothing'    ! Treatment of obc only
   CALL Agrif_Bc_Variable(ssh_interp_id,procname=interpssh)

   xeadv(limin:limax,ljmin:ljmax)=ssh(limin:limax,ljmin:ljmax)
   CALL init_height(ssh)

   ! TEMP & SAL
   IF (.NOT. ((l_initfromfile .AND. l_init_rtemp .AND. l_init_rsal) .OR. (l_initfromfile .AND. l_init_restart))) THEN
     ! None restart file, interpolation of data over the full grid using Max Search set to 15
     CALL Agrif_Set_MaskMaxSearch(15)
     c_agrif_interp='traceur_init'
     CALL Agrif_Init_Variable(traceur_init_id,procname=interp3dtraceur)
     CALL Agrif_Set_MaskMaxSearch(MaxSearchAfterInit)
   END IF
   c_agrif_interp='do_nothing'    ! Treatment of obc only
   CALL Agrif_Bc_Variable(traceur_id,procname=interp3dtraceur)

#ifdef key_substance
   ! SUBS
   IF (.NOT. l_restart_subs) THEN
     ! None restart file, interpolation of data over the full grid using Max Search set to 15
     CALL Agrif_Set_MaskMaxSearch(15)
     c_agrif_interp='cv_wat_init'
     CALL Agrif_Init_Variable(cv_wat_init_id,procname=interp3dsub)
     CALL Agrif_Set_MaskMaxSearch(MaxSearchAfterInit)
   END IF
   c_agrif_interp = 'do_nothing'    ! Treatment of obc only
   CALL Agrif_Bc_Variable(cv_wat_id,procname=interp3dsub)
#endif

   ! As velocities are set to 0 when there is no water, we could use them to interpolate spatially 
   Agrif_UseSpecialValue=.FALSE.

   ! U & V
   IF ((l_initfromfile .AND. l_init_rbtvel) .OR. (l_initfromfile .AND. l_init_restart)) THEN
     c_agrif_interp = 'do_nothing'
     CALL Agrif_Bc_Variable(u_id,procname=interpu_stepu)

     c_agrif_interp = 'do_nothing'
     CALL Agrif_Bc_Variable(v_id,procname=interpv_stepu)
   ELSE
     c_agrif_interp='u_init'
     CALL Agrif_Init_Variable(u_id,procname=interpu_stepu)
     u(:,:) = nesting_tabtemp_2d(:,:,1)
     ! Despite the profile, these MPI exchanges are required due to the loop indexes used in interpv_stepu
     CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
     CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))

     c_agrif_interp='v_init'
     CALL Agrif_Init_Variable(v_id,procname=interpv_stepu)
     v(:,:) = nesting_tabtemp_2d(:,:,1)
     ! Despite the profile, these MPI exchanges are required due to the loop indexes used in interpv_stepu
     CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))
     CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))
   END IF

   ! UZ & VZ
   IF ((l_initfromfile .AND. l_init_r3dvel) .OR. (l_initfromfile .AND. l_init_restart)) THEN
     c_agrif_interp = 'do_nothing'
     CALL Agrif_Bc_Variable(uz_id,procname=interpuz_stepu)

     c_agrif_interp = 'do_nothing'
     CALL Agrif_Bc_Variable(vz_id,procname=interpvz_stepu)
   ELSE
     c_agrif_interp='uz_init'
     CALL Agrif_Init_Variable(uz_id,procname=interpuz_stepu)
     uz(:,:,:) = nesting_tabtemp_3d(:,:,:)
     ! Despite the profile, these MPI exchanges are required due to the loop indexes used in interpuz_stepu
     CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
     CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))

     c_agrif_interp='vz_init'
     CALL Agrif_Init_Variable(vz_id,procname=interpvz_stepu)
     vz(:,:,:) = nesting_tabtemp_3d(:,:,:)
     ! Despite the profile, these MPI exchanges are required due to the loop indexes used in interpvz_stepu
     CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
     CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   END IF

#if ! defined key_agrif_1way
   ! Update of all variables in case of :
   !    - coarsening
   !    - restart configuration
   !    - start configuration (not important except for update surf)

   ! We fill daughter cells if they differ from Agrif_SpecialValueFineGrid to update spatially
   Agrif_UseSpecialValueInUpdate=.TRUE.
   Agrif_SpecialValueFineGrid = valmanq

   first_agrif_update = .TRUE.
   IF (l_init_restart .AND. l_initfromfile) first_agrif_update = .FALSE.

   IF (.NOT. Agrif_Root() .AND. icon >= 1) THEN
     CALL Agrif_Update_Variable(sal_id ,procname = updatesal)
     CALL Agrif_Update_Variable(temp_id,procname = updatetemp)
#ifdef key_substance
     CALL Agrif_Update_Variable(cv_wat_id,procname = updatecv_wat)
#endif
   END IF

   IF (.NOT. Agrif_Root()) THEN
     CALL Agrif_Update_Variable(ssh_id,procname = updatessh)
   END IF

   IF (Agrif_Root()) THEN
     CALL Agrif_Update_Variable(kz_id,procname = updatekz)
   END IF

   ! As velocities are set to 0 when there is no water, we could use them to update spatially 
   Agrif_UseSpecialValueInUpdate=.FALSE.

   ! Update surf for land mask only at first time step only
   CALL Agrif_Update_Variable(surf_id, procname = updatesurf)

   ! Update fwet
   IF (fwetmin /= 1.0_rsh) THEN
     CALL Agrif_Update_Variable(fwet_id, procname = updatefwet)
   END IF 

   IF (l_init_restart .AND. l_initfromfile) THEN

     IF (Agrif_Root()) THEN
       CALL Agrif_Update_Variable(hzex_id, locupdate2=(/-1,-2/) , procname = updatehzex)
       CALL Agrif_Update_Variable(hzey_id, locupdate3=(/-1,-2/) , procname = updatehzey)
     END IF
     CALL Agrif_Update_Variable(uz_update_id, locupdate2=(/-1,-2/), procname = updateuz)
     CALL Agrif_Update_Variable(vz_update_id, locupdate3=(/-1,-2/), procname = updatevz)

#if defined key_MPI_2D
     IF (.NOT. Agrif_GrandMother()) THEN
       CALL Agrif_ChildGrid_To_ParentGrid()
       CALL exchange_uv() ! MPI exchanges of uz,vz,u,v
       CALL Agrif_ParentGrid_To_ChildGrid()
     END IF
#endif

     IF (Agrif_Root()) THEN
       ! Coarsening : Update SSH and estimate Wz everywhere over grandmother domain 
       CALL Agrif_ChildGrid_To_ParentGrid()
       CALL update_wz_coarse()
       CALL Agrif_ParentGrid_To_ChildGrid()
     ELSE
       ! TODO : for this part, there must be a better way without using procname
       ! Estimate Wz in parent grid over the child domain with updated fluxes (fxeau & fyeau)
       CALL Agrif_Update_Variable(wz_id ,procname = updatewz)
#if defined key_MPI_2D
       CALL Agrif_ChildGrid_To_ParentGrid()
       CALL exchange_wz() ! MPI exchanges of wz
       CALL Agrif_ParentGrid_To_ChildGrid()
#endif
       ! TODO : add "Correct Boz (& Turb ?)" in parent grid over the child domain with updated temperature & salinity
       ! CALL Agrif_ChildGrid_To_ParentGrid()
       ! CALL Correct_density_and_turb()
       ! CALL Agrif_ParentGrid_To_ChildGrid()
     END IF

   END IF

   first_agrif_update = .FALSE.
#endif /* key_agrif_1way */

   ! On met a True quand besoin
   Agrif_UseSpecialValue=.FALSE.
   Agrif_UseSpecialValueInUpdate=.FALSE.

   PRINT_DBG*, 'exit Agrif_Init_mainVars routine'

   RETURN

  END SUBROUTINE Agrif_Init_mainVars

   !!======================================================================

  SUBROUTINE Compute_Heights
   USE parameters
   USE comvars2d
   USE agrifcommonid
   INTEGER :: test
   REAL(KIND=rsh),DIMENSION(limin:limax)                     :: tfh,tfb
   REAL(KIND=rsh),DIMENSION(ljmin:ljmax)                     :: tfg,tfd

   RETURN

   test = agrif_indic
   CALL dyn3d_height_mask(ssh,tfg,tfd,tfb,tfh,test)

  END SUBROUTINE compute_heights

   !!======================================================================

#if defined key_MPI_2D
  SUBROUTINE exchange_uv()
   USE parameters
   use comvarp3d
   use comvars2d
   use comvarp2d
   use toolmpi, ONLY : ex_i_rsh,ex_j_rsh

OMPMPI barrier
OMPMPI master

   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))

   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))

   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))

!cval pourquoi ces echanges n etaient ils pas presents ? Oubli ou volontoire ? Si volontaire pourquoi ?
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))

OMPMPI end master
OMPMPI barrier
OMPMPI flush(uz,vz,u,v)

  END SUBROUTINE exchange_uv
#endif

   !!======================================================================

#if defined key_MPI_2D
  SUBROUTINE exchange_wz()
   USE parameters
   USE comvars3d, ONLY : wz
   USE toolmpi,   ONLY : ex_i_rsh,ex_j_rsh

OMPMPI barrier
OMPMPI master

   CALL_MPI ex_i_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL_MPI ex_j_rsh(-1,1,kmax+1,liminm1,limaxp1,ljminm1,ljmaxp1,wz(:,liminm1:limaxp1,ljminm1:ljmaxp1))

OMPMPI end master
OMPMPI barrier
OMPMPI flush(wz)

  END SUBROUTINE exchange_wz
#endif

   !!======================================================================


  SUBROUTINE Correct_density_and_turb
   USE comvars2d
   USE comvars3d
   USE comvarp2d
   USE comvarp3d

! Estimate buoyancy from sal and temp
   IF (l_modele2d) THEN
     IF (icon == 2) CALL ts_stateeq2d(sal)
   ELSE
     IF (icon >= 1 .AND. l_stateeq_lin) CALL ts_stateeq_linear (sal,temp,boz,pn2)
     IF (icon >= 1 .AND. .NOT.l_stateeq_lin) CALL ts_stateeq_blumberg(sal,temp,boz,pn2)
   END IF

   IF (turb_nbeq == 2) THEN
     CALL turblepskznz(ssh,pn2,psi,ect)
   END IF

  END SUBROUTINE Correct_density_and_turb

   !!======================================================================

  SUBROUTINE Correct_h0_and_derived()
#include "toolcpp.h"
   USE comvarp2d
   USE comvarp3d
   USE comvars2d

#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hc_sig,hcx_sig,hcy_sig,hc_sig_g,hcx_sig_g,hcy_sig_g, &
                         dcwsds,hzdcusds,hzdcwsds,hzdcusds_x,hzdcwsds_x,hzdcusds_y,hzdcwsds_y
#endif
   USE_MPI toolmpi, ONLY : exchange_rsh
   integer :: i,j,k
   REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: hx_t,hy_t

   PRINT_DBG*, 'enter Correct_h0_and_derived'

   IF (.not.Agrif_Root()) THEN
     DO j=MAX(ljmin,imin+1),ljmax
       DO i=MAX(limin,imin+1),limax
         h0(i,j)=MAX(hx(i,j),hx(i-1,j),hy(i,j),hy(i,j-1))+h0fond
       END DO
     END DO
   END IF
   IF_MPI (limin==imin) THEN
   DO j=ljminp1,ljmax
      h0(imin,j)=MAX(hx(imin,j),hy(imin,j),hy(imin,j-1))+h0fond
   END DO
   ENDIF_MPI
   IF_MPI (ljmin==jmin) THEN
   DO i=liminp1,limax
      h0(i,jmin)=MAX(hx(i,jmin),hx(i-1,jmin),hy(i,jmin))+h0fond
   END DO
   ENDIF_MPI
   IF_MPI (limin==imin .AND. ljmin==jmin) THEN
   h0(imin,jmin)=MAX(hx(imin,jmin),hy(imin,jmin))+h0fond 
   ENDIF_MPI

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,h0(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,h0(liminm2:limaxp2,ljminm2:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(h0)
   WHERE (h0(:,:) < -valmanq+20.0_rsh) h0(:,:) = -valmanq


   IF (.not.Agrif_Root()) THEN
     DO j=ljmin,MIN(ljmax,jmax-1)
       DO i=limin,MIN(limax,imax-1)
         IF (hx(i,j) >= hrdef) hx(i,j)=(h0(i,j)+h0(i+1,j))/2.0_rsh
         IF (hy(i,j) >= hrdef) hy(i,j)=(h0(i,j)+h0(i,j+1))/2.0_rsh
       END DO
     END DO
   END IF

!cval pas indispensable
!OMPMPI barrier
!OMPMPI master
!   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hx(liminm2:limaxp2,ljminm2:ljmaxp2))
!   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hx(liminm2:limaxp2,ljminm2:ljmaxp2))
!   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hy(liminm2:limaxp2,ljminm2:ljmaxp2))
!   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hy(liminm2:limaxp2,ljminm2:ljmaxp2))
!OMPMPI end master
!OMPMPI barrier
!OMPMPI flush(hx,hy)

   hx_g(limin:limax,ljmin:ljmax) = hx(limin:limax,ljmin:ljmax)
   hy_g(limin:limax,ljmin:ljmax) = hy(limin:limax,ljmin:ljmax)
   h0_g(limin:limax,ljmin:ljmax) = h0(limin:limax,ljmin:ljmax)

   ! All the cpus need to know the content of the whole array
   OMPMPI barrier
   OMPMPI master
   CALL_MPI exchange_rsh(hx_g,1)
   CALL_MPI exchange_rsh(hy_g,1)
   CALL_MPI exchange_rsh(h0_g,1)
   OMPMPI end master
   OMPMPI barrier
   OMPMPI flush(hx_g)
   OMPMPI flush(hy_g)
   OMPMPI flush(h0_g)

   hx_t(:,:)=100000.0_rsh
   hy_t(:,:)=100000.0_rsh
   WHERE (hx_g(:,:) /= -valmanq -h0fond) hx_t(:,:)=hx_g(:,:)
   WHERE (hy_g(:,:) /= -valmanq -h0fond) hy_t(:,:)=hy_g(:,:)
   dhxy_g(:,:)=1.0_rsh
   DO j=jmin+1,jmax
     DO i=imin+1,imax
       dhxy_g(i,j)=MAX(hx_g(i,j),hx_g(i-1,j),hy_g(i,j),hy_g(i,j-1))- &
       MIN(hx_t(i,j),hx_t(i-1,j),hy_t(i,j),hy_t(i,j-1))+epsilon
     END DO
   END DO

   dhx_g(:,:)=0.0_rsh
   dhy_g(:,:)=0.0_rsh
   WHERE (hx_g(:,:) < hrdef) dhx_g(:,:)=dhj
   WHERE (hy_g(:,:) < hrdef) dhy_g(:,:)=dhj

#if defined key_siggen || defined key_gencoord
   hc_sig_g(:,:)=-valmanq
   hcx_sig_g(:,:)=-valmanq
   hcy_sig_g(:,:)=-valmanq
   IF (l_bathy_meanlev) THEN
     WHERE (h0_g(:,:)/=-valmanq) hc_sig_g(:,:)=MIN(h0_g(:,:),MAX(hnivmoy_g(:,:),hc))
     WHERE (hx_g(:,:)/=-valmanq-h0fond) hcx_sig_g(:,:)=MIN(hx_g(:,:),MAX(hnivmoy_g(:,:),hc))
     WHERE (hy_g(:,:)/=-valmanq-h0fond) hcy_sig_g(:,:)=MIN(hy_g(:,:),MAX(hnivmoy_g(:,:),hc))
   ELSE
     WHERE (h0_g(:,:)/=-valmanq) hc_sig_g(:,:)=MIN(h0_g(:,:),hc)
     WHERE (hx_g(:,:)/=-valmanq-h0fond) hcx_sig_g(:,:)=MIN(hx_g(:,:),hc)
     WHERE (hy_g(:,:)/=-valmanq-h0fond) hcy_sig_g(:,:)=MIN(hy_g(:,:),hc)
   END IF
   hc_sig(liminm2:limaxp2,ljminm2:ljmaxp2)  = hc_sig_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hcx_sig(liminm1:limaxp1,ljminm1:ljmaxp1) = hcx_sig_g(liminm1:limaxp1,ljminm1:ljmaxp1)
   hcy_sig(liminm1:limaxp1,ljminm1:ljmaxp1) = hcy_sig_g(liminm1:limaxp1,ljminm1:ljmaxp1)
#endif

   h0(:,:) = h0_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hx(:,:) = hx_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hy(:,:) = hy_g(liminm2:limaxp2,ljminm2:ljmaxp2)

   dhx(:,:) = dhx_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhy(:,:) = dhy_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhxy(:,:) = dhxy_g(limin:limax,ljmin:ljmax)

   CALL grid_mask

#if defined key_siggen || defined key_gencoord
   DO j=ljminm1,ljmaxp1
     DO i=liminm1,limaxp1
       DO k=1,kmax
         hzdcusds(k,i,j)=hc_sig(i,j)+(h0(i,j)-hc_sig(i,j))*dcusds(k)
         hzdcusds_x(k,i,j)=hcx_sig(i,j)+(hx(i,j)-hcx_sig(i,j))*dcusds(k)
         hzdcusds_y(k,i,j)=hcy_sig(i,j)+(hy(i,j)-hcy_sig(i,j))*dcusds(k)
       END DO
     END DO
   END DO
   DO j=ljmin,ljmaxp1
     DO i=limin,limaxp1
       DO k=1,kmax
         hzdcwsds(k,i,j)=hc_sig(i,j)+(h0(i,j)-hc_sig(i,j))*dcwsds(k)
         hzdcwsds_x(k,i,j)=hcx_sig(i,j)+(hx(i,j)-hcx_sig(i,j))*dcwsds(k)
         hzdcwsds_y(k,i,j)=hcy_sig(i,j)+(hy(i,j)-hcy_sig(i,j))*dcwsds(k)
       END DO
     END DO
   END DO
#endif

   PRINT_DBG*, 'exit Correct_h0_and_derived'

  END SUBROUTINE Correct_h0_and_derived

   !!======================================================================

  SUBROUTINE Agrif_detect(g,size)
!
!     Modules used:
!
   USE Agrif_Util
!
!
!     Declarations:
!
   IMPLICIT NONE
!
!     Variables
!
   Type(Agrif_Grid),pointer   :: g        ! Pointer on the current grid
   INTEGER, DIMENSION(3)      :: size
!
!     Begin
!
!

!
   RETURN

  END SUBROUTINE Agrif_detect

   !!======================================================================

  SUBROUTINE  agrif_before_regridding()
  END SUBROUTINE agrif_before_regridding

   !!======================================================================

#ifdef key_MPI_2D
  SUBROUTINE Agrif_Invloc(indloc,proc,dir,indglob)
   implicit none
   integer indloc, proc, dir, indglob
!# include "param.h"
   IF (dir == 1) THEN
!           indglob = indloc + iminmpi-1
     indglob = indloc
   ELSE IF (dir == 2) THEN
!           indglob = indloc + jminmpi-1
     indglob = indloc
   ELSE
     indglob = indloc
   END IF

   RETURN

   END SUBROUTINE Agrif_Invloc

   !!======================================================================

  SUBROUTINE Agrif_get_proc_info ( i1, i2, j1, j2 )
   use parameters
   integer, intent(out) :: i1,i2
   integer, intent(out) :: j1,j2

   i1 = limin
   i2 = limax
   j1 = ljmin
   j2 = ljmax

  END SUBROUTINE Agrif_get_proc_info

   !!======================================================================

  SUBROUTINE agrif_estimate_parallel_cost(i1,i2,j1,j2,nbprocs,grid_cost)
   use agrif_util
   use parameters
   use comvars2d, ONLY : h0_g
   integer :: i1,i2,j1,j2,nbprocs
   real :: grid_cost

   grid_cost = COUNT(h0_g(i1:i2,j1:j2) /= -valmanq)/real(nbprocs)

   IF (agrif_procrank == 0) THEN
     print *,'COST = ',i1,i2,j1,j2,nbprocs,grid_cost
   END IF
  END SUBROUTINE agrif_estimate_parallel_cost
#endif

!
#else
  SUBROUTINE Subcalledbyagrif
   write(*,*) 'Impossible to bet here'
  END SUBROUTINE Subcalledbyagrif
#endif

   !!======================================================================

  SUBROUTINE agrif_dyn3d_height_mask(xe,igd,hexy, &
#if defined key_siggen || defined key_gencoord
   hzexy,       &
#endif
   ii1, ii2, jj1, jj2, test)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dyn2d_height_mask  ***
   !&E
   !&E ** Purpose : Evaluation of water column height, of the kind
   !&E              of cell we deal with and of the type of bound
   !&E              for each row and column
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : dyn3dxy.F90,dyn3dyx.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : hx,hy,xe,dhx,goulet,govlet,id,ig,jb,jh
   !&E
   !&E ** Modified variables : htx,hex,igdu,tfg,tfd,mn,ms,me,mo,
   !&E                         hty,hey,igdv,tfb,tfh,whx,why,
   !&E                         hzutx,hzwtx,hzuty,hzwty,hzex,hzey
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-07  (T. Duhaut)  Original code
   !&E       !  2008-07  (V. Garnier)  Introduction of generalized sigma
   !&E       !  2008-08  (T. Odaka) MPI
   !&E       !  2009-03  (V. Garnier) masks advection terms
   !&E       !  2009-04  (A. Mano) MPI2D
   !&E       !  2009-04  (F. Dumas) generalized sigma : wet & drying management
   !&E       !  2011-03  (F. Dumas)  bottom stress parametrisation + remove case 2 and 3
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2016-01  (S. Petton)  introduction of ip0,jp0 to estimate hzx/y in a symmetric way 
   !&E                                as we now define Agrif_Set_bc(ssh_interp_id,(/0,0/))
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

#if defined key_agrif

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : htx,hty,hx,hy,hey,hex,goulet,govlet,    &
                         h0fond,jb,jh,ig,id,l_obc_cycl_x,l_obc_cycl_y, &
                         dhx,dhy,hm,dhxy,h0
   USE comcolig2d,ONLY : igdu,igdv,mn,ms,me,mo,whx,why
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : dcusds,dcwsds,hzutx,hzwtx,hzuty,hzwty,hzex,hzey
   USE toolgeom,  ONLY : f_hzu_x,f_hzu_y,f_hzw_x,f_hzw_y
#endif
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh,ex_i_int,ex_j_int

   IMPLICIT NONE

   !! * Arguments
   INTEGER, INTENT(in)                                        :: test
   INTEGER, INTENT(in)                                        :: ii1, ii2, jj1, jj2
   INTEGER, DIMENSION(ii1:ii2, jj1:jj2), INTENT(out)          :: igd
   REAL(KIND=rsh), DIMENSION(ii1:ii2, jj1:jj2), INTENT(out)   :: hexy
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rsh), DIMENSION(1:kmax,ii1:ii2, jj1:jj2), INTENT(out)   :: hzexy
#endif
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xe

   !! * Local declarations
   INTEGER                                                     :: i,j,k,i1,i2,j1,j2
   REAL(KIND=rsh), DIMENSION(ii1:ii2,jj1:jj2)                 :: htxy
   INTEGER                                                    :: ip1, jp1, ip0, jp0


   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter AGRIF_DYN3D_HEIGHT_MASK routine'

   !-----------------------------------------------------------
   ! Initialize htx and hty :
   ! water depth at time n  (1st call from ligne : xev)
   ! water depth at time nprime (2nd call from ligne : xeint)
   ! water depth at time n+1/2  (1st call from colonne : xeu)
   ! water depth at time (n+1/2)prime (2nd call from colonne : xeint)
   ! htx and hty are used in momentum equations
   !-----------------------------------------------------------

!$OMP PARALLEL PRIVATE(i,j,k,ip1,jp1,ip0,jp0)
   IF (test == 2) THEN
!$OMP DO SCHEDULE(RUNTIME)
     DO j=MAX(jj1,ljminm1),MIN(jj2,ljmaxp1)
       DO i=MAX(ii1,liminm1),MIN(ii2,limaxp1)

         jp1 = j + 1
         jp0 = j
         IF (j==jmax) jp1 = j
         IF (j==jmin) jp0 = j + 1


         htxy(i,j)=hy(i,j)+0.5_rsh*(xe(i,jp0)+xe(i,jp1))

!       IF ( htxy(i,j)<=0.0 ) With such a test mn and ms are set to 1 more often as htx is smaller when negative

         IF ( htxy(i,j)<=0.0_rsh .AND. ( hy(i,j)+xe(i,jp0) > 0.0_rsh .OR. hy(i,j)+xe(i,jp1) > 0.0_rsh )) &
         htxy(i,j)=hy(i,j)+MAX(xe(i,jp0),xe(i,jp1))

         IF (htxy(i,j) > 0.0_rsh) THEN
           igd(i,j)=1
           IF (htxy(i,j) < hm) igd(i,j)=5
           IF (htxy(i,j) < MIN(epsilon,REAL(hm,rlg))) igd(i,j)=4
#ifdef key_dyn_adv_linear
           hexy(i,j)=MAX(0.0_rsh,hy(i,j))*govlet(i,j)
#if defined key_siggen || defined key_gencoord
           DO k=1,kmax
             hzexy(k,i,j)=f_hzu_y(hy(i,j),dhy(i,j),k,i,j)*govlet(i,j)
           END DO
#endif
#else
           hexy(i,j)=(htxy(i,j)+dhy(i,j))*govlet(i,j)
#if defined key_siggen || defined key_gencoord
           IF ( hy(i,j)+0.5_rsh*(xe(i,jp0)+xe(i,jp1)) < 0.0_rsh ) THEN
             DO k=1,kmax
               hzexy(k,i,j)=f_hzu_y(hy(i,j),MAX(xe(i,jp0),xe(i,jp1))+dhy(i,j),k,i,j)*govlet(i,j)
             END DO
           ELSE
             DO k=1,kmax
               hzexy(k,i,j)=f_hzu_y(hy(i,j),0.5_rsh*(xe(i,jp0)+xe(i,jp1))+dhy(i,j),k,i,j)*govlet(i,j)
             END DO
           END IF
#endif
#endif
         ELSE
           igd(i,j)=4
           hexy(i,j)=0.0_rsh
#if defined key_siggen || defined key_gencoord
           hzexy(:,i,j)=0.0_rsh
#endif
         END IF
       END DO
     END DO
!$OMP END DO
   ELSE! test = 1
!$OMP DO
     DO j=MAX(jj1,ljminm1),MIN(jj2,ljmaxp1)
       DO i=MAX(ii1,liminm1),MIN(ii2,limaxp1)

         ip1 = i + 1
         ip0 = i
         IF (i==imax) ip1 = i
         IF (i==imin) ip0 = i + 1
         htxy(i,j)=hx(i,j)+0.5_rsh*(xe(ip0,j)+xe(ip1,j))

!       IF ( htxy(i,j)<=0.0 ) With such a test mn and ms are set to 1 more often as htx is smaller when negative
         IF ( htxy(i,j)<=0.0_rsh .AND. ( hx(i,j)+xe(ip0,j) > 0.0_rsh .OR. hx(i,j)+xe(ip1,j) > 0.0_rsh) ) &
         htxy(i,j)=hx(i,j)+MAX(xe(ip0,j),xe(ip1,j))


         IF (htxy(i,j) > 0.0_rsh) THEN
           igd(i,j)=1
           IF (htxy(i,j) < hm) igd(i,j)=5
           IF (htxy(i,j) < MIN(epsilon,REAL(hm,rlg))) igd(i,j)=4
#ifdef key_dyn_adv_linear
           hexy(i,j)=MAX(0.0_rsh,hx(i,j))*goulet(i,j)
#if defined key_siggen || defined key_gencoord
           DO k=1,kmax
             hzexy(k,i,j)=f_hzu_x(hx(i,j),dhx(i,j),k,i,j)*goulet(i,j)
           END DO
#endif
#else
           hexy(i,j)=(htxy(i,j)+dhx(i,j))*goulet(i,j)
#if defined key_siggen || defined key_gencoord
           IF ( hx(i,j)+0.5_rsh*(xe(ip0,j)+xe(ip1,j)) < 0.0_rsh ) THEN
             DO k=1,kmax
               hzexy(k,i,j)=f_hzu_x(hx(i,j),MAX(xe(ip0,j),xe(ip1,j))+dhx(i,j),k,i,j)*goulet(i,j)
             END DO
           ELSE
             DO k=1,kmax
               hzexy(k,i,j)=f_hzu_x(hx(i,j),0.5_rsh*(xe(ip0,j)+xe(ip1,j))+dhx(i,j),k,i,j)*goulet(i,j)
             END DO
           END IF
#endif
#endif
         ELSE
           igd(i,j)=4
           hexy(i,j)=0.0_rsh
#if defined key_siggen || defined key_gencoord
           hzexy(:,i,j)=0.0_rsh
#endif
         END IF
       END DO  ! ends loop on i

     END DO
!$OMP END DO
   END IF
!$OMP END PARALLEL

   PRINT_DBG*, 'exit AGRIF_DYN3D_HEIGHT_MASK routine'


#endif

  END SUBROUTINE agrif_dyn3d_height_mask

   !!======================================================================

  SUBROUTINE agrif_dyn3d_wz(xe3, ii1, ii2, jj1, jj2)

#include "toolcpp.h"

#ifdef key_agrif

   USE parameters
   USE comvars2d
   USE comvars3d
   USE comvarp2d
#if defined key_siggen || defined key_gencoord
   USE toolgeom,  ONLY : f_hzw
#endif
   
   IMPLICIT NONE

   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT ( in ) :: xe3
   INTEGER, INTENT(IN) :: ii1,ii2,jj1,jj2 

   INTEGER             :: k,i,j
   REAL(KIND=rsh)      :: d,xet,div,divt
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rsh)      :: hz3_w
#endif

   PRINT_DBG*, 'enter AGRIF_DYN3D_WZ routine'

  ! Estimate wz a cell outside the domain
   DO j=MAX(jj1-1,ljminp1),MIN(jj2+1,ljmax)
     DO i=MAX(ii1-1,liminp1),MIN(ii2+1,limax)

       d=h0(i,j)+xe3(i,j)

       IF (d > 0.0_rsh) THEN
         divt=0.0_rsh
         xet=-( SUM(dsigu(1:kmax)*fxeau(1:kmax,i  ,j))*dy_u(i  ,j) &
               -SUM(dsigu(1:kmax)*fxeau(1:kmax,i-1,j))*dy_u(i-1,j) &
               +SUM(dsigu(1:kmax)*fyeau(1:kmax,i  ,j))*dx_v(i  ,j) &
               -SUM(dsigu(1:kmax)*fyeau(1:kmax,i,j-1))*dx_v(i,j-1) &
               -SUM(phieau(1:kmax,i,j))*halfdt )/(surf(i,j)*halfdt)

         DO  k=1,kmax-1
         ! estimates vertical velocity (sigma)
           div=(fxeau(k,i,j)*dy_u(i,j)-fxeau(k,i-1,j)*dy_u(i-1,j))/surf(i,j)/halfdt &
              +(fyeau(k,i,j)*dx_v(i,j)-fyeau(k,i,j-1)*dx_v(i,j-1))/surf(i,j)/halfdt &
              -phieau(k,i,j)/(surf(i,j)*fwet(i,j))/dsigu(k)
           divt=divt+div*dsigu(k)
#if defined key_siggen || defined key_gencoord
           hz3_w=f_hzw(h0(i,j),xe3(i,j),k,i,j)
           IF (h0(i,j)+xe3(i,j) <= hm) hz3_w=hm/REAL(kmax,rsh)
           wz(k,i,j)=-(xet*(sigw(k)+1.0_rsh)+divt)/hz3_w
#else
           wz(k,i,j)=-(xet*(sigw(k)+1.0_rsh)+divt)/d
#endif
         END DO   ! ends loop on k

       ELSE
         wz(:,i,j)=0.0_rsh
       END IF   ! ends test on water depth

     END DO 
   END DO

   PRINT_DBG*, 'exit AGRIF_DYN3D_WZ routine'
#endif /* key_agrif */

  END SUBROUTINE agrif_dyn3d_wz

   !!======================================================================

