Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 5950)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 5951)
@@ -25,6 +25,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-   USE dynadv          ! dynamics: vector invariant versus flux form
-   USE dynspg_oce, ONLY: lk_dynspg_ts
+   USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form
 
    IMPLICIT NONE
@@ -87,9 +86,24 @@
       ENDIF
 
-      ! 0. Local constant initialization
-      ! --------------------------------
-      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
-
-      ! 1. Apply semi-implicit bottom friction
+      ! 1. Time step dynamics
+      ! ---------------------
+
+      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
+
+      ! 2. Apply semi-implicit bottom friction
       ! --------------------------------------
       ! Only needed for semi-implicit bottom friction setup. The explicit
@@ -97,5 +111,4 @@
       ! column vector of the tri-diagonal matrix equation
       !
-
       IF( ln_bfrimp ) THEN
          DO jj = 2, jpjm1
@@ -119,22 +132,9 @@
       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
+      ! With split-explicit free surface, barotropic stress is treated explicitly
+      ! Update velocities at the bottom.
+      ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 
+      !            not lead to the effective stress seen over the whole barotropic loop. 
+      IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN
          ! remove barotropic velocities:
          DO jk = 1, jpkm1
@@ -166,7 +166,6 @@
          END IF
       ENDIF
-#endif
-
-      ! 2. Vertical diffusion on u
+
+      ! 3. Vertical diffusion on u
       ! ---------------------------
       ! Matrix and second member construction
@@ -209,6 +208,5 @@
       !-----------------------------------------------------------------------
       !
-      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-      DO jk = 2, jpkm1
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -220,13 +218,7 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_dynspg_ts
             ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
             ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
                &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
-#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     ) * umask(ji,jj,1) ) 
-#endif
          END DO
       END DO
@@ -234,9 +226,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-#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
@@ -257,16 +245,5 @@
       END DO
 
-#if ! defined key_dynspg_ts
-      ! Normalization to obtain the general momentum trend ua
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
-            END DO
-         END DO
-      END DO
-#endif
-
-      ! 3. Vertical diffusion on v
+      ! 4. Vertical diffusion on v
       ! ---------------------------
       ! Matrix and second member construction
@@ -309,6 +286,5 @@
       !-----------------------------------------------------------------------
       !
-      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-      DO jk = 2, jpkm1        
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -319,14 +295,8 @@
       !
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_dynspg_ts            
+         DO ji = fs_2, fs_jpim1   ! vector opt.          
             ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
             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
@@ -334,9 +304,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-#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
@@ -357,17 +323,6 @@
       END DO
 
-      ! Normalization to obtain the general momentum trend va
-#if ! defined key_dynspg_ts
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
-            END DO
-         END DO
-      END DO
-#endif
-
       ! J. Chanut: Lines below are useless ?
-      !! restore bottom layer avmu(v) 
+      !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces)
       IF( ln_bfrimp ) THEN
         DO jj = 2, jpjm1
