Index: /branches/2013/dev_MERGE_2013/NEMOGCM/ARCH/OLD/arch-ifort_MERCATOR_CLUSTER.fcm
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/ARCH/OLD/arch-ifort_MERCATOR_CLUSTER.fcm	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/ARCH/OLD/arch-ifort_MERCATOR_CLUSTER.fcm	(revision 4370)
@@ -18,5 +18,5 @@
 %NCDF_LIB            -L $NETCDF_LIB -lnetcdf
 %FC                  mpif90
-%FCFLAGS             -assume byterecl  -convert big_endian -i4 -r8 -automatic -align all -O0
+%FCFLAGS             -assume byterecl -convert big_endian -i4 -r8 -O2 -fp-model precise
 %FFLAGS              %FCFLAGS
 %LD                  mpif90
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 4370)
@@ -50,5 +50,4 @@
 !-----------------------------------------------------------------------
    rn_rdt      =  300.     !  time step for the dynamics (and tracer if nn_acc=0)
-   nn_baro     =   30      !  number of barotropic time step            ("key_dynspg_ts")
    rn_rdtmin   =   300.          !  minimum time step on tracers (used if nn_acc=1)
    rn_rdtmax   =   300.          !  maximum time step on tracers (used if nn_acc=1)
@@ -73,4 +72,11 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+   ln_bt_nn_auto =    .FALSE.          !  Set nn_baro automatically to be just below
+                                       !  a user defined maximum courant number (rn_bt_cmax)
+   nn_baro       =    30               !  Number of iterations of barotropic mode
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 4370)
@@ -41,5 +41,4 @@
    nn_msh      =    3      !  create (=1) a mesh file or not (=0)
    rn_rdt      =  400.     !  time step for the dynamics (and tracer if nn_acc=0)
-   nn_baro     =   60      !  number of barotropic time step            ("key_dynspg_ts")
    rn_rdtmin   =  360.           !  minimum time step on tracers (used if nn_acc=1)
    rn_rdtmax   =  360.           !  maximum time step on tracers (used if nn_acc=1)
@@ -64,4 +63,8 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg	(revision 4370)
@@ -41,5 +41,4 @@
    nn_bathy    =    0      !  compute (=0) or read (=1) the bathymetry file
    rn_rdt      = 7200.     !  time step for the dynamics (and tracer if nn_acc=0)
-   nn_baro     =   60      !  number of barotropic time step            ("key_dynspg_ts")
    rn_rdtmin   = 7200.           !  minimum time step on tracers (used if nn_acc=1)
    rn_rdtmax   = 7200.           !  maximum time step on tracers (used if nn_acc=1)
@@ -64,4 +63,8 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 4370)
@@ -47,5 +47,4 @@
    nn_bathy    =    0      !  compute (=0) or read (=1) the bathymetry file
    rn_rdt      = 7200.     !  time step for the dynamics (and tracer if nn_acc=0)
-   nn_baro     =   60      !  number of barotropic time step            ("key_dynspg_ts")
    rn_rdtmin   = 7200.           !  minimum time step on tracers (used if nn_acc=1)
    rn_rdtmax   = 7200.           !  maximum time step on tracers (used if nn_acc=1)
@@ -70,4 +69,8 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 4370)
@@ -37,5 +37,4 @@
    nn_bathy    =    0      !  compute (=0) or read (=1) the bathymetry file
    rn_rdt      = 7200.     !  time step for the dynamics (and tracer if nn_acc=0)
-   nn_baro     =   60      !  number of barotropic time step            ("key_dynspg_ts")
    rn_rdtmin   = 7200.           !  minimum time step on tracers (used if nn_acc=1)
    rn_rdtmax   = 7200.           !  maximum time step on tracers (used if nn_acc=1)
@@ -60,4 +59,8 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 4370)
@@ -46,4 +46,8 @@
    ppkth2      =  999999.              !
    ppacr2      =  999999.              !
+/
+!-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_cfg	(revision 4370)
@@ -46,4 +46,8 @@
    ppkth2      =  999999.              !
    ppacr2      =  999999.              !
+/
+!-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/namelist_cfg	(revision 4370)
@@ -67,4 +67,8 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist_cfg	(revision 4370)
@@ -46,4 +46,8 @@
    ppkth2      =  999999.              !
    ppacr2      =  999999.              !
+/
+!-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist_cfg	(revision 4370)
@@ -53,4 +53,8 @@
    ppkth2      =  999999.              !
    ppacr2      =  999999.              !
+/
+!-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_SAS_LIM/EXP00/namelist_cfg
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_SAS_LIM/EXP00/namelist_cfg	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_SAS_LIM/EXP00/namelist_cfg	(revision 4370)
@@ -50,4 +50,8 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4370)
@@ -31,4 +31,5 @@
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
    ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
+   nn_euler    =       1   !  = 0 : start with forward time step if ln_rstart=.true.
    nn_rstctl   =       0   !  restart control => activated only if ln_rstart = T
                            !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
@@ -119,5 +120,4 @@
                            !
    rn_rdt      = 5760.     !  time step for the dynamics (and tracer if nn_acc=0)
-   nn_baro     =   64      !  number of barotropic time step            ("key_dynspg_ts")
    rn_atfp     =    0.1    !  asselin time filter parameter
    nn_acc      =    0      !  acceleration of convergence : =1      used, rdt < rdttra(k)
@@ -152,4 +152,19 @@
 /
 !-----------------------------------------------------------------------
+&namsplit      !   time splitting parameters                            ("key_dynspg_ts")
+!-----------------------------------------------------------------------
+   ln_bt_fw      =    .TRUE.           !  Forward integration of barotropic equations
+   ln_bt_av      =    .TRUE.           !  Time filtering of barotropic variables
+   ln_bt_nn_auto =    .TRUE.           !  Set nn_baro automatically to be just below
+                                       !  a user defined maximum courant number (rn_bt_cmax)
+   nn_baro       =    30               !  Number of iterations of barotropic mode
+                                       !  during rn_rdt seconds. Only used if ln_bt_nn_auto=F
+   rn_bt_cmax    =    0.8              !  Maximum courant number allowed if ln_bt_nn_auto=T 
+   nn_bt_flt     =    1                !  Time filter choice
+                                       !  = 0 None
+                                       !  = 1 Boxcar over   nn_baro barotropic steps
+                                       !  = 2 Boxcar over 2*nn_baro     "        "  
+/
+!-----------------------------------------------------------------------
 &namcrs        !   Grid coarsening for dynamics output and/or
                !   passive tracer coarsened online simulations
@@ -724,4 +739,5 @@
    ln_vvl_layer  = .false.          !  full layer vertical coordinate
    ln_vvl_ztilde_as_zstar = .false. !  ztilde vertical coordinate emulating zstar
+   ln_vvl_zstar_at_eqtor = .false.  !  ztilde near the equator
    rn_ahe3       = 0.0e0            !  thickness diffusion coefficient
    rn_rst_e3t    = 30.e0            !  ztilde to zstar restoration timescale [days]
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/cfg.txt	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/cfg.txt	(revision 4370)
@@ -5,7 +5,7 @@
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
 ORCA2_SAS_LIM OPA_SRC SAS_SRC LIM_SRC_2 NST_SRC
-AMM12 OPA_SRC
 C1D_PAPA OPA_SRC
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
+AMM12 OPA_SRC
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 4370)
@@ -61,5 +61,4 @@
       !!
       REAL(wp), POINTER, DIMENSION(:,:) :: pua2d, pva2d     ! after barotropic velocities
-      REAL(wp), POINTER, DIMENSION(:,:) :: phura, phvra     ! after inverse depth at u and v points
 
       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn')
@@ -82,5 +81,5 @@
       !-------------------------------------------------------
 
-      CALL wrk_alloc(jpi,jpj,pua2d,pva2d,phura,phvra) 
+      CALL wrk_alloc(jpi,jpj,pua2d,pva2d) 
 
       !-------------------------------------------------------
@@ -91,27 +90,12 @@
 
       pua2d(:,:) = 0.e0
-      pva2d(:,:) = 0.e0
-      
-      IF (lk_vvl) THEN
-         phura(:,:) = 0.
-         phvra(:,:) = 0.
-         DO jk = 1, jpkm1
-            phura(:,:) = phura(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
-            phvra(:,:) = phvra(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
-            pua2d(:,:) = pua2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-            pva2d(:,:) = pva2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
-         END DO
-         phura(:,:) = umask(:,:,1) / ( phura(:,:) + 1. - umask(:,:,1) )
-         phvra(:,:) = vmask(:,:,1) / ( phvra(:,:) + 1. - vmask(:,:,1) )
-      ELSE
-         phura(:,:) = hur(:,:)
-         phvra(:,:) = hvr(:,:)
-         DO jk = 1, jpkm1
-            pua2d(:,:) = pua2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-            pva2d(:,:) = pva2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
-         END DO
-      ENDIF
-      pua2d(:,:) = pua2d(:,:) * phura(:,:)
-      pva2d(:,:) = pva2d(:,:) * phvra(:,:)
+      pva2d(:,:) = 0.e0      
+      DO jk = 1, jpkm1
+         pua2d(:,:) = pua2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+         pva2d(:,:) = pva2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+      END DO
+
+      pua2d(:,:) = pua2d(:,:) * hur_a(:,:)
+      pva2d(:,:) = pva2d(:,:) * hvr_a(:,:)
 
       DO jk = 1 , jpkm1
@@ -134,5 +118,5 @@
       !-------------------------------------------------------
 
-      IF( ll_dyn2d ) CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, phura, phvra, ssha )
+      IF( ll_dyn2d ) CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, hur_a(:,:), hvr_a(:,:), ssha )
 
       IF( ll_dyn3d ) CALL bdy_dyn3d( kt )
@@ -154,5 +138,5 @@
       END IF
 
-      CALL wrk_dealloc(jpi,jpj,pua2d,pva2d,phura,phvra) 
+      CALL wrk_dealloc(jpi,jpj,pua2d,pva2d) 
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn')
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 4370)
@@ -41,7 +41,16 @@
    REAL(wp), PUBLIC ::   rn_rdtmax       !: maximum time step on tracers
    REAL(wp), PUBLIC ::   rn_rdth         !: depth variation of tracer step
-   INTEGER , PUBLIC ::   nn_baro         !: number of barotropic time steps (key_dynspg_ts)
    INTEGER , PUBLIC ::   nn_closea       !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
+   INTEGER , PUBLIC ::   nn_euler        !: =0 start with forward time step or not (=1)
    LOGICAL , PUBLIC ::   ln_crs          !: Apply grid coarsening to dynamical model output or online passive tracers
+
+   !! Time splitting parameters
+   !! =========================
+   LOGICAL,  PUBLIC :: ln_bt_fw          !: Forward integration of barotropic sub-stepping
+   LOGICAL,  PUBLIC :: ln_bt_av          !: Time averaging of barotropic variables
+   LOGICAL,  PUBLIC :: ln_bt_nn_auto     !: Set number of barotropic iterations automatically
+   INTEGER,  PUBLIC :: nn_bt_flt         !: Filter choice
+   INTEGER,  PUBLIC :: nn_baro           !: Number of barotropic iterations during one baroclinic step (rdt)
+   REAL(wp), PUBLIC :: rn_bt_cmax        !: Maximum allowed courant number (used if ln_bt_nn_auto=T)
 
    !! Horizontal grid parameters for domhgr
@@ -93,5 +102,5 @@
 
    !                                                  !!! associated variables
-   INTEGER , PUBLIC                 ::   neuler  = 0   !: restart euler forward option (0=Euler)
+   INTEGER , PUBLIC                 ::   neuler        !: restart euler forward option (0=Euler)
    REAL(wp), PUBLIC                 ::   atfp1         !: asselin time filter coeff. (atfp1= 1-2*atfp)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdttra  !: vertical profile of tracer time step
@@ -194,14 +203,19 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag
 #endif
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur  , hvr     !: inverse of u and v-points ocean depth (1/m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu   , hv      !: depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ht_0           !: reference depth at t-       points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_0 , hv_0    !: reference depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   re2u_e1u       !: scale factor coeffs at u points (e2u/e1u)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   re1v_e2v       !: scale factor coeffs at v points (e1v/e2v)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12t , r1_e12t !: horizontal cell surface and inverse at t points
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12u , r1_e12u !: horizontal cell surface and inverse at u points
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12v , r1_e12v !: horizontal cell surface and inverse at v points
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   e12f , r1_e12f !: horizontal cell surface and inverse at f points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur  , hvr     !: Now    inverse of u and v-points ocean depth (1/m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu   , hv      !:        depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht             !:        depth at t-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehur_a, ehvr_a !: After  inverse of u and v-points ocean depth (1/m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehu_a , ehv_a  !:        depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehur_b, ehvr_b !: Before inverse of u and v-points ocean depth (1/m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ehu_b , ehv_b  !:        depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0           !: reference depth at t-       points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0    !: reference depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   re2u_e1u       !: scale factor coeffs at u points (e2u/e1u)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   re1v_e2v       !: scale factor coeffs at v points (e1v/e2v)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12t , r1_e12t !: horizontal cell surface and inverse at t points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12u , r1_e12u !: horizontal cell surface and inverse at u points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12v , r1_e12v !: horizontal cell surface and inverse at v points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12f , r1_e12f !: horizontal cell surface and inverse at f points
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
@@ -336,9 +350,13 @@
          &      e3t_b   (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) ,                           &
          &      e3uw_b  (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,                                                 &
-         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk)                       , STAT=ierr(5) )
+         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk),                            &
+         &      ehu_a    (jpi,jpj)    , ehv_a  (jpi,jpj),                                                     &
+         &      ehur_a   (jpi,jpj)    , ehvr_a (jpi,jpj),                                                     &
+         &      ehu_b    (jpi,jpj)    , ehv_b  (jpi,jpj),                                                     &
+         &      ehur_b   (jpi,jpj)    , ehvr_b (jpi,jpj),                                  STAT=ierr(5) )                          
 #endif
          !
       ALLOCATE( hu      (jpi,jpj) , hur     (jpi,jpj) , hu_0(jpi,jpj) , ht_0  (jpi,jpj) ,     &
-         &      hv      (jpi,jpj) , hvr     (jpi,jpj) , hv_0(jpi,jpj) ,                       &
+         &      hv      (jpi,jpj) , hvr     (jpi,jpj) , hv_0(jpi,jpj) , ht    (jpi,jpj) ,     &
          &      re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) ,                                       &
          &      e12t    (jpi,jpj) , r1_e12t (jpi,jpj) ,                                       &
@@ -375,2 +393,3 @@
    !!======================================================================
 END MODULE dom_oce
+
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 4370)
@@ -91,9 +91,12 @@
       IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
       !
-      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
-      hv(:,:) = 0._wp
-      DO jk = 1, jpk
+      !
+      hu(:,:) = 0._wp                          ! Ocean depth at U-points
+      hv(:,:) = 0._wp                          ! Ocean depth at V-points
+      ht(:,:) = 0._wp                          ! Ocean depth at T-points
+      DO jk = 1, jpkm1
          hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
          hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
@@ -124,8 +127,8 @@
       NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
          &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
-         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz
+         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz, nn_euler
       NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
          &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
-         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs,    &
+         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
          &             jphgr_msh, &
          &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
@@ -158,4 +161,5 @@
          WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
          WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
+         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
          WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
          WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
@@ -182,5 +186,10 @@
       nstock = nn_stock
       nwrite = nn_write
-
+      neuler = nn_euler
+      IF ( neuler == 1 .AND. .NOT.ln_rstart ) THEN
+         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
+         CALL ctl_warn( ctmp1 )
+         neuler = 0
+      ENDIF
 
       !                             ! control of output frequency
@@ -240,5 +249,4 @@
          WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
          WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
-         WRITE(numout,*) '      time-splitting: nb of sub time-step   nn_baro   = ', nn_baro
          WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
          WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
@@ -447,3 +455,2 @@
    !!======================================================================
 END MODULE domain
-
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 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4370)
@@ -182,4 +182,15 @@
       END DO
 
+      ! Before depth and Inverse of the local depth of the water column at u- and v- points
+      ! -----------------------------------------------------------------------------------
+      hu_b(:,:) = 0.
+      hv_b(:,:) = 0.
+      DO jk = 1, jpkm1
+         hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+      END DO
+      hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
+      hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
+
       ! Restoring frequencies for z_tilde coordinate
       ! ============================================
@@ -514,4 +525,18 @@
       CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
 
+      ! *********************************** !
+      ! After depths at u- v points         !
+      ! *********************************** !
+
+      hu_a(:,:) = 0._wp                        ! Ocean depth at U-points
+      hv_a(:,:) = 0._wp                        ! Ocean depth at V-points
+      DO jk = 1, jpkm1
+         hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
+         hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
+      END DO
+      !                                        ! Inverse of the local depth
+      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
+      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
+
       CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
       CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
@@ -584,4 +609,5 @@
       ! --------------------------------------
       ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
+      ! - JC - hu_b, hv_b, hur_b, hvr_b also
       CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
       ! Vertical scale factor interpolations
@@ -604,13 +630,17 @@
       ! Local depth and Inverse of the local depth of the water column at u- and v- points
       ! ----------------------------------------------------------------------------------
-      hu(:,:) = 0.
-      hv(:,:) = 0.
-      DO jk = 1, jpk
-         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+      hu (:,:) = hu_a (:,:)
+      hv (:,:) = hv_a (:,:)
+
+      ! Inverse of the local depth
+      hur(:,:) = hur_a(:,:)
+      hvr(:,:) = hvr_a(:,:)
+
+      ! Local depth of the water column at t- points
+      ! --------------------------------------------
+      ht(:,:) = 0.
+      DO jk = 1, jpkm1
+         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
       END DO
-      ! Inverse of the local depth
-      hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) )
-      hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - vmask(:,:,1) )
 
       ! Write outputs
@@ -1322,2 +1352,5 @@
    !!======================================================================
 END MODULE domvvl
+
+
+
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90	(revision 4370)
@@ -34,4 +34,13 @@
 
 #   define  fse3t_m(i,j)     e3t_m(i,j)
+
+#   define  hu_a(i,j)        ehu_a(i,j)
+#   define  hur_a(i,j)       ehur_a(i,j)
+#   define  hv_a(i,j)        ehv_a(i,j)
+#   define  hvr_a(i,j)       ehvr_a(i,j)
+#   define  hu_b(i,j)        ehu_b(i,j)
+#   define  hur_b(i,j)       ehur_b(i,j)
+#   define  hv_b(i,j)        ehv_b(i,j)
+#   define  hvr_b(i,j)       ehvr_b(i,j)
 
 ! This part should be removed one day ...
@@ -75,4 +84,13 @@
 #   define  fse3t_m(i,j)     e3t_0(i,j,1)
 
+#   define  hu_a(i,j)        hu(i,j)
+#   define  hur_a(i,j)       hur(i,j)
+#   define  hv_a(i,j)        hv(i,j)
+#   define  hvr_a(i,j)       hvr(i,j)
+#   define  hu_b(i,j)        hu(i,j)
+#   define  hur_b(i,j)       hur(i,j)
+#   define  hv_b(i,j)        hv(i,j)
+#   define  hvr_b(i,j)       hvr(i,j)
+
 ! This part should be removed one day ...
 ! ... In that case all occurence of the above statement functions
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 4370)
@@ -40,7 +40,5 @@
    USE domvvl          ! varying vertical mesh
    USE dynspg_oce      ! pressure gradient schemes
-   USE dynspg_flt      ! pressure gradient schemes
-   USE dynspg_exp      ! pressure gradient schemes
-   USE dynspg_ts       ! pressure gradient schemes
+   USE dynspg_flt      ! filtered free surface
    USE sol_oce         ! ocean solver variables
    USE lib_mpp         ! MPP library
@@ -72,5 +70,4 @@
       ! - ML - needed for initialization of e3t_b
       INTEGER  ::  ji,jj,jk     ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zhur_b, zhvr_b ! U & Inverse of before depths
       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::  zuvd    ! U & V data workspace
       !!----------------------------------------------------------------------
@@ -93,5 +90,4 @@
       IF( ln_rstart ) THEN                    ! Restart from a file
          !                                    ! -------------------
-         neuler = 1                              ! Set time-step indicator at nit000 (leap-frog)
          CALL rst_read                           ! Read the restart file
          CALL day_init                           ! model calendar (using both namelist and restart infos)
@@ -161,15 +157,4 @@
       ! being eventually used
       !
-      IF (lk_vvl) THEN 
-         CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b  )
-         zhur_b(:,:) = 0._wp
-         zhvr_b(:,:) = 0._wp
-         DO jk = 1, jpk
-            zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
-            zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
-         END DO
-         zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
-         zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
-      ENDIF
       !
       un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
@@ -187,22 +172,16 @@
                vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
                !
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      un_b(:,:) = un_b(:,:) * hur(:,:)
-      vn_b(:,:) = vn_b(:,:) * hvr(:,:)
-      !
-      IF( lk_vvl ) THEN
-         ub_b(:,:) = ub_b(:,:) * zhur_b(:,:)
-         vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:)
-      ELSE
-         ub_b(:,:) = ub_b(:,:) * hur(:,:)
-         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-      ENDIF
-      !
-      IF (lk_vvl) CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b  ) 
+      un_b(:,:) = un_b(:,:) * hur  (:,:)
+      vn_b(:,:) = vn_b(:,:) * hvr  (:,:)
+      !
+      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
+      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
+      !
       !
       IF( nn_timing == 1 )  CALL timing_stop('istate_init')
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 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4370)
@@ -39,5 +39,4 @@
    USE wrk_nemo        ! Memory Allocation
    USE prtctl          ! Print control
-   USE dynspg_ts       ! Barotropic velocities
 
 #if defined key_agrif
@@ -102,5 +101,5 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva, zhura, zhvra, zhurb, zhvrb
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
       !!----------------------------------------------------------------------
@@ -109,5 +108,5 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva, zhura, zhvra )
+      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
       !
       IF( kt == nit000 ) THEN
@@ -157,29 +156,12 @@
       zua(:,:) = 0._wp
       zva(:,:) = 0._wp
-      IF (lk_vvl) THEN
-         zhura(:,:) = 0._wp
-         zhvra(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
-            zhura(:,:) = zhura(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
-            zhvra(:,:) = zhvra(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)  
-         END DO
-         zhura(:,:) = umask(:,:,1) / ( zhura(:,:) + 1._wp - umask(:,:,1) ) 
-         zhvra(:,:) = vmask(:,:,1) / ( zhvra(:,:) + 1._wp - vmask(:,:,1) )  
-         DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * zhura(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * zhvra(:,:) + 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
+      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_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+      END DO
 
       IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
@@ -321,15 +303,8 @@
             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
+               zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
+               zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+            END DO
             DO jk = 1, jpkm1
                ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
@@ -344,14 +319,14 @@
       ! integration
       !
-      IF (lk_vvl) THEN 
-         CALL wrk_alloc( jpi, jpj, zhurb, zhvrb  )
-         zhurb(:,:) = 0._wp
-         zhvrb(:,:) = 0._wp
-         DO jk = 1, jpk
-            zhurb(:,:) = zhurb(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
-            zhvrb(:,:) = zhvrb(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
-         END DO
-         zhurb(:,:) = umask(:,:,1) / ( zhurb(:,:) + 1._wp - umask(:,:,1) )
-         zhvrb(:,:) = vmask(:,:,1) / ( zhvrb(:,:) + 1._wp - vmask(:,:,1) )
+      !
+      IF (lk_vvl) THEN
+         hu_b(:,:) = 0.
+         hv_b(:,:) = 0.
+         DO jk = 1, jpkm1
+            hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+         END DO
+         hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
+         hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
       ENDIF
       !
@@ -367,25 +342,19 @@
             DO ji = 1, jpi
 #endif                  
-               un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
                !
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      un_b(:,:) = un_b(:,:) * hur(:,:)
-      vn_b(:,:) = vn_b(:,:) * hvr(:,:)
-      !
-      IF( lk_vvl ) THEN
-         ub_b(:,:) = ub_b(:,:) * zhurb(:,:)
-         vb_b(:,:) = vb_b(:,:) * zhvrb(:,:)
-      ELSE
-         ub_b(:,:) = ub_b(:,:) * hur(:,:)
-         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-      ENDIF
-      !
-      IF (lk_vvl) CALL wrk_dealloc( jpi, jpj, zhurb, zhvrb  ) 
+      !
+      un_b(:,:) = un_b(:,:) * hur_a(:,:)
+      vn_b(:,:) = vn_b(:,:) * hvr_a(:,:)
+      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
+      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
+      !
       !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
@@ -393,5 +362,5 @@
       ! 
       CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva, zhura, zhvra )
+      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 4370)
@@ -34,11 +34,10 @@
 
   !                                                                         !!! Time splitting scheme (key_dynspg_ts) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_e, ssha_e   ! sea surface heigth (now, after, average)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   ua_e  , va_e     ! barotropic velocities (after)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_e  , hv_e     ! now ocean depth ( = Ho+sshn_e )
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur_e , hvr_e    ! inverse of hu_e and hv_e
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_b           ! before field without time-filter
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv ! Advection vel. at "now" barocl. step
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b,  vb2_b  ! Advection vel. at "now-0.5" barocl. step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshn_e, ssha_e   ! sea surface heigth (now, after)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_e  , va_e     ! barotropic velocities (after)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_e  , hv_e     ! now ocean depth ( = Ho+sshn_e )
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur_e , hvr_e    ! inverse of hu_e and hv_e
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv   ! Advection vel. at "now" barocl. step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b,  vb2_b    ! Advection vel. at "now-0.5" barocl. step
 
    !!----------------------------------------------------------------------
@@ -56,6 +55,5 @@
          &      ssha_e(jpi,jpj) , va_e(jpi,jpj) , hv_e(jpi,jpj) , hvr_e(jpi,jpj) ,      &
          &      ub2_b(jpi,jpj)  , vb2_b(jpi,jpj)                                 ,      &
-         &      un_adv(jpi,jpj) , vn_adv(jpi,jpj)                                ,      &
-         &      sshn_b(jpi,jpj)                                                  , STAT = dynspg_oce_alloc )
+         &      un_adv(jpi,jpj) , vn_adv(jpi,jpj)                                , STAT = dynspg_oce_alloc )
          !
       IF( lk_mpp                )   CALL mpp_sum ( dynspg_oce_alloc )
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4370)
@@ -14,5 +14,5 @@
 #if defined key_dynspg_ts   ||   defined key_esopa
    !!----------------------------------------------------------------------
-   !!   'key_dynspg_ts'         free surface cst volume with time splitting
+   !!   'key_dynspg_ts'         split explicit free surface
    !!----------------------------------------------------------------------
    !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
@@ -55,12 +55,4 @@
    PUBLIC dyn_spg_ts_init   !    "      "     "    "
 
-   ! Potential namelist parameters below to be read in dyn_spg_ts_init
-   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.TRUE.        !: Forward integration of barotropic sub-stepping
-   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_av=.TRUE.        !: Time averaging of barotropic variables
-   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE.  !: Set number of iterations automatically
-   INTEGER,  PRIVATE, PARAMETER :: nn_bt_flt=1            !: Filter choice
-   REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp      !: Max. courant number (used if ln_bt_nn_auto=T)
-   ! End namelist parameters
-
    INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
    REAL(wp),SAVE :: rdtbt   ! Barotropic time step
@@ -160,7 +152,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv
       REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
-      REAL(wp), POINTER, DIMENSION(:,:) :: zhur_b, zhvr_b
       REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
-      REAL(wp), POINTER, DIMENSION(:,:) :: zht, zhf
+      REAL(wp), POINTER, DIMENSION(:,:) :: zhf
       !!----------------------------------------------------------------------
       !
@@ -172,7 +163,6 @@
       CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
       CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
-      CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b                                     )
       CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   )
-      CALL wrk_alloc( jpi, jpj, zht, zhf )
+      CALL wrk_alloc( jpi, jpj, zhf )
       !
       !                                         !* Local constant initialization
@@ -228,18 +218,8 @@
       IF ( kt == nit000 .OR. lk_vvl ) THEN
          IF ( ln_dynvor_een ) THEN
-            ! JC: Simplification needed below: define ht_0 even when volume is fixed
-            IF (lk_vvl) THEN
-               zht(:,:) = (ht_0(:,:) + sshn(:,:)) * tmask(:,:,1) 
-            ELSE
-               zht(:,:) = 0.
-               DO jk = 1, jpkm1
-                  zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
-               END DO
-            ENDIF
-
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-                  zwz(ji,jj) =   ( zht(ji  ,jj+1) + zht(ji+1,jj+1) +                     &
-                        &          zht(ji  ,jj  ) + zht(ji+1,jj  )   )                   &
+                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     &
+                        &          ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   &
                         &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
                         &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
@@ -261,28 +241,28 @@
          ELSE
             zwz(:,:) = 0._wp
-            zht(:,:) = 0.
+            zhf(:,:) = 0.
             IF ( .not. ln_sco ) THEN
 !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
 !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
 !              ENDIF
-!              zht(:,:) = gdepw_0(:,:,jk+1)
+!              zhf(:,:) = gdepw_0(:,:,jk+1)
             ELSE
-               zht(:,:) = hbatf(:,:)
+               zhf(:,:) = hbatf(:,:)
             END IF
 
             DO jj = 1, jpjm1
-               zht(:,jj) = zht(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+               zhf(:,jj) = zhf(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
             END DO
 
             DO jk = 1, jpkm1
                DO jj = 1, jpjm1
-                  zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
-               END DO
-            END DO
-            CALL lbc_lnk( zht, 'F', 1._wp )
+                  zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
+               END DO
+            END DO
+            CALL lbc_lnk( zhf, 'F', 1._wp )
             ! JC: TBC. hf should be greater than 0 
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF( zht(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array
+                  IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
                END DO
             END DO
@@ -296,19 +276,4 @@
          ll_fw_start=.FALSE.
          CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
-      ENDIF
-
-      ! before inverse water column height at u- and v- points
-      IF( lk_vvl ) THEN
-         zhur_b(:,:) = 0.
-         zhvr_b(:,:) = 0.
-         DO jk = 1, jpk
-            zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
-            zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
-         END DO
-         zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) )
-         zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) )
-      ELSE
-         zhur_b(:,:) = hur(:,:)
-         zhvr_b(:,:) = hvr(:,:)
       ENDIF
                           
@@ -331,6 +296,6 @@
             DO ji = 1, jpi
 #endif                                                                   
-               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
-               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)         
+               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
+               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)         
             END DO
          END DO
@@ -492,32 +457,26 @@
       !                                             !    Initialisations   !
       !                                             ! ==================== !  
-      ! Initialize barotropic variables:    
-      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                             
-         sshn_e (:,:) = sshn (:,:)            
-         zun_e  (:,:) = un_b (:,:)            
-         zvn_e  (:,:) = vn_b (:,:)
-      ELSE                                ! CENTRED integration: start from BEFORE fields
-         sshn_e (:,:) = sshb (:,:)
-         zun_e  (:,:) = ub_b (:,:)         
-         zvn_e  (:,:) = vb_b (:,:)
-      ENDIF
-      !
-      ! Initialize depths:
-      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN
-         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
-         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
-         hur_e (:,:) = zhur_b(:,:)
-         hvr_e (:,:) = zhvr_b(:,:)
-      ELSE
+      ! Initialize barotropic variables:      
+      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                    
+         sshn_e(:,:) = sshn (:,:)            
+         zun_e (:,:) = un_b (:,:)            
+         zvn_e (:,:) = vn_b (:,:)
+         !
          hu_e  (:,:) = hu   (:,:)       
          hv_e  (:,:) = hv   (:,:) 
          hur_e (:,:) = hur  (:,:)    
          hvr_e (:,:) = hvr  (:,:)
-      ENDIF
-      !
-      IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5:
-         zhup2_e (:,:) = hu(:,:)
-         zhvp2_e (:,:) = hv(:,:)
-      ENDIF
+      ELSE                                ! CENTRED integration: start from BEFORE fields
+         sshn_e(:,:) = sshb (:,:)
+         zun_e (:,:) = ub_b (:,:)         
+         zvn_e (:,:) = vb_b (:,:)
+         !
+         hu_e  (:,:) = hu_b (:,:)       
+         hv_e  (:,:) = hv_b (:,:) 
+         hur_e (:,:) = hur_b(:,:)    
+         hvr_e (:,:) = hvr_b(:,:)
+      ENDIF
+      !
+      !
       !
       ! Initialize sums:
@@ -560,10 +519,10 @@
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
                DO ji = 2, fs_jpim1   ! Vector opt.
-                  zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
-                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
-                     &              +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
-                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
-                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
-                     &              +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)     &
+                     &              * ( e12t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
+                     &              +   e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
+                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)     &
+                     &              * ( e12t(ji,jj  ) * zsshp2_e(ji,jj  )  &
+                     &              +   e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
                END DO
             END DO
@@ -572,4 +531,7 @@
             zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)               ! Ocean depth at U- and V-points
             zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
+         ELSE
+            zhup2_e (:,:) = hu(:,:)
+            zhvp2_e (:,:) = hv(:,:)
          ENDIF
          !                                                !* after ssh
@@ -583,6 +545,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
-                  &             + zwy(ji,jj) - zwy(ji,jj-1)   &
-                  &           ) / ( e1t(ji,jj) * e2t(ji,jj) )
+                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e12t(ji,jj)
             END DO
          END DO
@@ -609,10 +570,10 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1      ! NO Vector Opt.
-                  zsshu_a(ji,jj) = z1_2  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                 &
-                     &                                   * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha_e(ji  ,jj) &
-                     &                                     + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) )
-                  zsshv_a(ji,jj) = z1_2  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                 &
-                     &                                   * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha_e(ji,jj  ) &
-                     &                                     + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) )
+                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)  &
+                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
+                     &              +   e12t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
+                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)  &
+                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
+                     &              +   e12t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
                END DO
             END DO
@@ -648,10 +609,10 @@
             DO jj = 2, jpjm1                            
                DO ji = 2, jpim1
-                  zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
-                    &        * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
-                    &        +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )                 
-                  zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
-                     &       * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
-                     &       +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
+                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    &
+                     &      * ( e12t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
+                     &      +   e12t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
+                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e12v(ji  ,jj  )  &
+                     &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
+                     &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
                   zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
                   zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
@@ -833,10 +794,10 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1      ! NO Vector Opt.
-               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
-                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
-               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
-                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj) &
+                  &              * ( e12t(ji  ,jj) * ssha_e(ji  ,jj)    &
+                  &              +   e12t(ji+1,jj) * ssha_e(ji+1,jj) )
+               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj) &
+                  &              * ( e12t(ji,jj  ) * ssha_e(ji,jj  )    &
+                  &              +   e12t(ji,jj+1) * ssha_e(ji,jj+1) )
             END DO
          END DO
@@ -865,9 +826,7 @@
          END DO
       ELSE
-         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
-         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
          DO jk=1,jpkm1
-            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b
-            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b
+            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
+            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
          END DO
          ! Save barotropic velocities not transport:
@@ -890,7 +849,6 @@
       CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
       CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
-      CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b                                     )
       CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
-      CALL wrk_dealloc( jpi, jpj, zht, zhf )
+      CALL wrk_dealloc( jpi, jpj, zhf )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
@@ -989,5 +947,5 @@
          CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
          CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
-         IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN
+         IF( .NOT.ln_bt_av ) THEN
             CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
             CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
@@ -996,11 +954,4 @@
             CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
             CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
-         ELSE
-            sshbb_e = sshn_b                                                ! ACC GUESS WORK
-            ubb_e   = ub_b
-            vbb_e   = vb_b
-            sshb_e  = sshn_b
-            ub_e    = ub_b
-            vb_e    = vb_b
          ENDIF
       !
@@ -1029,39 +980,48 @@
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
       !
-      INTEGER  :: ji ,jj, jk
+      INTEGER  :: ji ,jj
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
       REAL(wp) :: zxr2, zyr2, zcmax
-      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zht
+      REAL(wp), POINTER, DIMENSION(:,:) :: zcu
       !!
-!      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
-!      &                  nn_baro, rn_bt_cmax, nn_bt_flt
+      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
+      &                  nn_baro, rn_bt_cmax, nn_bt_flt
       !!----------------------------------------------------------------------
-!      REWIND( numnam )              !* Namelist namsplit: split-explicit free surface
-!      READ  ( numnam, namsplit )
+      !
+      REWIND( numnam_ref )              ! Namelist namsplit in reference namelist : time splitting parameters
+      READ  ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsplit in configuration namelist : time splitting parameters
+      READ  ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
+      WRITE ( numond, namsplit )
+      !
       !         ! Max courant number for ext. grav. waves
       !
-      CALL wrk_alloc( jpi, jpj, zcu, zht )
-      !
-      ! JC: Simplification needed below: define ht_0 even when volume is fixed
+      CALL wrk_alloc( jpi, jpj, zcu )
+      !
       IF (lk_vvl) THEN 
-         zht(:,:) = ht_0(:,:) * tmask(:,:,1)
+         DO jj = 1, jpj
+            DO ji =1, jpi
+               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
+               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
+               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
+            END DO
+         END DO
       ELSE
-         zht(:,:) = 0.
-         DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-      ENDIF
-
-      DO jj = 1, jpj
-         DO ji =1, jpi
-            zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
-            zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
-            zcu(ji,jj) = sqrt(grav*zht(ji,jj)*(zxr2 + zyr2) )
-         END DO
-      END DO
+         DO jj = 1, jpj
+            DO ji =1, jpi
+               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
+               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
+               zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) )
+            END DO
+         END DO
+      ENDIF
 
       zcmax = MAXVAL(zcu(:,:))
       IF( lk_mpp )   CALL mpp_max( zcmax )
 
-      ! Estimate number of iterations to satisfy a max courant number=0.8 
+      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
       IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
       
@@ -1073,34 +1033,36 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
       IF( ln_bt_nn_auto ) THEN
-         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
-         IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
+         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.true. Automatically set nn_baro '
+         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
       ELSE
-         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
-      ENDIF
-      IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
-      IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
-      IF(lwp) WRITE(numout,*) ' Maximum Courant number is   :', zcmax
+         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.false.: Use nn_baro in namelist '
+      ENDIF
 
       IF(ln_bt_av) THEN
-         IF(lwp) WRITE(numout,*) ' ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
+         IF(lwp) WRITE(numout,*) '     ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
       ELSE
-         IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
+         IF(lwp) WRITE(numout,*) '     ln_bt_av=.false. => No time averaging of barotropic variables '
       ENDIF
       !
       !
       IF(ln_bt_fw) THEN
-         IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true.  => Forward integration of barotropic variables '
+         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
       ELSE
-         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
-      ENDIF
-      !
-      IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
+         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
+      ENDIF
+      !
+      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
       SELECT CASE ( nn_bt_flt )
-           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '      Dirac'
-           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = nn_baro'
-           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = 2*nn_baro' 
+           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac'
+           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
+           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
       END SELECT
       !
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
+      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
+      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
+      !
       IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
          CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
@@ -1110,5 +1072,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, zcu, zht )
+      CALL wrk_dealloc( jpi, jpj, zcu )
       !
    END SUBROUTINE dyn_spg_ts_init
@@ -1116,9 +1078,6 @@
 #else
    !!---------------------------------------------------------------------------
-   !!   Default case :   Empty module   No standard free surface constant volume
+   !!   Default case :   Empty module   No split explicit free surface
    !!---------------------------------------------------------------------------
-
-   USE par_kind
-   LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping
 CONTAINS
    INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4370)
@@ -27,5 +27,4 @@
    USE dynadv          ! dynamics: vector invariant versus flux form
    USE dynspg_oce, ONLY: lk_dynspg_ts
-   USE dynspg_ts
 
    IMPLICIT NONE
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 4369)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4370)
@@ -33,7 +33,4 @@
    USE diaar5, ONLY:   lk_diaar5
    USE iom
-   USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div   ! River runoff 
-   USE dynspg_ts,   ONLY: ln_bt_fw
-   USE dynspg_oce, ONLY: lk_dynspg_ts
 #if defined key_agrif
    USE agrif_opa_update
