Index: /trunk/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 1501)
+++ /trunk/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 1502)
@@ -7,4 +7,5 @@
    !!             -   !  2007-07  (D. Storkey) Move Flather implementation to separate routine.
    !!            3.0  !  2008-04  (NEMO team)  add in the reference version
+   !!            3.2  !  2008-04  (R. Benshila) consider velocity instead of transport 
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -86,6 +87,6 @@
          END DO 
          !
-         CALL lbc_lnk( ua, 'U', 1. )   ! Boundary points should be updated
-         CALL lbc_lnk( va, 'V', 1. )   !
+         CALL lbc_lnk( ua, 'U', -1. )   ! Boundary points should be updated
+         CALL lbc_lnk( va, 'V', -1. )   !
          !
       ENDIF ! ln_bdy_dyn_frs
@@ -96,5 +97,5 @@
 #if defined key_dynspg_exp || defined key_dynspg_ts
 !! Option to use Flather with dynspg_flt not coded yet...
-   SUBROUTINE bdy_dyn_fla
+   SUBROUTINE bdy_dyn_fla( pssh )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn_fla  ***
@@ -116,7 +117,10 @@
       !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
       !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh
+
       INTEGER  ::   ib, igrd                         ! dummy loop indices
       INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
       REAL(wp) ::   zcorr                            ! Flather correction
+      REAL(wp) ::   zforc                            ! temporary scalar
       !!----------------------------------------------------------------------
 
@@ -127,46 +131,45 @@
       IF(ln_bdy_dyn_fla .OR. ln_bdy_tides) THEN ! If these are both false, then this routine does nothing. 
 
-      ! Fill temporary array with ssh data (here spgu):
-      igrd = 1
-      spgu(:,:) = 0.0
-      DO ib = 1, nblenrim(igrd)
-         ii = nbi(ib,igrd)
-         ij = nbj(ib,igrd)
-         IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib)
-         IF( ln_bdy_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(ib)
-      END DO
+         ! Fill temporary array with ssh data (here spgu):
+         igrd = 1
+         spgu(:,:) = 0.0
+         DO ib = 1, nblenrim(igrd)
+            ii = nbi(ib,igrd)
+            ij = nbj(ib,igrd)
+            IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib)
+            IF( ln_bdy_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(ib)
+         END DO
+         !
+         igrd = 2      ! Flather bc on u-velocity; 
+         !             ! remember that flagu=-1 if normal velocity direction is outward
+         !             ! I think we should rather use after ssh ?
+         DO ib = 1, nblenrim(igrd)
+            ii  = nbi(ib,igrd)
+            ij  = nbj(ib,igrd) 
+            iim1 = ii + MAX( 0, INT( flagu(ib) ) )   ! T pts i-indice inside the boundary
+            iip1 = ii - MIN( 0, INT( flagu(ib) ) )   ! T pts i-indice outside the boundary 
+            !
+            zcorr = - flagu(ib) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
+            zforc = ubtbdy(ib) + utide(ib)
+            ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
+         END DO
+         !
+         igrd = 3      ! Flather bc on v-velocity
+         !             ! remember that flagv=-1 if normal velocity direction is outward
+         DO ib = 1, nblenrim(igrd)
+            ii  = nbi(ib,igrd)
+            ij  = nbj(ib,igrd) 
+            ijm1 = ij + MAX( 0, INT( flagv(ib) ) )   ! T pts j-indice inside the boundary
+            ijp1 = ij - MIN( 0, INT( flagv(ib) ) )   ! T pts j-indice outside the boundary 
+            !
+            zcorr = - flagv(ib) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
+            zforc = vbtbdy(ib) + vtide(ib)
+            va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
+         END DO
+         CALL lbc_lnk( ua_e, 'U', 1. )   ! Boundary points should be updated
+         CALL lbc_lnk( va_e, 'V', 1. )   !
+         !
+      ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides
       !
-      igrd = 2      ! Flather bc on u-velocity; 
-      !             ! remember that flagu=-1 if normal velocity direction is outward
-      !             ! I think we should rather use after ssh ?
-      DO ib = 1, nblenrim(igrd)
-         ii  = nbi(ib,igrd)
-         ij  = nbj(ib,igrd) 
-         iim1 = ii + MAX( 0, INT( flagu(ib) ) )   ! T pts i-indice inside the boundary
-         iip1 = ii - MIN( 0, INT( flagu(ib) ) )   ! T pts i-indice outside the boundary 
-         !
-         zcorr = - flagu(ib) * SQRT( grav / (hu_e(ii, ij) + 1.e-20) ) * ( sshn_e(iim1, ij) - spgu(iip1,ij) )
-         ua_e(ii, ij) = ( ubtbdy(ib) + utide(ib) ) * hu_e(ii,ij)
-         ua_e(ii,ij) = ua_e(ii,ij) + zcorr * umask(ii,ij,1) * hu_e(ii,ij)
-      END DO
-      !
-      igrd = 3      ! Flather bc on v-velocity
-      !             ! remember that flagv=-1 if normal velocity direction is outward
-      DO ib = 1, nblenrim(igrd)
-         ii  = nbi(ib,igrd)
-         ij  = nbj(ib,igrd) 
-         ijm1 = ij + MAX( 0, INT( flagv(ib) ) )   ! T pts j-indice inside the boundary
-         ijp1 = ij - MIN( 0, INT( flagv(ib) ) )   ! T pts j-indice outside the boundary 
-         !
-         zcorr = - flagv(ib) * SQRT( grav / (hv_e(ii, ij) + 1.e-20) ) * ( sshn_e(ii, ijm1) - spgu(ii,ijp1) )
-         va_e(ii, ij) = ( vbtbdy(ib) + vtide(ib) ) * hv_e(ii,ij)
-         va_e(ii,ij) = va_e(ii,ij) + zcorr * vmask(ii,ij,1) * hv_e(ii,ij)
-      END DO
-      !
-      CALL lbc_lnk( ua_e, 'U', 1. ) ! Boundary points should be updated
-      CALL lbc_lnk( va_e, 'V', 1. ) !
-      !
-      ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides
-
    END SUBROUTINE bdy_dyn_fla
 #endif
Index: /trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1501)
+++ /trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 1502)
@@ -1,8 +1,7 @@
 MODULE dynnxt
-   !!======================================================================
+   !!=========================================================================
    !!                       ***  MODULE  dynnxt  ***
    !! Ocean dynamics: time stepping
-   !!======================================================================
-   !!======================================================================
+   !!=========================================================================
    !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
    !!                 !  1990-10  (C. Levy, G. Madec)
@@ -15,13 +14,15 @@
    !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
    !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines. 
-   !!            3.2  !  2009-04  (G. Madec, R.Benshila))  re-introduce the vvl option
-   !!----------------------------------------------------------------------
+   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
+   !!-------------------------------------------------------------------------
   
-   !!----------------------------------------------------------------------
-   !!   dyn_nxt      : update the horizontal velocity from the momentum trend
-   !!----------------------------------------------------------------------
+   !!-------------------------------------------------------------------------
+   !!   dyn_nxt      : obtain the next (after) horizontal velocity
+   !!-------------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O manager
+   USE dynspg_oce      ! type of surface pressure gradient
+   USE dynadv          ! dynamics: vector invariant versus flux form
+   USE domvvl          ! variable volume
    USE obc_oce         ! ocean open boundary conditions
    USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
@@ -31,10 +32,9 @@
    USE bdydta          ! unstructured open boundary conditions
    USE bdydyn          ! unstructured open boundary conditions
-   USE dynspg_oce      ! type of surface pressure gradient
+   USE agrif_opa_update
+   USE agrif_opa_interp
+   USE in_out_manager  ! I/O manager
    USE lbclnk          ! lateral boundary condition (or mpp link)
    USE prtctl          ! Print control
-   USE agrif_opa_update
-   USE agrif_opa_interp
-   USE domvvl          ! variable volume
 
    IMPLICIT NONE
@@ -57,22 +57,33 @@
       !!                  ***  ROUTINE dyn_nxt  ***
       !!                   
-      !! ** Purpose :   Compute the after horizontal velocity from the 
-      !!      momentum trend.
-      !!
-      !! ** Method  :   Apply lateral boundary conditions on the trends (ua,va) 
-      !!      through calls to routine lbc_lnk.
-      !!      After velocity is compute using a leap-frog scheme environment:
-      !!         (ua,va) = (ub,vb) + 2 rdt (ua,va)
-      !!      Note that if lk_dynspg_flt=T, the time stepping has already been
-      !!      performed in dynspg module
-      !!      Time filter applied on now horizontal velocity to avoid the
-      !!      divergence of two consecutive time-steps and swap of dynamics
-      !!      arrays to start the next time step:
-      !!         (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
-      !!         (un,vn) = (ua,va) 
-      !!
-      !! ** Action : - Update ub,vb arrays, the before horizontal velocity
-      !!             - Update un,vn arrays, the now horizontal velocity
-      !!
+      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary 
+      !!             condition on the after velocity, achieved the time stepping 
+      !!             by applying the Asselin filter on now fields and swapping 
+      !!             the fields.
+      !!
+      !! ** Method  : * After velocity is compute using a leap-frog scheme:
+      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
+      !!             Note that with flux form advection and variable volume layer
+      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
+      !!             velocity.
+      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
+      !!             the time stepping has already been done in dynspg module
+      !!
+      !!              * Apply lateral boundary conditions on after velocity 
+      !!             at the local domain boundaries through lbc_lnk call,
+      !!             at the radiative open boundaries (lk_obc=T),
+      !!             at the relaxed   open boundaries (lk_bdy=T), and
+      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
+      !!
+      !!              * Apply the time filter applied and swap of the dynamics
+      !!             arrays to start the next time step:
+      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
+      !!                (un,vn) = (ua,va).
+      !!             Note that with flux form advection and variable volume layer
+      !!             (lk_vvl=T), the time filter is applied on thickness weighted
+      !!             velocity.
+      !!
+      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
+      !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
@@ -85,5 +96,5 @@
       REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         - 
       REAL(wp) ::   zuf   , zvf              !    -         - 
-     !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
 
       IF( kt == nit000 ) THEN
@@ -93,49 +104,55 @@
       ENDIF
 
-      ! Local constant initialization
-      z2dt = 2. * rdt
+#if defined key_dynspg_flt
+      !
+      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
+      ! -------------
+
+      ! Update after velocity on domain lateral boundaries      (only local domain required)
+      ! --------------------------------------------------
+      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
+      CALL lbc_lnk( va, 'V', -1. ) 
+      !
+#else
+      ! Next velocity :   Leap-frog time stepping
+      ! -------------
+      z2dt = 2. * rdt                                 ! Euler or leap-frog time step 
       IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
-
-      !! Explicit physics with thickness weighted updates
-
-      ! Lateral boundary conditions on ( ua, va )
-      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. ) 
-
-      ! Next velocity
-      ! -------------
-#if defined key_dynspg_flt
-      ! Leap-frog time stepping already done in dynspg_flt.F routine
-#else
-      IF( lk_vvl ) THEN                                  ! Varying levels
+      !
+      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
          DO jk = 1, jpkm1
-            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)     &
-               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) )   &
-               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
-            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)     &
-               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) )   &
-               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
-         END DO
-      ELSE
-         DO jk = 1, jpkm1
-                  ! Leap-frog time stepping
             ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
             va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
          END DO
-      ENDIF
-
+      ELSE                                            ! applied on thickness weighted velocity
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
+               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
+               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
+            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
+               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
+               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
+         END DO
+      ENDIF
+
+
+      ! Update after velocity on domain lateral boundaries
+      ! --------------------------------------------------      
+      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
+      CALL lbc_lnk( va, 'V', -1. ) 
+      !
 # if defined key_obc
-      ! Update (ua,va) along open boundaries (only in the rigid-lid case)
+      !                                !* OBC open boundaries
       CALL obc_dyn( kt )
-
+      !
       IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
-         ! Flather boundary condition :
-         !        - Update sea surface height on each open boundary
-         !                 sshn (= after ssh) for explicit case
-         !                 sshn_b (= after ssha_b) for time-splitting case
-         !        - Correct the barotropic velocities
+         ! Flather boundary condition : - Update sea surface height on each open boundary
+         !                                       sshn   (= after ssh   ) for explicit case
+         !                                       sshn_b (= after ssha_b) for time-splitting case
+         !                              - Correct the barotropic velocities
          CALL obc_dyn_bt( kt )
          !
-         ! Boundary conditions on sshn ( after ssh)
-         CALL lbc_lnk( sshn, 'T', 1. )
+!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
+         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
          !
          IF( ln_vol_cst )   CALL obc_vol( kt )
@@ -143,38 +160,40 @@
          IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
       ENDIF
-
+      !
 # elif defined key_bdy 
-      ! Update (ua,va) along open boundaries (for exp or ts options).
-      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
+      !                                !* BDY open boundaries
+      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option
          !
          CALL bdy_dyn_frs( kt )
          !
-         IF( ln_bdy_fla ) THEN
+         IF( ln_bdy_dyn_fla ) THEN
             ua_e(:,:) = 0.e0
             va_e(:,:) = 0.e0
             ! Set these variables for use in bdy_dyn_fla
-            hu_e(:,:) = hu(:,:)
-            hv_e(:,:) = hv(:,:)
+            hur_e(:,:) = hur(:,:)
+            hvr_e(:,:) = hvr(:,:)
             DO jk = 1, jpkm1   !! Vertically integrated momentum trends
                ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
                va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
             END DO
+            ua_e(:,:) = ua_e(:,:) * hur(:,:)
+            va_e(:,:) = va_e(:,:) * hvr(:,:)
             DO jk = 1 , jpkm1
-               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)
-               va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)
+               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
+               va(:,:,jk) = va(:,:,jk) - va_e(:,:)
             END DO
             CALL bdy_dta_bt( kt+1, 0)
-            CALL bdy_dyn_fla
+            CALL bdy_dyn_fla( sshn_b )
+            DO jk = 1 , jpkm1
+               ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
+               va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
+            END DO
          ENDIF
          !
-         DO jk = 1 , jpkm1
-            ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)
-            va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)
-         END DO
       ENDIF
 # endif
-
+      !
 # if defined key_agrif
-      CALL Agrif_dyn( kt )
+      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
 # endif
 #endif
@@ -182,12 +201,26 @@
       ! Time filter and swap of dynamics arrays
       ! ------------------------------------------
-      IF( neuler == 0 .AND. kt == nit000 ) THEN
-         DO jk = 1, jpkm1                                 
-            un(:,:,jk) = ua(:,:,jk)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
+         DO jk = 1, jpkm1
+            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
             vn(:,:,jk) = va(:,:,jk)
          END DO
-      ELSE
-         IF( lk_vvl ) THEN                                ! Varying levels
-            DO jk = 1, jpkm1                              ! filter applied on thickness weighted velocities
+      ELSE                                             !* Leap-Frog : Asselin filter and swap
+         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity
+            DO jk = 1, jpkm1                              
+               DO jj = 1, jpj
+                  DO ji = 1, jpi    
+                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
+                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
+                     !
+                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
+                     vb(ji,jj,jk) = zvf
+                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+                     vn(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ELSE                                                ! applied on thickness weighted velocity
+            DO jk = 1, jpkm1
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -206,22 +239,12 @@
                      zve3b = vb(ji,jj,jk) * ze3v_b
                      !
-                     ub(ji,jj,jk) = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
-                        &         / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
-                     vb(ji,jj,jk) = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
-                        &         / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
-                     un(ji,jj,jk) = ua(ji,jj,jk) * umask(ji,jj,jk)
-                     vn(ji,jj,jk) = va(ji,jj,jk) * vmask(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ELSE                                             ! Fixed levels
-            DO jk = 1, jpkm1                              ! filter applied on velocities
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
-                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
-                     ub(ji,jj,jk) = zuf
+                     zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
+                        & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
+                     zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
+                        & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
+                     !
+                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
                      vb(ji,jj,jk) = zvf
-                     un(ji,jj,jk) = ua(ji,jj,jk)
+                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
                      vn(ji,jj,jk) = va(ji,jj,jk)
                   END DO
@@ -232,4 +255,5 @@
 
 #if defined key_agrif
+      ! Update velocity at AGRIF zoom boundaries
       IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
 #endif      
@@ -240,4 +264,4 @@
    END SUBROUTINE dyn_nxt
 
-   !!======================================================================
+   !!=========================================================================
 END MODULE dynnxt
Index: /trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 1501)
+++ /trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 1502)
@@ -40,12 +40,10 @@
 
 #if   defined key_dynspg_ts   ||  defined key_vvl   ||  defined key_esopa
-   !! Time splitting variables
-   !! ------------------------
-      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop
-         un_e  , vn_e,                        & ! vertically integrated horizontal velocities (now)
-         ua_e  , va_e                           ! vertically integrated horizontal velocities (after)
-      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop
-         sshn_e, ssha_e,                      & ! sea surface heigth (now, after)
-         hu_e  , hv_e                           ! depth arrays for the barotropic solution
+  !! Time splitting variables   ! variables of the explicit barotropic loop
+  !! ------------------------
+     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ua_e  , va_e             ! barotropic velocities (after)
+     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshn_e, ssha_e, sshn_b   ! sea surface heigth (now, after, average)
+     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hu_e  , hv_e             ! depth arrays for the barotropic solution
+     REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hur_e , hvr_e            ! inverse of depth arrays 
 #endif
 
Index: /trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1501)
+++ /trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1502)
@@ -1,14 +1,14 @@
 MODULE dynspg_ts
    !!======================================================================
-   !! History :   1.0  !  04-12  (L. Bessieres, G. Madec)  Original code
-   !!              -   !  05-11  (V. Garnier, G. Madec)  optimization
-   !!              -   !  06-08  (S. Masson)  distributed restart using iom
-   !!             2.0  !  07-07  (D. Storkey) calls to BDY routines
-   !!              -   !  08-01  (R. Benshila)  change averaging method
+   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
+   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
+   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
+   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
+   !!              -   ! 2008-01  (R. Benshila)  change averaging method
+   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
   !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
    !!----------------------------------------------------------------------
    !!   'key_dynspg_ts'         free surface cst volume with time splitting
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
@@ -45,6 +45,9 @@
    PUBLIC ts_rst      ! routine called by istate.F90
 
+
    REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter
    REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme)
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   un_b, vn_b   ! averaged velocity
 
    !! * Substitutions
@@ -66,5 +69,4 @@
       !!      gradient in case of free surface formulation with time-splitting.
       !!      Add it to the general trend of momentum equation.
-      !!      Compute the free surface.
       !!
       !! ** Method  :   Free surface formulation with time-splitting
@@ -75,11 +77,11 @@
       !!          the barotropic integration.
       !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and 
-      !!          barotropic transports (ua_e and va_e) through barotropic 
+      !!          barotropic velocity (ua_e and va_e) through barotropic 
       !!          momentum and continuity integration. Barotropic former 
       !!          variables are time averaging over the full barotropic cycle
-      !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b 
+      !!          (= 2 * baroclinic time step) and saved in zuX_b 
       !!          and zvX_b (X specifying after, now or before).
       !!      -3- The new general trend becomes :
-      !!          ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H
+      !!          ua = ua - sum_k(ua)/H + ( ua_e - sum_k(ub) )
       !!
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
@@ -87,29 +89,23 @@
       !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
       !!---------------------------------------------------------------------
-      INTEGER, INTENT( in )  ::   kt           ! ocean time-step index
-      !!
-      INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices
-      INTEGER  ::  icycle                      ! temporary scalar
-      REAL(wp) ::                           &
-         zraur, zcoef, z2dt_e, z2dt_b, zfac25,   &  ! temporary scalars
-         zfact1, zspgu, zcubt, zx1, zy1,    &  !     "        "
-         zfact2, zspgv, zcvbt, zx2, zy2        !     "        "
-      REAL(wp), DIMENSION(jpi,jpj) ::       &
-         zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays
-         zua, zva, zub, zvb,                &  !     "        "
-         zua_e, zva_e, zssha_e,             &  !     "        "
-         zub_e, zvb_e, zsshb_e,             &  !     "        "
-         zu_sum, zv_sum, zssh_sum
-      !! Variable volume
-      REAL(wp), DIMENSION(jpi,jpj)     ::   &
-         zspgu_1, zspgv_1, zsshun_e, zsshvn_e                     ! 2D workspace
+      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
+      INTEGER  ::   icycle                      ! temporary scalar
+      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b,   &  ! temporary scalars
+         z1_8, zspgu, zcubt, zx1, zy1,    &  !     "        "
+         z1_4, zspgv, zcvbt, zx2, zy2        !     "        "
+      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e
+      !! 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcu, zwx, zua, zun, zub   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcv, zwy, zva, zvn, zvb   !  -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zun_e, zub_e, zu_sum      ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zvn_e, zvb_e, zv_sum      !  -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssh_sum                  !  -      -
       !!----------------------------------------------------------------------
 
-      ! Arrays initialization
-      ! ---------------------
-      zhdiv(:,:) = 0.e0
-
-
-      IF( kt == nit000 ) THEN
+      IF( kt == nit000 ) THEN             !* initialisation
          !
          IF(lwp) WRITE(numout,*)
@@ -117,13 +113,14 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
          IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
-
+         !
          CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields:
-         !                               ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b
-
-         zssha_e(:,:) = sshn(:,:)
-         zua_e  (:,:) = un_e(:,:)
-         zva_e  (:,:) = vn_e(:,:)
-         hu_e   (:,:) = hu  (:,:)
-         hv_e   (:,:) = hv  (:,:)
+         !                               ! un_b, vn_b  
+         !
+         ua_e  (:,:) = un_b (:,:)
+         va_e  (:,:) = vn_b (:,:)
+         hu_e  (:,:) = hu   (:,:)
+         hv_e  (:,:) = hv   (:,:)
+         hur_e (:,:) = hur  (:,:)
+         hvr_e (:,:) = hvr  (:,:)
          IF( ln_dynvor_een ) THEN
             ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
@@ -140,87 +137,57 @@
       ENDIF
 
-      ! Substract the surface pressure gradient from the momentum
-      ! ---------------------------------------------------------
-      IF( lk_vvl ) THEN    ! Variable volume
-
-         ! 0. Local initialization
-         ! -----------------------
-         zspgu_1(:,:) = 0.e0
-         zspgv_1(:,:) = 0.e0
-
-         ! 1. Surface pressure gradient (now)
-         ! ----------------------------
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  ) &
-                  &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e1u(ji,jj)
-               zspgv_1(ji,jj) = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1) &
-                  &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e2v(ji,jj)
-            END DO
-         END DO
-
-         ! 2. Substract the surface pressure trend from the general trend
-         ! ------------------------------------------------------
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj)
-                  va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj)
-               END DO
-            END DO
-         END DO
-
-      ENDIF
-
-      ! Local constant initialization
-      ! --------------------------------
+      !                                   !* Local constant initialization
       z2dt_b = 2.0 * rdt                                    ! baroclinic time step
-      IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt
-      zfact1 = 0.5 * 0.25                                   ! coefficient for vorticity estimates
-      zfact2 = 0.5 * 0.5
+      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
+      z1_4 = 0.5 * 0.5
       zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water
+      !
+      zhdiv(:,:) = 0.e0                                     ! barotropic divergence
 
       ! -----------------------------------------------------------------------------
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
       ! -----------------------------------------------------------------------------
-
-      ! Vertically integrated quantities
-      ! --------------------------------
-      zua(:,:) = 0.e0
-      zva(:,:) = 0.e0
-      zub(:,:) = 0.e0
-      zvb(:,:) = 0.e0
-      zwx(:,:) = 0.e0
-      zwy(:,:) = 0.e0
-
-      ! vertical sum
-      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-         DO jk = 1, jpkm1
+      !      
+      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
+      !                                   ! --------------------------
+      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   zub(:,:) = 0.e0
+      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   zvb(:,:) = 0.e0
+      !
+      DO jk = 1, jpkm1
+#if defined key_vectopt_loop
+         DO jj = 1, 1         !Vector opt. => forced unrolling
             DO ji = 1, jpij
-               !                                                           ! Vertically integrated momentum trends
-               zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk)
-               zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk)
-               !                                                           ! Vertically integrated transports (before)
-               zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
-               zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
-               !                                                           ! Planetary vorticity transport fluxes (now)
-               zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk)
-               zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk)
-            END DO
-         END DO
-      ELSE                             ! No  vector opt.
-         DO jk = 1, jpkm1
-            !                                                           ! Vertically integrated momentum trends
-            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
-            !                                                           ! Vertically integrated transports (before)
-            zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
-            zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
-            !                                                           ! Planetary vorticity (now)
-            zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
-            zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         END DO
-      ENDIF
-
+#else 
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+#endif
+               !                                                                              ! now trend
+               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
+               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
+               !                                                                              ! now velocity 
+               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
+               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
+               !                                                                              ! before velocity
+               zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) 
+               zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
+      DO jk = 1, jpkm1                    ! --------------------------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
+               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
+            END DO
+         END DO
+      END DO
+
+      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
+      !                                   ! ---------------------------====
+      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport 
+      zwy(:,:) = zvn(:,:) * e1v(:,:)
+      !
       IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
          DO jj = 2, jpjm1
@@ -231,6 +198,6 @@
                zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
                ! energy conserving formulation for planetary vorticity term
-               zcu(ji,jj) = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
-               zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
+               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
+               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
             END DO
          END DO
@@ -239,8 +206,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
-                              + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-               zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
-                              + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
+               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
                zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
                zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
@@ -249,13 +214,10 @@
          !
       ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
-         zfac25 = 0.25
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcu(ji,jj) = + zfac25 / e1u(ji,jj)   &
-                  &       * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                  &          + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zcv(ji,jj) = - zfac25 / e2v(ji,jj)   &
-                  &       * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                  &          + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
+               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
+                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
+               zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
+                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
             END DO
          END DO
@@ -263,19 +225,18 @@
       ENDIF
 
-
-      ! Remove barotropic trend from general momentum trend
-      ! ---------------------------------------------------
-      DO jk = 1 , jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
-               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
-            END DO
-         END DO
-      END DO
-
-      ! Remove coriolis term from barotropic trend
-      ! ------------------------------------------
-      DO jj = 2, jpjm1
+      !                                   !* Right-Hand-Side of the barotropic momentum equation
+      !                                   ! ----------------------------------------------------
+      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
+         DO jj = 2, jpjm1 
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    &
+                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
+               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    &
+                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
+            END DO
+         END DO
+      ENDIF
+
+      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
          DO ji = fs_2, fs_jpim1
             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
@@ -284,43 +245,46 @@
       END DO
 
+      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
+      !                                   ! -------------------------- 
+      zua(:,:) = zua(:,:) * hur(:,:)
+      zva(:,:) = zva(:,:) * hvr(:,:)
+      !
+      IF( lk_vvl ) THEN
+         zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
+         zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
+      ELSE
+         zub(:,:) = zub(:,:) * hur(:,:)
+         zvb(:,:) = zvb(:,:) * hvr(:,:)
+      ENDIF
+
       ! -----------------------------------------------------------------------
       !  Phase 2 : Integration of the barotropic equations with time splitting
       ! -----------------------------------------------------------------------
-
-      ! Initialisations
-      !----------------
-      ! Number of iteration of the barotropic loop
-      icycle = 2  * nn_baro + 1
-
-      ! variables for the barotropic equations
-      zu_sum  (:,:) = 0.e0
-      zv_sum  (:,:) = 0.e0
-      zssh_sum(:,:) = 0.e0
-      hu_e    (:,:) = hu    (:,:)      ! (barotropic) ocean depth at u-point
-      hv_e    (:,:) = hv    (:,:)      ! (barotropic) ocean depth at v-point
-      zsshb_e (:,:) = sshn_e(:,:)      ! (barotropic) sea surface height (before and now)
-      ! vertical sum
-      un_e  (:,:) = 0.e0
-      vn_e  (:,:) = 0.e0
-      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-         DO jk = 1, jpkm1
-            DO ji = 1, jpij
-               un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
-               vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
-            END DO
-         END DO
-      ELSE                             ! No  vector opt.
-         DO jk = 1, jpkm1
-            un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)
-            vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
-         END DO
-      ENDIF
-      zub_e (:,:) = un_e(:,:)
-      zvb_e (:,:) = vn_e(:,:)
-
-
+      !
+      !                                             ! ==================== !
+      !                                             !    Initialisations   !
+      !                                             ! ==================== !
+      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
+      
+      !                                ! Start from NOW field
+      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
+      hv_e   (:,:) = hv   (:,:) 
+      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
+      hvr_e  (:,:) = hvr  (:,:)
+      zsshb_e(:,:) = sshn (:,:)            ! sea surface height (before and now)
+      sshn_e (:,:) = sshn (:,:)
+      
+      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
+      zvn_e  (:,:) = vn_b (:,:)
+      zub_e  (:,:) = un_b (:,:)
+      zvb_e  (:,:) = vn_b (:,:)
+
+      zu_sum  (:,:) = un_b (:,:)           ! summation
+      zv_sum  (:,:) = vn_b (:,:)
+      zssh_sum(:,:) = sshn (:,:)
+
+#if defined key_obc
       ! set ssh corrections to 0
       ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
-#if defined key_obc
       IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
       IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
@@ -329,73 +293,61 @@
 #endif
 
-      ! Barotropic integration over 2 baroclinic time steps
-      ! ---------------------------------------------------
-
-      !                                                    ! ==================== !
-      DO jit = 1, icycle                                   !  sub-time-step loop  !
-         !                                                 ! ==================== !
+      !                                             ! ==================== !
+      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
+         !                                          ! ==================== !
          z2dt_e = 2. * ( rdt / nn_baro )
-         IF ( jit == 1 )   z2dt_e = rdt / nn_baro
-
-         ! Time interpolation of open boundary condition data
-         IF( lk_obc )   CALL obc_dta_bt( kt, jit )
-         IF( lk_bdy .OR. ln_bdy_tides)   CALL bdy_dta_bt( kt, jit+1 )
-
-         ! Horizontal divergence of barotropic transports
-         !--------------------------------------------------
-         zhdiv(:,:) = 0.e0
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * un_e(ji  ,jj)              &
-                  &            -e2u(ji-1,jj  ) * un_e(ji-1,jj)              &
-                  &            +e1v(ji  ,jj  ) * vn_e(ji  ,jj)              &
-                  &            -e1v(ji  ,jj-1) * vn_e(ji  ,jj-1) )          &
-                  &           / (e1t(ji,jj)*e2t(ji,jj))
-            END DO
-         END DO
-
+         IF( jn == 1 )   z2dt_e = rdt / nn_baro
+
+         !                                                !* Update the forcing (OBC, BDY and tides)
+         !                                                !  ------------------
+         IF( lk_obc                     )   CALL obc_dta_bt( kt, jn   )
+         IF( lk_bdy  .OR.  ln_bdy_tides )   CALL bdy_dta_bt( kt, jn+1 )
+
+         !                                                !* after ssh_e
+         !                                                !  -----------
+         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
+                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
+                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
+                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
+            END DO
+         END DO
+         !
 #if defined key_obc
-         ! open boundaries (div must be zero behind the open boundary)
-         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
-         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0      ! east
-         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0      ! west
+         !                                                     ! OBC : zhdiv must be zero behind the open boundary
+!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
+         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
+         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
-         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0      ! south
+         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
 #endif
-
 #if defined key_bdy
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
-            END DO
-         END DO
+         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
 #endif
-
-         ! Sea surface height from the barotropic system
-         !----------------------------------------------
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
-                  &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
-            END DO
-         END DO
-
-         ! evolution of the barotropic transport ( following the vorticity scheme used)
-         ! ----------------------------------------------------------------------------
-         zwx(:,:) = e2u(:,:) * un_e(:,:)
-         zwy(:,:) = e1v(:,:) * vn_e(:,:)
-
-         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
+         !
+         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
+            END DO
+         END DO
+
+         !                                                !* after barotropic velocities (vorticity scheme dependent)
+         !                                                !  ---------------------------  
+         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport
+         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
+         !
+         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
                   IF( lk_vvl) THEN
-                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) &
-                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj)
-                     zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) &
-                        &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj)
+                     zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
+                        &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
+                     zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
+                        &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
                   ELSE
-                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
+                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
                   ENDIF
                   ! energy conserving formulation for planetary vorticity term
@@ -404,128 +356,115 @@
                   zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
                   zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-                  zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
-                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
-                  ! after transports
-                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  zcubt = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
+                  zcvbt =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
+                  ! after barotropic velocity
+                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
             !
-         ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
+         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                   ! surface pressure gradient
+                  IF( lk_vvl) THEN
+                     zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
+                        &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
+                     zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
+                        &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
+                  ELSE
+                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
+                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
+                  ENDIF
+                  ! enstrophy conserving formulation for planetary vorticity term
+                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
+                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
+                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
+                  ! after barotropic velocity
+                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+               END DO
+            END DO
+            !
+         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
                   IF( lk_vvl) THEN
-                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) &
-                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj)
-                     zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) &
-                        &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj)
+                     zspgu = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
+                        &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
+                     zspgv = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
+                        &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
                   ELSE
-                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
-                  ENDIF
-                  ! enstrophy conserving formulation for planetary vorticity term
-                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
-                     + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
-                     + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
-                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
-                  ! after transports
-                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
-               END DO
-            END DO
-            !
-         ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
-            zfac25 = 0.25
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     zspgu = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  ) &
-                        &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj)
-                     zspgv = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) &
-                        &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj)
-                  ELSE
-                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
+                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
                   ENDIF
                   ! energy/enstrophy conserving formulation for planetary vorticity term
-                  zcubt = + zfac25 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                     &                             + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-                  zcvbt = - zfac25 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                     &                             + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
-                  ! after transports
-                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  zcubt = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
+                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
+                  zcvbt = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
+                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
+                  ! after barotropic velocity
+                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
             ! 
          ENDIF
-
-         ! Flather's boundary condition for the barotropic loop :
-         !         - Update sea surface height on each open boundary
-         !         - Correct the barotropic transports
-         IF( lk_obc )   CALL obc_fla_ts
-         IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla
-
-         ! ... Boundary conditions on ua_e, va_e, ssha_e
-         CALL lbc_lnk( zua_e  , 'U', -1. )
-         CALL lbc_lnk( zva_e  , 'V', -1. )
-         CALL lbc_lnk( zssha_e, 'T',  1. )
-
-         ! temporal sum
-         !-------------
-         zu_sum  (:,:) = zu_sum  (:,:) + zua_e  (:,:)
-         zv_sum  (:,:) = zv_sum  (:,:) + zva_e  (:,:) 
-         zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 
-
-         ! Time filter and swap of dynamics arrays
-         ! ---------------------------------------
-         IF( jit == 1 ) THEN   ! Euler (forward) time stepping
-            zsshb_e(:,:) = sshn_e (:,:)
-            zub_e  (:,:) = un_e   (:,:)
-            zvb_e  (:,:) = vn_e   (:,:)
-            sshn_e (:,:) = zssha_e(:,:)
-            un_e   (:,:) = zua_e  (:,:)
-            vn_e   (:,:) = zva_e  (:,:)
-         ELSE                                        ! Asselin filtering
-            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:)
-            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * un_e  (:,:)
-            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * vn_e  (:,:)
-            sshn_e (:,:) = zssha_e(:,:)
-            un_e   (:,:) = zua_e  (:,:)
-            vn_e   (:,:) = zva_e  (:,:)
+         !                                                !* domain lateral boundary
+         !                                                !  -----------------------
+         !                                                      ! Flather's boundary condition for the barotropic loop :
+         !                                                      !         - Update sea surface height on each open boundary
+         !                                                      !         - Correct the velocity
+
+         IF( lk_obc                   )   CALL obc_fla_ts
+         IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla( sshn_e ) 
+         !
+         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries 
+         CALL lbc_lnk( va_e  , 'V', -1. )
+         CALL lbc_lnk( ssha_e, 'T',  1. )
+
+         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
+         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
+         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
+
+         !                                                !* Time filter and swap
+         !                                                !  --------------------
+         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
+            zsshb_e(:,:) = sshn_e(:,:)
+            zub_e  (:,:) = zun_e (:,:)
+            zvb_e  (:,:) = zvn_e (:,:)
+            sshn_e (:,:) = ssha_e(:,:)
+            zun_e  (:,:) = ua_e  (:,:)
+            zvn_e  (:,:) = va_e  (:,:)
+         ELSE                                                   ! Swap + Filter
+            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
+            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
+            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
+            sshn_e (:,:) = ssha_e(:,:)
+            zun_e  (:,:) = ua_e  (:,:)
+            zvn_e  (:,:) = va_e  (:,:)
          ENDIF
 
-         IF( lk_vvl ) THEN ! Variable volume   ! needed for BDY ???
-
-            ! Sea surface elevation
-            ! ---------------------
-            DO jj = 1, jpjm1
-               DO ji = 1,jpim1
-
-                  ! Sea Surface Height at u-point before
-                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )   &
-                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
-                     &                + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
-
-                  ! Sea Surface Height at v-point before
-                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )   &
-                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
-                     &                + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
-
-               END DO
-            END DO
-
-            ! Boundaries conditions
-            CALL lbc_lnk( zsshun_e, 'U', 1. )    ;  CALL lbc_lnk( zsshvn_e, 'V', 1. )
-
-            ! Ocean depth at U- and V-points
-            ! -------------------------------------------
-            hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:)
-            hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:)
-
+         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
+            !                                             !  ------------------
+            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
+               DO ji = 1, fs_jpim1   ! Vector opt.
+                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
+                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
+                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
+                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
+                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
+                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
+            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
+            !
+            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
+            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
+            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
+            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
             !
          ENDIF
@@ -534,12 +473,6 @@
       !                                                    ! ==================== !
 
-      ! Time average of after barotropic variables
-      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
-      un_e  (:,:) = zcoef *  zu_sum  (:,:) 
-      vn_e  (:,:) = zcoef *  zv_sum  (:,:) 
-      sshn_e(:,:) = zcoef *  zssh_sum(:,:) 
-
 #if defined key_obc
-      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
+      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
       IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
       IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
@@ -548,21 +481,21 @@
 
       ! -----------------------------------------------------------------------------
-      ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step)
+      ! Phase 3. update the general trend with the barotropic trend
       ! -----------------------------------------------------------------------------
-
-
-
-      ! add time averaged barotropic coriolis and surface pressure gradient
-      ! terms to the general momentum trend
-      ! --------------------------------------------------------------------
+      !
+      !                                   !* Time average ==> after barotropic u, v, ssh
+      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
+      un_b  (:,:) = zcoef * zu_sum  (:,:) 
+      vn_b  (:,:) = zcoef * zv_sum  (:,:) 
+      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
+      ! 
+      !                                   !* update the general momentum trend
       DO jk=1,jpkm1
-         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b
-         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b
+         ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b
+         va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b
       END DO
-
-      ! write filtered free surface arrays in restart file
-      ! --------------------------------------------------
+      !
+      !                                   !* write time-spliting arrays in the restart
       IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
-
       !
    END SUBROUTINE dyn_spg_ts
@@ -582,31 +515,30 @@
       !
       IF( TRIM(cdrw) == 'READ' ) THEN
-         IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) )   ! free surface and
-            CALL iom_get( numror, jpdom_autoglo, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued
-            CALL iom_get( numror, jpdom_autoglo, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop
+         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
+            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
          ELSE
-            sshn_e(:,:) = sshn(:,:)
-            un_e  (:,:) = 0.e0
-            vn_e  (:,:) = 0.e0
+            un_b (:,:) = 0.e0
+            vn_b (:,:) = 0.e0
             ! vertical sum
             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
                DO jk = 1, jpkm1
                   DO ji = 1, jpij
-                     un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
-                     vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
+                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
+                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
                   END DO
                END DO
             ELSE                             ! No  vector opt.
                DO jk = 1, jpkm1
-                  un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)
-                  vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
+                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
+                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
                END DO
             ENDIF
+            un_b (:,:) = un_b(:,:) * hur(:,:)
+            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
          ENDIF
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) )   ! free surface and
-         CALL iom_rstput( kt, nitrst, numrow, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued
-         CALL iom_rstput( kt, nitrst, numrow, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop
+         CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! external velocity issued
+         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
       ENDIF
       !
