Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4292)
@@ -767,6 +767,6 @@
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
-          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad)
-          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad)
+          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation
+          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation
         END DO
       END DO
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4292)
@@ -17,4 +17,5 @@
    !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
    !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!-------------------------------------------------------------------------
   
@@ -42,4 +43,6 @@
    USE wrk_nemo        ! Memory Allocation
    USE prtctl          ! Print control
+   USE dynspg_ts       ! Barotropic velocities
+
 #if defined key_agrif
    USE agrif_opa_interp
@@ -103,4 +106,5 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
       !!----------------------------------------------------------------------
@@ -109,4 +113,5 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
+      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
       !
       IF( kt == nit000 ) THEN
@@ -127,4 +132,6 @@
       !
 #else
+
+# if defined key_dynspg_exp
       ! Next velocity :   Leap-frog time stepping
       ! -------------
@@ -147,5 +154,41 @@
          END DO
       ENDIF
-
+# endif
+
+# if defined key_dynspg_ts
+      ! Ensure below that barotropic velocities match time splitting estimate
+      ! Compute actual transport and replace it with ts estimate at "after" time step
+      zua(:,:) = 0._wp
+      zva(:,:) = 0._wp
+      IF (lk_vvl) THEN
+         DO jk = 1, jpkm1
+            zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)    
+         END DO
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_e(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_e(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+         END DO
+      ELSE
+         DO jk = 1, jpkm1
+            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)    
+         END DO
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk)
+         END DO
+      ENDIF
+
+      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
+         ! Remove advective velocity from "now velocities" 
+         ! prior to asselin filtering     
+         ! In the forward case, this is done below after asselin filtering    
+         DO jk = 1, jpkm1
+            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
+            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+         END DO  
+      ENDIF
+# endif
 
       ! Update after velocity on domain lateral boundaries
@@ -194,4 +237,11 @@
             vn(:,:,jk) = va(:,:,jk)
          END DO
+         IF (lk_vvl) THEN
+            DO jk = 1, jpkm1
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
+               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
+            ENDDO
+         ENDIF
       ELSE                                             !* Leap-Frog : Asselin filter and swap
          !                                ! =============!
@@ -201,6 +251,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi    
-                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
-                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
+                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
+                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -214,23 +264,30 @@
          ELSE                             ! Variable volume !
             !                             ! ================!
+            ! Before scale factor at t-points
+            ! (used as a now filtered scale factor until the swap)
+            ! ----------------------------------------------------
+            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
+               ! Remove asselin filtering on thicknesses if forward time splitting
+                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
+            ELSE
+               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
+               ! Add volume filter correction: compatibility with tracer advection scheme
+               ! => time filter + conservation correction (only at the first level)
+               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
             !
-            DO jk = 1, jpkm1                 ! Before scale factor at t-points
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
-                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
-                  &                         - 2._wp * fse3t_n(:,:,jk)            )
-            END DO
-            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
-            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            ENDIF
             !
-            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
-               !
-               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b)
-               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
-               !
-               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity
-                  DO jj = 1, jpj                      !                                                 --------
+            IF( ln_dynadv_vec ) THEN
+               ! Before scale factor at (u/v)-points
+               ! -----------------------------------
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
+               ! Leap-Frog - Asselin filter and swap: applied on velocity
+               ! -----------------------------------
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
-                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
+                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
+                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                         !
                         ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -242,10 +299,13 @@
                END DO
                !
-            ELSE                             ! flux form (thickness weighted calulation)
-               !
-               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b)
-               !
-               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: 
-                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity
+            ELSE
+               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
+               !------------------------------------------------
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
+               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
+               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
+               ! -----------------------------------             ===========================
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
                      DO ji = 1, jpi                   !                              ---------------------------
                         zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
@@ -272,5 +332,27 @@
          ENDIF
          !
-      ENDIF
+         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
+         ! Remove asselin filtering of barotropic velocities if forward time splitting
+         ! note that we replace barotropic velocities by advective velocities       
+            zua(:,:) = 0._wp
+            zva(:,:) = 0._wp
+            IF (lk_vvl) THEN
+               DO jk = 1, jpkm1
+                  zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
+                  zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+               END DO
+            ELSE
+               DO jk = 1, jpkm1
+                  zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
+                  zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+               END DO
+            ENDIF
+            DO jk = 1, jpkm1
+               ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
+               vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
+            END DO
+         ENDIF
+         !
+      ENDIF ! neuler =/0
 
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
@@ -278,4 +360,5 @@
       ! 
       CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
+      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 4292)
@@ -23,4 +23,7 @@
    USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine)
    USE dynadv         ! dynamics: vector invariant versus flux form
+   USE dynhpg, ONLY: ln_dynhpg_imp
+   USE sbctide
+   USE updtide
    USE trdmod         ! ocean dynamics trends
    USE trdmod_oce     ! ocean variables trends
@@ -101,15 +104,57 @@
       ENDIF
 
-      IF( ln_apr_dyn ) THEN                   !==  Atmospheric pressure gradient  ==!
-         zg_2 = grav * 0.5
-         DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
+      IF(      ln_apr_dyn                                                &   ! atmos. pressure
+         .OR.  ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting)
+         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice
+         !
+         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               spgu(ji,jj) =  zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
-                  &                   + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
-               spgv(ji,jj) =  zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
-                  &                   + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
-            END DO
-         END DO
-         DO jk = 1, jpkm1                          ! Add the apg to the general trend
+               spgu(ji,jj) = 0._wp
+               spgv(ji,jj) = 0._wp
+            END DO
+         END DO         
+         !
+         IF( ln_apr_dyn .AND. (.NOT. lk_dynspg_ts) ) THEN                    !==  Atmospheric pressure gradient (added later in time-split case) ==!
+            zg_2 = grav * 0.5
+            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
+                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
+                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
+                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         !
+         !                                    !==  tide potential forcing term  ==!
+         IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
+            !
+            CALL upd_tide( kt )                      ! update tide potential
+            !
+            DO jj = 2, jpjm1                         ! add tide potential forcing
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  spgv(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
+                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
+               END DO 
+            END DO
+         ENDIF
+         !
+         IF( nn_ice_embd == 2 ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
+            CALL wrk_alloc( jpi, jpj, zpice )
+            !                                            
+            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
+            zgrau0r     = - grav * r1_rau0
+            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
+                  spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
+               END DO
+            END DO
+            !
+            CALL wrk_dealloc( jpi, jpj, zpice )         
+         ENDIF
+         !
+         DO jk = 1, jpkm1                     !== Add all terms to the general trend
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -118,31 +163,6 @@
                END DO
             END DO
-         END DO
-      ENDIF
-
-      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
-         CALL wrk_alloc( jpi, jpj, zpice )
-         !                                            
-         zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
-         zgrau0r     = - grav * r1_rau0
-         zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               spgu(ji,jj) = ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
-               spgv(ji,jj) = ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
-            END DO
-         END DO
-         DO jk = 1, jpkm1                             ! Add the surface pressure trend to the general trend
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
-                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
-               END DO
-            END DO
-         END DO
-         !
-         CALL wrk_dealloc( jpi, jpj, zpice )
-      ENDIF
-
+         END DO         
+      ENDIF
 
       SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend
@@ -209,4 +229,6 @@
       ENDIF
 
+      IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 )
+      ! (do it now, to set nn_baro, used to allocate some arrays later on)
       !                        ! allocate dyn_spg arrays
       IF( lk_dynspg_ts ) THEN
@@ -248,7 +270,7 @@
       ENDIF
 
-      !                        ! Control of momentum formulation
-      IF( lk_dynspg_ts .AND. lk_vvl ) THEN
-         IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' )
+      !			       ! Control of hydrostatic pressure choice
+      IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN
+         CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )
       ENDIF
       !
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 4292)
@@ -91,5 +91,6 @@
                spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
             END DO 
-         END DO 
+         END DO
+         !
          DO jk = 1, jpkm1                    ! Add it to the general trend
             DO jj = 2, jpjm1
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 4292)
@@ -39,4 +39,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur_e , hvr_e    ! inverse of hu_e and hv_e
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_b           ! before field without time-filter
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b     ! after  averaged velocities
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b     ! now    averaged velocities
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b     ! before averaged velocities
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv ! Advection vel. at "now" barocl. step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b,  vb2_b  ! Advection vel. at "now-0.5" barocl. step
 
    !!----------------------------------------------------------------------
@@ -53,4 +58,8 @@
       ALLOCATE( sshn_e(jpi,jpj) , ua_e(jpi,jpj) , hu_e(jpi,jpj) , hur_e(jpi,jpj) ,      &
          &      ssha_e(jpi,jpj) , va_e(jpi,jpj) , hv_e(jpi,jpj) , hvr_e(jpi,jpj) ,      &
+         &      ub_b(jpi,jpj)   , vb_b(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj)  ,      &
+         &      ua_b(jpi,jpj)   , va_b(jpi,jpj)                                  ,      & 
+         &      ub2_b(jpi,jpj)  , vb2_b(jpi,jpj)                                 ,      &
+         &      un_adv(jpi,jpj) , vn_adv(jpi,jpj)                                ,      &
          &      sshn_b(jpi,jpj)                                                  , STAT = dynspg_oce_alloc )
          !
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4292)
@@ -9,4 +9,6 @@
    !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
    !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
+   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
+   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
    !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
@@ -16,5 +18,4 @@
    !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
    !!                 splitting scheme and add to the general trend 
-   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -24,16 +25,11 @@
    USE phycst          ! physical constants
    USE domvvl          ! variable volume
-   USE zdfbfr          ! bottom friction
    USE dynvor          ! vorticity term
-   USE obc_oce         ! Lateral open boundary condition
-   USE obc_par         ! open boundary condition parameters
-   USE obcdta          ! open boundary condition data     
-   USE obcfla          ! Flather open boundary condition  
    USE bdy_par         ! for lk_bdy
    USE bdy_oce         ! Lateral open boundary condition
-   USE bdydta          ! open boundary condition data     
+   USE bdytides        ! open boundary condition data     
    USE bdydyn2d        ! open boundary conditions on barotropic variables
-   USE sbctide
-   USE updtide
+   USE sbctide         ! tides
+   USE updtide         ! tide potential
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -41,7 +37,13 @@
    USE in_out_manager  ! I/O manager
    USE iom             ! IOM library
+   USE restart         ! only for lrst_oce
    USE zdf_oce         ! Vertical diffusion
    USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE timing          ! Timing    
+   USE sbcapr          ! surface boundary condition: atmospheric pressure
+   USE dynadv, ONLY: ln_dynadv_vec
+#if defined key_agrif
+   USE agrif_opa_interp ! agrif
+#endif
 
 
@@ -49,14 +51,32 @@
    PRIVATE
 
-   PUBLIC dyn_spg_ts        ! routine called by step.F90
-   PUBLIC ts_rst            ! routine called by istate.F90
-   PUBLIC dyn_spg_ts_alloc  ! routine called by dynspg.F90
-
-
+   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90 
+   PUBLIC dyn_spg_ts_alloc  !    "      "     "    "
+   PUBLIC dyn_spg_ts_init   !    "      "     "    "
+
+   ! Potential namelist parameters below to be read in dyn_spg_ts_init
+   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.TRUE.        !: Forward integration of barotropic sub-stepping
+   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_av=.TRUE.        !: Time averaging of barotropic variables
+   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE.  !: Set number of iterations automatically
+   INTEGER,  PRIVATE, PARAMETER :: nn_bt_flt=1            !: Filter choice
+   REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp      !: Max. courant number (used if ln_bt_nn_auto=T)
+   ! End namelist parameters
+
+   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
+   REAL(wp),SAVE :: rdtbt   ! Barotropic time step
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: &
+                    wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables
+                    wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b   ! now    averaged velocity
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b   ! before averaged velocity
+   ! Arrays below are saved to allow testing of the "no time averaging" option
+   ! If this option is not retained, these could be replaced by temporary arrays
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sshbb_e, sshb_e, & ! Instantaneous barotropic arrays
+                                                   ubb_e, ub_e,     &
+                                                   vbb_e, vb_e
 
    !! * Substitutions
@@ -64,6 +84,6 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
-   !! $Id$
+   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
+   !! $Id: dynspg_ts.F90
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -74,7 +94,19 @@
       !!                  ***  routine dyn_spg_ts_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( ftnw  (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) ,     &
-         &      ftsw  (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc )
-         !
+      INTEGER :: ierr(3)
+      !!----------------------------------------------------------------------
+      ierr(:) = 0
+
+      ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
+         &      ub_e(jpi,jpj)  , vb_e(jpi,jpj)   , &
+         &      ubb_e(jpi,jpj) , vbb_e(jpi,jpj)  , STAT= ierr(1) )
+
+      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
+
+      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
+                             &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
+
+      dyn_spg_ts_alloc = MAXVAL(ierr(:))
+
       IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc )
       IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
@@ -82,57 +114,87 @@
    END FUNCTION dyn_spg_ts_alloc
 
-
    SUBROUTINE dyn_spg_ts( kt )
       !!----------------------------------------------------------------------
-      !!                  ***  routine dyn_spg_ts  ***
       !!
-      !! ** Purpose :   Compute the now trend due to the surface pressure
-      !!      gradient in case of free surface formulation with time-splitting.
-      !!      Add it to the general trend of momentum equation.
+      !! ** Purpose :   
+      !!      -Compute the now trend due to the explicit time stepping
+      !!      of the quasi-linear barotropic system. Barotropic variables are
+      !!      advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F)
+      !!      or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a
+      !!      generalized forward-backward (see ref. below) time stepping.
+      !!      -Update the free surface at step "n+1" (ssha, zsshu_a, zsshv_a).
+      !!      -Compute barotropic advective velocities at step "n" to be used 
+      !!      to advect tracers latter on. These are compliant with discrete
+      !!      continuity equation taken at the baroclinic time steps, thus 
+      !!      ensuring tracers conservation.
       !!
-      !! ** Method  :   Free surface formulation with time-splitting
-      !!      -1- Save the vertically integrated trend. This general trend is
-      !!          held constant over the barotropic integration.
-      !!          The Coriolis force is removed from the general trend as the
-      !!          surface gradient and the Coriolis force are updated within
-      !!          the barotropic integration.
-      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and 
-      !!          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 uX_b 
-      !!          and vX_b (X specifying after, now or before).
-      !!      -3- The new general trend becomes :
-      !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b )
+      !! ** Method  :  
       !!
-      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
+      !! ** Action : - Update barotropic velocities: ua_b, va_b
+      !!             - Update trend (ua,va) with barotropic component
+      !!             - Update ssha, zsshu_a, zsshv_a
+      !!             - Update barotropic advective velocity at kt=now
       !!
-      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
+      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 
+      !!              The regional oceanic modeling system (ROMS): 
+      !!              a split-explicit, free-surface,
+      !!              topography-following-coordinate oceanic model. 
+      !!              Ocean Modelling, 9, 347-404. 
       !!---------------------------------------------------------------------
       !
       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::   icycle           ! local scalar
-      INTEGER  ::   ikbu, ikbv       ! local scalar
-      REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars
-      REAL(wp) ::   z1_8, zx1, zy1                            !   -      -
-      REAL(wp) ::   z1_4, zx2, zy2                            !   -      -
-      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      -
-      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      -
-      REAL(wp) ::   ua_btm, va_btm                            !   -      -
-      !
-      REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv 
-      REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 
-      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum
+      LOGICAL  ::   ll_fw_start			! if true, forward integration 
+      LOGICAL  ::   ll_init			! if true, special startup of 2d equations
+      INTEGER  ::   ji, jj, jk, jn   		! dummy loop indices
+      INTEGER  ::   ikbu, ikbv, noffset    	! local integers
+      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf  	! local scalars
+      REAL(wp) ::   zx1, zy1, zx2, zy2         	!   -      -
+      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2 	!   -      -
+      REAL(wp) ::   zu_spg, zv_spg		!   -      -
+      REAL(wp) ::   zhura, zhvra         	!   -      -
+      REAL(wp) ::   za0, za1, za2, za3		!   -      -
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
+      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
+      REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv
+      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
+      REAL(wp), POINTER, DIMENSION(:,:) :: zhur_b, zhvr_b
+      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
+      REAL(wp), POINTER, DIMENSION(:,:) :: zht, zhf
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
       !
-      CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     )
-      CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
-      CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
-      !
-      IF( kt == nit000 ) THEN             !* initialisation
+      !                                         !* Allocate temporay arrays
+      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv )
+      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  )
+      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
+      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
+      CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b                                     )
+      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   )
+      CALL wrk_alloc( jpi, jpj, zht, zhf )
+      !
+      !                                         !* Local constant initialization
+      z1_12 = 1._wp / 12._wp 
+      z1_8  = 0.125_wp                                   
+      z1_4  = 0.25_wp
+      z1_2  = 0.5_wp     
+      zraur = 1._wp / rau0
+      !
+      IF( kt == nit000 .AND. neuler == 0 ) THEN  	! reciprocal of baroclinic time step 
+        z2dt_bf = rdt
+      ELSE
+        z2dt_bf = 2.0_wp * rdt
+      ENDIF
+      z1_2dt_b = 1.0_wp / z2dt_bf 
+      !
+      ll_init = ln_bt_av                       		! if no time averaging, then no specific restart 
+      ll_fw_start = .FALSE.
+      !
+                                     	                ! time offset in steps for bdy data update
+      IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF
+      !
+      IF( kt == nit000 ) THEN             	!* initialisation
          !
          IF(lwp) WRITE(numout,*)
@@ -141,52 +203,125 @@
          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: 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._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
+         IF (neuler==0) ll_init=.TRUE.
+         !
+         IF (ln_bt_fw.OR.(neuler==0)) THEN
+           ll_fw_start=.TRUE.
+           noffset = 0
+         ELSE
+           ll_fw_start=.FALSE.
+         ENDIF
+         !
+         ! Set averaging weights and cycle length:
+         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
+         !
+         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) 
+         !
+      ENDIF
+      !
+      ! Set arrays to remove/compute coriolis trend.
+      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
+      ! Note that these arrays are also used during barotropic loop. These are however frozen
+      ! although they should be updated in variable volume case. Not a big approximation.
+      ! To remove this approximation, copy lines below inside barotropic loop
+      ! and update depths at T-F points (ht and hf resp.) at each barotropic time step
+      !
+      IF ( kt == nit000 .OR. lk_vvl ) THEN
+         IF ( ln_dynvor_een ) THEN
+            ! JC: Simplification needed below: define ht_0 even when volume is fixed
+            IF (lk_vvl) THEN
+               zht(:,:) = (ht_0(:,:) + sshn(:,:)) * tmask(:,:,1) 
+            ELSE
+               zht(:,:) = 0.
+               DO jk = 1, jpkm1
+                  zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+               END DO
+            ENDIF
+
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zwz(ji,jj) =   ( zht(ji  ,jj+1) + zht(ji+1,jj+1) +                     &
+                        &          zht(ji  ,jj  ) + zht(ji+1,jj  )   )                   &
+                        &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
+                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
+                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( zwz, 'F', 1._wp )
+            zwz(:,:) = ff(:,:) * zwz(:,:)
+
+            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
             DO jj = 2, jpj
                DO ji = fs_2, jpi   ! vector opt.
-                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp
-                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp
-                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp
-                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp
-               END DO
-            END DO
-         ENDIF
-         !
-      ENDIF
-
-      !                                                     !* Local constant initialization
-      z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                   ! reciprocal of baroclinic time step
-      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! reciprocal of baroclinic 
-                                                                        ! time step (euler timestep)
-      z1_8     = 0.125_wp                                   ! coefficient for vorticity estimates
-      z1_4     = 0.25_wp        
-      zraur    = 1._wp / rau0                               ! 1 / volumic mass
-      !
-      zhdiv(:,:) = 0._wp                                    ! barotropic divergence
-      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F)
-      zv_sld = 0._wp   ;   zv_asp = 0._wp
-
-      IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction
-        z2dt_bf = rdt
-      ELSE
-        z2dt_bf = 2.0_wp * rdt
-      ENDIF
-
+                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
+                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
+                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
+                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
+               END DO
+            END DO
+         ELSE
+            zwz(:,:) = 0._wp
+            zht(:,:) = 0.
+            IF ( .not. ln_sco ) THEN
+!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
+!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
+!              ENDIF
+!              zht(:,:) = gdepw_0(:,:,jk+1)
+            ELSE
+               zht(:,:) = hbatf(:,:)
+            END IF
+
+            DO jj = 1, jpjm1
+               zht(:,jj) = zht(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+            END DO
+
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1
+                  zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
+               END DO
+            END DO
+            CALL lbc_lnk( zht, 'F', 1._wp )
+            ! JC: TBC. hf should be greater than 0 
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( zht(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array
+               END DO
+            END DO
+            zwz(:,:) = ff(:,:) * zwz(:,:)
+         ENDIF
+      ENDIF
+      !
+      ! If forward start at previous time step, and centered integration, 
+      ! then update averaging weights:
+      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
+         ll_fw_start=.FALSE.
+         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
+      ENDIF
+
+      ! before inverse water column height at u- and v- points
+      IF( lk_vvl ) THEN
+         zhur_b(:,:) = 0.
+         zhvr_b(:,:) = 0.
+         DO jk = 1, jpk
+            zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+            zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+         END DO
+         zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) )
+         zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) )
+      ELSE
+         zhur_b(:,:) = hur(:,:)
+         zhvr_b(:,:) = hvr(:,:)
+      ENDIF
+                          
       ! -----------------------------------------------------------------------------
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
       ! -----------------------------------------------------------------------------
       !      
+      ! Some vertical sums (at now and before time steps) below could be suppressed 
+      ! if one swap barotropic arrays somewhere
+      !
       !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
-      !                                   ! --------------------------
-      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp
-      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
+      !                                   ! --------------------------------------------------
+      zu_frc(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp  ;  un_b(:,:) = 0._wp
+      zv_frc(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp  ;  vn_b(:,:) = 0._wp
       !
       DO jk = 1, jpkm1
@@ -198,36 +333,41 @@
             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)               
-               !
-#if defined key_vvl
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk) 
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk)
-#else
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
-#endif
+               !        ! now trend:                                                                   
+               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
+               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)         
+               !        ! now bt transp:                   
+               un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)        
+               vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               !	! before bt transp:
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-
+      !
+      zu_frc(:,:) = zu_frc(:,:) * hur(:,:)
+      zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
+      !
+      IF( lk_vvl ) THEN
+          ub_b(:,:) = ub_b(:,:) * zhur_b(:,:)
+          vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:)
+      ELSE
+          ub_b(:,:) = ub_b(:,:) * hur(:,:)
+          vb_b(:,:) = vb_b(:,:) * hvr(:,:)
+      ENDIF
+      !
       !                                   !* baroclinic momentum trend (remove the vertical mean trend)
-      DO jk = 1, jpkm1                    ! --------------------------
+      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)
+               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-
-      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
-      !                                   ! ---------------------------====
-      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport 
-      zwy(:,:) = zvn(:,:) * e1v(:,:)
+      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
+      !                                   ! --------------------------------------------------------
+      zwx(:,:) = un_b(:,:) * e2u(:,:)           ! now transport 
+      zwy(:,:) = vn_b(:,:) * e1v(:,:)
       !
       IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
@@ -239,250 +379,319 @@
                zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
                ! energy conserving formulation for planetary vorticity term
-               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
-         !
-      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
+               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
+            END DO
+         END DO
+         !
+      ELSEIF ( ln_dynvor_ens ) THEN             ! enstrophy conserving scheme
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               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) )
-            END DO
-         END DO
-         !
-      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
+               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)
+               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
+               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
+            END DO
+         END DO
+         !
+      ELSEIF ( ln_dynvor_een ) THEN             ! enstrophy and energy conserving scheme
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               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
-         !
-      ENDIF
-
+               zu_trd(ji,jj) = + z1_12 / 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) )
+               zv_trd(ji,jj) = - z1_12 / 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
+         !
+      ENDIF 
+      !
+      un_b (:,:) = un_b(:,:) * hur(:,:)         ! Revert now transport to barotropic velocities
+      vn_b (:,:) = vn_b(:,:) * hvr(:,:)  
       !                                   !* Right-Hand-Side of the barotropic momentum equation
       !                                   ! ----------------------------------------------------
-      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
+      IF( lk_vvl ) THEN                         ! Variable volume : remove 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
+               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
+               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(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)
-             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
+             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
+             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
           END DO
-      END DO
-
-                    
-      !                                             ! Remove barotropic contribution of bottom friction 
-      !                                             ! from the barotropic transport trend
-      zcoef = -1._wp * z1_2dt_b
-
-      IF(ln_bfrimp) THEN
-      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient
-      !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction 
-        DO jj = 2, jpjm1         
-           DO ji = fs_2, fs_jpim1
-              ikbu = mbku(ji,jj)
-              ikbv = mbkv(ji,jj)
-              ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu)
-              va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv)
-
-              zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm
-              zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm
-           END DO
-        END DO
-
-      ELSE
-
-# if defined key_vectopt_loop
-        DO jj = 1, 1
-           DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-# else
-        DO jj = 2, jpjm1
-           DO ji = 2, jpim1
-# endif
-            ! Apply stability criteria for bottom friction
-            !RBbug for vvl and external mode we may need to use varying fse3
-            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
-              zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
-              zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
-           END DO
-        END DO
-
-        IF( lk_vvl ) THEN
-           DO jj = 2, jpjm1
-              DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
-                    &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) )
-                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
-                    &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) )
-              END DO
-           END DO
-        ELSE
-           DO jj = 2, jpjm1
-              DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
-                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
-              END DO
-           END DO
-        ENDIF
-      END IF    ! end (ln_bfrimp)
-
-                    
-      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
-      !                                   ! -------------------------- 
-      zua(:,:) = zua(:,:) * hur(:,:)
-      zva(:,:) = zva(:,:) * hvr(:,:)
-      !
-      IF( lk_vvl ) THEN
-         ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
-         vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
-      ELSE
-         ub_b(:,:) = ub_b(:,:) * hur(:,:)
-         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-      ENDIF
-
+      END DO 
+      !
+      !						! Add bottom stress contribution from baroclinic velocities:      
+      IF (ln_bt_fw) THEN
+         DO jj = 2, jpjm1                          
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ikbu = mbku(ji,jj)       
+               ikbv = mbkv(ji,jj)    
+               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
+               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
+            END DO
+         END DO
+      ELSE
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ikbu = mbku(ji,jj)       
+               ikbv = mbkv(ji,jj)    
+               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
+               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
+      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
+      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
+      !       
+      IF (ln_bt_fw) THEN                        ! Add wind forcing
+         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
+         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
+      ELSE
+         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
+         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
+      ENDIF  
+      !
+      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
+         IF (ln_bt_fw) THEN
+            DO jj = 2, jpjm1              
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
+                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
+                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
+                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
+               END DO
+            END DO
+         ELSE
+            DO jj = 2, jpjm1              
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
+                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
+                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
+                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
+                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
+                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
+               END DO
+            END DO
+         ENDIF 
+      ENDIF
+      !                                   !* Right-Hand-Side of the barotropic ssh equation
+      !                                   ! -----------------------------------------------
+      !                                         ! Surface net water flux and rivers
+      IF (ln_bt_fw) THEN
+         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) )
+      ELSE
+         zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:))
+      ENDIF
+#if defined key_asminc
+      !                                         ! Include the IAU weighted SSH increment
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
+         zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:)
+      ENDIF
+#endif
+      !
       ! -----------------------------------------------------------------------
-      !  Phase 2 : Integration of the barotropic equations with time splitting
+      !  Phase 2 : Integration of the barotropic equations 
       ! -----------------------------------------------------------------------
       !
       !                                             ! ==================== !
       !                                             !    Initialisations   !
+      !                                             ! ==================== !  
+      ! Initialize barotropic variables:    
+      IF (ln_bt_fw) THEN                  ! FORWARD integration:  start from NOW fields                             
+         sshn_e (:,:) = sshn (:,:)            
+         zun_e  (:,:) = un_b (:,:)            
+         zvn_e  (:,:) = vn_b (:,:)
+      ELSE                                ! CENTERED integration: start from BEFORE fields
+         sshn_e (:,:) = sshb (:,:)
+         zun_e  (:,:) = ub_b (:,:)         
+         zvn_e  (:,:) = vb_b (:,:)
+      ENDIF
+      !
+      ! Initialize depths:
+      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN
+         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
+         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
+         hur_e (:,:) = zhur_b(:,:)
+         hvr_e (:,:) = zhvr_b(:,:)
+      ELSE
+         hu_e  (:,:) = hu   (:,:)       
+         hv_e  (:,:) = hv   (:,:) 
+         hur_e (:,:) = hur  (:,:)    
+         hvr_e (:,:) = hvr  (:,:)
+      ENDIF
+      !
+      IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5:
+         zhup2_e (:,:) = hu(:,:)
+         zhvp2_e (:,:) = hv(:,:)
+      ENDIF
+      !
+      ! Initialize sums:
+      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)          
+      va_b  (:,:) = 0._wp
+      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
+      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop
+      zv_sum(:,:) = 0._wp
       !                                             ! ==================== !
-      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  (:,:)
-!RBbug     zsshb_e(:,:) = sshn (:,:)  
-      zsshb_e(:,:) = sshn_b(:,:)           ! 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( lp_obc_east  )   sshfoe_b(:,:) = 0._wp
-      IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp
-      IF( lp_obc_south )   sshfos_b(:,:) = 0._wp
-      IF( lp_obc_north )   sshfon_b(:,:) = 0._wp
-#endif
-
-      !                                             ! ==================== !
-      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
+      DO jn = 1, icycle                             !  sub-time-step loop  !
          !                                          ! ==================== !
-         z2dt_e = 2. * ( rdt / nn_baro )
-         IF( jn == 1 )   z2dt_e = rdt / nn_baro
-
          !                                                !* Update the forcing (BDY and tides)
          !                                                !  ------------------
-         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
-         IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 )
-         IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn )
-
-         !                                                !* after ssh_e
+         ! Update only tidal forcing at open boundaries
+#if defined key_tide
+         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
+         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )
+#endif
+         !
+         ! Set extrapolation coefficients for predictor step:
+         IF ((jn<3).AND.ll_init) THEN      ! Forward           
+           za1 = 1._wp                                          
+           za2 = 0._wp                        
+           za3 = 0._wp                        
+         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105 
+           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
+           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
+           za3 =  0.281105_wp              ! za3 = bet
+         ENDIF
+
+         ! Extrapolate barotropic velocities at step jit+0.5:
+         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
+         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
+
+         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
+            !                                             !  ------------------
+            ! Extrapolate Sea Level at step jit+0.5:
+            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
+            !
+            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
+               DO ji = 2, fs_jpim1   ! Vector opt.
+                  zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
+                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
+                     &              +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
+                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
+                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
+                     &              +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( zwx, 'U', 1._wp )    ;   CALL lbc_lnk( zwy, 'V', 1._wp )
+            !
+            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)               ! Ocean depth at U- and V-points
+            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
+         ENDIF
+         !                                                !* after ssh
          !                                                !  -----------
-         DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports
+         ! One should enforce volume conservation at open boundaries here
+         ! considering fluxes below:
+         !
+         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
+         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
+         DO jj = 2, jpjm1                                 
             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
-         !                                                     ! 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._wp      ! east
-         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west
-         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north
-         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south
+               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
+                  &             + zwy(ji,jj) - zwy(ji,jj-1)   &
+                  &           ) / ( e1t(ji,jj) * e2t(ji,jj) )
+            END DO
+         END DO
+         !
+         ! Sum over sub-time-steps to compute advective velocities
+         za2 = wgtbtp2(jn)
+         zu_sum  (:,:) = zu_sum  (:,:) + za2 * ua_e  (:,:) * zhup2_e  (:,:)
+         zv_sum  (:,:) = zv_sum  (:,:) + za2 * va_e  (:,:) * zhvp2_e  (:,:)
+         !
+         ! Set next sea level:
+         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
+         CALL lbc_lnk( ssha_e, 'T',  1._wp )
+
+#if defined key_bdy
+         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
+         IF (lk_bdy) CALL bdy_ssh( ssha_e )
 #endif
-#if defined key_bdy
-         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
+#if defined key_agrif
+         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
 #endif
-         !
-         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)-rnf(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(:,:)
+         !  
+         ! Sea Surface Height at u-,v-points (vvl case only)
+         IF ( lk_vvl ) THEN                                
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1      ! NO Vector Opt.
+                  zsshu_a(ji,jj) = z1_2  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                 &
+                     &                                   * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha_e(ji  ,jj) &
+                     &                                     + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) )
+                  zsshv_a(ji,jj) = z1_2  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                 &
+                     &                                   * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha_e(ji,jj  ) &
+                     &                                     + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp )
+         ENDIF   
+         !                                 
+         ! Half-step back interpolation of SSH for surface pressure computation:
+         !----------------------------------------------------------------------
+         IF ((jn==1).AND.ll_init) THEN
+           za0=1._wp                        ! Forward-backward
+           za1=0._wp                           
+           za2=0._wp
+           za3=0._wp
+         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
+           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
+           za1=-0.1666666666666_wp          ! za1 = gam
+           za2= 0.0833333333333_wp          ! za2 = eps
+           za3= 0._wp              
+         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
+           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps    
+           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps 
+           za2=0.088_wp                     ! za2 = gam
+           za3=0.013_wp                     ! za3 = eps
+         ENDIF
+
+         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
+          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
+
+         !
+         ! Compute associated depths at U and V points:
+         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN      
+            !                                        
+            DO jj = 2, jpjm1                            
+               DO ji = 2, jpim1
+                  zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
+                    &        * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
+                    &        +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )                 
+                  zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
+                     &       * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
+                     &       +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
+                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
+               END DO
+            END DO
+         ENDIF
+         !
+         ! Add Coriolis trend:
+         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated
+         ! at each time step. We however keep them constant here for optimization.
+         ! Recall that zwx and zwy arrays hold fluxes at this stage:
+         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
+         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_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
-                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zv_spg = -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
-                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! add tidal astronomical forcing
-                  IF ( ln_tide_pot .AND. lk_tide ) THEN 
-                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! energy conserving formulation for planetary vorticity term
                   zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
                   zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
                   zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
                   zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
-                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
-                  ! after velocities with implicit bottom friction
-
-                  IF( ln_bfrimp ) THEN      ! implicit bottom friction
-                     !   A new method to implement the implicit bottom friction. 
-                     !   H. Liu
-                     !   Sept 2011
-                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
-                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                     !
-                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
-                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
-                     !
-                  ELSE
-                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
-                      &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
-                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
-                      &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
-                  ENDIF
+                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
                END DO
             END DO
@@ -491,45 +700,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                   ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zv_spg = -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
-                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! add tidal astronomical forcing
-                  IF ( ln_tide_pot .AND. lk_tide ) THEN
-                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(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)
-                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
-                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
-                  ! after velocities with implicit bottom friction
-                  IF( ln_bfrimp ) THEN
-                     !   A new method to implement the implicit bottom friction. 
-                     !   H. Liu
-                     !   Sept 2011
-                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
-                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                     !
-                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
-                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
-                     !
-                  ELSE
-                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
-                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
-                  ENDIF
+                  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)
+                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
+                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
                END DO
             END DO
@@ -538,149 +712,211 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! surface pressure gradient
-                  IF( lk_vvl) THEN
-                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
-                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
-                     zv_spg = -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
-                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
-                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! add tidal astronomical forcing
-                  IF ( ln_tide_pot .AND. lk_tide ) THEN
-                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
-                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
-                  ENDIF
-                  ! energy/enstrophy conserving formulation for planetary vorticity term
-                  zu_cor = + 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)
-                  zv_cor = - 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 velocities with implicit bottom friction
-                  IF( ln_bfrimp ) THEN
-                     !   A new method to implement the implicit bottom friction. 
-                     !   H. Liu
-                     !   Sept 2011
-                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
-                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
-                     !
-                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
-                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
-                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
-                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
-                     !
-                  ELSE
-                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
-                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
-                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
-                  ENDIF
+                  zu_trd(ji,jj) = + z1_12 / 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) )
+                  zv_trd(ji,jj) = - z1_12 / 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
             ! 
          ENDIF
-         !                                                     !* domain lateral boundary
-         !                                                     !  -----------------------
-
-                                                               ! OBC open boundaries
-         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
-
-                                                               ! BDY open boundaries
-#if defined key_bdy
-         pssh => sshn_e
-         phur => hur_e
-         phvr => hvr_e
-         pu2d => ua_e
-         pv2d => va_e
-
-         IF( lk_bdy )   CALL bdy_dyn2d( kt ) 
-#endif
-
-         !
-         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                                !* 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_wp * 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_wp * 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(:,:)
+         !
+         ! Add tidal astronomical forcing if defined
+         IF ( lk_tide.AND.ln_tide_pot ) THEN
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
+                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
+                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
+                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
+               END DO
+            END DO
+         ENDIF
+         !
+         ! Add bottom stresses:
+         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
+         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
+         !
+         ! Surface pressure trend:
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Add surface pressure gradient
+               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
+               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
+               zwx(ji,jj) = zu_spg
+               zwy(ji,jj) = zv_spg
+            END DO
+         END DO
+         !
+         ! Set next velocities:
+         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN 	! Vector form
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
+                            &     + rdtbt * (                      zwx(ji,jj)   &
+                            &                                 + zu_trd(ji,jj)   &
+                            &                                 + zu_frc(ji,jj) ) & 
+                            &   ) * umask(ji,jj,1)
+
+                  va_e(ji,jj) = (                                zvn_e(ji,jj)   &
+                            &     + rdtbt * (                      zwy(ji,jj)   &
+                            &                                 + zv_trd(ji,jj)   &
+                            &                                 + zv_frc(ji,jj) ) &
+                            &   ) * vmask(ji,jj,1)
+               END DO
+            END DO
+
+         ELSE					 	! Flux form
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+
+                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
+                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
+
+                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
+                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
+                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
+                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
+                            &   ) * zhura
+
+                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   &
+                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
+                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
+                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
+                            &   ) * zhvra
+               END DO
+            END DO
+         ENDIF
+         !
+         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
+            !                                          !  ----------------------------------------------        
+            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
+            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
             !
          ENDIF
+         !                                                 !* domain lateral boundary
+         !                                                 !  -----------------------
+         !
+         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries 
+         CALL lbc_lnk( va_e  , 'V', -1._wp )
+
+#if defined key_bdy  
+
+         pssh => ssha_e
+         phur => hur_e
+         phvr => hvr_e
+         pua2d => ua_e
+         pva2d => va_e
+         pub2d => zun_e
+         pvb2d => zvn_e
+                                      
+         IF( lk_bdy )   CALL bdy_dyn2d( kt )               ! open boundaries
+#endif
+#if defined key_agrif
+         IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( kt, jn ) ! Agrif
+#endif
+         !                                             !* Swap
+         !                                             !  ----
+         ubb_e  (:,:) = ub_e  (:,:)
+         ub_e   (:,:) = zun_e (:,:)
+         zun_e  (:,:) = ua_e  (:,:)
+         !
+         vbb_e  (:,:) = vb_e  (:,:)
+         vb_e   (:,:) = zvn_e (:,:)
+         zvn_e  (:,:) = va_e  (:,:)
+         !
+         sshbb_e(:,:) = sshb_e(:,:)
+         sshb_e (:,:) = sshn_e(:,:)
+         sshn_e (:,:) = ssha_e(:,:)
+
+         !                                             !* Sum over whole bt loop
+         !                                             !  ----------------------
+         za1 = wgtbtp1(jn)                                    
+         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities
+            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
+            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
+         ELSE                     			   ! Sum transports
+            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
+            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
+         ENDIF
+         !                       			   ! Sum sea level
+         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
          !                                                 ! ==================== !
       END DO                                               !        end loop      !
       !                                                    ! ==================== !
-
-#if defined key_obc
-      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(:,:)
-      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
-#endif
-
       ! -----------------------------------------------------------------------------
       ! Phase 3. update the general trend with the barotropic trend
       ! -----------------------------------------------------------------------------
       !
-      !                                   !* Time average ==> after barotropic u, v, ssh
-      zcoef =  1._wp / ( 2 * nn_baro  + 1 ) 
-      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
-      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
-      ! 
-      !                                   !* update the general momentum trend
-      DO jk=1,jpkm1
-         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b
-         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b
+      ! At this stage ssha holds a time averaged value
+      !                                                ! Sea Surface Height at u-,v- and f-points
+      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1      ! NO Vector Opt.
+               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
+                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
+                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
+               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
+                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
+                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions
+      ENDIF
+      !
+      ! Set advection velocity correction:
+      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
+         un_adv(:,:) = zu_sum(:,:)*hur(:,:)
+         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
+      ELSE
+         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
+         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
+      END IF
+
+      IF (ln_bt_fw) THEN ! Save integrated transport for next computation
+         ub2_b(:,:) = zu_sum(:,:)
+         vb2_b(:,:) = zv_sum(:,:)
+      ENDIF
+      !
+      ! Update barotropic trend:
+      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
+         DO jk=1,jpkm1
+            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
+            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
+         END DO
+      ELSE
+         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
+         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
+         DO jk=1,jpkm1
+            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b
+            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b
+         END DO
+         ! Save barotropic velocities not transport:
+         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
+         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
+      ENDIF
+      !
+      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)
+         !
       END DO
-      un_b  (:,:) =  zu_sum(:,:) 
-      vn_b  (:,:) =  zv_sum(:,:) 
-      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
       !
       !                                   !* write time-spliting arrays in the restart
-      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
-      !
-      CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     )
-      CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
-      CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
+      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
+      !
+      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
+      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
+      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
+      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
+      CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b                                     )
+      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
+      CALL wrk_dealloc( jpi, jpj, zht, zhf )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
@@ -688,4 +924,80 @@
    END SUBROUTINE dyn_spg_ts
 
+   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ts_wgt  ***
+      !!
+      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
+      !!----------------------------------------------------------------------
+      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
+      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
+      INTEGER, INTENT(inout) :: jpit      ! cycle length    
+      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
+                                                         zwgt2    ! Secondary weights
+      
+      INTEGER ::  jic, jn, ji                      ! temporary integers
+      REAL(wp) :: za1, za2
+      !!----------------------------------------------------------------------
+
+      zwgt1(:) = 0._wp
+      zwgt2(:) = 0._wp
+
+      ! Set time index when averaged value is requested
+      IF (ll_fw) THEN 
+         jic = nn_baro
+      ELSE
+         jic = 2 * nn_baro
+      ENDIF
+
+      ! Set primary weights:
+      IF (ll_av) THEN
+           ! Define simple boxcar window for primary weights 
+           ! (width = nn_baro, centered around jic)     
+         SELECT CASE ( nn_bt_flt )
+              CASE( 0 )  ! No averaging
+                 zwgt1(jic) = 1._wp
+                 jpit = jic
+
+              CASE( 1 )  ! Boxcar, width = nn_baro
+                 DO jn = 1, 3*nn_baro
+                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+                    IF (za1 < 0.5_wp) THEN
+                      zwgt1(jn) = 1._wp
+                      jpit = jn
+                    ENDIF
+                 ENDDO
+
+              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
+                 DO jn = 1, 3*nn_baro
+                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+                    IF (za1 < 1._wp) THEN
+                      zwgt1(jn) = 1._wp
+                      jpit = jn
+                    ENDIF
+                 ENDDO
+              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
+         END SELECT
+
+      ELSE ! No time averaging
+         zwgt1(jic) = 1._wp
+         jpit = jic
+      ENDIF
+    
+      ! Set secondary weights
+      DO jn = 1, jpit
+        DO ji = jn, jpit
+             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
+        END DO
+      END DO
+
+      ! Normalize weigths:
+      za1 = 1._wp / SUM(zwgt1(1:jpit))
+      za2 = 1._wp / SUM(zwgt2(1:jpit))
+      DO jn = 1, jpit
+        zwgt1(jn) = zwgt1(jn) * za1
+        zwgt2(jn) = zwgt2(jn) * za2
+      END DO
+      !
+   END SUBROUTINE ts_wgt
 
    SUBROUTINE ts_rst( kt, cdrw )
@@ -698,81 +1010,141 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
       !
-      INTEGER ::  ji, jk        ! dummy loop indices
       !!----------------------------------------------------------------------
       !
       IF( TRIM(cdrw) == 'READ' ) THEN
-         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
+         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
+         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
+         IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN
+            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
+            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
+            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
          ELSE
-            un_b (:,:) = 0._wp
-            vn_b (:,:) = 0._wp
-            ! vertical sum
-            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-               DO jk = 1, jpkm1
-                  DO ji = 1, jpij
-                     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_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
-
-         ! Vertically integrated velocity (before)
-         IF (neuler/=0) THEN
-            ub_b (:,:) = 0._wp
-            vb_b (:,:) = 0._wp
-
-            ! vertical sum
-            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-               DO jk = 1, jpkm1
-                  DO ji = 1, jpij
-                     ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
-                     vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
-                  END DO
-               END DO
-            ELSE                             ! No  vector opt.
-               DO jk = 1, jpkm1
-                  ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
-                  vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
-               END DO
-            ENDIF
-
-            IF( lk_vvl ) THEN
-               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
-               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
-            ELSE
-               ub_b(:,:) = ub_b(:,:) * hur(:,:)
-               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-            ENDIF
-         ELSE                                 ! neuler==0
-            ub_b (:,:) = un_b (:,:)
-            vb_b (:,:) = vn_b (:,:)
-         ENDIF
-
-         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered ssh
-         ELSE
-            sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
-         ENDIF 
+            sshbb_e = sshn_b                                                ! ACC GUESS WORK
+            ubb_e   = ub_b
+            vbb_e   = vb_b
+            sshb_e  = sshn_b
+            ub_e    = ub_b
+            vb_e    = vb_b
+         ENDIF
+      !
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
-         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
-         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   ! 
+         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
+         !
+         IF (.NOT.ln_bt_av) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
+            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
+         ENDIF
       ENDIF
       !
    END SUBROUTINE ts_rst
 
+   SUBROUTINE dyn_spg_ts_init( kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dyn_spg_ts_init  ***
+      !!
+      !! ** Purpose : Set time splitting options
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      !
+      INTEGER  :: ji ,jj, jk
+      REAL(wp) :: zxr2, zyr2, zcmax
+      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zht
+      !!
+!      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
+!      &                  nn_baro, rn_bt_cmax, nn_bt_flt
+      !!----------------------------------------------------------------------
+!      REWIND( numnam )              !* Namelist namsplit: split-explicit free surface
+!      READ  ( numnam, namsplit )
+      !         ! Max courant number for ext. grav. waves
+      !
+      CALL wrk_alloc( jpi, jpj, zcu, zht )
+      !
+      ! JC: Simplification needed below: define ht_0 even when volume is fixed
+      IF (lk_vvl) THEN 
+         zht(:,:) = ht_0(:,:) * tmask(:,:,1)
+      ELSE
+         zht(:,:) = 0.
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+      ENDIF
+
+      DO jj = 1, jpj
+         DO ji =1, jpi
+            zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
+            zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
+            zcu(ji,jj) = sqrt(grav*zht(ji,jj)*(zxr2 + zyr2) )
+         END DO
+      END DO
+
+      zcmax = MAXVAL(zcu(:,:))
+      IF( lk_mpp )   CALL mpp_max( zcmax )
+
+      ! Estimate number of iterations to satisfy a max courant number=0.8 
+      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
+      
+      rdtbt = rdt / FLOAT(nn_baro)
+      zcmax = zcmax * rdtbt
+							! Print results
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+      IF( ln_bt_nn_auto ) THEN
+         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
+         IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
+      ELSE
+         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
+      ENDIF
+      IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
+      IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
+      IF(lwp) WRITE(numout,*) ' Maximum Courant number is   :', zcmax
+
+      IF(ln_bt_av) THEN
+         IF(lwp) WRITE(numout,*) ' ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
+      ELSE
+         IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
+      ENDIF
+      !
+      !
+      IF(ln_bt_fw) THEN
+         IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true.  => Forward integration of barotropic variables '
+      ELSE
+         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables '
+      ENDIF
+      !
+      IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
+      SELECT CASE ( nn_bt_flt )
+           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '      Dirac'
+           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = nn_baro'
+           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = 2*nn_baro' 
+           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
+      END SELECT
+      !
+      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
+         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
+      ENDIF
+      IF ( zcmax>0.9_wp ) THEN
+         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )          
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, zcu, zht )
+      !
+   END SUBROUTINE dyn_spg_ts_init
+
 #else
-   !!----------------------------------------------------------------------
-   !!   Default case :   Empty module   No standart free surface cst volume
-   !!----------------------------------------------------------------------
+   !!---------------------------------------------------------------------------
+   !!   Default case :   Empty module   No standard free surface constant volume
+   !!---------------------------------------------------------------------------
+
+   USE par_kind
+   LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping
 CONTAINS
    INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
@@ -787,5 +1159,9 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
       WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
-   END SUBROUTINE ts_rst    
+   END SUBROUTINE ts_rst  
+   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt
+   END SUBROUTINE dyn_spg_ts_init
 #endif
    
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 4292)
@@ -572,4 +572,5 @@
       INTEGER  ::   ierr                                          ! local integer
       REAL(wp) ::   zfac12, zua, zva                              ! local scalars
+      REAL(wp) ::   zmsk, ze3                                     ! local scalars
       !                                                           !  3D workspace 
       REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
@@ -577,6 +578,5 @@
 #if defined key_vvl
       REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
-#endif
-#if ! defined key_vvl
+#else
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
 #endif
@@ -604,11 +604,13 @@
       ENDIF
 
-      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
+      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
          DO jk = 1, jpk
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-                  ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25
-                  IF( ze3f(ji,jj,jk) /= 0._wp )   ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk)
+                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
+                  zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
+                     &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
+                  IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
                END DO
             END DO
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4292)
@@ -16,4 +16,5 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE domvvl          ! variable volume
    USE sbc_oce         ! surface boundary condition: ocean
    USE zdf_oce         ! ocean vertical physics
@@ -24,4 +25,7 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE dynadv          ! dynamics: vector invariant versus flux form
+   USE dynspg_oce, ONLY: lk_dynspg_ts, ua_b, va_b
+   USE dynspg_ts
 
    IMPLICIT NONE
@@ -29,4 +33,6 @@
 
    PUBLIC   dyn_zdf_imp   ! called by step.F90
+
+   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise 
 
    !! * Substitutions
@@ -64,8 +70,8 @@
       INTEGER  ::   ikbu, ikbv   ! local integers
       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
+      REAL(wp) ::   ze3ua, ze3va
       !!----------------------------------------------------------------------
 
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zavmu, zavmv
       !!----------------------------------------------------------------------
       !
@@ -73,5 +79,4 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
-      CALL wrk_alloc( jpi,jpj, zavmu, zavmv ) 
       !
       IF( kt == nit000 ) THEN
@@ -79,4 +84,8 @@
          IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+         !
+         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
+         ELSE                ;    r_vvl = 0._wp       
+         ENDIF
       ENDIF
 
@@ -94,19 +103,54 @@
       IF( ln_bfrimp ) THEN
 # if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+         DO jj = 1, 1
+            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
 # endif
-            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            zavmu(ji,jj) = avmu(ji,jj,ikbu+1)
-            zavmv(ji,jj) = avmv(ji,jj,ikbv+1)
-            avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
-            avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
-         END DO
-      END DO
-      ENDIF
+               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
+               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
+               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
+               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
+            END DO
+         END DO
+      ENDIF
+
+#if defined key_dynspg_ts
+      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
+            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
+         END DO
+      ELSE                                            ! applied on thickness weighted velocity
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
+               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
+               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
+            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
+               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
+               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
+         END DO
+      ENDIF
+
+      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
+         ! remove barotropic velocities:
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
+            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
+         ENDDO
+         ! Add bottom stress due to barotropic component only:
+         DO jj = 2, jpjm1        
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
+               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
+               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
+            END DO
+         END DO
+      ENDIF
+#endif
 
       ! 2. Vertical diffusion on u
@@ -119,5 +163,6 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef = - p2dt / fse3u(ji,jj,jk)
+               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
+               zcoef = - p2dt / ze3ua      
                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
@@ -128,5 +173,5 @@
          END DO
       END DO
-      DO jj = 2, jpjm1        ! Surface boudary conditions
+      DO jj = 2, jpjm1        ! Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
@@ -160,6 +205,12 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-               &                                                       * r1_rau0 / fse3u(ji,jj,1)       )
+            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
+#if defined key_dynspg_ts
+            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+               &                                      / ( ze3ua * rau0 ) 
+#else
+            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+               &                                                     / ( fse3u(ji,jj,1) * rau0     ) ) 
+#endif
          END DO
       END DO
@@ -167,5 +218,9 @@
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
+#if defined key_dynspg_ts
+               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
+#else
+               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
+#endif
                ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
             END DO
@@ -186,4 +241,5 @@
       END DO
 
+#if ! defined key_dynspg_ts
       ! Normalization to obtain the general momentum trend ua
       DO jk = 1, jpkm1
@@ -194,5 +250,5 @@
          END DO
       END DO
-
+#endif
 
       ! 3. Vertical diffusion on v
@@ -205,5 +261,6 @@
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef = -p2dt / fse3v(ji,jj,jk)
+               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
+               zcoef = - p2dt / ze3va
                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
@@ -214,5 +271,5 @@
          END DO
       END DO
-      DO jj = 2, jpjm1        ! Surface boudary conditions
+      DO jj = 2, jpjm1        ! Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
@@ -246,6 +303,12 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                                       * r1_rau0 / fse3v(ji,jj,1)       )
+            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
+#if defined key_dynspg_ts            
+            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+               &                                      / ( ze3va * rau0 ) 
+#else
+            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
+#endif
          END DO
       END DO
@@ -253,5 +316,9 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
+#if defined key_dynspg_ts
+               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
+#else
+               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
+#endif
                va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
             END DO
@@ -259,5 +326,5 @@
       END DO
       !
-      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
+      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
@@ -273,4 +340,5 @@
 
       ! Normalization to obtain the general momentum trend va
+#if ! defined key_dynspg_ts
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1   
@@ -280,5 +348,7 @@
          END DO
       END DO
-
+#endif
+
+      ! J. Chanut: Lines below are useless ?
       !! restore bottom layer avmu(v) 
       IF( ln_bfrimp ) THEN
@@ -292,6 +362,6 @@
             ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
             ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            avmu(ji,jj,ikbu+1) = zavmu(ji,jj)
-            avmv(ji,jj,ikbv+1) = zavmv(ji,jj)
+            avmu(ji,jj,ikbu+1) = 0.e0
+            avmv(ji,jj,ikbv+1) = 0.e0
          END DO
       END DO
@@ -299,5 +369,4 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
-      CALL wrk_dealloc( jpi,jpj, zavmu, zavmv) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4292)
@@ -8,9 +8,11 @@
    !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   ssh_wzv        : after ssh & now vertical velocity
-   !!   ssh_nxt        : filter ans swap the ssh arrays
+   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   ssh_nxt        : after ssh
+   !!   ssh_swp        : filter ans swap the ssh arrays
+   !!   wzv            : compute now vertical velocity
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -20,4 +22,5 @@
    USE divcur          ! hor. divergence and curl      (div & cur routines)
    USE iom             ! I/O library
+   USE restart         ! only for lrst_oce
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
@@ -28,7 +31,11 @@
    USE obc_oce
    USE bdy_oce
+   USE bdy_par         
+   USE bdydyn2d        ! bdy_ssh routine
    USE diaar5, ONLY:   lk_diaar5
    USE iom
-   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff 
+   USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div   ! River runoff 
+   USE dynspg_ts,   ONLY: ln_bt_fw
+   USE dynspg_oce, ONLY: lk_dynspg_ts
 #if defined key_agrif
    USE agrif_opa_update
@@ -44,6 +51,7 @@
    PRIVATE
 
-   PUBLIC   ssh_wzv    ! called by step.F90
    PUBLIC   ssh_nxt    ! called by step.F90
+   PUBLIC   wzv        ! called by step.F90
+   PUBLIC   ssh_swp    ! called by step.F90
 
    !! * Substitutions
@@ -57,97 +65,36 @@
 CONTAINS
 
-   SUBROUTINE ssh_wzv( kt ) 
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ssh_wzv  ***
+   SUBROUTINE ssh_nxt( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE ssh_nxt  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
-      !!              and update the now vertical coordinate (lk_vvl=T).
-      !!
-      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
-      !!      velocity is computed by integrating the horizontal divergence  
-      !!      from the bottom to the surface minus the scale factor evolution.
-      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !! ** Purpose :   compute the after ssh (ssha)
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
+      !!      is computed by integrating the horizontal divergence and multiply by
+      !!      by the time step.
       !!
       !! ** action  :   ssha    : after sea surface height
-      !!                wn      : now vertical velocity
-      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
-      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! time step
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv')
-      !
-      CALL wrk_alloc( jpi, jpj, z2d, zhdiv ) 
+      !
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
+      INTEGER, INTENT(in) ::   kt                      ! time step
+      ! 
+      INTEGER             ::   jk                      ! dummy loop indice
+      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
+      !
+      CALL wrk_alloc( jpi, jpj, zhdiv ) 
       !
       IF( kt == nit000 ) THEN
          !
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
+         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
          !
-         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
-         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
-                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
-                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
-                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
-                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
-                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
-                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )
-         ENDIF
-         !
-      ENDIF
-
-      !                                           !------------------------------------------!
-      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
-         !                                        !------------------------------------------!
-         DO jk = 1, jpkm1
-            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
-            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
-            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
-            !
-            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
-            fse3u (:,:,jk) = fse3u_n (:,:,jk)
-            fse3v (:,:,jk) = fse3v_n (:,:,jk)
-            fse3f (:,:,jk) = fse3f_n (:,:,jk)
-            fse3w (:,:,jk) = fse3w_n (:,:,jk)
-            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
-            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
-         END DO
-         !
-         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
-         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
-         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
-         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
-         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
-         ! 
       ENDIF
       !
@@ -162,10 +109,10 @@
       zhdiv(:,:) = 0._wp
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
+        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
       END DO
       !                                                ! Sea surface elevation time stepping
       ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
       ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
-      z1_rau0 = 0.5 / rau0
+      z1_rau0 = 0.5_wp * r1_rau0
       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
 
@@ -180,7 +127,13 @@
 #endif
 #if defined key_bdy
-      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
-      CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm)
-#endif
+      ! bg jchanut tschanges
+      ! These lines are not necessary with time splitting since
+      ! boundary condition on sea level is set during ts loop
+      IF (lk_bdy) THEN
+         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
+         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
+      ENDIF
+#endif
+      ! end jchanut tschanges
 #if defined key_asminc
       !                                                ! Include the IAU weighted SSH increment
@@ -190,43 +143,104 @@
       ENDIF
 #endif
-      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
-      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
-                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
-               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
-                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+
+      !                                           !------------------------------!
+      !                                           !           outputs            !
+      !                                           !------------------------------!
+      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
+      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
+      !
+      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
+      !
+   END SUBROUTINE ssh_nxt
+
+   
+   SUBROUTINE wzv( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE wzv  ***
+      !!                   
+      !! ** Purpose :   compute the now vertical velocity
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
+      !!      velocity is computed by integrating the horizontal divergence  
+      !!      from the bottom to the surface minus the scale factor evolution.
+      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !!
+      !! ** action  :   wn      : now vertical velocity
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt           ! time step
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
+      !
+      INTEGER             ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp)            ::   z1_2dt       ! local scalars
+      !!----------------------------------------------------------------------
+      
+      IF( nn_timing == 1 )  CALL timing_start('wzv')
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
+         IF(lwp) WRITE(numout,*) '~~~~~ '
+         !
+         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         !
+      ENDIF
+      !                                           !------------------------------!
+      !                                           !     Now Vertical Velocity    !
+      !                                           !------------------------------!
+      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
+      !
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
+         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
+         !
+         DO jk = 1, jpkm1
+            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
+            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+               END DO
             END DO
          END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
-      ENDIF
-
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-      z1_2dt = 1.e0 / z2dt
-      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
-         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-            &                         * tmask(:,:,jk) * z1_2dt
+         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
+         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
+         !                             ! Same question holds for hdivn. Perhaps just for security
+         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+            ! computation of w
+            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
+               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+         END DO
+         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
+         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
+      ELSE   ! z_star and linear free surface cases
+         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+            ! computation of w
+            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
+               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+         END DO
+      ENDIF
+
 #if defined key_bdy
          wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
 #endif
-      END DO
-
+      !
       !                                           !------------------------------!
       !                                           !           outputs            !
       !                                           !------------------------------!
-      CALL iom_put( "woce", wn                    )   ! vertical velocity
-      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
-      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
+      CALL iom_put( "woce", wn )   ! vertical velocity
       IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
+         CALL wrk_alloc( jpi, jpj, z2d ) 
+         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         CALL wrk_alloc( jpi,jpj,jpk, z3d )
-         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
+         z2d(:,:) = rau0 * e12t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
@@ -234,17 +248,14 @@
          CALL iom_put( "w_masstr" , z3d                     )  
          CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-         CALL wrk_dealloc( jpi,jpj,jpk, z3d )
-      ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
-      !
-      CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv')
-      !
-   END SUBROUTINE ssh_wzv
-
-
-   SUBROUTINE ssh_nxt( kt )
+         CALL wrk_dealloc( jpi, jpj, z2d  ) 
+         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('wzv')
+
+
+   END SUBROUTINE wzv
+
+   SUBROUTINE ssh_swp( kt )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ssh_nxt  ***
@@ -252,5 +263,5 @@
       !! ** Purpose :   achieve the sea surface  height time stepping by 
       !!              applying Asselin time filter and swapping the arrays
-      !!              ssha  already computed in ssh_wzv  
+      !!              ssha  already computed in ssh_nxt  
       !!
       !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
@@ -266,86 +277,25 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zec      ! temporary scalar
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
       !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
+         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
 
-      !                       !--------------------------!
-      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
-         !                    !--------------------------!
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
-            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
-            sshu_n(:,:) = sshu_a(:,:)
-            sshv_n(:,:) = sshv_a(:,:)
-            DO jj = 1, jpjm1                                ! ssh now at f-point
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
-                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
-            !
-         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
-            zec = atfp * rdt / rau0
-            DO jj = 1, jpj
-               DO ji = 1, jpi                               ! before <-- now filtered
-                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
-                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
-                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
-                  sshu_n(ji,jj) = sshu_a(ji,jj)
-                  sshv_n(ji,jj) = sshv_a(ji,jj)
-               END DO
-            END DO
-            DO jj = 1, jpjm1                                ! ssh now at f-point
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
-                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
-            !
-            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-               END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
-            !
-         ENDIF
-         !                    !--------------------------!
-      ELSE                    !        fixed levels      !     (ssh at t-point only)
-         !                    !--------------------------!
-         !
-         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
-            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
-            !
-         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
-            DO jj = 1, jpj
-               DO ji = 1, jpi                               ! before <-- now filtered
-                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
-                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
-               END DO
-            END DO
-         ENDIF
-         !
+# if defined key_dynspg_ts
+      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
+# else
+      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
+#endif
+         sshb(:,:) = sshn(:,:)                           ! before <-- now
+         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
+      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
+         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
+         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+         sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
       !
@@ -357,7 +307,7 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
-      !
-   END SUBROUTINE ssh_nxt
+      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
+      !
+   END SUBROUTINE ssh_swp
 
    !!======================================================================
