Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5866)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6004)
@@ -18,32 +18,33 @@
    !!            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
-   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
+   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
+   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
    !!-------------------------------------------------------------------------
   
    !!-------------------------------------------------------------------------
-   !!   dyn_nxt      : obtain the next (after) horizontal velocity
+   !!   dyn_nxt       : obtain the next (after) horizontal velocity
    !!-------------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE phycst          ! physical constants
-   USE dynspg_oce      ! type of surface pressure gradient
-   USE dynadv          ! dynamics: vector invariant versus flux form
-   USE domvvl          ! variable volume
-   USE bdy_oce         ! ocean open boundary conditions
-   USE bdydta          ! ocean open boundary conditions
-   USE bdydyn          ! ocean open boundary conditions
-   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
-   USE trd_oce         ! trends: ocean variables
-   USE trddyn          ! trend manager: dynamics
-   USE trdken          ! trend manager: kinetic energy
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE phycst         ! physical constants
+   USE dynadv         ! dynamics: vector invariant versus flux form
+   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
+   USE domvvl         ! variable volume
+   USE bdy_oce        ! ocean open boundary conditions
+   USE bdydta         ! ocean open boundary conditions
+   USE bdydyn         ! ocean open boundary conditions
+   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   USE trdken         ! trend manager: kinetic energy
    !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O manager library
-   USE lbclnk          ! lateral boundary condition (or mpp link)
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory Allocation
-   USE prtctl          ! Print control
-   USE timing          ! Timing
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lbclnk         ! lateral boundary condition (or mpp link)
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! Memory Allocation
+   USE prtctl         ! Print control
+   USE timing         ! Timing
 #if defined key_agrif
    USE agrif_opa_interp
@@ -66,15 +67,11 @@
       !!                  ***  ROUTINE dyn_nxt  ***
       !!                   
-      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary 
-      !!             condition on the after velocity, achieved the time stepping 
+      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary 
+      !!             condition on the after velocity, achieve the time stepping 
       !!             by applying the Asselin filter on now fields and swapping 
       !!             the fields.
       !!
-      !! ** Method  : * After velocity is compute using a leap-frog scheme:
-      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
-      !!             Note that with flux form advection and non linear free surface,
-      !!             the leap-frog is applied on thickness weighted velocity.
-      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
-      !!             the time stepping has already been done in dynspg module
+      !! ** Method  : * Ensure after velocities transport matches time splitting
+      !!             estimate (ln_dynspg_ts=T)
       !!
       !!              * Apply lateral boundary conditions on after velocity 
@@ -89,4 +86,5 @@
       !!             Note that with flux form advection and non linear free surface,
       !!             the time filter is applied on thickness weighted velocity.
+      !!             As a result, dyn_nxt MUST be called after tra_nxt.
       !!
       !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
@@ -97,8 +95,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   iku, ikv     ! local integers
-#if ! defined key_dynspg_flt
-      REAL(wp) ::   z2dt         ! temporary scalar
-#endif
-      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
+      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
@@ -108,6 +103,6 @@
       IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
       !
-      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
-      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
+      IF( ln_dynspg_ts       )   CALL wrk_alloc( jpi,jpj,       zue, zve)
+      IF( l_trddyn           )   CALL wrk_alloc( jpi,jpj,jpk,   zua, zva)
       !
       IF( kt == nit000 ) THEN
@@ -117,84 +112,48 @@
       ENDIF
 
-#if defined key_dynspg_flt
-      !
-      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
-      ! -------------
-
-      ! Update after velocity on domain lateral boundaries      (only local domain required)
-      ! --------------------------------------------------
-      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
-      CALL lbc_lnk( va, 'V', -1. ) 
-      !
-#else
-
-# if defined key_dynspg_exp
-      ! Next velocity :   Leap-frog time stepping
-      ! -------------
-      z2dt = 2. * rdt                                 ! Euler or leap-frog time step 
-      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
-      !
-      IF( ln_dynadv_vec .OR. ln_linssh ) THEN         !==  applied on velocity  ==!
+      IF ( ln_dynspg_ts ) THEN
+         ! Ensure below that barotropic velocities match time splitting estimate
+         ! Compute actual transport and replace it with ts estimate at "after" time step
+         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
+         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+         END DO
          DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
-            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
+            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
          END DO
-      ELSE                                            !==  applied on thickness weighted velocity  ==!
-         DO jk = 1, jpkm1
-            ua(:,:,jk) = (          ub(:,:,jk) * e3u_b(:,:,jk)    &
-               &           + z2dt * ua(:,:,jk) * e3u_n(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
-            va(:,:,jk) = (          vb(:,:,jk) * e3v_b(:,:,jk)    &
-               &           + z2dt * va(:,:,jk) * e3v_n(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
-         END DO
-      ENDIF
-# endif
-
-# if defined key_dynspg_ts
-!!gm IF ( lk_dynspg_ts ) THEN ....
-      ! Ensure below that barotropic velocities match time splitting estimate
-      ! Compute actual transport and replace it with ts estimate at "after" time step
-      zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
-      zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
-      DO jk = 2, jpkm1
-         zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
-      END DO
-      DO jk = 1, jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
-      END DO
-
-      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   
-         ! so that asselin contribution is removed at the same time 
-         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
-!!gm ENDIF
-# endif
+         !
+         IF( .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   
+            ! so that asselin contribution is removed at the same time 
+            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
       ! --------------------------------------------------      
-      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
-      CALL lbc_lnk( va, 'V', -1. ) 
-      !
-# if defined key_bdy
-      !                                !* BDY open boundaries
-      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
-      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
-
-!!$   Do we need a call to bdy_vol here??
-      !
-# endif
-      !
 # if defined key_agrif
       CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
 # endif
-#endif
-
+      !
+      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
+      CALL lbc_lnk( va, 'V', -1. ) 
+      !
+# if defined key_bdy
+      !                                !* BDY open boundaries
+      IF( lk_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
+      IF( lk_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
+
+!!$   Do we need a call to bdy_vol here??
+      !
+# endif
+      !
       IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
          z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
@@ -253,22 +212,24 @@
             ! (used as a now filtered scale factor until the swap)
             ! ----------------------------------------------------
-            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
-               ! No asselin filtering on thicknesses if forward time splitting
-                  e3t_b(:,:,:) = e3t_n(:,:,:)
+            IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting
+               e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1)
             ELSE
-               e3t_b(:,:,:) = e3t_n(:,:,:) + atfp * ( e3t_b(:,:,:) - 2._wp * e3t_n(:,:,:) + e3t_a(:,:,:) )
+               DO jk = 1, jpkm1
+                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
+               END DO
                ! Add volume filter correction: compatibility with tracer advection scheme
                ! => time filter + conservation correction (only at the first level)
-               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
-                  e3t_b(:,:,1) = e3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
-                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
+               IF( nn_isf == 0) THEN   ! if no ice shelf melting
+                  zcoef = atfp * rdt * r1_rau0
+                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * (  emp_b(:,:) - emp(:,:)  &
+                     &                                   - rnf_b(:,:) + rnf(:,:)  ) * tmask(:,:,1)
                ELSE                     ! if ice shelf melting
-                  DO jj = 1,jpj
-                     DO ji = 1,jpi
+                  zcoef = atfp * rdt * r1_rau0
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
                         jk = mikt(ji,jj)
-                        e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
-                                          &                      * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
-                                          &                        - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
-                                          &                        + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
+                        e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
+                           &                                         - rnf_b   (ji,jj) + rnf   (ji,jj)  &
+                           &                                         + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,jk)
                      END DO
                   END DO
@@ -276,11 +237,8 @@
             ENDIF
             !
-            IF( ln_dynadv_vec ) THEN
-               ! Before scale factor at (u/v)-points
-               ! -----------------------------------
+            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
+               ! Before filtered scale factor at (u/v)-points
                CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
                CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
-               ! Leap-Frog - Asselin filter and swap: applied on velocity
-               ! -----------------------------------
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
@@ -297,20 +255,19 @@
                END DO
                !
-            ELSE
-               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
-               !------------------------------------------------
+            ELSE                          ! Asselin filter applied on thickness weighted velocity
+               !
+               CALL wrk_alloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
+               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
                CALL dom_vvl_interpol( e3t_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) * e3u_a(ji,jj,jk)
-                        zve3a = va(ji,jj,jk) * e3v_a(ji,jj,jk)
-                        zue3n = un(ji,jj,jk) * e3u_n(ji,jj,jk)
-                        zve3n = vn(ji,jj,jk) * e3v_n(ji,jj,jk)
-                        zue3b = ub(ji,jj,jk) * e3u_b(ji,jj,jk)
-                        zve3b = vb(ji,jj,jk) * e3v_b(ji,jj,jk)
+                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
+                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
+                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
+                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
+                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
                         !
                         zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
@@ -324,11 +281,13 @@
                   END DO
                END DO
-               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
+               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
                e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
+               !
+               CALL wrk_dealloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
             ENDIF
             !
          ENDIF
          !
-         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
+         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
             ! Revert "before" velocities to time split estimate
             ! Doing it here also means that asselin filter contribution is removed  
@@ -364,16 +323,13 @@
       ENDIF
       !
-      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
-      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               un_b(ji,jj) = un_b(ji,jj) + e3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-               vn_b(ji,jj) = vn_b(ji,jj) + e3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
-               !
-               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
-            END DO
-         END DO
+      un_b(:,:) = e3u_a(:,:,jk) * un(:,:,1) * umask(:,:,1)
+      ub_b(:,:) = e3u_b(:,:,jk) * ub(:,:,1) * umask(:,:,1)
+      vn_b(:,:) = e3v_a(:,:,jk) * vn(:,:,1) * vmask(:,:,1)
+      vb_b(:,:) = e3v_b(:,:,jk) * vb(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(ji,jj,jk)
+         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(ji,jj,jk)
+         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(ji,jj,jk)
+         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(ji,jj,jk)
       END DO
       un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
@@ -382,4 +338,8 @@
       vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
       !
+      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
+         CALL iom_put(  "ubar", un_b(:,:) )
+         CALL iom_put(  "vbar", vn_b(:,:) )
+      ENDIF
       IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
          zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
@@ -391,6 +351,6 @@
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
       ! 
-      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
-      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
+      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj,       zue, zve )
+      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk,   zua, zva )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
