Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4245)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4246)
@@ -901,4 +901,11 @@
       ENDIF
       !
+      DO jk = 1, jpkm1
+         ! Correct velocities:
+         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
+         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
+         !
+      END DO
+      !
       !                                   !* write time-spliting arrays in the restart
       IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4245)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4246)
@@ -50,6 +50,5 @@
 
    PUBLIC   ssh_nxt    ! called by step.F90
-   PUBLIC   wzv_1      ! called by step.F90
-   PUBLIC   wzv_2      ! called by step.F90
+   PUBLIC   wzv        ! called by step.F90
    PUBLIC   ssh_swp    ! called by step.F90
 
@@ -152,7 +151,7 @@
 
    
-   SUBROUTINE wzv_1( kt )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE wzv_1  ***
+   SUBROUTINE wzv( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE wzv  ***
       !!                   
       !! ** Purpose :   compute the now vertical velocity
@@ -169,5 +168,6 @@
       !
       INTEGER, INTENT(in) ::   kt           ! time step
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhdiv
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
       !
       INTEGER             ::   ji, jj, jk   ! dummy loop indices
@@ -175,12 +175,10 @@
       !!----------------------------------------------------------------------
       
-      IF( nn_timing == 1 )  CALL timing_start('wzv_1')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
+      IF( nn_timing == 1 )  CALL timing_start('wzv')
       !
       IF( kt == nit000 ) THEN
          !
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'wzv_1 : now vertical velocity '
+         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
          IF(lwp) WRITE(numout,*) '~~~~~ '
          !
@@ -195,4 +193,6 @@
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
+         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
+         !
          DO jk = 1, jpkm1
             ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
@@ -213,4 +213,5 @@
          END DO
          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
+         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
       ELSE   ! z_star and linear free surface cases
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
@@ -224,119 +225,5 @@
          wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
 #endif
-
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('wzv_1')
-
-
-   END SUBROUTINE wzv_1
-
-   SUBROUTINE wzv_2( kt )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE wzv_2  ***
-      !!                   
-      !! ** Purpose :   compute the now vertical velocity
-      !!
-      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
-      !!      velocity is computed by integrating the horizontal divergence  
-      !!      from the bottom to the surface minus the scale factor evolution.
-      !!        The boundary conditions are w=0 at the bottom (no flux) and.
-      !!
-      !! ** action  :   wn      : now vertical velocity
-      !!
-      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
-      !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt           ! time step
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      !
-      INTEGER             ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp)            ::   z1_2dt       ! local scalars
-      !!----------------------------------------------------------------------
-      
-      IF( nn_timing == 1 )  CALL timing_start('wzv_2')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 
-      !
-      IF( kt == nit000 ) THEN
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'wzv_2 : now vertical velocity '
-         IF(lwp) WRITE(numout,*) '~~~~~ '
-         !
-      ENDIF
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
-
-      IF(  lk_dynspg_ts ) THEN
-         ! At this stage: 
-         ! 1) vertical scale factors have been updated. 
-         ! 2) Time averaged barotropic velocities around step kt+1 (ua_b, va_b) as well as 
-         !   "advective" barotropic velocities (un_adv, vn_adv) at step kt in agreement 
-         !   with continuity equation are available.
-         ! 3) 3d velocities (ua, va) hold ua_b, va_b issued from time splitting as barotropic components.
-         ! Below => Update "now" velocities, divergence, then vertical velocity
-         ! NB: hdivn is NOT updated such that hdivb is not updated also. This means that horizontal
-         !     momentum diffusion is still performed on Instantaneous barotropic arrays. I experencied
-         !     some issues with UBS with the current method. See "divup" comments in development branch to 
-         !     update the divergence fully if necessary (2013/dev_r3867_MERCATOR1_DYN).
-         !
-         DO jk = 1, jpkm1
-            ! Correct velocities:                            
-            un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
-            vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
-            !
-            ! Compute divergence:
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z3d(ji,jj,jk) =   &
-                     (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)       &
-                      + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)  )    &
-                     / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-      !
-         IF( ln_rnf      )   CALL sbc_rnf_div( z3d )      ! runoffs (update divergence)
-      ELSE
-         z3d(:,:,:) = hdivn(:,:,:)
-      ENDIF
-      !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
-         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)
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
-         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
-         !                             ! Same question holds for hdivn. Perhaps just for security
-         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
-            ! computation of w
-            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * z3d(:,:,jk) + zhdiv(:,:,jk)                    &
-               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
-         END DO
-         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
-      ELSE   ! z_star and linear free surface cases
-         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
-            ! computation of w
-            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * z3d(:,:,jk)                                    &
-               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
-         END DO
-      ENDIF
-
-#if defined key_bdy
-         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
-#endif
-
+      !
       !                                           !------------------------------!
       !                                           !           outputs            !
@@ -345,4 +232,5 @@
       IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
          CALL wrk_alloc( jpi, jpj, z2d ) 
+         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
          z2d(:,:) = rau0 * e12t(:,:)
@@ -353,12 +241,11 @@
          CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
          CALL wrk_dealloc( jpi, jpj, z2d  ) 
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('wzv_2')
-
-
-   END SUBROUTINE wzv_2
+         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('wzv')
+
+
+   END SUBROUTINE wzv
 
    SUBROUTINE ssh_swp( kt )
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4245)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4246)
@@ -98,7 +98,7 @@
       !  Ocean dynamics : hdiv, rot, ssh, e3, wn
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL ssh_nxt       ( kstp )  ! after ssh
+                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
       IF( lk_dynspg_ts ) THEN
-                                  CALL wzv_1         ( kstp )  ! now cross-level velocity 
+                                  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
@@ -124,4 +124,5 @@
                                   CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
                                   CALL dyn_spg( kstp, indic )  ! surface pressure gradient
+                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
 
                                   ua_sv(:,:,:) = ua(:,:,:)     ! save next velocities (not trends !)
@@ -129,5 +130,5 @@
       ENDIF
       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
-                         CALL wzv_2         ( kstp )  ! now cross-level velocity (original)
+                         CALL wzv           ( kstp )  ! now cross-level velocity (original)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 4245)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 4246)
@@ -11,4 +11,5 @@
    USE ldftra_oce       ! ocean tracer   - trends
    USE ldfdyn_oce       ! ocean dynamics - trends
+   USE divcur           ! hor. divergence and curl      (div & cur routines)
    USE in_out_manager   ! I/O manager
    USE iom              !
