Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/1_namelist_cfg
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/1_namelist_cfg	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/1_namelist_cfg	(revision 8011)
@@ -0,0 +1,372 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OPA  Configuration namelist : used to overwrite defaults values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+   cn_exp      =  "VORTEX"    !  experience name
+   nn_date0    = 19500101     !  initial date
+   nn_it000    =        1     !  first time step
+   nn_itend    =    16200     !  last  time step (std 5475)
+   nn_stock    =        0     !  frequency of creation of a restart file (modulo referenced to 1)
+   nn_istate   =        0     !  output the initial state (1) or not (0)
+/
+!-----------------------------------------------------------------------
+&namagrif      !  AGRIF zoom                                            ("key_agrif")
+!-----------------------------------------------------------------------
+   nn_cln_update =    1    !  baroclinic update frequency
+   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
+   rn_sponge_tra =  300.   !  coefficient for tracer   sponge layer [m2/s]
+   rn_sponge_dyn =  300.   !  coefficient for dynamics sponge layer [m2/s]
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration
+!-----------------------------------------------------------------------
+   cp_cfg      =  "vortex"             !  name of the configuration
+   cp_cfz      =      ''               !  name of the zoom of configuration
+   jp_cfg      =       0               !  resolution of the configuration
+   jpidta      =      61               !  1st lateral dimension ( >= jpi )
+   jpjdta      =      61               !  2nd    "         "    ( >= jpj )
+   jpkdta      =      11               !  number of levels      ( >= jpk )
+   jpiglo      =      61               !  1st dimension of global domain --> i =jpidta
+   jpjglo      =      61               !  2nd    -                  -    --> j  =jpjdta
+   jpizoom     =       1               !  left bottom (i,j) indices of the zoom
+   jpjzoom     =       1               !  in data domain indices
+   jperio      =       0               !  lateral cond. type (between 0 and 6)
+/
+!-----------------------------------------------------------------------
+&namzgr        !   vertical coordinate
+!-----------------------------------------------------------------------
+   ln_zco      = .true.    !  z-coordinate - full    steps   (T/F)      ("key_zco" may also be defined)
+   ln_zps      = .false.   !  z-coordinate - partial steps   (T/F)
+/
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdom        !   space and time domain (bathymetry, mesh, timestep)
+!-----------------------------------------------------------------------
+   nn_msh      =    0      !  create (=1) a mesh file or not (=0)
+   nn_bathy    =    0      !  compute (=0) or read (=1) the bathymetry file
+   rn_bathy    =  5000.    !  value of the bathymetry. if (=0) bottom flat at jpkm1                       !
+   rn_rdt      =   480.    !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_atfp     =    0.05   !  asselin time filter parameter
+   jphgr_msh   =    3                  !  type of horizontal mesh
+                                       !  = 0 curvilinear coordinate on the sphere read in coordinate.nc
+                                       !  = 1 geographical mesh on the sphere with regular grid-spacing
+                                       !  = 2 f-plane with regular grid-spacing
+                                       !  = 3 beta-plane with regular grid-spacing
+                                       !  = 4 Mercator grid with T/U point at the equator
+   ppgphi0     =    38.5               ! latitude  of first raw and column T-point (jphgr_msh = 1)
+   ppe1_m      =   30000.0             !  zonal      grid-spacing (meters)
+   ppe2_m      =   30000.0             !  meridional grid-spacing (meters)
+   ppsur       =  999999.              !  ORCA r4, r2 and r05 coefficients
+   ppa0        =  999999.              ! (default coefficients)
+   ppa1        =  999999.              !
+   ppkth       =  0.                   ! =0 means uniform vertical grid
+   pphmax      =  5000.                !  Maximum depth
+   ldbletanh   =  .FALSE.              !  Use/do not use double tanf function for vertical coordinates
+
+/
+!-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+   ln_bt_nn_auto =    .FALSE.          !  Set nn_baro automatically to be just below
+                                       !  a user defined maximum courant number (rn_bt_cmax)
+   nn_baro       =    20               !  Number of iterations of barotropic mode
+                                       !  during rn_rdt seconds. Only used if ln_bt_nn_auto=F
+/
+!-----------------------------------------------------------------------
+&namcrs        !   Grid coarsening for dynamics output and/or
+               !   passive tracer coarsened online simulations
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtsd    !   data : Temperature  & Salinity
+!-----------------------------------------------------------------------
+   ln_tsd_init   = .false.   !  Initialisation of ocean T & S with T &S input data (T) or not (F)
+   ln_tsd_tradmp = .false.   !  damping of ocean T & S toward T &S input data (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namsbc        !   Surface Boundary Condition (surface module)
+!-----------------------------------------------------------------------
+   ln_ana      = .true.     !  analytical formulation                   (T => fill namsbc_ana )
+   ln_blk_core = .false.    !  CORE bulk formulation                    (T => fill namsbc_core)
+   nn_ice      = 0          !  =0 no ice boundary condition   ,
+   ln_rnf      = .false.    !  runoffs                                  (T   => fill namsbc_rnf)
+   ln_ssr      = .false.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
+   nn_fwb      = 0          !  FreshWater Budget: =0 unchecked
+/
+!-----------------------------------------------------------------------
+&namsbc_ana    !   analytical surface boundary condition
+!-----------------------------------------------------------------------
+   rn_utau0    =   0.e0    !  uniform value for the i-stress
+/
+!-----------------------------------------------------------------------
+&namsbc_core   !   namsbc_core  CORE bulk formulae
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtra_qsr    !   penetrative solar radiation
+!-----------------------------------------------------------------------
+   ln_traqsr   = .false.    !  Light penetration (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namsbc_rnf    !   runoffs namelist surface boundary condition
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_ssr    !   surface boundary condition : sea surface restoring
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_alb    !   albedo parameters
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namberg       !   iceberg parameters
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition
+!-----------------------------------------------------------------------
+   rn_shlat    =     0     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat
+                           !  free slip  !   partial slip  !   no slip   ! strong slip
+/
+!-----------------------------------------------------------------------
+&namcla        !   cross land advection
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nambfr        !   bottom friction
+!-----------------------------------------------------------------------
+   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
+                           !                              = 2 : nonlinear friction
+/
+!-----------------------------------------------------------------------
+&nambbc        !   bottom temperature boundary condition
+!-----------------------------------------------------------------------
+   ln_trabbc   = .false.    !  Apply a geothermal heating at the ocean bottom
+/
+!-----------------------------------------------------------------------
+&nambbl        !   bottom boundary layer scheme
+!-----------------------------------------------------------------------
+   nn_bbl_ldf  =  0      !  diffusive bbl (=1)   or not (=0)
+/
+!-----------------------------------------------------------------------
+&nameos        !   ocean physical parameters
+!-----------------------------------------------------------------------
+   nn_eos      =  1     !  type of equation of state and Brunt-Vaisala frequency
+                                 !  =-1, TEOS-10
+                                 !  = 0, EOS-80
+                                 !  = 1, S-EOS   (simplified eos)
+   ln_useCT    = .false. ! use of Conservative Temp. ==> surface CT converted in Pot. Temp. in sbcssm
+   !                             !
+   !                     ! S-EOS coefficients :
+   !                             !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
+   rn_a0       =  0.28           !  thermal expension coefficient (nn_eos= 1)
+   rn_b0       =  0.             !  saline  expension coefficient (nn_eos= 1)
+   rn_lambda1  =  0.             !  cabbeling coeff in T^2  (=0 for linear eos)
+   rn_lambda2  =  0.             !  cabbeling coeff in S^2  (=0 for linear eos)
+   rn_mu1      =  0.             !  thermobaric coeff. in T (=0 for linear eos)
+   rn_mu2      =  0.             !  thermobaric coeff. in S (=0 for linear eos)
+   rn_nu       =  0.             !  cabbeling coeff in T*S  (=0 for linear eos)
+/
+!-----------------------------------------------------------------------
+&namtra_adv    !   advection scheme for tracer
+!-----------------------------------------------------------------------
+   ln_traadv_cen2   =  .false.   !  2nd order centered scheme
+   ln_traadv_tvd    =  .true.    !  TVD scheme
+   ln_traadv_ubs    =  .false.   !  UBS scheme
+/
+!-----------------------------------------------------------------------
+&namtra_adv_mle !  mixed layer eddy parametrisation (Fox-Kemper param)
+!-----------------------------------------------------------------------
+   ln_mle    = .false.      ! (T) use the Mixed Layer Eddy (MLE) parameterisation
+/
+!-----------------------------------------------------------------------
+&namtra_adv_mle !  mixed layer eddy parametrisation (Fox-Kemper param)
+!-----------------------------------------------------------------------
+/
+!----------------------------------------------------------------------------------
+&namtra_ldf    !   lateral diffusion scheme for tracers
+!----------------------------------------------------------------------------------
+   !                       !  Operator type:
+   ln_traldf_lap    =  .false.  !  laplacian operator
+   ln_traldf_bilap  =  .true.   !  bilaplacian operator
+   !                       !  Direction of action:
+   ln_traldf_level  =  .false.  !  iso-level
+   ln_traldf_hor    =  .true.   !  horizontal (geopotential)   (needs "key_ldfslp" when ln_sco=T)
+   ln_traldf_iso    =  .false.  !  iso-neutral                 (needs "key_ldfslp")
+   rn_aht_0         =   -2.e9  !  horizontal eddy diffusivity for tracers [m2/s]
+/
+!-----------------------------------------------------------------------
+&namtra_dmp    !   tracer: T & S newtonian damping
+!-----------------------------------------------------------------------
+   ln_tradmp   =  .false.   !  add a damping termn (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namdyn_adv    !   formulation of the momentum advection
+!-----------------------------------------------------------------------
+   ln_dynadv_vec = .false.  !  vector form (T) or flux form (F)
+   nn_dynkeg     = 0       ! scheme for grad(KE): =0   C2  ;  =1   Hollingsworth correction
+   ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
+   ln_dynadv_ubs = .true. !  flux form - 3rd order UBS      scheme
+   ln_dynzad_zts = .false. !  Use (T) sub timestepping for vertical momentum advection
+/
+!-----------------------------------------------------------------------
+&namdyn_vor    !   option of physics/algorithm (not control by CPP keys)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdyn_hpg    !   Hydrostatic pressure gradient option
+!-----------------------------------------------------------------------
+   ln_hpg_zps  = .false.   !  z-coordinate - partial steps (interpolation)
+   ln_hpg_sco  = .true.    !  s-coordinate (standard jacobian formulation)
+/
+!-----------------------------------------------------------------------
+&namdyn_ldf    !   lateral diffusion on momentum
+!-----------------------------------------------------------------------
+   !                       !  Type of the operator :
+   ln_dynldf_lap    =  .false.  !  laplacian operator
+   ln_dynldf_bilap  =  .false.  !  bilaplacian operator
+   !                       !  Direction of action  :
+   ln_dynldf_level  =  .false.  !  iso-level
+   ln_dynldf_hor    =  .true.   !  horizontal (geopotential)            (require "key_ldfslp" in s-coord.)
+   ln_dynldf_iso    =  .false.  !  iso-neutral                          (require "key_ldfslp")
+   !                       !  Coefficient
+/
+!-----------------------------------------------------------------------
+&namzdf        !   vertical physics
+!-----------------------------------------------------------------------
+   rn_avm0     =    1.e-4  !  vertical eddy viscosity   [m2/s]          (background Kz if not "key_zdfcst")
+   rn_avt0     =    1.e-5  !  vertical eddy diffusivity [m2/s]          (background Kz if not "key_zdfcst")
+   ln_zdfevd   = .true.    !  enhanced vertical diffusion (evd) (T) or not (F)
+   nn_evdm     =    0      !  evd apply on tracer (=0) or on tracer and momentum (=1)
+   rn_avevd    =    0.     !  evd mixing coefficient [m2/s]
+/
+!-----------------------------------------------------------------------
+&nambdy        !  unstructured open boundaries                          ("key_bdy")
+    nb_bdy         = 0                    !  number of open boundary sets
+    ln_coords_file = .false.,.false.,.false.,.false.      !  =T : read bdy coordinates from file
+    cn_coords_file = '','','',''          !  bdy coordinates files
+    ln_mask_file   = .false.              !  =T : read mask from file
+    cn_mask_file   = ''                   !  name of mask file (if ln_mask_file=.TRUE.)
+    cn_dyn2d       = 'flather','flather','flather','flather' !
+    nn_dyn2d_dta   =  0,0,0,0             !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+                                          !  = 2, use tidal harmonic forcing data from files
+                                          !  = 3, use external data AND tidal harmonic forcing
+    cn_dyn3d      =  'specified','specified','specified','specified'      !
+    nn_dyn3d_dta  =  0,0,0,0              !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    cn_tra        =  'specified','specified','specified','specified'      !
+    nn_tra_dta    =  0,0,0,0              !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    ln_tra_dmp    =.false.,.false.,.false.,.false.      !  open boudaries conditions for tracers
+    ln_dyn3d_dmp  =.false.,.false.,.false.,.false.      !  open boundary condition for baroclinic velocities
+    rn_time_dmp   =  1., 1., 1., 1.             ! Damping time scale in days
+    rn_time_dmp_out =  30., 30., 30., 30.       ! Outflow damping time scale
+    nn_rimwidth   = 1, 1, 1, 1                  !  width of the relaxation zone
+    ln_vol        = .false.               !  total volume correction (see nn_volctl parameter)
+    nn_volctl     = 1                     !  = 0, the total water flux across open boundaries is zero
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='S'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='N'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='W'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='E'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  ("key_zdftke")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_ddm    !   double diffusive mixing parameterization             ("key_zdfddm")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_tmx    !   tidal mixing parameterization                        ("key_zdftmx")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsol        !   elliptic solver / island / free surface
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl        !   Control prints & Benchmark
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namptr       !   Poleward Transport Diagnostic
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namhsb       !  Heat and salt budgets
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namobs       !  observation usage                                      ('key_diaobs')
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nam_asminc   !   assimilation increments                               ('key_asminc')
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsto       ! Stochastic parametrization of EOS
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrd        !   diagnostics on dynamics and/or tracer trends
+!              !       and/or mixed-layer trends and/or barotropic vorticity
+!-----------------------------------------------------------------------
+/
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/1_namelist_ref
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/1_namelist_ref	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/1_namelist_ref	(revision 8011)
@@ -0,0 +1,1 @@
+link ../../SHARED/1_namelist_ref
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/AGRIF_FixedGrids.in
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/AGRIF_FixedGrids.in	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/AGRIF_FixedGrids.in	(revision 8011)
@@ -0,0 +1,5 @@
+1
+20 39 20 39 3 3 3 
+0
+10 49 10 49 3 3 3
+0
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/domain_def.xml
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/domain_def.xml	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/domain_def.xml	(revision 8011)
@@ -0,0 +1,1 @@
+link ../../SHARED/domain_def.xml
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/field_def.xml
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/field_def.xml	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/field_def.xml	(revision 8011)
@@ -0,0 +1,1 @@
+link ../../SHARED/field_def.xml
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/iodef.xml
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/iodef.xml	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/iodef.xml	(revision 8011)
@@ -0,0 +1,310 @@
+<?xml version="1.0"?>
+<simulation> 
+
+ <context id="nemo" time_origin="1950-01-01 00:00:00" >
+    
+    <!-- $id$ -->
+    
+    <!-- 
+============================================================================================================
+=                                  definition of all existing variables                                    =
+=                                            DO NOT CHANGE                                                 =
+============================================================================================================
+    -->
+    <field_definition src="./field_def.xml"/>
+    <!-- 
+============================================================================================================
+=                                           output files definition                                        =
+=                                            Define your own files                                         =
+=                                         put the variables you want...                                    =
+============================================================================================================
+    -->
+    
+    <file_definition type="multiple_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="10d" min_digits="4">
+    
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE." > <!-- 1 time step files -->
+        <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" >
+          <field field_ref="ssh"          name="zos"    />
+          <field field_ref="toce"         name="thetao" /> 
+          <field field_ref="e3t"         name="e3t" />
+        </file>
+        <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
+          <field field_ref="ubar"         name="ubar"   />
+          <field field_ref="uoce"         name="uoce"   />
+        </file>
+        <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
+          <field field_ref="vbar"         name="vbar"   />
+          <field field_ref="voce"         name="voce"   />
+        </file>
+        <file id="file4" name_suffix="_grid_W" description="ocean V grid variables" >
+          <field field_ref="woce"         name="woce"   />
+        </file>
+      </file_group>
+
+      <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1h files -->
+      <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->     
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files -->
+      <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->    
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/> <!-- 5d files -->  
+      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files -->
+      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m" output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m" output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+
+      <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."/> <!-- real yearly files -->
+      <file_group id="2y"  output_freq="2y" output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y" output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
+    
+    <!-- 
+============================================================================================================
+= grid definition = = DO NOT CHANGE =
+============================================================================================================
+    -->
+    
+   <axis_definition>  
+      <axis id="deptht"  long_name="Vertical T levels"  unit="m" positive="down" />
+      <axis id="depthu"  long_name="Vertical U levels"  unit="m" positive="down" />
+      <axis id="depthv"  long_name="Vertical V levels"  unit="m" positive="down" />
+      <axis id="depthw"  long_name="Vertical W levels"  unit="m" positive="down" />
+      <axis id="nfloat"  long_name="Float number"       unit="1"                 />
+      <axis id="icbcla"  long_name="Iceberg class"      unit="1"                 />
+      <axis id="ncatice" long_name="Ice category"       unit="1"                 />
+      <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"              />
+      <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"              />
+   </axis_definition> 
+   
+   <domain_definition src="./domain_def.xml"/>
+   
+   <grid_definition>    
+     <grid id="grid_T_2D" domain_ref="grid_T"/>
+     <grid id="grid_T_3D" domain_ref="grid_T" axis_ref="deptht"/>
+     <grid id="grid_U_2D" domain_ref="grid_U"/>
+     <grid id="grid_U_3D" domain_ref="grid_U" axis_ref="depthu"/>
+     <grid id="grid_V_2D" domain_ref="grid_V"/>
+     <grid id="grid_V_3D" domain_ref="grid_V" axis_ref="depthv"/>
+     <grid id="grid_W_2D" domain_ref="grid_W"/>
+     <grid id="grid_W_3D" domain_ref="grid_W" axis_ref="depthw"/>
+     <grid id="gznl_T_2D" domain_ref="gznl"/>
+     <grid id="gznl_T_3D" domain_ref="gznl" axis_ref="deptht"/>
+     <grid id="gznl_W_3D" domain_ref="gznl" axis_ref="depthw"/>
+    </grid_definition>   
+  </context>
+  
+ <context id="1_nemo" time_origin="1950-01-01 00:00:00" >
+    
+    <!-- $id$ -->
+    
+    <!-- 
+============================================================================================================
+=                                  definition of all existing variables                                    =
+=                                            DO NOT CHANGE                                                 =
+============================================================================================================
+    -->
+    <field_definition src="./field_def.xml"/>
+    <!-- 
+============================================================================================================
+=                                           output files definition                                        =
+=                                            Define your own files                                         =
+=                                         put the variables you want...                                    =
+============================================================================================================
+    -->
+    
+    <file_definition type="multiple_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="10d" min_digits="4">
+    
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE." > <!-- 1 time step files -->
+        <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" >
+          <field field_ref="ssh"          name="zos"    />
+          <field field_ref="toce"         name="thetao" /> 
+          <field field_ref="e3t"          name="e3t"    />
+        </file>
+        <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
+          <field field_ref="ubar"         name="ubar"   />
+          <field field_ref="uoce"         name="uoce"   />
+        </file>
+        <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
+          <field field_ref="vbar"         name="vbar"   />
+          <field field_ref="voce"         name="voce"   />
+        </file>
+        <file id="file4" name_suffix="_grid_W" description="ocean V grid variables" >
+          <field field_ref="woce"         name="woce"   />
+        </file>
+      </file_group>
+
+      <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1h files -->
+      <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->     
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files -->
+      <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->    
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/> <!-- 5d files -->  
+      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files -->
+      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m" output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m" output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+
+      <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."/> <!-- real yearly files -->
+      <file_group id="2y"  output_freq="2y" output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y" output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
+    
+    <!-- 
+============================================================================================================
+= grid definition = = DO NOT CHANGE =
+============================================================================================================
+    -->
+    
+   <axis_definition>  
+      <axis id="deptht"  long_name="Vertical T levels"  unit="m" positive="down" />
+      <axis id="depthu"  long_name="Vertical U levels"  unit="m" positive="down" />
+      <axis id="depthv"  long_name="Vertical V levels"  unit="m" positive="down" />
+      <axis id="depthw"  long_name="Vertical W levels"  unit="m" positive="down" />
+      <axis id="nfloat"  long_name="Float number"       unit="1"                 />
+      <axis id="icbcla"  long_name="Iceberg class"      unit="1"                 />
+      <axis id="ncatice" long_name="Ice category"       unit="1"                 />
+      <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"              />
+      <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"              />
+   </axis_definition> 
+   
+   <domain_definition src="./domain_def.xml"/>
+   
+   <grid_definition>    
+     <grid id="grid_T_2D" domain_ref="grid_T"/>
+     <grid id="grid_T_3D" domain_ref="grid_T" axis_ref="deptht"/>
+     <grid id="grid_U_2D" domain_ref="grid_U"/>
+     <grid id="grid_U_3D" domain_ref="grid_U" axis_ref="depthu"/>
+     <grid id="grid_V_2D" domain_ref="grid_V"/>
+     <grid id="grid_V_3D" domain_ref="grid_V" axis_ref="depthv"/>
+     <grid id="grid_W_2D" domain_ref="grid_W"/>
+     <grid id="grid_W_3D" domain_ref="grid_W" axis_ref="depthw"/>
+     <grid id="gznl_T_2D" domain_ref="gznl"/>
+     <grid id="gznl_T_3D" domain_ref="gznl" axis_ref="deptht"/>
+     <grid id="gznl_W_3D" domain_ref="gznl" axis_ref="depthw"/>
+    </grid_definition>   
+  </context> 
+
+ <context id="2_nemo" time_origin="1950-01-01 00:00:00" >
+    
+    <!-- $id$ -->
+    
+    <!-- 
+============================================================================================================
+=                                  definition of all existing variables                                    =
+=                                            DO NOT CHANGE                                                 =
+============================================================================================================
+    -->
+    <field_definition src="./field_def.xml"/>
+    <!-- 
+============================================================================================================
+=                                           output files definition                                        =
+=                                            Define your own files                                         =
+=                                         put the variables you want...                                    =
+============================================================================================================
+    -->
+    
+    <file_definition type="multiple_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="10d" min_digits="4">
+    
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE." > <!-- 1 time step files -->
+        <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" >
+          <field field_ref="ssh"          name="zos"    />
+          <field field_ref="toce"         name="thetao" /> 
+          <field field_ref="e3t"          name="e3t"    />
+        </file>
+        <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
+          <field field_ref="ubar"         name="ubar"   />
+          <field field_ref="uoce"         name="uoce"   />
+        </file>
+        <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
+          <field field_ref="vbar"         name="vbar"   />
+          <field field_ref="voce"         name="voce"   />
+        </file>
+        <file id="file4" name_suffix="_grid_W" description="ocean V grid variables" >
+          <field field_ref="woce"         name="woce"   />
+        </file>
+      </file_group>
+
+      <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1h files -->
+      <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->     
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files -->
+      <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->    
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/> <!-- 5d files -->  
+      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files -->
+      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m" output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m" output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+
+      <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."/> <!-- real yearly files -->
+      <file_group id="2y"  output_freq="2y" output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y" output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
+    
+    <!-- 
+============================================================================================================
+= grid definition = = DO NOT CHANGE =
+============================================================================================================
+    -->
+    
+   <axis_definition>  
+      <axis id="deptht"  long_name="Vertical T levels"  unit="m" positive="down" />
+      <axis id="depthu"  long_name="Vertical U levels"  unit="m" positive="down" />
+      <axis id="depthv"  long_name="Vertical V levels"  unit="m" positive="down" />
+      <axis id="depthw"  long_name="Vertical W levels"  unit="m" positive="down" />
+      <axis id="nfloat"  long_name="Float number"       unit="1"                 />
+      <axis id="icbcla"  long_name="Iceberg class"      unit="1"                 />
+      <axis id="ncatice" long_name="Ice category"       unit="1"                 />
+      <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"              />
+      <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"              />
+   </axis_definition> 
+   
+   <domain_definition src="./domain_def.xml"/>
+   
+   <grid_definition>    
+     <grid id="grid_T_2D" domain_ref="grid_T"/>
+     <grid id="grid_T_3D" domain_ref="grid_T" axis_ref="deptht"/>
+     <grid id="grid_U_2D" domain_ref="grid_U"/>
+     <grid id="grid_U_3D" domain_ref="grid_U" axis_ref="depthu"/>
+     <grid id="grid_V_2D" domain_ref="grid_V"/>
+     <grid id="grid_V_3D" domain_ref="grid_V" axis_ref="depthv"/>
+     <grid id="grid_W_2D" domain_ref="grid_W"/>
+     <grid id="grid_W_3D" domain_ref="grid_W" axis_ref="depthw"/>
+     <grid id="gznl_T_2D" domain_ref="gznl"/>
+     <grid id="gznl_T_3D" domain_ref="gznl" axis_ref="deptht"/>
+     <grid id="gznl_W_3D" domain_ref="gznl" axis_ref="depthw"/>
+    </grid_definition>   
+  </context> 
+
+  <context id="xios">
+
+      <variable_definition>
+	
+     <!-- 
+        We must have buffer_size > jpi*jpj*jpk*8 (with jpi and jpj the subdomain size)
+-->
+	  <variable id="buffer_size"               type="integer">50000000</variable>
+	  <variable id="buffer_server_factor_size" type="integer">2</variable>
+	  <variable id="info_level"                type="integer">0</variable>
+	  <variable id="using_server"              type="boolean">false</variable>
+	  <variable id="using_oasis"               type="boolean">false</variable>
+	  <variable id="oasis_codes_id"            type="string" >oceanx</variable>
+	
+      </variable_definition>
+               
+  </context>
+  
+</simulation>
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/namelist_cfg	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/namelist_cfg	(revision 8011)
@@ -0,0 +1,364 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OPA  Configuration namelist : used to overwrite defaults values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+   cn_exp      =  "VORTEX"    !  experience name
+   nn_date0    = 19500101     !  initial date
+   nn_it000    =        1     !  first time step
+   nn_itend    =     5400     !  last  time step (std 5475)
+   nn_stock    =        0     !  frequency of creation of a restart file (modulo referenced to 1)
+   nn_istate   =        0     !  output the initial state (1) or not (0)
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration
+!-----------------------------------------------------------------------
+   cp_cfg      =  "vortex"             !  name of the configuration
+   cp_cfz      =      ''               !  name of the zoom of configuration
+   jp_cfg      =       0               !  resolution of the configuration
+   jpidta      =      61               !  1st lateral dimension ( >= jpi )
+   jpjdta      =      61               !  2nd    "         "    ( >= jpj )
+   jpkdta      =      11               !  number of levels      ( >= jpk )
+   jpiglo      =      61               !  1st dimension of global domain --> i =jpidta
+   jpjglo      =      61               !  2nd    -                  -    --> j  =jpjdta
+   jpizoom     =       1               !  left bottom (i,j) indices of the zoom
+   jpjzoom     =       1               !  in data domain indices
+   jperio      =       0               !  lateral cond. type (between 0 and 6)
+/
+!-----------------------------------------------------------------------
+&namzgr        !   vertical coordinate
+!-----------------------------------------------------------------------
+   ln_zco      = .true.    !  z-coordinate - full    steps   (T/F)      ("key_zco" may also be defined)
+   ln_zps      = .false.   !  z-coordinate - partial steps   (T/F)
+/
+!-----------------------------------------------------------------------
+&nam_vvl    !   vertical coordinate options
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdom        !   space and time domain (bathymetry, mesh, timestep)
+!-----------------------------------------------------------------------
+   nn_msh      =    0      !  create (=1) a mesh file or not (=0)
+   nn_bathy    =    0      !  compute (=0) or read (=1) the bathymetry file
+   rn_bathy    =  5000.    !  value of the bathymetry. if (=0) bottom flat at jpkm1                       !
+   rn_rdt      =  1440.    !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_atfp     =    0.05   !  asselin time filter parameter
+   jphgr_msh   =    3                  !  type of horizontal mesh
+                                       !  = 0 curvilinear coordinate on the sphere read in coordinate.nc
+                                       !  = 1 geographical mesh on the sphere with regular grid-spacing
+                                       !  = 2 f-plane with regular grid-spacing
+                                       !  = 3 beta-plane with regular grid-spacing
+                                       !  = 4 Mercator grid with T/U point at the equator
+   ppgphi0     =    38.5               ! latitude  of first raw and column T-point (jphgr_msh = 1)
+   ppe1_m      =   30000.0             !  zonal      grid-spacing (meters)
+   ppe2_m      =   30000.0             !  meridional grid-spacing (meters)
+   ppsur       =  999999.              !  ORCA r4, r2 and r05 coefficients
+   ppa0        =  999999.              ! (default coefficients)
+   ppa1        =  999999.              !
+   ppkth       =  0.                   ! =0 means uniform vertical grid
+   pphmax      =  5000.                !  Maximum depth
+   ldbletanh   =  .FALSE.              !  Use/do not use double tanf function for vertical coordinates
+
+/
+!-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+   ln_bt_nn_auto =    .FALSE.          !  Set nn_baro automatically to be just below
+                                       !  a user defined maximum courant number (rn_bt_cmax)
+   nn_baro       =    20               !  Number of iterations of barotropic mode
+                                       !  during rn_rdt seconds. Only used if ln_bt_nn_auto=F
+/
+!-----------------------------------------------------------------------
+&namcrs        !   Grid coarsening for dynamics output and/or
+               !   passive tracer coarsened online simulations
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtsd    !   data : Temperature  & Salinity
+!-----------------------------------------------------------------------
+   ln_tsd_init   = .false.   !  Initialisation of ocean T & S with T &S input data (T) or not (F)
+   ln_tsd_tradmp = .false.   !  damping of ocean T & S toward T &S input data (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namsbc        !   Surface Boundary Condition (surface module)
+!-----------------------------------------------------------------------
+   ln_ana      = .true.     !  analytical formulation                   (T => fill namsbc_ana )
+   ln_blk_core = .false.    !  CORE bulk formulation                    (T => fill namsbc_core)
+   nn_ice      = 0          !  =0 no ice boundary condition   ,
+   ln_rnf      = .false.    !  runoffs                                  (T   => fill namsbc_rnf)
+   ln_ssr      = .false.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
+   nn_fwb      = 0          !  FreshWater Budget: =0 unchecked
+/
+!-----------------------------------------------------------------------
+&namsbc_ana    !   analytical surface boundary condition
+!-----------------------------------------------------------------------
+   rn_utau0    =   0.e0    !  uniform value for the i-stress
+/
+!-----------------------------------------------------------------------
+&namsbc_core   !   namsbc_core  CORE bulk formulae
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtra_qsr    !   penetrative solar radiation
+!-----------------------------------------------------------------------
+   ln_traqsr   = .false.    !  Light penetration (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namsbc_rnf    !   runoffs namelist surface boundary condition
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_ssr    !   surface boundary condition : sea surface restoring
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_alb    !   albedo parameters
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namberg       !   iceberg parameters
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition
+!-----------------------------------------------------------------------
+   rn_shlat    =     0     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat
+                           !  free slip  !   partial slip  !   no slip   ! strong slip
+/
+!-----------------------------------------------------------------------
+&namcla        !   cross land advection
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nambfr        !   bottom friction
+!-----------------------------------------------------------------------
+   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
+                           !                              = 2 : nonlinear friction
+/
+!-----------------------------------------------------------------------
+&nambbc        !   bottom temperature boundary condition
+!-----------------------------------------------------------------------
+   ln_trabbc   = .false.    !  Apply a geothermal heating at the ocean bottom
+/
+!-----------------------------------------------------------------------
+&nambbl        !   bottom boundary layer scheme
+!-----------------------------------------------------------------------
+   nn_bbl_ldf  =  0      !  diffusive bbl (=1)   or not (=0)
+/
+!-----------------------------------------------------------------------
+&nameos        !   ocean physical parameters
+!-----------------------------------------------------------------------
+   nn_eos      =  1     !  type of equation of state and Brunt-Vaisala frequency
+                                 !  =-1, TEOS-10
+                                 !  = 0, EOS-80
+                                 !  = 1, S-EOS   (simplified eos)
+   ln_useCT    = .false. ! use of Conservative Temp. ==> surface CT converted in Pot. Temp. in sbcssm
+   !                             !
+   !                     ! S-EOS coefficients :
+   !                             !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
+   rn_a0       =  0.28           !  thermal expension coefficient (nn_eos= 1)
+   rn_b0       =  0.             !  saline  expension coefficient (nn_eos= 1)
+   rn_lambda1  =  0.             !  cabbeling coeff in T^2  (=0 for linear eos)
+   rn_lambda2  =  0.             !  cabbeling coeff in S^2  (=0 for linear eos)
+   rn_mu1      =  0.             !  thermobaric coeff. in T (=0 for linear eos)
+   rn_mu2      =  0.             !  thermobaric coeff. in S (=0 for linear eos)
+   rn_nu       =  0.             !  cabbeling coeff in T*S  (=0 for linear eos)
+/
+!-----------------------------------------------------------------------
+&namtra_adv    !   advection scheme for tracer
+!-----------------------------------------------------------------------
+   ln_traadv_cen2   =  .false.   !  2nd order centered scheme
+   ln_traadv_tvd    =  .true.    !  TVD scheme
+   ln_traadv_ubs    =  .false.   !  UBS scheme
+/
+!-----------------------------------------------------------------------
+&namtra_adv_mle !  mixed layer eddy parametrisation (Fox-Kemper param)
+!-----------------------------------------------------------------------
+   ln_mle    = .false.      ! (T) use the Mixed Layer Eddy (MLE) parameterisation
+/
+!-----------------------------------------------------------------------
+&namtra_adv_mle !  mixed layer eddy parametrisation (Fox-Kemper param)
+!-----------------------------------------------------------------------
+/
+!----------------------------------------------------------------------------------
+&namtra_ldf    !   lateral diffusion scheme for tracers
+!----------------------------------------------------------------------------------
+   !                       !  Operator type:
+   ln_traldf_lap    =  .false.  !  laplacian operator
+   ln_traldf_bilap  =  .true.   !  bilaplacian operator
+   !                       !  Direction of action:
+   ln_traldf_level  =  .false.  !  iso-level
+   ln_traldf_hor    =  .true.   !  horizontal (geopotential)   (needs "key_ldfslp" when ln_sco=T)
+   ln_traldf_iso    =  .false.  !  iso-neutral                 (needs "key_ldfslp")
+   rn_aht_0         =   -5.e10  !  horizontal eddy diffusivity for tracers [m2/s]
+/
+!-----------------------------------------------------------------------
+&namtra_dmp    !   tracer: T & S newtonian damping
+!-----------------------------------------------------------------------
+   ln_tradmp   =  .false.   !  add a damping termn (T) or not (F)
+/
+!-----------------------------------------------------------------------
+&namdyn_adv    !   formulation of the momentum advection
+!-----------------------------------------------------------------------
+   ln_dynadv_vec = .false.  !  vector form (T) or flux form (F)
+   nn_dynkeg     = 0       ! scheme for grad(KE): =0   C2  ;  =1   Hollingsworth correction
+   ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
+   ln_dynadv_ubs = .true. !  flux form - 3rd order UBS      scheme
+   ln_dynzad_zts = .false. !  Use (T) sub timestepping for vertical momentum advection
+/
+!-----------------------------------------------------------------------
+&namdyn_vor    !   option of physics/algorithm (not control by CPP keys)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdyn_hpg    !   Hydrostatic pressure gradient option
+!-----------------------------------------------------------------------
+   ln_hpg_zps  = .false.   !  z-coordinate - partial steps (interpolation)
+   ln_hpg_sco  = .true.    !  s-coordinate (standard jacobian formulation)
+/
+!-----------------------------------------------------------------------
+&namdyn_ldf    !   lateral diffusion on momentum
+!-----------------------------------------------------------------------
+   !                       !  Type of the operator :
+   ln_dynldf_lap    =  .false.  !  laplacian operator
+   ln_dynldf_bilap  =  .false.  !  bilaplacian operator
+   !                       !  Direction of action  :
+   ln_dynldf_level  =  .false.  !  iso-level
+   ln_dynldf_hor    =  .true.   !  horizontal (geopotential)            (require "key_ldfslp" in s-coord.)
+   ln_dynldf_iso    =  .false.  !  iso-neutral                          (require "key_ldfslp")
+   !                       !  Coefficient
+/
+!-----------------------------------------------------------------------
+&namzdf        !   vertical physics
+!-----------------------------------------------------------------------
+   rn_avm0     =    1.e-4  !  vertical eddy viscosity   [m2/s]          (background Kz if not "key_zdfcst")
+   rn_avt0     =    1.e-5  !  vertical eddy diffusivity [m2/s]          (background Kz if not "key_zdfcst")
+   ln_zdfevd   = .true.    !  enhanced vertical diffusion (evd) (T) or not (F)
+   nn_evdm     =    0      !  evd apply on tracer (=0) or on tracer and momentum (=1)
+   rn_avevd    =    0.     !  evd mixing coefficient [m2/s]
+/
+!-----------------------------------------------------------------------
+&nambdy        !  unstructured open boundaries                          ("key_bdy")
+    nb_bdy         = 4                    !  number of open boundary sets
+    ln_coords_file = .false.,.false.,.false.,.false.      !  =T : read bdy coordinates from file
+    cn_coords_file = '','','',''          !  bdy coordinates files
+    ln_mask_file   = .false.              !  =T : read mask from file
+    cn_mask_file   = ''                   !  name of mask file (if ln_mask_file=.TRUE.)
+    cn_dyn2d       = 'flather','flather','flather','flather' !
+    nn_dyn2d_dta   =  0,0,0,0             !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+                                          !  = 2, use tidal harmonic forcing data from files
+                                          !  = 3, use external data AND tidal harmonic forcing
+    cn_dyn3d      =  'specified','specified','specified','specified'      !
+    nn_dyn3d_dta  =  0,0,0,0              !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    cn_tra        =  'specified','specified','specified','specified'      !
+    nn_tra_dta    =  0,0,0,0              !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    ln_tra_dmp    =.false.,.false.,.false.,.false.      !  open boudaries conditions for tracers
+    ln_dyn3d_dmp  =.false.,.false.,.false.,.false.      !  open boundary condition for baroclinic velocities
+    rn_time_dmp   =  100., 100., 100., 100.             ! Damping time scale in days
+    rn_time_dmp_out =  100., 100., 100., 100.       ! Outflow damping time scale
+    nn_rimwidth   = 1, 1, 1, 1                  !  width of the relaxation zone
+    ln_vol        = .false.               !  total volume correction (see nn_volctl parameter)
+    nn_volctl     = 1                     !  = 0, the total water flux across open boundaries is zero
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='S'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='N'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='W'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&nambdy_index  !  structured open boundaries definition     ("key_bdy")
+!-----------------------------------------------------------------------
+    ctypebdy ='E'                   ! Open boundary type (W,E,S or N)
+    nbdyind  = -1                   ! indice of velocity row or column
+                                            ! if ==-1, set obc at the domain boundary
+                                            !        , discard start and end indices
+    nbdybeg  =  -99                 ! indice of segment start
+    nbdyend  =  -99                 ! indice of segment end
+/
+!-----------------------------------------------------------------------
+&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  ("key_zdftke")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_ddm    !   double diffusive mixing parameterization             ("key_zdfddm")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_tmx    !   tidal mixing parameterization                        ("key_zdftmx")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsol        !   elliptic solver / island / free surface
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl        !   Control prints & Benchmark
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namptr       !   Poleward Transport Diagnostic
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namhsb       !  Heat and salt budgets
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namobs       !  observation usage                                      ('key_diaobs')
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nam_asminc   !   assimilation increments                               ('key_asminc')
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsto       ! Stochastic parametrization of EOS
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrd        !   diagnostics on dynamics and/or tracer trends
+!              !       and/or mixed-layer trends and/or barotropic vorticity
+!-----------------------------------------------------------------------
+/
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/namelist_ref
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/namelist_ref	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/EXP00/namelist_ref	(revision 8011)
@@ -0,0 +1,1 @@
+link ../../SHARED/namelist_ref
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/domhgr.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/domhgr.F90	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/domhgr.F90	(revision 8011)
@@ -0,0 +1,704 @@
+MODULE domhgr
+   !!==============================================================================
+   !!                       ***  MODULE domhgr   ***
+   !! Ocean initialization : domain initialization
+   !!==============================================================================
+   !! History :  OPA  ! 1988-03  (G. Madec) Original code
+   !!            7.0  ! 1996-01  (G. Madec)  terrain following coordinates
+   !!            8.0  ! 1997-02  (G. Madec)  print mesh informations
+   !!            8.1  ! 1999-11  (M. Imbard) NetCDF format with IO-IPSL
+   !!            8.2  ! 2000-08  (D. Ludicone) Reduced section at Bab el Mandeb
+   !!             -   ! 2001-09  (M. Levy)  eel config: grid in km, beta-plane
+   !!  NEMO      1.0  ! 2002-08  (G. Madec)  F90: Free form and module, namelist
+   !!             -   ! 2004-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
+   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist
+   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration 
+   !!            4.0  ! 2011-02  (G. Madec) add cell surface (e1e2t)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_hgr       : initialize the horizontal mesh 
+   !!   hgr_read      : read "coordinate" NetCDF file 
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
+
+   PUBLIC   dom_hgr   ! called by domain.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: domhgr.F90 6204 2016-01-04 13:47:06Z cetlod $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_hgr
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_hgr  ***
+      !!
+      !! ** Purpose :   Compute the geographical position (in degre) of the 
+      !!      model grid-points,  the horizontal scale factors (in meters) and 
+      !!      the Coriolis factor (in s-1).
+      !!
+      !! ** Method  :   The geographical position of the model grid-points is
+      !!      defined from analytical functions, fslam and fsphi, the deriva-
+      !!      tives of which gives the horizontal scale factors e1,e2.
+      !!      Defining two function fslam and fsphi and their derivatives in 
+      !!      the two horizontal directions (fse1 and fse2), the model grid-
+      !!      point position and scale factors are given by:
+      !!         t-point:
+      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
+      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
+      !!         u-point:
+      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
+      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
+      !!         v-point:
+      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
+      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
+      !!            f-point:
+      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
+      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
+      !!      Where fse1 and fse2 are defined by:
+      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
+      !!                                     +          di(fsphi) **2 )(i,j)
+      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
+      !!                                     +          dj(fsphi) **2 )(i,j)
+      !!
+      !!        The coriolis factor is given at z-point by:
+      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
+      !!
+      !!        This routine is given as an example, it must be modified
+      !!      following the user s desiderata. nevertheless, the output as
+      !!      well as the way to compute the model grid-point position and
+      !!      horizontal scale factors must be respected in order to insure
+      !!      second order accuracy schemes.
+      !!
+      !! N.B. If the domain is periodic, verify that scale factors are also
+      !!      periodic, and the coriolis term again.
+      !!
+      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-, 
+      !!                u-, v- and f-points (in degre)
+      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
+      !!               u-, v-  and f-points (in degre)
+      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
+      !!      scale factors (in meters) at t-, u-, v-, and f-points.
+      !!        define ff: coriolis factor at f-point
+      !!
+      !! References :   Marti, Madec and Delecluse, 1992, JGR
+      !!                Madec, Imbard, 1996, Clim. Dyn.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
+      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4)
+      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
+      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
+      REAL(wp) ::   zphi0, zbeta, znorme !
+      REAL(wp) ::   zarg, zf0, zminff, zmaxff
+      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
+      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
+      INTEGER  ::   isrow                ! index for ORCA1 starting row
+
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dom_hgr')
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
+         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
+         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
+         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
+         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
+         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
+         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m  
+         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m  
+      ENDIF
+
+
+      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
+
+      CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
+
+         CALL hgr_read           ! Defaultl option  :   NetCDF file
+
+         !                                                ! =====================
+         IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
+            !                                             ! =====================
+            IF( nn_cla == 0 ) THEN
+               !
+               ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km)
+               ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar    : e2u reduced to 20 km'
+               !
+               ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km)
+               ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3
+                                               e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km'
+               IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km'
+            ENDIF
+
+            ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u = 10 km)
+            ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km'
+            !
+         ENDIF
+
+            !                                             ! =====================
+         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
+            !                                             ! =====================
+            ! This dirty section will be suppressed by simplification process: all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            ! which had a grid-size of 362x292.
+            ! 
+            isrow = 332 - jpjglo
+            !
+            ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km)
+            ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km'
+
+            ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km)
+            ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km'
+
+            ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km)
+            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km'
+
+            ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
+            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km'
+
+            ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km)
+            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km'
+
+            ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km)
+            ij0 = 164 - isrow   ;   ij1 = 145 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km'
+
+            ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km)
+            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km'
+
+            ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km)
+            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km'
+            !
+            !
+         ENDIF
+
+         !                                                ! ======================
+         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
+            !                                             ! ======================
+            ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km)
+            ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait'
+            !
+            ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km)
+            ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait'
+            !
+            ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km)
+            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait'
+            !
+            ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km)
+            ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait'
+            !
+            ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km)
+            ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait'
+            !
+            ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km)
+            ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait'
+            !
+            !
+            ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km)
+            ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb'
+            !
+         ENDIF
+
+
+         ! N.B. :  General case, lat and long function of both i and j indices:
+         !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   &
+         !                                  + (                           fsdiph( zti, ztj ) )**2  )
+         !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   &
+         !                                  + (                           fsdiph( zui, zuj ) )**2  )
+         !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   &
+         !                                  + (                           fsdiph( zvi, zvj ) )**2  )
+         !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   &
+         !                                  + (                           fsdiph( zfi, zfj ) )**2  )
+         !
+         !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   &
+         !                                  + (                           fsdjph( zti, ztj ) )**2  )
+         !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   &
+         !                                  + (                           fsdjph( zui, zuj ) )**2  )
+         !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   &
+         !                                  + (                           fsdjph( zvi, zvj ) )**2  )
+         !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   &
+         !                                  + (                           fsdjph( zfi, zfj ) )**2  )
+
+
+      CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
+         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 )
+               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 )
+               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
+               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
+         ! Longitude
+               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
+               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
+               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
+               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
+         ! Latitude
+               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
+               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
+               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
+               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
+         ! e1
+               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
+               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
+               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
+               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
+         ! e2
+               e2t(ji,jj) = ra * rad * ppe2_deg
+               e2u(ji,jj) = ra * rad * ppe2_deg
+               e2v(ji,jj) = ra * rad * ppe2_deg
+               e2f(ji,jj) = ra * rad * ppe2_deg
+            END DO
+         END DO
+
+
+      CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
+         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
+
+         ! Position coordinates (in kilometers)
+         !                          ==========
+         glam0 = 0.e0
+         gphi0 = - ppe2_m * 1.e-3
+         
+#if defined key_agrif 
+         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
+            IF( .NOT. Agrif_Root() ) THEN
+              glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
+              gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
+              ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
+              ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()          
+            ENDIF
+         ENDIF
+#endif         
+         IF ( cp_cfg == 'vortex' ) THEN
+#if defined key_agrif 
+            IF( .NOT. Agrif_Root() ) THEN
+               ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
+               ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()
+            ENDIF
+#endif
+            glam0 = -(jpiglo-1)/2 * ppe1_m * 1.e-3
+            gphi0 = -(jpjglo-1)/2 * ppe2_m * 1.e-3
+         ENDIF 
+        
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       )
+               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
+               glamv(ji,jj) = glamt(ji,jj)
+               glamf(ji,jj) = glamu(ji,jj)
+   
+               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       )
+               gphiu(ji,jj) = gphit(ji,jj)
+               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
+               gphif(ji,jj) = gphiv(ji,jj)
+            END DO
+         END DO
+
+         ! Horizontal scale factors (in meters)
+         !                              ======
+         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
+         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
+         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
+         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
+
+      CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
+         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
+         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
+
+         !  Find index corresponding to the equator, given the grid spacing e1_deg
+         !  and the (approximate) southern latitude ppgphi0.
+         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
+         !  The formula should work even if the equator is outside the domain.
+         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
+         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
+         IF(  ppgphi0 > 0 )  ijeq = -ijeq
+
+         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 )
+               zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 )
+               zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
+               zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
+         ! Longitude
+               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
+               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
+               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
+               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
+         ! Latitude
+               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
+               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
+               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
+               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
+         ! e1
+               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
+               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
+               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
+               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
+         ! e2
+               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
+               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
+               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
+               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
+            END DO
+         END DO
+
+      CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
+         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
+
+         ! Position coordinates (in kilometers)
+         !                          ==========
+
+         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
+         zlam1 = -85
+         zphi1 = 29
+         ! resolution in meters
+         ze1 = 106000. / FLOAT(jp_cfg)            
+         ! benchmark: forced the resolution to be about 100 km
+         IF( nbench /= 0 )   ze1 = 106000.e0     
+         zsin_alpha = - SQRT( 2. ) / 2.
+         zcos_alpha =   SQRT( 2. ) / 2.
+         ze1deg = ze1 / (ra * rad)
+         IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon
+         !                                                          ! at the right jp_cfg resolution
+         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
+         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
+
+         IF( nprint==1 .AND. lwp )   THEN
+            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
+            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
+         ENDIF
+
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+             zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
+             zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
+
+             glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+             gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
+
+             glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
+             gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
+
+             glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
+             gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
+
+             glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+             gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
+           END DO
+          END DO
+
+         ! Horizontal scale factors (in meters)
+         !                              ======
+         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
+         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
+         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
+         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
+
+      CASE DEFAULT
+         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
+         CALL ctl_stop( ctmp1 )
+
+      END SELECT
+      
+      ! T-cell surface
+      ! --------------
+      e1e2t(:,:) = e1t(:,:) * e2t(:,:)
+    
+      ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning)
+      ! ---------------------------------------------------------------------------
+      e12t    (:,:) = e1t(:,:) * e2t(:,:)
+      e12u    (:,:) = e1u(:,:) * e2u(:,:)
+      e12v    (:,:) = e1v(:,:) * e2v(:,:)
+      e12f    (:,:) = e1f(:,:) * e2f(:,:)
+      r1_e12t (:,:) = 1._wp    / e12t(:,:)
+      r1_e12u (:,:) = 1._wp    / e12u(:,:)
+      r1_e12v (:,:) = 1._wp    / e12v(:,:)
+      r1_e12f (:,:) = 1._wp    / e12f(:,:)
+      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
+      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
+      r1_e1t  (:,:) = 1._wp    / e1t(:,:)
+      r1_e1u  (:,:) = 1._wp    / e1u(:,:)
+      r1_e1v  (:,:) = 1._wp    / e1v(:,:)
+      r1_e1f  (:,:) = 1._wp    / e1f(:,:)
+      r1_e2t  (:,:) = 1._wp    / e2t(:,:)
+      r1_e2u  (:,:) = 1._wp    / e2u(:,:)
+      r1_e2v  (:,:) = 1._wp    / e2v(:,:)
+      r1_e2f  (:,:) = 1._wp    / e2f(:,:)
+
+      ! Control printing : Grid informations (if not restart)
+      ! ----------------
+
+      IF( lwp .AND. .NOT.ln_rstart ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '          longitude and e1 scale factors'
+         WRITE(numout,*) '          ------------------------------'
+         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
+            glamv(ji,1), glamf(ji,1),   &
+            e1t(ji,1), e1u(ji,1),   &
+            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
+9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
+            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
+         
+         WRITE(numout,*)
+         WRITE(numout,*) '          latitude and e2 scale factors'
+         WRITE(numout,*) '          -----------------------------'
+         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
+            &                     gphiv(1,jj), gphif(1,jj),   &
+            &                     e2t  (1,jj), e2u  (1,jj),   &
+            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
+      ENDIF
+
+      
+      IF( nprint == 1 .AND. lwp ) THEN
+         WRITE(numout,*) '          e1u e2u '
+         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
+         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
+         WRITE(numout,*) '          e1v e2v  '
+         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
+         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
+         WRITE(numout,*) '          e1f e2f  '
+         CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
+         CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
+      ENDIF
+
+
+      ! ================= !
+      !  Coriolis factor  !
+      ! ================= !
+
+      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
+
+      CASE ( 0, 1, 4 )               ! mesh on the sphere
+
+         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
+
+      CASE ( 2 )                     ! f-plane at ppgphi0 
+
+         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
+
+         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
+
+      CASE ( 3 )                     ! beta-plane
+
+         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
+         zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
+         
+#if defined key_agrif
+         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
+            IF( .NOT. Agrif_Root() ) THEN
+              zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m)   & 
+                    &           / (ra * rad)
+            ENDIF
+         ENDIF
+#endif    
+
+         IF ( cp_cfg == 'vortex' ) THEN
+#if defined key_agrif 
+            IF( .NOT. Agrif_Root() ) THEN
+               ppgphi0=Agrif_Parent(ppgphi0)
+            ENDIF
+#endif
+            zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra
+            zphi0 = ppgphi0  
+         ENDIF 
+     
+         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
+
+         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
+         
+         IF(lwp) THEN
+            WRITE(numout,*) 
+            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
+            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
+         ENDIF
+         IF( lk_mpp ) THEN 
+            zminff=ff(nldi,nldj)
+            zmaxff=ff(nldi,nlej)
+            CALL mpp_min( zminff )   ! min over the global domain
+            CALL mpp_max( zmaxff )   ! max over the global domain
+            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
+         END IF
+
+      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
+
+         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
+         zphi0 = 15.e0                                                      ! latitude of the first row F-points
+         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
+
+         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
+
+         IF(lwp) THEN
+            WRITE(numout,*) 
+            WRITE(numout,*) '          Beta-plane and rotated domain : '
+            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
+         ENDIF
+
+         IF( lk_mpp ) THEN 
+            zminff=ff(nldi,nldj)
+            zmaxff=ff(nldi,nlej)
+            CALL mpp_min( zminff )   ! min over the global domain
+            CALL mpp_max( zmaxff )   ! max over the global domain
+            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
+         END IF
+
+      END SELECT
+
+
+      ! Control of domain for symetrical condition
+      ! ------------------------------------------
+      ! The equator line must be the latitude coordinate axe
+
+      IF( nperio == 2 ) THEN
+         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
+         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
+      !
+   END SUBROUTINE dom_hgr
+
+
+   SUBROUTINE hgr_read
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE hgr_read  ***
+      !!
+      !! ** Purpose :   Read a coordinate file in NetCDF format 
+      !!
+      !! ** Method  :   The mesh file has been defined trough a analytical 
+      !!      or semi-analytical method. It is read in a NetCDF file. 
+      !!     
+      !!----------------------------------------------------------------------
+      USE iom
+
+      INTEGER ::   inum   ! temporary logical unit
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
+         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
+      ENDIF
+      
+      CALL iom_open( 'coordinates', inum )
+      
+      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
+      
+      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
+      
+      CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr )
+      
+      CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr )
+      
+      CALL iom_close( inum )
+      
+! need to be define for the extended grid south of -80S
+! some point are undefined but you need to have e1 and e2 .NE. 0
+      WHERE (e1t==0.0_wp)
+         e1t=1.0e2
+      END WHERE
+      WHERE (e1v==0.0_wp)
+         e1v=1.0e2
+      END WHERE
+      WHERE (e1u==0.0_wp)
+         e1u=1.0e2
+      END WHERE
+      WHERE (e1f==0.0_wp)
+         e1f=1.0e2
+      END WHERE
+      WHERE (e2t==0.0_wp)
+         e2t=1.0e2
+      END WHERE
+      WHERE (e2v==0.0_wp)
+         e2v=1.0e2
+      END WHERE
+      WHERE (e2u==0.0_wp)
+         e2u=1.0e2
+      END WHERE
+      WHERE (e2f==0.0_wp)
+         e2f=1.0e2
+      END WHERE
+       
+    END SUBROUTINE hgr_read
+    
+   !!======================================================================
+END MODULE domhgr
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/domvvl.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/domvvl.F90	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/domvvl.F90	(revision 8011)
@@ -0,0 +1,1448 @@
+MODULE domvvl
+   !!======================================================================
+   !!                       ***  MODULE domvvl   ***
+   !! Ocean : 
+   !!======================================================================
+   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
+   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
+   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
+   !!                                          vvl option includes z_star and z_tilde coordinates
+   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
+   !!----------------------------------------------------------------------
+   !!   'key_vvl'                              variable volume
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
+   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
+   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
+   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
+   !!   dom_vvl_rst      : read/write restart file
+   !!   dom_vvl_ctl      : Check the vvl options
+   !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors 
+   !!                    : to account for manual changes to e[1,2][u,v] in some Straits 
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! ocean surface boundary condition
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+   USE restart         ! ocean restart
+   USE lib_mpp         ! distributed memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE wrk_nemo        ! Memory allocation
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC  dom_vvl_init       ! called by domain.F90
+   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
+   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
+   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
+   PRIVATE dom_vvl_orca_fix   ! called by dom_vvl_interpol
+
+   !!* Namelist nam_vvl
+   LOGICAL , PUBLIC                                      :: ln_vvl_zstar = .FALSE.              ! zstar  vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.             ! ztilde vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_layer = .FALSE.              ! level  vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_zstar_at_eqtor = .FALSE.     ! ztilde vertical coordinate
+   LOGICAL , PUBLIC                                      :: ln_vvl_kepe = .FALSE.               ! kinetic/potential energy transfer
+   !                                                                                            ! conservation: not used yet
+   REAL(wp)                                              :: rn_ahe3                   ! thickness diffusion coefficient
+   REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
+   REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
+   REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
+   LOGICAL , PUBLIC                                      :: ln_vvl_dbg = .FALSE.      ! debug control prints
+
+   !! * Module variables
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                       ! thickness diffusion transport
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 
+   !! $Id: domvvl.F90 7971 2017-04-26 10:03:26Z jchanut $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   INTEGER FUNCTION dom_vvl_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION dom_vvl_alloc  ***
+      !!----------------------------------------------------------------------
+      IF( ln_vvl_zstar ) dom_vvl_alloc = 0
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
+            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
+            &      STAT = dom_vvl_alloc        )
+         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
+         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
+         un_td = 0.0_wp
+         vn_td = 0.0_wp
+      ENDIF
+      IF( ln_vvl_ztilde ) THEN
+         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
+         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
+         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
+      ENDIF
+
+   END FUNCTION dom_vvl_alloc
+
+
+   SUBROUTINE dom_vvl_init
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_init  ***
+      !!                   
+      !! ** Purpose :  Initialization of all scale factors, depths
+      !!               and water column heights
+      !!
+      !! ** Method  :  - use restart file and/or initialize
+      !!               - interpolate scale factors
+      !!
+      !! ** Action  : - fse3t_(n/b) and tilde_e3t_(n/b)
+      !!              - Regrid: fse3(u/v)_n
+      !!                        fse3(u/v)_b       
+      !!                        fse3w_n           
+      !!                        fse3(u/v)w_b      
+      !!                        fse3(u/v)w_n      
+      !!                        fsdept_n, fsdepw_n and fsde3w_n
+      !!              - h(t/u/v)_0
+      !!              - frq_rst_e3t and frq_rst_hdv
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      USE phycst,  ONLY : rpi, rsmall, rad
+      !! * Local declarations
+      INTEGER ::   ji,jj,jk
+      INTEGER ::   ii0, ii1, ij0, ij1
+      REAL(wp)::   zcoef
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+
+      ! choose vertical coordinate (z_star, z_tilde or layer)
+      ! ==========================
+      CALL dom_vvl_ctl
+
+      ! Allocate module arrays
+      ! ======================
+      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
+
+      ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
+      ! ============================================================================
+      CALL dom_vvl_rst( nit000, 'READ' )
+      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
+
+      ! Reconstruction of all vertical scale factors at now and before time steps
+      ! =============================================================================
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+      ! Vertical scale factor interpolations
+      ! ------------------------------------
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
+      ! t- and w- points depth
+      ! ----------------------
+      ! set the isf depth as it is in the initial step
+      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+      fsdepw_n(:,:,1) = 0.0_wp
+      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
+      fsdepw_b(:,:,1) = 0.0_wp
+
+      DO jk = 2, jpk
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+                                                     ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
+                                                     ! 0.5 where jk = mikt  
+               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
+                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
+               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
+               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
+               fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  &
+                   &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk)) 
+            END DO
+         END DO
+      END DO
+
+      ! Before depth and Inverse of the local depth of the water column at u- and v- points
+      ! -----------------------------------------------------------------------------------
+      hu_b(:,:) = 0.
+      hv_b(:,:) = 0.
+      DO jk = 1, jpkm1
+         hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+      END DO
+      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
+      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
+
+      ! Restoring frequencies for z_tilde coordinate
+      ! ============================================
+      IF( ln_vvl_ztilde ) THEN
+         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
+         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
+         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
+         IF( ln_vvl_ztilde_as_zstar ) THEN
+            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
+            frq_rst_e3t(:,:) = 0.0_wp 
+            frq_rst_hdv(:,:) = 1.0_wp / rdt
+         ENDIF
+         IF ( ln_vvl_zstar_at_eqtor ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
+                     ! values outside the equatorial band and transition zone (ztilde)
+                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
+                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
+                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
+                     ! values inside the equatorial band (ztilde as zstar)
+                     frq_rst_e3t(ji,jj) =  0.0_wp
+                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
+                  ELSE
+                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
+                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
+                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+                        &                                          * 180._wp / 3.5_wp ) )
+                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
+                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
+                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+                        &                                          * 180._wp / 3.5_wp ) )
+                  ENDIF
+               END DO
+            END DO
+            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
+               ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
+               ij0 = 128   ;   ij1 = 135   ;   
+               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
+               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
+            ENDIF
+         ENDIF
+      ENDIF
+
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
+
+   END SUBROUTINE dom_vvl_init
+
+
+   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
+      !!                   
+      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
+      !!                 tranxt and dynspg routines
+      !!
+      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
+      !!               - z_tilde_case: after scale factor increment = 
+      !!                                    high frequency part of horizontal divergence
+      !!                                  + retsoring towards the background grid
+      !!                                  + thickness difusion
+      !!                               Then repartition of ssh INCREMENT proportionnaly
+      !!                               to the "baroclinic" level thickness.
+      !!
+      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
+      !!               - tilde_e3t_a: after increment of vertical scale factor 
+      !!                              in z_tilde case
+      !!               - fse3(t/u/v)_a
+      !!
+      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
+      !! * Arguments
+      INTEGER, INTENT( in )                  :: kt                    ! time step
+      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence
+      !! * Local declarations
+      INTEGER                                :: ji, jj, jk            ! dummy loop indices
+      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
+      REAL(wp)                               :: z2dt                  ! temporary scalars
+      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
+      LOGICAL                                :: ll_do_bclinic         ! temporary logical
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
+      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
+      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
+
+      IF(kt == nit000)   THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+      ENDIF
+
+      ll_do_bclinic = .TRUE.
+      IF( PRESENT(kcall) ) THEN
+         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
+      ENDIF
+
+      ! ******************************* !
+      ! After acale factors at t-points !
+      ! ******************************* !
+
+      !                                                ! --------------------------------------------- !
+                                                       ! z_star coordinate and barotropic z-tilde part !
+      !                                                ! --------------------------------------------- !
+
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+      DO jk = 1, jpkm1
+         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
+         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+      END DO
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
+         !                                                            ! ------baroclinic part------ !
+
+         ! I - initialization
+         ! ==================
+
+         ! 1 - barotropic divergence
+         ! -------------------------
+         zhdiv(:,:) = 0.
+         zht(:,:)   = 0.
+         DO jk = 1, jpkm1
+            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
+            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
+
+         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
+         ! --------------------------------------------------
+         IF( ln_vvl_ztilde ) THEN
+            IF( kt .GT. nit000 ) THEN
+               DO jk = 1, jpkm1
+                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
+                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+               END DO
+            ENDIF
+         END IF
+
+         ! II - after z_tilde increments of vertical scale factors
+         ! =======================================================
+         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
+
+         ! 1 - High frequency divergence term
+         ! ----------------------------------
+         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
+            DO jk = 1, jpkm1
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+            END DO
+         ELSE                         ! layer case
+            DO jk = 1, jpkm1
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+            END DO
+         END IF
+
+         ! 2 - Restoring term (z-tilde case only)
+         ! ------------------
+         IF( ln_vvl_ztilde ) THEN
+            DO jk = 1, jpk
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
+            END DO
+         END IF
+
+         ! 3 - Thickness diffusion term
+         ! ----------------------------
+         zwu(:,:) = 0.0_wp
+         zwv(:,:) = 0.0_wp
+         ! a - first derivative: diffusive fluxes
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
+                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         ! b - correction for last oceanic u-v points
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
+               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
+            END DO
+         END DO
+         ! c - second derivative: divergence of diffusive fluxes
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
+                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
+                     &                                            ) * r1_e12t(ji,jj)
+               END DO
+            END DO
+         END DO
+         ! d - thickness diffusion transport: boundary conditions
+         !     (stored for tracer advction and continuity equation)
+         CALL lbc_lnk( un_td , 'U' , -1._wp)
+         CALL lbc_lnk( vn_td , 'V' , -1._wp)
+
+         ! 4 - Time stepping of baroclinic scale factors
+         ! ---------------------------------------------
+         ! Leapfrog time stepping
+         ! ~~~~~~~~~~~~~~~~~~~~~~
+         IF( neuler == 0 .AND. kt == nit000 ) THEN
+            z2dt =  rdt
+         ELSE
+            z2dt = 2.0_wp * rdt
+         ENDIF
+         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
+         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
+
+         ! Maximum deformation control
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         ze3t(:,:,jpk) = 0.0_wp
+         DO jk = 1, jpkm1
+            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
+         END DO
+         z_tmax = MAXVAL( ze3t(:,:,:) )
+         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
+         z_tmin = MINVAL( ze3t(:,:,:) )
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         ! - ML - test: for the moment, stop simulation for too large e3_t variations
+         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
+            IF( lk_mpp ) THEN
+               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
+               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
+            ELSE
+               ijk_max = MAXLOC( ze3t(:,:,:) )
+               ijk_max(1) = ijk_max(1) + nimpp - 1
+               ijk_max(2) = ijk_max(2) + njmpp - 1
+               ijk_min = MINLOC( ze3t(:,:,:) )
+               ijk_min(1) = ijk_min(1) + nimpp - 1
+               ijk_min(2) = ijk_min(2) + njmpp - 1
+            ENDIF
+            IF (lwp) THEN
+               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
+               WRITE(numout, *) 'at i, j, k=', ijk_max
+               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
+               WRITE(numout, *) 'at i, j, k=', ijk_min            
+               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
+            ENDIF
+         ENDIF
+         ! - ML - end test
+         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
+         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
+         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
+
+         !
+         ! "tilda" change in the after scale factor
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         DO jk = 1, jpkm1
+            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
+         END DO
+         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
+         ! ==================================================================================
+         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
+         ! - ML - baroclinicity error should be better treated in the future
+         !        i.e. locally and not spread over the water column.
+         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
+         zht(:,:) = 0.
+         DO jk = 1, jpkm1
+            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+         DO jk = 1, jpkm1
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         END DO
+
+      ENDIF
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
+      !                                           ! ---baroclinic part--------- !
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
+      ENDIF
+
+      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
+         !
+         IF( lwp ) WRITE(numout, *) 'kt =', kt
+         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
+            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
+            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
+         END IF
+         !
+         zht(:,:) = 0.0_wp
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
+         !
+         zht(:,:) = 0.0_wp
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
+         !
+         zht(:,:) = 0.0_wp
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
+      END IF
+
+      ! *********************************** !
+      ! After scale factors at u- v- points !
+      ! *********************************** !
+
+      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
+      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
+
+      ! *********************************** !
+      ! After depths at u- v points         !
+      ! *********************************** !
+
+      hu_a(:,:) = 0._wp                        ! Ocean depth at U-points
+      hv_a(:,:) = 0._wp                        ! Ocean depth at V-points
+      DO jk = 1, jpkm1
+         hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
+         hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
+      END DO
+      !                                        ! Inverse of the local depth
+      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+
+      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
+
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
+
+   END SUBROUTINE dom_vvl_sf_nxt
+
+
+   SUBROUTINE dom_vvl_sf_swp( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_swp  ***
+      !!                   
+      !! ** Purpose :  compute time filter and swap of scale factors 
+      !!               compute all depths and related variables for next time step
+      !!               write outputs and restart file
+      !!
+      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
+      !!               - reconstruct scale factor at other grid points (interpolate)
+      !!               - recompute depths and water height fields
+      !!
+      !! ** Action  :  - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
+      !!               - Recompute:
+      !!                    fse3(u/v)_b       
+      !!                    fse3w_n           
+      !!                    fse3(u/v)w_b      
+      !!                    fse3(u/v)w_n      
+      !!                    fsdept_n, fsdepw_n  and fsde3w_n
+      !!                    h(u/v) and h(u/v)r
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in )               :: kt       ! time step
+      !! * Local declarations
+      INTEGER                             :: ji,jj,jk       ! dummy loop indices
+      REAL(wp)                            :: zcoef
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
+      !
+      IF( kt == nit000 )   THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
+      ENDIF
+      !
+      ! Time filter and swap of scale factors
+      ! =====================================
+      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         IF( neuler == 0 .AND. kt == nit000 ) THEN
+            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
+         ELSE
+            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
+            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+         ENDIF
+         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
+      ENDIF
+      fsdept_b(:,:,:) = fsdept_n(:,:,:)
+      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
+
+      fse3t_n(:,:,:) = fse3t_a(:,:,:)
+      fse3u_n(:,:,:) = fse3u_a(:,:,:)
+      fse3v_n(:,:,:) = fse3v_a(:,:,:)
+
+      ! Compute all missing vertical scale factor and depths
+      ! ====================================================
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
+      ! - JC - hu_b, hv_b, hur_b, hvr_b also
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
+      ! Vertical scale factor interpolations
+      ! ------------------------------------
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
+      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
+      ! t- and w- points depth
+      ! ----------------------
+      ! set the isf depth as it is in the initial step
+      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+      fsdepw_n(:,:,1) = 0.0_wp
+      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+
+      DO jk = 2, jpk
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+                                                                 ! 1 for jk = mikt
+               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
+                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
+               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
+            END DO
+         END DO
+      END DO
+
+      ! Local depth and Inverse of the local depth of the water column at u- and v- points
+      ! ----------------------------------------------------------------------------------
+      hu (:,:) = hu_a (:,:)
+      hv (:,:) = hv_a (:,:)
+
+      ! Inverse of the local depth
+      hur(:,:) = hur_a(:,:)
+      hvr(:,:) = hvr_a(:,:)
+
+      ! Local depth of the water column at t- points
+      ! --------------------------------------------
+      ht(:,:) = 0.
+      DO jk = 1, jpkm1
+         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+      END DO
+
+      ! write restart file
+      ! ==================
+      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
+
+   END SUBROUTINE dom_vvl_sf_swp
+
+
+   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl__interpol  ***
+      !!
+      !! ** Purpose :   interpolate scale factors from one grid point to another
+      !!
+      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
+      !!                - horizontal interpolation: grid cell surface averaging
+      !!                - vertical interpolation: simple averaging
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
+      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
+      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
+      LOGICAL ::   l_is_orca                                           ! local logical
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
+         !
+      l_is_orca = .FALSE.
+      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
+
+      SELECT CASE ( pout )
+         !               ! ------------------------------------- !
+      CASE( 'U' )        ! interpolation from T-point to U-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
+                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+               END DO
+            END DO
+         END DO
+         !
+         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'V' )        ! interpolation from T-point to V-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
+                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+               END DO
+            END DO
+         END DO
+         !
+         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'F' )        ! interpolation from U-point to F-point !
+         !               ! ------------------------------------- !
+         ! horizontal surface weighted interpolation
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
+                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+               END DO
+            END DO
+         END DO
+         !
+         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
+         ! boundary conditions
+         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
+         !               ! ------------------------------------- !
+      CASE( 'W' )        ! interpolation from T-point to W-point !
+         !               ! ------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
+               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
+         END DO
+         !               ! -------------------------------------- !
+      CASE( 'UW' )       ! interpolation from U-point to UW-point !
+         !               ! -------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
+               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
+         END DO
+         !               ! -------------------------------------- !
+      CASE( 'VW' )       ! interpolation from V-point to VW-point !
+         !               ! -------------------------------------- !
+         ! vertical simple interpolation
+         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
+         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+         DO jk = 2, jpk
+            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
+               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
+         END DO
+      END SELECT
+      !
+
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
+
+   END SUBROUTINE dom_vvl_interpol
+
+   SUBROUTINE dom_vvl_rst( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dom_vvl_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write VVL file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!                if the restart does not contain vertical scale factors,
+      !!                they are set to the _0 values
+      !!                if the restart does not contain vertical scale factors increments (z_tilde),
+      !!                they are set to 0.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !! * Local declarations
+      INTEGER ::   jk
+      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+         !                                   ! ===============
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            CALL rst_read_open                  !  open the restart file if necessary
+            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
+            !
+            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
+            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
+            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
+            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
+            !                             ! --------- !
+            !                             ! all cases !
+            !                             ! --------- !
+            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
+               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+               ! needed to restart if land processor not computed 
+               IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
+               WHERE ( tmask(:,:,:) == 0.0_wp ) 
+                  fse3t_n(:,:,:) = e3t_0(:,:,:)
+                  fse3t_b(:,:,:) = e3t_0(:,:,:)
+               END WHERE
+               IF( neuler == 0 ) THEN
+                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               ENDIF
+            ELSE IF( id1 > 0 ) THEN
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
+               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
+               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
+               fse3t_n(:,:,:) = fse3t_b(:,:,:)
+               neuler = 0
+            ELSE IF( id2 > 0 ) THEN
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
+               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
+               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+               fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               neuler = 0
+            ELSE
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
+               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
+               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               DO jk=1,jpk
+                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
+                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
+               END DO
+               fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               neuler = 0
+            ENDIF
+            !                             ! ----------- !
+            IF( ln_vvl_zstar ) THEN       ! z_star case !
+               !                          ! ----------- !
+               IF( MIN( id3, id4 ) > 0 ) THEN
+                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
+               ENDIF
+               !                          ! ----------------------- !
+            ELSE                          ! z_tilde and layer cases !
+               !                          ! ----------------------- !
+               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
+                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
+                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
+               ELSE                            ! one at least array is missing
+                  tilde_e3t_b(:,:,:) = 0.0_wp
+                  tilde_e3t_n(:,:,:) = 0.0_wp
+               ENDIF
+               !                          ! ------------ !
+               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
+                  !                       ! ------------ !
+                  IF( id5 > 0 ) THEN  ! required array exists
+                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
+                  ELSE                ! array is missing
+                     hdiv_lf(:,:,:) = 0.0_wp
+                  ENDIF
+               ENDIF
+            ENDIF
+            !
+         ELSE                                   !* Initialize at "rest" or from analytical solution
+            sshn(:,:) = 0.0_wp
+            IF ( cp_cfg == 'vortex' ) THEN
+               CALL istate_vortex_ssh
+            ENDIF
+            DO jk=1,jpk
+               fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+                   &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
+                   &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
+            END DO
+            fse3t_b(:,:,:) = fse3t_n(:,:,:)
+
+            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
+               tilde_e3t_b(:,:,:) = 0.0_wp
+               tilde_e3t_n(:,:,:) = 0.0_wp
+               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
+            END IF
+         ENDIF
+
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! ===================
+         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
+         !                                           ! --------- !
+         !                                           ! all cases !
+         !                                           ! --------- !
+         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
+         !                                           ! ----------------------- !
+         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
+            !                                        ! ----------------------- !
+            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
+         END IF
+         !                                           ! -------------!    
+         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
+            !                                        ! ------------ !
+            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
+         ENDIF
+
+      ENDIF
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
+
+   END SUBROUTINE dom_vvl_rst
+
+
+   SUBROUTINE dom_vvl_ctl
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_ctl  ***
+      !!                
+      !! ** Purpose :   Control the consistency between namelist options
+      !!                for vertical coordinate
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio
+      INTEGER ::   ios
+
+      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
+                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
+                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
+      !!---------------------------------------------------------------------- 
+
+      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist : 
+      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_vvl )
+
+      IF(lwp) THEN                    ! Namelist print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
+         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
+         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
+         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
+         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
+         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
+         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
+         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
+         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
+         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
+         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
+         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
+         IF( ln_vvl_ztilde_as_zstar ) THEN
+            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
+            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
+            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
+            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
+            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
+            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
+         ELSE
+            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
+            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
+            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
+            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
+         ENDIF
+         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
+         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
+      ENDIF
+
+      ioptio = 0                      ! Parameter control
+      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
+      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
+      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
+      IF( ln_vvl_layer           )        ioptio = ioptio + 1
+
+      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
+      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
+
+      IF(lwp) THEN                   ! Print the choice
+         WRITE(numout,*)
+         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
+         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
+         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
+         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
+         ! - ML - Option not developed yet
+         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
+         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
+      ENDIF
+
+#if defined key_agrif
+      IF ((.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar)) CALL ctl_stop( 'AGRIF implemented with zstar coordinate only (key_vvl)' )
+#endif
+
+   END SUBROUTINE dom_vvl_ctl
+
+   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
+      !!                     
+      !! ** Purpose :   Correct surface weighted, horizontally interpolated, 
+      !!                scale factors at locations that have been individually
+      !!                modified in domhgr. Such modifications break the
+      !!                relationship between e12t and e1u*e2u etc.
+      !!                Recompute some scale factors ignoring the modified metric.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
+      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
+      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
+      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
+      INTEGER ::   isrow                                               ! index for ORCA1 starting row
+      !! acc
+      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
+      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
+      !! 
+      !                                                ! =====================
+      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
+         !                                             ! =====================
+      !! acc
+         IF( nn_cla == 0 ) THEN
+            !
+            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
+            ij0 = 102   ;   ij1 = 102
+            DO jk = 1, jpkm1
+               DO jj = mj0(ij0), mj1(ij1)
+                  DO ji = mi0(ii0), mi1(ii1)
+                     SELECT CASE ( pout )
+                     CASE( 'U' )
+                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                     CASE( 'F' )
+                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                     END SELECT
+                  END DO
+               END DO
+            END DO
+            !
+            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
+            ij0 =  88   ;   ij1 =  88
+            DO jk = 1, jpkm1
+               DO jj = mj0(ij0), mj1(ij1)
+                  DO ji = mi0(ii0), mi1(ii1)
+                     SELECT CASE ( pout )
+                     CASE( 'U' )
+                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                     CASE( 'V' )
+                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                     CASE( 'F' )
+                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                     END SELECT
+                  END DO
+               END DO
+            END DO
+         ENDIF
+
+         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
+         ij0 = 116   ;   ij1 = 116
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+         !                                             ! =====================
+      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
+         !                                             ! =====================
+         ! This dirty section will be suppressed by simplification process:
+         ! all this will come back in input files
+         ! Currently these hard-wired indices relate to configuration with
+         ! extend grid (jpjglo=332)
+         ! which had a grid-size of 362x292.
+         isrow = 332 - jpjglo
+         !
+         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified)
+         ij0 = 241 - isrow   ;   ij1 = 241 - isrow
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
+         ij0 = 248 - isrow   ;   ij1 = 248 - isrow
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &  
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &  
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified)
+         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
+         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified)
+         ij0 = 164 - isrow  ;   ij1 = 165  - isrow  
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified)
+         ij0 = 164 - isrow    ;   ij1 = 165  - isrow  
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
+         ij0 = 181 - isrow    ;   ij1 = 182 - isrow  
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
+         ij0 = 181 - isrow    ;   ij1 = 182 - isrow  
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+      ENDIF
+         !                                             ! =====================
+      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
+         !                                             ! =====================
+         !
+         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
+         ij0 = 327   ;   ij1 = 327
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
+         ij0 = 343   ;   ij1 = 343
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &  
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &  
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
+         ij0 = 232   ;   ij1 = 232
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
+         ij0 = 232   ;   ij1 = 232
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
+         ij0 = 270   ;   ij1 = 270
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'U' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
+                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
+                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
+                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
+                  CASE( 'F' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
+                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
+                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
+                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
+         ij0 = 232   ;   ij1 = 233
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+         !
+         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
+         ij0 = 276   ;   ij1 = 276
+         DO jk = 1, jpkm1
+            DO jj = mj0(ij0), mj1(ij1)
+               DO ji = mi0(ii0), mi1(ii1)
+                  SELECT CASE ( pout )
+                  CASE( 'V' )
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
+                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
+                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
+                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
+                  END SELECT
+               END DO
+            END DO
+         END DO
+      ENDIF
+   END SUBROUTINE dom_vvl_orca_fix
+
+   SUBROUTINE istate_vortex_ssh
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE istate_vortex_ssh  ***
+      !! 
+      !! ** Purpose :   Initialization of baroclinic vortex ssh
+      !!
+      !! ** Method  : - set ssh
+      !!----------------------------------------------------------------------
+      INTEGER :: ji, jj, jk  ! dummy loop indices
+      REAL(wp) :: zx, zy, zP0, zumax, zlambda, zf0, zH, zrho1, za
+      !!----------------------------------------------------------------------
+      !
+      !
+      zf0   = 2._wp * omega * SIN( rad * ppgphi0 )
+      zumax = 1._wp
+      zlambda = SQRT(2._wp)*60.e3
+      zH = 0.5_wp * rn_bathy
+      !
+      zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
+      ! Sea level:
+      za = -zP0 * (1._wp-EXP(-zH)) / (grav*(zH-1._wp + exp(-zH)))
+      DO ji=1, jpi
+         DO jj=1, jpj
+            zx = glamt(ji,jj) * 1.e3
+            zy = gphit(ji,jj) * 1.e3
+            zrho1 = rau0 + za * exp(-(zx**2+zy**2)/zlambda**2)
+            sshn(ji,jj) = zP0 * exp(-(zx**2+zy**2)/zlambda**2)/(zrho1*grav) * tmask(ji,jj,1)
+         END DO
+      END DO
+      sshb(:,:) = sshn(:,:)
+      !
+   END SUBROUTINE istate_vortex_ssh
+
+   !!======================================================================
+END MODULE domvvl
+
+
+
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/istate.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/istate.F90	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/istate.F90	(revision 8011)
@@ -0,0 +1,676 @@
+MODULE istate
+   !!======================================================================
+   !!                     ***  MODULE  istate  ***
+   !! Ocean state   :  initial state setting
+   !!=====================================================================
+   !! History :  OPA  !  1989-12  (P. Andrich)  Original code
+   !!            5.0  !  1991-11  (G. Madec)  rewritting
+   !!            6.0  !  1996-01  (G. Madec)  terrain following coordinates
+   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_eel
+   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_uvg
+   !!   NEMO     1.0  !  2003-08  (G. Madec, C. Talandier)  F90: Free form, modules + EEL R5
+   !!             -   !  2004-05  (A. Koch-Larrouy)  istate_gyre 
+   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
+   !!            3.3  !  2010-10  (C. Ethe) merge TRC-TRA
+   !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   istate_init   : initial state setting
+   !!   istate_tem    : analytical profile for initial Temperature
+   !!   istate_sal    : analytical profile for initial Salinity
+   !!   istate_eel    : initial state setting of EEL R5 configuration
+   !!   istate_gyre   : initial state setting of GYRE configuration
+   !!   istate_uvg    : initial velocity in geostropic balance
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers 
+   USE dom_oce         ! ocean space and time domain 
+   USE c1d             ! 1D vertical configuration
+   USE daymod          ! calendar
+   USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
+   USE ldftra_oce      ! ocean active tracers: lateral physics
+   USE zdf_oce         ! ocean vertical physics
+   USE phycst          ! physical constants
+   USE dtatsd          ! data temperature and salinity   (dta_tsd routine)
+   USE dtauvd          ! data: U & V current             (dta_uvd routine)
+   USE zpshde          ! partial step: hor. derivative (zps_hde routine)
+   USE eosbn2          ! equation of state            (eos bn2 routine)
+   USE domvvl          ! varying vertical mesh
+   USE dynspg_oce      ! pressure gradient schemes
+   USE dynspg_flt      ! filtered free surface
+   USE sol_oce         ! ocean solver variables
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O library
+   USE lib_mpp         ! MPP library
+   USE restart         ! restart
+   USE wrk_nemo        ! Memory allocation
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   istate_init ! routine called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! $Id: istate.F90 5332 2015-06-01 16:59:09Z mathiot $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE istate_init
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE istate_init  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracer fields.
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('istate_init')
+      !
+
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+
+      CALL dta_tsd_init                       ! Initialisation of T & S input data
+      IF( lk_c1d ) CALL dta_uvd_init          ! Initialization of U & V input data
+
+      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
+      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
+      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
+      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
+
+      IF( ln_rstart ) THEN                    ! Restart from a file
+         !                                    ! -------------------
+         CALL rst_read                           ! Read the restart file
+         CALL day_init                           ! model calendar (using both namelist and restart infos)
+      ELSE
+         !                                    ! Start from rest
+         !                                    ! ---------------
+         numror = 0                              ! define numror = 0 -> no restart file to read
+         neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
+         CALL day_init                           ! model calendar (using both namelist and restart infos)
+         !                                       ! Initialization of ocean to zero
+         !   before fields      !       now fields     
+         sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
+         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
+         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp  
+         rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp
+         hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp
+         !
+         IF( cp_cfg == 'eel' ) THEN
+            CALL istate_eel                      ! EEL   configuration : start from pre-defined U,V T-S fields
+         ELSEIF( cp_cfg == 'gyre' ) THEN         
+            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
+         ELSEIF( cp_cfg == 'vortex' ) THEN         
+            CALL istate_vortex                   ! GYRE  configuration : start from pre-defined T-S fields
+         ELSE                                    ! Initial T-S, U-V fields read in files
+            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
+               CALL dta_tsd( nit000, tsb )  
+               tsn(:,:,:,:) = tsb(:,:,:,:)
+               !
+            ELSE                                 ! Initial T-S fields defined analytically
+               CALL istate_t_s
+            ENDIF
+            IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
+               CALL wrk_alloc( jpi, jpj, jpk, 2, zuvd )
+               CALL dta_uvd( nit000, zuvd )
+               ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
+               vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
+               CALL wrk_dealloc( jpi, jpj, jpk, 2, zuvd )
+            ENDIF
+         ENDIF
+         !
+         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities
+#if ! defined key_c1d
+         IF( ln_zps .AND. .NOT. ln_isfcav)                                 &
+            &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+         IF( ln_zps .AND.       ln_isfcav)                                 &
+            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
+            &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+#endif
+         !   
+         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpk
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+            ENDDO
+         ENDIF
+         ! 
+      ENDIF
+      !
+      IF( lk_agrif ) THEN                  ! read free surface arrays in restart file
+         IF( ln_rstart ) THEN
+            IF( lk_dynspg_flt )  THEN      ! read or initialize the following fields
+               !                           ! gcx, gcxb for agrif_opa_init
+               IF( sol_oce_alloc()  > 0 )   CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')
+               CALL flt_rst( nit000, 'READ' )
+            ENDIF
+         ENDIF                             ! explicit case not coded yet with AGRIF
+      ENDIF
+      !
+      ! 
+      ! Initialize "now" and "before" barotropic velocities:
+      ! Do it whatever the free surface method, these arrays
+      ! being eventually used
+      !
+      !
+      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
+      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
+      !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               !
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      un_b(:,:) = un_b(:,:) * hur  (:,:)
+      vn_b(:,:) = vn_b(:,:) * hvr  (:,:)
+      !
+      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
+      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
+      !
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('istate_init')
+      !
+   END SUBROUTINE istate_init
+
+
+   SUBROUTINE istate_t_s
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE istate_t_s  ***
+      !!   
+      !! ** Purpose :   Intialization of the temperature field with an 
+      !!      analytical profile or a file (i.e. in EEL configuration)
+      !!
+      !! ** Method  : - temperature: use Philander analytic profile
+      !!              - salinity   : use to a constant value 35.5
+      !!
+      !! References :  Philander ???
+      !!----------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk
+      REAL(wp) ::   zsal = 35.50
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~   and constant salinity (',zsal,' psu)'
+      !
+      DO jk = 1, jpk
+         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   &
+            &                + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.)  ) * tmask(:,:,jk)
+         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
+      END DO
+      tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
+      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
+      !
+   END SUBROUTINE istate_t_s
+
+
+   SUBROUTINE istate_eel
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE istate_eel  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracers for EEL R5
+      !!      configuration (channel with or without a topographic bump)
+      !!
+      !! ** Method  : - set temprature field
+      !!              - set salinity field
+      !!              - set velocity field including horizontal divergence
+      !!                and relative vorticity fields
+      !!----------------------------------------------------------------------
+      USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine)
+      USE iom
+      !
+      INTEGER  ::   inum              ! temporary logical unit
+      INTEGER  ::   ji, jj, jk        ! dummy loop indices
+      INTEGER  ::   ijloc
+      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
+      REAL(wp) ::   zt1  = 15._wp                   ! surface temperature value (EEL R5)
+      REAL(wp) ::   zt2  =  5._wp                   ! bottom  temperature value (EEL R5)
+      REAL(wp) ::   zsal = 35.0_wp                  ! constant salinity (EEL R2, R5 and R6)
+      REAL(wp) ::   zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
+      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( jp_cfg ) 
+         !                                              ! ====================
+         CASE ( 5 )                                     ! EEL R5 configuration
+            !                                           ! ====================
+            !
+            ! set temperature field with a linear profile
+            ! -------------------------------------------
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+            !
+            zh1 = gdept_1d(  1  )
+            zh2 = gdept_1d(jpkm1)
+            !
+            zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
+            zcst   = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
+            !
+            DO jk = 1, jpk
+               tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
+               tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
+            END DO
+            !
+            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
+               &                             1     , jpi   , 5     , 1     , jpk   ,   &
+               &                             1     , 1.    , numout                  )
+            !
+            ! set salinity field to a constant value
+            ! --------------------------------------
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+            !
+            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
+            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
+            !
+            ! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
+            ! ----------------
+            ! Start EEL5 configuration with barotropic geostrophic velocities 
+            ! according the sshb and sshn SSH imposed.
+            ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
+            ! we use the Coriolis frequency at mid-channel.   
+            ub(:,:,:) = zueel * umask(:,:,:)
+            un(:,:,:) = ub(:,:,:)
+            ijloc = mj0(INT(jpjglo-1)/2)
+            zfcor = ff(1,ijloc)
+            !
+            DO jj = 1, jpjglo
+               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 
+            END DO
+            !
+            IF(lwp) THEN
+               WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
+               WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
+               WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
+            ENDIF
+            !
+            DO jj = 1, nlcj
+               DO ji = 1, nlci
+                  sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
+               END DO
+            END DO
+            sshb(nlci+1:jpi,      :   ) = 0.e0      ! set to zero extra mpp columns
+            sshb(      :   ,nlcj+1:jpj) = 0.e0      ! set to zero extra mpp rows
+            !
+            sshn(:,:) = sshb(:,:)                   ! set now ssh to the before value
+            !
+            IF( nn_rstssh /= 0 ) THEN  
+               nn_rstssh = 0                        ! hand-made initilization of ssh 
+               CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
+            ENDIF
+            !
+            CALL div_cur( nit000 )                  ! horizontal divergence and relative vorticity (curl)
+            ! N.B. the vertical velocity will be computed from the horizontal divergence field
+            ! in istate by a call to wzv routine
+
+
+            !                                     ! ==========================
+         CASE ( 2 , 6 )                           ! EEL R2 or R6 configuration
+            !                                     ! ==========================
+            !
+            ! set temperature field with a NetCDF file
+            ! ----------------------------------------
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+            !
+            CALL iom_open ( 'eel.initemp', inum )
+            CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
+            CALL iom_close( inum )
+            !
+            tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem)                            ! set nox temperature to tb
+            !
+            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
+               &                            1     , jpi   , 5     , 1     , jpk   ,   &
+               &                            1     , 1.    , numout                  )
+            !
+            ! set salinity field to a constant value
+            ! --------------------------------------
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+            !
+            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
+            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
+            !
+            !                                    ! ===========================
+         CASE DEFAULT                            ! NONE existing configuration
+            !                                    ! ===========================
+            WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
+            CALL ctl_stop( ctmp1 )
+            !
+      END SELECT
+      !
+   END SUBROUTINE istate_eel
+
+
+   SUBROUTINE istate_gyre
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE istate_gyre  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracers for GYRE
+      !!      configuration (double gyre with rotated domain)
+      !!
+      !! ** Method  : - set temprature field
+      !!              - set salinity field
+      !!----------------------------------------------------------------------
+      INTEGER :: ji, jj, jk  ! dummy loop indices
+      INTEGER            ::   inum          ! temporary logical unit
+      INTEGER, PARAMETER ::   ntsinit = 0   ! (0/1) (analytical/input data files) T&S initialization
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( ntsinit)
+      !
+      CASE ( 0 )                  ! analytical T/S profil deduced from LEVITUS
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+         !
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  tsn(ji,jj,jk,jp_tem) = (  16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 )         )   &
+                       &           * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2               &
+                       &       + (      15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) )       &
+                       &                - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.)               &    
+                       &                + 7.  * (1500. - fsdept(ji,jj,jk)) / 1500.             )   & 
+                       &           * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
+                  tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
+                  tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
+
+                  tsn(ji,jj,jk,jp_sal) =  (  36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 )  )  &
+                     &              * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2          &
+                     &          + (  35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000.         &
+                     &                - 1.62 * TANH( (fsdept(ji,jj,jk) - 60.  ) / 650. )       &
+                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 35.  ) / 100. )       &
+                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.)    )  &
+                     &              * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 
+                  tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
+                  tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 1 )                  ! T/S data fields read in dta_tem.nc/data_sal.nc files
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+         IF(lwp) WRITE(numout,*) '              NetCDF FORMAT'
+
+         ! Read temperature field
+         ! ----------------------
+         CALL iom_open ( 'data_tem', inum )
+         CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 
+         CALL iom_close( inum )
+
+         tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 
+         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
+
+         ! Read salinity field
+         ! -------------------
+         CALL iom_open ( 'data_sal', inum )
+         CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 
+         CALL iom_close( inum )
+
+         tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
+         tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
+         !
+      END SELECT
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '              Initial temperature and salinity profiles:'
+         WRITE(numout, "(9x,' level   gdept_1d   temperature   salinity   ')" )
+         WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
+      ENDIF
+      !
+   END SUBROUTINE istate_gyre
+
+   SUBROUTINE istate_vortex
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE istate_vortex  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracers for 
+      !!                baroclinic vortex
+      !!
+      !! ** Method  : - set ssh
+      !!              - set temperature/salinity fields
+      !!              - set velocities
+      !!
+      !! ** Ref     : Penven et al. (2006) : Evaluation and application of 
+      !!              the ROMS 1-way embedding procedure to the central 
+      !!              california upwelling system. Ocean Mod., 12, 1–2, 157–187.
+      !!----------------------------------------------------------------------
+      INTEGER :: ji, jj, jk  ! dummy loop indices
+      REAL(wp) :: zx, zy, zP0, zumax, zlambda, zn2, zf0, zH, zrho1, za, zf
+      REAL(wp) :: zdt, zdu, zdv
+      !!----------------------------------------------------------------------
+      !
+      !
+      zf0   = 2._wp * omega * SIN( rad * ppgphi0 )
+      zumax = 1._wp ! Anticyclonic: set zumax=-1 for cyclonic
+      zlambda = SQRT(2._wp)*60.e3
+      zn2 = 3.e-3**2
+      zH = 0.5_wp * rn_bathy
+      !
+      zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
+      ! Sea level:
+      za = -zP0 * (1._wp-EXP(-zH)) / (grav*(zH-1._wp + EXP(-zH)))
+      DO ji=1, jpi
+         DO jj=1, jpj
+            zx = glamt(ji,jj) * 1.e3
+            zy = gphit(ji,jj) * 1.e3
+            zrho1 = rau0 + za * EXP(-(zx**2+zy**2)/zlambda**2)
+            sshn(ji,jj) = zP0 * EXP(-(zx**2+zy**2)/zlambda**2)/(zrho1*grav) * tmask(ji,jj,1)
+         END DO
+      END DO
+      sshb(:,:) = sshn(:,:)
+      !
+      ! temperature:         
+      DO ji=1, jpi
+         DO jj=1, jpj
+            zx = glamt(ji,jj) * 1.e3
+            zy = gphit(ji,jj) * 1.e3
+            DO jk=1,jpk
+               zdt =  fsdept(ji,jj,jk) 
+               zrho1 = rau0 * (1._wp + zn2*zdt/grav)
+               IF (zdt < zH) THEN
+                  zrho1 = zrho1 - zP0 * (1._wp-EXP(zdt-zH)) &
+                          & * EXP(-(zx**2+zy**2)/zlambda**2) / (grav*(zH -1._wp + exp(-zH)));
+               ENDIF
+               tsn(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
+      !
+      ! salinity:  
+      tsn(:,:,:,jp_sal) = 35._wp 
+      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
+      !
+      ! velocities:
+      za = 2._wp * zP0 / (zf0 * rau0 * zlambda**2)
+      DO ji=1, jpim1
+         DO jj=1, jpj
+            zx = glamu(ji,jj) * 1.e3
+            zy = gphiu(ji,jj) * 1.e3
+            DO jk=1, jpk
+               zdu = 0.5_wp * (fsdept(ji  ,jj,jk) + fsdept(ji+1,jj,jk))
+               zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
+               IF (zdu < zH) THEN
+                  un(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * umask(ji,jj,jk)
+               ELSE
+                  un(ji,jj,jk) = 0._wp
+               ENDIF
+            END DO
+         END DO
+      END DO
+      !
+      DO ji=1, jpi
+         DO jj=1, jpjm1
+            zx = glamv(ji,jj) * 1.e3
+            zy = gphiv(ji,jj) * 1.e3
+            DO jk=1, jpk
+               zdv = 0.5_wp * (fsdept(ji  ,jj,jk) + fsdept(ji,jj+1,jk))
+               zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
+               IF (zdv < zH) THEN
+                  vn(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * umask(ji,jj,jk)
+               ELSE
+                  vn(ji,jj,jk) = 0._wp
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk( un, 'U', -1. )
+      CALL lbc_lnk( vn, 'V', -1. )
+
+      ub(:,:,:) = un(:,:,:)
+      vb(:,:,:) = vn(:,:,:)
+
+   END SUBROUTINE istate_vortex
+
+   SUBROUTINE istate_uvg
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE istate_uvg  ***
+      !!
+      !! ** Purpose :   Compute the geostrophic velocities from (tn,sn) fields
+      !!
+      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic 
+      !!      pressure is computed by integrating the in-situ density from the
+      !!      surface to the bottom.
+      !!                 p=integral [ rau*g dz ]
+      !!----------------------------------------------------------------------
+      USE dynspg          ! surface pressure gradient             (dyn_spg routine)
+      USE divcur          ! hor. divergence & rel. vorticity      (div_cur routine)
+      USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+      !
+      INTEGER ::   ji, jj, jk        ! dummy loop indices
+      INTEGER ::   indic             ! ???
+      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zprn)
+      !
+      IF(lwp) WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+
+      ! Compute the now hydrostatic pressure
+      ! ------------------------------------
+
+      zalfg = 0.5 * grav * rau0
+      
+      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
+
+      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
+         zprn(:,:,jk) = zprn(:,:,jk-1)   &
+            &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
+      END DO  
+
+      ! Compute geostrophic balance
+      ! ---------------------------
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vertor opt.
+               zmsv = 1. / MAX(  umask(ji-1,jj+1,jk) + umask(ji  ,jj+1,jk)   &
+                               + umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) , 1.  )
+               zphv = ( zprn(ji  ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1)   &
+                    + ( zprn(ji+1,jj+1,jk) - zprn(ji  ,jj+1,jk) ) * umask(ji  ,jj+1,jk) / e1u(ji  ,jj+1)   &
+                    + ( zprn(ji  ,jj  ,jk) - zprn(ji-1,jj  ,jk) ) * umask(ji-1,jj  ,jk) / e1u(ji-1,jj  )   &
+                    + ( zprn(ji+1,jj  ,jk) - zprn(ji  ,jj  ,jk) ) * umask(ji  ,jj  ,jk) / e1u(ji  ,jj  )
+               zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
+
+               zmsu = 1. / MAX(  vmask(ji+1,jj  ,jk) + vmask(ji  ,jj  ,jk)   &
+                               + vmask(ji+1,jj-1,jk) + vmask(ji  ,jj-1,jk) , 1.  )
+               zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj  ,jk) ) * vmask(ji+1,jj  ,jk) / e2v(ji+1,jj  )   &
+                    + ( zprn(ji  ,jj+1,jk) - zprn(ji  ,jj  ,jk) ) * vmask(ji  ,jj  ,jk) / e2v(ji  ,jj  )   &
+                    + ( zprn(ji+1,jj  ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1)   &
+                    + ( zprn(ji  ,jj  ,jk) - zprn(ji  ,jj-1,jk) ) * vmask(ji  ,jj-1,jk) / e2v(ji  ,jj-1)
+               zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
+
+               ! Compute the geostrophic velocities
+               un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji  ,jj-1) )
+               vn(ji,jj,jk) =  2. * zphv / ( ff(ji,jj) + ff(ji-1,jj  ) )
+            END DO
+         END DO
+      END DO
+
+      IF(lwp) WRITE(numout,*) '         we force to zero bottom velocity'
+
+      ! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
+      ! to have a zero bottom velocity
+
+      DO jk = 1, jpkm1
+         un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
+         vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
+      END DO
+
+      CALL lbc_lnk( un, 'U', -1. )
+      CALL lbc_lnk( vn, 'V', -1. )
+      
+      ub(:,:,:) = un(:,:,:)
+      vb(:,:,:) = vn(:,:,:)
+      
+      ! WARNING !!!!!
+      ! after initializing u and v, we need to calculate the initial streamfunction bsf.
+      ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
+      ! to do that, we call dyn_spg with a special trick:
+      ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
+      ! right value assuming the velocities have been set up in one time step.
+      ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
+      !  sets up s false trend to calculate the barotropic streamfunction.
+
+      ua(:,:,:) = ub(:,:,:) / rdt
+      va(:,:,:) = vb(:,:,:) / rdt
+
+      ! calls dyn_spg. we assume euler time step, starting from rest.
+      indic = 0
+      CALL dyn_spg( nit000, indic )       ! surface pressure gradient
+
+      ! the new velocity is ua*rdt
+
+      CALL lbc_lnk( ua, 'U', -1. )
+      CALL lbc_lnk( va, 'V', -1. )
+
+      ub(:,:,:) = ua(:,:,:) * rdt
+      vb(:,:,:) = va(:,:,:) * rdt
+      ua(:,:,:) = 0.e0
+      va(:,:,:) = 0.e0
+      un(:,:,:) = ub(:,:,:)
+      vn(:,:,:) = vb(:,:,:)
+       
+      ! Compute the divergence and curl
+
+      CALL div_cur( nit000 )            ! now horizontal divergence and curl
+
+      hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value
+      rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zprn)
+      !
+   END SUBROUTINE istate_uvg
+
+   !!=====================================================================
+END MODULE istate
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/cpp_VORTEX.fcm
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/cpp_VORTEX.fcm	(revision 8011)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/cpp_VORTEX.fcm	(revision 8011)
@@ -0,0 +1,1 @@
+bld::tool::fppkeys key_dynspg_ts key_bdy key_zdfcst key_iomput key_mpp_mpi key_vvl key_agrif
