Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg	(revision 10001)
@@ -23,5 +23,5 @@
    nn_it000    =       1   !  first time step
    nn_itend    =      75   !  last  time step (std 5475)
-   nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_istate   = .false.   !  output the initial state
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg	(revision 10001)
@@ -23,5 +23,5 @@
    nn_it000    =    1      !  first time step
    nn_itend    =  300      !  last  time step
-   nn_istate   =      0    !  output the initial state (1) or not (0)
+   ln_istate   = .false.   !  output the initial state
    ln_clobber  = .true.    !  clobber (overwrite) an existing file
 /
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg	(revision 10001)
@@ -23,5 +23,5 @@
    nn_it000    =    1      !  first time step
    nn_itend    =  900      !  last  time step
-   nn_istate   =      0    !  output the initial state (1) or not (0)
+   ln_istate   = .false.   !  output the initial state
    ln_clobber  = .true.    !  clobber (overwrite) an existing file
 /
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg	(revision 10001)
@@ -23,5 +23,5 @@
    nn_it000    =       1   !  first time step
    nn_itend    =      75   !  last  time step (std 5475)
-   nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_istate   = .false.   !  output the initial state
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_BFM/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_BFM/EXPREF/namelist_cfg	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_BFM/EXPREF/namelist_cfg	(revision 10001)
@@ -26,5 +26,5 @@
    nn_stock    =    4320   !  frequency of creation of a restart file (modulo referenced to 1)
    nn_write    =      60   !  frequency of write in the output file   (modulo referenced to nn_it000)
-   nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_istate   = .false.   !  output the initial state
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_PISCES/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_PISCES/EXPREF/namelist_cfg	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_PISCES/EXPREF/namelist_cfg	(revision 10001)
@@ -157,11 +157,11 @@
 &namdyn_adv    !   formulation of the momentum advection                (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynadv_vec = .true.  !  vector form - 2nd centered scheme
-     nn_dynkeg     = 0        ! grad(KE) scheme: =0   C2  ;  =1   Hollingsworth correction
+   ln_dynadv_cen2= .true.  !  flux form - 2nd order centered scheme
+   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
 /
 !-----------------------------------------------------------------------
 &namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynvor_ene = .true.  !  energy conserving scheme
+   ln_dynvor_enT = .true.  !  energy conserving scheme (T-point)
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 10001)
@@ -23,5 +23,5 @@
    nn_it000    =       1   !  first time step
    nn_itend    =    5475   !  last  time step (std 5475)
-   nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_istate   = .false.   !  output the initial state
 /
 !-----------------------------------------------------------------------
@@ -295,12 +295,11 @@
 &namdyn_adv    !   formulation of the momentum advection                (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynadv_vec = .true.  !  vector form - 2nd centered scheme
-     nn_dynkeg     = 0        ! grad(KE) scheme: =0   C2  ;  =1   Hollingsworth correction
+   ln_dynadv_cen2= .true.  !  flux form - 2nd order centered scheme
+   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
 /
 !-----------------------------------------------------------------------
 &namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_dynvor_een = .true.  !  energy & enstrophy scheme
-      nn_een_e3f = 0          ! =0   e3f = mean masked e3t divided by 4
+   ln_dynvor_enT = .true.  !  energy conserving scheme (T-point)
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/SHARED/namelist_ref	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/SHARED/namelist_ref	(revision 10001)
@@ -52,5 +52,5 @@
       cn_ocerst_outdir= "."         !  directory in which to write output ocean restarts
    ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model
-   nn_istate   =       0   !  output the initial state (1) or not (0)
+   ln_istate   = .false.   !  output the initial state
    ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
    nn_stock    =    5475   !  frequency of creation of a restart file (modulo referenced to 1)
@@ -67,9 +67,13 @@
 &namdom        !   time and space domain
 !-----------------------------------------------------------------------
+   ln_RK3      = .false.   !  3rd order Runge-Kutta time stepping scheme
+   ln_MLF      = .false.   !  Modified Leap-Frog    time stepping scheme
+      rn_atfp     =    0.1    !  asselin time filter parameter
+   !
+   rn_Dt       = 5760.     !  ocean time step (dynamics and T-S)
+   !
    ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time
+   !
    rn_isfhmin  =    1.00   !  treshold [m] to discriminate grounding ice from floating ice
-   !
-   rn_Dt       = 5760.     !  time step for the dynamics and tracer
-   rn_atfp     =    0.1    !  asselin time filter parameter
    !
    ln_crs      = .false.   !  Logical switch for coarsening module      (T => fill namcrs)
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceistate.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceistate.F90	(revision 10001)
@@ -94,9 +94,9 @@
       REAL(wp) ::   zarg, zV, zconv, zdv, zfac
       INTEGER , DIMENSION(4)           ::   itest
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
       REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
       REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini , za_i_ini                        !data by cattegories to fill
+      REAL(wp), DIMENSION(jpi,jpj)     ::   z_ssh_h0, zsshu, zsshv, zsshf
       !--------------------------------------------------------------------
 
@@ -413,47 +413,19 @@
       snwice_mass_b(:,:) = snwice_mass(:,:)
       !
-      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
+      IF( ln_ice_embd ) THEN           ! embedded sea-ice: deplete the initial ssh below sea-ice area
          !
          sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rho0
          sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rho0
          !
-         IF( .NOT.ln_linssh ) THEN
-            !
-            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
-            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
-            !
-            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                
-               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
-               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-               e3t_a(:,:,jk) = e3t_n(:,:,jk)
-            END DO
-            !
-            ! Reconstruction of all vertical scale factors at now and before time-steps
-            ! =========================================================================
-            ! Horizontal scale factor interpolations
-            ! --------------------------------------
-            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
-            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
-            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
-            ! Vertical scale factor interpolations
-            ! ------------------------------------
-            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
-            ! t- and w- points depth
-            ! ----------------------
-            !!gm not sure of that....
-            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-            gdepw_n(:,:,1) = 0.0_wp
-            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-            DO jk = 2, jpk
-               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
-               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
-               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
-            END DO
+         IF( .NOT.ln_linssh ) THEN     ! modified the now and before vertical mesh and scale factors 
+            !
+            !                             !* BEFORE fields : 
+            CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
+            !                                    !  e3t, e3w, e3u, e3uw, e3v, e3vw
+            !
+            !                             !* NOW fields : 
+            CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
+            !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f
+            !                                !      gdept_n, gdepw_n, gde3w_n
          ENDIF
       ENDIF
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asminc.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asminc.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asminc.F90	(revision 10001)
@@ -58,15 +58,18 @@
     LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
 #endif
-   LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields
-   LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment
-   LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization
-   LOGICAL, PUBLIC :: ln_trainc     !: No tracer (T and S) assimilation increments
-   LOGICAL, PUBLIC :: ln_dyninc     !: No dynamics (u and v) assimilation increments
-   LOGICAL, PUBLIC :: ln_sshinc     !: No sea surface height assimilation increment
-   LOGICAL, PUBLIC :: ln_seaiceinc  !: No sea ice concentration increment
-   LOGICAL, PUBLIC :: ln_salfix     !: Apply minimum salinity check
-   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
-   INTEGER, PUBLIC :: nn_divdmp     !: Apply divergence damping filter nn_divdmp times
-
+   LOGICAL, PUBLIC ::   ln_bkgwri     !: No output of the background state fields
+   LOGICAL, PUBLIC ::   ln_asmiau     !: No applying forcing with an assimilation increment
+   LOGICAL, PUBLIC ::   ln_asmdin     !: No direct initialization
+   LOGICAL, PUBLIC ::   ln_trainc     !: No tracer (T and S) assimilation increments
+   LOGICAL, PUBLIC ::   ln_dyninc     !: No dynamics (u and v) assimilation increments
+   LOGICAL, PUBLIC ::   ln_sshinc     !: No sea surface height assimilation increment
+   LOGICAL, PUBLIC ::   ln_seaiceinc  !: No sea ice concentration increment
+   LOGICAL, PUBLIC ::   ln_salfix     !: Apply minimum salinity check
+   LOGICAL, PUBLIC ::   ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
+   INTEGER, PUBLIC ::   nn_divdmp     !: Apply divergence damping filter nn_divdmp times
+
+   LOGICAL, PUBLIC ::   l_dynasm      !: flag derived from ln_asm... flags for calling 
+   LOGICAL, PUBLIC ::   l_traasm      !: dyn_asm_inc and tra_asm_inc routines in step.F90
+   
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diawri.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diawri.F90	(revision 10001)
@@ -118,7 +118,7 @@
       ! 
       ! Output the initial state and forcings
-      IF( ninist == 1 ) THEN                       
+      IF( ln_istate ) THEN                       
          CALL dia_wri_state( 'output.init', kt )
-         ninist = 0
+         ln_istate = .FALSE.
       ENDIF
 
@@ -444,7 +444,7 @@
       IF( ln_timing )   CALL timing_start('dia_wri')
       !
-      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
+      IF( ln_istate ) THEN     !==  Output the initial state and forcings  ==!
          CALL dia_wri_state( 'output.init', kt )
-         ninist = 0
+         ln_istate = .FALSE.
       ENDIF
       !
@@ -878,5 +878,5 @@
       !!
       !! ** Method  :   NetCDF files using ioipsl
-      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
+      !!      File 'output.init.nc'  is created if ln_istate = T (namelist)
       !!      File 'output.abort.nc' is created in case of abnormal job end
       !!----------------------------------------------------------------------
@@ -1029,5 +1029,5 @@
       CALL histclo( id_i )
 #if ! defined key_iomput
-      IF( ninist /= 1  ) THEN
+      IF( ln_istate  ) THEN
          CALL histclo( nid_T )
          CALL histclo( nid_U )
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90	(revision 10001)
@@ -30,9 +30,11 @@
    !! ----------------------------
    !                                   !!* Namelist namdom : time & space domain *
+   LOGICAL , PUBLIC ::   ln_RK3         !: 2rd order Runge-Kutta scheme 
+   LOGICAL , PUBLIC ::   ln_MLF         !: Modified Leap-Frog    scheme (Leclair & Madec OM2010)
+   REAL(wp), PUBLIC ::      rn_atfp        !: asselin time filter parameter
+   REAL(wp), PUBLIC ::   rn_Dt          !: time step for the dynamics and tracer (for both RK3 & MLF)
    LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time
    LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc)
    REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice
-   REAL(wp), PUBLIC ::   rn_dt          !: time step for the dynamics and tracer
-   REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter
    LOGICAL , PUBLIC ::   ln_1st_euler   !: =0 start with forward time step or not (=1)
    LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet
@@ -143,12 +145,15 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hf_0                                  !: v-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ht_0                                 !: inverse of u-depth [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hf_0                                 !: inverse of v-depth [1/m]
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
    INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1) 
 
-   !! 1D reference  vertical coordinate
-   !! =-----------------====------
+   !! 1D reference vertical coordinate
+   !! ==------------------------------
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_1d  , e3w_1d   !: reference vertical scale factors at T- and W-pts (m)
@@ -170,5 +175,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                    (ISF)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts
@@ -265,9 +270,9 @@
          &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
          !
-      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           &
-         &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     &
-         &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     &
-         &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  )
-         !
+      ALLOCATE( ht_0(jpi,jpj)    , hu_0(jpi,jpj)   , hv_0(jpi,jpj)    , hf_0(jpi,jpj)    ,                     &
+         &                         hu_b(jpi,jpj)   , hv_b(jpi,jpj)    , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,  &
+         &      ht_n(jpi,jpj)    , hu_n(jpi,jpj)   , hv_n(jpi,jpj)    , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,  &
+         &                         hu_a(jpi,jpj)   , hv_a(jpi,jpj)    , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) ,  &
+         &      r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj), r1_hv_0(jpi,jpj) , r1_hf_0(jpi,jpj) , STAT=ierr(6)  )
          !
       ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) )
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90	(revision 10001)
@@ -78,5 +78,4 @@
       CHARACTER (len=*), INTENT(IN) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables
       INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
-      REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0
       !!----------------------------------------------------------------------
       !
@@ -138,48 +137,51 @@
       CALL dom_msk( ik_top, ik_bot )   ! Masks
       IF( ln_closea )   CALL dom_clo   ! ln_closea=T : closed seas included in the simulation
-                                       ! Read in masks to define closed seas and lakes 
+      !                                ! Read in masks to define closed seas and lakes 
       !
       DO jj = 1, jpj                   ! depth of the iceshelves
          DO ji = 1, jpi
             ik = mikt(ji,jj)
-            risfdep(ji,jj) = gdepw_0(ji,jj,ik)
+            risfdep(ji,jj) = gdepw_0(ji,jj,ik)        !!gm  RENAME it as h_isf(:,:)  better no?
          END DO
       END DO
       !
-      ht_0(:,:) = 0._wp  ! Reference ocean thickness
-      hu_0(:,:) = 0._wp
+      ht_0(:,:) = 0._wp  ! Reference water column thickness  (wet point only, i.e. without ISF thickness
+      hu_0(:,:) = 0._wp  !           ------------                                  and zero over land    )    
       hv_0(:,:) = 0._wp
+      hf_0(:,:) = 0._wp
       DO jk = 1, jpk
          ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
          hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
          hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
+         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * fmask(:,:,jk)
       END DO
       !
+      r1_ht_0(:,:) =  ssmask(:,:) / ( ht_0(:,:) + 1._wp -  ssmask(:,:) )    ! =0 on lands
+      r1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )    ! not in ISF areas
+      r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp - ssfmask(:,:) )
+      !
       !           !==  time varying part of coordinate system  ==!
       !
       IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all
-      !
-         !       before        !          now          !       after         !
-            gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
-            gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
-                                   gde3w_n = gde3w_0   !        ---          !
-         !                                                                  
-              e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
-              e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
-              e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
-                                     e3f_n =   e3f_0   !        ---          !
-              e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
-             e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
-             e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
-         !
-         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
-         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
-         !
-         !        before       !          now          !       after         !
-                                      ht_n =    ht_0   !                     ! water column thickness
-               hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   ! 
-               hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
-            r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
-            r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
+         !
+         !    before         !          now          !       after         !
+         gdept_b = gdept_0   ;   gdept_n = gdept_0   !        ---          !    depth 
+         gdepw_b = gdepw_0   ;   gdepw_n = gdepw_0   !        ---          !     of
+                                 gde3w_n = gde3w_0   !        ---          ! grid-points
+         !                   !                       !                     !                                              
+          e3t_b  =   e3t_0   ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    !  scale 
+          e3u_b  =   e3u_0   ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! factors
+          e3v_b  =   e3v_0   ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
+                                   e3f_n =   e3f_0   !        ---          !
+          e3w_b  =   e3w_0   ;     e3w_n =   e3w_0   !        ---          !
+         e3uw_b  =  e3uw_0   ;    e3uw_n =  e3uw_0   !        ---          !
+         e3vw_b  =  e3vw_0   ;    e3vw_n =  e3vw_0   !        ---          !
+         !                   !
+                                    ht_n =    ht_0   !                     ! water column
+            hu_b =    hu_0   ;      hu_n =    hu_0   ;    hu_a =    hu_0   !  thickness
+            hv_b =    hv_0   ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
+         r1_hu_b = r1_hu_0   ;   r1_hu_n = r1_hu_0   ; r1_hu_a = r1_hu_0   ! 1 / water column
+         r1_hv_b = r1_hv_0   ;   r1_hv_n = r1_hv_0   ; r1_hv_a = r1_hv_0   !       thickness
          !
          !
@@ -192,9 +194,9 @@
       IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
       !
-      IF( ln_meshmask .AND. .NOT.ln_iscpl )                        CALL dom_wri     ! Create a domain file
+      IF( ln_meshmask .AND. .NOT.ln_iscpl                      )   CALL dom_wri     ! Create a domain file
       IF( ln_meshmask .AND.      ln_iscpl .AND. .NOT.ln_rstart )   CALL dom_wri     ! Create a domain file
       IF(                                       .NOT.ln_rstart )   CALL dom_ctl     ! Domain control
       !
-      IF( ln_write_cfg )   CALL cfg_write         ! create the configuration file
+      IF( ln_write_cfg )   CALL cfg_write       ! create the configuration file
       !
       IF(lwp) THEN
@@ -204,5 +206,5 @@
          WRITE(numout,*) 
       ENDIF
-      !
+       !
    END SUBROUTINE dom_init
 
@@ -290,8 +292,8 @@
       NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                   &
          &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl   ,     &
-         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate   ,     &
+         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , ln_istate   ,     &
          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, ln_1st_euler,     &
          &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios
-      NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_Dt, rn_atfp, ln_crs, ln_meshmask
+      NAMELIST/namdom/ ln_RK3, ln_MLF, rn_Dt, rn_atfp, ln_linssh, rn_isfhmin, ln_crs, ln_meshmask
 #if defined key_netcdf4
       NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
@@ -305,4 +307,33 @@
       ENDIF
       !
+      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
+      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
+      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
+904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namdom )
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain'
+         WRITE(numout,*) '      3rd order Runge-Kutta scheme            ln_RK3      = ', ln_RK3
+         WRITE(numout,*) '      Modified Leap-Frog    scheme            ln_MLF      = ', ln_MLF
+         WRITE(numout,*) '         asselin time filter parameter           rn_atfp  = ', rn_atfp
+         WRITE(numout,*) '      ocean time step (RK3 & MLF)             rn_Dt       = ', rn_Dt
+         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
+         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
+         WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]'
+         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
+         WRITE(numout,*)  
+      ENDIF
+      !
+      IF ( (      ln_RK3 .AND. .NOT.ln_MLF ) .OR.   &
+         & ( .NOT.ln_RK3 .AND.      ln_MLF ) ) THEN
+         IF( ln_RK3 .AND. lwp )   WRITE(numout,*)'   ==>>>   a RK3 time-stepping scheme is used'
+         IF( ln_MLF .AND. lwp )   WRITE(numout,*)'   ==>>>   a MLF time-stepping scheme is used'
+      ELSE
+         CALL ctl_stop( 'dom_nam: select one and only one time stepping scheme')
+      ENDIF
       !
       REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
@@ -315,4 +346,5 @@
       !
       IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
          WRITE(numout,*) '   Namelist : namrun   ---   run parameters'
          WRITE(numout,*) '      Assimilation cycle              nn_no           = ', nn_no
@@ -330,5 +362,5 @@
          WRITE(numout,*) '      initial time of day in hhmm     nn_time0        = ', nn_time0
          WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy        = ', nn_leapy
-         WRITE(numout,*) '      initial state output            nn_istate       = ', nn_istate
+         WRITE(numout,*) '      output the initial state        ln_istate       = ', ln_istate
          IF( ln_rst_list ) THEN
             WRITE(numout,*) '      list of restart dump times      nn_stocklist    =', nn_stocklist
@@ -357,15 +389,17 @@
       ndate0 = nn_date0
       nleapy = nn_leapy
-      ninist = nn_istate
       nstock = nn_stock
       nstocklist = nn_stocklist
       nwrite = nn_write
-      IF( ln_rstart ) THEN       ! restart : set 1st time-step scheme (LF or forward) 
-         l_1st_euler = ln_1st_euler
-      ELSE                       ! start from rest : always an Euler scheme for the 1st time-step
-         IF(lwp) WRITE(numout,*)  
-         IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)'
-         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used '   
-         l_1st_euler = .TRUE.
+      
+      IF( ln_MLF ) THEN       ! Leap-Frog only
+         IF( ln_rstart ) THEN       ! restart : set 1st time-step scheme (LF or forward) 
+            l_1st_euler = ln_1st_euler
+         ELSE                       ! start from rest : always an Euler scheme for the 1st time-step
+            IF(lwp) WRITE(numout,*)  
+            IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)'
+            IF(lwp) WRITE(numout,*)'           an Euler initial time step is used '   
+            l_1st_euler = .TRUE.
+         ENDIF
       ENDIF
       !                             ! control of output frequency
@@ -399,27 +433,7 @@
       ENDIF
 #endif
-
-      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
-      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
-903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
-      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
-      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
-904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namdom )
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain'
-         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
-         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
-         WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]'
-         WRITE(numout,*) '      ocean time step                         rn_dt       = ', rn_dt
-         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
-         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
-      ENDIF
-      !
-      IF( TRIM(Agrif_CFixed()) == '0' ) THEN
+      !
+      IF( TRIM(Agrif_CFixed()) == '0' ) THEN    !set output file type for XIOS based on NEMO namelist 
          lrxios = ln_xios_read.AND.ln_rstart
-!set output file type for XIOS based on NEMO namelist 
          IF (nn_wxios > 0) lwxios = .TRUE. 
          nxioso = nn_wxios
@@ -553,6 +567,6 @@
             CALL iom_getatt( inum, 'cn_cfg', cd_cfg )  ! returns   !  if not found
             CALL iom_getatt( inum, 'nn_cfg', kk_cfg )  ! returns -999 if not found
-            IF( TRIM(cd_cfg) == '!') cd_cfg = 'UNKNOWN'
-            IF( kk_cfg == -999     ) kk_cfg = -9999999
+            IF( TRIM(cd_cfg) == '!')   cd_cfg = 'UNKNOWN'
+            IF( kk_cfg == -999     )   kk_cfg = -9999999
          ENDIF
          !
@@ -701,5 +715,5 @@
       !
    END SUBROUTINE cfg_write
-
+   
    !!======================================================================
 END MODULE domain
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dommsk.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dommsk.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dommsk.F90	(revision 10001)
@@ -202,4 +202,5 @@
       ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
       ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
+      ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 )
 
 
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90	(revision 10001)
@@ -4,8 +4,10 @@
    !! 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
+   !! History :  5.0  !  2018-07  (G. Madec)  Flux Form with Kinetic Energy conservation
+   !!                             ==>>>  here z* and s* only (no z-tilde) 
+   
+   ! 1- remove   z-tilde          ==>>>  pure z-star (or s-star)
+   ! 2- remove   dom_vvl_interpol  
+   
    !!----------------------------------------------------------------------
 
@@ -14,5 +16,4 @@
    !!   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
@@ -38,5 +39,6 @@
    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
+   PUBLIC  ssh2e3_before      ! ...
+   PUBLIC  ssh2e3_now         ! ...
 
    !                                                      !!* Namelist nam_vvl
@@ -61,4 +63,10 @@
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv       ! retoring period for low freq. divergence
 
+
+!!gm add
+!!gm 
+
+
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -74,4 +82,5 @@
       !!                ***  FUNCTION dom_vvl_alloc  ***
       !!----------------------------------------------------------------------
+      !
       IF( ln_vvl_zstar )   dom_vvl_alloc = 0
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
@@ -115,7 +124,4 @@
       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk
-      INTEGER ::   ii0, ii1, ij0, ij1
-      REAL(wp)::   zcoef, z1_Dt
       !!----------------------------------------------------------------------
       !
@@ -129,132 +135,49 @@
       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
       !
-      !                    ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf
+      !                    ! Read or initialize e3t_(b/n), ssh(b/n)
       CALL dom_vvl_rst( nit000, 'READ' )
-      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
-      !
-      !                    !== Set of all other vertical scale factors  ==!  (now and before)
-      !                                ! Horizontal interpolation of e3t
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
-      !                                ! Vertical interpolation of e3t,u,v 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
-
-      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
-      e3t_a(:,:,:) = e3t_n(:,:,:)
-      e3u_a(:,:,:) = e3u_n(:,:,:)
-      e3v_a(:,:,:) = e3v_n(:,:,:)
-      !
-      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
-      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
-      gdepw_b(:,:,1) = 0.0_wp
-      DO jk = 2, jpk                               ! vertical sum
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               !    zcoef = tmask - wmask    ! 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     
-!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
-               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
-               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
-               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
-            END DO
-         END DO
-      END DO
-      !
-      !                    !==  thickness of the water column  !!   (ocean portion only)
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
-      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
-      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
-      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
-      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
-      DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-      END DO
-      !
-      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
-      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
-      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
-      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
-      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
-
-      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
-      IF( ln_vvl_ztilde ) THEN
-!!gm : idea: add here a READ in a file of custumized restoring frequency
-         !                                   ! Values in days provided via the namelist
-         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
-         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
-         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
-         !
-         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
-            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
-            frq_rst_hdv(:,:) = 1._wp / rn_Dt
-         ENDIF
-         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
-            z1_Dt = 1._wp / rn_Dt
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
-                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
-                     ! values outside the equatorial band and transition zone (ztilde)
-                     frq_rst_e3t(ji,jj) =  2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400._wp )
-                     frq_rst_hdv(ji,jj) =  2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400._wp )
-                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
-                     ! values inside the equatorial band (ztilde as zstar)
-                     frq_rst_e3t(ji,jj) =  0._wp
-                     frq_rst_hdv(ji,jj) =  z1_Dt
-                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
-                     !                                      ! (linearly transition from z-tilde to z-star)
-                     frq_rst_e3t(ji,jj) = 0._wp + ( frq_rst_e3t(ji,jj) - 0._wp  ) * 0.5_wp                             &
-                        &                       * (  1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp )  )
-                     frq_rst_hdv(ji,jj) = z1_Dt + (  frq_rst_hdv(ji,jj) - z1_Dt ) * 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( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
-               ii0 = 103   ;   ii1 = 111       
-               ij0 = 128   ;   ij1 = 135   ;   
-               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0._wp
-               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  z1_Dt
-            ENDIF
-         ENDIF
-      ENDIF
-      !
-      IF(lwxios) THEN
-! define variables in restart file when writing with XIOS
+      !
+      !                    !== Set of all other vertical mesh fields  ==!  (now and before)      
+      !
+      !                          !* BEFORE fields : 
+      CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
+      !                                    !  e3t, e3w, e3u, e3uw, e3v, e3vw
+      !
+      !                                ! set one for all last level to the e3._0 value
+      e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
+      e3w_b(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk)
+      !
+      !                          !* NOW fields : 
+      CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
+      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f
+      !                                !      gdept_n, gdepw_n, gde3w_n
+      !
+      !                                ! set one for all last level to the e3._0 value
+      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
+      e3w_n(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk)
+      e3f_n(:,:,jpk) = e3f_0(:,:,jpk)
+      !
+      !                          !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation)
+      e3t_a(:,:,:) = e3t_n(:,:,:)   ;   e3u_a(:,:,:) = e3u_n(:,:,:)   ;   e3v_a(:,:,:) = e3v_n(:,:,:)
+      
+!!gm           
+!!gm        ===>>>>   below:  some issues to think about !!!
+!!gm  
+!!gm   fmask definition checked (0 or 1 nothing else)
+!!        in z-tilde or ALE    e3._0  should be the time varying fields step forward with an euler scheme
+!!gm   e3w_b  & gdept_b are not used   except its update in agrif  
+!!     and used to compute before slope of surface in dynldf_iso ==>>>  remove it !!!!
+!!         NB: in triads on TRA, gdept_n  is used !!!!   BUG?
+!!gm   e3f_n  almost not used  ===>>>>  verify whether it can be removed or not...
+!!       verify the use of wumask & wvmask   mau be replaced by a multiplication by umask(k)*umask(k+1)  ???
+!!
+!!gm  ISF case to be checked by Pierre Mathiot
+!!
+!!gm  setting of e3._a  for agrif....  TO BE CHECKED
+
+      !
+      IF(lwxios) THEN      ! define variables in restart file when writing with XIOS
          CALL iom_set_rstw_var_active('e3t_b')
          CALL iom_set_rstw_var_active('e3t_n')
-         !                                           ! ----------------------- !
-         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
-            !                                        ! ----------------------- !
-            CALL iom_set_rstw_var_active('tilde_e3t_b')
-            CALL iom_set_rstw_var_active('tilde_e3t_n')
-         END IF
-         !                                           ! -------------!    
-         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
-            !                                        ! ------------ !
-            CALL iom_set_rstw_var_active('hdiv_lf')
-         ENDIF
-         !
       ENDIF
       !
@@ -287,10 +210,6 @@
       INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
       !
-      INTEGER                ::   ji, jj, jk            ! dummy loop indices
-      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
-      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
-      LOGICAL                ::   ll_do_bclinic         ! local logical
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
       !!----------------------------------------------------------------------
       !
@@ -305,247 +224,30 @@
       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 e3t_a = e3t_0*(1+ssha/ht_0)
-         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_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._wp
-         zht(:,:)   = 0._wp
-         DO jk = 1, jpkm1
-            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + e3t_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 > nit000 ) THEN
-               DO jk = 1, jpkm1
-                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
-               END DO
-            ENDIF
-         ENDIF
-
-         ! II - after z_tilde increments of vertical scale factors
-         ! =======================================================
-         te3t_a(:,:,:) = 0._wp  ! te3t_a used to store tendency terms
-
-         ! 1 - High frequency divergence term
-         ! ----------------------------------
-         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
-            DO jk = 1, jpkm1
-               te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
-            END DO
-         ELSE                         ! layer case
-            DO jk = 1, jpkm1
-               te3t_a(:,:,jk) = te3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
-            END DO
-         ENDIF
-
-         ! 2 - Restoring term (z-tilde case only)
-         ! ------------------
-         IF( ln_vvl_ztilde ) THEN
-            DO jk = 1, jpk
-               te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk)
-            END DO
-         ENDIF
-
-         ! 3 - Thickness diffusion term
-         ! ----------------------------
-         zwu(:,:) = 0._wp
-         zwv(:,:) = 0._wp
-         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
-                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj  ,jk) )
-                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
-                     &            * ( te3t_b(ji,jj,jk) - te3t_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
-         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
-            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
-         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  te3t_a(ji,jj,jk) = te3t_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_e1e2t(ji,jj)
-               END DO
-            END DO
-         END DO
-         !                       ! d - thickness diffusion transport: boundary conditions
-         !                             (stored for tracer advction and continuity equation)
-         CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
-
-         ! 4 - Time stepping of baroclinic scale factors
-         ! ---------------------------------------------
-         ! Leapfrog time stepping
-         ! ~~~~~~~~~~~~~~~~~~~~~~
-         CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp )
-         te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:)
-
-         ! Maximum deformation control
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ze3t(:,:,jpk) = 0._wp
-         DO jk = 1, jpkm1
-            ze3t(:,:,jk) = te3t_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 >  rn_zdef_max ) .OR. ( z_tmin < - 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( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
-               WRITE(numout, *) 'at i, j, k=', ijk_max
-               WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
-               WRITE(numout, *) 'at i, j, k=', ijk_min            
-               CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
-            ENDIF
-         ENDIF
-         ! - ML - end test
-         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
-         te3t_a(:,:,:) = MIN( te3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
-         te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
-
-         !
-         ! "tilda" change in the after scale factor
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         DO jk = 1, jpkm1
-            dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_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._wp
-         DO jk = 1, jpkm1
-            zht(:,:)  = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
-         DO jk = 1, jpkm1
-            dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_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
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_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(te3t_a))) =', z_tmax
-         END IF
-         !
-         zht(:,:) = 0.0_wp
-         DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_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(e3t_n))) =', z_tmax
-         !
-         zht(:,:) = 0.0_wp
-         DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_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(e3t_a))) =', z_tmax
-         !
-         zht(:,:) = 0.0_wp
-         DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_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(e3t_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( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
-
-      ! *********************************** !
-      ! After depths at u- v points         !
-      ! *********************************** !
-
-      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
-      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
-      DO jk = 2, jpkm1
-         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
-         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
-      END DO
-      !                                        ! Inverse of the local depth
-!!gm BUG ?  don't understand the use of umask_i here .....
+      !                             ! ------------------!
+      !                             ! z_star coordinate !     (and barotropic z-tilde part)
+      !                             ! ------------------!
+      !
+      !                                   !==  after ssh  ==!  (u- and v-points)
+      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
+         zsshu_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji+1,jj) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji,jj+1) ) * ssvmask(ji,jj)
+      END DO             ;   END DO      
+      CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
+      !
+      !                                   !==  after depths and its inverse  ==! 
+         hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:)
+         hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:)
       r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
       r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
+      !
+      !                                   !==  after scale factors  ==!  (e3t , e3u , e3v)
+      zssht_h(:,:) = ssha   (:,:) * r1_ht_0(:,:)           ! t-point
+      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! u-point
+      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! v-point
+      DO jk = 1, jpkm1
+         e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+         e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
+         e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
+      END DO
       !
       IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
@@ -581,5 +283,5 @@
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      REAL(wp) ::   zcoef, ze3f   ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
       !!----------------------------------------------------------------------
       !
@@ -594,75 +296,87 @@
       ENDIF
       !
-      ! Time filter and swap of scale factors
-      ! =====================================
-      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         IF( l_1st_euler ) THEN
-            te3t_n(:,:,:) = te3t_a(:,:,:)
-         ELSE
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ze3f = te3t_n(ji,jj,jk)   &
-                        & + rn_atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) )
-                     te3t_b(ji,jj,jk) = ze3f
-                     te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ENDIF
-      ENDIF
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-
-      e3t_n(:,:,:) = e3t_a(:,:,:)
-      e3u_n(:,:,:) = e3u_a(:,:,:)
-      e3v_n(:,:,:) = e3v_a(:,:,:)
-
-      ! Compute all missing vertical scale factor and depths
-      ! ====================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
+      ! Swap and Compute all missing vertical scale factor and depths
+      ! =============================================================
       ! - ML - e3u_b and e3v_b are allready computed in dynnxt
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
-      
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
-      
-      ! Vertical scale factor interpolations
-      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
-
-      ! t- and w- points depth (set the isf depth as it is in the initial step)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_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))
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
-                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
-            END DO
-         END DO
-      END DO
-
+      !
+      ! - GM - to be updated :   e3f_n,  e3w_n , e3uw_n , e3vw_n 
+      !                                  e3w_b , e3uw_b , e3vw_b
+      !                          gdept_n , gdepw_n , gde3w_n
+      !                          ht_n
+      !        to be swap    :   hu_a , hv_a  , r1_hu_a , r1_hv_a
+      !
       ! Local depth and Inverse of the local depth of the water
       ! -------------------------------------------------------
+      !                       ! swap of depth and scale factors
+      !                       ! ===============================
+      DO jk = 1, jpkm1                          ! swap n--> a  
+         gdept_b(:,:,jk) = gdept_n(:,:,jk)         ! depth at t and w
+         gdepw_b(:,:,jk) = gdepw_n(:,:,jk)
+         e3t_n  (:,:,jk) = e3t_a  (:,:,jk)         ! e3t, e3u, e3v
+         e3u_n  (:,:,jk) = e3u_a  (:,:,jk)
+         e3v_n  (:,:,jk) = e3v_a  (:,:,jk)
+      END DO
+      ht_n(:,:) = ht_0(:,:) + sshn(:,:)            ! ocean thickness
+      !
       hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
       hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
       !
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
-      DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+      !                    !==  before : 
+      !                                            !* ssh at u- and v-points)
+      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
+         zsshu_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+      END DO             ;   END DO      
+      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
+      !
+      !                                            !*  e3w_b , e3uw_b , e3vw_b
+      zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:)           ! w-point
+      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! uw-point
+      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! vw-point
+      DO jk = 1, jpkm1
+          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
+         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
       END DO
-
+      !      
+      !                    !==  now    : 
+      !                                            !* ssh at u- and v-points)
+      DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 for f-point
+         zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj)
+         zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
+            &                       + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj)
+      END DO             ;   END DO      
+      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )      
+      !
+      !                                            !* e3w_n , e3uw_n , e3vw_n, e3f_n 
+      zssht_h(:,:) = sshn   (:,:) * r1_ht_0(:,:)           ! t- & w-point
+      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! uw-point
+      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! vw-point
+      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)           ! f-point
+      DO jk = 1, jpkm1
+          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,jk) )
+         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
+         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
+          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) *  fmask(:,:,jk) )
+      END DO      
+      ! 
+      zssht_h(:,:) = 1._wp + sshn (:,:) * r1_ht_0(:,:)   ! t-point
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpkm1
+            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
+            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
+            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - sshn   (:,:)
+         END DO
+      ELSE                    ! no ISF cavities 
+         DO jk = 1, jpkm1
+            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+         END DO
+      ENDIF
+      !
       ! write restart file
       ! ==================
@@ -672,111 +386,4 @@
       !
    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
-      !!----------------------------------------------------------------------
-      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'
-      !
-      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
-      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
-      !!----------------------------------------------------------------------
-      !
-      IF(ln_wd_il) THEN
-        zlnwd = 1.0_wp
-      ELSE
-        zlnwd = 0.0_wp
-      END IF
-      !
-      SELECT CASE ( pout )    !==  type of interpolation  ==!
-         !
-      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-         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) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
-                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
-         !
-      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-         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)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
-                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
-         !
-      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         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) * (1.0_wp - zlnwd) + zlnwd ) &
-                     &                       *    r1_e1e2f(ji,jj)                                                  &
-                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
-         !
-      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
-         !
-         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
-         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
-!!gm BUG? use here wmask in case of ISF ?  to be checked
-         DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
-               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
-               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
-               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
-         END DO
-         !
-      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
-         !
-         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
-!!gm BUG? use here wumask in case of ISF ?  to be checked
-         DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
-               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
-               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
-               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
-         END DO
-         !
-      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
-         !
-         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
-!!gm BUG? use here wvmask in case of ISF ?  to be checked
-         DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
-               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
-               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
-               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
-         END DO
-      END SELECT
-      !
-   END SUBROUTINE dom_vvl_interpol
 
 
@@ -804,119 +411,70 @@
          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, ldxios = lrxios    )
             !
-            id1 = iom_varid( numror, 'e3t_b'      , ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'e3t_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 !
-            !                             ! --------- !
+            id1 = iom_varid( numror, 'sshb'      , ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'sshn'      , ldstop = .FALSE. )
+            !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
-               ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
-               WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  e3t_n(:,:,:) = e3t_0(:,:,:)
-                  e3t_b(:,:,:) = e3t_0(:,:,:)
-               END WHERE
+               IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files'
+               !
+!!gm  Question: use jpdom_data above to read data over jpi x jpj    (like is dom_hgr_read and dom_zgr_read)
+!!              so that it will work with land processor suppression
+!               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
+!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios    )
+!!gm 
+               CALL iom_get( numror, jpdom_data, 'sshn'   , sshn, ldxios = lrxios    )
+               CALL iom_get( numror, jpdom_data, 'sshb'   , sshb, ldxios = lrxios    )
+!!gm end
                IF( l_1st_euler ) THEN
-                  e3t_b(:,:,:) = e3t_n(:,:,:)
+                  sshb(:,:) = sshn(:,:)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
-               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               e3t_n(:,:,:) = e3t_b(:,:,:)
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files'
+               IF(lwp) write(numout,*) '   set sshn = sshb  and force l_1st_euler = true'
+!!gm               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_data, 'sshb', sshb, ldxios = lrxios )
+               sshn(:,:) = sshb(:,:)
                l_1st_euler = .TRUE.
             ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
-               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files'
+               IF(lwp) write(numout,*) 'set sshb = sshn  and force l_1st_euler = true'
+               CALL iom_get( numror, jpdom_data, 'sshn', sshb, ldxios = lrxios )
+               sshb(:,:) = sshn(:,:)
                l_1st_euler = .TRUE.
             ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
-               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
-               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
-               DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
-                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
-               END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file'
+               IF(lwp) write(numout,*) 'set sshb = sshn = 0  and force l_1st_euler = true'
+               sshb(:,:) = 0._wp
+               sshn(:,:) = 0._wp
                l_1st_euler = .TRUE.
             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', te3t_b(:,:,:), ldxios = lrxios )
-                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios )
-               ELSE                            ! one at least array is missing
-                  te3t_b(:,:,:) = 0.0_wp
-                  te3t_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(:,:,:), ldxios = lrxios )
-                  ELSE                ! array is missing
-                     hdiv_lf(:,:,:) = 0.0_wp
-                  ENDIF
-               ENDIF
-            ENDIF
-            !
          ELSE                                   !* Initialize at "rest"
             !
-
             IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential 
                !
-               IF( cn_cfg == 'wad' ) THEN
-                  ! Wetting and drying test case
+               IF( cn_cfg == 'wad' ) THEN             ! Wetting and drying test case
                   CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
-                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
+                  tsn  (:,:,:,:) = tsb (:,:,:,:)            ! set now values from to before ones
                   sshn (:,:)     = sshb(:,:)
                   un   (:,:,:)   = ub  (:,:,:)
                   vn   (:,:,:)   = vb  (:,:,:)
-               ELSE
-                  ! if not test case
+               ELSE                                   ! Not the test case
                   sshn(:,:) = -ssh_ref
                   sshb(:,:) = -ssh_ref
-
+                  !
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
-
+                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN  ! if total depth is less than min depth
                            sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
                            sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
                            ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
                         ENDIF
-                     ENDDO
-                  ENDDO
-               ENDIF !If test case else
-
-               ! Adjust vertical metrics for all wad
-               DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
-                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
-               END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
-
+                     END DO
+                  END DO
+               ENDIF     ! If test case else
+               !
                DO ji = 1, jpi
                   DO jj = 1, jpj
-                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
+                     IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN
                        CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
                      ENDIF
@@ -926,48 +484,18 @@
             ELSE
                !
-               ! Just to read set ssh in fact, called latter once vertical grid
-               ! is set up:
+               ! Just to read set ssh in fact, called latter once vertical grid is set up:
 !               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
-!               !
-!               DO jk=1,jpk
-!                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
-!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
-!               END DO
-!               e3t_n(:,:,:) = e3t_b(:,:,:)
-                sshn(:,:)=0._wp
-                e3t_n(:,:,:)=e3t_0(:,:,:)
-                e3t_b(:,:,:)=e3t_0(:,:,:)
+               sshn(:,:) = 0._wp
+               sshb(:,:) = 0._wp
                !
-            END IF           ! end of ll_wd edits
-
-            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
-               te3t_b(:,:,:) = 0._wp
-               te3t_n(:,:,:) = 0._wp
-               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
             END IF
+            !
          ENDIF
          !
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
          !                                   ! ===================
-         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
-         IF( lwxios ) CALL iom_swap(      cwxios_context          )
-         !                                           ! --------- !
-         !                                           ! all cases !
-         !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
-         !                                           ! ----------------------- !
-         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
-            !                                        ! ----------------------- !
-            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios)
-            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios)
-         END IF
-         !                                           ! -------------!    
-         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
-            !                                        ! ------------ !
-            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
-         ENDIF
-         !
-         IF( lwxios ) CALL iom_swap(      cxios_context          )
+
+!!gm      DO NOTHING,   sshb and sshn  are written in restart.F90
+
       ENDIF
       !
@@ -1024,4 +552,10 @@
       ENDIF
       !
+
+!!gm 
+      IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar )   CALL ctl_stop( 'z-tilde NOT available in this branch' )
+!!gm
+
+      !
       ioptio = 0                      ! Parameter control
       IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
@@ -1047,4 +581,119 @@
    END SUBROUTINE dom_vvl_ctl
 
+   
+   SUBROUTINE ssh2e3_now
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ssh2e3_now  ***
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
+      !!----------------------------------------------------------------------
+      !
+      !                             !==  ssh at u- and v-points  ==!
+      !
+      DO jj = 1, jpjm1                    ! start from 1 due to f-point
+         DO ji = 1, jpim1
+            zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj)
+            zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
+               &                       + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj)
+         END DO
+      END DO      
+      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )
+      !
+      !                             !==  ht, hu and hv  == !   (and their inverse)
+      !
+      ht_n   (:,:) = ht_0(:,:) +  sshn  (:,:)
+      hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:)
+      hv_n   (:,:) = hv_0(:,:) + zsshv_h(:,:)
+      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF
+      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+      !      
+      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
+      !
+      zssht_h(:,:) =  sshn  (:,:) * r1_ht_0(:,:)
+      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)
+      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)
+      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)
+      !
+      !                             !==  e3t, e3w  ,  e3u, e3uw ,  e3v, e3vw  , and e3f  ==!
+      !      
+      DO jk = 1, jpkm1
+          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          !
+          e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) *  umask(:,:,jk) )
+         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
+         !
+          e3v_n(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) *  vmask(:,:,jk) )
+         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
+          !
+          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
+      END DO
+      !      
+      !                             !== depth of t- and w-points  ==!
+      !
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + sshn / ht_0
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpkm1
+            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
+            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
+            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - sshn(:,:)
+         END DO
+      ELSE                    ! no ISF cavities
+         DO jk = 1, jpkm1
+            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE ssh2e3_now
+   
+   
+   SUBROUTINE ssh2e3_before
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ssh2e3_before  ***
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
+      !!----------------------------------------------------------------------
+      !
+      !                             !==  ssh at u- and v-points  ==!
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            zsshu_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+         END DO
+      END DO      
+      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
+      !
+      !                             !==  ht, hu and hv  == !   (and their inverse)
+         hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:)
+         hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:)
+      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
+      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
+      !
+      !      
+      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
+      zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:)
+      zsshu_h   (:,:) = zsshu_h(:,:) * r1_hu_0(:,:)
+      zsshv_h   (:,:) = zsshv_h(:,:) * r1_hv_0(:,:)
+      !
+      !                             !==  e3t, e3w  ,  e3u, e3uw , and  e3v, e3vw  ==!
+      DO jk = 1, jpkm1
+          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          !
+          e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h  (:,:) *  umask(:,:,jk) )
+         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h  (:,:) * wumask(:,:,jk) )
+          !
+          e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h  (:,:) *  vmask(:,:,jk) )
+         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h  (:,:) * wvmask(:,:,jk) )
+      END DO
+      !   
+   END SUBROUTINE ssh2e3_before
+   
    !!======================================================================
 END MODULE domvvl
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90	(revision 10001)
@@ -14,5 +14,4 @@
    USE dom_oce         ! ocean space and time domain
    USE domwri          ! ocean space and time domain
-   USE domvvl   , ONLY : dom_vvl_interpol
    USE phycst          ! physical constants
    USE sbc_oce         ! surface boundary condition variables
@@ -137,7 +136,8 @@
       REAL(wp), DIMENSION(jpi,jpj)          :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t
       REAL(wp), DIMENSION(jpi,jpj)          :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1
-      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn, ztrp
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn
       REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
+      REAL(wp), DIMENSION(jpi,jpj) ::   z_ssh_h0, zsshu, zsshv, zsshf
       !!----------------------------------------------------------------------
       !
@@ -163,13 +163,11 @@
          DO jj = 2,jpj-1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               jip1=ji+1; jim1=ji-1;
-               jjp1=jj+1; jjm1=jj-1;
-               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
+               summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1)
                IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
-                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
-                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
-                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
-                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
-                  zsmask1(ji,jj)=1._wp
+                  sshn(ji,jj)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj)     &
+                  &           + zssh0(ji-1,jj)*zsmask0(ji-1,jj)     &
+                  &           + zssh0(ji,jj+1)*zsmask0(ji,jj+1)     &
+                  &           + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk
+                  zsmask1(ji,jj) = 1._wp
                ENDIF
             END DO
@@ -185,76 +183,95 @@
       IF ( .NOT.ln_linssh ) THEN
       ! Reconstruction of all vertical scale factors at now time steps
-      ! =============================================================================
+      ! ======================================================================
+      
+!!gm Question : bug ????
+      ! at this stage  is ht_0, r1_ht0 already computed ????
+      ! I have the feeling that this should be done in a different manner...
+      ! that is iscpl_rst_interpol  should be call directly in restart !!!
+      
+      ! in my changes here I assume that ht_0 AND  r1_ht_0  have been updated 
+      ! Note that the former calculation were using ht_0  so if it as not been updated ===>>> BUG
+!!gm
+      
+      
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-         DO jk = 1,jpk
-            DO jj=1,jpj
-               DO ji=1,jpi
-                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
-                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) /   &
-                     &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
-
-         ! Vertical scale factor interpolations
-         ! ------------------------------------
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp ) THEN
+                  DO jk = 1, jpk
+                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
+                  END DO
+               ENDIF
+            END DO
+         END DO
+         !
+!!gm  Note that if this routine is called in dom_vvl_init then all the lines below are uselss !!!
+!!        they are a duplication of dom_vvl_init lines
+
+         !                                   !==  now fields  ==!
+         !
+         !                                            !* ssh at u- and v-points)
+         DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 due to f-point
+            zsshu(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj)
+            zsshv(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj)
+            zsshf(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
+               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj)
+         END DO             ;   END DO      
+         CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp , zsshf(:,:),'F', 1._wp )
+         !
+         !                                            !* hu and hv (and their inverse) 
+         ht_n   (:,:) = ht_0(:,:) +  sshn(:,:)
+         hu_n   (:,:) = hu_0(:,:) + zsshu(:,:)
+         hv_n   (:,:) = hv_0(:,:) + zsshv(:,:)
+         r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF
+         r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+         !
+         !                                            !* e3u, e3uw  and  e3v, e3vw
+         z_ssh_h0(:,:) = sshn (:,:) * r1_ht_0(:,:)           ! t-point
+         DO jk = 1, jpkm1
+            e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * tmask(:,:,jk) )
+         END DO
+         z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:)           ! u-point
+         DO jk = 1, jpkm1
+            e3u_n (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) *  umask(:,:,jk) )
+            e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wumask(:,:,jk) )
+         END DO
+         z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:)           ! v-point
+         DO jk = 1, jpkm1
+            e3v_n (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) *  vmask(:,:,jk) )
+            e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wvmask(:,:,jk) )
+         END DO
+         z_ssh_h0(:,:) = zsshf(:,:) * r1_hf_0(:,:)           ! f-point
+         DO jk = 1, jpkm1
+            e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * fmask(:,:,jk) )
+         END DO
+
+         z_ssh_h0(:,:) = 1._wp + sshn(:,:) * r1_ht_0(:,:)    ! t-point
+         !
+         IF( ln_isfcav ) THEN    ! iceshelf cavities : ssh scaling not applied over the iceshelf thickness 
+            DO jk = 1, jpkm1
+               gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:)
+               gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:)
+               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+            END DO
+         ELSE
+            DO jk = 1, jpkm1
+               gdept_n(:,:,jk) = gdept_0(:,:,jk) * z_ssh_h0(:,:)
+               gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * z_ssh_h0(:,:)
+               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+            END DO
+         ENDIF
          
-         ! t- and w- points depth
-         ! ----------------------
-         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-         gdepw_n(:,:,1) = 0.0_wp
-         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               DO jk = 2,mikt(ji,jj)-1
-                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
-               END DO
-               IF (mikt(ji,jj) > 1) THEN
-                  jk = mikt(ji,jj)
-                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
-               END IF
-               DO jk = mikt(ji,jj)+1, jpk
-                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
-               END DO
-            END DO
-         END DO
-
-      ! t-, u- and v- water column thickness
-      ! ------------------------------------
-         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         !                                        ! Inverse of the local depth
-         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
-         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
-
-      END IF
+      ENDIF
 
 !=============================================================================
 ! compute velocity
 ! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
-      ub(:,:,:)=un(:,:,:)
-      vb(:,:,:)=vn(:,:,:)
-      DO jk = 1,jpk
-         DO jj = 1,jpj
-            DO ji = 1,jpi
+      ub(:,:,:) = un(:,:,:)
+      vb(:,:,:) = vn(:,:,:)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
                un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk);
                vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk);
@@ -265,12 +282,8 @@
 ! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
 ! compute barotropic velocity now and after 
-      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
-      zbub(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
-      zbvb(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
-      zbun(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
-      zbvn(:,:)   = SUM(ztrp,DIM=3)
+      zbub(:,:)   = SUM( ub(:,:,:)*pe3u_b(:,:,:) , DIM=3 )
+      zbvb(:,:)   = SUM( vb(:,:,:)*pe3v_b(:,:,:) , DIM=3 )
+      zbun(:,:)   = SUM( un(:,:,:)* e3u_n(:,:,:) , DIM=3 )
+      zbvn(:,:)   = SUM( vn(:,:,:)* e3v_n(:,:,:) , DIM=3 )
 
       ! new water column
@@ -285,6 +298,6 @@
       zucorr = 0._wp
       zvcorr = 0._wp
-      DO jj = 1,jpj
-         DO ji = 1,jpi
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
                zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
@@ -297,7 +310,7 @@
       
       ! update velocity
-      DO jk  = 1,jpk
-         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
-         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
+      DO jk  = 1, jpk
+         un(:,:,jk) = ( un(:,:,jk) - zucorr(:,:) ) * umask(:,:,jk)
+         vn(:,:,jk) = ( vn(:,:,jk) - zvcorr(:,:) ) * vmask(:,:,jk)
       END DO
 
@@ -305,49 +318,49 @@
       ! compute temp and salt
       ! compute new tn and sn if we open a new cell
-      tsb (:,:,:,:) = tsn(:,:,:,:)
-      zts0(:,:,:,:) = tsn(:,:,:,:)
+      tsb  (:,:,:,:) = tsn   (:,:,:,:)
+      zts0 (:,:,:,:) = tsn   (:,:,:,:)
       ztmask1(:,:,:) = ptmask_b(:,:,:)
       ztmask0(:,:,:) = ptmask_b(:,:,:)
       DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
-          DO jk = 1,jpk-1
-              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
-              DO jj = 2,jpj-1
-                 DO ji = fs_2,fs_jpim1
-                      jip1=ji+1; jim1=ji-1;
-                      jjp1=jj+1; jjm1=jj-1;
-                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
-                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
-                      !! horizontal basic extrapolation
-                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
-                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
-                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
-                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
-                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
-                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
-                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
-                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
-                         ztmask1(ji,jj,jk)=1
-                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
-                      !! vertical extrapolation if horizontal extrapolation failed
-                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
-                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
-                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
-                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
-                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
-                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
-                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
-                            ztmask1(ji,jj,jk)=1._wp
-                         END IF
-                      END IF
-                  END DO
-              END DO
-          END DO
-          
-          CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
-
-          ! update
-          zts0(:,:,:,:) = tsn(:,:,:,:)
-          ztmask0 = ztmask1
-
+         DO jk = 1,jpk-1
+            zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
+            DO jj = 2,jpj-1
+               DO ji = fs_2,fs_jpim1
+                  jip1=ji+1; jim1=ji-1;
+                  jjp1=jj+1; jjm1=jj-1;
+                  summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
+                  IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
+                     !! horizontal basic extrapolation
+                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
+                         &            +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
+                         &            +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
+                         &            +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
+                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
+                         &            +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
+                         &            +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
+                         &            +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
+                     ztmask1(ji,jj,jk)=1
+                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
+                     !! vertical extrapolation if horizontal extrapolation failed
+                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
+                     summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
+                     IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
+                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
+                            &            +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / summsk
+                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
+                            &            +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / summsk
+                        ztmask1(ji,jj,jk)=1._wp
+                     ENDIF
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
+         !
+         ! update
+         zts0(:,:,:,:) = tsn(:,:,:,:)
+         ztmask0 = ztmask1
+         !
       END DO
 
@@ -395,6 +408,6 @@
       ! set mbkt and mikt to 1 in thiese location
       WHERE (SUM(tmask,dim=3) == 0)
-         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
-         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
+         mbkt(:,:) = 1   ;   mbku(:,:)=1   ;   mbkv(:,:) = 1
+         mikt(:,:) = 1   ;   miku(:,:)=1   ;   mikv(:,:) = 1
       END WHERE
       ! -------------------------------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynadv.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynadv.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynadv.F90	(revision 10001)
@@ -121,4 +121,8 @@
       ENDIF
 
+!!gm 
+      IF ( ln_dynadv_vec )   CALL ctl_stop( 'Vector form not available in this branch' )
+!!gm
+
       ioptio = 0                      ! parameter control and set n_dynadv
       IF( ln_dynadv_OFF  ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_LIN_dyn   ;   ENDIF
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90	(revision 10001)
@@ -20,4 +20,6 @@
    !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
    !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
+   !!            4.0  !  2018-07  (G. Madec)   1- z-star (s-star) only  
+   !!                                          2- remove   dom_vvl_interpol  
    !!-------------------------------------------------------------------------
   
@@ -34,5 +36,5 @@
    USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
    USE domvvl         ! variable volume
-   USE bdy_oce   , ONLY: ln_bdy
+   USE bdy_oce , ONLY : ln_bdy
    USE bdydta         ! ocean open boundary conditions
    USE bdydyn         ! ocean open boundary conditions
@@ -92,5 +94,5 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -100,4 +102,5 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
+      REAL(wp), DIMENSION(jpi,jpj) ::   z_ssh_h0, zsshu, zsshv
       !!----------------------------------------------------------------------
       !
@@ -212,7 +215,7 @@
             !    => time filter + conservation correction (only at the first level)
             zcoef = rn_atfp * rn_Dt * r1_rho0
-
+            !
             e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
-
+            !
             IF ( ln_rnf ) THEN
                IF( ln_rnf_depth ) THEN
@@ -228,8 +231,8 @@
                   END DO
                ELSE
-                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
+                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:) )*tmask(:,:,1)
                ENDIF
             ENDIF
-
+            !
             IF ( ln_isf ) THEN   ! if ice shelf melting
                DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
@@ -245,8 +248,24 @@
             ENDIF
             !
+            !
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
-               ! Before filtered scale factor at (u/v)-points
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
+               !                                            !* ssh at u- and v-points)
+               DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
+                  zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
+                  zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+               END DO             ;   END DO      
+               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )
+               !
+               !
+               !                                            !* e3u and e3v 
+               z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:)           ! u-point
+               DO jk = 1, jpkm1
+                  e3u_b (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) )
+               END DO
+               z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:)           ! v-point
+               DO jk = 1, jpkm1
+                  e3v_b (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) )
+               END DO
+               !
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
@@ -266,7 +285,23 @@
                !
                ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
-               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
-               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
+               !
+               !                                            !* ssh at u- and v-points)
+               DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
+                  zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
+                  zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+               END DO             ;   END DO      
+               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )
+               !
+               !
+               !                                            !* e3u and e3v 
+               z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:)           ! u-point
+               DO jk = 1, jpkm1
+                  ze3u_f(:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) )
+               END DO
+               z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:)           ! v-point
+               DO jk = 1, jpkm1
+                  ze3u_f(:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) )
+               END DO
+               !
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
@@ -319,14 +354,17 @@
       !
       IF(.NOT.ln_linssh ) THEN
-         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
-         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
-         DO jk = 2, jpkm1
-            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
-         END DO
-         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
+         DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
+            zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
+            zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+         END DO             ;   END DO      
+         CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )
+         !
+         hu_b   (:,:) = hu_0(:,:) + zsshu(:,:)
+         hv_b   (:,:) = hv_0(:,:) + zsshv(:,:)
+         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
          r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
       ENDIF
       !
+      !                                   ! 
       un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
       ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90	(revision 10001)
@@ -458,5 +458,5 @@
          !                                ! Note also that now e3 scale factors are used as after one are not computed !
          !
-         CALL wrk_alloc(jpi,jpj,   z2d )
+         ALLOCATE( z2d(jpi,jpj) )
          z2d(:,:) = 0._wp
          DO jk = 1, jpkm1
@@ -498,4 +498,5 @@
          CALL lbc_lnk( z2d,'T', 1. )
          CALL iom_put( 'dispkevfo', z2d )
+         DEALLOCATE( z2d )
       ENDIF
       !
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/IOM/in_out_manager.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/IOM/in_out_manager.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/IOM/in_out_manager.F90	(revision 10001)
@@ -35,5 +35,5 @@
    INTEGER       ::   nn_time0         !: initial time of day in hhmm
    INTEGER       ::   nn_leapy         !: Leap year calendar flag (0/1 or 30)
-   INTEGER       ::   nn_istate        !: initial state output flag (0/1)
+   LOGICAL       ::   ln_istate        !: initial state output flag
    INTEGER       ::   nn_write         !: model standard output frequency
    INTEGER       ::   nn_stock         !: restart file frequency
@@ -79,5 +79,4 @@
    INTEGER       ::   ndate0                      !: initial calendar date aammjj
    INTEGER       ::   nleapy                      !: Leap year calendar flag (0/1 or 30)
-   INTEGER       ::   ninist                      !: initial state output flag (0/1)
    INTEGER       ::   nwrite                      !: model standard output frequency
    INTEGER       ::   nstock                      !: restart file frequency
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcice_cice.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcice_cice.F90	(revision 10001)
@@ -205,5 +205,5 @@
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
-         ENDDO
+         END DO
       ENDIF
 
@@ -214,6 +214,6 @@
             fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
             fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
       CALL lbc_lnk_multi( fr_iu , 'U', 1.,  fr_iv , 'V', 1. )
@@ -235,7 +235,7 @@
                !
                DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-                  e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-               ENDDO
+                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * r1_ht_0(:,:) )
+                  e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshb(:,:) * r1_ht_0(:,:) )
+               END DO
                e3t_a(:,:,:) = e3t_b(:,:,:)
                ! Reconstruction of all vertical scale factors at now and before time-steps
@@ -307,6 +307,6 @@
                ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
                                     + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
          CALL nemo2cice(ztmp,strax,'F', -1. )
 
@@ -317,6 +317,6 @@
                ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
                                     + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
          CALL nemo2cice(ztmp,stray,'F', -1. )
 
@@ -325,5 +325,5 @@
             DO jl=1,ncat
                ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
-            ENDDO
+            END DO
          ELSE
 ! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
@@ -332,8 +332,8 @@
             DO jj=1,jpj
                DO ji=1,jpi
-                  IF (fr_i(ji,jj).eq.0.0) THEN
+                  IF (fr_i(ji,jj) == 0._wp ) THEN
                      DO jl=1,ncat
-                        ztmpn(ji,jj,jl)=0.0
-                     ENDDO
+                        ztmpn(ji,jj,jl) = 0._wp
+                     END DO
                      ! This will then be conserved in CICE
                      ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
@@ -341,8 +341,8 @@
                      DO jl=1,ncat
                         ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
-                     ENDDO
+                     END DO
                   ENDIF
-               ENDDO
-            ENDDO
+               END DO
+            END DO
          ENDIF
          DO jl=1,ncat
@@ -366,5 +366,5 @@
             ENDIF
             CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
-         ENDDO
+         END DO
 
       ELSE IF (ksbc == jp_blk) THEN
@@ -437,6 +437,6 @@
          DO ji=1,jpi
             ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
       CALL nemo2cice(ztmp,uocn,'F', -1. )
 
@@ -445,6 +445,6 @@
          DO ji=1,jpim1
             ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
       CALL nemo2cice(ztmp,vocn,'F', -1. )
 
@@ -511,6 +511,6 @@
          DO ji=2,jpim1
             ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
       CALL lbc_lnk( ss_iou , 'U', -1. )
 
@@ -523,6 +523,6 @@
          DO ji=2,jpim1
             ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
       CALL lbc_lnk( ss_iov , 'V', -1. )
 
@@ -609,6 +609,6 @@
          DO ji=1,jpi
             nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
 #if defined key_cice4
@@ -627,5 +627,5 @@
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
-         ENDDO
+         END DO
       ENDIF
 
@@ -636,6 +636,6 @@
             fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
             fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
       CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. )
@@ -872,6 +872,6 @@
          DO ji=2,nx_block-1
             pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
 #else
@@ -898,12 +898,12 @@
                DO ji=nldit(jn),nleit(jn)
                   png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
-               ENDDO
-            ENDDO
-         ENDDO
+               END DO
+            END DO
+         END DO
          DO jj=1,ny_global
             DO ji=1,nx_global
                pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
 
@@ -999,6 +999,6 @@
          DO ji=1,jpim1
             pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
 #else
@@ -1021,7 +1021,7 @@
                DO ji=nldit(jn),nleit(jn)
                   png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
-               ENDDO
-            ENDDO
-         ENDDO
+               END DO
+            END DO
+         END DO
       ENDIF
 
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90	(revision 10001)
@@ -175,4 +175,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
+!!gm BUG   here e3w_a  should be used !!!!!   but then should be added in the system 
                      zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
                      zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/nemogcm.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/nemogcm.F90	(revision 10001)
@@ -36,4 +36,6 @@
    !!----------------------------------------------------------------------
    !!   nemo_gcm      : solve ocean dynamics, tracer, biogeochemistry and/or sea-ice
+   !!   nemo_MLF      : Modified Leap-Frog time stepping loop
+   !!   nemo_RK3      : 3rd order Runge-Kutta time stepping loop
    !!   nemo_init     : initialization of the NEMO system
    !!   nemo_ctl      : initialisation of the contol print
@@ -60,4 +62,5 @@
    USE diacfl         ! CFL diagnostics               (dia_cfl_init routine)
    USE step           ! NEMO time-stepping                 (stp     routine)
+   USE stpRK3         ! NEMO time-stepping                 (stp     routine)
    USE icbini         ! handle bergs, initialisation
    USE icbstp         ! handle bergs, calving, themodynamics and transport
@@ -123,6 +126,4 @@
       !!              Madec, 2008, internal report, IPSL.
       !!----------------------------------------------------------------------
-      INTEGER ::   istp   ! time step index
-      !!----------------------------------------------------------------------
       !
 #if defined key_agrif
@@ -152,5 +153,49 @@
       !                            !-----------------------!
       !
-      !                                               !== set the model time-step  ==!
+      IF( ln_MLF )   CALL nemo_MLF        ! Modified Leap-Frog
+      !
+      IF( ln_RK3 )   CALL nemo_RK3        ! 3rd order Runge-Kutta
+      !
+      !
+      IF( ln_diaobs   )   CALL dia_obs_wri
+      !
+      IF( ln_icebergs )   CALL icb_end( nitend )
+
+      !                            !------------------------!
+      !                            !==  finalize the run  ==!
+      !                            !------------------------!
+      IF(lwp) WRITE(numout,cform_aaa)        ! Flag AAAAAAA
+      !
+      IF( nstop /= 0 .AND. lwp ) THEN        ! error print
+         WRITE(numout,cform_err)
+         WRITE(numout,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found'
+         WRITE(numout,*)
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_finalize
+      !
+      CALL nemo_closefile
+      !
+#if defined key_iomput
+                                    CALL xios_finalize  ! end mpp communications with xios
+      IF( lk_oasis     )            CALL cpl_finalize   ! end coupling and mpp communications with OASIS
+#else
+      IF    ( lk_oasis ) THEN   ;   CALL cpl_finalize   ! end coupling and mpp communications with OASIS
+      ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop        ! end mpp communications
+      ENDIF
+#endif
+      !
+   END SUBROUTINE nemo_gcm
+
+
+   SUBROUTINE nemo_MLF
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_MLF  ***
+      !!
+      !! ** Purpose :   Modified Leap-Frog time stepping loop
+      !!----------------------------------------------------------------------
+      INTEGER ::   istp   ! time step index
+      !!----------------------------------------------------------------------
+      !                                               !== set the MLF time-step  ==!
       !
       IF( l_1st_euler ) THEN   ;   rDt =         rn_Dt   ;   l_1st_euler = .TRUE.    ! start or restart with Euler 1st time-step
@@ -212,33 +257,50 @@
 #endif
       !
-      IF( ln_diaobs   )   CALL dia_obs_wri
-      !
-      IF( ln_icebergs )   CALL icb_end( nitend )
-
-      !                            !------------------------!
-      !                            !==  finalize the run  ==!
-      !                            !------------------------!
-      IF(lwp) WRITE(numout,cform_aaa)        ! Flag AAAAAAA
-      !
-      IF( nstop /= 0 .AND. lwp ) THEN        ! error print
-         WRITE(numout,cform_err)
-         WRITE(numout,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found'
-         WRITE(numout,*)
-      ENDIF
-      !
-      IF( ln_timing )   CALL timing_finalize
-      !
-      CALL nemo_closefile
-      !
-#if defined key_iomput
-                                    CALL xios_finalize  ! end mpp communications with xios
-      IF( lk_oasis     )            CALL cpl_finalize   ! end coupling and mpp communications with OASIS
+   END SUBROUTINE nemo_MLF
+
+
+   SUBROUTINE nemo_RK3
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_RK3  ***
+      !!
+      !! ** Purpose :   3rd order Runge-Kutta time stepping loop
+      !!----------------------------------------------------------------------
+      INTEGER ::   istp   ! time step index
+      !!----------------------------------------------------------------------
+      !
+      !                                         !== set the model time-step  ==!
+      rDt = rn_Dt
+      r1_Dt = 1._wp / rDt
+      istp = nit000
+      !
+#if defined key_agrif
+      !                                         !==  AGRIF time-stepping  ==!
+      CALL Agrif_Regrid()
+      !
+      CALL Agrif_step_child_adj(Agrif_Update_All)        ! Recursive update from highest to lowest nested levels
+      !
+      DO WHILE( istp <= nitend .AND. nstop == 0 )
+         CALL stp_RK3
+         istp = istp + 1
+      END DO
+      !
+      IF( .NOT. Agrif_Root() ) THEN
+         CALL Agrif_ParentGrid_To_ChildGrid()
+         IF( ln_diaobs )   CALL dia_obs_wri
+         IF( ln_timing )   CALL timing_finalize
+         CALL Agrif_ChildGrid_To_ParentGrid()
+      ENDIF
+      !
 #else
-      IF    ( lk_oasis ) THEN   ;   CALL cpl_finalize   ! end coupling and mpp communications with OASIS
-      ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop        ! end mpp communications
-      ENDIF
-#endif
-      !
-   END SUBROUTINE nemo_gcm
+      !                                         !==  Standard time-stepping  ==!
+      !
+      DO WHILE( istp <= nitend .AND. nstop == 0 )
+         CALL stp_RK3( istp ) 
+         istp = istp + 1
+      END DO
+      !
+#endif
+      !
+   END SUBROUTINE nemo_RK3
 
 
@@ -459,5 +521,12 @@
 
       !                                      ! Assimilation increments
-      IF( lk_asminc    )   CALL asm_inc_init    ! Initialize assimilation increments
+      IF( lk_asminc    ) THEN
+                           CALL asm_inc_init    ! Initialize assimilation increments
+                           l_dynasm = ln_asmiau .AND. ln_dyninc
+                           l_traasm = ln_asmiau .AND. ln_trainc
+      ELSE
+                           l_dynasm = .FALSE.
+                           l_traasm = .FALSE.
+      ENDIF
       !
       IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90	(revision 10001)
@@ -34,5 +34,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b   ,  un_b  ,  ua_b  !: Barotropic velocities at u-point [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vb_b   ,  vn_b  ,  va_b  !: Barotropic velocities at v-point [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   ,  sshn  ,  ssha  !: sea surface height at t-point [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   ,  sshn  , ssha   !: sea surface height at t-point [m]
+
+!!gm??   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ssh  , sshu   , sshv   !: sea surface height at t-, u- v-points [m]
+
 
    !! Arrays at barotropic time step:                   ! befbefore! before !  now   ! after  !
@@ -90,5 +93,5 @@
          &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) )
          !
-      ALLOCATE( sshb(jpi,jpj)     , sshn(jpi,jpj)   , ssha(jpi,jpj)   ,     &
+      ALLOCATE( sshb  (jpi,jpj)  , sshn   (jpi,jpj) , ssha(jpi,jpj)   ,     &
          &      ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     &
          &      vb_b(jpi,jpj)     , vn_b(jpi,jpj)   , va_b(jpi,jpj)   ,     &
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/par_oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/par_oce.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/par_oce.F90	(revision 10001)
@@ -12,4 +12,12 @@
    IMPLICIT NONE
    PUBLIC
+
+
+   !!
+   INTEGER, PUBLIC ::   Nnn, Np1, Nm1   ! =now, before, after
+
+   INTEGER, PUBLIC ::   Nb, Nn, Na      ! before, now, after  index
+
+
 
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step.F90	(revision 10001)
@@ -179,6 +179,5 @@
                          va(:,:,:) = 0._wp
 
-      IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
-               &         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment
+      IF( l_dynasm   )   CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment
       IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends
 #if defined key_agrif
@@ -232,6 +231,5 @@
                          tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero
 
-      IF(  lk_asminc .AND. ln_asmiau .AND. &
-         & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
+      IF(  l_traasm  )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
                          CALL tra_sbc       ( kstp )  ! surface boundary condition
       IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step_oce.F90	(revision 10000)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step_oce.F90	(revision 10001)
@@ -108,4 +108,5 @@
 #if defined key_top
    USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
+   USE trcadv           ! passive tracer time-stepping      (trc_stp routine)
 #endif
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/stpRK3.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/stpRK3.F90	(revision 10001)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/stpRK3.F90	(revision 10001)
@@ -0,0 +1,376 @@
+MODULE stpRK3
+   !!======================================================================
+   !!                       ***  MODULE stpRK3  ***
+   !! RK3 Time-stepping: manager of the ocean, tracer and ice time stepping
+   !!======================================================================
+   !! History :  5.0  !  2018-07  (G. Madec)  original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_RK3       : NEMO system RK3 time-stepping scheme
+   !!   stp_RK3_init  : initialize the RK3 scheme
+   !!----------------------------------------------------------------------
+   USE step_oce       ! time stepping definition modules
+   !
+   USE iom            ! xIOs server
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   stp_RK3   ! called by nemogcm.F90
+   
+   LOGICAL ::   l_1st_stg = .TRUE.     ! 1st stage only flag
+   LOGICAL ::   l_2nd_stg = .TRUE.     ! 2nd stage only flag
+   LOGICAL ::   l_3rd_stg = .TRUE.     ! 3rd stage only flag
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 5.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if defined key_agrif
+   RECURSIVE SUBROUTINE stp_RK3( )
+      INTEGER             ::   kstp   ! ocean time-step index
+#else
+   SUBROUTINE stp_RK3( kstp )
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
+#endif
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp  ***
+      !!
+      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.)
+      !!              - Time stepping of SI3 (dynamic and thermodynamic eqs.)
+      !!              - Time stepping of TRC  (passive tracer eqs.)
+      !!
+      !! ** Method  : -1- Update forcings and data
+      !!              -2- Update ocean physics
+      !!              -3- Compute the t and s trends
+      !!              -4- Update t and s
+      !!              -5- Compute the momentum trends
+      !!              -6- Update the horizontal velocity
+      !!              -7- Compute the diagnostics variables (rd,N2, hdiv,w)
+      !!              -8- Outputs and diagnostics
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk, jstg   ! dummy loop indice
+      INTEGER ::   indic              ! error indicator if < 0
+!!gm kcall can be removed, I guess
+      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
+      !! ---------------------------------------------------------------------
+#if defined key_agrif
+      kstp = nit000 + Agrif_Nb_Step()
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
+      ENDIF
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
+# if defined key_iomput
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
+# endif
+#endif
+      !
+      IF( ln_timing )   CALL timing_start('stp')
+
+!!!======================!!!
+!!!   First STAGE only   !!!
+!!!======================!!!
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! update I/O and calendar 
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                             indic = 0                ! reset to no error condition
+                             
+      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
+                             CALL iom_init(      cxios_context          )  ! for model grid (including passible AGRIF zoom)
+         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
+      ENDIF
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
+      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_tide    )   CALL sbc_tide( kstp )                   ! update tide potential
+      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
+      IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
+      IF ( ln_diurnal )  CALL stp_diurnal( kstp )                ! diagnose cool skin
+      !
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update stochastic parameters and random T/S fluctuations
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      IF( ln_sto_eos )   CALL sto_par( kstp )          ! Stochastic parameters
+      IF( ln_sto_eos )   CALL sto_pts( tsn  )          ! Random T/S fluctuations
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !  THERMODYNAMICS
+                         CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
+                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+
+      !  VERTICAL PHYSICS
+                         CALL zdf_phy( kstp )         ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
+
+      !  LATERAL  PHYSICS
+      !
+      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
+                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
+
+         IF( ln_zps .AND. .NOT. ln_isfcav)                               &
+            &            CALL zps_hde    ( kstp, 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( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
+         IF( ln_traldf_triad ) THEN 
+                         CALL ldf_slp_triad( kstp )                       ! before slope for triad operator
+         ELSE     
+                         CALL ldf_slp     ( kstp, rhd, rn2b )             ! before slope for standard operator
+         ENDIF
+      ENDIF
+      !                                                                   ! eddy diffusivity coeff.
+      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp )       !       and/or eiv coeff.
+      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp )       ! eddy viscosity coeff. 
+
+
+!!!======================!!!
+!!!   Loop over stages   !!!
+!!!======================!!!
+
+      DO jstg = 1, 3
+
+         SELECT CASE( jstg )
+         CASE( 1 )   ;   rDt = rn_Dt / 3._wp
+         CASE( 2 )   ;   rDt = rn_Dt / 2._wp
+         CASE( 3 )   ;   rDt = rn_Dt
+         END SELECT
+         
+         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+         !  Ocean dynamics : hdiv, ssh, e3, u, v, w
+         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+                               CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
+         IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
+                               CALL wzv           ( kstp )  ! now cross-level velocity 
+                               CALL eos    ( tsb, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
+                            
+         !!jc: fs simplification
+         !!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636) 
+         !!                                         but ensures reproductible results
+         !!                                         with previous versions using split-explicit free surface          
+         IF( ln_zps .AND. .NOT. ln_isfcav )                               &
+            &            CALL zps_hde    ( kstp, jpts, tsn, 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( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
+         !!jc: fs simplification
+
+                        ua (:,:,:)  = 0._wp          ! set the RHS of dyn Eq. to zero
+                        va (:,:,:)  = 0._wp
+                        tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero
+         !
+         !                          ! ================ !    
+         IF ( jstg <= 3 ) THEN      !   stages 1 & 2   :   ADV, COR, HPG and SPG trends only
+            !                       ! ================ !
+            !
+            !                             !==  dynamics  ==!
+                        CALL dyn_adv( kstp )    ! advection (vector or flux form)
+                        CALL dyn_vor( kstp )    ! vorticity term including Coriolis
+                        CALL dyn_hpg( kstp )    ! horizontal gradient of Hydrostatic pressure
+                        CALL dyn_spg( kstp )    ! surface pressure gradient
+                                                ! With split-explicit free surface, since now transports have been updated and ssha as well
+            IF( ln_dynspg_ts ) THEN             ! vertical scale factors and vertical velocity need to be updated
+                        CALL div_hor( kstp )          ! Horizontal divergence  (2nd call in time-split case)
+               IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
+                        CALL wzv    ( kstp )    ! now cross-level velocity 
+         ENDIF
+!!gm  to be added here :  time stepping ==>>>   un & vn
+
+
+            !                             !==  tracers  ==!
+#if defined key_top
+                        CALL trc_adv( kstp )    ! horizontal & vertical advection
+#endif
+                        CALL tra_adv( kstp )    ! horizontal & vertical advection
+
+!!gm  to be added here :  time stepping ==>>>   tsn
+                        
+            !
+            !                       ! ================ !
+         ELSE                       !     stage 3      :   add all dynamical trends 
+            !                       ! ================ !
+            !
+            !                             !==  dynamics  ==!
+                        CALL dyn_adv( kstp )  ! advection (vector or flux form)
+                        CALL dyn_vor( kstp )  ! vorticity term including Coriolis
+                        CALL dyn_hpg( kstp )  ! horizontal gradient of Hydrostatic pressure
+            !
+            IF(l_dynasm)   CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment
+            IF( ln_bdy )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends
+#if defined key_agrif
+            IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momentum sponge
+#endif
+            IF( ln_zdfosm )   CALL dyn_osm( kstp )  ! OSMOSIS non-local velocity fluxes
+                        CALL dyn_ldf( kstp )  ! lateral mixing
+                        CALL dyn_spg( kstp )    ! surface pressure gradient
+                                                ! With split-explicit free surface, since now transports have been updated and ssha as well
+            IF( ln_dynspg_ts ) THEN             ! vertical scale factors and vertical velocity need to be updated
+                        CALL div_hor( kstp )          ! Horizontal divergence  (2nd call in time-split case)
+               IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
+                        CALL wzv    ( kstp )    ! now cross-level velocity 
+            ENDIF
+                        CALL dyn_zdf( kstp )    ! vertical diffusion  & time-stepping
+
+            !                             !==  tracers  ==!
+            IF( ln_crs )   CALL crs_fld( kstp )   ! ocean model: online field coarsening & output
+#if defined key_top
+                         CALL trc_stp       ( kstp )  ! time-stepping
+#endif
+
+                         tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero
+
+                         CALL tra_adv       ( kstp )  ! horizontal & vertical advection
+      IF(  l_traasm  )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
+                         CALL tra_sbc       ( kstp )  ! surface boundary condition
+      IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
+      IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux
+      IF( ln_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme
+      IF( ln_tradmp  )   CALL tra_dmp       ( kstp )  ! internal damping trends
+      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends
+#if defined key_agrif
+      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_tra        ! tracers sponge
+#endif
+      IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes
+      IF( lrst_oce .AND. ln_zdfosm ) &
+           &             CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
+                         CALL tra_ldf       ( kstp )  ! lateral mixing
+
+!!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
+      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics
+!!gm
+                         CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields
+      IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection
+
+         
+         ENDIF
+
+
+
+
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Set boundary conditions and Swap
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap 
+!!    (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields. 
+!!    If so: 
+!!    (i) no need to call agrif update at initialization time
+!!    (ii) no need to update "before" fields 
+!!
+!!    Apart from creating new tra_swp/dyn_swp routines, this however: 
+!!    (i) makes boundary conditions at initialization time computed from updated fields which is not the case between 
+!!    two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation", 
+!!    e.g. a shift of the feedback interface inside child domain. 
+!!    (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
+!!    place.
+!! 
+!!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
+                         CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap
+                         CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
+                         CALL ssh_swp       ( kstp )  ! swap of sea surface height
+      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
+      !
+
+
+
+
+!!!==========================!!!
+!!!   end Loop over stages   !!!
+!!!==========================!!!
+
+      END DO
+
+
+      IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics
+
+!!gm : This does not only concern the dynamics ==>>> add a new title
+!!gm2: why ouput restart before AGRIF update?
+!!
+!!jc: That would be better, but see comment above
+!!
+      IF( lrst_oce   )   CALL rst_write    ( kstp )   ! write output ocean restart file
+      IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
+
+#if defined key_agrif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! AGRIF
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
+                         CALL Agrif_Integrate_ChildGrids( stp )  ! allows to finish all the Child Grids before updating
+
+                         IF( Agrif_NbStepint() == 0 ) CALL Agrif_update_all( ) ! Update all components
+#endif
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diagnostics and outputs
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_diaobs  )   CALL dia_obs ( kstp )        ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+      IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats
+      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics
+      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth)
+      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports
+                         CALL dia_ar5 ( kstp )        ! ar5 diag
+      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
+                         CALL dia_wri ( kstp )        ! ocean model: outputs
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Control
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL stp_ctl      ( kstp, indic )
+                         
+      IF( kstp == nit000 ) THEN                          ! 1st time step only
+                                        CALL iom_close( numror )   ! close input  ocean restart file
+         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
+         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Coupled mode
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!gm why lk_oasis and not lk_cpl ????
+      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
+      !
+#if defined key_iomput
+      IF( kstp == nitend .OR. indic < 0 ) THEN 
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
+      ENDIF
+#endif
+      !
+      IF( ln_timing )   CALL timing_stop('stp')
+      !
+   END SUBROUTINE stp_RK3
+
+
+   SUBROUTINE stp_RK3_init
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE stp_RK3_init  ***
+      !!
+      !! ** Purpose :   RK3 time stepping initialization
+      !!
+      !! ** Method  : 
+      !!----------------------------------------------------------------------
+      
+      
+   END SUBROUTINE stp_RK3_init
+   
+   !!======================================================================
+END MODULE stpRK3
