Index: NEMO/releases/r4.0/r4.0-HEAD/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/cfgs/SHARED/namelist_ref	(revision 15609)
+++ NEMO/releases/r4.0/r4.0-HEAD/cfgs/SHARED/namelist_ref	(revision 15610)
@@ -877,4 +877,8 @@
    rn_zdef_max   =  0.9             !  maximum fractional e3t deformation
    ln_vvl_dbg    = .true.           !  debug prints    (T/F)
+   nn_vvl_interp =  2               !  interpolation method of scale factor anomalies at U/V/F points
+                                    !  =0 linear even at the bottom (old)
+                                    !  =1 linear with bottom correction
+                                    !  =2 proportionnal to scale factors at rest ("qco" like) 
 /
 !-----------------------------------------------------------------------
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/dom_oce.F90	(revision 15609)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/dom_oce.F90	(revision 15610)
@@ -149,4 +149,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hf_0                                  !: f-depth              [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
@@ -275,5 +276,5 @@
          &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
          !
-      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           &
+      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , hf_0(jpi,jpj)    ,                        &
          &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     &
          &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     &
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/domain.F90	(revision 15609)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/domain.F90	(revision 15610)
@@ -150,4 +150,5 @@
       hu_0(:,:) = 0._wp
       hv_0(:,:) = 0._wp
+      hf_0(:,:) = 0._wp
       DO jk = 1, jpk
          ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
@@ -155,4 +156,13 @@
          hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpim1  
+               hf_0(ji,jj) = hf_0(ji,jj) + e3f_0(ji,jj,jk)*vmask(ji,jj,jk)*vmask(ji+1,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk('domain', hf_0, 'F', 1._wp)
       !
       !           !==  time varying part of coordinate system  ==!
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/domvvl.F90	(revision 15609)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/domvvl.F90	(revision 15610)
@@ -49,4 +49,10 @@
    LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
    LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
+   !
+   INTEGER          :: nn_vvl_interp                       ! scale factors anomaly interpolation method at U-V-F points
+                                                           ! =0 linear with no bottom correction over steps (old)
+                                                           ! =1 linear with bottom correction over steps
+                                                           ! =2 "qco like", i.e. proportional to thicknesses at rest
+   !
    !                                                       ! conservation: not used yet
    REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
@@ -720,5 +726,6 @@
       INTEGER ::   ji, jj, jk                                       ! dummy loop indices
       INTEGER ::   iku, ikum1, ikv, ikvm1, ikf, ikfm1               ! 
-      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
+      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/
+      REAL(wp), DIMENSION(jpi,jpj) :: zssh                          ! work array to retrieve ssh (nn_vvl_interp > 1)
       !!----------------------------------------------------------------------
       !
@@ -732,26 +739,55 @@
          !
       CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )   
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE ( 1 ) 
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
-                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               iku    = mbku(ji  ,jj)
-               ikum1  = iku - 1
-               pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              & 
-                  &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            &
-                  &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   &
-                  &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikum1)))
-            END DO
-         END DO
+                  iku    = mbku(ji  ,jj)
+                  ikum1  = iku - 1
+                  pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              & 
+                     &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            &
+                     &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   &
+                     &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikum1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)    &
+                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji+1,jj) * zssh(ji+1,jj)) &
+                        &                       * e3u_0(ji,jj,jk) / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !   
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
@@ -759,26 +795,55 @@
          !
       CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )   
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !                     !
+         CASE ( 1 ) 
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
-                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ikv    = mbkv(ji  ,jj)
-               ikvm1  = ikv - 1
-               pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              & 
-                  &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            &
-                  &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   &
-                  &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikvm1)))
-            END DO
-         END DO
+                  ikv    = mbkv(ji  ,jj)
+                  ikvm1  = ikv - 1
+                  pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              & 
+                     &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            &
+                     &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   &
+                     &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikvm1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)    &
+                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji,jj+1) * zssh(ji,jj+1)) &
+                        &                       * e3v_0(ji,jj,jk) / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !   
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
@@ -786,27 +851,57 @@
          !
       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )  
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                        &                       *    r1_e1e2f(ji,jj)                                                  &
+                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE ( 1 ) 
+            !
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                        &                       *    r1_e1e2f(ji,jj)                                                  &
+                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
-                     &                       *    r1_e1e2f(ji,jj)                                                  &
-                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1))
-               ikfm1  = ikf - 1
-               pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         &
-                  &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            &
-                  &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   &
-                  &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikfm1)))
-            END DO
-         END DO
+                  ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1))
+                  ikfm1  = ikf - 1
+                  pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         &
+                     &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            &
+                     &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   &
+                     &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikfm1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(umask(:,:,:)*(pe3_in(:,:,:)-e3u_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) =  (  umask(ji,jj,jk)* umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )   &
+                        &                 * 0.5_wp * r1_e1e2f(ji,jj)                                           &
+                        &                 * (e1e2u(ji  ,jj) * zssh(ji  ,jj) + e1e2u(ji,jj+1) * zssh(ji,jj+1))  &
+                        &                 * e3f_0(ji,jj,jk) / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) )
+                  END DO
+               END DO
+            END DO      
+            !
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
@@ -1059,5 +1154,6 @@
       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
+         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg            , &! not yet implemented: ln_vvl_kepe
+         &              nn_vvl_interp
       !!---------------------------------------------------------------------- 
       !
@@ -1095,4 +1191,5 @@
          ENDIF
          WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
+         WRITE(numout,*) '         Method to compute scale factors anomaly at U/V/F points  nn_vvl_interp   = ', nn_vvl_interp
       ENDIF
       !
@@ -1106,4 +1203,5 @@
       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
+      IF( .NOT. ln_vvl_zstar .AND. (nn_vvl_interp==2 ) )  CALL ctl_stop( 'nn_vvl_interp must be < 2 if ln_vvl_zstar=F' )
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
Index: NEMO/releases/r4.0/r4.0-HEAD/tests/CANAL/MY_SRC/domvvl.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/tests/CANAL/MY_SRC/domvvl.F90	(revision 15609)
+++ NEMO/releases/r4.0/r4.0-HEAD/tests/CANAL/MY_SRC/domvvl.F90	(revision 15610)
@@ -49,4 +49,10 @@
    LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
    LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
+   !
+   INTEGER          :: nn_vvl_interp                       ! scale factors anomaly interpolation method at U-V-F points
+                                                           ! =0 linear with no bottom correction over steps (old)
+                                                           ! =1 linear with bottom correction over steps
+                                                           ! =2 "qco like", i.e. proportional to thicknesses at rest
+   !
    !                                                       ! conservation: not used yet
    REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
@@ -720,5 +726,6 @@
       INTEGER ::   ji, jj, jk                                       ! dummy loop indices
       INTEGER ::   iku, ikum1, ikv, ikvm1, ikf, ikfm1               ! 
-      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
+      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/
+      REAL(wp), DIMENSION(jpi,jpj) :: zssh                          ! work array to retrieve ssh (nn_vvl_interp > 1)
       !!----------------------------------------------------------------------
       !
@@ -732,26 +739,55 @@
          !
       CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )   
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE ( 1 ) 
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
-                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               iku    = mbku(ji  ,jj)
-               ikum1  = iku - 1
-               pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              & 
-                  &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            &
-                  &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   &
-                  &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikum1)))
-            END DO
-         END DO
+                  iku    = mbku(ji  ,jj)
+                  ikum1  = iku - 1
+                  pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              & 
+                     &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            &
+                     &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   &
+                     &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikum1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)    &
+                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji+1,jj) * zssh(ji+1,jj)) &
+                        &                       * e3u_0(ji,jj,jk) / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !   
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
@@ -759,26 +795,55 @@
          !
       CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )   
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !                     !
+         CASE ( 1 ) 
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
-                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ikv    = mbkv(ji  ,jj)
-               ikvm1  = ikv - 1
-               pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              & 
-                  &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            &
-                  &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   &
-                  &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikvm1)))
-            END DO
-         END DO
+                  ikv    = mbkv(ji  ,jj)
+                  ikvm1  = ikv - 1
+                  pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              & 
+                     &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            &
+                     &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   &
+                     &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikvm1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)    &
+                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji,jj+1) * zssh(ji,jj+1)) &
+                        &                       * e3v_0(ji,jj,jk) / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !   
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
@@ -786,27 +851,57 @@
          !
       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )  
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                        &                       *    r1_e1e2f(ji,jj)                                                  &
+                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE ( 1 ) 
+            !
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                        &                       *    r1_e1e2f(ji,jj)                                                  &
+                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
-                     &                       *    r1_e1e2f(ji,jj)                                                  &
-                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1))
-               ikfm1  = ikf - 1
-               pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         &
-                  &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            &
-                  &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   &
-                  &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikfm1)))
-            END DO
-         END DO
+                  ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1))
+                  ikfm1  = ikf - 1
+                  pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         &
+                     &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            &
+                     &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   &
+                     &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikfm1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(umask(:,:,:)*(pe3_in(:,:,:)-e3u_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) =  (  umask(ji,jj,jk)* umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )   &
+                        &                 * 0.5_wp * r1_e1e2f(ji,jj)                                           &
+                        &                 * (e1e2u(ji  ,jj) * zssh(ji  ,jj) + e1e2u(ji,jj+1) * zssh(ji,jj+1))  &
+                        &                 * e3f_0(ji,jj,jk) / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) )
+                  END DO
+               END DO
+            END DO      
+            !
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
@@ -997,5 +1092,5 @@
                !
             ELSE
-               !
+              !
                ! usr_def_istate called here only to get sshb, that is needed to
                ! initialize e3t_b and e3t_n 
@@ -1061,5 +1156,6 @@
       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
+         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg            , &! not yet implemented: ln_vvl_kepe
+         &              nn_vvl_interp
       !!---------------------------------------------------------------------- 
       !
@@ -1097,4 +1193,5 @@
          ENDIF
          WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
+         WRITE(numout,*) '         Method to compute scale factors anomaly at U/V/F points  nn_vvl_interp   = ', nn_vvl_interp
       ENDIF
       !
@@ -1108,4 +1205,5 @@
       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
+      IF( .NOT. ln_vvl_zstar .AND. (nn_vvl_interp==2 ) )  CALL ctl_stop( 'nn_vvl_interp must be < 2 if ln_vvl_zstar=F' )
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
Index: NEMO/releases/r4.0/r4.0-HEAD/tests/VORTEX/MY_SRC/domvvl.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/tests/VORTEX/MY_SRC/domvvl.F90	(revision 15609)
+++ NEMO/releases/r4.0/r4.0-HEAD/tests/VORTEX/MY_SRC/domvvl.F90	(revision 15610)
@@ -49,4 +49,10 @@
    LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
    LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
+   !
+   INTEGER          :: nn_vvl_interp                       ! scale factors anomaly interpolation method at U-V-F points
+                                                           ! =0 linear with no bottom correction over steps (old)
+                                                           ! =1 linear with bottom correction over steps
+                                                           ! =2 "qco like", i.e. proportional to thicknesses at rest
+   !
    !                                                       ! conservation: not used yet
    REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
@@ -720,5 +726,6 @@
       INTEGER ::   ji, jj, jk                                       ! dummy loop indices
       INTEGER ::   iku, ikum1, ikv, ikvm1, ikf, ikfm1               ! 
-      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
+      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/
+      REAL(wp), DIMENSION(jpi,jpj) :: zssh                          ! work array to retrieve ssh (nn_vvl_interp > 1)
       !!----------------------------------------------------------------------
       !
@@ -732,26 +739,55 @@
          !
       CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )   
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE ( 1 ) 
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
-                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               iku    = mbku(ji  ,jj)
-               ikum1  = iku - 1
-               pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              & 
-                  &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            &
-                  &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   &
-                  &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikum1)))
-            END DO
-         END DO
+                  iku    = mbku(ji  ,jj)
+                  ikum1  = iku - 1
+                  pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd )                              & 
+                     &      *(  0.5_wp * r1_e1e2u(ji,jj)                                                            &
+                     &      *(  e1e2t(ji  ,jj) * ( SUM(tmask(ji  ,jj,:)*(pe3_in(ji  ,jj,:) - e3t_0(ji  ,jj,:))) )   &
+                     &        + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikum1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)    &
+                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji+1,jj) * zssh(ji+1,jj)) &
+                        &                       * e3u_0(ji,jj,jk) / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !   
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
@@ -759,26 +795,55 @@
          !
       CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )   
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !                     !
+         CASE ( 1 ) 
+            !   
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
-                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ikv    = mbkv(ji  ,jj)
-               ikvm1  = ikv - 1
-               pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              & 
-                  &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            &
-                  &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   &
-                  &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikvm1)))
-            END DO
-         END DO
+                  ikv    = mbkv(ji  ,jj)
+                  ikvm1  = ikv - 1
+                  pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd )                              & 
+                     &      *(  0.5_wp * r1_e1e2v(ji,jj)                                                            &
+                     &      *(  e1e2t(ji,jj  ) * ( SUM(tmask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3t_0(ji,jj  ,:))) )   &
+                     &        + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikvm1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)    &
+                        &                       * (   e1e2t(ji  ,jj) * zssh(ji  ,jj) + e1e2t(ji,jj+1) * zssh(ji,jj+1)) &
+                        &                       * e3v_0(ji,jj,jk) / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )
+                  END DO
+               END DO
+            END DO
+            !   
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
@@ -786,27 +851,57 @@
          !
       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         SELECT CASE ( nn_vvl_interp )
+         CASE ( 0 )  
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                        &                       *    r1_e1e2f(ji,jj)                                                  &
+                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE ( 1 ) 
+            !
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                        &                       *    r1_e1e2f(ji,jj)                                                  &
+                        &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                        &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+            !
+            ! Bottom correction:
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
-                     &                       *    r1_e1e2f(ji,jj)                                                  &
-                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! Bottom correction:
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1))
-               ikfm1  = ikf - 1
-               pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         &
-                  &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            &
-                  &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   &
-                  &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
-                  &     - SUM(pe3_out(ji,jj,1:ikfm1)))
-            END DO
-         END DO
+                  ikf    = MIN(mbku(ji  ,jj),mbku(ji,jj+1))
+                  ikfm1  = ikf - 1
+                  pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd )         &
+                     &     * ( 0.5_wp *  r1_e1e2f(ji,jj)                                                            &
+                     &     * (  e1e2u(ji,jj  ) * ( SUM(umask(ji,jj  ,:)*(pe3_in(ji,jj  ,:) - e3u_0(ji,jj  ,:))) )   &
+                     &        + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
+                     &     - SUM(pe3_out(ji,jj,1:ikfm1)))
+               END DO
+            END DO
+            !
+         CASE ( 2 ) 
+            zssh(:,:) = SUM(umask(:,:,:)*(pe3_in(:,:,:)-e3u_0(:,:,:)), DIM=3)
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) =  (  umask(ji,jj,jk)* umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )   &
+                        &                 * 0.5_wp * r1_e1e2f(ji,jj)                                           &
+                        &                 * (e1e2u(ji  ,jj) * zssh(ji  ,jj) + e1e2u(ji,jj+1) * zssh(ji,jj+1))  &
+                        &                 * e3f_0(ji,jj,jk) / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) )
+                  END DO
+               END DO
+            END DO      
+            !
+         END SELECT
          !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
@@ -997,5 +1092,5 @@
                !
             ELSE
-               !
+              !
                ! usr_def_istate called here only to get sshb, that is needed to
                ! initialize e3t_b and e3t_n 
@@ -1061,5 +1156,6 @@
       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
+         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg            , &! not yet implemented: ln_vvl_kepe
+         &              nn_vvl_interp
       !!---------------------------------------------------------------------- 
       !
@@ -1097,4 +1193,5 @@
          ENDIF
          WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
+         WRITE(numout,*) '         Method to compute scale factors anomaly at U/V/F points  nn_vvl_interp   = ', nn_vvl_interp
       ENDIF
       !
@@ -1108,4 +1205,5 @@
       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
+      IF( .NOT. ln_vvl_zstar .AND. (nn_vvl_interp==2 ) )  CALL ctl_stop( 'nn_vvl_interp must be < 2 if ln_vvl_zstar=F' )
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
