Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 7970)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 7971)
@@ -137,5 +137,5 @@
          DO jk=1,jpkm1
             DO jj=1,jpj
-               spgu(2,jj)=spgu(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
+               spgu(2,jj)=spgu(2,jj)+fse3u_a(2,jj,jk)*ua(2,jj,jk)
             END DO
          END DO
@@ -143,5 +143,5 @@
          DO jj=1,jpj
             IF (umask(2,jj,1).NE.0.) THEN
-               spgu(2,jj)=spgu(2,jj)/hu(2,jj)
+               spgu(2,jj)=spgu(2,jj)*hur_a(2,jj)
             ENDIF
          END DO
@@ -161,5 +161,5 @@
          DO jk=1,jpkm1
             DO jj=1,jpj
-               spgu1(2,jj)=spgu1(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
+               spgu1(2,jj)=spgu1(2,jj)+fse3u_a(2,jj,jk)*ua(2,jj,jk)
             END DO
          END DO
@@ -167,5 +167,5 @@
          DO jj=1,jpj
             IF (umask(2,jj,1).NE.0.) THEN
-               spgu1(2,jj)=spgu1(2,jj)/hu(2,jj)
+               spgu1(2,jj)=spgu1(2,jj)*hur_a(2,jj)
             ENDIF
          END DO
@@ -207,10 +207,10 @@
          DO jk=1,jpkm1
             DO jj=1,jpj
-               spgu(nlci-2,jj)=spgu(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)
+               spgu(nlci-2,jj)=spgu(nlci-2,jj)+fse3u_a(nlci-2,jj,jk)*ua(nlci-2,jj,jk)
             ENDDO
          ENDDO
          DO jj=1,jpj
             IF (umask(nlci-2,jj,1).NE.0.) THEN
-               spgu(nlci-2,jj)=spgu(nlci-2,jj)/hu(nlci-2,jj)
+               spgu(nlci-2,jj)=spgu(nlci-2,jj)*hur_a(nlci-2,jj)
             ENDIF
          END DO
@@ -229,10 +229,10 @@
          DO jk=1,jpkm1
             DO jj=1,jpj
-               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)*umask(nlci-2,jj,jk)
+               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+fse3u_a(nlci-2,jj,jk)*ua(nlci-2,jj,jk)*umask(nlci-2,jj,jk)
             END DO
          END DO
          DO jj=1,jpj
             IF (umask(nlci-2,jj,1).NE.0.) THEN
-               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)/hu(nlci-2,jj)
+               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)*hur_a(nlci-2,jj)
             ENDIF
          END DO
@@ -278,5 +278,5 @@
          DO jk=1,jpkm1
             DO ji=1,jpi
-               spgv(ji,2)=spgv(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)
+               spgv(ji,2)=spgv(ji,2)+fse3v_a(ji,2,jk)*va(ji,2,jk)
             END DO
          END DO
@@ -284,5 +284,5 @@
          DO ji=1,jpi
             IF (vmask(ji,2,1).NE.0.) THEN
-               spgv(ji,2)=spgv(ji,2)/hv(ji,2)
+               spgv(ji,2)=spgv(ji,2)*hvr_a(ji,2)
             ENDIF
          END DO
@@ -302,5 +302,5 @@
          DO jk=1,jpkm1
             DO ji=1,jpi
-               spgv1(ji,2)=spgv1(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)*vmask(ji,2,jk)
+               spgv1(ji,2)=spgv1(ji,2)+fse3v_a(ji,2,jk)*va(ji,2,jk)*vmask(ji,2,jk)
             END DO
          END DO
@@ -308,5 +308,5 @@
          DO ji=1,jpi
             IF (vmask(ji,2,1).NE.0.) THEN
-               spgv1(ji,2)=spgv1(ji,2)/hv(ji,2)
+               spgv1(ji,2)=spgv1(ji,2)*hvr_a(ji,2)
             ENDIF
          END DO
@@ -353,5 +353,5 @@
          DO jk=1,jpkm1
             DO ji=1,jpi
-               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
+               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+fse3v_a(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
             END DO
          END DO
@@ -359,5 +359,5 @@
          DO ji=1,jpi
             IF (vmask(ji,nlcj-2,1).NE.0.) THEN
-               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)/hv(ji,nlcj-2)
+               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)*hvr_a(ji,nlcj-2)
             ENDIF
          END DO
@@ -378,5 +378,5 @@
          DO jk=1,jpkm1
             DO ji=1,jpi
-               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
+               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+fse3v_a(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
             END DO
          END DO
@@ -384,5 +384,5 @@
          DO ji=1,jpi
             IF (vmask(ji,nlcj-2,1).NE.0.) THEN
-               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)/hv(ji,nlcj-2)
+               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)*hvr_a(ji,nlcj-2)
             ENDIF
          END DO
@@ -503,10 +503,4 @@
          zt = REAL(Agrif_NbStepint(),wp) / zrhot
       ENDIF
-
-      ! Linear interpolation of sea level
-      Agrif_SpecialValue    = 0.e0
-      Agrif_UseSpecialValue = .TRUE.
-      CALL Agrif_Bc_variable(sshn_id,calledweight=zt, procname=interpsshn )
-      Agrif_UseSpecialValue = .FALSE.
 
       ! Interpolate barotropic fluxes
@@ -539,30 +533,41 @@
    SUBROUTINE Agrif_ssh( kt )
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_DYN  ***
+      !!                  ***  ROUTINE Agrif_ssh  ***
       !!----------------------------------------------------------------------  
       INTEGER, INTENT(in) ::   kt
       !!
+      INTEGER :: ji, jj
       !!----------------------------------------------------------------------  
 
       IF( Agrif_Root() )   RETURN
 
+      ! Linear interpolation of sea level
+      Agrif_SpecialValue    = 0.e0
+      Agrif_UseSpecialValue = .TRUE.
+      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
+      Agrif_UseSpecialValue = .FALSE.
+
       IF((nbondi == -1).OR.(nbondi == 2)) THEN
-         ssha(2,:)=ssha(3,:)
-         sshn(2,:)=sshn(3,:)
+         DO jj=1,jpj
+            ssha(2,jj) = hbdy_w(jj)
+         END DO
       ENDIF
 
       IF((nbondi == 1).OR.(nbondi == 2)) THEN
-         ssha(nlci-1,:)=ssha(nlci-2,:)
-         sshn(nlci-1,:)=sshn(nlci-2,:)
+         DO jj=1,jpj
+            ssha(nlci-1,jj) = hbdy_e(jj)
+         END DO
       ENDIF
 
       IF((nbondj == -1).OR.(nbondj == 2)) THEN
-         ssha(:,2)=ssha(:,3)
-         sshn(:,2)=sshn(:,3)
+         DO ji=1,jpi
+            ssha(ji,2) = hbdy_s(ji)
+         END DO
       ENDIF
 
       IF((nbondj == 1).OR.(nbondj == 2)) THEN
-         ssha(:,nlcj-1)=ssha(:,nlcj-2)
-         sshn(:,nlcj-1)=sshn(:,nlcj-2)
+         DO ji=1,jpi
+            ssha(ji,nlcj-1) = hbdy_n(ji)
+         END DO
       ENDIF
 
@@ -812,5 +817,5 @@
                DO ji=i1,i2
                   ptab(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk)
-                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u(ji,jj,jk)
+                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u_n(ji,jj,jk)
                END DO
             END DO
@@ -821,5 +826,5 @@
             DO jj=j1,j2
                ua(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj)))
-               ua(i1:i2,jj,jk) = ua(i1:i2,jj,jk) / fse3u(i1:i2,jj,jk)
+               ua(i1:i2,jj,jk) = ua(i1:i2,jj,jk) / fse3u_a(i1:i2,jj,jk)
             END DO
          END DO
@@ -880,5 +885,5 @@
                DO ji=i1,i2
                   ptab(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk)
-                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v(ji,jj,jk)
+                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v_n(ji,jj,jk)
                END DO
             END DO
@@ -889,5 +894,5 @@
             DO jj=j1,j2
                va(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj)))
-               va(i1:i2,jj,jk) = va(i1:i2,jj,jk) / fse3v(i1:i2,jj,jk)
+               va(i1:i2,jj,jk) = va(i1:i2,jj,jk) / fse3v_a(i1:i2,jj,jk)
             END DO
          END DO
@@ -1110,5 +1115,5 @@
          ! Polynomial interpolation coefficients:
          zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
-               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
+                 &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
          ! 
          IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
@@ -1151,5 +1156,5 @@
          ! Polynomial interpolation coefficients:
          zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
-               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
+                 &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
          !
          IF(western_side ) vbdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 7970)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 7971)
@@ -13,9 +13,11 @@
    USE dynspg_oce
    USE zdf_oce        ! vertical physics: ocean variables 
+   USE domvvl         ! Need interpolation routines 
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn,Update_Scales
+   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl
+
 # if defined key_zdftke
    PUBLIC Agrif_Update_Tke
@@ -177,9 +179,100 @@
 # endif /* key_zdftke */
 
+   RECURSIVE SUBROUTINE Agrif_Update_vvl( )
+      !!---------------------------------------------
+      !!   *** ROUTINE Agrif_Update_vvl ***
+      !!---------------------------------------------
+      !
+      IF (Agrif_Root()) RETURN
+      !
+#if defined TWO_WAY  
+      !
+      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
+      !
+      Agrif_UseSpecialValueInUpdate = .FALSE.
+      Agrif_SpecialValueFineGrid = 0.
+      ! 
+# if ! defined DECAL_FEEDBACK
+      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)
+# else
+      CALL Agrif_Update_Variable(e3t_id, locupdate=(/1,0/), procname=updatee3t)
+# endif 
+      !
+      Agrif_UseSpecialValueInUpdate = .FALSE.
+      !
+      CALL Agrif_ChildGrid_To_ParentGrid()
+      CALL dom_vvl_update_UVF
+      CALL Agrif_ParentGrid_To_ChildGrid()
+      !
+      IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:
+         CALL Agrif_ChildGrid_To_ParentGrid()
+         CALL Agrif_Update_vvl()
+         CALL Agrif_ParentGrid_To_ChildGrid()
+      ENDIF
+      !
+#endif
+      !
+   END SUBROUTINE Agrif_Update_vvl
+
+   SUBROUTINE dom_vvl_update_UVF
+      !!---------------------------------------------
+      !!       *** ROUTINE dom_vvl_update_UVF ***
+      !!---------------------------------------------
+#  include "domzgr_substitute.h90"
+      !!
+      INTEGER :: jk
+      REAL(wp):: zcoef
+      !!---------------------------------------------
+
+      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', &
+                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step()
+
+      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( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+
+      ! Update total depths:
+      hu(:,:) = 0._wp                        ! Ocean depth at U-points
+      hv(:,:) = 0._wp                        ! Ocean depth at V-points
+      DO jk = 1, jpkm1
+         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+      END DO
+      !                                      ! Inverse of the local depth
+      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+
+#if ! defined key_dynspg_ts
+      IF  (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+#else
+      IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN
+#endif
+         CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:),  'U'  )
+         CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:),  'V'  )
+
+         CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
+         CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
+
+         hu_b(:,:) = 0._wp                        ! Ocean depth at U-points
+         hv_b(:,:) = 0._wp                        ! Ocean depth at V-points
+         DO jk = 1, jpkm1
+            hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+         END DO
+         !                                      ! Inverse of the local depth
+         hur_b(:,:) = 1._wp / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+         hvr_b(:,:) = 1._wp / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+      ENDIF
+
+      !
+   END SUBROUTINE dom_vvl_update_UVF
+
    SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!---------------------------------------------
       !!           *** ROUTINE updateT ***
       !!---------------------------------------------
-#  include "domzgr_substitute.h90"
       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
@@ -188,4 +281,5 @@
       INTEGER :: ji,jj,jk,jn
       !!---------------------------------------------
+      !
       !
       IF (before) THEN
@@ -194,5 +288,8 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
+!> jc tmp
+                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * fse3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
+!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * fse3t_n(ji,jj,jk)
+!< jc tmp
                   END DO
                END DO
@@ -200,4 +297,10 @@
          END DO
       ELSE
+!> jc tmp
+         DO jn = n1,n2
+            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2)
+         ENDDO
+!< jc tmp
+
          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
             ! Add asselin part
@@ -209,5 +312,6 @@
                            tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
                                  & + atfp * ( tabres(ji,jj,jk,jn) &
-                                 &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                                 &             - tsn(ji,jj,jk,jn)*fse3t_a(ji,jj,jk) ) & 
+                                 & * tmask(ji,jj,jk) / fse3t_b(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -221,5 +325,9 @@
                   DO ji=i1,i2
                      IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk)
+#if defined key_vvl
+                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / fse3t_n(ji,jj,jk)
+#else
+                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn)
+#endif
                      END IF
                   END DO
@@ -260,12 +368,14 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e2u(ji,jj) / fse3u_n(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e2u(ji,jj) 
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                      ub(ji,jj,jk) = ub(ji,jj,jk) & 
-                           & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
+                           & + atfp * ( tabres(ji,jj,jk)                     & 
+                           &              - un(ji,jj,jk)*fse3u_a(ji,jj,jk) ) & 
+                           & * umask(ji,jj,jk) / fse3u_b(ji,jj,jk)
                   ENDIF
                   !
-                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk)
+                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / fse3u_n(ji,jj,jk)
                END DO
             END DO
@@ -304,12 +414,14 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e1v(ji,jj) / fse3v_n(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e1v(ji,jj)
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                      vb(ji,jj,jk) = vb(ji,jj,jk) & 
-                           & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
+                           & + atfp * ( tabres(ji,jj,jk)                     & 
+                           &              - vn(ji,jj,jk)*fse3v_a(ji,jj,jk) ) & 
+                           & * vmask(ji,jj,jk) / fse3v_b(ji,jj,jk)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk)
+                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / fse3v_n(ji,jj,jk)
                END DO
             END DO
@@ -345,5 +457,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) =  tabres(ji,jj) * hur(ji,jj) / e2u(ji,jj)  
+               tabres(ji,jj) =  tabres(ji,jj) / e2u(ji,jj)  
                !    
                ! Update "now" 3d velocities:
@@ -352,7 +464,6 @@
                   spgu(ji,jj) = spgu(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk)
                END DO
-               spgu(ji,jj) = spgu(ji,jj) * hur(ji,jj)
                !
-               zcorr = tabres(ji,jj) - spgu(ji,jj)
+               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * hur(ji,jj)
                DO jk=1,jpkm1              
                   un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
@@ -360,11 +471,13 @@
                !
                ! Update barotropic velocities:
-#if defined key_dynspg_ts
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                  zcorr = tabres(ji,jj) - un_b(ji,jj)
+#if ! defined key_dynspg_ts
+               IF  (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+#else
+               IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN
+#endif
+                  zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * hur_b(ji,jj)
                   ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
-               END IF
-#endif               
-               un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1)
+               END IF             
+               un_b(ji,jj) = tabres(ji,jj) * hur(ji,jj) * umask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
@@ -373,7 +486,6 @@
                   spgu(ji,jj) = spgu(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)
                END DO
-               spgu(ji,jj) = spgu(ji,jj) * hur_b(ji,jj)
                !
-               zcorr = ub_b(ji,jj) - spgu(ji,jj)
+               zcorr = ub_b(ji,jj) - spgu(ji,jj) * hur_b(ji,jj)
                DO jk=1,jpkm1              
                   ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
@@ -410,5 +522,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) =  tabres(ji,jj) * hvr(ji,jj) / e1v(ji,jj)  
+               tabres(ji,jj) =  tabres(ji,jj) / e1v(ji,jj)  
                !    
                ! Update "now" 3d velocities:
@@ -417,7 +529,6 @@
                   spgv(ji,jj) = spgv(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk)
                END DO
-               spgv(ji,jj) = spgv(ji,jj) * hvr(ji,jj)
                !
-               zcorr = tabres(ji,jj) - spgv(ji,jj)
+               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * hvr(ji,jj)
                DO jk=1,jpkm1              
                   vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
@@ -425,11 +536,13 @@
                !
                ! Update barotropic velocities:
-#if defined key_dynspg_ts
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                  zcorr = tabres(ji,jj) - vn_b(ji,jj)
+#if ! defined key_dynspg_ts
+               IF  (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+#else
+               IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN
+#endif
+                  zcorr = (tabres(ji,jj) - vn_b(ji,jj)*hv_a(ji,jj)) * hvr_b(ji,jj)
                   vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
-               END IF
-#endif               
-               vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1)
+               END IF              
+               vn_b(ji,jj) = tabres(ji,jj) * hvr(ji,jj) * vmask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
@@ -438,7 +551,6 @@
                   spgv(ji,jj) = spgv(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
                END DO
-               spgv(ji,jj) = spgv(ji,jj) * hvr_b(ji,jj)
                !
-               zcorr = vb_b(ji,jj) - spgv(ji,jj)
+               zcorr = vb_b(ji,jj) - spgv(ji,jj) * hvr_b(ji,jj)
                DO jk=1,jpkm1              
                   vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
@@ -471,5 +583,8 @@
       ELSE
 #if ! defined key_dynspg_ts
-         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+         IF  (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+#else
+         IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN
+#endif
             DO jj=j1,j2
                DO ji=i1,i2
@@ -479,5 +594,5 @@
             END DO
          ENDIF
-#endif
+         !
          DO jj=j1,j2
             DO ji=i1,i2
@@ -655,4 +770,99 @@
 # endif /* key_zdftke */ 
 
+   SUBROUTINE updatee3t( ptab, i1, i2, j1, j2, k1, k2, before )
+      !!---------------------------------------------
+      !!           *** ROUTINE updatee3t ***
+      !!---------------------------------------------
+      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
+      LOGICAL, INTENT(in) :: before
+      INTEGER :: ji,jj,jk
+      REAL(wp) :: zcoef
+      !!---------------------------------------------
+      !
+      IF (before) THEN
+!> jc tmp:
+!         ptab(i1:i2,j1:j2,k1:k2) = fse3t_n(i1:i2,j1:j2,k1:k2)
+         ptab(i1:i2,j1:j2,k1:k2) = fse3t_n(i1:i2,j1:j2,k1:k2) / e3t_0(i1:i2,j1:j2,k1:k2)
+!< jc tmp:
+      ELSE
+         !
+         ! 1) Updates at before time step:
+         ! -------------------------------
+         !
+!> jc tmp:
+         ptab(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
+!< jc tmp:
+
+#if ! defined key_dynspg_ts
+         IF  (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+#else
+         IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN
+#endif
+            DO jk = 1, jpkm1
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     fse3t_b(ji,jj,jk) =   fse3t_b(ji,jj,jk) &
+                           & + atfp * ( ptab(ji,jj,jk) - fse3t_n(ji,jj,jk) )
+                  END DO
+               END DO
+            END DO
+            !
+            fse3w_b (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + fse3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
+            fsdepw_b(i1:i2,j1:j2,1) = 0.0_wp
+            fsdept_b(i1:i2,j1:j2,1) = 0.5_wp * fse3w_b(i1:i2,j1:j2,1)
+            !
+            DO jk = 2, jpkm1
+               DO jj = j1,j2
+                  DO ji = i1,i2            
+                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+                     fse3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
+                     &                                        ( fse3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
+                     &                                    +            0.5_wp * tmask(ji,jj,jk)   *        &
+                     &                                        ( fse3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
+                     fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
+                     fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  &
+                         &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk)) 
+                  END DO
+               END DO
+            END DO
+         ENDIF        
+         !
+         ! 2) Updates at now time step:
+         ! ----------------------------
+         !
+         ! Update vertical scale factor at T-points:
+         fse3t_n(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2)
+         !
+         ! Update total depth:
+         ht(i1:i2,j1:j2) = 0._wp
+         DO jk = 1, jpkm1
+            ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + fse3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)
+         END DO
+         !
+         ! Update vertical scale factor at W-points and depths:
+         fse3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + fse3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
+         fsdept_n(i1:i2,j1:j2,1) = 0.5_wp * fse3w_n(i1:i2,j1:j2,1)
+         fsdepw_n(i1:i2,j1:j2,1) = 0.0_wp
+         fsde3w_n(i1:i2,j1:j2,1) = fsdept_n(i1:i2,j1:j2,1) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
+         !
+         DO jk = 2, jpkm1
+            DO jj = j1,j2
+               DO ji = i1,i2            
+               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+               fse3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( fse3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   &
+               &                                    +            0.5_wp * tmask(ji,jj,jk)   * ( fse3t_n(ji,jj,jk  ) - e3t_0(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) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE updatee3t
+
 #else
 CONTAINS
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 7970)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 7971)
@@ -309,5 +309,6 @@
 !     or the absolute maximum nesting level...TBC                        
    IF ( Agrif_Level().EQ.Agrif_MaxLevel() ) THEN 
-      ! NB: Do tracers first, dynamics after because nbcline incremented in dynamics
+      ! NB: Order matters below: 
+      CALL Agrif_Update_vvl() 
       CALL Agrif_Update_tra()
       CALL Agrif_Update_dyn()
@@ -445,4 +446,6 @@
    CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Average)
 # endif
+
+   CALL Agrif_Set_Updatetype(e3t_id, update = AGRIF_Update_Average)
 
 ! High order updates
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 7970)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 7971)
@@ -1006,5 +1006,5 @@
 
 #if defined key_agrif
-      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
+      IF ((.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar)) CALL ctl_stop( 'AGRIF implemented with zstar coordinate only (key_vvl)' )
 #endif
 
Index: /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7970)
+++ /branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7971)
@@ -350,4 +350,5 @@
 
       IF ( Agrif_NbStepint().EQ.0 ) THEN
+         IF( lk_vvl )          CALL Agrif_Update_vvl()      ! Update vert. grid 
                                CALL Agrif_Update_Tra()      ! Update active tracers
                                CALL Agrif_Update_Dyn()      ! Update momentum
