Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 2337)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 2338)
@@ -19,4 +19,5 @@
    USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
    USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine)
+   USE dynadv          ! dynamics: vector invariant versus flux form
    USE trdmod         ! ocean dynamics trends
    USE trdmod_oce     ! ocean variables trends
@@ -167,4 +168,10 @@
       ENDIF
 
+      !                        ! Control of momentum forulation
+      IF( lk_dynspg_ts .AND. lk_vvl ) THEN
+         IF( .NOT. ln_dynadv_vec )   &
+           &   CALL ctl_stop( ' Flux formulae not implemented for this free surface formulation ' )
+      ENDIF
+
 #if defined key_obc
       !                        ! Conservation of ocean volume (key_dynspg_flt)
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2337)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2338)
@@ -51,5 +51,7 @@
    REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme)
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   un_b, vn_b   ! averaged velocity
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   un_b, vn_b   ! now    averaged velocity
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ub_b, vb_b   ! before averaged velocity
+
 
    !! * Substitutions
@@ -82,8 +84,8 @@
       !!          momentum and continuity integration. Barotropic former 
       !!          variables are time averaging over the full barotropic cycle
-      !!          (= 2 * baroclinic time step) and saved in zuX_b 
-      !!          and zvX_b (X specifying after, now or before).
+      !!          (= 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 + ( ua_e - sum_k(ub) )
+      !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b )
       !!
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
@@ -109,6 +111,6 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace
       !!
-      REAL(wp), DIMENSION(jpi,jpj) ::   zcu, zwx, zua, zun, zub   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj) ::   zcv, zwy, zva, zvn, zvb   !  -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcu, zwx, zua, zun   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcv, zwy, zva, zvn   !  -      -
       REAL(wp), DIMENSION(jpi,jpj) ::   zun_e, zub_e, zu_sum      ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zvn_e, zvb_e, zv_sum      !  -      -
@@ -161,6 +163,6 @@
       !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
       !                                   ! --------------------------
-      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   zub(:,:) = 0.e0
-      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   zvb(:,:) = 0.e0
+      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0 
+      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0 
       !
       DO jk = 1, jpkm1
@@ -178,7 +180,4 @@
                zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
                zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
-               !                                                                              ! before velocity
-               zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) 
-               zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
             END DO
          END DO
@@ -280,7 +279,7 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj)   &
+               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
                   &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
-               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj)   &
+               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
                   &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
             END DO
@@ -289,6 +288,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj) * hur(ji,jj)
-               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj) * hvr(ji,jj)
+               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
@@ -300,11 +299,4 @@
       zva(:,:) = zva(:,:) * hvr(:,:)
       !
-      IF( lk_vvl ) THEN
-         zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
-         zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
-      ELSE
-         zub(:,:) = zub(:,:) * hur(:,:)
-         zvb(:,:) = zvb(:,:) * hvr(:,:)
-      ENDIF
 
       ! -----------------------------------------------------------------------
@@ -352,5 +344,5 @@
          !                                                !* Update the forcing (OBC, BDY and tides)
          !                                                !  ------------------
-         IF( lk_obc )   CALL obc_dta_bt( kt, jn   )
+         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
          IF( lk_bdy )   CALL bdy_dta_fla( kt, jn+1, icycle )
 
@@ -476,6 +468,6 @@
          !                                                      !         - Correct the velocity
 
-         IF( lk_obc                   )   CALL obc_fla_ts
-         IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla( sshn_e ) 
+         IF( lk_obc               )   CALL obc_fla_ts
+         IF( lk_bdy .OR. ln_tides )   CALL bdy_dyn_fla( sshn_e ) 
          !
          CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries 
@@ -543,13 +535,17 @@
       !                                   !* Time average ==> after barotropic u, v, ssh
       zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
-      un_b  (:,:) = zcoef * zu_sum  (:,:) 
-      vn_b  (:,:) = zcoef * zv_sum  (:,:) 
-      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
+      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
+      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
       ! 
       !                                   !* update the general momentum trend
       DO jk=1,jpkm1
-         ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b
-         va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b
+         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
+         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
       END DO
+      ub_b (:,:) = un_b(:,:)
+      vb_b (:,:) = vn_b(:,:)
+      un_b  (:,:) =  zu_sum(:,:) 
+      vn_b  (:,:) =  zv_sum(:,:) 
+      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
       !
       !                                   !* write time-spliting arrays in the restart
@@ -596,4 +592,37 @@
             vn_b (:,:) = vn_b(:,:) * hvr(:,:)
          ENDIF
+
+         ! Vertically integrated velocity (before)
+         IF (neuler/=0) THEN
+            ub_b (:,:) = 0.e0
+            vb_b (:,:) = 0.e0
+
+            ! 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.e0 - umask(:,:,1) )
+               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - 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 extrenal ssh
