Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4334)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4338)
@@ -61,5 +61,5 @@
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a                        ! baroclinic scale factors
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence
@@ -82,6 +82,7 @@
       IF( ln_vvl_zstar ) dom_vvl_alloc = 0
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
-            &      un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     , STAT = dom_vvl_alloc        )
+         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
+            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
+            &      STAT = dom_vvl_alloc        )
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
@@ -215,5 +216,5 @@
                END DO
             END DO
-            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
+            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
                ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
                ij0 = 128   ;   ij1 = 135   ;   
@@ -229,5 +230,5 @@
 
 
-   SUBROUTINE dom_vvl_sf_nxt( kt ) 
+   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_sf_nxt  ***
@@ -255,4 +256,5 @@
       !! * Arguments
       INTEGER, INTENT( in )                  :: kt                    ! time step
+      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence
       !! * Local declarations
       INTEGER                                :: ji, jj, jk            ! dummy loop indices
@@ -260,4 +262,5 @@
       REAL(wp)                               :: z2dt                  ! temporary scalars
       REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
+      LOGICAL                                :: ll_do_bclinic         ! temporary logical
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
@@ -271,20 +274,25 @@
       ENDIF
 
+      ll_do_bclinic = .TRUE.
+      IF( PRESENT(kcall) ) THEN
+         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
+      ENDIF
+
       ! ******************************* !
       ! After acale factors at t-points !
       ! ******************************* !
 
-      !                                                ! ----------------- !
-      IF( ln_vvl_zstar ) THEN                          ! z_star coordinate !
-         !                                             ! ----------------- !
-
-         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
-         DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-         END DO
-
-      !                                                ! --------------------------- !
-      ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN   ! z_tilde or layer coordinate !
-         !                                             ! --------------------------- !
+      !                                                ! --------------------------------------------- !
+                                                       ! z_star coordinate and barotropic z-tilde part !
+      !                                                ! --------------------------------------------- !
+
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+      DO jk = 1, jpkm1
+         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
+         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+      END DO
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
+         !                                                            ! ------baroclinic part------ !
 
          ! I - initialization
@@ -423,15 +431,13 @@
          tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
 
-         ! Add "tilda" part to the after scale factor
+         !
+         ! "tilda" change in the after scale factor
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         fse3t_a(:,:,:) = e3t_0(:,:,:) + tilde_e3t_a(:,:,:)
-
+         DO jk = 1, jpkm1
+            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
+         END DO
          ! III - Barotropic repartition of the sea surface height over the baroclinic profile
          ! ==================================================================================
-         ! add e3t(n-1) "star" Asselin-filtered
-         DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - tilde_e3t_b(:,:,jk)
-         END DO
-         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n)
+         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
          ! - ML - baroclinicity error should be better treated in the future
          !        i.e. locally and not spread over the water column.
@@ -441,12 +447,19 @@
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
-         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
-         DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-         END DO
-
-      ENDIF
-
-      IF( ln_vvl_dbg ) THEN   ! - ML - test: control prints for debuging
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+         DO jk = 1, jpkm1
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         END DO
+
+      ENDIF
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
+      !                                           ! ---baroclinic part--------- !
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk)
+         END DO
+      ENDIF
+
+      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
          !
          IF( lwp ) WRITE(numout, *) 'kt =', kt
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 4334)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4338)
@@ -252,5 +252,5 @@
             ! ----------------------------------------------------
             IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
-               ! Remove asselin filtering on thicknesses if forward time splitting
+               ! No asselin filtering on thicknesses if forward time splitting
                   fse3t_b(:,:,:) = fse3t_n(:,:,:)
             ELSE
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 4334)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4338)
@@ -110,6 +110,7 @@
       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
+      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
+      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
+      ! 
       z1_rau0 = 0.5_wp * r1_rau0
       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
@@ -197,5 +198,5 @@
          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)
+            ! - ML - note: computation already 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.
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4334)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4338)
@@ -72,4 +72,5 @@
       INTEGER ::   jk       ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
+      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
       !! ---------------------------------------------------------------------
 
@@ -106,11 +107,12 @@
                          CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
       IF( lk_dynspg_ts ) THEN
+          IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
                                   CALL wzv           ( kstp )  ! now cross-level velocity 
           ! In case the time splitting case, update almost all momentum trends here:
           ! Note that the computation of vertical velocity above, hence "after" sea level
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
-                                  CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )                 ! now in situ density for hpg computation
-          IF( ln_zps      )       CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &  ! zps: now hor. derivative
-                &                                          rhd, gru , grv  )     ! of t, s, rd at the last ocean level
+                                  CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
+          IF( ln_zps      )       CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &   ! zps: now hor. derivative
+                &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
 
                                   ua(:,:,:) = 0.e0             ! set dynamics trends to zero
@@ -135,5 +137,5 @@
                                   CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
       ENDIF
-      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
+      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors
                          CALL wzv           ( kstp )  ! now cross-level velocity 
 
