Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/EXP00/1_namelist_cfg
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/EXP00/1_namelist_cfg	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/EXP00/1_namelist_cfg	(revision 8741)
@@ -6,4 +6,8 @@
 &namusr_def    !   User defined :   VORTEX configuration: Flat bottom, beta-plane
 !-----------------------------------------------------------------------
+   rn_dx       =  30000.   !  x horizontal resolution   [meters]
+   rn_dy       =  30000.   !  y horizontal resolution   [meters]
+   rn_dz       =    500.   !  z vertical resolution [meters]
+   rn_ppgphi0  =   38.5    !  Reference latitude [degrees]
 /
 !
@@ -37,5 +41,6 @@
    rn_sponge_tra = 2850.   !  coefficient for tracer   sponge layer [m2/s]
    rn_sponge_dyn = 2850.   !  coefficient for dynamics sponge layer [m2/s]
-   ln_chk_bathy  = .false. !
+   ln_chk_bathy  = .FALSE. !
+/
 /
 !-----------------------------------------------------------------------
@@ -160,4 +165,5 @@
 &namdyn_spg    !   Surface pressure gradient
 !-----------------------------------------------------------------------
+   ln_dynspg_exp  = .false. 
    ln_dynspg_ts  = .true.   ! split-explicit free surface
       ln_bt_fw      = .true.     ! Forward integration of barotropic Eqs.
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/EXP00/namelist_cfg	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/EXP00/namelist_cfg	(revision 8741)
@@ -37,4 +37,10 @@
 &namagrif      !  AGRIF zoom                                            ("key_agrif")
 !-----------------------------------------------------------------------
+   nn_cln_update =    1    !  baroclinic update frequency
+   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
+   rn_sponge_tra = 2850.   !  coefficient for tracer   sponge layer [m2/s]
+   rn_sponge_dyn = 2850.   !  coefficient for dynamics sponge layer [m2/s]
+   ln_chk_bathy  = .FALSE. !
+/
 /
 !-----------------------------------------------------------------------
@@ -159,4 +165,5 @@
 &namdyn_spg    !   Surface pressure gradient
 !-----------------------------------------------------------------------
+   ln_dynspg_exp  = .false. 
    ln_dynspg_ts  = .true.   ! split-explicit free surface
       ln_bt_fw      = .true.     ! Forward integration of barotropic Eqs.
Index: anches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/agrif_opa_interp.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/agrif_opa_interp.F90	(revision 8740)
+++ 	(revision )
@@ -1,1209 +1,0 @@
-MODULE agrif_opa_interp
-   !!======================================================================
-   !!                   ***  MODULE  agrif_opa_interp  ***
-   !! AGRIF: interpolation package
-   !!======================================================================
-   !! History :  2.0  !  2002-06  (XXX)  Original cade
-   !!             -   !  2005-11  (XXX) 
-   !!            3.2  !  2009-04  (R. Benshila) 
-   !!            3.6  !  2014-09  (R. Benshila) 
-   !!----------------------------------------------------------------------
-#if defined key_agrif
-   !!----------------------------------------------------------------------
-   !!   'key_agrif'                                              AGRIF zoom
-   !!----------------------------------------------------------------------
-   !!   Agrif_tra     :
-   !!   Agrif_dyn     : 
-   !!   interpu       :
-   !!   interpv       :
-   !!----------------------------------------------------------------------
-   USE par_oce
-   USE oce
-   USE dom_oce      
-   USE zdf_oce
-   USE agrif_oce
-   USE phycst
-   !
-   USE in_out_manager
-   USE agrif_opa_sponge
-   USE lib_mpp
-   USE wrk_nemo
- 
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
-   PUBLIC   interpun, interpvn
-   PUBLIC   interptsn,  interpsshn
-   PUBLIC   interpunb, interpvnb, interpub2b, interpvb2b
-   PUBLIC   interpe3t, interpumsk, interpvmsk
-# if defined key_zdftke
-   PUBLIC   Agrif_tke, interpavm
-# endif
-
-   INTEGER ::   bdy_tinterp = 0
-
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/NST 3.7 , NEMO Consortium (2015)
-   !! $Id: agrif_opa_interp.F90 7646 2017-02-06 09:25:03Z timgraham $
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE Agrif_tra
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_tra  ***
-      !!----------------------------------------------------------------------
-      !
-      IF( Agrif_Root() )   RETURN
-      !
-      Agrif_SpecialValue    = 0._wp
-      Agrif_UseSpecialValue = .TRUE.
-      !
-      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
-      !
-      Agrif_UseSpecialValue = .FALSE.
-      !
-   END SUBROUTINE Agrif_tra
-
-
-   SUBROUTINE Agrif_dyn( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_DYN  ***
-      !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) ::   kt
-      !
-      INTEGER ::   ji, jj, jk       ! dummy loop indices
-      INTEGER ::   j1, j2, i1, i2
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zub, zvb
-      !!----------------------------------------------------------------------  
-      !
-      IF( Agrif_Root() )   RETURN
-      !
-      CALL wrk_alloc( jpi,jpj,   zub, zvb )
-      !
-      Agrif_SpecialValue    = 0._wp
-      Agrif_UseSpecialValue = ln_spc_dyn
-      !
-      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
-      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
-      !
-      Agrif_UseSpecialValue = .FALSE.
-      !
-      ! prevent smoothing in ghost cells
-      i1 =  1   ;   i2 = jpi
-      j1 =  1   ;   j2 = jpj
-      IF( nbondj == -1 .OR. nbondj == 2 )   j1 = 3
-      IF( nbondj == +1 .OR. nbondj == 2 )   j2 = nlcj-2
-      IF( nbondi == -1 .OR. nbondi == 2 )   i1 = 3
-      IF( nbondi == +1 .OR. nbondi == 2 )   i2 = nlci-2
-
-      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
-         !
-         ! Smoothing
-         ! ---------
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            ua_b(2,:) = 0._wp
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  ua_b(2,jj) = ua_b(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
-               END DO
-            END DO
-            DO jj = 1, jpj
-               ua_b(2,jj) = ua_b(2,jj) * r1_hu_a(2,jj)            
-            END DO
-         ENDIF
-         !
-         DO jk=1,jpkm1                 ! Smooth
-            DO jj=j1,j2
-               ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk))
-               ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk)
-            END DO
-         END DO
-         !
-         zub(2,:) = 0._wp              ! Correct transport
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               zub(2,jj) = zub(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
-            END DO
-         END DO
-         DO jj=1,jpj
-            zub(2,jj) = zub(2,jj) * r1_hu_a(2,jj)
-         END DO
-
-         DO jk=1,jpkm1
-            DO jj=1,jpj
-               ua(2,jj,jk) = (ua(2,jj,jk)+ua_b(2,jj)-zub(2,jj))*umask(2,jj,jk)
-            END DO
-         END DO
-
-         ! Set tangential velocities to time splitting estimate
-         !-----------------------------------------------------
-         IF( ln_dynspg_ts ) THEN
-            zvb(2,:) = 0._wp
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  zvb(2,jj) = zvb(2,jj) + e3v_a(2,jj,jk) * va(2,jj,jk)
-               END DO
-            END DO
-            DO jj = 1, jpj
-               zvb(2,jj) = zvb(2,jj) * r1_hv_a(2,jj)
-            END DO
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-zvb(2,jj)) * vmask(2,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         ! Mask domain edges:
-         !-------------------
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               ua(1,jj,jk) = 0._wp
-               va(1,jj,jk) = 0._wp
-            END DO
-         END DO         
-         !
-      ENDIF
-
-      IF( nbondi == 1 .OR. nbondi == 2 ) THEN
-
-         ! Smoothing
-         ! ---------
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            ua_b(nlci-2,:) = 0._wp
-            DO jk=1,jpkm1
-               DO jj=1,jpj
-                  ua_b(nlci-2,jj) = ua_b(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
-               END DO
-            END DO
-            DO jj=1,jpj
-               ua_b(nlci-2,jj) = ua_b(nlci-2,jj) * r1_hu_a(nlci-2,jj)            
-            END DO
-         ENDIF
-
-         DO jk = 1, jpkm1              ! Smooth
-            DO jj = j1, j2
-               ua(nlci-2,jj,jk) = 0.25_wp * umask(nlci-2,jj,jk)      &
-                  &             * ( ua(nlci-3,jj,jk) + 2._wp*ua(nlci-2,jj,jk) + ua(nlci-1,jj,jk) )
-            END DO
-         END DO
-
-         zub(nlci-2,:) = 0._wp        ! Correct transport
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               zub(nlci-2,jj) = zub(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
-            END DO
-         END DO
-         DO jj = 1, jpj
-            zub(nlci-2,jj) = zub(nlci-2,jj) * r1_hu_a(nlci-2,jj)
-         END DO
-
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               ua(nlci-2,jj,jk) = ( ua(nlci-2,jj,jk) + ua_b(nlci-2,jj) - zub(nlci-2,jj) ) * umask(nlci-2,jj,jk)
-            END DO
-         END DO
-         !
-         ! Set tangential velocities to time splitting estimate
-         !-----------------------------------------------------
-         IF( ln_dynspg_ts ) THEN
-            zvb(nlci-1,:) = 0._wp
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  zvb(nlci-1,jj) = zvb(nlci-1,jj) + e3v_a(nlci-1,jj,jk) * va(nlci-1,jj,jk)
-               END DO
-            END DO
-            DO jj=1,jpj
-               zvb(nlci-1,jj) = zvb(nlci-1,jj) * r1_hv_a(nlci-1,jj)
-            END DO
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  va(nlci-1,jj,jk) = ( va(nlci-1,jj,jk) + va_b(nlci-1,jj) - zvb(nlci-1,jj) ) * vmask(nlci-1,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         ! Mask domain edges:
-         !-------------------
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               ua(nlci-1,jj,jk) = 0._wp
-               va(nlci  ,jj,jk) = 0._wp
-            END DO
-         END DO 
-         !
-      ENDIF
-
-      IF( nbondj == -1 .OR. nbondj == 2 ) THEN
-
-         ! Smoothing
-         ! ---------
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            va_b(:,2) = 0._wp
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  va_b(ji,2) = va_b(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk)
-               END DO
-            END DO
-            DO ji=1,jpi
-               va_b(ji,2) = va_b(ji,2) * r1_hv_a(ji,2)            
-            END DO
-         ENDIF
-         !
-         DO jk = 1, jpkm1              ! Smooth
-            DO ji = i1, i2
-               va(ji,2,jk) = 0.25_wp * vmask(ji,2,jk)    &
-                  &        * ( va(ji,1,jk) + 2._wp*va(ji,2,jk) + va(ji,3,jk) )
-            END DO
-         END DO
-         !
-         zvb(:,2) = 0._wp              ! Correct transport
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               zvb(ji,2) = zvb(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk)
-            END DO
-         END DO
-         DO ji = 1, jpi
-            zvb(ji,2) = zvb(ji,2) * r1_hv_a(ji,2)
-         END DO
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               va(ji,2,jk) = ( va(ji,2,jk) + va_b(ji,2) - zvb(ji,2) ) * vmask(ji,2,jk)
-            END DO
-         END DO
-
-         ! Set tangential velocities to time splitting estimate
-         !-----------------------------------------------------
-         IF( ln_dynspg_ts ) THEN
-            zub(:,2) = 0._wp
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  zub(ji,2) = zub(ji,2) + e3u_a(ji,2,jk) * ua(ji,2,jk) * umask(ji,2,jk)
-               END DO
-            END DO
-            DO ji = 1, jpi
-               zub(ji,2) = zub(ji,2) * r1_hu_a(ji,2)
-            END DO
-
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  ua(ji,2,jk) = ( ua(ji,2,jk) + ua_b(ji,2) - zub(ji,2) ) * umask(ji,2,jk)
-               END DO
-            END DO
-         ENDIF
-
-         ! Mask domain edges:
-         !-------------------
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               ua(ji,1,jk) = 0._wp
-               va(ji,1,jk) = 0._wp
-            END DO
-         END DO 
-
-      ENDIF
-
-      IF( nbondj == 1 .OR. nbondj == 2 ) THEN
-         !
-         ! Smoothing
-         ! ---------
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            va_b(:,nlcj-2) = 0._wp
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  va_b(ji,nlcj-2) = va_b(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk)
-               END DO
-            END DO
-            DO ji = 1, jpi
-               va_b(ji,nlcj-2) = va_b(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)            
-            END DO
-         ENDIF
-         !
-         DO jk = 1, jpkm1              ! Smooth
-            DO ji = i1, i2
-               va(ji,nlcj-2,jk) = 0.25_wp * vmask(ji,nlcj-2,jk)   &
-                  &             * ( va(ji,nlcj-3,jk) + 2._wp * va(ji,nlcj-2,jk) + va(ji,nlcj-1,jk) )
-            END DO
-         END DO
-         !
-         zvb(:,nlcj-2) = 0._wp         ! Correct transport
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk)
-            END DO
-         END DO
-         DO ji = 1, jpi
-            zvb(ji,nlcj-2) = zvb(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)
-         END DO
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               va(ji,nlcj-2,jk) = ( va(ji,nlcj-2,jk) + va_b(ji,nlcj-2) - zvb(ji,nlcj-2) ) * vmask(ji,nlcj-2,jk)
-            END DO
-         END DO
-         !
-         ! Set tangential velocities to time splitting estimate
-         !-----------------------------------------------------
-         IF( ln_dynspg_ts ) THEN
-            zub(:,nlcj-1) = 0._wp
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  zub(ji,nlcj-1) = zub(ji,nlcj-1) + e3u_a(ji,nlcj-1,jk) * ua(ji,nlcj-1,jk) * umask(ji,nlcj-1,jk)
-               END DO
-            END DO
-            DO ji = 1, jpi
-               zub(ji,nlcj-1) = zub(ji,nlcj-1) * r1_hu_a(ji,nlcj-1)
-            END DO
-            !
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  ua(ji,nlcj-1,jk) = ( ua(ji,nlcj-1,jk) + ua_b(ji,nlcj-1) - zub(ji,nlcj-1) ) * umask(ji,nlcj-1,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         ! Mask domain edges:
-         !-------------------
-         DO jk = 1, jpkm1
-            DO ji = 1, jpi
-               ua(ji,nlcj  ,jk) = 0._wp
-               va(ji,nlcj-1,jk) = 0._wp
-            END DO
-         END DO 
-         !
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj,   zub, zvb )
-      !
-   END SUBROUTINE Agrif_dyn
-
-
-   SUBROUTINE Agrif_dyn_ts( jn )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_dyn_ts  ***
-      !!----------------------------------------------------------------------  
-      !! 
-      INTEGER, INTENT(in) ::   jn
-      !!
-      INTEGER :: ji, jj
-      !!----------------------------------------------------------------------  
-      !
-      IF( Agrif_Root() )   RETURN
-      !
-      IF((nbondi == -1).OR.(nbondi == 2)) THEN
-         DO jj=1,jpj
-            va_e(2,jj) = vbdy_w(jj) * hvr_e(2,jj)
-            ! Specified fluxes:
-            ua_e(2,jj) = ubdy_w(jj) * hur_e(2,jj)
-            ! Characteristics method:
-            !alt            ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) &
-            !alt                       &           - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) )
-         END DO
-      ENDIF
-      !
-      IF((nbondi == 1).OR.(nbondi == 2)) THEN
-         DO jj=1,jpj
-            va_e(nlci-1,jj) = vbdy_e(jj) * hvr_e(nlci-1,jj)
-            ! Specified fluxes:
-            ua_e(nlci-2,jj) = ubdy_e(jj) * hur_e(nlci-2,jj)
-            ! Characteristics method:
-            !alt            ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) &
-            !alt                            &           + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) )
-         END DO
-      ENDIF
-      !
-      IF((nbondj == -1).OR.(nbondj == 2)) THEN
-         DO ji=1,jpi
-            ua_e(ji,2) = ubdy_s(ji) * hur_e(ji,2)
-            ! Specified fluxes:
-            va_e(ji,2) = vbdy_s(ji) * hvr_e(ji,2)
-            ! Characteristics method:
-            !alt            va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) &
-            !alt                       &           - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) )
-         END DO
-      ENDIF
-      !
-      IF((nbondj == 1).OR.(nbondj == 2)) THEN
-         DO ji=1,jpi
-            ua_e(ji,nlcj-1) = ubdy_n(ji) * hur_e(ji,nlcj-1)
-            ! Specified fluxes:
-            va_e(ji,nlcj-2) = vbdy_n(ji) * hvr_e(ji,nlcj-2)
-            ! Characteristics method:
-            !alt            va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2)  + va_e(ji,nlcj-3) &
-            !alt                            &           + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) )
-         END DO
-      ENDIF
-      !
-   END SUBROUTINE Agrif_dyn_ts
-
-
-   SUBROUTINE Agrif_dta_ts( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_dta_ts  ***
-      !!----------------------------------------------------------------------  
-      !! 
-      INTEGER, INTENT(in) ::   kt
-      !!
-      INTEGER :: ji, jj
-      LOGICAL :: ll_int_cons
-      REAL(wp) :: zrhot, zt
-      !!----------------------------------------------------------------------  
-      !
-      IF( Agrif_Root() )   RETURN
-      !
-      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
-      !
-      zrhot = Agrif_rhot()
-      !
-      ! "Central" time index for interpolation:
-      IF( ln_bt_fw ) THEN
-         zt = REAL( Agrif_NbStepint()+0.5_wp, wp ) / zrhot
-      ELSE
-         zt = REAL( Agrif_NbStepint()       , wp ) / zrhot
-      ENDIF
-      !
-      ! Linear interpolation of sea level
-      Agrif_SpecialValue    = 0._wp
-      Agrif_UseSpecialValue = .TRUE.
-      CALL Agrif_Bc_variable( sshn_id, calledweight=zt, procname=interpsshn )
-      Agrif_UseSpecialValue = .FALSE.
-      !
-      ! Interpolate barotropic fluxes
-      Agrif_SpecialValue=0.
-      Agrif_UseSpecialValue = ln_spc_dyn
-      !
-      IF( ll_int_cons ) THEN  ! Conservative interpolation
-         ! orders matters here !!!!!!
-         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
-         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )
-         bdy_tinterp = 1
-         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
-         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )
-         bdy_tinterp = 2
-         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
-         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )         
-      ELSE ! Linear interpolation
-         bdy_tinterp = 0
-         ubdy_w(:) = 0._wp   ;   vbdy_w(:) = 0._wp 
-         ubdy_e(:) = 0._wp   ;   vbdy_e(:) = 0._wp 
-         ubdy_n(:) = 0._wp   ;   vbdy_n(:) = 0._wp 
-         ubdy_s(:) = 0._wp   ;   vbdy_s(:) = 0._wp
-         CALL Agrif_Bc_variable( unb_id, calledweight=zt, procname=interpunb )
-         CALL Agrif_Bc_variable( vnb_id, calledweight=zt, procname=interpvnb )
-      ENDIF
-      Agrif_UseSpecialValue = .FALSE.
-      ! 
-   END SUBROUTINE Agrif_dta_ts
-
-
-   SUBROUTINE Agrif_ssh( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_DYN  ***
-      !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) ::   kt
-      !!
-      !!----------------------------------------------------------------------  
-      !
-      IF( Agrif_Root() )   RETURN
-      !
-      IF((nbondi == -1).OR.(nbondi == 2)) THEN
-         ssha(2,:)=ssha(3,:)
-         sshn(2,:)=sshn(3,:)
-      ENDIF
-      !
-      IF((nbondi == 1).OR.(nbondi == 2)) THEN
-         ssha(nlci-1,:)=ssha(nlci-2,:)
-         sshn(nlci-1,:)=sshn(nlci-2,:)
-      ENDIF
-      !
-      IF((nbondj == -1).OR.(nbondj == 2)) THEN
-         ssha(:,2)=ssha(:,3)
-         sshn(:,2)=sshn(:,3)
-      ENDIF
-      !
-      IF((nbondj == 1).OR.(nbondj == 2)) THEN
-         ssha(:,nlcj-1)=ssha(:,nlcj-2)
-         sshn(:,nlcj-1)=sshn(:,nlcj-2)
-      ENDIF
-      !
-   END SUBROUTINE Agrif_ssh
-
-
-   SUBROUTINE Agrif_ssh_ts( jn )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_ssh_ts  ***
-      !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) ::   jn
-      !!
-      INTEGER :: ji,jj
-      !!----------------------------------------------------------------------  
-      !
-      IF((nbondi == -1).OR.(nbondi == 2)) THEN
-         DO jj = 1, jpj
-            ssha_e(2,jj) = hbdy_w(jj)
-         END DO
-      ENDIF
-      !
-      IF((nbondi == 1).OR.(nbondi == 2)) THEN
-         DO jj = 1, jpj
-            ssha_e(nlci-1,jj) = hbdy_e(jj)
-         END DO
-      ENDIF
-      !
-      IF((nbondj == -1).OR.(nbondj == 2)) THEN
-         DO ji = 1, jpi
-            ssha_e(ji,2) = hbdy_s(ji)
-         END DO
-      ENDIF
-      !
-      IF((nbondj == 1).OR.(nbondj == 2)) THEN
-         DO ji = 1, jpi
-            ssha_e(ji,nlcj-1) = hbdy_n(ji)
-         END DO
-      ENDIF
-      !
-   END SUBROUTINE Agrif_ssh_ts
-
-# if defined key_zdftke
-
-   SUBROUTINE Agrif_tke
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_tke  ***
-      !!----------------------------------------------------------------------  
-      REAL(wp) ::   zalpha
-      !!----------------------------------------------------------------------  
-      !
-      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
-      IF( zalpha > 1. )   zalpha = 1.
-      !
-      Agrif_SpecialValue    = 0.e0
-      Agrif_UseSpecialValue = .TRUE.
-      !
-      CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
-      !
-      Agrif_UseSpecialValue = .FALSE.
-      !
-   END SUBROUTINE Agrif_tke
-   
-# endif
-
-   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!   *** ROUTINE interptsn ***
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
-      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
-      LOGICAL                                     , INTENT(in   ) ::   before
-      INTEGER                                     , INTENT(in   ) ::   nb , ndir
-      !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::   imin, imax, jmin, jmax
-      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
-      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------
-      !
-      IF (before) THEN         
-         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
-      ELSE
-         !
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         !
-         zrhox = Agrif_Rhox()
-         ! 
-         zalpha1 = ( zrhox - 1. ) * 0.5
-         zalpha2 = 1. - zalpha1
-         ! 
-         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
-         zalpha4 = 1. - zalpha3
-         ! 
-         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
-         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
-         zalpha5 = 1. - zalpha6 - zalpha7
-         !
-         imin = i1
-         imax = i2
-         jmin = j1
-         jmax = j2
-         ! 
-         ! Remove CORNERS
-         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
-         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
-         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
-         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2        
-         !
-         IF( eastern_side ) THEN
-            DO jn = 1, jpts
-               tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)
-               DO jk = 1, jpkm1
-                  DO jj = jmin,jmax
-                     IF( umask(nlci-2,jj,jk) == 0._wp ) THEN
-                        tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
-                     ELSE
-                        tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
-                        IF( un(nlci-2,jj,jk) > 0._wp ) THEN
-                           tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
-                                 + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
-                        ENDIF
-                     ENDIF
-                  END DO
-               END DO
-               tsa(nlci,j1:j2,k1:k2,jn) = 0._wp
-            END DO
-         ENDIF
-         ! 
-         IF( northern_side ) THEN            
-            DO jn = 1, jpts
-               tsa(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)
-               DO jk = 1, jpkm1
-                  DO ji = imin,imax
-                     IF( vmask(ji,nlcj-2,jk) == 0._wp ) THEN
-                        tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
-                     ELSE
-                        tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)        
-                        IF (vn(ji,nlcj-2,jk) > 0._wp ) THEN
-                           tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
-                                 + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
-                        ENDIF
-                     ENDIF
-                  END DO
-               END DO
-               tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp
-            END DO
-         ENDIF
-         !
-         IF( western_side ) THEN            
-            DO jn = 1, jpts
-               tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)
-               DO jk = 1, jpkm1
-                  DO jj = jmin,jmax
-                     IF( umask(2,jj,jk) == 0._wp ) THEN
-                        tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
-                     ELSE
-                        tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)        
-                        IF( un(2,jj,jk) < 0._wp ) THEN
-                           tsa(2,jj,jk,jn)=(zalpha6*tsa(3,jj,jk,jn)+zalpha5*tsa(1,jj,jk,jn)+zalpha7*tsa(4,jj,jk,jn))*tmask(2,jj,jk)
-                        ENDIF
-                     ENDIF
-                  END DO
-               END DO
-               tsa(1,j1:j2,k1:k2,jn) = 0._wp
-            END DO
-         ENDIF
-         !
-         IF( southern_side ) THEN           
-            DO jn = 1, jpts
-               tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)
-               DO jk = 1, jpk      
-                  DO ji=imin,imax
-                     IF( vmask(ji,2,jk) == 0._wp ) THEN
-                        tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
-                     ELSE
-                        tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
-                        IF( vn(ji,2,jk) < 0._wp ) THEN
-                           tsa(ji,2,jk,jn)=(zalpha6*tsa(ji,3,jk,jn)+zalpha5*tsa(ji,1,jk,jn)+zalpha7*tsa(ji,4,jk,jn))*tmask(ji,2,jk)
-                        ENDIF
-                     ENDIF
-                  END DO
-               END DO
-               tsa(i1:i2,1,k1:k2,jn) = 0._wp
-            END DO
-         ENDIF
-         !
-         ! Treatment of corners
-         ! 
-         ! East south
-         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:)
-         ENDIF
-         ! East north
-         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:)
-         ENDIF
-         ! West south
-         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tsa(2,2,:,:) = ptab(2,2,:,:)
-         ENDIF
-         ! West north
-         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:)
-         ENDIF
-         !
-      ENDIF
-      !
-   END SUBROUTINE interptsn
-
-
-   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpsshn  ***
-      !!----------------------------------------------------------------------  
-      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
-      LOGICAL                         , INTENT(in   ) ::   before
-      INTEGER                         , INTENT(in   ) ::   nb , ndir
-      !
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------  
-      !
-      IF( before) THEN
-         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
-      ELSE
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         IF(western_side)  hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
-         IF(eastern_side)  hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
-         IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
-         IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
-      ENDIF
-      !
-   END SUBROUTINE interpsshn
-
-
-   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before )
-      !!----------------------------------------------------------------------
-      !!   *** ROUTINE interpun ***
-      !!----------------------------------------------------------------------
-      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) ::   zrhoy	
-      !!----------------------------------------------------------------------
-      !
-      IF( before ) THEN 
-         DO jk = k1, jpk
-            ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
-         END DO
-      ELSE
-         zrhoy = Agrif_Rhoy()
-         DO jk = 1, jpkm1
-            DO jj=j1,j2
-               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) )
-            END DO
-         END DO
-      ENDIF
-      ! 
-   END SUBROUTINE interpun
-
-
-   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before )
-      !!----------------------------------------------------------------------
-      !!   *** ROUTINE interpvn ***
-      !!----------------------------------------------------------------------
-      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) ::   zrhox	
-      !!----------------------------------------------------------------------
-      !      
-      IF( before ) THEN       !interpv entre 1 et k2 et interpv2d en jpkp1
-         DO jk = k1, jpk
-            ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk)
-         END DO
-      ELSE          
-         zrhox= Agrif_Rhox()
-         DO jk = 1, jpkm1
-            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) )
-         END DO
-      ENDIF
-      !        
-   END SUBROUTINE interpvn
-   
-
-   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpunb  ***
-      !!----------------------------------------------------------------------  
-      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
-      LOGICAL                         , INTENT(in   ) ::   before
-      INTEGER                         , INTENT(in   ) ::   nb , ndir
-      !
-      INTEGER  ::   ji, jj
-      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------  
-      !
-      IF( before ) THEN 
-         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
-      ELSE
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         zrhoy = Agrif_Rhoy()
-         zrhot = Agrif_rhot()
-         ! Time indexes bounds for integration
-         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
-         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot      
-         ! Polynomial interpolation coefficients:
-         IF( bdy_tinterp == 1 ) THEN
-            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-               &               - zt0**2._wp * (       zt0 - 1._wp)        )
-         ELSEIF( bdy_tinterp == 2 ) THEN
-            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
-               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
-
-         ELSE
-            ztcoeff = 1
-         ENDIF
-         !   
-         IF(western_side) THEN
-            ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2)  
-         ENDIF
-         IF(eastern_side) THEN
-            ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2)  
-         ENDIF
-         IF(southern_side) THEN
-            ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
-         ENDIF
-         IF(northern_side) THEN
-            ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
-         ENDIF
-         !            
-         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
-            IF(western_side) THEN
-               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
-            ENDIF
-            IF(eastern_side) THEN
-               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
-            ENDIF
-            IF(southern_side) THEN
-               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
-            ENDIF
-            IF(northern_side) THEN
-               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
-            ENDIF
-         ENDIF
-      ENDIF
-      ! 
-   END SUBROUTINE interpunb
-
-
-   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpvnb  ***
-      !!----------------------------------------------------------------------  
-      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
-      LOGICAL                         , INTENT(in   ) ::   before
-      INTEGER                         , INTENT(in   ) ::   nb , ndir
-      !
-      INTEGER  ::   ji,jj
-      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------  
-      ! 
-      IF( before ) THEN 
-         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
-      ELSE
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         zrhox = Agrif_Rhox()
-         zrhot = Agrif_rhot()
-         ! Time indexes bounds for integration
-         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
-         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot      
-         IF( bdy_tinterp == 1 ) THEN
-            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-               &               - zt0**2._wp * (       zt0 - 1._wp)        )
-         ELSEIF( bdy_tinterp == 2 ) THEN
-            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
-               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
-         ELSE
-            ztcoeff = 1
-         ENDIF
-         !
-         IF(western_side) THEN
-            vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2)  
-         ENDIF
-         IF(eastern_side) THEN
-            vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2)  
-         ENDIF
-         IF(southern_side) THEN
-            vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
-         ENDIF
-         IF(northern_side) THEN
-            vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
-         ENDIF
-         !            
-         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
-            IF(western_side) THEN
-               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
-                     &                                  * vmask(i1,j1:j2,1)
-            ENDIF
-            IF(eastern_side) THEN
-               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
-                     &                                  * vmask(i1,j1:j2,1)
-            ENDIF
-            IF(southern_side) THEN
-               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
-                     &                                  * vmask(i1:i2,j1,1)
-            ENDIF
-            IF(northern_side) THEN
-               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
-                     &                                  * vmask(i1:i2,j1,1)
-            ENDIF
-         ENDIF
-      ENDIF
-      !
-   END SUBROUTINE interpvnb
-
-
-   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpub2b  ***
-      !!----------------------------------------------------------------------  
-      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
-      LOGICAL                         , INTENT(in   ) ::   before
-      INTEGER                         , INTENT(in   ) ::   nb , ndir
-      !
-      INTEGER  ::   ji,jj
-      REAL(wp) ::   zrhot, zt0, zt1,zat
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------  
-      IF( before ) THEN
-         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
-      ELSE
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         zrhot = Agrif_rhot()
-         ! Time indexes bounds for integration
-         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
-         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
-         ! Polynomial interpolation coefficients:
-         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
-            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
-         ! 
-         IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
-         IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
-      ENDIF
-      ! 
-   END SUBROUTINE interpub2b
-   
-
-   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpvb2b  ***
-      !!----------------------------------------------------------------------  
-      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
-      LOGICAL                         , INTENT(in   ) ::   before
-      INTEGER                         , INTENT(in   ) ::   nb , ndir
-      !
-      INTEGER ::   ji,jj
-      REAL(wp) ::   zrhot, zt0, zt1,zat
-      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------  
-      !
-      IF( before ) THEN
-         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
-      ELSE      
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         zrhot = Agrif_rhot()
-         ! Time indexes bounds for integration
-         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
-         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
-         ! Polynomial interpolation coefficients:
-         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
-            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
-         !
-         IF(western_side )   vbdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(eastern_side )   vbdy_e(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(southern_side)   vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
-         IF(northern_side)   vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
-      ENDIF
-      !      
-   END SUBROUTINE interpvb2b
-
-
-   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpe3t  ***
-      !!----------------------------------------------------------------------  
-      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                              , INTENT(in   ) :: nb , ndir
-      !
-      INTEGER :: ji, jj, jk
-      LOGICAL :: western_side, eastern_side, northern_side, southern_side
-      REAL(wp) :: ztmpmsk      
-      !!----------------------------------------------------------------------  
-      !    
-      IF( before ) THEN
-         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
-      ELSE
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-
-         DO jk = k1, k2
-            DO jj = j1, j2
-               DO ji = i1, i2
-                  ! Get velocity mask at boundary edge points:
-                  IF( western_side )   ztmpmsk = umask(ji    ,jj    ,1)
-                  IF( eastern_side )   ztmpmsk = umask(nlci-2,jj    ,1)
-                  IF( northern_side)   ztmpmsk = vmask(ji    ,nlcj-2,1)
-                  IF( southern_side)   ztmpmsk = vmask(ji    ,2     ,1)
-                  !
-                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) )*ztmpmsk > 1.D-2) THEN
-                     IF (western_side) THEN
-                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
-                     ELSEIF (eastern_side) THEN
-                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
-                     ELSEIF (southern_side) THEN
-                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
-                     ELSEIF (northern_side) THEN
-                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
-                     ENDIF
-                     WRITE(numout,*) '      ptab(ji,jj,jk), e3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
-                     kindic_agr = kindic_agr + 1
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE interpe3t
-
-
-   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpumsk  ***
-      !!----------------------------------------------------------------------  
-      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                              , INTENT(in   ) ::   nb , ndir
-      !
-      INTEGER ::   ji, jj, jk
-      LOGICAL ::   western_side, eastern_side   
-      !!----------------------------------------------------------------------  
-      !    
-      IF( before ) THEN
-         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
-      ELSE
-         western_side = (nb == 1).AND.(ndir == 1)
-         eastern_side = (nb == 1).AND.(ndir == 2)
-         DO jk = k1, k2
-            DO jj = j1, j2
-               DO ji = i1, i2
-                   ! Velocity mask at boundary edge points:
-                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
-                     IF (western_side) THEN
-                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
-                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
-                        kindic_agr = kindic_agr + 1
-                     ELSEIF (eastern_side) THEN
-                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
-                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
-                        kindic_agr = kindic_agr + 1
-                     ENDIF
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE interpumsk
-
-
-   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interpvmsk  ***
-      !!----------------------------------------------------------------------  
-      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                              , INTENT(in   ) :: nb , ndir
-      !
-      INTEGER ::   ji, jj, jk
-      LOGICAL ::   northern_side, southern_side     
-      !!----------------------------------------------------------------------  
-      !    
-      IF( before ) THEN
-         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
-      ELSE
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         DO jk = k1, k2
-            DO jj = j1, j2
-               DO ji = i1, i2
-                   ! Velocity mask at boundary edge points:
-                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
-                     IF (southern_side) THEN
-                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
-                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
-                        kindic_agr = kindic_agr + 1
-                     ELSEIF (northern_side) THEN
-                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
-                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
-                        kindic_agr = kindic_agr + 1
-                     ENDIF
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE interpvmsk
-
-# if defined key_zdftke
-
-   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE interavm  ***
-      !!----------------------------------------------------------------------  
-      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
-      !!----------------------------------------------------------------------  
-      !      
-      IF( before ) THEN
-         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
-      ELSE
-         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
-      ENDIF
-      !
-   END SUBROUTINE interpavm
-
-# endif /* key_zdftke */
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Empty module                                          no AGRIF zoom
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE Agrif_OPA_Interp_empty
-      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
-   END SUBROUTINE Agrif_OPA_Interp_empty
-#endif
-
-   !!======================================================================
-END MODULE agrif_opa_interp
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/domvvl.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/domvvl.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/domvvl.F90	(revision 8741)
@@ -1034,5 +1034,5 @@
       !
 #if defined key_agrif
-      IF(.NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' )
+      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
 #endif
       !
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8741)
@@ -35,5 +35,4 @@
    LOGICAL , PUBLIC :: spongedoneU = .FALSE.       !: dynamics sponge layer indicator
    LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE.     !: if true: first step
-   LOGICAL , PUBLIC :: lk_agrif_doupd = .TRUE.     !: if true: send update from current grid
    LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE.    !: if true: print debugging info
 
@@ -65,5 +64,5 @@
    INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id
    INTEGER :: scales_t_id
-# if defined key_zdftke
+# if defined key_zdftke || defined key_zdfgls
    INTEGER :: avt_id, avm_id, en_id
 # endif  
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8741)
@@ -24,4 +24,5 @@
    USE agrif_oce
    USE phycst
+   USE dynspg_ts, ONLY: un_adv, vn_adv
    !
    USE in_out_manager
@@ -38,6 +39,6 @@
    PUBLIC   interpunb, interpvnb, interpub2b, interpvb2b
    PUBLIC   interpe3t, interpumsk, interpvmsk
-# if defined key_zdftke
-   PUBLIC   Agrif_tke, interpavm
+# if defined key_zdftke || defined key_zdfgls
+   PUBLIC   Agrif_avm, interpavm
 # endif
 
@@ -449,5 +450,4 @@
       INTEGER :: ji, jj
       LOGICAL :: ll_int_cons
-      REAL(wp) :: zrhot, zt
       !!----------------------------------------------------------------------  
       !
@@ -456,21 +456,9 @@
       ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
       !
-      zrhot = Agrif_rhot()
-      !
-      ! "Central" time index for interpolation:
-      IF( ln_bt_fw ) THEN
-         zt = REAL( Agrif_NbStepint()+0.5_wp, wp ) / zrhot
-      ELSE
-         zt = REAL( Agrif_NbStepint()       , wp ) / zrhot
-      ENDIF
-      !
-      ! Linear interpolation of sea level
-      Agrif_SpecialValue    = 0._wp
-      Agrif_UseSpecialValue = .TRUE.
-      CALL Agrif_Bc_variable( sshn_id, calledweight=zt, procname=interpsshn )
-      Agrif_UseSpecialValue = .FALSE.
+      ! Enforce volume conservation if no time refinement:  
+      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE.  
       !
       ! Interpolate barotropic fluxes
-      Agrif_SpecialValue=0.
+      Agrif_SpecialValue=0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
       !
@@ -491,6 +479,6 @@
          ubdy_n(:) = 0._wp   ;   vbdy_n(:) = 0._wp 
          ubdy_s(:) = 0._wp   ;   vbdy_s(:) = 0._wp
-         CALL Agrif_Bc_variable( unb_id, calledweight=zt, procname=interpunb )
-         CALL Agrif_Bc_variable( vnb_id, calledweight=zt, procname=interpvnb )
+         CALL Agrif_Bc_variable( unb_id, procname=interpunb )
+         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb )
       ENDIF
       Agrif_UseSpecialValue = .FALSE.
@@ -501,30 +489,42 @@
    SUBROUTINE Agrif_ssh( kt )
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_DYN  ***
+      !!                  ***  ROUTINE Agrif_ssh  ***
       !!----------------------------------------------------------------------  
       INTEGER, INTENT(in) ::   kt
       !!
+      INTEGER :: ji, jj
       !!----------------------------------------------------------------------  
       !
       IF( Agrif_Root() )   RETURN
+      !      
+      ! Linear interpolation in time of sea level
+      !
+      Agrif_SpecialValue    = 0._wp
+      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
       !
@@ -541,4 +541,7 @@
       !!----------------------------------------------------------------------  
       !
+      !
+      IF( Agrif_Root() )   RETURN
+      !
       IF((nbondi == -1).OR.(nbondi == 2)) THEN
          DO jj = 1, jpj
@@ -567,15 +570,19 @@
    END SUBROUTINE Agrif_ssh_ts
 
-# if defined key_zdftke
-
-   SUBROUTINE Agrif_tke
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_tke  ***
+# if defined key_zdftke || defined key_zdfgls
+
+   SUBROUTINE Agrif_avm
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE Agrif_avm  ***
       !!----------------------------------------------------------------------  
       REAL(wp) ::   zalpha
       !!----------------------------------------------------------------------  
       !
-      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
-      IF( zalpha > 1. )   zalpha = 1.
+      IF( Agrif_Root() )   RETURN
+      !
+!      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
+!      IF( zalpha > 1. )   zalpha = 1.
+      zalpha = 1._wp ! JC: proper time interpolation impossible  
+                     ! => use last available value from parent 
       !
       Agrif_SpecialValue    = 0.e0
@@ -586,5 +593,5 @@
       Agrif_UseSpecialValue = .FALSE.
       !
-   END SUBROUTINE Agrif_tke
+   END SUBROUTINE Agrif_avm
    
 # endif
@@ -781,5 +788,5 @@
       !
       IF( before ) THEN 
-         DO jk = k1, jpk
+         DO jk = 1, jpkm1
             ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
          END DO
@@ -788,5 +795,5 @@
          DO jk = 1, jpkm1
             DO jj=j1,j2
-               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) )
+               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
             END DO
          END DO
@@ -808,6 +815,6 @@
       !!----------------------------------------------------------------------
       !      
-      IF( before ) THEN       !interpv entre 1 et k2 et interpv2d en jpkp1
-         DO jk = k1, jpk
+      IF( before ) THEN  
+         DO jk = 1, jpkm1
             ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk)
          END DO
@@ -815,5 +822,5 @@
          zrhox= Agrif_Rhox()
          DO jk = 1, jpkm1
-            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) )
+            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
          END DO
       ENDIF
@@ -978,5 +985,9 @@
       !!----------------------------------------------------------------------  
       IF( before ) THEN
-         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
+         IF ( ln_bt_fw ) THEN
+            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
+         ELSE
+            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
+         ENDIF
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1016,5 +1027,9 @@
       !
       IF( before ) THEN
-         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
+         IF ( ln_bt_fw ) THEN
+            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
+         ELSE
+            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
+         ENDIF
       ELSE      
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1175,5 +1190,5 @@
    END SUBROUTINE interpvmsk
 
-# if defined key_zdftke
+# if defined key_zdftke || defined key_zdfgls
 
    SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
@@ -1189,10 +1204,10 @@
          ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
       ELSE
-         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
+         avm  (i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
       ENDIF
       !
    END SUBROUTINE interpavm
 
-# endif /* key_zdftke */
+# endif /* key_zdftke || key_zdfgls */
 
 #else
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8741)
@@ -12,9 +12,11 @@
    USE wrk_nemo  
    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
@@ -27,5 +29,5 @@
 CONTAINS
 
-   RECURSIVE SUBROUTINE Agrif_Update_Tra( )
+   SUBROUTINE Agrif_Update_Tra( )
       !!---------------------------------------------
       !!   *** ROUTINE Agrif_Update_Tra ***
@@ -56,16 +58,9 @@
       Agrif_UseSpecialValueInUpdate = .FALSE.
       !
-      IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:
-         CALL Agrif_ChildGrid_To_ParentGrid()
-         CALL Agrif_Update_Tra()
-         CALL Agrif_ParentGrid_To_ChildGrid()
-      ENDIF
-      !
 #endif
       !
    END SUBROUTINE Agrif_Update_Tra
 
-
-   RECURSIVE SUBROUTINE Agrif_Update_Dyn( )
+   SUBROUTINE Agrif_Update_Dyn( )
       !!---------------------------------------------
       !!   *** ROUTINE Agrif_Update_Dyn ***
@@ -140,21 +135,15 @@
 #endif
       !
-      ! Do recursive update:
-      IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:
-         CALL Agrif_ChildGrid_To_ParentGrid()
-         CALL Agrif_Update_Dyn()
-         CALL Agrif_ParentGrid_To_ChildGrid()
-      ENDIF
-      !
    END SUBROUTINE Agrif_Update_Dyn
 
 # if defined key_zdftke
 
-   SUBROUTINE Agrif_Update_Tke( kt )
+   SUBROUTINE Agrif_Update_Tke( )
       !!---------------------------------------------
       !!   *** ROUTINE Agrif_Update_Tke ***
       !!---------------------------------------------
       !!
-      INTEGER, INTENT(in) :: kt
+      ! 
+      IF (Agrif_Root()) RETURN
       !       
       IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
@@ -176,5 +165,112 @@
 # endif /* key_zdftke */
 
-   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
+   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 = .TRUE.
+      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()
+      !
+#endif
+      !
+   END SUBROUTINE Agrif_Update_vvl
+
+   SUBROUTINE dom_vvl_update_UVF
+      !!---------------------------------------------
+      !!       *** ROUTINE dom_vvl_update_UVF ***
+      !!---------------------------------------------
+      !!
+      INTEGER :: jk
+      REAL(wp):: zcoef
+      !!---------------------------------------------
+
+      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', &
+                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step()
+
+      ! Save "old" scale factor (prior update) for subsequent asselin correction
+      ! of prognostic variables (needed to update initial state only)
+      ! -------------------------------------------------------------
+      !
+      e3u_a(:,:,:) = e3u_n(:,:,:)
+      e3v_a(:,:,:) = e3v_n(:,:,:)
+!      ua(:,:,:) = e3u_b(:,:,:)
+!      va(:,:,:) = e3v_b(:,:,:)
+      hu_a(:,:) = hu_n(:,:)
+      hv_a(:,:) = hv_n(:,:)
+
+      ! 1) NOW fields
+      !--------------
+      
+         ! Vertical scale factor interpolations
+         ! ------------------------------------
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' )
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' )
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' )
+
+      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+
+         ! Update total depths:
+         ! --------------------
+      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points
+      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
+      DO jk = 1, jpkm1
+         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
+         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+      END DO
+      !                                        ! Inverse of the local depth
+      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
+      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+
+
+      ! 2) BEFORE fields:
+      !------------------
+      IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
+         & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
+         & .AND.(.NOT.ln_bt_fw)))) THEN
+         !
+         ! Vertical scale factor interpolations
+         ! ------------------------------------
+         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  )
+         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  )
+
+         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
+         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+
+         ! Update total depths:
+         ! --------------------
+         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(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
+            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+         END DO
+         !                                     ! Inverse of the local depth
+         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
+         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
+      ENDIF
+      !
+   END SUBROUTINE dom_vvl_update_UVF
+
+   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
       !!---------------------------------------------
       !!           *** ROUTINE updateT ***
@@ -183,6 +279,9 @@
       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
       LOGICAL, INTENT(in) :: before
+      INTEGER, INTENT(in) :: nb, ndir
       !!
+      LOGICAL :: western_side, eastern_side, southern_side, northern_side 
       INTEGER :: ji,jj,jk,jn
+      REAL(wp) :: ztb, ztnu, ztno
       !!---------------------------------------------
       !
@@ -192,5 +291,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)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
+!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
+!< jc tmp
                   END DO
                END DO
@@ -198,4 +300,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) &
+                                         & * tmask(i1:i2,j1:j2,k1:k2)
+         ENDDO
+!< jc tmp
          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
             ! Add asselin part
@@ -205,7 +313,9 @@
                      DO ji=i1,i2
                         IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
-                           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)
+                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
+                           ztnu = tabres(ji,jj,jk,jn)
+                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
+                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -219,5 +329,5 @@
                   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)
+                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
                      END IF
                   END DO
@@ -225,4 +335,111 @@
             END DO
          END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            tsb(i1:i2,j1:j2,k1:k2,n1:n2)  = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
+         ENDIF
+         !
+         !
+# if defined DECAL_FEEDBACK
+         IF (.NOT.ln_linssh) THEN 
+            western_side  = (nb == 1).AND.(ndir == 1)
+            eastern_side  = (nb == 1).AND.(ndir == 2)
+            southern_side = (nb == 2).AND.(ndir == 1)
+            northern_side = (nb == 2).AND.(ndir == 2)
+            !
+            ! Asselin correction 
+            IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+               IF (southern_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO ji=i1,i2
+                           ztb  = tsb(ji,j1-1,jk,jn) * e3t_b(ji,j1-1,jk) ! fse3t_b prior update should be used
+                           ztnu = tsn(ji,j1-1,jk,jn) * e3t_n(ji,j1-1,jk)
+                           ztno = tsn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk)
+                           tsb(ji,j1-1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,j1-1,jk) / e3t_b(ji,j1-1,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+               IF (northern_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO ji=i1,i2
+                           ztb  = tsb(ji,j2+1,jk,jn) * e3t_b(ji,j2+1,jk) ! fse3t_b prior update should be used
+                           ztnu = tsn(ji,j2+1,jk,jn) * e3t_n(ji,j2+1,jk)
+                           ztno = tsn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk)
+                           tsb(ji,j2+1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,j2+1,jk) / e3t_b(ji,j2+1,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+               IF (western_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO jj=j1,j2
+                           ztb  = tsb(i1-1,jj,jk,jn) * e3t_b(i1-1,jj,jk) ! fse3t_b prior update should be used
+                           ztnu = tsn(i1-1,jj,jk,jn) * e3t_n(i1-1,jj,jk)
+                           ztno = tsn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk)
+                           tsb(i1-1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(i1-1,jj,jk) / e3t_b(i1-1,jj,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+               IF (eastern_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO jj=j1,j2
+                           ztb  = tsb(i2+1,jj,jk,jn) * e3t_b(i2+1,jj,jk) ! fse3t_b prior update should be used
+                           ztnu = tsn(i2+1,jj,jk,jn) * e3t_n(i2+1,jj,jk)
+                           ztno = tsn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk)
+                           tsb(i2+1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(i2+1,jj,jk) / e3t_b(i2+1,jj,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+            ENDIF ! Asselin correction
+
+            IF (southern_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO ji=i1,i2
+                        tsn(ji,j1-1,jk,jn) = tsn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk) / e3t_n(ji,j1-1,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+            IF (northern_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO ji=i1,i2
+                        tsn(ji,j2+1,jk,jn) = tsn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk) / e3t_n(ji,j2+1,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+            IF (western_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO jj=j1,j2
+                        tsn(i1-1,jj,jk,jn) = tsn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk) / e3t_n(i1-1,jj,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+            IF (eastern_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO jj=j1,j2
+                        tsn(i2+1,jj,jk,jn) = tsn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk) / e3t_n(i2+1,jj,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+         ENDIF
+#endif
       ENDIF
       ! 
@@ -238,6 +455,6 @@
       LOGICAL                               , INTENT(in   ) :: before
       !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhoy
+      INTEGER  :: ji, jj, jk
+      REAL(wp) :: zrhoy, zub, zunu, zuno
       !!---------------------------------------------
       ! 
@@ -251,15 +468,23 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_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)
+                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)
+                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
+                     zunu = tabres(ji,jj,jk)
+                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &      
+                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
                   ENDIF
                   !
-                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
+         ENDIF
+         !
       ENDIF
       ! 
@@ -267,14 +492,14 @@
 
 
-   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before )
+   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before)
       !!---------------------------------------------
       !!           *** ROUTINE updatev ***
       !!---------------------------------------------
-      INTEGER :: i1,i2,j1,j2,k1,k2
-      INTEGER :: ji,jj,jk
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres
-      LOGICAL :: before
-      !!
-      REAL(wp) :: zrhox
+      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
+      LOGICAL                               , INTENT(in   ) :: before
+      !
+      INTEGER  :: ji, jj, jk
+      REAL(wp) :: zrhox, zvb, zvnu, zvno
       !!---------------------------------------------      
       !
@@ -292,15 +517,23 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_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)
+                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk)
+                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
+                     zvnu = tabres(ji,jj,jk)
+                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &      
+                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
+         ENDIF
+         !
       ENDIF
       ! 
@@ -316,5 +549,5 @@
       LOGICAL, INTENT(in) :: before
       !! 
-      INTEGER :: ji, jj, jk
+      INTEGER  :: ji, jj, jk
       REAL(wp) :: zrhoy
       REAL(wp) :: zcorr
@@ -331,5 +564,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) =  tabres(ji,jj) * r1_hu_n(ji,jj) * r1_e2u(ji,jj)  
+               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)  
                !    
                ! Update "now" 3d velocities:
@@ -338,7 +571,6 @@
                   spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
                END DO
-               spgu(ji,jj) = spgu(ji,jj) * r1_hu_n(ji,jj)
                !
-               zcorr = tabres(ji,jj) - spgu(ji,jj)
+               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
                DO jk=1,jpkm1              
                   un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
@@ -348,9 +580,9 @@
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zcorr = tabres(ji,jj) - un_b(ji,jj)
+                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_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)
+               ENDIF    
+               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
@@ -359,7 +591,6 @@
                   spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
                END DO
-               spgu(ji,jj) = spgu(ji,jj) * r1_hu_b(ji,jj)
                !
-               zcorr = ub_b(ji,jj) - spgu(ji,jj)
+               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
                DO jk=1,jpkm1              
                   ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
@@ -368,4 +599,8 @@
             END DO
          END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
+         ENDIF
       ENDIF
       !
@@ -381,5 +616,5 @@
       LOGICAL, INTENT(in) :: before
       !! 
-      INTEGER :: ji, jj, jk
+      INTEGER  :: ji, jj, jk
       REAL(wp) :: zrhox
       REAL(wp) :: zcorr
@@ -396,5 +631,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) =  tabres(ji,jj) * r1_hv_n(ji,jj) * r1_e1v(ji,jj)  
+               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)  
                !    
                ! Update "now" 3d velocities:
@@ -403,7 +638,6 @@
                   spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
                END DO
-               spgv(ji,jj) = spgv(ji,jj) * r1_hv_n(ji,jj)
                !
-               zcorr = tabres(ji,jj) - spgv(ji,jj)
+               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
                DO jk=1,jpkm1              
                   vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
@@ -413,9 +647,9 @@
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zcorr = tabres(ji,jj) - vn_b(ji,jj)
+                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_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)
+               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
@@ -424,7 +658,6 @@
                   spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
                END DO
-               spgv(ji,jj) = spgv(ji,jj) * r1_hv_b(ji,jj)
                !
-               zcorr = vb_b(ji,jj) - spgv(ji,jj)
+               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
                DO jk=1,jpkm1              
                   vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
@@ -433,4 +666,9 @@
             END DO
          END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
+         ENDIF
+         !
       ENDIF
       ! 
@@ -438,5 +676,5 @@
 
 
-   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
+   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before, nb, ndir )
       !!---------------------------------------------
       !!          *** ROUTINE updateSSH ***
@@ -445,5 +683,7 @@
       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
       LOGICAL, INTENT(in) :: before
+      INTEGER, INTENT(in) :: nb, ndir
       !!
+      LOGICAL :: western_side, eastern_side, southern_side, northern_side 
       INTEGER :: ji, jj
       !!---------------------------------------------
@@ -472,4 +712,40 @@
             END DO
          END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
+         ENDIF
+         !
+# if defined DECAL_FEEDBACK
+!         western_side  = (nb == 1).AND.(ndir == 1)
+!         eastern_side  = (nb == 1).AND.(ndir == 2)
+!         southern_side = (nb == 2).AND.(ndir == 1)
+!         northern_side = (nb == 2).AND.(ndir == 2)
+!         !
+!         ! Asselin correction 
+!         IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
+!            IF (southern_side) THEN
+!               DO ji=i1,i2
+!                  sshn(ji,j1-1) = sshn(ji,j1-1) - rdt * r1_e2t(ji,j1-1) * (vb2_b_s(ji,j1-1)-vb2_b(ji,j1-1))
+!               END DO
+!            ENDIF
+!            IF (northern_side) THEN
+!               DO ji=i1,i2
+!                  sshn(ji,j1+1) = sshn(ji,j1+1) + rdt * r1_e2t(ji,j1+1) * (vb2_b_s(ji,j1)-vb2_b(ji,j1))
+!               END DO
+!            ENDIF
+!            IF (western_side) THEN
+!               DO jj=j1,j2
+!                  sshn(i1-1,jj) = sshn(i1-1,jj) - rdt * r1_e2t(i1-1,jj) * (ub2_b_s(i1-1,jj)-ub2_b(i1-1,jj))
+!               END DO
+!            ENDIF
+!            IF (eastern_side) THEN
+!               DO jj=j1,j2
+!                  sshn(i1+1,jj) = sshn(i1+1,jj) + rdt * r1_e2t(i1+1,jj) * (ub2_b_s(i1,jj)-ub2_b(i1,jj))
+!               END DO
+!            ENDIF
+!            ! 
+!         ENDIF
+#endif
       ENDIF
       !
@@ -486,5 +762,5 @@
       !!
       INTEGER :: ji, jj
-      REAL(wp) :: zrhoy
+      REAL(wp) :: zrhoy, za1
       !!---------------------------------------------
       !
@@ -498,7 +774,12 @@
          tabres = zrhoy * tabres
       ELSE
+         za1 = 1._wp / REAL(Agrif_rhot(), wp)
+         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
          DO jj=j1,j2
             DO ji=i1,i2
-               ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj)
+               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) & 
+                & + za1 * (tabres(ji,jj) - ub2_b(ji,jj))
+!               ub2_b_s(ji,jj) = ub2_b(ji,jj)
+               ub2_b(ji,jj) = tabres(ji,jj)
             END DO
          END DO
@@ -517,5 +798,5 @@
       !!
       INTEGER :: ji, jj
-      REAL(wp) :: zrhox
+      REAL(wp) :: zrhox, za1
       !!---------------------------------------------
       !
@@ -529,7 +810,12 @@
          tabres = zrhox * tabres
       ELSE
+         za1 = 1._wp / REAL(Agrif_rhot(), wp)
+         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
          DO jj=j1,j2
             DO ji=i1,i2
-               vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj)
+               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) & 
+                & + za1 * (tabres(ji,jj) - vb2_b(ji,jj))
+!               vb2_b_s(ji,jj) = vb2_b(ji,jj)
+               vb2_b(ji,jj) = tabres(ji,jj)
             END DO
          END DO
@@ -644,4 +930,122 @@
 # 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) = e3t_n(i1:i2,j1:j2,k1:k2)
+         ptab(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2) / e3t_0(i1:i2,j1:j2,k1:k2) * tmask(i1:i2,j1:j2,k1:k2)
+!< jc tmp:
+      ELSE
+         !
+         ! 1) Updates at BEFORE time step:
+         ! -------------------------------
+         !
+!> jc tmp:
+!         DO jk = 1, jpkm1
+!            DO jj=j1,j2
+!               DO ji=i1,i2
+!                  IF (tmask(ji,jj,jk)==1) THEN
+!                     ptab(ji,jj,jk) = ptab(ji,jj,jk) * e3t_0(ji,jj,jk)
+!                  ELSE
+!                     ptab(ji,jj,jk) = e3t_0(ji,jj,jk)
+!                  ENDIF
+!               END DO
+!            END DO
+!         END DO
+         ptab(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
+!< jc tmp:
+
+         ! Save "old" scale factor (prior update) for subsequent asselin correction
+         ! of prognostic variables (needed to update initial state only)
+         e3t_a(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2)
+!         hdivb(i1:i2,j1:j2,k1:k2)   = e3t_b(i1:i2,j1:j2,k1:k2)
+
+         IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
+            & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
+            & .AND.(.NOT.ln_bt_fw)))) THEN
+
+            DO jk = 1, jpkm1
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
+                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
+                  END DO
+               END DO
+            END DO
+            !
+            e3w_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
+            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
+            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
+            !
+            DO jk = 2, jpk
+               DO jj = j1,j2
+                  DO ji = i1,i2            
+                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
+                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
+                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
+                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
+                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
+                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
+                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF        
+         !
+         ! 2) Updates at NOW time step:
+         ! ----------------------------
+         !
+         ! Update vertical scale factor at T-points:
+         e3t_n(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2)
+         !
+         ! Update total depth:
+         ht_n(i1:i2,j1:j2) = 0._wp
+         DO jk = 1, jpkm1
+            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)
+         END DO
+         !
+         ! Update vertical scale factor at W-points and depths:
+         e3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
+         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
+         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
+         gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
+         !
+         DO jk = 2, jpk
+            DO jj = j1,j2
+               DO ji = i1,i2            
+               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+               e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   &
+               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
+               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
+               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
+                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
+               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
+               END DO
+            END DO
+         END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
+            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
+            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
+            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE updatee3t
+
 #else
 CONTAINS
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90	(revision 8741)
@@ -33,12 +33,11 @@
 CONTAINS
 
-   SUBROUTINE Agrif_Update_Trc( kt )
+   SUBROUTINE Agrif_Update_Trc( )
       !!----------------------------------------------------------------------
       !!                   *** ROUTINE Agrif_Update_Trc ***
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt
-      !!----------------------------------------------------------------------
       ! 
-      IF((Agrif_NbStepint() .NE. (Agrif_irhot()-1)).AND.(kt /= 0)) RETURN
+      IF (Agrif_Root()) RETURN 
+      !
 #if defined TWO_WAY   
       Agrif_UseSpecialValueInUpdate = .TRUE.
@@ -66,5 +65,5 @@
 
 
-   SUBROUTINE updateTRC( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
+   SUBROUTINE updateTRC( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                      *** ROUTINE updateT ***
@@ -73,33 +72,56 @@
       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
       LOGICAL                                    , INTENT(in   ) ::   before
+      INTEGER, INTENT(in) :: nb, ndir
       !!
-      INTEGER ::   ji, jj, jk, jn
-      !!----------------------------------------------------------------------
-      !
-      IF( before ) THEN
-         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2)
-      ELSE
-         IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN
-            ! Add asselin part
-            DO jn = n1,n2
-               DO jk = k1, k2
-                  DO jj = j1, j2
-                     DO ji = i1, i2
-                        IF( ptab(ji,jj,jk,jn) /= 0._wp ) THEN
-                           trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn)             & 
-                              &             + atfp * ( ptab(ji,jj,jk,jn)   &
-                                 &                    - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
-                        ENDIF
-                     END DO
+      LOGICAL :: western_side, eastern_side, southern_side, northern_side 
+      INTEGER :: ji,jj,jk,jn
+      REAL(wp) :: ztb, ztnu, ztno
+      !!----------------------------------------------------------------------
+      !
+      !
+      IF (before) THEN
+         DO jn = n1,n2
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+!> jc tmp
+                     tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
+!                     tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
+!< jc tmp
                   END DO
                END DO
             END DO
+         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) &
+                                         & * tmask(i1:i2,j1:j2,k1:k2)
+         ENDDO
+!< jc tmp
+         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+            ! Add asselin part
+            DO jn = n1,n2
+               DO jk=k1,k2
+                  DO jj=j1,j2
+                     DO ji=i1,i2
+                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
+                           ztb  = trb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
+                           ztnu = tabres(ji,jj,jk,jn)
+                           ztno = trn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
+                           trb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
+                        ENDIF
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDDO
          ENDIF
-         DO jn = n1, n2
-            DO jk = k1, k2
-               DO jj = j1, j2
-                  DO ji = i1, i2
-                     IF( ptab(ji,jj,jk,jn) /= 0._wp ) THEN 
-                        trn(ji,jj,jk,jn) = ptab(ji,jj,jk,jn) * tmask(ji,jj,jk)
+         DO jn = n1,n2
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
+                        trn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
                      END IF
                   END DO
@@ -107,4 +129,111 @@
             END DO
          END DO
+         !
+         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
+            trb(i1:i2,j1:j2,k1:k2,n1:n2)  = trn(i1:i2,j1:j2,k1:k2,n1:n2)
+         ENDIF
+         !
+         !
+# if defined DECAL_FEEDBACK
+         IF (.NOT.ln_linssh) THEN 
+            western_side  = (nb == 1).AND.(ndir == 1)
+            eastern_side  = (nb == 1).AND.(ndir == 2)
+            southern_side = (nb == 2).AND.(ndir == 1)
+            northern_side = (nb == 2).AND.(ndir == 2)
+            !
+            ! Asselin correction 
+            IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+               IF (southern_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO ji=i1,i2
+                           ztb  = trb(ji,j1-1,jk,jn) * e3t_b(ji,j1-1,jk) ! fse3t_b prior update should be used
+                           ztnu = trn(ji,j1-1,jk,jn) * e3t_n(ji,j1-1,jk)
+                           ztno = trn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk)
+                           trb(ji,j1-1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,j1-1,jk) / e3t_b(ji,j1-1,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+               IF (northern_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO ji=i1,i2
+                           ztb  = trb(ji,j2+1,jk,jn) * e3t_b(ji,j2+1,jk) ! fse3t_b prior update should be used
+                           ztnu = trn(ji,j2+1,jk,jn) * e3t_n(ji,j2+1,jk)
+                           ztno = trn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk)
+                           trb(ji,j2+1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,j2+1,jk) / e3t_b(ji,j2+1,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+               IF (western_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO jj=j1,j2
+                           ztb  = trb(i1-1,jj,jk,jn) * e3t_b(i1-1,jj,jk) ! fse3t_b prior update should be used
+                           ztnu = trn(i1-1,jj,jk,jn) * e3t_n(i1-1,jj,jk)
+                           ztno = trn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk)
+                           trb(i1-1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(i1-1,jj,jk) / e3t_b(i1-1,jj,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+               IF (eastern_side) THEN
+                  DO jn = n1,n2
+                     DO jk=k1,k2
+                        DO jj=j1,j2
+                           ztb  = trb(i2+1,jj,jk,jn) * e3t_b(i2+1,jj,jk) ! fse3t_b prior update should be used
+                           ztnu = trn(i2+1,jj,jk,jn) * e3t_n(i2+1,jj,jk)
+                           ztno = trn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk)
+                           trb(i2+1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(i2+1,jj,jk) / e3t_b(i2+1,jj,jk)
+                        END DO
+                     ENDDO
+                  ENDDO
+               ENDIF
+            ENDIF ! Asselin correction
+
+            IF (southern_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO ji=i1,i2
+                        trn(ji,j1-1,jk,jn) = trn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk) / e3t_n(ji,j1-1,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+            IF (northern_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO ji=i1,i2
+                        trn(ji,j2+1,jk,jn) = trn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk) / e3t_n(ji,j2+1,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+            IF (western_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO jj=j1,j2
+                        trn(i1-1,jj,jk,jn) = trn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk) / e3t_n(i1-1,jj,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+            IF (eastern_side) THEN
+               DO jn = n1,n2
+                  DO jk=k1,k2
+                     DO jj=j1,j2
+                        trn(i2+1,jj,jk,jn) = trn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk) / e3t_n(i2+1,jj,jk)
+                     END DO
+                  ENDDO
+               ENDDO
+            ENDIF
+         ENDIF
+#endif
       ENDIF
       ! 
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 8741)
@@ -1,2 +1,3 @@
+#undef UPD_HIGH   /* MIX HIGH UPDATE */
 #if defined key_agrif
 !!----------------------------------------------------------------------
@@ -88,4 +89,13 @@
 # endif
    !
+   nbcline     = 0
+#if defined key_top
+   nbcline_trc = 0
+#endif
+   !
+   IF ( Agrif_Level().EQ.Agrif_MaxLevel() ) CALL agrif_Update_ini()
+
+   Agrif_UseSpecialValueInUpdate = .FALSE.     
+
 END SUBROUTINE Agrif_initvalues
 
@@ -144,13 +154,14 @@
    CALL Agrif_Set_bc(e2v_id,(/0,0/))
 
-   ! 5. Update type
+   ! 4. Update type
    !--------------- 
+# if defined UPD_HIGH
+   CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting)
+   CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average)
+#else
    CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average)
    CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy)
-
-! High order updates
-!   CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average,            update2=Agrif_Update_Full_Weighting)
-!   CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting,     update2=Agrif_Update_Average)
-    !
+#endif
+
 END SUBROUTINE agrif_declare_var_dom
 
@@ -175,5 +186,5 @@
    !
    LOGICAL :: check_namelist
-   CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3
+   CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 
    !!----------------------------------------------------------------------
 
@@ -205,4 +216,6 @@
    Agrif_UseSpecialValue = .TRUE.
    CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn )
+   hbdy_w(:) = 0.e0 ; hbdy_e(:) = 0.e0 ; hbdy_n(:) = 0.e0 ; hbdy_s(:) = 0.e0
+   ssha(:,:) = 0.e0
 
    IF ( ln_dynspg_ts ) THEN
@@ -212,8 +225,8 @@
       CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b)
       CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b)
-      ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 ; hbdy_w(:) =0.e0
-      ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 ; hbdy_e(:) =0.e0 
-      ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 ; hbdy_n(:) =0.e0 
-      ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 ; hbdy_s(:) =0.e0
+      ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0
+      ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0
+      ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0
+      ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0
    ENDIF
 
@@ -234,5 +247,5 @@
          WRITE(cl_check2,*)  NINT(rdt)
          WRITE(cl_check3,*)  NINT(Agrif_Parent(rdt)/Agrif_Rhot())
-         CALL ctl_stop( 'incompatible time step between ocean grids',   &
+         CALL ctl_stop( 'Incompatible time step between ocean grids',   &
                &               'parent grid value : '//cl_check1    ,   & 
                &               'child  grid value : '//cl_check2    ,   & 
@@ -245,39 +258,36 @@
          WRITE(cl_check1,*)  (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1
          WRITE(cl_check2,*)   Agrif_Parent(nitend)   *Agrif_IRhot()
-         CALL ctl_warn( 'incompatible run length between grids'               ,   &
-               &              ' nit000 on fine grid will be change to : '//cl_check1,   &
-               &              ' nitend on fine grid will be change to : '//cl_check2    )
+         CALL ctl_warn( 'Incompatible run length between grids'                      ,   &
+               &               'nit000 on fine grid will be changed to : '//cl_check1,   &
+               &               'nitend on fine grid will be changed to : '//cl_check2    )
          nit000 = (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1
          nitend =  Agrif_Parent(nitend)   *Agrif_IRhot()
       ENDIF
 
-      ! Check coordinates
-     !SF  IF( ln_zps ) THEN
-     !SF     ! check parameters for partial steps 
-     !SF     IF( Agrif_Parent(e3zps_min) .NE. e3zps_min ) THEN
-     !SF        WRITE(*,*) 'incompatible e3zps_min between grids'
-     !SF        WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_min)
-     !SF        WRITE(*,*) 'child grid  :',e3zps_min
-     !SF        WRITE(*,*) 'those values should be identical'
-     !SF        STOP
-     !SF     ENDIF
-     !SF     IF( Agrif_Parent(e3zps_rat) /= e3zps_rat ) THEN
-     !SF        WRITE(*,*) 'incompatible e3zps_rat between grids'
-     !SF        WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_rat)
-     !SF        WRITE(*,*) 'child grid  :',e3zps_rat
-     !SF        WRITE(*,*) 'those values should be identical'                  
-     !SF        STOP
-     !SF     ENDIF
-     !SF  ENDIF
-
       ! Check free surface scheme
       IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.&
          & ( Agrif_Parent(ln_dynspg_exp).AND.ln_dynspg_ts ) ) THEN
-         WRITE(*,*) 'incompatible free surface scheme between grids'
-         WRITE(*,*) 'parent grid ln_dynspg_ts  :', Agrif_Parent(ln_dynspg_ts )
-         WRITE(*,*) 'parent grid ln_dynspg_exp :', Agrif_Parent(ln_dynspg_exp)
-         WRITE(*,*) 'child grid  ln_dynspg_ts  :', ln_dynspg_ts
-         WRITE(*,*) 'child grid  ln_dynspg_exp :', ln_dynspg_exp
-         WRITE(*,*) 'those logicals should be identical'                  
+         WRITE(cl_check1,*)  Agrif_Parent( ln_dynspg_ts )
+         WRITE(cl_check2,*)  ln_dynspg_ts
+         WRITE(cl_check3,*)  Agrif_Parent( ln_dynspg_exp )
+         WRITE(cl_check4,*)  ln_dynspg_exp
+         CALL ctl_stop( 'Incompatible free surface scheme between grids' ,  &
+               &               'parent grid ln_dynspg_ts  :'//cl_check1  ,  & 
+               &               'child  grid ln_dynspg_ts  :'//cl_check2  ,  &
+               &               'parent grid ln_dynspg_exp :'//cl_check3  ,  &
+               &               'child  grid ln_dynspg_exp :'//cl_check4  ,  &
+               &               'those logicals should be identical' )                 
+         STOP
+      ENDIF
+
+      ! Check if identical linear free surface option
+      IF ( ( Agrif_Parent(ln_linssh ).AND.(.NOT.ln_linssh )).OR.&
+         & ( (.NOT.Agrif_Parent(ln_linssh)).AND.ln_linssh ) ) THEN
+         WRITE(cl_check1,*)  Agrif_Parent(ln_linssh )
+         WRITE(cl_check2,*)  ln_linssh
+         CALL ctl_stop( 'Incompatible linearized fs option between grids',  &
+               &               'parent grid ln_linssh  :'//cl_check1     ,  &
+               &               'child  grid ln_linssh  :'//cl_check2     ,  &
+               &               'those logicals should be identical' )                  
          STOP
       ENDIF
@@ -306,27 +316,39 @@
    ENDIF
    ! 
-   ! Do update at initialisation because not done before writing restarts
-   ! This would indeed change boundary conditions values at initial time
-   ! hence produce restartability issues.
-   ! Note that update below is recursive (with lk_agrif_doupd=T):
-   ! 
-! JC: I am not sure if Agrif_MaxLevel() is the "relative"
-!     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
-      CALL Agrif_Update_tra()
-      CALL Agrif_Update_dyn()
-   ENDIF
-   !
+END SUBROUTINE Agrif_InitValues_cont
+
+RECURSIVE SUBROUTINE Agrif_Update_ini( )
+   !!----------------------------------------------------------------------
+   !!                 *** ROUTINE agrif_Update_ini ***
+   !!
+   !! ** Purpose :: Recursive update done at initialization
+   !!----------------------------------------------------------------------
+   USE dom_oce
+   USE agrif_opa_update
+#if defined key_top
+   USE agrif_top_update
+#endif
+   !
+   IMPLICIT NONE
+   !!----------------------------------------------------------------------
+   !
+   IF (Agrif_Root()) RETURN
+   !
+   IF (.NOT.ln_linssh) CALL Agrif_Update_vvl()
+   CALL Agrif_Update_tra()
+#if defined key_top
+   CALL Agrif_Update_Trc()
+#endif
+   CALL Agrif_Update_dyn()
 # if defined key_zdftke
-   CALL Agrif_Update_tke(0)
-# endif
-   !
-   Agrif_UseSpecialValueInUpdate = .FALSE.
-   nbcline = 0
-   lk_agrif_doupd = .FALSE.
-   !
-END SUBROUTINE Agrif_InitValues_cont
-
+! JC remove update because this precludes from perfect restartability
+!!   CALL Agrif_Update_tke()
+# endif
+
+   CALL Agrif_ChildGrid_To_ParentGrid()
+   CALL Agrif_Update_ini()
+   CALL Agrif_ParentGrid_To_ChildGrid()
+
+END SUBROUTINE agrif_update_ini
 
 SUBROUTINE agrif_declare_var
@@ -371,8 +393,8 @@
    CALL agrif_declare_variable((/2,2/),(/3,3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id)
 
-# if defined key_zdftke
-   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id)
-   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avt_id)
-   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avm_id)
+# if defined key_zdftke || defined key_zdfgls
+   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/), en_id)
+   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avt_id)
+   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avm_id)
 # endif
 
@@ -400,5 +422,5 @@
    CALL Agrif_Set_bcinterp(vmsk_id,interp=AGRIF_constant)
 
-# if defined key_zdftke
+# if defined key_zdftke || defined key_zdfgls
    CALL Agrif_Set_bcinterp(avm_id ,interp=AGRIF_linear)
 # endif
@@ -411,7 +433,4 @@
    CALL Agrif_Set_bc(vn_interp_id,(/0,1/))
 
-!   CALL Agrif_Set_bc(tsn_sponge_id,(/-3*Agrif_irhox(),0/))
-!   CALL Agrif_Set_bc(un_sponge_id,(/-2*Agrif_irhox()-1,0/))
-!   CALL Agrif_Set_bc(vn_sponge_id,(/-2*Agrif_irhox()-1,0/))
    CALL Agrif_Set_bc(tsn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/))
    CALL Agrif_Set_bc(un_sponge_id ,(/-nn_sponge_len*Agrif_irhox()-1,0/))
@@ -428,21 +447,37 @@
    CALL Agrif_Set_bc(vmsk_id,(/0,0/))
 
+# if defined key_zdftke || defined key_zdfgls
+   CALL Agrif_Set_bc(avm_id ,(/0,1/))
+# endif
+
+   ! 4. Update type
+   !--------------- 
+   CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average)
+
+# if defined UPD_HIGH
+   CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting)
+   CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting)
+   CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average)
+
+   CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting)
+   CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average)
+   CALL Agrif_Set_Updatetype(sshn_id, update = Agrif_Update_Full_Weighting)
+   CALL Agrif_Set_Updatetype(e3t_id, update = Agrif_Update_Full_Weighting)
+
 # if defined key_zdftke
-   CALL Agrif_Set_bc(avm_id ,(/0,1/))
-# endif
-
-   ! 5. Update type
-   !--------------- 
+   CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Full_Weighting)
+   CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Full_Weighting)
+   CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Full_Weighting)
+# endif
+
+#else
    CALL Agrif_Set_Updatetype(tsn_id, update = AGRIF_Update_Average)
-
-   CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average)
-
    CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average)
    CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy)
 
-   CALL Agrif_Set_Updatetype(sshn_id, update = AGRIF_Update_Average)
-
    CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average)
    CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy)
+   CALL Agrif_Set_Updatetype(sshn_id, update = AGRIF_Update_Average)
+   CALL Agrif_Set_Updatetype(e3t_id, update = AGRIF_Update_Average)
 
 # if defined key_zdftke
@@ -452,13 +487,5 @@
 # endif
 
-! High order updates
-!   CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting)
-!   CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting)
-!   CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average)
-!
-!   CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting)
-!   CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average)
-!   CALL Agrif_Set_Updatetype(sshn_id, update = Agrif_Update_Full_Weighting)
-
+#endif
    !
 END SUBROUTINE agrif_declare_var
@@ -733,21 +760,4 @@
       ENDIF
 
-      ! Check coordinates
-      IF( ln_zps ) THEN
-         ! check parameters for partial steps 
-         IF( Agrif_Parent(e3zps_min) .NE. e3zps_min ) THEN
-            WRITE(*,*) 'incompatible e3zps_min between grids'
-            WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_min)
-            WRITE(*,*) 'child grid  :',e3zps_min
-            WRITE(*,*) 'those values should be identical'
-            STOP
-         ENDIF
-         IF( Agrif_Parent(e3zps_rat) .NE. e3zps_rat ) THEN
-            WRITE(*,*) 'incompatible e3zps_rat between grids'
-            WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_rat)
-            WRITE(*,*) 'child grid  :',e3zps_rat
-            WRITE(*,*) 'those values should be identical'                  
-            STOP
-         ENDIF
       ENDIF
       ! Check passive tracer cell
@@ -756,9 +766,4 @@
       ENDIF
    ENDIF
-
-   CALL Agrif_Update_trc(0)
-   !
-   Agrif_UseSpecialValueInUpdate = .FALSE.
-   nbcline_trc = 0
    !
 END SUBROUTINE Agrif_InitValues_cont_top
@@ -792,14 +797,13 @@
    !-----------------------------
    CALL Agrif_Set_bc(trn_id,(/0,1/))
-!   CALL Agrif_Set_bc(trn_sponge_id,(/-3*Agrif_irhox(),0/))
    CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/))
 
-   ! 5. Update type
+   ! 4. Update type
    !--------------- 
+# if defined UPD_HIGH
+   CALL Agrif_Set_Updatetype(trn_id, update = Agrif_Update_Full_Weighting)
+#else
    CALL Agrif_Set_Updatetype(trn_id, update = AGRIF_Update_Average)
-
-!   Higher order update
-!   CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting)
-
+#endif
    !
 END SUBROUTINE agrif_declare_var_top
@@ -866,4 +870,6 @@
    IF (lk_mpp) iminspon = MIN(iminspon,FLOOR(REAL(jpi-2)/REAL(Agrif_irhox())), FLOOR(REAL(jpj-2)/REAL(Agrif_irhox())) )
    IF (nn_sponge_len > iminspon)  CALL ctl_stop('agrif sponge length is too large')
+   ! Check update frequency
+   IF (MOD((nitend-nit000+1), nbclineupdate).NE.0 ) CALL ctl_stop('number of time steps should be a multiple of nn_cln_update')
    !
    IF( agrif_oce_alloc()  > 0 )   CALL ctl_warn('agrif agrif_oce_alloc: allocation of arrays failed')
@@ -878,5 +884,5 @@
 SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob )
    !!----------------------------------------------------------------------
-   !!                     *** ROUTINE Agrif_detect ***
+   !!                     *** ROUTINE Agrif_InvLoc ***
    !!----------------------------------------------------------------------
    USE dom_oce
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 8741)
@@ -276,4 +276,12 @@
                ENDIF
             END DO
+#if defined key_agrif 
+	    IF( .NOT. AGRIF_Root() ) THEN 
+	       IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east 
+	       IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west 
+ 	       IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north 
+	       IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south 
+	    ENDIF 
+#endif 
          END DO
          !
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 8741)
@@ -1029,5 +1029,5 @@
       !
 #if defined key_agrif
-      IF(.NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' )
+      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
 #endif
       !
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 8741)
@@ -132,6 +132,6 @@
             ! so that asselin contribution is removed at the same time 
             DO jk = 1, jpkm1
-               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
-               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
+               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
             END DO  
          ENDIF
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 8741)
@@ -52,4 +52,5 @@
 #if defined key_agrif
    USE agrif_opa_interp ! agrif
+   USE agrif_oce
 #endif
 #if defined key_asminc   
@@ -76,5 +77,5 @@
 
    !! Time filtered arrays at baroclinic time step:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv     !: Advection vel. at "now" barocl. step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv     !: Advection fluxes at "now" barocl. step
 
    !! * Substitutions
@@ -127,5 +128,5 @@
       !!      -Update the filtered free surface at step "n+1"      : ssha
       !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
-      !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
+      !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
       !!      These are used to advect tracers and are compliant with discrete
       !!      continuity equation taken at the baroclinic time steps. This 
@@ -1019,15 +1020,12 @@
       !
       ! Set advection velocity correction:
-      zwx(:,:) = un_adv(:,:)
-      zwy(:,:) = vn_adv(:,:)
-      IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN     
-         un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:)
-         vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
-      ELSE
-         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)
-         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)
-      END IF
-
-      IF( ln_bt_fw ) THEN ! Save integrated transport for next computation
+      IF (ln_bt_fw) THEN
+         zwx(:,:) = un_adv(:,:)
+         zwy(:,:) = vn_adv(:,:)
+         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
+            un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) 
+            vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) )
+         END IF
+         ! Save integrated transport for next computation
          ub2_b(:,:) = zwx(:,:)
          vb2_b(:,:) = zwy(:,:)
@@ -1065,11 +1063,11 @@
       DO jk = 1, jpkm1
          ! Correct velocities:
-         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)
-         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
+         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)
+         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
          !
       END DO
       !
-      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
-      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic i-current
+      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu_n(:,:) )    ! barotropic i-current
+      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv_n(:,:) )    ! barotropic i-current
       !
 #if defined key_agrif
@@ -1297,5 +1295,5 @@
 #if defined key_agrif
       ! Restrict the use of Agrif to the forward case only
-      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
+!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
 #endif
       !
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 8741)
@@ -109,11 +109,10 @@
       ! 
       ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
-
+      !
+#if defined key_agrif
+      CALL agrif_ssh( kt )
+#endif
+      !
       IF ( .NOT.ln_dynspg_ts ) THEN
-         ! These lines are not necessary with time splitting since
-         ! boundary condition on sea level is set during ts loop
-# if defined key_agrif
-         CALL agrif_ssh( kt )
-# endif
          IF( ln_bdy ) THEN
             CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
@@ -214,4 +213,13 @@
          END DO
       ENDIF
+      !
+#if defined key_agrif 
+      IF( .NOT. AGRIF_Root() ) THEN 
+         IF ((nbondi ==  1).OR.(nbondi == 2)) wn(nlci-1 , :     ,:) = 0.e0      ! east 
+         IF ((nbondi == -1).OR.(nbondi == 2)) wn(2      , :     ,:) = 0.e0      ! west 
+         IF ((nbondj ==  1).OR.(nbondj == 2)) wn(:      ,nlcj-1 ,:) = 0.e0      ! north 
+         IF ((nbondj == -1).OR.(nbondj == 2)) wn(:      ,2      ,:) = 0.e0      ! south 
+      ENDIF 
+#endif 
       !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8741)
@@ -33,5 +33,8 @@
    USE iom            ! I/O manager library
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
+#if defined key_agrif
+   USE agrif_opa_interp ! Set bc on avm
+#endif 
 
    IMPLICIT NONE
@@ -204,10 +207,12 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
-                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &
-                  &                            / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
-                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
-                  &                            / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+               avmu(ji,jj,jk) = 0.5_wp * ( avm(ji,jj,jk) + avm(ji+1,jj,jk) )   &
+                  &                    * ( un(ji,jj,jk-1) - un(ji,jj,jk)   )   &
+                  &                    * ( ub(ji,jj,jk-1) - ub(ji,jj,jk)   )   &
+                  &                    / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
+               avmv(ji,jj,jk) = 0.5_wp * ( avm(ji,jj,jk) + avm(ji,jj+1,jk) )   & 
+                  &                    * ( vn(ji,jj,jk-1) - vn(ji,jj,jk)   )   &
+                  &                    * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )     &
+                  &                    / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
                eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
             END DO
@@ -800,4 +805,9 @@
       ! Lateral boundary conditions (sign unchanged)
       avt(:,:,1)  = 0._wp
+      !
+#if defined key_agrif 
+      CALL Agrif_avm 
+#endif
+      !
       CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
 
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 8741)
@@ -173,9 +173,4 @@
       !!----------------------------------------------------------------------
       !
-#if defined key_agrif 
-      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
-      IF( .NOT.Agrif_Root() )   CALL Agrif_Tke
-#endif
-      !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
          avt (:,:,:) = avt_k (:,:,:) 
@@ -196,5 +191,6 @@
 #if defined key_agrif
       ! Update child grid f => parent grid 
-      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
+!!! JC: suppress update since restartability is not possible in that case
+!!!      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( )      ! children only
 #endif      
      ! 
@@ -666,4 +662,9 @@
          END DO
       END DO
+      !
+#if defined key_agrif 
+      CALL Agrif_avm 
+#endif
+      !
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8741)
@@ -315,6 +315,10 @@
       IF( Agrif_NbStepint() == 0 ) THEN               ! AGRIF Update 
 !!jc in fact update is useless at last time step, but do it for global diagnostics
-                         CALL Agrif_Update_Tra()      ! Update active tracers
-                         CALL Agrif_Update_Dyn()      ! Update momentum
+         IF(.NOT.ln_linssh)  CALL Agrif_Update_vvl()      ! Update vertical scale factors  
+                             CALL Agrif_Update_Tra()      ! Update active tracers
+                             CALL Agrif_Update_Dyn()      ! Update momentum
+#if defined key_top
+                             CALL Agrif_Update_Trc()      ! Update passive tracers
+#endif
       ENDIF
 #endif
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8741)
@@ -113,4 +113,7 @@
    USE agrif_opa_sponge ! Momemtum and tracers sponges
    USE agrif_opa_update ! Update (2-way nesting)
+#if defined key_top
+   USE agrif_top_update ! passive tracers update (2-way nesting)
+#endif
 #endif
 #if defined key_top
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 8740)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 8741)
@@ -30,5 +30,4 @@
 #if defined key_agrif
    USE agrif_top_sponge ! tracers sponges
-   USE agrif_top_update ! tracers updates
 #endif
 
@@ -83,7 +82,4 @@
          IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kt )      ! internal damping trends on closed seas only
 
-#if defined key_agrif
-         IF( .NOT.Agrif_Root()) CALL Agrif_Update_Trc( kt ) ! Update tracer at AGRIF zoom boundaries : children only
-#endif
          !
       ELSE                                               ! 1D vertical configuration
