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 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90	(revision 10023)
@@ -127,12 +127,15 @@
    LOGICAL, PUBLIC ::   ln_sco       !: s-coordinate or hybrid z-s coordinate
    LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF 
+   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e3t, e3u, e3v, e3w        !: t-, u-, v-, w-points vertical scale factors [m]
+   
    !                                                        !  ref.   ! before  !   now   ! after  !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n            !: w- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
 
    !                                                        !  ref.   ! before  !   now   !
@@ -171,7 +174,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book)
 
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level             (ISF)
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level (ISF)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                    (ISF)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_isf                  !: Iceshelf draft                    (ISF)
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts
@@ -264,5 +266,7 @@
          &      e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) ,                      e3w_b(jpi,jpj,jpk) ,   & 
          &      e3t_n(jpi,jpj,jpk) , e3u_n(jpi,jpj,jpk) , e3v_n(jpi,jpj,jpk) , e3f_n(jpi,jpj,jpk) , e3w_n(jpi,jpj,jpk) ,   & 
-         &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) ,                                             &
+         &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk)                      , e3w_a(jpi,jpj,jpk) ,   &
+         !
+         &      e3t(jpi,jpj,jpk,Nt) , e3u(jpi,jpj,jpk,Nt) , e3v(jpi,jpj,jpk,Nt) ,                                          &
          !                                                          !
          &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         &
@@ -282,6 +286,6 @@
          &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) )
          !
-      ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) ,     &
-         &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) )
+      ALLOCATE( ht_isf(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) ,     &
+         &                        mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) )
          !
       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     & 
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 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90	(revision 10023)
@@ -141,5 +141,5 @@
          DO ji = 1, jpi
             ik = mikt(ji,jj)
-            risfdep(ji,jj) = gdepw_0(ji,jj,ik)        !!gm  RENAME it as h_isf(:,:)  better no?
+            ht_isf(ji,jj) = gdepw_0(ji,jj,ik)
          END DO
       END DO
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 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90	(revision 10023)
@@ -9,4 +9,5 @@
    ! 1- remove   z-tilde          ==>>>  pure z-star (or s-star)
    ! 2- remove   dom_vvl_interpol  
+   ! 3- 
    
    !!----------------------------------------------------------------------
@@ -142,23 +143,24 @@
       !                          !* BEFORE fields : 
       CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
-      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw        (from 1 to jpkm1)
+      !                                !      e3t, e3u , e3v              (from 1 to jpkm1)
+      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
+      !                                !      gdept, gdepw                (from 1 to jpk  )
       !
       !                                ! set jpk level one to the e3._0 values
-      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)
+      e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3u_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
       !
       !                          !* NOW fields : 
       CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
-      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f   (from 1 to jpkm1)
-      !                                !      gdept_n, gdepw_n, gde3w_n
-!!gm issue?   gdept_n, gdepw_n, gde3w_n never defined at jpk
+      !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
+      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
+      !                                !      gdept, gdepw, gde3w         (from 1 to jpk  )
       !
       !                                ! 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)
+      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3u_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_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(:,:,:)
+      e3t_a(:,:,:) = e3t_n(:,:,:)   ;   e3u_a(:,:,:) = e3u_n(:,:,:)
+      e3w_a(:,:,:) = e3w_n(:,:,:)   ;   e3v_a(:,:,:) = e3v_n(:,:,:)
       
 !!gm           
@@ -207,4 +209,29 @@
       !!
       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
+      !!
+      ! !  ref.   ! before  !   now   ! after  !
+      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
+      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
+      !
+      ! !  ref.   ! before  !   now   !
+      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
+      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
+      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
+      ! 
+      ! !  ref. ! before  !   now   !  after  !
+      !   ht_0            ,    ht_n             !: t-depth              [m]
+      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+      !   hf_0                                  !: v-depth              [m]
+      ! r1_ht_0                                 !: inverse of u-depth [1/m]
+      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+      ! r1_hf_0                                 !: inverse of v-depth [1/m]
+      !
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in )           ::   kt      ! time step
@@ -230,8 +257,10 @@
       !
       !                                   !==  after ssh  ==!  (u- and v-points)
-      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-         zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa) ) * ssvmask(ji,jj)
-      END DO             ;   END DO      
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa) ) * ssvmask(ji,jj)
+         END DO
+      END DO      
       CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
       !
@@ -247,7 +276,8 @@
       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) )
+         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)                     )
+         e3w_a(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
       END DO
       !
@@ -315,7 +345,8 @@
          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
+         e3t_n  (:,:,jk) = e3t_a  (:,:,jk)         ! e3t, e3u, e3v, e3w
          e3u_n  (:,:,jk) = e3u_a  (:,:,jk)
          e3v_n  (:,:,jk) = e3v_a  (:,:,jk)
+         e3w_n  (:,:,jk) = e3w_a  (:,:,jk)
       END DO
       ht_n(:,:) = ht_0(:,:) + ssh(:,:,Nnn)            ! ocean thickness
@@ -324,5 +355,5 @@
       hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
       !
-      !                    !==  before : 
+      !                    !==  before  ==!
       !                                            !* ssh at u- and v-points)
       DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
@@ -341,6 +372,20 @@
          e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
       END DO
+      ! 
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)          !* gdept , gdepw
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpkm1
+            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+         END DO
+      ELSE                    ! no ISF cavities 
+         DO jk = 1, jpkm1
+            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+         END DO
+      ENDIF
       !      
-      !                    !==  now    : 
+      !                    !==   now    ==!
       !                                            !* ssh at u- and v-points)
       DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 for f-point
@@ -358,5 +403,4 @@
       zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)     ! f-point
       DO jk = 1, jpkm1
-          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
          e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
          e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
@@ -364,10 +408,10 @@
       END DO      
       ! 
-      zssht_h(:,:) = 1._wp + zssht_h(:,:)               ! t-point
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)          !* gdept , gdepw , gde3w
       !
       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(:,:)
+            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
             gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
          END DO
@@ -584,4 +628,28 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ssh2e3_now  ***
+      ! !  ref.   ! before  !   now   ! after  !
+      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
+      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
+      !
+      ! !  ref.   ! before  !   now   !
+      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
+      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
+      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
+      ! 
+      ! !  ref. ! before  !   now   !  after  !
+      !   ht_0            ,    ht_n             !: t-depth              [m]
+      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+      !   hf_0                                  !: v-depth              [m]
+      ! r1_ht_0                                 !: inverse of u-depth [1/m]
+      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+      ! r1_hf_0                                 !: inverse of v-depth [1/m]
+      !
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk
@@ -616,19 +684,24 @@
       zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)
       !
-      !                             !==  e3t, e3w  ,  e3u, e3uw ,  e3v, e3vw  , and e3f  ==!
+      !                             !==  e3t  ,  e3u  ,  e3v  ,  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(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
-          !
-          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) )
-          !
+          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
+          e3v_n(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
           e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
       END DO
       !      
+      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
+      !
+       e3w_n(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
+      e3uw_n(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
+      e3vw_n(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      DO jk = 2, jpk
+          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
+         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
+         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
+      END DO
+      !
       !                             !== depth of t- and w-points  ==!
       !
@@ -636,12 +709,11 @@
       !
       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(:,:)
+         DO jk = 1, jpk
+            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
             gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh(:,:,Nnn)
          END DO
       ELSE                    ! no ISF cavities
-!!gm BUG ???    gdept should be updated down to the ocean floor !  ===>> jpk NOT jpkm1 !!!
-         DO jk = 1, jpkm1
+         DO jk = 1, jpk
             gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
             gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
@@ -656,4 +728,28 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ssh2e3_before  ***
+      ! !  ref.   ! before  !   now   ! after  !
+      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
+      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
+      !
+      ! !  ref.   ! before  !   now   !
+      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
+      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
+      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
+      ! 
+      ! !  ref. ! before  !   now   !  after  !
+      !   ht_0            ,    ht_n             !: t-depth              [m]
+      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+      !   hf_0                                  !: v-depth              [m]
+      ! r1_ht_0                                 !: inverse of u-depth [1/m]
+      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+      ! r1_hf_0                                 !: inverse of v-depth [1/m]
+      !
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk
@@ -677,21 +773,44 @@
       !
       !      
-      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
+      !                             !==  ssh / h  factor at t-, u- ,v-points  ==!
       zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)
       zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
       zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
       !
-      !                             !==  e3t, e3w  ,  e3u, e3uw , and  e3v, e3vw  ==!
+      !                             !==  e3t  ,  e3u  ,  e3v  ==!
+      !      
       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(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
-          !
-          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) )
+          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
+          e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
+      END DO
+      !      
+      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
+      !
+       e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
+      e3uw_b(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
+      e3vw_b(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      DO jk = 2, jpk
+          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
+         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
+         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
       END DO
       !   
+      !                             !== depth of t- and w-points  ==!
+      !
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpk
+            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+         END DO
+      ELSE                    ! no ISF cavities
+         DO jk = 1, jpk
+            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+         END DO
+      ENDIF
+      !
    END SUBROUTINE ssh2e3_before
    
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl_RK3.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl_RK3.F90	(revision 10023)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl_RK3.F90	(revision 10023)
@@ -0,0 +1,818 @@
+MODULE domvvl_RK3
+   !!======================================================================
+   !!                       ***  MODULE domvvl_RK3   ***
+   !! Ocean : 
+   !!======================================================================
+   !! 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  
+   ! 3- 
+   
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
+   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
+   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
+   !!   dom_vvl_rst      : read/write restart file
+   !!   dom_vvl_ctl      : Check the vvl options
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE phycst          ! physical constant
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! ocean surface boundary condition
+   USE wet_dry         ! wetting and drying
+   USE usrdef_istate   ! user defined initial state (wad only)
+   USE restart         ! ocean restart
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+   USE lib_mpp         ! distributed memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC  dom_vvl_init       ! called by domain.F90
+   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
+   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
+   PUBLIC  ssh2e3_before      ! ...
+   PUBLIC  ssh2e3_now         ! ...
+
+   !                                                      !!* Namelist nam_vvl
+   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
+   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
+   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
+   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
+   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
+   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
+   !                                                       ! conservation: not used yet
+   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
+   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
+   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
+   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
+   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td      ! thickness diffusion transport
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf           ! low frequency part of hz divergence
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n    ! baroclinic scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a   ! baroclinic scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t       ! retoring period for scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv       ! retoring period for low freq. divergence
+
+
+!!gm add
+!!gm 
+
+
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
+   !! $Id$
+   !! Software governed by the CeCILL licence     (./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION dom_vvl_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION dom_vvl_alloc  ***
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         ALLOCATE( te3t_b(jpi,jpj,jpk)  , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) ,   &
+            &      dte3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
+            &      STAT = dom_vvl_alloc        )
+         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
+         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
+         un_td = 0._wp
+         vn_td = 0._wp
+      ENDIF
+      IF( ln_vvl_ztilde ) THEN
+         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
+         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
+         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
+      ENDIF
+      !
+   END FUNCTION dom_vvl_alloc
+
+
+   SUBROUTINE dom_vvl_init
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_init  ***
+      !!                   
+      !! ** Purpose :  Initialization of all scale factors, depths
+      !!               and water column heights
+      !!
+      !! ** Method  :  - use restart file and/or initialize
+      !!               - interpolate scale factors
+      !!
+      !! ** Action  : - e3t_(n/b) and te3t_(n/b)
+      !!              - Regrid: e3(u/v)_n
+      !!                        e3(u/v)_b       
+      !!                        e3w_n           
+      !!                        e3(u/v)w_b      
+      !!                        e3(u/v)w_n      
+      !!                        gdept_n, gdepw_n and gde3w_n
+      !!              - h(t/u/v)_0
+      !!              - frq_rst_e3t and frq_rst_hdv
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      !
+      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
+      !
+      !                    ! Allocate module arrays
+      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
+      !
+      !                    ! Read or initialize ssh(Nbb) & ssh(Nnn)
+      CALL dom_vvl_rst( nit000, 'READ' )
+      !
+      !                    !== Set of all other vertical mesh fields  ==!  (now and before)      
+      !
+      !                          !* BEFORE fields : 
+      CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
+      !                                !      e3t, e3u , e3v              (from 1 to jpkm1)
+      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
+      !                                !      gdept, gdepw                (from 1 to jpk  )
+      !
+      !                                ! set jpk level one to the e3._0 values
+      e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3u_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
+      !
+      !                          !* NOW fields : 
+      CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
+      !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
+      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
+      !                                !      gdept, gdepw, gde3w         (from 1 to jpk  )
+      !
+      !                                ! set one for all last level to the e3._0 value
+      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3u_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_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(:,:,:)
+      e3w_a(:,:,:) = e3w_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')
+      ENDIF
+      !
+   END SUBROUTINE dom_vvl_init
+
+
+   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
+      !!                   
+      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
+      !!                 tranxt and dynspg routines
+      !!
+      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
+      !!               - z_tilde_case: after scale factor increment = 
+      !!                                    high frequency part of horizontal divergence
+      !!                                  + retsoring towards the background grid
+      !!                                  + thickness difusion
+      !!                               Then repartition of ssh INCREMENT proportionnaly
+      !!                               to the "baroclinic" level thickness.
+      !!
+      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
+      !!               - te3t_a: after increment of vertical scale factor 
+      !!                              in z_tilde case
+      !!               - e3(t/u/v)_a
+      !!
+      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
+      !!
+      ! !  ref.   ! before  !   now   ! after  !
+      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
+      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
+      !
+      ! !  ref.   ! before  !   now   !
+      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
+      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
+      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
+      ! 
+      ! !  ref. ! before  !   now   !  after  !
+      !   ht_0            ,    ht_n             !: t-depth              [m]
+      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+      !   hf_0                                  !: v-depth              [m]
+      ! r1_ht_0                                 !: inverse of u-depth [1/m]
+      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+      ! r1_hf_0                                 !: inverse of v-depth [1/m]
+      !
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )           ::   kt      ! time step
+      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
+      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+      ENDIF
+
+      !                             ! ------------------!
+      !                             ! 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 * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa) ) * 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(:,:) = ssh    (:,:,Naa) * 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)                     )
+         e3w_a(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
+      END DO
+      !
+      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
+      !
+   END SUBROUTINE dom_vvl_sf_nxt
+
+
+   SUBROUTINE dom_vvl_sf_swp( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_swp  ***
+      !!                   
+      !! ** Purpose :  compute time filter and swap of scale factors 
+      !!               compute all depths and related variables for next time step
+      !!               write outputs and restart file
+      !!
+      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
+      !!               - reconstruct scale factor at other grid points (interpolate)
+      !!               - recompute depths and water height fields
+      !!
+      !! ** Action  :  - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step
+      !!               - Recompute:
+      !!                    e3(u/v)_b       
+      !!                    e3w_n           
+      !!                    e3(u/v)w_b      
+      !!                    e3(u/v)w_n      
+      !!                    gdept_n, gdepw_n  and gde3w_n
+      !!                    h(u/v) and h(u/v)r
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
+      !
+      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
+      !
+      IF( kt == nit000 )   THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
+      ENDIF
+      !
+      ! 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
+      !
+      ! - 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, e3w
+         e3u_n  (:,:,jk) = e3u_a  (:,:,jk)
+         e3v_n  (:,:,jk) = e3v_a  (:,:,jk)
+         e3w_n  (:,:,jk) = e3w_a  (:,:,jk)
+      END DO
+      ht_n(:,:) = ht_0(:,:) + ssh(:,:,Nnn)            ! ocean thickness
+      !
+      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
+      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
+      !
+      !                    !==  before  ==!
+      !                                            !* ssh at u- and v-points)
+      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
+         zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * 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(:,:) = ssh    (:,:,Nbb) * 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(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
+         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
+      ! 
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)          !* gdept , gdepw
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpkm1
+            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+         END DO
+      ELSE                    ! no ISF cavities 
+         DO jk = 1, jpkm1
+            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+         END DO
+      ENDIF
+      !      
+      !                    !==   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 * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
+         zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
+            &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * 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(:,:) = ssh    (:,:,Nnn) * 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
+         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 + zssht_h(:,:)          !* gdept , gdepw , gde3w
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpkm1
+            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
+         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) - ssh    (:,:,Nnn)
+         END DO
+      ENDIF
+      !
+      ! write restart file
+      ! ==================
+      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
+      !
+      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
+      !
+   END SUBROUTINE dom_vvl_sf_swp
+
+
+   SUBROUTINE dom_vvl_rst( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dom_vvl_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write VVL file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!                if the restart does not contain vertical scale factors,
+      !!                they are set to the _0 values
+      !!                if the restart does not contain vertical scale factors increments (z_tilde),
+      !!                they are set to 0.
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !
+      INTEGER ::   ji, jj, jk
+      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+         !                                   ! ===============
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            CALL rst_read_open                  !  open the restart file if necessary
+            !
+            id1 = iom_varid( numror, 'sshb'      , ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'sshn'      , ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+               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'   , ssh(:,:,Nnn), ldxios = lrxios    )
+!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
+!!gm 
+               CALL iom_get( numror, jpdom_data, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
+               CALL iom_get( numror, jpdom_data, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
+!!gm end
+               IF( l_1st_euler )   ssh(:,:,Nbb) = ssh(:,:,Nnn)
+            ELSE IF( id1 > 0 ) THEN
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files'
+               IF(lwp) write(numout,*) '   set ssh(Nnn) = ssh(Nbb)  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', ssh(:,:,Nbb), ldxios = lrxios )
+               ssh(:,:,Nnn) = ssh(:,:,Nbb)
+               l_1st_euler = .TRUE.
+            ELSE IF( id2 > 0 ) THEN
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files'
+               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn)  and force l_1st_euler = true'
+               CALL iom_get( numror, jpdom_data, 'sshn', ssh(:,:,Nnn), ldxios = lrxios )
+               ssh(:,:,Nbb) = ssh(:,:,Nnn)
+               l_1st_euler = .TRUE.
+            ELSE
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file'
+               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn) = 0  and force l_1st_euler = true'
+               ssh(:,:,Nbb) = 0._wp
+               ssh(:,:,Nnn) = 0._wp
+               l_1st_euler = .TRUE.
+            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
+                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
+                  tsn(:,:,:,:) = tsb(:,:,:,:)            ! set now values from to before ones
+                  ssh(:,:,Nnn) = ssh(:,:,Nbb)
+                  un (:,:,:)   = ub (:,:,:)
+                  vn (:,:,:)   = vb (:,:,:)
+               ELSE                                   ! Not the test case
+                  ssh(:,:,Nnn) = -ssh_ref
+                  ssh(:,:,Nbb) = -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
+                           ssh(ji,jj,Nbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Nnn) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Naa) = rn_wdmin1 - (ht_0(ji,jj) )
+                        ENDIF
+                     END DO
+                  END DO
+               ENDIF     ! If test case else
+               !
+               DO ji = 1, jpi
+                  DO jj = 1, jpj
+                     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
+                  END DO 
+               END DO 
+               !
+            ELSE
+               !
+               ! 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, ssh(:,:,Nbb)  )
+               ssh(:,:,Nnn) = 0._wp
+               ssh(:,:,Nbb) = 0._wp
+               !
+            END IF
+            !
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! ===================
+
+!!gm      DO NOTHING,   ssh(Nbb) and ssh(Nnn)  are written in restart.F90
+
+      ENDIF
+      !
+   END SUBROUTINE dom_vvl_rst
+
+
+   SUBROUTINE dom_vvl_ctl
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_ctl  ***
+      !!                
+      !! ** Purpose :   Control the consistency between namelist options
+      !!                for vertical coordinate
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio, ios
+      !!
+      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
+         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
+         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
+      !!---------------------------------------------------------------------- 
+      !
+      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist : 
+      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
+      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_vvl )
+      !
+      IF(lwp) THEN                    ! Namelist print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
+         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
+         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
+         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
+         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
+         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
+         WRITE(numout,*) '      !'
+         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
+         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
+         IF( ln_vvl_ztilde_as_zstar ) THEN
+            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
+            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
+            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
+            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
+            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
+            WRITE(numout,*) '                         rn_lf_cutoff   = 1/rn_Dt'
+         ELSE
+            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
+            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
+         ENDIF
+         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
+      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.
+      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
+      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
+      IF( ln_vvl_layer           )   ioptio = ioptio + 1
+      !
+      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
+      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
+      !
+      IF(lwp) THEN                   ! Print the choice
+         WRITE(numout,*)
+         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
+         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
+         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
+         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
+      ENDIF
+      !
+#if defined key_agrif
+      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
+#endif
+      !
+   END SUBROUTINE dom_vvl_ctl
+
+   
+   SUBROUTINE ssh2e3_now
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ssh2e3_now  ***
+      ! !  ref.   ! before  !   now   ! after  !
+      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
+      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
+      !
+      ! !  ref.   ! before  !   now   !
+      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
+      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
+      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
+      ! 
+      ! !  ref. ! before  !   now   !  after  !
+      !   ht_0            ,    ht_n             !: t-depth              [m]
+      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+      !   hf_0                                  !: v-depth              [m]
+      ! r1_ht_0                                 !: inverse of u-depth [1/m]
+      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+      ! r1_hf_0                                 !: inverse of v-depth [1/m]
+      !
+      !!----------------------------------------------------------------------
+      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 * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
+            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
+               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * 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(:,:) +  ssh   (:,:,Nnn)
+      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(:,:) =  ssh   (:,:,Nnn) * 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  ,  e3u  ,  e3v  ,  e3f  ==!
+      !      
+      DO jk = 1, jpkm1
+          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
+          e3v_n(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
+          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
+      END DO
+      !      
+      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
+      !
+       e3w_n(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
+      e3uw_n(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
+      e3vw_n(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      DO jk = 2, jpk
+          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
+         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
+         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
+      END DO
+      !
+      !                             !== depth of t- and w-points  ==!
+      !
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpk
+            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh(:,:,Nnn)
+         END DO
+      ELSE                    ! no ISF cavities
+         DO jk = 1, jpk
+            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh(:,:,Nnn)
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE ssh2e3_now
+   
+   
+   SUBROUTINE ssh2e3_before
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ssh2e3_before  ***
+      ! !  ref.   ! before  !   now   ! after  !
+      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
+      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
+      !
+      ! !  ref.   ! before  !   now   !
+      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
+      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
+      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
+      ! 
+      ! !  ref. ! before  !   now   !  after  !
+      !   ht_0            ,    ht_n             !: t-depth              [m]
+      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
+      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+      !   hf_0                                  !: v-depth              [m]
+      ! r1_ht_0                                 !: inverse of u-depth [1/m]
+      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
+      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+      ! r1_hf_0                                 !: inverse of v-depth [1/m]
+      !
+      !!----------------------------------------------------------------------
+      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  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * 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-points  ==!
+      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)
+      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
+      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
+      !
+      !                             !==  e3t  ,  e3u  ,  e3v  ==!
+      !      
+      DO jk = 1, jpkm1
+          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
+          e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
+      END DO
+      !      
+      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
+      !
+       e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
+      e3uw_b(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
+      e3vw_b(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      DO jk = 2, jpk
+          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
+         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
+         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
+      END DO
+      !   
+      !                             !== depth of t- and w-points  ==!
+      !
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpk
+            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
+         END DO
+      ELSE                    ! no ISF cavities
+         DO jk = 1, jpk
+            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE ssh2e3_before
+   
+   !!======================================================================
+END MODULE domvvl_RK3
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domwri.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domwri.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domwri.F90	(revision 10023)
@@ -159,5 +159,5 @@
       zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
       CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
-      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
+      zprt(:,:) = ssmask(:,:) * REAL( ht_isf(:,:) , wp )
       CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
       !															             ! vertical mesh
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 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90	(revision 10023)
@@ -13,4 +13,5 @@
    USE oce             ! global tra/dyn variable
    USE dom_oce         ! ocean space and time domain
+   USE domvvl          !
    USE domwri          ! ocean space and time domain
    USE phycst          ! physical constants
@@ -149,5 +150,5 @@
       zwmaskn(:,:,1) = tmask   (:,:,1)
       zwmaskb(:,:,1) = ptmask_b(:,:,1)
-      DO jk = 2,jpk
+      DO jk = 2, jpk
          zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
          zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
@@ -179,9 +180,35 @@
       ssh(:,:,Nnn) = ssh(:,:,Nnn) * ssmask(:,:)
 
+!!gm BUGs....
+!
+!     for me ht_0, hu_0, hv_0 and hf_0 should be recomputed whatever the value of ln_linssh
+!     further more mask at all grid-point should be recomputed
+!     and mikt, u, v, f also...
+!
+!     perhaps, not, if dom_cfg.nc file has been modified....
+!
+!     Pierre we should discuss of that !
+
+
+
+
 !=============================================================================
 !PM: Is this needed since introduction of VVL by default?
       IF ( .NOT.ln_linssh ) THEN
-      ! Reconstruction of all vertical scale factors at now time steps
-      ! ======================================================================
+         ! Reconstruction of all vertical scale factors at now time steps
+         ! ======================================================================
+
+
+         CALL ctl_stop( 'iscplrst : gm: here there is a BUG: not all required fields are defined')
+
+         !
+         !                          !* NOW fields : 
+         CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
+         !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
+         !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
+         !                                !      gdept, gdepw, gde3w         (from 1 to jpk  )
+         !
+
+
       
 !!gm Question : bug ????
@@ -193,75 +220,4 @@
       ! Note that the former calculation were using ht_0  so if it as not been updated ===>>> BUG
 !!gm
-      
-      
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-         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 + ssh(ji,jj,Nnn) * 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  * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
-            zsshv(ji,jj) = 0.5_wp  * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
-            zsshf(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
-               &                     + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * 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(:,:) +  ssh (:,:,Nnn)
-         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(:,:) = ssh(:,:,Nnn) * 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 + ssh(:,:,Nnn) * 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) - ssh(:,:,Nnn)
-            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) - ssh(:,:,Nnn)
-            END DO
-         ENDIF
-         
       ENDIF
 
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/istate.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/istate.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/istate.F90	(revision 10023)
@@ -104,5 +104,9 @@
             !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb) )         
+!!gm            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb) )    
+!!gm  when ln_linssh=.FALSE. (non linear free surface, gdept_b  is NOT initialized !
+!!gm     
+            CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, ssh(:,:,Nbb) ) 
+!!gm        
          ENDIF
          tsn  (:,:,:,:) = tsb(:,:,:,:)       ! set now values from to before ones
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynhpg.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynhpg.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynhpg.F90	(revision 10023)
@@ -232,5 +232,5 @@
          !
          !                                !- compute rhd at the ice/oce interface (ice shelf side)
-         CALL eos( zts_top , risfdep, zrhdtop_isf )
+         CALL eos( zts_top , ht_isf, zrhdtop_isf )
          !
          !                                !- Surface value + ice shelf gradient
@@ -245,5 +245,5 @@
                END DO
                IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
-                  &                                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
+                  &                                              * ( ht_isf(ji,jj) - gdept_1d(ikt-1) )
             END DO
          END DO
@@ -599,5 +599,5 @@
         END DO
       END DO
-      CALL eos( zts_top, risfdep, zrhdtop_oce )
+      CALL eos( zts_top, ht_isf, zrhdtop_oce )
 
 !==================================================================================     
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 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90	(revision 10023)
@@ -255,5 +255,5 @@
                   zsshv(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
                END DO             ;   END DO      
-               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )
+               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshv(:,:),'V', 1._wp )
                !
                !
@@ -291,5 +291,5 @@
                   zsshv(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
                END DO             ;   END DO      
-               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )
+               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshv(:,:),'V', 1._wp )
                !
                !
@@ -301,5 +301,5 @@
                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) )
+                  ze3v_f(:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) )
                END DO
                !
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/sshwzv.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/sshwzv.F90	(revision 10023)
@@ -245,5 +245,5 @@
          Naa   = isave                 ! after <-- previously now index
          !
-      ELSE                       !==  Leap-Frog time-stepping  ==!   Asselin filter + swap
+      ELSEIF ( ln_MLF  ) THEN    !==  Leap-Frog time-stepping  ==!   Asselin filter + swap
          !
          !                                   ! before <-- now filtered
@@ -258,4 +258,10 @@
          Nnn   = Naa                   ! now    <-- after
          Naa   = isave                 ! after  <-- previously now index
+         !
+      ELSEIF ( ln_MLF  ) THEN    !==  RK3 time-stepping  ==!  swap
+         isave = Nnn
+         Nnn   = Naa                   ! now    <-- after
+         Naa   = isave                 ! after  <-- previously now index
+
       ENDIF
       !
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/LDF/ldfslp.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/LDF/ldfslp.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/LDF/ldfslp.F90	(revision 10023)
@@ -180,8 +180,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt  (ji,jj), hmlpt  (ji+1,jj  ), 5._wp) &
-                  &                                  - MAX(risfdep(ji,jj), risfdep(ji+1,jj  )       ) ) 
-               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt  (ji,jj), hmlpt  (ji  ,jj+1), 5._wp) &
-                  &                                  - MAX(risfdep(ji,jj), risfdep(ji  ,jj+1)       ) )
+               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX( hmlpt (ji,jj) , hmlpt (ji+1,jj  ) , 5._wp )  &
+                  &                                  - MAX( ht_isf(ji,jj) , ht_isf(ji+1,jj  )         )  ) 
+               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX( hmlpt (ji,jj) , hmlpt (ji  ,jj+1) , 5._wp )  &
+                  &                                  - MAX( ht_isf(ji,jj) , ht_isf(ji  ,jj+1)         )  )
             END DO
          END DO
@@ -189,6 +189,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
-               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
+               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX( hmlpt(ji,jj) , hmlpt(ji+1,jj  ) , 5._wp )
+               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX( hmlpt(ji,jj) , hmlpt(ji  ,jj+1) , 5._wp )
             END DO
          END DO
@@ -211,8 +211,8 @@
                zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
                ! thickness of water column between surface and level k at u/v point
-               zdepu = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji+1,jj,jk) )                            &
-                                - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj))   )
-               zdepv = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji,jj+1,jk) )                            &
-                                - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj))   )
+               zdepu = 0.5_wp * ( ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )                            &
+                                - 2 * MAX( ht_isf(ji,jj), ht_isf(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj))   )
+               zdepv = 0.5_wp * ( ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )                            &
+                                - 2 * MAX( ht_isf(ji,jj), ht_isf(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj))   )
                !
                zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps )                                     &
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcisf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcisf.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcisf.F90	(revision 10023)
@@ -517,5 +517,5 @@
          CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
             ! Calculate freezing temperature
-            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
+            CALL eos_fzp( stbl(:,:), zfrz(:,:), ht_isf(:,:) )
 
             ! compute gammat every where (2d)
@@ -543,15 +543,15 @@
                DO ji = 1, jpi
                   ! compute coeficient to solve the 2nd order equation
-                  zeps1 = rcp*rho0*zgammat(ji,jj)
-                  zeps2 = rLfus*rho0*zgammas(ji,jj)
-                  zeps3 = rho_isf*rcp_isf*rkappa/MAX(risfdep(ji,jj),zeps)
-                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
-                  zeps6 = zeps4-ttbl(ji,jj)
-                  zeps7 = zeps4-tsurf
-                  zaqe  = zlamb1 * (zeps1 + zeps3)
-                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
-                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
-                  zcqe  = zeps2*stbl(ji,jj)
-                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
+                  zeps1 = rho0_rcp   * zgammat(ji,jj)
+                  zeps2 = rho0*rLfus * zgammas(ji,jj)
+                  zeps3 = rho_isf*rcp_isf*rkappa / MAX( ht_isf(ji,jj) , zeps )
+                  zeps4 = zlamb2 + zlamb3 * ht_isf(ji,jj)
+                  zeps6 = zeps4 - ttbl(ji,jj)
+                  zeps7 = zeps4 - tsurf
+                  zaqe  = zlamb1 * ( zeps1 + zeps3 )
+                  zaqer = 0.5_wp / MIN( zaqe , -zeps )
+                  zbqe  = zeps1 * zeps6 + zeps3 * zeps7 - zeps2
+                  zcqe  = zeps2 * stbl(ji,jj)
+                  zdis  = zbqe * zbqe - 4.0_wp * zaqe * zcqe               
 
                   ! Presumably zdis can never be negative because gammas is very small compared to gammat
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trasbc.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trasbc.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trasbc.F90	(revision 10023)
@@ -122,10 +122,11 @@
       ENDIF
       !                             !==  Now sbc tracer content fields  ==!
-      DO jj = 2, jpj
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
-               IF ( ssh(ji,jj,Nnn) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN
+      !
+      IF ( ll_wd ) THEN                   !* WAD case:  If near WAD point limit the flux for now
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               IF     ( ssh(ji,jj,Nnn) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN
                   sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
-               ELSE IF ( ssh(ji,jj,Nnn) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
+               ELSEIF ( ssh(ji,jj,Nnn) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
                   sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) &
                      &                  * tanh ( 5._wp * ( ( ssh(ji,jj,Nnn) + ht_0(ji,jj) - rn_wdmin1 ) * r_rn_wdmin1 ) )
@@ -133,11 +134,16 @@
                   sbc_tsc(ji,jj,jp_tem) = 0._wp
                ENDIF
-            ELSE 
+               sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
+            END DO
+         END DO
+      ELSE                                !* standard case
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
                sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
-            ENDIF
-
-            sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
-         END DO
-      END DO
+               sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
+            END DO
+         END DO
+      ENDIF
+      !
       IF( ln_linssh ) THEN                !* linear free surface  
          DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfmxl.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfmxl.F90	(revision 10023)
@@ -126,11 +126,11 @@
       IF( .NOT.l_offline ) THEN
          IF( iom_use("mldr10_1") ) THEN
-            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
-            ELSE                  ;  CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
+            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldr10_1", hmlp - ht_isf)   ! mixed layer thickness
+            ELSE                  ;  CALL iom_put( "mldr10_1", hmlp )           ! mixed layer depth
             END IF
          END IF
          IF( iom_use("mldkz5") ) THEN
-            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
-            ELSE                  ;  CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
+            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldkz5"  , hmld - ht_isf )   ! turbocline thickness
+            ELSE                  ;  CALL iom_put( "mldkz5"  , hmld )            ! turbocline depth
             END IF
          ENDIF
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdftke.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdftke.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdftke.F90	(revision 10023)
@@ -392,5 +392,5 @@
       !                            !  TKE due to surface and internal wave breaking
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!gm BUG : in the exp  remove the depth of ssh !!!
+!!gm BUG : in the exp  remove the depth of ssh !!!   ===>>>  not sure of that
 !!gm       i.e. use gde3w in argument (pdepw)
       
@@ -518,4 +518,8 @@
  !!gm Not sure of that coding for ISF....
       ! where wmask = 0 set zmxlm == p_e3w
+      
+!!gm  pdepw(ji,jj,mikt(ji,jj)+1) = ht_n ....
+!!gm  pdepw(ji,jj,mikt(ji,jj)) = ht_isf(:,:) + ssh(:,:,Nnn)
+      
       CASE ( 0 )           ! bounded by the distance to surface and bottom
          DO jk = 2, jpkm1
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90	(revision 10023)
@@ -81,5 +81,5 @@
       !!                   ***  FUNCTION oce_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(6)
+      INTEGER :: ierr(7)
       !!----------------------------------------------------------------------
       !
@@ -102,20 +102,16 @@
          &      riceload(jpi,jpj)                                     , STAT=ierr(2) )
          !
-         
-         
-      ALLOCATE( ssh(jpi,jpj,Nt) , STAT=ierr(5) )
-
-
+      ALLOCATE( ssh(jpi,jpj,Nt) , STAT=ierr(3) )
          !
-      ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(3) )
+      ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(4) )
          !
       ALLOCATE( ssha_e(jpi,jpj),  sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
          &        ua_e(jpi,jpj),    un_e(jpi,jpj),   ub_e(jpi,jpj),   ubb_e(jpi,jpj), &
          &        va_e(jpi,jpj),    vn_e(jpi,jpj),   vb_e(jpi,jpj),   vbb_e(jpi,jpj), &
-         &        hu_e(jpi,jpj),   hur_e(jpi,jpj),   hv_e(jpi,jpj),   hvr_e(jpi,jpj), STAT=ierr(4) )
+         &        hu_e(jpi,jpj),   hur_e(jpi,jpj),   hv_e(jpi,jpj),   hvr_e(jpi,jpj), STAT=ierr(5) )
          !
       ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_bf(jpi,jpj), vn_bf(jpi,jpj)      , STAT=ierr(6) )
 #if defined key_agrif
-      ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj)                                  , STAT=ierr(6) )
+      ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj)                                  , STAT=ierr(7) )
 #endif
          !
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 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/par_oce.F90	(revision 10023)
@@ -64,4 +64,6 @@
    INTEGER, PUBLIC, PARAMETER ::   nn_hls = 1   !: halo width (applies to both rows and columns)
 
+!!gm  thsi should be move in dom_oce 
+
    !!----------------------------------------------------------------------
    !!                   namcfg namelist parameters
@@ -82,4 +84,6 @@
    INTEGER       ::   nn_cfg           !: resolution of the configuration 
 
+!!gm  end
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step.F90	(revision 10022)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step.F90	(revision 10023)
@@ -94,4 +94,5 @@
       IF( ln_timing )   CALL timing_start('stp')
       !
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! update I/O and calendar 
@@ -273,4 +274,5 @@
                          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
