Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 9353)
@@ -96,4 +96,8 @@
          <field id="tpt_dep"      long_name="T-point depth"                  standard_name="depth_below_geoid"   unit="m"   grid_ref="grid_T_3D" />
          <field id="e3tdef"       long_name="T-cell thickness deformation"                                       unit="%"   grid_ref="grid_T_3D" />
+         <field id="stiff_tilde"  long_name="Stiffness"                                                          unit=" "   grid_ref="grid_T_3D" />
+         <field id="depw_tilde"   long_name="Interface displacement"                                             unit="m"   grid_ref="grid_W_3D" />
+         <field id="dz_tilde"     long_name="Vertical grid deformation"                                          unit=" "   grid_ref="grid_W_3D" />
+         <field id="dh_tilde"     long_name="Horizontal grid deformation"                                        unit=" "   grid_ref="grid_W_3D" />
       </field_group>
 
@@ -403,4 +407,6 @@
          <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_U_3D" />
          <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_U_3D" />
+         <field id="un_lf_tilde"  long_name="Sea Water X Velocity (low pass)"                        standard_name="low_pass_sea_water_x_velocity"     unit="m/s"        grid_ref="grid_U_3D" />
+
 
          <!-- variables available with MLE -->
@@ -445,4 +451,5 @@
          <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_V_3D" />
          <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_V_3D" />
+         <field id="vn_lf_tilde"  long_name="Sea Water Y Velocity (low pass)"                  standard_name="low_pass_sea_water_y_velocity"     unit="m/s"        grid_ref="grid_U_3D" />
 
          <!-- variables available with MLE -->
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9353)
@@ -820,8 +820,15 @@
    ln_vvl_layer  = .false.          !  full layer vertical coordinate
    ln_vvl_ztilde_as_zstar = .false. !  ztilde vertical coordinate emulating zstar
-   ln_vvl_zstar_at_eqtor = .false.  !  ztilde near the equator
-   rn_ahe3       = 0.0e0            !  thickness diffusion coefficient
+   ln_vvl_zstar_at_eqtor  = .false. !  revert to zstart at equator
+   ln_vvl_zstar_on_shelf  = .false. !  revert to zstar on shelves
+   ln_vvl_adv_cn2 = .false.         !  Centred thickness advection scheme
+   ln_vvl_adv_fct = .true.          !  FCT thickness advection scheme
+   ln_vvl_lap = .false.             !  Laplacian thickness diffusion
+      rn_ahe3_lap =  0.e0
+   ln_vvl_blp = .false.             !  Bilaplacian thickness diffusion
+      rn_ahe3_blp = -2.e8
    rn_rst_e3t    = 30.e0            !  ztilde to zstar restoration timescale [days]
    rn_lf_cutoff  = 5.0e0            !  cutoff frequency for low-pass filter  [days]
+   ln_vvl_regrid = .true.           !  Vertical regriding to ensure layer separation
    rn_zdef_max   = 0.9e0            !  maximum fractional e3t deformation
    ln_vvl_dbg    = .true.           !  debug prints    (T/F)
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt	(revision 9353)
@@ -11,2 +11,5 @@
 GYRE OPA_SRC
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
+ORCA2LIM3_16 OPA_SRC LIM_SRC_3 NST_SRC
+COMODO_IW OPA_SRC
+ORCA2LIM3_LONG OPA_SRC LIM_SRC_3 NST_SRC
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 9353)
@@ -487,5 +487,5 @@
          CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
          CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
-         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' )
          ! Vertical scale factor interpolations
          ! ------------------------------------
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 9353)
@@ -354,5 +354,5 @@
             CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
             CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
-            CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' )
             ! Vertical scale factor interpolations
             ! ------------------------------------
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9353)
@@ -214,5 +214,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehur_b, ehvr_b !: Before inverse of u and v-points ocean depth (1/m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehu_b , ehv_b  !:        depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0           !: reference depth at t-       points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0 , hf_0    !: reference depth at t- and f-points (meters)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0    !: reference depth at u- and v-points (meters)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   re2u_e1u       !: scale factor coeffs at u points (e2u/e1u)
@@ -370,5 +370,5 @@
 #endif
          !
-      ALLOCATE( hu      (jpi,jpj) , hur     (jpi,jpj) , hu_0(jpi,jpj) , ht_0  (jpi,jpj) ,     &
+      ALLOCATE( hu      (jpi,jpj) , hur     (jpi,jpj) , hu_0(jpi,jpj) , ht_0  (jpi,jpj) , hf_0  (jpi,jpj),   &
          &      hv      (jpi,jpj) , hvr     (jpi,jpj) , hv_0(jpi,jpj) , ht    (jpi,jpj) ,     &
          &      re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) ,                                       &
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9353)
@@ -70,5 +70,5 @@
       !!              - 1D configuration, move Coriolis, u and v at T-point
       !!----------------------------------------------------------------------
-      INTEGER ::   jk          ! dummy loop argument
+      INTEGER ::   jj, jk      ! dummy loop argument
       INTEGER ::   iconf = 0   ! local integers
       !!----------------------------------------------------------------------
@@ -92,4 +92,5 @@
       hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points
       hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points
+      hf_0(:,:) = 0.0_wp                       ! Reference ocean depth at F-points
       DO jk = 1, jpk
          ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
@@ -97,4 +98,11 @@
          hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
+      !   
+      DO jj=1, jpjm1
+         DO jk = 1, jpk
+            hf_0(:,jj) = hf_0(:,jj) + e3f_0(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
+         END DO
+      END DO
+      CALL lbc_lnk( hf_0, 'F', 1._wp )
       !
       IF( lk_c1d )           CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9353)
@@ -9,4 +9,5 @@
    !!                                          vvl option includes z_star and z_tilde coordinates
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
+   !!                 !  2018-01  (J. Chanut) improve ztilde robustness
    !!----------------------------------------------------------------------
    !!   'key_vvl'                              variable volume
@@ -33,4 +34,5 @@
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE bdy_oce         ! ocean open boundary conditions
 
    IMPLICIT NONE
@@ -48,21 +50,31 @@
    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
+   LOGICAL                                               :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
+   LOGICAL                                               :: ln_vvl_zstar_at_eqtor = .FALSE.     ! revert to zstar at equator
+   LOGICAL                                               :: ln_vvl_zstar_on_shelf =.FALSE.      ! revert to zstar on shelves
+   LOGICAL                                               :: ln_vvl_kepe =.FALSE.                ! kinetic/potential energy transfer
+   LOGICAL                                               :: ln_vvl_adv_fct =.FALSE.             ! Centred thickness advection
+   LOGICAL                                               :: ln_vvl_adv_cn2 =.TRUE.              ! FCT thickness advection
+   LOGICAL                                               :: ln_vvl_dbg = .FALSE.                ! debug control prints
+   LOGICAL                                               :: ln_vvl_lap                          ! Laplacian thickness diffusion
+   LOGICAL                                               :: ln_vvl_blp                          ! Bilaplacian thickness diffusion 
+   LOGICAL                                               :: ln_vvl_regrid                       ! ensure layer separation 
    !                                                                                            ! conservation: not used yet
-   REAL(wp)                                              :: rn_ahe3                   ! thickness diffusion coefficient
+   REAL(wp)                                              :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian)
+   REAL(wp)                                              :: rn_ahe3_blp               ! thickness diffusion coefficient (Bilaplacian)
    REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
    REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
    REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
-   LOGICAL , PUBLIC                                      :: ln_vvl_dbg = .FALSE.      ! debug control prints
 
    !! * Module variables
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                       ! thickness diffusion transport
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn_lf,  un_lf,  vn_lf  ! 1st order filtered arrays         
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n  ! baroclinic scale factors
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t               ! restoring period for scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv               ! restoring period for low freq. divergence
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: tildemask                 ! mask tilde tendency
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: hsm, dsm                  ! 
+   INTEGER         , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: i_int_bot
 
    !! * Substitutions
@@ -85,4 +97,5 @@
          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
             &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
+            &      i_int_bot(jpi,jpj), tildemask(jpi,jpj), hsm(jpi,jpj), dsm(jpi,jpj), &
             &      STAT = dom_vvl_alloc        )
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
@@ -92,5 +105,8 @@
       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 )
+         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj),  &
+            &        hdivn_lf(jpi,jpj,jpk),  &
+            &           un_lf(jpi,jpj,jpk),  &
+            &           vn_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')
@@ -124,7 +140,7 @@
       USE phycst,  ONLY : rpi, rsmall, rad
       !! * Local declarations
-      INTEGER ::   ji,jj,jk
+      INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
-      REAL(wp)::   zcoef
+      REAL(wp)::   zcoef, zwgt, ztmp, zhmin, zhmax
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
@@ -155,5 +171,5 @@
       CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
       CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' )
       ! Vertical scale factor interpolations
       ! ------------------------------------
@@ -199,18 +215,24 @@
          hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
       END DO
-      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
-      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
+      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
+      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
 
       ! Restoring frequencies for z_tilde coordinate
       ! ============================================
+      tildemask(:,:) = 1._wp
+
       IF( ln_vvl_ztilde ) THEN
          ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
          frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
          frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
+         ! tendency mask:
+         !
          IF( ln_vvl_ztilde_as_zstar ) THEN
             ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
             frq_rst_e3t(:,:) = 0.0_wp 
             frq_rst_hdv(:,:) = 1.0_wp / rdt
+            tildemask(:,:) = 0._wp
          ENDIF
+         
          IF ( ln_vvl_zstar_at_eqtor ) THEN
             DO jj = 1, jpj
@@ -219,9 +241,11 @@
                      ! values outside the equatorial band and transition zone (ztilde)
                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
-                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
+!                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
+              
                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
                      ! values inside the equatorial band (ztilde as zstar)
                      frq_rst_e3t(ji,jj) =  0.0_wp
-                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
+!                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
+                     tildemask(ji,jj) = 0._wp
                   ELSE
                      ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
@@ -229,8 +253,10 @@
                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
                         &                                          * 180._wp / 3.5_wp ) )
-                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
-                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
-                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
-                        &                                          * 180._wp / 3.5_wp ) )
+!                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
+!                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
+!                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+!                        &                                          * 180._wp / 3.5_wp ) )
+                     tildemask(ji,jj) = 0.5_wp * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
+                        &                                                 * 180._wp / 3.5_wp ) )
                   ENDIF
                END DO
@@ -240,7 +266,31 @@
                ij0 = 128   ;   ij1 = 135   ;   
                frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
-               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
+               tildemask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )   = 0.0_wp
+!               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
             ENDIF
          ENDIF
+         !
+         IF ( ln_vvl_zstar_on_shelf ) THEN
+            zhmin =  50._wp
+            zhmax = 100._wp
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zwgt = 1._wp
+                  IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
+                     zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
+                  ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
+                     zwgt = 0._wp
+                  ENDIF
+                  frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt)
+                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt)
+               END DO
+            END DO
+         ENDIF
+         !
+         ztmp = MAXVAL( frq_rst_hdv(:,:) )
+         IF( lk_mpp )   CALL mpp_max( ztmp )                 ! max over the global domain
+         !
+         IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' )
+         !
       ENDIF
 
@@ -272,5 +322,5 @@
       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t, ztu, ztv
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
       !! * Arguments
@@ -279,12 +329,15 @@
       !! * Local declarations
       INTEGER                                :: ji, jj, jk            ! dummy loop indices
+      INTEGER                                :: ib, ib_bdy, ip, jp    !   "     "     "
       INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
       REAL(wp)                               :: z2dt                  ! temporary scalars
       REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
+      REAL(wp)                               :: zalpha, zwgt          ! temporary scalars
+      REAL(wp)                               :: zdu, zdv
       LOGICAL                                :: ll_do_bclinic         ! temporary logical
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
       CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
-      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
+      CALL wrk_alloc( jpi, jpj, jpk, ze3t, ztu, ztv )
 
       IF(kt == nit000)   THEN
@@ -296,5 +349,5 @@
       ll_do_bclinic = .TRUE.
       IF( PRESENT(kcall) ) THEN
-         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
+         IF ( kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE.
       ENDIF
 
@@ -307,102 +360,168 @@
       !                                                ! --------------------------------------------- !
 
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - ssmask(:,:) )
       DO jk = 1, jpkm1
          ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
          fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
       END DO
-
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
-         !                                                            ! ------baroclinic part------ !
-
-         ! I - initialization
-         ! ==================
-
-         ! 1 - barotropic divergence
-         ! -------------------------
+      
+      IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
+         
+         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
+         un_td(:,:,:) = 0.0_wp        ! Transport corrections
+         vn_td(:,:,:) = 0.0_wp
+
          zhdiv(:,:) = 0.
-         zht(:,:)   = 0.
          DO jk = 1, jpkm1
             zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
-
-         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
-         ! --------------------------------------------------
-         IF( ln_vvl_ztilde ) THEN
-            IF( kt .GT. nit000 ) THEN
+         END DO
+         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )      
+
+         ! Thickness advection:
+         ! --------------------
+         ! Set advection velocities and source term
+         IF ( ln_vvl_ztilde ) THEN
+            !
+            IF ((kt==nit000).AND.(neuler==0)) THEN
                DO jk = 1, jpkm1
-                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+                  ztu(:,:,jk) = un(:,:,jk)
+                  ztv(:,:,jk) = vn(:,:,jk)
+               END DO
+            ELSE
+               DO jk = 1, jpkm1
+                  ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk)
+                  ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk)
+                  tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk)
                END DO
             ENDIF
-         END IF
-
-         ! II - after z_tilde increments of vertical scale factors
-         ! =======================================================
-         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
-
-         ! 1 - High frequency divergence term
-         ! ----------------------------------
-         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
+            !
+         ELSEIF ( ln_vvl_layer ) THEN
+            !
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
-            END DO
-         ELSE                         ! layer case
+               ztu(:,:,jk) = un(:,:,jk)
+               ztv(:,:,jk) = vn(:,:,jk)
+            END DO
+            !
+         ENDIF
+         !
+         ! Do advection
+         IF     (ln_vvl_adv_fct) THEN
+            CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv )
+            !
+         ELSEIF (ln_vvl_adv_cn2) THEN
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
-            END DO
-         END IF
-
-         ! 2 - Restoring term (z-tilde case only)
-         ! ------------------
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * hdivn(:,:,jk)
+            END DO
+         ENDIF
+         ! 
+         ! Thickness anomlaly diffusion:
+         ! -----------------------------
+         zwu(:,:) = 0.0_wp
+         zwv(:,:) = 0.0_wp
+         ztu(:,:,:) = 0.0_wp
+         ztv(:,:,:) = 0.0_wp
+
+         IF ( ln_vvl_blp ) THEN  ! Bilaplacian
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1                 ! First derivative (gradient)
+                  DO ji = 1, fs_jpim1   ! vector opt.                  
+                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+                  END DO
+               END DO
+
+               DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
+                  DO ji = fs_2, fs_jpim1 ! vector opt.
+                     zht(ji,jj) =  rn_ahe3_blp * r1_e12t(ji,jj) * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
+                        &                                           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
+                  END DO
+               END DO
+
+#if defined key_bdy
+               DO ib_bdy=1, nb_bdy
+                  DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
+                     ji = idx_bdy(ib_bdy)%nbi(ib,1)
+                     jj = idx_bdy(ib_bdy)%nbj(ib,1)
+                     zht(ji,jj) = 0._wp
+                  END DO 
+               END DO
+#endif
+               CALL lbc_lnk( zht, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
+
+               DO jj = 1, jpjm1                 ! third derivative (gradient)
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) )
+                     ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( zht(ji  ,jj+1) - zht(ji,jj) )
+                     zwu(ji,jj) = zwu(ji,jj) + ztu(ji,jj,jk)
+                     zwv(ji,jj) = zwv(ji,jj) + ztv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+
+         IF ( ln_vvl_lap ) THEN    ! Laplacian
+            DO jk = 1, jpkm1                    ! First derivative (gradient)
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.                  
+                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+                     zwu(ji,jj) = zwu(ji,jj) + zdu
+                     zwv(ji,jj) = zwv(ji,jj) + zdv
+                     ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu
+                     ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv
+                  END DO
+               END DO
+            END DO
+         ENDIF 
+
+         ! Ensure barotropic fluxes are null:
+!         DO jj = 1, jpj
+!            DO ji = 1, jpi
+!               DO jk = 1, jpkm1
+!                  ztu(ji,jj,jk) = ztu(ji,jj,jk) - zwu(ji,jj)*fse3u_b(ji,jj,jk)*hur_b(ji,jj)*umask(ji,jj,jk)
+!                  ztv(ji,jj,jk) = ztv(ji,jj,jk) - zwv(ji,jj)*fse3v_b(ji,jj,jk)*hvr_b(ji,jj)*vmask(ji,jj,jk)
+!               END DO
+!            END DO
+!         END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ztu(ji,jj,mbku(ji,jj)) = ztu(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
+               ztv(ji,jj,mbkv(ji,jj)) = ztv(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
+            END DO
+         END DO
+
+         ! divergence of diffusive fluxes
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk) - ztu(ji,jj,jk)    &
+                     &                                          +     ztv(ji  ,jj-1,jk) - ztv(ji,jj,jk)    &
+                     &                                            ) * r1_e12t(ji,jj)
+               END DO
+            END DO
+         END DO
+ 
+         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:)
+         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:)
+         CALL lbc_lnk( un_td , 'U' , -1.)
+         CALL lbc_lnk( vn_td , 'V' , -1.)  
+         !
+         !
+         ! Restoring:
          IF( ln_vvl_ztilde ) THEN
             DO jk = 1, jpk
                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
             END DO
-         END IF
-
-         ! 3 - Thickness diffusion term
-         ! ----------------------------
-         zwu(:,:) = 0.0_wp
-         zwv(:,:) = 0.0_wp
-         ! a - first derivative: diffusive fluxes
+         ENDIF
+
+         ! Remove "external thickness" tendency:
          DO jk = 1, jpkm1
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
-                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
-                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
-                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
-                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
-                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         ! b - correction for last oceanic u-v points
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
-               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
-            END DO
-         END DO
-         ! c - second derivative: divergence of diffusive fluxes
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
-                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
-                     &                                            ) * r1_e12t(ji,jj)
-               END DO
-            END DO
-         END DO
-         ! d - thickness diffusion transport: boundary conditions
-         !     (stored for tracer advction and continuity equation)
-         CALL lbc_lnk( un_td , 'U' , -1._wp)
-         CALL lbc_lnk( vn_td , 'V' , -1._wp)
-
-         ! 4 - Time stepping of baroclinic scale factors
-         ! ---------------------------------------------
+            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   fse3t_n(:,:,jk) * zhdiv(:,:)
+         END DO 
+                   
          ! Leapfrog time stepping
          ! ~~~~~~~~~~~~~~~~~~~~~~
@@ -412,6 +531,37 @@
             z2dt = 2.0_wp * rdt
          ENDIF
-         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
+
          tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
+
+         ! Revert to zstar locally:
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~
+         DO jk=1,jpkm1
+            tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) * tildemask(:,:) 
+         END DO
+
+         ! Ensure layer separation:
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~
+         IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt )
+  
+         ! Boundary conditions:
+         ! ~~~~~~~~~~~~~~~~~~~~
+#if defined key_bdy
+         DO ib_bdy=1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
+!!            DO ib = 1, idx_bdy(ib_bdy)%nblen(1)
+               ji = idx_bdy(ib_bdy)%nbi(ib,1)
+               jj = idx_bdy(ib_bdy)%nbj(ib,1)
+               zwgt = idx_bdy(ib_bdy)%nbw(ib,1)
+               ip = bdytmask(ji+1,jj  ) - bdytmask(ji-1,jj  )
+               jp = bdytmask(ji  ,jj+1) - bdytmask(ji  ,jj-1)
+               DO jk = 1, jpkm1  
+                  tilde_e3t_a(ji,jj,jk) = 0.e0
+!!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt)
+!!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk)
+               END DO
+            END DO 
+         END DO
+#endif
+         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )        
 
          ! Maximum deformation control
@@ -446,48 +596,85 @@
             ENDIF
          ENDIF
-         ! - ML - end test
-         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
-         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
-         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
-
-         !
-         ! "tilda" change in the after scale factor
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         IF ( ln_vvl_ztilde ) THEN
+            zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
+            DO jk = 1, jpkm1
+               ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:) + un_td(:,:,jk)
+               ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk)
+               ze3t(:,:,jk) = -fse3t_n(:,:,jk) * zhdiv(:,:)
+!
+!               DO jj = 2, jpjm1
+!                  DO ji = fs_2, fs_jpim1   ! vector opt.
+!                     
+!                     ze3t(ji,jj,jk) =   -fse3t_n(ji,jj,jk) * zhdiv(ji,jj)   &
+!                        & + (  un_td(ji,jj,jk) - un_td(ji-1,jj  ,jk)        &
+!                        &    + vn_td(ji,jj,jk) - vn_td(ji  ,jj-1,jk)  )     &
+!                        & / ( e1t(ji,jj) * e2t(ji,jj) )
+!                  END DO
+!               END DO 
+            END DO
+            !
+            un_lf(:,:,:)    = un_lf(:,:,:)    * (1._wp - zalpha) + zalpha * ztu(:,:,:)
+            vn_lf(:,:,:)    = vn_lf(:,:,:)    * (1._wp - zalpha) + zalpha * ztv(:,:,:)
+            hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
+         ENDIF
+
+      ENDIF
+
+      IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.(.NOT.ll_do_bclinic) ) THEN
+         zhdiv(:,:) = 0.
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
-         END DO
-         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
-         ! ==================================================================================
-         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
-         ! - ML - baroclinicity error should be better treated in the future
-         !        i.e. locally and not spread over the water column.
-         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
-         zht(:,:) = 0.
+            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * (hdivn(:,:,jk) - hdivb(:,:,jk))
+         END DO
+         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+
+         IF( neuler == 0 .AND. kt == nit000 ) THEN
+            z2dt =  rdt
+         ELSE
+            z2dt = 2.0_wp * rdt
+         ENDIF
+
          DO jk = 1, jpkm1
-            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * & 
+                                  & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:))
+         END DO
+      ENDIF
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
+      !                                             ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-         END DO
-
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                     
+         END DO
       ENDIF
 
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
-      !                                           ! ---baroclinic part--------- !
+      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
+         !
+         zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-      ENDIF
-
-      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
-         !
+            zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
          IF( lwp ) WRITE(numout, *) 'kt =', kt
          IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
+            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 )
             IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
             IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
          END IF
          !
+         z_tmin = MINVAL( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0  )
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin         
+         !
+         z_tmin = MINVAL( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 )
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3u_n) =', z_tmin
+         !
+         z_tmin = MINVAL( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 )
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3v_n) =', z_tmin
+         !
+         z_tmin = MINVAL( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 )
+         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3f_n) =', z_tmin
+         !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
@@ -496,33 +683,58 @@
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn  -SUM(fse3t_n))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+         END DO
+         zwu(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e12u(ji,jj)                             &
+                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji+1,jj) * sshn(ji+1,jj) )
+            END DO
+         END DO
+         CALL lbc_lnk( zwu(:,:), 'U', 1._wp )
+         z_tmax = MAXVAL( umask(:,:,1) * umask_i(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(fse3u_n))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+         END DO
+         zwv(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e12v(ji,jj)                             &
+                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji,jj+1) * sshn(ji,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( zwv(:,:), 'V', 1._wp )
+         z_tmax = MAXVAL( vmask(:,:,1) * vmask_i(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(fse3v_n))) =', z_tmax
+         !
+         zht(:,:) = 0.0_wp
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
+            END DO
+         END DO
+         zwu(:,:) = 0._wp
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e12f(ji,jj)  &
+                        &          * (  e12t(ji  ,jj) * sshn(ji  ,jj) + e12t(ji  ,jj+1) * sshn(ji  ,jj+1) &
+                        &             + e12t(ji+1,jj) * sshn(ji+1,jj) + e12t(ji+1,jj+1) * sshn(ji+1,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( zht(:,:), 'F', 1._wp )
+         CALL lbc_lnk( zwu(:,:), 'F', 1._wp )
+         z_tmax = MAXVAL( fmask(:,:,1) * fmask_i(:,:) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(fse3f_n))) =', z_tmax
+         !
       END IF
 
@@ -549,5 +761,5 @@
 
       CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
-      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze3t, ztu, ztv           )
 
       IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
@@ -583,6 +795,6 @@
       INTEGER, INTENT( in )               :: kt       ! time step
       !! * Local declarations
-      INTEGER                             :: ji,jj,jk       ! dummy loop indices
-      REAL(wp)                            :: zcoef
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwrk
+      INTEGER                             :: jk       ! dummy loop indices
       !!----------------------------------------------------------------------
 
@@ -607,4 +819,5 @@
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
+
       fsdept_b(:,:,:) = fsdept_n(:,:,:)
       fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
@@ -620,5 +833,5 @@
       ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
-      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n (:,:,:), 'F'  )
       ! Vertical scale factor interpolations
       ! ------------------------------------
@@ -631,23 +844,12 @@
       ! t- and w- points depth
       ! ----------------------
-      ! set the isf depth as it is in the initial step
       fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
       fsdepw_n(:,:,1) = 0.0_wp
       fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-
       DO jk = 2, jpk
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-                                                                 ! 1 for jk = mikt
-               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
-               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
-                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
-               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
-            END DO
-         END DO
+         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
+         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
+         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
       END DO
-
       ! Local depth and Inverse of the local depth of the water column at u- and v- points
       ! ----------------------------------------------------------------------------------
@@ -666,4 +868,8 @@
       END DO
 
+      ! Write additional diagnostics 
+      ! ============================
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) CALL dom_vvl_dia( kt) 
+
       ! write restart file
       ! ==================
@@ -674,4 +880,168 @@
    END SUBROUTINE dom_vvl_sf_swp
 
+   SUBROUTINE dom_vvl_dia( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_dia  ***
+      !!                   
+      !! ** Purpose :  Output some diagnostics in ztilde/zlayer cases
+      !!
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in )               :: kt       ! time step
+      !! * Local declarations
+      INTEGER                             :: ji,jj,jk ! dummy loop indices
+      REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, ztmp1, z2dt
+      REAL(wp), DIMENSION(4)              :: zr1
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zout
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_dia')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zwdw, zout )
+      !
+      ! Compute internal interfaces depths:
+      !------------------------------------
+      IF ( iom_use("dh_tilde").OR.iom_use("depw_tilde").OR.iom_use("stiff_tilde")) THEN
+         zwdw(:,:,1) = 0.e0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 2, jpkm1
+                  zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + & 
+                                 & (tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      ! Output interface depth anomaly:
+      ! -------------------------------
+      IF ( iom_use("depw_tilde") ) CALL iom_put( "depw_tilde", (zwdw(:,:,:)-gdepw_0(:,:,:))*tmask(:,:,:) )
+      !
+      ! Output grid stiffness (T-points):
+      ! ---------------------------------
+      IF ( iom_use("stiff_tilde"  ) ) THEN
+         zr1(:)   = 0.e0
+         zout(:,:,jpk) = 0.e0      
+         DO ji = 2, jpim1
+            DO jj = 2, jpjm1
+               DO jk = 1, jpkm1
+                  zr1(1) = umask(ji-1,jj  ,jk) *abs( (zwdw(ji  ,jj  ,jk  )-zwdw(ji-1,jj  ,jk  )  & 
+                       &                             +zwdw(ji  ,jj  ,jk+1)-zwdw(ji-1,jj  ,jk+1)) &
+                       &                            /(zwdw(ji  ,jj  ,jk  )+zwdw(ji-1,jj  ,jk  )  &
+                       &                             -zwdw(ji  ,jj  ,jk+1)-zwdw(ji-1,jj  ,jk+1) + rsmall) )
+                  zr1(2) = umask(ji  ,jj  ,jk) *abs( (zwdw(ji+1,jj  ,jk  )-zwdw(ji  ,jj  ,jk  )  &
+                       &                             +zwdw(ji+1,jj  ,jk+1)-zwdw(ji  ,jj  ,jk+1)) &
+                       &                            /(zwdw(ji+1,jj  ,jk  )+zwdw(ji  ,jj  ,jk  )  &
+                       &                             -zwdw(ji+1,jj  ,jk+1)-zwdw(ji  ,jj  ,jk+1) + rsmall) )
+                  zr1(3) = vmask(ji  ,jj  ,jk) *abs( (zwdw(ji  ,jj+1,jk  )-zwdw(ji  ,jj  ,jk  )  &
+                       &                             +zwdw(ji  ,jj+1,jk+1)-zwdw(ji  ,jj  ,jk+1)) &
+                       &                            /(zwdw(ji  ,jj+1,jk  )+zwdw(ji  ,jj  ,jk  )  &
+                       &                             -zwdw(ji  ,jj+1,jk+1)-zwdw(ji  ,jj  ,jk+1) + rsmall) )
+                  zr1(4) = vmask(ji  ,jj-1,jk) *abs( (zwdw(ji  ,jj  ,jk  )-zwdw(ji  ,jj-1,jk  )  &
+                       &                             +zwdw(ji  ,jj  ,jk+1)-zwdw(ji  ,jj-1,jk+1)) &
+                       &                            /(zwdw(ji  ,jj  ,jk  )+zwdw(ji  ,jj-1,jk  )  &
+                       &                             -zwdw(ji,  jj  ,jk+1)-zwdw(ji  ,jj-1,jk+1) + rsmall) )
+                  zout(ji,jj,jk) = MAXVAL(zr1(1:4))
+               END DO
+            END DO
+         END DO
+
+         CALL lbc_lnk( zout, 'T', 1. )
+         CALL iom_put( "stiff_tilde", zout(:,:,:) ) 
+      END IF
+      ! Output Horizontal Laplacian of interfaces depths (W-points):
+      ! ------------------------------------------------------------
+      IF ( iom_use("dh_tilde")   ) THEN
+         !
+         zout(:,:,1  )=0._wp
+         zout(:,:,jpk)=0._wp
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) )
+                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) )
+                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) )
+                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) )
+!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj))  
+                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 
+!                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)
+                  zout(ji,jj,jk) = ztmp1
+               END DO
+            END DO            
+         END DO
+         ! Mask open boundaries:
+#if defined key_bdy
+         IF (lk_bdy) THEN
+            DO jk = 1, jpkm1
+               zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:)
+            END DO
+         ENDIF
+#endif
+         CALL lbc_lnk( zout(:,:,:), 'T', 1. )
+         zwdw(:,:,:) = zout(:,:,:)
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) )
+                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) )
+                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) )
+                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) )
+!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj))  
+                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 
+                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)
+               END DO
+            END DO            
+         END DO
+         ! Mask open boundaries:
+#if defined key_bdy
+         IF (lk_bdy) THEN
+            DO jk = 1, jpkm1
+               zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:)
+            END DO
+         ENDIF
+#endif
+         CALL lbc_lnk( zout(:,:,:), 'T', 1. )
+         !
+         CALL iom_put( "dh_tilde", zout(:,:,:) )
+      ENDIF
+      !
+      ! Output vertical Laplacian of interfaces depths (W-points):
+      ! ----------------------------------------------------------
+      IF ( iom_use("dz_tilde"  ) ) THEN
+         zout(:,:,1  ) = 0._wp
+         zout(:,:,jpk) = 0._wp
+         DO jk=2,jpkm1
+            zout(:,:,jk) = 2._wp*ABS(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)-tilde_e3t_n(:,:,jk-1)-e3t_0(:,:,jk-1)) & 
+                               &   /(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)+tilde_e3t_n(:,:,jk-1)+e3t_0(:,:,jk-1)) &
+                               &   * tmask(:,:,jk)
+         END DO
+         CALL iom_put( "dz_tilde", zout(:,:,:) ) 
+      END IF
+      !
+      !
+      ! Output low pass U-velocity:
+      ! ---------------------------
+      IF ( iom_use("un_lf_tilde"  ).AND.ln_vvl_ztilde ) THEN
+         zout(:,:,jpk) = 0.e0  
+         DO jk=1,jpkm1
+            zout(:,:,jk) = un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:)
+         END DO
+         CALL iom_put( "un_lf_tilde", zout(:,:,:) )
+      END IF
+      !
+      ! Output low pass V-velocity:
+      ! ---------------------------
+      IF ( iom_use("vn_lf_tilde"  ).AND.ln_vvl_ztilde ) THEN
+         zout(:,:,jpk) = 0.e0  
+         DO jk=1,jpkm1
+            zout(:,:,jk) = vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:)
+         END DO
+         CALL iom_put( "vn_lf_tilde", zout(:,:,:) )
+      END IF
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zwdw, zout )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_dia')
+      !
+   END SUBROUTINE dom_vvl_dia
 
    SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
@@ -691,6 +1061,11 @@
       !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
       !! * Local declarations
+      INTEGER ::   nmet                                                ! horizontal interpolation method
       INTEGER ::   ji, jj, jk                                          ! dummy loop indices
+      INTEGER ::   jkbot                                               ! bottom level
       LOGICAL ::   l_is_orca                                           ! local logical
+      REAL(wp) :: zmin, zdo, zup, ztap
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zs                        ! Surface interface depth anomaly
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw                        ! Interface depth anomaly
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
@@ -698,19 +1073,118 @@
       l_is_orca = .FALSE.
       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
+
+ 
+      nmet=1 ! Original method (Surely wrong)
+!      nmet=1 ! Interface interpolation
+!      nmet=2 ! Internal interfaces interpolation only, spread barotropic increment
+!      nmet=1
+      ztap = 0.1_wp ! Minimum fraction of T-point thickness at cell interfaces
+
+      IF ( (nmet==1).OR.(nmet==2) ) THEN
+         SELECT CASE ( pout )
+            !
+         CASE( 'U', 'V', 'F' )
+            ! Compute interface depth anomaly at T-points
+            CALL wrk_alloc( jpi, jpj, jpk, zw )
+            CALL wrk_alloc( jpi, jpj, zs )
+            !
+            zw(:,:,:) =  0._wp    
+            !
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  jkbot = mbkt(ji,jj)
+                  DO jk=jkbot,1,-1
+                     zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk)
+                  END DO
+               END DO
+            END DO 
+            !
+            !
+            IF (nmet==2) THEN        ! Consider "internal" interfaces only
+               zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh)
+               !
+               zw(:,:,:) = 0._wp
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     jkbot = mbkt(ji,jj)
+                     DO jk=jkbot,1,-1
+                        zw(ji,jj,jk) = zw(ji,jj,jk+1) - ( pe3_in(ji,jj,jk)                     &
+                                     & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1))  &
+                                     & - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO 
+            ENDIF 
+            !
+         END SELECT
+      END IF
+
+      pe3_out(:,:,:) = 0._wp
 
       SELECT CASE ( pout )
          !               ! ------------------------------------- !
       CASE( 'U' )        ! interpolation from T-point to U-point !
-         !               ! ------------------------------------- !
+         !               ! ------------------------------------- !     
          ! horizontal surface weighted interpolation
-         DO jk = 1, jpk
+         IF (nmet==0) THEN
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
+                        &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         !
+         IF ( (nmet==1).OR.(nmet==2) ) THEN
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
-                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
+                  ! Correction at last level:
+                  jkbot = mbku(ji,jj)
+                  pe3_out(ji,jj,jkbot) = -0.5_wp * umask(ji,jj,jkbot) * r1_e12u(ji,jj) &
+                        &                    * (   e12t(ji  ,jj) * zw(ji  ,jj,jkbot)   &
+                        &                        + e12t(ji+1,jj) * zw(ji+1,jj,jkbot)   )
+                  !
+                  ! If there is a step, taper bottom interface:
+                  IF (hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ) THEN
+                     IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
+!                        zmin = ztap * pe3_in(ji+1,jj,jkbot)
+                        zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
+                     ELSE
+!                        zmin = ztap * pe3_in(ji  ,jj,jkbot)
+                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
+                     ENDIF
+                     zmin = -e3u_0(ji,jj,jkbot) + zmin
+                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
+                  ENDIF
+                  !
+                  zdo =  -pe3_out(ji,jj,jkbot)
+                  DO jk=jkbot-1,1,-1
+                     zup = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji  ,jj) &
+                              &                     *(   e12t(ji  ,jj) * zw(ji  ,jj,jk)  &
+                              &                         +e12t(ji+1,jj) * zw(ji+1,jj,jk)  )
+                     pe3_out(ji,jj,jk) = zdo - zup  
+                     zdo =  zdo - pe3_out(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF
+         ! 
+         IF (nmet==2) THEN           ! Spread sea level anomaly
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  DO jk=1,jpk
+                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                 &
+                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )             & 
+                           &               / ( hu_0(ji,jj) + 1._wp - umask_i(ji,jj) )            &
+                           &               * 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)           &
+                           &               * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji+1,jj)*zs(ji+1,jj) )       
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF
          !
          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
@@ -718,17 +1192,74 @@
          CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
+         ! 
+         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
+         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
+         !
          !               ! ------------------------------------- !
       CASE( 'V' )        ! interpolation from T-point to V-point !
          !               ! ------------------------------------- !
          ! horizontal surface weighted interpolation
-         DO jk = 1, jpk
+         IF (nmet==0) THEN
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
+                        &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         !
+         IF ( (nmet==1).OR.(nmet==2) ) THEN
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
-                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
+                  ! Correction at last level:
+                  jkbot = mbkv(ji,jj)
+                  pe3_out(ji,jj,jkbot) = -0.5_wp * vmask(ji,jj,jkbot) * r1_e12v(ji,jj) &
+                        &                    * (   e12t(ji,jj  ) * zw(ji,jj  ,jkbot)   &
+                        &                        + e12t(ji,jj+1) * zw(ji,jj+1,jkbot)   )
+                  !
+                  ! If there is a step, taper bottom interface:
+                  IF (hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ) THEN
+                     IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
+!                        zmin = ztap * pe3_in(ji,jj+1,jkbot)
+                        zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
+                     ELSE
+!                        zmin = ztap * pe3_in(ji,jj  ,jkbot)
+                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
+                     ENDIF
+                     zmin = -e3v_0(ji,jj,jkbot) + zmin
+                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
+                  ENDIF
+                  !
+                  zdo =  -pe3_out(ji,jj,jkbot)
+                  DO jk=jkbot-1,1,-1
+                     zup = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj  ) &
+                              &                     *  ( e12t(ji,jj  ) * zw(ji,jj  ,jk) &
+                              &                         +e12t(ji,jj+1) * zw(ji,jj+1,jk) )
+                     !
+                     pe3_out(ji,jj,jk) = zdo - zup 
+                     zdo =  zdo - pe3_out(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF
+         ! 
+         IF (nmet==2) THEN           ! Spread sea level anomaly
+            !
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  DO jk=1,jpk
+                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                 &
+                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )             & 
+                           &               / ( hv_0(ji,jj) + 1._wp - vmask_i(ji,jj) )            &
+                           &               * 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)           &
+                           &               * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji,jj+1)*zs(ji,jj+1) )       
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF
          !
          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
@@ -736,17 +1267,91 @@
          CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
+         ! 
+         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
+         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
+         !
          !               ! ------------------------------------- !
-      CASE( 'F' )        ! interpolation from U-point to F-point !
+      CASE( 'F' )        ! interpolation from T-point to F-point !
          !               ! ------------------------------------- !
          ! horizontal surface weighted interpolation
-         DO jk = 1, jpk
+         IF (nmet==0) THEN
+            DO jk=1,jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
+                           &      * (  e12t(ji  ,jj  ) * ( pe3_in(ji  ,jj  ,jk)-e3t_0(ji  ,jj  ,jk) )   & 
+                           &         + e12t(ji  ,jj+1) * ( pe3_in(ji  ,jj+1,jk)-e3t_0(ji  ,jj+1,jk) )   &
+                           &         + e12t(ji+1,jj  ) * ( pe3_in(ji+1,jj  ,jk)-e3t_0(ji+1,jj  ,jk) )   & 
+                           &         + e12t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ))                                                      
+                  END DO
+               END DO
+            END DO 
+         ENDIF
+         !
+         IF ( (nmet==1).OR.(nmet==2) ) THEN
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
-                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
+                  ! bottom correction:
+                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 
+                  pe3_out(ji,jj,jkbot) = -0.25_wp * umask(ji,jj,jkbot) * umask(ji,jj+1,jkbot) * r1_e12f(ji,jj) &
+                        &                         * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jkbot)   &
+                        &                            + e12t(ji+1,jj  ) * zw(ji+1,jj  ,jkbot)   &
+                        &                            + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jkbot)   &
+                        &                            + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jkbot)   )
+                  !
+                  ! If there is a step, taper bottom interface:
+                  IF (hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ) THEN
+                     IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
+                        IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN
+!                           zmin = ztap * pe3_in(ji+1,jj+1,jkbot)
+                           zmin = ztap * (-zw(ji+1,jj+1,jkbot)+e3t_0(ji+1,jj+1,jkbot))
+                        ELSE
+!                           zmin = ztap * pe3_in(ji  ,jj+1,jkbot)
+                           zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
+                        ENDIF
+                     ELSE
+                        IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
+!                           zmin = ztap * pe3_in(ji+1,jj  ,jkbot)
+                           zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
+                        ELSE
+!                           zmin = ztap * pe3_in(ji  ,jj  ,jkbot)
+                           zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
+                        ENDIF
+                     ENDIF
+                     zmin = -e3f_0(ji,jj,jkbot) + zmin
+                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
+                  ENDIF
+                  !
+                  zdo =  -pe3_out(ji,jj,jkbot)
+                  DO jk=jkbot-1,1,-1
+                     zup =  0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
+                           &        * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  & 
+                           &           + e12t(ji+1,jj  ) * zw(ji+1,jj  ,jk)  & 
+                           &           + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  & 
+                           &           + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  )
+                     pe3_out(ji,jj,jk) = zdo - zup 
+                     !
+                     zdo =  zdo - pe3_out(ji,jj,jk)                                    
+                  END DO
+                  !
+               END DO
+            END DO 
+         ENDIF
+         !
+         IF (nmet==2) THEN           ! Spread sea level anomaly
+            !
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  DO jk=1,jpk
+                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                         &
+                           &               + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) )                     & 
+                           &               / ( hf_0(ji,jj) + 1._wp - umask_i(ji,jj)*umask_i(ji,jj+1) )   &
+                           &               * 0.25_wp * umask(ji,jj,jk)*umask(ji,jj+1,jk)*r1_e12f(ji,jj)  &
+                           &               * ( e12t(ji  ,jj)*zs(ji  ,jj) + e12t(ji  ,jj+1)*zs(ji  ,jj+1) &
+                           &                  +e12t(ji+1,jj)*zs(ji+1,jj) + e12t(ji+1,jj+1)*zs(ji+1,jj+1) )               
+                  END DO
+               END DO
+            END DO
+         END IF
          !
          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
@@ -754,4 +1359,8 @@
          CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
+         ! 
+         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
+         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
+         !
          !               ! ------------------------------------- !
       CASE( 'W' )        ! interpolation from T-point to W-point !
@@ -808,5 +1417,6 @@
       !! * Local declarations
       INTEGER ::   jk
-      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
+      INTEGER, DIMENSION(7) ::   id            ! local integers
+      REAL(wp), POINTER, DIMENSION(:,:) :: zhdiv
       !!----------------------------------------------------------------------
       !
@@ -818,13 +1428,16 @@
             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
             !
-            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
-            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
-            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
-            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
+            id(1) = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
+            id(2) = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
+            id(3) = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
+            id(4) = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
+            id(5) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
+            id(6) = iom_varid( numror, 'un_lf'   , ldstop = .FALSE. )
+            id(7) = iom_varid( numror, 'vn_lf'   , ldstop = .FALSE. )
+
             !                             ! --------- !
             !                             ! all cases !
             !                             ! --------- !
-            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
+            IF( MIN( id(1), id(2) ) > 0 ) THEN       ! all required arrays exist
                CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
                CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
@@ -838,12 +1451,12 @@
                   fse3t_b(:,:,:) = fse3t_n(:,:,:)
                ENDIF
-            ELSE IF( id1 > 0 ) THEN
+            ELSE IF( id(1) > 0 ) THEN
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
                IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
                CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
-               fse3t_n(:,:,:) = fse3t_b(:,:,:)
+               fse3t_b(:,:,:) = fse3t_n(:,:,:)
                neuler = 0
-            ELSE IF( id2 > 0 ) THEN
+            ELSE IF( id(2) > 0 ) THEN
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
                IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
@@ -867,5 +1480,5 @@
             IF( ln_vvl_zstar ) THEN       ! z_star case !
                !                          ! ----------- !
-               IF( MIN( id3, id4 ) > 0 ) THEN
+               IF( MIN( id(3), id(4) ) > 0 ) THEN
                   CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
                ENDIF
@@ -873,5 +1486,5 @@
             ELSE                          ! z_tilde and layer cases !
                !                          ! ----------------------- !
-               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
+               IF( MIN( id(3), id(4) ) > 0 ) THEN  ! all required arrays exist
                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
@@ -879,14 +1492,21 @@
                   tilde_e3t_b(:,:,:) = 0.0_wp
                   tilde_e3t_n(:,:,:) = 0.0_wp
-               ENDIF
-               !                          ! ------------ !
+                  !
+                  neuler = 0
+               ENDIF                      ! ------------ !
                IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
                   !                       ! ------------ !
-                  IF( id5 > 0 ) THEN  ! required array exists
-                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
-                  ELSE                ! array is missing
-                     hdiv_lf(:,:,:) = 0.0_wp
+                  IF( MINVAL(id(5:7) ) > 0 ) THEN  ! all required arrays exist
+                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) )
+                     CALL iom_get( numror, jpdom_autoglo, 'un_lf'   ,    un_lf(:,:,:) )
+                     CALL iom_get( numror, jpdom_autoglo, 'vn_lf'   ,    vn_lf(:,:,:) )
+                  ELSE                            ! one at least array is missing
+                     hdivn_lf(:,:,:) = 0.0_wp
+                     un_lf(:,:,:)    = 0.0_wp
+                     vn_lf(:,:,:)    = 0.0_wp
+                     neuler = 0
                   ENDIF
                ENDIF
+               !
             ENDIF
             !
@@ -898,6 +1518,10 @@
                tilde_e3t_b(:,:,:) = 0.0_wp
                tilde_e3t_n(:,:,:) = 0.0_wp
-               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
             END IF
+            IF( ln_vvl_ztilde ) THEN
+               hdivn_lf(:,:,:) = 0.0_wp
+               un_lf(:,:,:)    = 0.0_wp
+               vn_lf(:,:,:)    = 0.0_wp
+            ENDIF
          ENDIF
 
@@ -919,5 +1543,7 @@
          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
             !                                        ! ------------ !
-            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'un_lf'   , un_lf(:,:,:)    )
+            CALL iom_rstput( kt, nitrst, numrow, 'vn_lf'   , vn_lf(:,:,:)    )
          ENDIF
 
@@ -926,5 +1552,4 @@
 
    END SUBROUTINE dom_vvl_rst
-
 
    SUBROUTINE dom_vvl_ctl
@@ -938,7 +1563,13 @@
       INTEGER ::   ios
 
-      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
-                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
-                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
+      NAMELIST/nam_vvl/ ln_vvl_zstar               , ln_vvl_ztilde                       , &
+                      & ln_vvl_layer               , ln_vvl_ztilde_as_zstar              , &
+                      & ln_vvl_zstar_at_eqtor      , ln_vvl_zstar_on_shelf               , &
+                      & ln_vvl_adv_cn2             , ln_vvl_adv_fct                      , &
+                      & ln_vvl_lap                 , ln_vvl_blp                          , &
+                      & rn_ahe3_lap                , rn_ahe3_blp                         , &
+                      & rn_rst_e3t                 , rn_lf_cutoff                        , &
+                      & ln_vvl_regrid              , rn_zdef_max                         , &
+                      & ln_vvl_dbg   ! not yet implemented: ln_vvl_kepe
       !!---------------------------------------------------------------------- 
 
@@ -961,9 +1592,18 @@
          WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
          WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
-         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
          ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
          ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
-         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
-         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
+         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
+         WRITE(numout,*) '      ztilde on shelves          ln_vvl_zstar_on_shelf  = ', ln_vvl_zstar_on_shelf
+         WRITE(numout,*) '           Namelist nam_vvl : thickness advection scheme'
+         WRITE(numout,*) '              2nd order                  ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2
+         WRITE(numout,*) '              2nd order FCT              ln_vvl_adv_fct = ', ln_vvl_adv_fct                    
+         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion scheme'
+         WRITE(numout,*) '              Laplacian                  ln_vvl_lap     = ', ln_vvl_lap
+         WRITE(numout,*) '              Bilaplacian                ln_vvl_blp     = ', ln_vvl_blp
+         WRITE(numout,*) '              Laplacian   coefficient    rn_ahe3_lap    = ', rn_ahe3_lap
+         WRITE(numout,*) '              Bilaplacian coefficient    rn_ahe3_blp    = ', rn_ahe3_blp
+         WRITE(numout,*) '           Namelist nam_vvl : layers regriding'
+         WRITE(numout,*) '              ln_vvl_regrid                              = ', ln_vvl_regrid
          WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
          WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
@@ -992,5 +1632,11 @@
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
-      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
+
+      IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN 
+         ioptio = 0                      ! Choose one advection scheme at most
+         IF( ln_vvl_adv_cn2         )        ioptio = ioptio + 1
+         IF( ln_vvl_adv_fct         )        ioptio = ioptio + 1
+         IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' )
+      ENDIF
 
       IF(lwp) THEN                   ! Print the choice
@@ -1401,7 +2047,1069 @@
    END SUBROUTINE dom_vvl_orca_fix
 
+   SUBROUTINE dom_vvl_zdf( kt, p2dt ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_zdf  ***
+      !!
+      !! ** Purpose :   Do vertical thicknesses anomaly diffusion
+      !!
+      !! ** Method  :  
+      !!
+      !! ** Action  : 
+      !!---------------------------------------------------------------------
+      USE oce     , ONLY:   zwd => ua       , zws => va         ! (ua,va) used as 3D workspace
+      !
+      INTEGER                         , INTENT(in) ::   kt     ! ocean time-step index
+      REAL(wp)                        , INTENT(in) ::   p2dt   ! time step
+      !
+      INTEGER  :: ji, jj, jk ! dummy loop indices
+      REAL(wp) :: zr_tscale
+      REAL(wp) :: za1, za2, za3, za4, zdiff
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_zdf')
+      !
+      zr_tscale = 0.5
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt) 
+      !
+      !
+
+      zwt(:,:,1:jpk) = 1.e-10
+
+      zwt(:,:,1) = 0.0_wp
+!      DO jk = 2, jpkm1
+!         DO jj = 1, jpj
+!            DO ji = 1, jpi
+!               zwt(ji,jj,jk) = zwt(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+!            END DO
+!         END DO
+!      END DO
+      !
+      !
+      ! Set diffusivity (homogeneous to an inverse time scale)
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+! Taper a little bit:
+                za1 = tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)
+                za2 = tilde_e3t_n(ji,jj,jk  )+e3t_0(ji,jj,jk  )
+                za4 = 0.5_wp * (e3t_0(ji,jj,jk-1) + e3t_0(ji,jj,jk  ))
+                za3 = 0.5_wp * (za1 + za2) 
+                zdiff = ABS(za3-za4)/za4
+!                IF (zdiff>=0.8) THEN
+!                   zwt(ji,jj,jk) =  zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk)
+                   zwt(ji,jj,jk) =  dsm(ji,jj)/ht_0(ji,jj)*(1._wp-tanh((mbkt(ji,jj)+1-jk)*0.2))*tmask(ji,jj,jk)
+
+!                ELSE
+!                   zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk)
+!                ENDIF
+            END DO
+         END DO
+      END DO
+      ! 
+      !
+      ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+                zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / fse3w(ji,jj,jk  )
+                zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / fse3w(ji,jj,jk+1)
+                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      !! Matrix inversion from the first level
+      !!----------------------------------------------------------------------
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            zwt(ji,jj,1) = zwd(ji,jj,1)
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !         
+      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * tilde_e3t_a(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            tilde_e3t_a(ji,jj,jpkm1) = tilde_e3t_a(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 1, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               tilde_e3t_a(ji,jj,jk) = ( tilde_e3t_a(ji,jj,jk) - zws(ji,jj,jk) * tilde_e3t_a(ji,jj,jk+1) )   &
+                  &            / zwt(ji,jj,jk) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_zdf')
+      !
+   END SUBROUTINE dom_vvl_zdf
+
+   SUBROUTINE dom_vvl_zdf2( kt, p2dt ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_zdf  ***
+      !!
+      !! ** Purpose :   Do vertical interface diffusion
+      !!
+      !! ** Method  :  
+      !!
+      !! ** Action  : 
+      !!---------------------------------------------------------------------
+      USE oce     , ONLY:   zwd => ua       , zws => va         ! (ua,va) used as 3D workspace
+      !
+      INTEGER                         , INTENT(in) ::   kt     ! ocean time-step index
+      REAL(wp)                        , INTENT(in) ::   p2dt   ! time step
+      !
+      INTEGER  :: ji, jj, jk, kbot, kbotm1 ! dummy loop indices
+      REAL(wp) :: zr_tscale
+      REAL(wp) :: za1, za2, za3, za4, zdiff
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt, zw
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_zdf')
+      !
+      zr_tscale = 0.5
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zw) 
+      !
+      ! Compute internal interfaces depths:
+      zw(:,:,1) = 0._wp
+      DO jk = 2, jpk
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zw(ji,jj,jk) = zw(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))  * tmask(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      ! Set diffusivities at interfaces
+      zwt(:,:,:) = 0.00000001_wp * tmask(:,:,:)      
+      zwt(:,:,1) = 0._wp
+      ! 
+      !
+      ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+                zwi(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk-1) + zwt(ji,jj,jk  )) & 
+                              &                   / fse3w(ji,jj,jk-1) / fse3t(ji,jj,jk-1) &
+                              &                   * ht(ji,jj) 
+                zws(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk  ) + zwt(ji,jj,jk+1)) & 
+                              &                   / fse3w(ji,jj,jk  ) / fse3t(ji,jj,jk  ) &
+                              &                   * ht(ji,jj) 
+
+                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      ! Boundary conditions (Neumann)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            zwi(ji,jj,1) = 0._wp
+            zws(ji,jj,1) = 0._wp
+            zwd(ji,jj,1) = 1._wp
+            zw (ji,jj,1) = 0._wp
+            !
+            zwd(ji,jj,2) = zwd(ji,jj,2) +  zwi(ji,jj,2)
+            zwi(ji,jj,2) = 0._wp
+!            zwi(ji,jj,2) = 0._wp
+!            zws(ji,jj,2) = 0._wp
+!            zwd(ji,jj,2) = 1._wp
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            kbot   = mbkt(ji,jj) + 1
+            kbotm1 = mbkt(ji,jj)
+            zwi(ji,jj,kbot  ) = 0._wp
+            zws(ji,jj,kbot  ) = 0._wp
+            zwd(ji,jj,kbot  ) = 1._wp
+            !
+            zwd(ji,jj,kbotm1) = zwd(ji,jj,kbotm1) +  zws(ji,jj,kbotm1)
+            zws(ji,jj,kbotm1) = 0._wp
+!            zwi(ji,jj,kbotm1) = 0._wp
+!            zws(ji,jj,kbotm1) = 0._wp
+!            zwd(ji,jj,kbotm1) = 1._wp
+         END DO
+      END DO
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jk = 2, jpk
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zwi(ji,jj,jk) = zw(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) *zwi(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jk = jpk-1, 2, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zw(ji,jj,jk) = ( zwi(ji,jj,jk) - zws(ji,jj,jk) * zw(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      ! Revert to thicknesses anomalies:
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               tilde_e3t_a(ji,jj,jk) = (zw(ji,jj,jk+1)-zw(ji,jj,jk)-e3t_0(ji,jj,jk))* tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zw) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_zdf')
+      !
+   END SUBROUTINE dom_vvl_zdf2
+
+   SUBROUTINE dom_vvl_regrid( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_regrid  ***
+      !!                   
+      !! ** Purpose :  Ensure "well-behaved" vertical grid
+      !!
+      !! ** Method  :  More or less adapted from references below.
+      !!
+      !! ** Action  :  Ensure that thickness are above a given value, spaced enough
+      !!               and revert to Eulerian coordinates near the bottom.      
+      !!
+      !! References :  Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction
+      !!               with a Model Combining Terrain-following and Isentropic
+      !!               coordinates. Part I: Model Description. Monthly Weather Rev.,
+      !!               121, 1770-1785.
+      !!               Toy, M., 2011: Incorporating Condensational Heating into a
+      !!               Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic-
+      !!               Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in )               :: kt         ! time step
+
+      !! * Local declarations
+      INTEGER                             :: ji, jj, jk ! dummy loop indices
+      LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond
+      LOGICAL                             :: ll_e3tsurf_const, ll_zdiff_cond, ll_blpdiff_cond
+      INTEGER                             :: jkbot
+      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd, dzmin_vvl
+      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj
+      REAL(wp)                            :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn
+      REAL(wp)                            :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot
+      REAL(wp)                            :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zdw
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zwdw_b
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_regrid')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zwdw) 
+      !
+      ! Some user defined parameters below:
+      ll_chk_bot2top   = .TRUE.
+      ll_chk_top2bot   = .TRUE.
+      ll_e3tsurf_const = .FALSE.
+      dzmin_vvl  = 1._wp   ! Absolute minimum depth (in meters)
+      zfrch_stp  = 0.5_wp   ! Maximum fractionnal thickness change in one time step (<= 1.)
+      zfrch_rel  = 0.0_wp   ! Maximum relative thickness change in the vertical (<= 1.)
+      zfrac_bot  = 0.01_wp  ! Fraction of bottom level allowed to change
+      zscal_bot  = 2._wp    ! Depth lengthscale
+      ll_zdiff_cond    = .FALSE.  ! Conditionnal vertical diffusion of interfaces
+         zvdiff  = 0.1_wp   ! m/s
+         zvlim   = 0.4_wp   ! max d2h/dh
+      ll_lapdiff_cond  = .FALSE.  ! Conditionnal Laplacian diffusion of interfaces
+         zhdiff  = 0.1_wp   ! ad.
+         zhlim   = 0.01_wp  ! max lap(z)/e1
+      ll_blpdiff_cond  = .FALSE.  ! Conditionnal Bilaplacian diffusion of interfaces
+         zhdiff2 = 0.8_wp   ! ad.
+         zhlim2  = 5.e-11_wp  ! max bilap(z)/e1**3
+      ! --------------------------------------------------------------------------------------- 
+      !
+      ! Set arrays determining maximum vertical displacement at the bottom:
+      !--------------------------------------------------------------------
+      IF ( kt==nit000 ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1))
+               i_int_bot(ji,jj) = jk
+            END DO
+         END DO
+         dsm(:,:) = REAL( i_int_bot(:,:), wp )   ;   CALL lbc_lnk(dsm(:,:),'T',1.)
+         i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 )
+
+         CALL wrk_alloc( jpi, jpj, zdw )
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji  ,jj,1), &
+                          &     ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), &
+                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj  ,1), &
+                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1)  )
+               zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall ) 
+            END DO
+         END DO
+         CALL lbc_lnk( zdw(:,:), 'T', 1. )
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               dsm(ji,jj) = 1._wp/16._wp * (            zdw(ji-1,jj-1) + zdw(ji+1,jj-1)      &
+                  &                           +         zdw(ji-1,jj+1) + zdw(ji+1,jj+1)      &
+                  &                           + 2._wp*( zdw(ji  ,jj-1) + zdw(ji-1,jj  )      &
+                  &                           +         zdw(ji+1,jj  ) + zdw(ji  ,jj+1) )    &
+                  &                           + 4._wp*  zdw(ji  ,jj  )                       )
+            END DO
+         END DO
+
+         CALL lbc_lnk( dsm(:,:), 'T', 1. )
+         CALL wrk_dealloc( jpi, jpj, zdw)         
+     
+         IF (ln_zps) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  jk = i_int_bot(ji,jj)
+                  hsm(ji,jj) = zfrac_bot * e3w_1d(jk)
+               END DO
+            END DO
+         ELSE
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  jk = i_int_bot(ji,jj)
+                  hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+      END IF
+
+      ! Provisionnal interface depths:
+      !-------------------------------
+      zwdw(:,:,1) = 0.e0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            DO jk = 2, jpk
+               zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + & 
+                              & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      ! Conditionnal horizontal Laplacian diffusion:
+      !---------------------------------------------
+      IF ( ll_lapdiff_cond ) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zwdw_b)
+         !
+         zwdw_b(:,:,1) = 0._wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk=2,jpk
+                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
+                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+               END DO
+            END DO
+         END DO
+         !
+         DO jk = 2, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.                  
+                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
+                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+            
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
+                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
+                  zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e12t(ji,jj)), 0._wp)
+                  ztmp = SIGN(zh2, ztmp1)
+                  zeu2 = zhdiff * e12t(ji,jj)*e12t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj))
+                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO         
+         !
+         CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b)
+         !
+      ENDIF
+
+      ! Conditionnal horizontal Bilaplacian diffusion:
+      !-----------------------------------------------
+      IF ( ll_blpdiff_cond ) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zwdw_b)
+         !
+         zwdw_b(:,:,1) = 0._wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk=2,jpk
+                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
+                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
+               END DO
+            END DO
+         END DO
+         !
+         DO jk = 2, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.                  
+                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
+                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+            
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zwdw_b(ji,jj,jk) = -( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
+                                &  +    (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
+               END DO
+            END DO
+         END DO   
+         !
+         CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )      
+         !
+         DO jk = 2, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.                  
+                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
+                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
+                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
+                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+         !   
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
+                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
+                  zh2 = MAX(abs(ztmp1)-zhlim2, 0._wp)
+                  ztmp = SIGN(zh2, ztmp1)
+                  zeu2 = zhdiff2 * e12t(ji,jj)*e12t(ji,jj) / 16._wp
+                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO   
+         !
+         CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b)
+         !
+      ENDIF
+
+      ! Conditionnal vertical diffusion:
+      !---------------------------------
+      IF ( ll_zdiff_cond ) THEN
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.    
+                  ztmp  = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) & 
+                            -(tilde_e3t_b(ji,jj,jk  )+e3t_0(ji,jj,jk  ))*tmask(ji,jj,jk  ) ) 
+                  ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) & 
+                        &           +tilde_e3t_b(ji,jj,jk  ) + e3t_0(ji,jj,jk  ) )
+                  zh2   = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp)
+                  ztmp  = SIGN(zh2, ztmp)
+                  IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0
+                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      ! Check grid from the bottom to the surface
+      !------------------------------------------
+      IF ( ll_chk_bot2top ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               jkbot = mbkt(ji,jj)                   
+               DO jk = jkbot,2,-1
+                  !
+                  zh_0   = e3t_0(ji,jj,jk)
+                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1))
+                  zh_old = tilde_e3t_a(ji,jj,jk  ) + zh_0
+!                  zh_dwn = tilde_e3t_a(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)
+                  zh_min = MIN(zh_0/3._wp, dzmin_vvl)
+                  ! 
+                  ! Set maximum and minimum vertical excursions   
+                  ztmph = hsm(ji,jj)
+                  ztmpd = dsm(ji,jj)
+                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)               
+                  zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 )
+                  zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff)
+                  zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 )
+                  zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff)
+                  !
+                  ! New layer thickness:
+                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
+                  !
+                  ! Ensure minimum layer thickness:                  
+!                  zh_new = MIN(zh_old, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) )
+                  zh_new = cush(zh_new, zh_min)
+                  !
+                  ! Final flux:
+                  zdiff = (zh_new - zh_old) * tmask(ji,jj,jk)
+                  !
+                  ! Limit flux: 
+                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
+                  zdiff = SIGN(ztmp, zh_new - zh_old)
+                  zh_new = zdiff + zh_old
+                  !
+                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
+                  zwdw(ji,jj,jk)          = zwdw(ji,jj,jk+1) - zh_new
+                  tilde_e3t_a(ji,jj,jk-1) = (-zdiff + tilde_e3t_a(ji,jj,jk-1) ) * tmask(ji,jj,jk-1)            
+               END DO    
+            END DO 
+         END DO
+      END IF    
+      !
+      ! Check grid from the surface to the bottom 
+      !------------------------------------------ 
+      IF ( ll_chk_top2bot ) THEN      
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               jkbot = mbkt(ji,jj)   
+               DO jk = 1, jkbot-1
+                  !
+                  zh_0   = e3t_0(ji,jj,jk)
+                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1))
+                  zh_old = tilde_e3t_a(ji,jj,jk  ) + zh_0   
+!                  zh_up  = tilde_e3t_a(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)
+                  zh_min = MIN(zh_0/3._wp, dzmin_vvl)
+                  IF ((jk<=5).AND.ll_e3tsurf_const) zh_min = MAX(e3t_0(ji,jj,1)/3._wp, dzmin_vvl)
+                  !
+                  ! Ensure minimum layer thickness:
+!                  zh_new=MIN(zh_old, zh_up * zfrch_rel / (2._wp-zfrch_rel) )
+                  zh_new = cush(zh_old, zh_min)
+                  !
+                  ! Final flux:
+                  zdiff = (zh_new -zh_old) * tmask(ji,jj,jk)
+                  ! 
+                  ! Limit flux:                 
+                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
+                  zdiff = SIGN(ztmp, zdiff)
+                  zh_new = zdiff + zh_old
+                  !
+                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
+                  zwdw(ji,jj,jk+1)        = zwdw(ji,jj,jk) + zh_new
+                  tilde_e3t_a(ji,jj,jk+1) = (-zdiff + tilde_e3t_a(ji,jj,jk+1) ) * tmask(ji,jj,jk+1)
+               END DO
+               !               
+            END DO
+         END DO
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zwdw ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_regrid')
+      !
+   END SUBROUTINE dom_vvl_regrid
+
+   FUNCTION cush(hin, hmin)  RESULT(hout)
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION cush  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  hin, hmin
+      REAL(wp)             ::  hout, zx, zh_cri
+      !!----------------------------------------------------------------------
+      zh_cri = 3._wp * hmin
+      !
+      IF ( hin<=0._wp ) THEN
+         hout = hmin
+      !
+      ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN
+         zx = hin/zh_cri
+         hout = hmin * (1._wp + zx + zx*zx)      
+      !
+      ELSEIF ( hin>zh_cri ) THEN
+         hout = hin
+      !
+      ENDIF
+      !
+   END FUNCTION cush
+
+   FUNCTION cush_max(hin, hmax)  RESULT(hout)
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION cush  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  hin, hmax
+      REAL(wp)             ::  hout, hmin, zx, zh_cri
+      !!----------------------------------------------------------------------
+      hmin = 0.1_wp * hmax
+      zh_cri = 3._wp * hmin
+      !
+      IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN
+         zx = (hmax-hin)/zh_cri
+         hout = hmax - hmin * (1._wp + zx + zx*zx)
+         !
+      ELSEIF ( hin>(hmax-zh_cri) ) THEN
+         hout = hmax - hmin
+         !
+      ELSE
+         hout = hin
+         !
+      ENDIF
+      !
+   END FUNCTION cush_max
+
+   SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_vvl_adv_fct  ***
+      !! 
+      !! **  Purpose :  Do thickness advection
+      !!
+      !! **  Method  :  FCT scheme to ensure positivity 
+      !!
+      !! **  Action  : - Update pta thickness tendency and diffusive fluxes
+      !!               - this is the total trend, hence it does include sea level motions
+      !!               - Upstream corrections to antidiffusive fluxes ensure
+      !!                 that barotropic transport matches what is contained in input fluxes
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta    ! thickness baroclinic trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)    ::   uin, vin  ! input velocities
+      !
+      INTEGER  ::   ji, jj, jk               ! dummy loop indices  
+      INTEGER  ::   ikbu, ikbv, ibot
+      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
+      REAL(wp) ::   zdi, zdj, zmin           !   -      -
+      REAL(wp) ::   zfp_ui, zfp_vj           !   -      -
+      REAL(wp) ::   zfm_ui, zfm_vj           !   -      -
+      REAL(wp) ::   zfp_hi, zfp_hj           !   -      -
+      REAL(wp) ::   zfm_hi, zfm_hj           !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zbu, zbv, zhu_b, zhv_b
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwi
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_adv_fct')
+      !
+      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
+      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwi)
+      !
+      ! 1. Initializations
+      ! ------------------
+      !
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+         z2dtt =  rdt
+      ELSE
+         z2dtt = 2.0_wp * rdt
+      ENDIF
+      !
+      zwi(:,:,:) = 0.e0
+      zwx(:,:,:) = 0.e0 
+      zwy(:,:,:) = 0.e0
+      !
+      !
+      ! 2. upstream advection with initial mass fluxes & intermediate update
+      ! --------------------------------------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               !
+               zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
+               zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
+               zfp_hi = fse3t_b(ji  ,jj  ,jk)
+               zfm_hi = fse3t_b(ji+1,jj  ,jk)
+               zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
+               !
+               zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
+               zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
+               zfp_hj = fse3t_b(ji  ,jj  ,jk)
+               zfm_hj = fse3t_b(ji  ,jj+1,jk)
+               zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      IF ( .NOT.ln_sco ) THEN
+         ! Correct bottom upstream fluxes
+         ! considering "shelf horizon depths"
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               ! upstream scheme
+               ikbu = mbku(ji,jj)
+               ikbv = mbkv(ji,jj)
+               zfp_ui = uin(ji,jj,ikbu) + ABS( uin(ji,jj,ikbu) )
+               zfm_ui = uin(ji,jj,ikbu) - ABS( uin(ji,jj,ikbu) )
+               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
+               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
+               zwx(ji,jj,ikbu) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,ikbu)
+               !
+               zfp_vj = vin(ji,jj,ikbv) + ABS( vin(ji,jj,ikbv) )
+               zfm_vj = vin(ji,jj,ikbv) - ABS( vin(ji,jj,ikbv) )
+               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
+               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
+               zwy(ji,jj,ikbv) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,ikbv)
+            END DO
+         END DO
+      ENDIF
+
+      ! total advective trend
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zbtr = r1_e12t(ji,jj)
+               ! total intermediate advective trends
+               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )  &
+                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )  )
+               !
+               ! update and guess with monotonic sheme
+               pta(ji,jj,jk) =   pta(ji,jj,jk) + ztra
+               zwi(ji,jj,jk) = (fse3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !                             ! Lateral boundary conditions on zwi  (unchanged sign)
+      CALL lbc_lnk( zwi, 'T', 1. )  
+
+      IF ( ln_vvl_dbg ) THEN
+         zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 ) 
+         IF( lk_mpp )   CALL mpp_min( zmin )
+         IF( zmin < 0._wp) THEN
+            IF(lwp) CALL ctl_stop('vvl_adv: CFL issue here')
+         ENDIF
+      ENDIF
+
+      ! 3. antidiffusive flux : high order minus low order
+      ! --------------------------------------------------
+      ! antidiffusive flux on i and j
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * fse3u_n(ji,jj,jk) &
+                                & - zwx(ji,jj,jk)) * umask(ji,jj,jk)
+               zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * fse3v_n(ji,jj,jk) & 
+                                & - zwy(ji,jj,jk)) * vmask(ji,jj,jk)
+               !
+               ! Update advective fluxes
+               un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk)
+               vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      
+      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
+
+      ! 4. monotonicity algorithm
+      ! -------------------------
+      CALL nonosc_2d( fse3t_b(:,:,:), zwx, zwy, zwi, z2dtt )
+
+      ! 5. final trend with corrected fluxes
+      ! ------------------------------------
+      !
+      ! Update advective fluxes
+      un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:)
+      vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:)
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.  
+               !
+               zbtr = r1_e12t(ji,jj)
+               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
+               ! add them to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
+            END DO
+         END DO
+      END DO
+
+      !
+      ! 6. Correct barotropic flux
+      ! --------------------------
+      ! Compute barotropic flux difference:
+      zbu(:,:) = 0.e0
+      zbv(:,:) = 0.e0
+      DO jj = 1, jpj
+         DO ji = 1, jpi   ! vector opt.
+            DO jk = 1, jpkm1
+               zbu(ji,jj) = zbu(ji,jj) - un_td(ji,jj,jk) * umask(ji,jj,jk)
+               zbv(ji,jj) = zbv(ji,jj) - vn_td(ji,jj,jk) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      ENDDO
+
+      ! Compute upstream depths:
+      zhu_b(:,:) = 0.e0
+      zhv_b(:,:) = 0.e0
+      IF ( ln_sco ) THEN; ibot=0 ; ELSE ; ibot=1 ; ENDIF
+
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1   ! vector opt.
+            zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
+            ikbu = mbku(ji,jj) - ibot 
+            DO jk = 1, ikbu
+               zfp_hi = fse3t_b(ji  ,jj  ,jk)
+               zfm_hi = fse3t_b(ji+1,jj  ,jk)
+               zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi &
+                            &                  + (1._wp-zdi) * zfm_hi &
+                            &                ) * umask(ji,jj,jk)
+            END DO
+            !
+            zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
+            ikbv = mbkv(ji,jj) - ibot
+            DO jk = 1, ikbv
+               zfp_hj = fse3t_b(ji  ,jj  ,jk)
+               zfm_hj = fse3t_b(ji  ,jj+1,jk)
+               zhv_b(ji,jj) = zhv_b(ji,jj) + (          zdj  * zfp_hj &
+                            &                  + (1._wp-zdj) * zfm_hj &
+                            &                ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      IF ( .NOT.ln_sco ) THEN
+         ! Correct bottom value
+         ! considering "shelf horizon depth"
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1   ! vector opt.
+               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
+               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
+               ikbu = mbku(ji,jj)
+               ikbv = mbkv(ji,jj)
+               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
+               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
+               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
+               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
+               zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi &
+                            &                 + (1._wp-zdi) * zfm_hi &
+                            &                ) * umask(ji,jj,ikbu)
+               zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj &
+                            &                 + (1._wp-zdj) * zfm_hj &
+                            &                ) * vmask(ji,jj,ikbv)
+            END DO
+         END DO
+      ENDIF
+
+      ! Corrective barotropic velocity (times hor. scale factor)
+      zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1))
+      zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1))
+
+      CALL lbc_lnk( zbu(:,:), 'U', -1. )
+      CALL lbc_lnk( zbv(:,:), 'V', -1. )
+      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
+      
+      ! Set corrective fluxes in upstream direction:
+      !
+      zwx(:,:,:) = 0.e0
+      zwy(:,:,:) = 0.e0
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            ! upstream scheme
+            zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
+            zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
+            zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
+            zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
+            DO jk = 1, jpkm1
+               zfp_hi = fse3t_b(ji  ,jj  ,jk)
+               zfm_hi = fse3t_b(ji+1,jj  ,jk)
+               !
+               zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
+
+               zfp_hj = fse3t_b(ji  ,jj  ,jk)
+               zfm_hj = fse3t_b(ji  ,jj+1,jk)
+               !
+               zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      IF ( .NOT.ln_sco ) THEN
+      ! Bottom correction:
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            ! upstream scheme
+            ikbu = mbku(ji,jj)
+            ikbv = mbkv(ji,jj)
+            !
+            zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
+            zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
+            zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
+            zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
+            !
+            zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
+            zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
+            !
+            zwx(ji,jj,ikbu) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi )
+            !
+            zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
+            zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
+            !
+            zwy(ji,jj,ikbv) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj )
+         END DO
+      END DO
+      ENDIF
+
+      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
+
+      un_td(:,:,:) = un_td(:,:,:) + zwx(:,:,:)
+      vn_td(:,:,:) = vn_td(:,:,:) + zwy(:,:,:)
+      !
+      ! Update trend with corrective fluxes:
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.  
+               !
+               zbtr = r1_e12t(ji,jj)
+               ! total advective trends
+               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
+               ! add them to the general tracer trends
+               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
+            END DO
+         END DO
+      END DO
+
+      CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
+      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi)
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_adv_fct')
+      !
+   END SUBROUTINE dom_vvl_adv_fct
+
+   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE nonosc_2d  ***
+      !!     
+      !! **  Purpose :   compute monotonic thickness fluxes from the upstream 
+      !!       scheme and the before field by a nonoscillatory algorithm 
+      !!
+      !! **  Method  :   ... ???
+      !!       warning : pbef and paft must be masked, but the boundaries
+      !!       conditions on the fluxes are not necessary zalezak (1979)
+      !!       drange (1995) multi-dimensional forward-in-time and upstream-
+      !!       in-space based differencing for fluid
+      !!----------------------------------------------------------------------
+      !
+      !!----------------------------------------------------------------------
+      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb        ! monotonic fluxes in the 3 directions
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
+      REAL(wp) ::   zupip1, zupim1, zupjp1, zupjm1, zupb, zupa
+      REAL(wp) ::   zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('nonosc2')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
+      !
+
+      zbig  = 1.e+40_wp
+      zrtrn = 1.e-15_wp
+      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp
+
+
+      ! Search local extrema
+      ! --------------------
+      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
+      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   &
+         &        paft * tmask - zbig * ( 1.e0 - tmask )  )
+      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   &
+         &        paft * tmask + zbig * ( 1.e0 - tmask )  )
+
+      DO jk = 1, jpkm1
+         z2dtt = p2dt
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+
+               ! search maximum in neighbourhood
+               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
+                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
+                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ))
+
+               ! search minimum in neighbourhood
+               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
+                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
+                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ))
+
+               ! positive part of the flux
+               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
+                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   
+
+               ! negative part of the flux
+               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
+                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )
+
+               ! up & down beta terms
+               zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt
+               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
+               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
+
+      ! 3. monotonic flux in the i & j direction (paa & pbb)
+      ! ----------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
+               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
+               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
+               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
+
+               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
+               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
+               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
+               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('nonosc2')
+      !
+   END SUBROUTINE nonosc_2d
+
    !!======================================================================
 END MODULE domvvl
-
-
-
Index: /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 9352)
+++ /branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 9353)
@@ -260,5 +260,5 @@
             CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
             CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
-            CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+            CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' )
             ! Vertical scale factor interpolations
             ! ------------------------------------
