Index: trunk/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- trunk/DOC/TexFiles/Biblio/Biblio.bib	(revision 4559)
+++ trunk/DOC/TexFiles/Biblio/Biblio.bib	(revision 4560)
@@ -2505,4 +2505,12 @@
   volume = {87},
   pages = {394--409}
+}
+
+@BOOK{Shchepetkin_McWilliams_Bk08,
+  author = {A. F. Shchepetkin and J. C. McWilliams},
+  title = {Handbook of Numerical Analysis, Vol. XIV: Computational Methods for the Ocean and the Atmosphere, pp 121-183},
+  publisher = {P. G. Ciarlet, editor, R. Temam and J. Tribbia, guest eds., Elsevier Science},
+  year = {2008},
+  pages = {784}
 }
 
Index: trunk/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 4559)
+++ trunk/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 4560)
@@ -762,8 +762,11 @@
 \subsection{Split-Explicit free surface (\key{dynspg\_ts})}
 \label{DYN_spg_ts}
+%------------------------------------------namsplit-----------------------------------------------------------
+\namdisplay{namsplit} 
+%-------------------------------------------------------------------------------------------------------------
 
 The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined),
 also called the time-splitting formulation, follows the one 
-proposed by \citet{Griffies_Bk04}. The general idea is to solve the free surface 
+proposed by \citet{Shchepetkin_McWilliams_OM05}. The general idea is to solve the free surface 
 equation and the associated barotropic velocity equations with a smaller time 
 step than $\rdt$, the time step used for the three dimensional prognostic 
@@ -771,35 +774,70 @@
 The size of the small time step, $\rdt_e$ (the external mode or barotropic time step)
  is provided through the \np{nn\_baro} namelist parameter as: 
-$\rdt_e = \rdt / nn\_baro$.
- 
+$\rdt_e = \rdt / nn\_baro$. This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}=true) 
+considering that the stability of the barotropic system is essentially controled by external waves propagation. 
+Maximum allowed Courant number is in that case time independent, and easily computed online from the input bathymetry.
+
+%%%
+The barotropic mode solves the following equations:
+\begin{subequations} \label{Eq_BT}
+  \begin{equation}     \label{Eq_BT_dyn}
+\frac{\partial {\rm \overline{{\bf U}}_h} }{\partial t}=
+ -f\;{\rm {\bf k}}\times {\rm \overline{{\bf U}}_h} 
+-g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \rm {\overline{{\bf U}}_h} + \rm {\overline{\bf G}}
+  \end{equation}
+
+  \begin{equation} \label{Eq_BT_ssh}
+\frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]+P-E
+  \end{equation}
+\end{subequations}
+where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. The third term on the right hand side of \eqref{Eq_BT_dyn} represents the bottom stress (see section \S\ref{ZDF_bfr}), explicitly accounted for at each barotropic iteration. Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm detailed in \citet{Shchepetkin_McWilliams_OM05}. AB3-AM4 coefficients used in \NEMO follow the second-order accurate, "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08} (see their figure 12, lower left). 
 
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
 \begin{figure}[!t]    \begin{center}
-\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_DYN_dynspg_ts.pdf}
+\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_time_split.pdf}
 \caption{  \label{Fig_DYN_dynspg_ts}
 Schematic of the split-explicit time stepping scheme for the external 
-and internal modes. Time increases to the right. 
+and internal modes. Time increases to the right. In this particular exemple, 
+a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_filt=1$) and $nn\_baro=5$.
 Internal mode time steps (which are also the model time steps) are denoted 
-by $t-\rdt$, $t, t+\rdt$, and $t+2\rdt$. 
-The curved line represents a leap-frog time step, and the smaller time 
-steps $N \rdt_e=\frac{3}{2}\rdt$ are denoted by the zig-zag line. 
-The vertically integrated forcing \textbf{M}(t) computed at the model time step $t$ 
-represents the interaction between the external and internal motions. 
-While keeping \textbf{M} and freshwater forcing field fixed, a leap-frog 
-integration carries the external mode variables (surface height and vertically 
-integrated velocity) from $t$ to $t+\frac{3}{2} \rdt$ using N external time 
-steps of length $\rdt_e$. Time averaging the external fields over the 
-$\frac{2}{3}N+1$ time steps (endpoints included) centers the vertically integrated 
-velocity and the sea surface height at the model timestep $t+\rdt$. 
-These averaged values are used to update \textbf{M}(t) with both the surface 
-pressure gradient and the Coriolis force, therefore providing the $t+\rdt$
-velocity.  The model time stepping scheme can then be achieved by a baroclinic 
-leap-frog time step that carries the surface height from $t-\rdt$ to $t+\rdt$.  }
+by $t-\rdt$, $t$ and $t+\rdt$. Variables with $k$ superscript refer to instantaneous barotropic variables, 
+$< >$ and $<< >>$ operator refer to time filtered variables using respectively primary (red vertical bars) and secondary weights (blue vertical bars). 
+The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged 
+transports to advect tracers.
+a) Forward time integration: \np{ln\_bt\_fw}=true, \np{ln\_bt\_ave}=true. 
+b) Centred time integration: \np{ln\_bt\_fw}=false, \np{ln\_bt\_ave}=true. 
+c) Forward time integration with no time filtering (POM-like scheme): \np{ln\_bt\_fw}=true, \np{ln\_bt\_ave}=false. }
 \end{center}    \end{figure}
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
 
-The split-explicit formulation has a damping effect on external gravity waves, 
-which is weaker damping than that for the filtered free surface but still significant, as 
-shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave. 
+In the default case (\np{ln\_bt\_fw}=true), the external mode is integrated 
+between \textit{now} and  \textit{after} baroclinic time-steps (Fig.~\ref{Fig_DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic 
+quantities (\np{ln\_bt\_ave}=true). In that case, the integration is extended slightly beyond  \textit{after} time step to provide time filtered quantities. 
+These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. 
+Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, 
+asselin filtering is not applied to barotropic quantities. \\
+Alternatively, one can choose to integrate barotropic equations starting 
+from \textit{before} time step (\np{ln\_bt\_fw}=false). Although more computationaly expensive ( \np{nn\_baro} additional iterations are indeed necessary), the baroclinic to barotropic forcing term given at \textit{now} time step 
+become centred in the middle of the integration window. It can easily be shown that this property 
+removes part of splitting errors between modes, which increases the overall numerical robustness.
+%references to Patrick Marsaleix' work here. Also work done by SHOM group.
+
+%%%
+
+As far as tracer conservation is concerned, barotropic velocities used to advect tracers must also be updated 
+at \textit{now} time step. This implies to change the traditional order of computations in \NEMO: most of momentum  
+trends (including the barotropic mode calculation) updated first, tracers' after. This \textit{de facto} makes semi-implicit hydrostatic 
+pressure gradient (see section \S\ref{DYN_hpg_imp}) and time splitting not compatible. 
+Advective barotropic velocities are obtained by using a secondary set of filtering weights, uniquely defined from the filter 
+coefficients used for the time averaging (\citet{Shchepetkin_McWilliams_OM05}). Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to obtain exact conservation.
+
+%%%
+
+One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_ave}=false). 
+In that case, external mode equations are continuous in time, ie they are not re-initialized when starting a new 
+sub-stepping sequence. This is the method used so far in the POM model, the stability being maintained by refreshing at (almost) 
+each barotropic time step advection and horizontal diffusion terms. Since the latter terms have not been added in \NEMO for 
+computational efficiency, removing time filtering is not recommended except for debugging purposes. 
+This may be used for instance to appreciate the damping effect of the standard formulation on external gravity waves in idealized or weakly non-linear cases. Although the damping is lower than for the filtered free surface, it is still significant as shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave.
 
 %>>>>>===============
Index: trunk/DOC/TexFiles/Namelist/namdom
===================================================================
--- trunk/DOC/TexFiles/Namelist/namdom	(revision 4559)
+++ trunk/DOC/TexFiles/Namelist/namdom	(revision 4560)
@@ -10,5 +10,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)
Index: trunk/DOC/TexFiles/Namelist/namsplit
===================================================================
--- trunk/DOC/TexFiles/Namelist/namsplit	(revision 4560)
+++ trunk/DOC/TexFiles/Namelist/namsplit	(revision 4560)
@@ -0,0 +1,15 @@
+!-----------------------------------------------------------------------
+&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     "        "
+/
