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')
