Index: branches/nemo_v3_3_beta/DOC/NEMO_book.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/NEMO_book.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/NEMO_book.tex	(revision 2376)
@@ -298,4 +298,6 @@
 \include{./TexFiles/Chapters/Chap_MISC}			% Miscellaneous topics
 
+\include{./TexFiles/Chapters/Chap_CFG}			% Predefined configurations
+
 % ================================================================
 % APPENDIX
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Biblio/Biblio.bib	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Biblio/Biblio.bib	(revision 2376)
@@ -28,4 +28,6 @@
 @STRING{E = {Eyrolles}}
 
+@STRING{EFM ={Environ. Fluid Mech.}}
+
 @STRING{GRL = {Geophys. Res. Let.}}
 
@@ -40,4 +42,6 @@
 @STRING{JCP = {J. Comput. Phys.}}
 
+@STRING{JFM = {J. Fluid Mech.}}
+
 @STRING{JGR = {J. Geophys. Res}}
 
@@ -66,9 +70,9 @@
 @STRING{Ocean = {Oceanology}}
 
-@STRING{OD = {Ocean Dynamics}}
+@STRING{OD = {Ocean Dyn.}}
 
 @STRING{OM = {Ocean Modelling}}
 
-@STRING{OS = {Ocean Science}}
+@STRING{OS = {Ocean Sc.}}
 
 @STRING{OUP = {Oxford University Press}}
@@ -93,4 +97,5 @@
 
 @STRING{Tellus = {Tellus}}
+
 
 @ARTICLE{Adcroft_Campin_OM04,
@@ -188,4 +193,15 @@
   year = {2002},
   volume = {107},
+  doi = {10.1029/2001JC000922},
+  url = {http://dx.doi.org/10.1029/2001JC000922}
+}
+
+@ARTICLE{Axell_Liungman_EFM01,
+  author = {L. B. Axell and O. Liungman},
+  title = {A one-equation turbulence model for geophysical applications: 
+  Comparison with data and the k-$\epsilon$ model},
+  journal = EFM,
+  year = {2001},
+  volume = {1}, pages = {71--106},
   doi = {10.1029/2001JC000922},
   url = {http://dx.doi.org/10.1029/2001JC000922}
@@ -609,6 +625,6 @@
 }
 
-@ARTICLE{Craig_Banner_1994,
-  author = {P. D. Banner and M. L. Banner},
+@ARTICLE{Craig_Banner_JPO94,
+  author = {P. D. Craig and M. L. Banner},
   title = {Modeling wave-enhanced turbulence in the ocean surface layer},
   journal = JPO,
@@ -616,4 +632,13 @@
   volume = {24},  number = {12},
   pages = {2546--2559}
+}
+
+@ARTICLE{Craik_Leibovich_JFM76,
+  author = {A. D. D. Craik and S. Leibovich},
+  title = {A rational model for Langmuir circulations},
+  journal = JFM,
+  year = {1976},
+  volume = {73}, 
+  pages = {401--426}
 }
 
@@ -634,4 +659,13 @@
   year = {2001},
   pages = {163pp}
+}
+
+@ARTICLE{D'Alessio_al_JPO98,
+  author = {S. J. D. D'Alessio and K. Abdella and N. A. McFarlane},
+  title = {A new second-order turbulence closure scheme for modeling the oceanic mixed layer},
+  journal = {JPO},
+  year = {1998},
+  volume = {28},
+  pages = {1624--1641}
 }
 
@@ -1685,4 +1719,13 @@
 }
 
+@ARTICLE{Li_Garrett_JMR93,
+  author = {M. Li and C. Garrett},
+  title = {Cell merging and the jet/downwelling ratio in Langmuir circulation},
+  journal = JMR,
+  year = {1993},
+  volume = {51},
+  pages = {737--769}
+}
+
 @TECHREPORT{Lott1989,
   author = {F. Lott and G. Madec},
@@ -1958,10 +2001,12 @@
 
 @ARTICLE{Mellor_Blumberg_JPO04,
-  author = {G. Mellor and A. Blumberg},
-  title = {Wave Breaking and Ocean Surface Layer Thermal Response},
-  journal = JPO,
-  year = {2004},
-  volume = {34},
-  pages = {693--698}
+author = {G. Mellor and A. Blumberg},
+title = {Wave Breaking and Ocean Surface Layer Thermal Response},
+journal = {JPO},
+volume = {34}, number = {3},
+pages = {693--698},
+year = {2004},
+doi = {10.1175/2517.1},
+URL = {http://journals.ametsoc.org/doi/abs/10.1175/2517.1}
 }
 
@@ -2379,4 +2424,14 @@
   doi = {10.1029/2002GL015633},
   url = {http://dx.doi.org/10.1029/2002GL015633}
+}
+
+@ARTICLE{Stacey_JPO99,
+  author = {M. W. Stacey},
+  title = {Simulations of the wind-forced near-surface circulation in Knight Inlet: 
+             A parameterization of the roughness length},
+  journal = JPO,
+  year = {1999},
+  volume = {29},
+  pages = {1363--1367},
 }
 
@@ -2558,4 +2613,15 @@
 }
 
+@ARTICLE{Warner_al_OM05,
+  author = {J. C. Warner and C. R. Sherwood and H. G. Arango and R. P. Signell},
+  title = {Performance of four turbulence closure models implemented using a generic length scale method},
+  journal = OM,
+  year = {2005},
+  volume = {8},
+  pages = {81--113},
+  doi = {10.1016/j.ocemod.2003.12.003},
+  url = {http://dx.doi.org/j.ocemod.2003.12.003}
+}
+
 @ARTICLE{Weatherly_JMR84,
   author = {G. L. Weatherly},
@@ -2564,5 +2630,5 @@
   year = {1984},
   volume = {42},  number = {2},
-  pages = {289-301}
+  pages = {289--301}
 }
 
@@ -2582,5 +2648,5 @@
   year = {1998},
   volume = {15},  number = {5},
-  pages = {1171-1187},
+  pages = {1171--1187},
 }
 
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_E.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_E.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_E.tex	(revision 2376)
@@ -300,5 +300,6 @@
 \begin{center}
 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_ISO_triad.pdf}
-\caption{Triads used in the Griffies's like iso-neutral diffision scheme for 
+\caption{  \label{Fig_ISO_triad}   
+Triads used in the Griffies's like iso-neutral diffision scheme for 
 $u$-component (upper panel) and $w$-component (lower panel).}
 \end{center}
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_ISO.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_ISO.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Annex_ISO.tex	(revision 2376)
@@ -100,8 +100,9 @@
 \begin{figure}[h] \begin{center}
     \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_triad_fluxes}
-    \caption{(a) Arrangement of triads $S_i$ and tracer gradients to
-      give lateral tracer flux from box $i,k$ to $i+1,k$ (b) Triads
-      $S'_i$ and tracer gradients to give vertical tracer flux from
-      box $i,k$ to $i,k+1$.}
+    \caption{  \label{Fig_ISO_triad}  
+      (a) Arrangement of triads $S_i$ and tracer gradients to
+           give lateral tracer flux from box $i,k$ to $i+1,k$ 
+      (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from
+            box $i,k$ to $i,k+1$.}
     \label{Fig_triad}
   \end{center} \end{figure}
@@ -168,5 +169,6 @@
 \begin{figure}[h] \begin{center}
     \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_qcells}
-    \caption{Triad notation for quarter cells.T-cells are inside
+    \caption{   \label{Fig_ISO_triad_notation}  
+    Triad notation for quarter cells.T-cells are inside
       boxes, while the  $i+\half,k$ u-cell is shaded in green and the
       $i,k+\half$ w-cell is shaded in pink.}
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_CFG.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_CFG.tex	(revision 2376)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_CFG.tex	(revision 2376)
@@ -0,0 +1,274 @@
+% ================================================================
+% Chapter Ñ Configurations
+% ================================================================
+\chapter{Configurations}
+\label{MISC}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new ligne
+
+% ================================================================
+% Introduction
+% ================================================================
+\section{Introduction}
+\label{CFG_intro}
+
+
+The purpose of this part of the manual is to introduce the \NEMO predefined configuration. 
+These configurations are offered as means to explore various numerical and physical options, 
+thus allowing the user to verify that the code is performing in a manner consistent with that 
+we are running. This form of verification is critical as one adopts the code for his or her particular 
+research purposes. The test cases also provide a sense for some of the options available 
+in the code, though by no means are all options exercised in the predefined configurations.
+
+
+%There is several predefined ocean configuration which use is controlled by a specific CPP key. 
+
+%The key set the domain sizes (\jp{jpiglo}, \jp{jpjglo}, \jp{jpk}), the mesh and the bathymetry, 
+%and, in some cases, add to the model physics some specific treatments.
+
+
+% ================================================================
+% 1D model functionality
+% ================================================================
+\section{Water column model: 1D model (\key{c1d})}
+\label{CFG_c1d}
+
+The 1D model option simulates a stand alone water column within the 3D \NEMO system. 
+It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers 
+or a biogeochemical model. It is set up by defining the \key{c1d} CPP key. 
+The 1D model is a very useful tool  
+\textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes ; 
+\textit{(b)} to investigate suitable parameterisations of unresolved turbulence (wind steering, 
+langmuir circulation, skin layers, ...) ; 
+\textit{(c)} to compare the behaviour of different vertical mixing schemes  ; 
+\textit{(d)} to perform sensitivity studies on the vertical diffusion at a particular point of an ocean domain ; 
+\textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model.
+
+The methodology is based on the use of the zoom functionality over the smallest possible 
+domain : a 3 x 3 domain centred on the grid point of interest (see \S\ref{MISC_zoom}), 
+with some extra routines. There is no need to define a new mesh, bathymetry, 
+initial state or forcing, since the 1D model will use those of the configuration it is a zoom of. 
+The chosen grid point is set in par\_oce.F90 module by setting the jpizoom and jpjzoom 
+parameters to the indices of the location of the chosen grid point.
+
+
+% ================================================================
+% ORCA family configurations
+% ================================================================
+\section{ORCA family: global ocean with tripolar grid}
+\label{CFG_orca}
+
+The ORCA family is a series of global ocean configurations that are run together with 
+the LIM sea-ice model (ORCA-LIM) and possibly with PISCES biogeochemical model 
+(ORCA-LIM-PISCES), using various resolutions.
+
+% -------------------------------------------------------------------------------------------------------------
+%       ORCA tripolar grid
+% -------------------------------------------------------------------------------------------------------------
+\subsection{ORCA tripolar grid}
+\label{CFG_orca_grid}
+
+The ORCA grid is a tripolar is based on the semi-analytical method of \citet{Madec_Imbard_CD96}. 
+It allows to construct a global orthogonal curvilinear ocean mesh which has no singularity point inside 
+the computational domain since two north mesh poles are introduced and placed on lands.
+The method involves defining an analytical set of mesh parallels in the stereographic polar plan, 
+computing the associated set of mesh meridians, and projecting the resulting mesh onto the sphere. 
+The set of mesh parallels used is a series of embedded ellipses which foci are the two mesh north 
+poles (Fig.~\ref{Fig_MISC_ORCA_msh}). The resulting mesh presents no loss of continuity in 
+either the mesh lines or the scale factors, or even the scale factor derivatives over the whole 
+ocean domain, as the mesh is not a composite mesh. 
+
+The method is applied to Mercator grid ($i.e.$ same zonal and meridional grid spacing) poleward 
+of $20\deg$N, so that the Equator is a mesh line, which provides a better numerical solution 
+for equatorial dynamics. The choice of the series of embedded ellipses (position of the foci and 
+variation of the ellipses) is a compromise between maintaining  the ratio of mesh anisotropy 
+($e_1 / e_2$) close to one in the ocean (especially in area of strong eddy activities such as 
+the Gulf Stream) and keeping the smallest scale factor in the northern hemisphere larger 
+than the smallest one in the southern hemisphere.
+The resulting mesh is shown in Fig.~\ref{Fig_MISC_ORCA_msh} and \ref{Fig_MISC_ORCA_e1e2} 
+for a half a degree grid (ORCA\_R05). The smallest ocean scale factor is found in along  
+Antarctica, while the ratio of anisotropy remains close to one except near the Victoria Island 
+in the Canadian Archipelago. 
+
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]   \begin{center}
+\includegraphics[width=0.98\textwidth]{./TexFiles/Figures/Fig_ORCA_NH_mesh.pdf}
+\caption{  \label{Fig_MISC_ORCA_msh}     
+ORCA mesh conception. The departure from an isotropic Mercator grid start poleward of 20\deg N.
+The two "north pole" are the foci of a series of embedded ellipses (blue curves) 
+which are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes). 
+Then, following \citet{Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed 
+which provide the j-lines of the mesh (pseudo longitudes).  }
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!tbp]  \begin{center}
+\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_ORCA_NH_msh05_e1_e2.pdf}
+\includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_ORCA_aniso.pdf}
+\caption {  \label{Fig_MISC_ORCA_e1e2}
+\textit{Top}: Horizontal scale factors ($e_1$, $e_2$) and 
+\textit{Bottom}: ratio of anisotropy ($e_1 / e_2$)
+for ORCA 0.5\deg ~mesh. South of 20\deg N a Mercator grid is used ($e_1 = e_2$) 
+so that the anisotropy ratio is 1. Poleward of 20\deg N, the two "north pole" 
+introduce a weak anisotropy over the ocean areas ($< 1.2$) except in vicinity of Victoria Island 
+(Canadian Arctic Archipelago). }
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+
+
+% -------------------------------------------------------------------------------------------------------------
+%       ORCA-LIM(-PISCES) configurations
+% -------------------------------------------------------------------------------------------------------------
+\subsection{ORCA-LIM(-PISCES) configurations}
+\label{CFG_orca_grid}
+
+
+The NEMO system is provided with four built-in ORCA configurations which differ in the 
+horizontal resolution
+\footnote{the value of the resolution is given by the resolution at the Equator expressed in degrees.} 
+used:
+\begin{description}
+\item[\key{orca\_r4}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~4
+\item[\key{orca\_r2}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~2
+%\item[\key{orca\_r1}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~1
+\item[\key{orca\_r05}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~05
+\item[\key{orca\_r025}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~025
+\end{description}
+
+The ORCA\_R2 configuration has the following specificity : starting from a 2\deg~ORCA mesh, 
+local mesh refinements were applied to the Mediterranean, Red, Black and Caspian Seas, 
+so that the resolution is $1\deg \time 1\deg$ there. A local transformation were also applied 
+with in the Tropics in order to refine the meridional resolution up to 0.5\deg at the Equator.
+
+The ORCA\_R1 configuration has only a local tropical transformation  to refine the meridional 
+resolution up to 1/3\deg~at the Equator. Note that the tropical mesh refinements in ORCA\_R2 
+and R1 strongly increases the mesh anisotropy there.
+
+The ORCA\_R05 configuration and higher ones does not incorporate any regional refinements.  
+
+Only the ORCA\_R2 is provided with all its input files in the \NEMO distribution. 
+It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC 
+assessment of climate change (Marti et al., 2009). It is also the basis for the \NEMO contribution 
+to the Coordinate Ocean-ice Reference Experiments (COREs) documented in \citet{Griffies_al_OM09}. 
+
+This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m) 
+in the upper 150m. The bottom topography and the coastlines are derived 
+from the global atlas of Smith and Sandwell (1997). The default forcing employ the boundary 
+forcing from \citet{Large_Yeager_Rep04} (see \S\ref{SBC_blk_core}), 
+which was developed for the purpose of running global coupled ocean-ice simulations 
+without an interactive atmosphere. This \citet{Large_Yeager_Rep04} dataset is available 
+through the GFDL web site \footnote{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}. 
+The "normal year" of \citet{Large_Yeager_Rep04} has been chosen of the \NEMO distribution 
+since release v3.3. 
+
+The vertical resolution can be increased by a factor of 10 by defining the \key{orca\_lev10} CPP key. 
+It can also be run with an AGRIF zoom over the Agulhas current area ( \key{agrif}  defined). 
+Also available are to keys, \key{arctic} and \key{antarctic}, which allows to run a regional Arctic 
+or peri-Antarctic configuration extracted from an ORCA configuration. (This does not work with ORCA\_R4 and R1).
+
+
+%--------------------------------------------------TABLE--------------------------------------------------
+\begin{table}[htbp]     \begin{center}
+\begin{tabular}{ccccc}
+key          	             & \jp{jp\_cfg} &  \jp{jpiglo} & \jp{jpiglo} &       \\  
+\hline  \hline
+\key{orca\_r4} 	    &        4         &         92     &      76      &       \\
+\key{orca\_r2}       &        2         &       182     &    149      &        \\
+%\key{orca\_r1}       &        1         &       362     &     511     &        \\
+\key{orca\_r05}     &        05       &       722     &     261     &        \\
+\key{orca\_r025}   &        025     &      1442    &   1021     &        \\
+%\key{orca\_r8}       &        8         &      2882    &   2042     &        \\
+%\key{orca\_r12}     &      12         &      4322    &   3062      &       \\
+\hline
+\hline
+\end{tabular}
+\caption{ \label{Tab_ORCA}   
+Set of predefined ORCA parameters. }
+\end{center}
+\end{table}
+%--------------------------------------------------------------------------------------------------------------
+
+
+% -------------------------------------------------------------------------------------------------------------
+%       GYRE family: double gyre basin
+% -------------------------------------------------------------------------------------------------------------
+\section{GYRE family: double gyre basin (\key{gyre})}
+\label{MISC_config_gyre}
+
+The GYRE configuration \citep{Levy_al_OM10} have been built to simulated 
+the seasonal cycle of a double-gyre box model. It consist in an idealized domain 
+similar to that used in the studies of \citet{Drijfhout_JPO94} and \citet{Hazeleger_Drijfhout_JPO98, 
+Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00}, 
+over which an analytical seasonal forcing is applied. This allows to investigate the 
+spontaneous generation of a large number of interacting, transient mesoscale eddies 
+and their contribution to the large scale circulation. 
+
+The domain geometry is a closed rectangular basin on the $\beta$-plane centred 
+at $\sim 30\deg$N and rotated by 45\deg, 3180~km long, 2120~km wide 
+and 4~km deep (Fig.~\ref{Fig_MISC_strait_hand}). 
+The domain is bounded by vertical walls and by a ßat bottom. The configuration is 
+meant to represent an idealized North Atlantic or North Pacific basin. 
+The circulation is forced by analytical profiles of wind and buoyancy ßuxes. 
+The applied forcings vary seasonally in a sinusoidal manner between winter 
+and summer extrema \citep{Levy_al_OM10}. 
+The wind stress is zonal and its curl changes sign at 22\deg N and 36\deg N. 
+It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain 
+and a small recirculation gyre in the southern corner. 
+The net heat ßux takes the form of a restoring toward a zonal apparent air 
+temperature profile. A portion of the net heat ßux which comes from the solar radiation
+is allowed to penetrate within the water column. 
+The fresh water ßux is also prescribed and varies zonally. 
+It is determined such as, at each time step, the basin-integrated ßux is zero. 
+The basin is initialised at rest with vertical profiles of temperature and salinity 
+uniformly applied to the whole domain.
+
+The GYRE configuration is set through the \key{gyre} CPP key. Its horizontal resolution 
+(and thus the size of the domain) is determined by setting \jp{jp\_cfg} in \hf{par\_GYRE} file: \\
+\jp{jpiglo} $= 30 \times$ \jp{jp\_cfg} + 2   \\
+\jp{jpjglo} $= 20 \times$ \jp{jp\_cfg} + 2   \\
+Obviously, the namelist parameters have to be adjusted to the chosen resolution.
+In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}=31) (Fig.~\ref{Fig_zgr}).
+
+The GYRE configuration is also used in benchmark test as it is very simple to increase 
+its resolution and as it does not requires any input file. For example, keeping a same model size 
+on each processor while increasing the number of processor used is very easy, even though the 
+physical integrity of the solution can be compromised.
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]   \begin{center}
+\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_GYRE.pdf}
+\caption{  \label{Fig_GYRE}   
+Snapshot of relative vorticity at the surface of the model domain 
+in GYRE R9, R27 and R54. From \citet{Levy_al_OM10}.}
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+% -------------------------------------------------------------------------------------------------------------
+%       EEL family configuration
+% -------------------------------------------------------------------------------------------------------------
+\section{EEL family: periodic channel}
+\label{MISC_config_EEL}
+
+\begin{description}
+\item[\key{eel\_r2}]  to be described....
+\item[\key{eel\_r5}]  
+\item[\key{eel\_r6}]  
+\end{description}
+
+% -------------------------------------------------------------------------------------------------------------
+%       POMME configuration
+% -------------------------------------------------------------------------------------------------------------
+\section{POMME: mid-latitude sub-domain}
+\label{MISC_config_POMME}
+
+
+\key{pomme\_r025} : to be described....
+
+
+
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 2376)
@@ -1,3 +1,2 @@
-
 % ================================================================
 % Chapter 2 Ñ Space and Time Domain (DOM)
@@ -40,7 +39,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tb] \label{Fig_cell}  \begin{center}
+\begin{figure}[!tb]    \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf}
-\caption{Arrangement of variables. $t$ indicates scalar points where temperature, 
+\caption{ \label{Fig_cell}    
+Arrangement of variables. $t$ indicates scalar points where temperature, 
 salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) 
 indicates vector points, and $f$ indicates vorticity points where both relative and 
@@ -80,4 +80,5 @@
 as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section). 
 
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{table}[!tb] \label{Tab_cell}
 \begin{center} \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|}
@@ -92,10 +93,12 @@
 fw	& $i+1/2$ 	& $j+1/2$	& $k+1/2$ 	\\ \hline
 \end{tabular}
-\caption{Location of grid-points as a function of integer or integer and a half value 
-of the column, line or level. This indexing is only used for the writing of the semi-
-discrete equation. In the code, the indexing uses integer values only and has a 
-reverse direction in the vertical (see \S\ref{DOM_Num_Index})}
+\caption{ \label{Tab_cell}   
+Location of grid-points as a function of integer or integer and a half value of the column, 
+line or level. This indexing is only used for the writing of the semi-discrete equation. 
+In the code, the indexing uses integer values only and has a reverse direction 
+in the vertical (see \S\ref{DOM_Num_Index})}
 \end{center}
 \end{table}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 % -------------------------------------------------------------------------------------------------------------
@@ -206,9 +209,9 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tb] \label{Fig_index_hor}  \begin{center}
+\begin{figure}[!tb]  \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_index_hor.pdf}
-\caption{Horizontal integer indexing used in the \textsc{Fortran} code. The dashed 
-area indicates the cell in which variables contained in arrays have the same 
-$i$- and $j$-indices}
+\caption{   \label{Fig_index_hor}    
+Horizontal integer indexing used in the \textsc{Fortran} code. The dashed area indicates 
+the cell in which variables contained in arrays have the same $i$- and $j$-indices}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -256,7 +259,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!pt] \label{Fig_index_vert}  \begin{center}
+\begin{figure}[!pt]    \begin{center}
 \includegraphics[width=.90\textwidth]{./TexFiles/Figures/Fig_index_vert.pdf}
-\caption{Vertical integer indexing used in the \textsc{Fortran } code. Note that 
+\caption{ \label{Fig_index_vert}     
+Vertical integer indexing used in the \textsc{Fortran } code. Note that 
 the $k$-axis is orientated downward. The dashed area indicates the cell in 
 which variables contained in arrays have the same $k$-index.}
@@ -364,9 +368,10 @@
 Fig.~\ref{Fig_zgr_e3}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_zgr_e3}  \begin{center}
+\begin{figure}[!t]     \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr_e3.pdf}
-\caption{Comparison of (a) traditional definitions of grid-point position and grid-size 
-in the vertical, and (b) analytically derived grid-point position and scale factors. For 
-both grids here,  the same $w$-point depth has been chosen but in (a) the 
+\caption{ \label{Fig_zgr_e3}    
+Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
+and (b) analytically derived grid-point position and scale factors. 
+For both grids here,  the same $w$-point depth has been chosen but in (a) the 
 $t$-points are set half way between $w$-points while in (b) they are defined from 
 an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$. 
@@ -471,7 +476,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tb] \label{Fig_z_zps_s_sps}  \begin{center}
+\begin{figure}[!tb]    \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zps_s_sps.pdf}
-\caption{The ocean bottom as seen by the model: 
+\caption{  \label{Fig_z_zps_s_sps}   
+The ocean bottom as seen by the model: 
 (a) $z$-coordinate with full step, 
 (b) $z$-coordinate with partial step, 
@@ -575,8 +581,9 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tb] \label{Fig_zgr}  \begin{center}
+\begin{figure}[!tb]    \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr.pdf}
-\caption{Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level 
-functions for (a) T-point depth and (b) the associated scale factor as computed 
+\caption{ \label{Fig_zgr}    
+Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level functions for 
+(a) T-point depth and (b) the associated scale factor as computed 
 from \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.}
 \end{center}   \end{figure}
@@ -651,6 +658,5 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{table} \label{Tab_orca_zgr}
-\begin{center} \begin{tabular}{c||r|r|r|r}
+\begin{table}     \begin{center} \begin{tabular}{c||r|r|r|r}
 \hline
 \textbf{LEVEL}& \textbf{gdept}& \textbf{gdepw}& \textbf{e3t }& \textbf{e3w  } \\ \hline
@@ -687,7 +693,7 @@
 31	& 	\textbf{5250.23}&	5000.00 &	\textbf{500.56} &	500.33 \\ \hline
 \end{tabular} \end{center} 
-\caption{Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration 
-as computed from \eqref{DOM_zgr_ana} using the coefficients given in 
-\eqref{DOM_zgr_coef}}
+\caption{ \label{Tab_orca_zgr}   
+Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed 
+from \eqref{DOM_zgr_ana} using the coefficients given in \eqref{DOM_zgr_coef}}
 \end{table}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -766,7 +772,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tb] \label{Fig_sco_function}  \begin{center}
+\begin{figure}[!tb]    \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf}
-\caption{Examples of the stretching function applied to a sea mont; from left to right: 
+\caption{  \label{Fig_sco_function}   
+Examples of the stretching function applied to a sea mont; from left to right: 
 surface, surface and bottom, and bottom intensified resolutions}
 \end{center}   \end{figure}
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2376)
@@ -291,11 +291,10 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \label{Fig_DYN_een_triad}
-\begin{center}
+\begin{figure}[!ht]    \begin{center}
 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_DYN_een_triad.pdf}
-\caption{Triads used in the energy and enstrophy conserving scheme (een) for 
+\caption{ \label{Fig_DYN_een_triad}  
+Triads used in the energy and enstrophy conserving scheme (een) for 
 $u$-component (upper panel) and $v$-component (lower panel).}
-\end{center}
-\end{figure}
+\end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -769,8 +768,8 @@
 
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
-\begin{figure}[!t] \label{Fig_DYN_dynspg_ts}
-\begin{center}
+\begin{figure}[!t]    \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_DYN_dynspg_ts.pdf}
-\caption{Schematic of the split-explicit time stepping scheme for the external 
+\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. 
 Internal mode time steps (which are also the model time steps) are denoted 
@@ -790,6 +789,5 @@
 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$.  }
-\end{center}
-\end{figure}
+\end{center}    \end{figure}
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
 
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_LBC.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_LBC.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_LBC.tex	(revision 2376)
@@ -51,8 +51,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_LBC_uv}  \begin{center}
+\begin{figure}[!t]     \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_uv.pdf}
-\caption {Lateral boundary (thick line) at T-level. The velocity normal to the 
-		 boundary is set to zero.}
+\caption{  \label{Fig_LBC_uv}
+Lateral boundary (thick line) at T-level. The velocity normal to the boundary is set to zero.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -74,10 +74,11 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!p] \label{Fig_LBC_shlat}  \begin{center}
+\begin{figure}[!p] \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_shlat.pdf}
-\caption {lateral boundary condition (a) free-slip ($rn\_shlat=0$) ; (b) no-slip ($rn\_shlat=2$) 
+\caption{     \label{Fig_LBC_shlat} 
+lateral boundary condition (a) free-slip ($rn\_shlat=0$) ; (b) no-slip ($rn\_shlat=2$) 
 ; (c) "partial" free-slip ($0<rn\_shlat<2$) and (d) "strong" no-slip ($2<rn\_shlat$). 
 Implied "ghost" velocity inside land area is display in grey. }
-\end{center}   \end{figure}
+\end{center}    \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -192,7 +193,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_LBC_jperio}  \begin{center}
+\begin{figure}[!t]     \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_jperio.pdf}
-\caption {setting of (a) east-west cyclic  (b) symmetric across the equator boundary conditions.}
+\caption{    \label{Fig_LBC_jperio}
+setting of (a) east-west cyclic  (b) symmetric across the equator boundary conditions.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -209,9 +211,10 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_North_Fold_T}  \begin{center}
+\begin{figure}[!t]    \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_North_Fold_T.pdf}
-\caption {North fold boundary with a $T$-point pivot and cyclic east-west boundary condition 
-($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds to the inner 
-domain mask (see text). }
+\caption{    \label{Fig_North_Fold_T} 
+North fold boundary with a $T$-point pivot and cyclic east-west boundary condition 
+($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds 
+to the inner domain mask (see text). }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -274,7 +277,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_mpp}  \begin{center}
+\begin{figure}[!t]    \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mpp.pdf}
-\caption {Positioning of a sub-domain when massively parallel processing is used. }
+\caption{   \label{Fig_mpp} 
+Positioning of a sub-domain when massively parallel processing is used. }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -363,11 +367,13 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \label{Fig_mppini2}  \begin{center}
+\begin{figure}[!ht]     \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mppini2.pdf}
-\caption {Example of Atlantic domain defined for the CLIPPER projet. Initial grid is 
-composed of 773 x 1236 horizontal points. (a) the domain is split onto 9 \time 20 
-subdomains (jpni=9, jpnj=20). 52 subdomains are land areas. (b) 52 subdomains 
-are eliminated (white rectangles) and the resulting number of processors really 
-used during the computation is jpnij=128.}
+\caption {    \label{Fig_mppini2}
+Example of Atlantic domain defined for the CLIPPER projet. Initial grid is 
+composed of 773 x 1236 horizontal points. 
+(a) the domain is split onto 9 \time 20 subdomains (jpni=9, jpnj=20). 
+52 subdomains are land areas. 
+(b) 52 subdomains are eliminated (white rectangles) and the resulting number 
+of processors really used during the computation is jpnij=128.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -444,8 +450,5 @@
 
 %--------------------------------------------------TABLE--------------------------------------------------
-
-\begin{table}[htbp]  \label{Tab_obc_param}
-\begin{center}
-\begin{tabular}{|l|c|c|c|}
+\begin{table}[htbp]     \begin{center}    \begin{tabular}{|l|c|c|c|}
 \hline
 Boundary and  & Constant index  & Starting index (d\'{e}but) & Ending index (fin) \\
@@ -464,10 +467,11 @@
 lp\_obc\_north  & $j$-index of a $v$ point      & $i$  of a $T$ point & $i$ of a $T$ point \\
 \hline
-\end{tabular}
-\end{center}
-\caption{Names of different indices relating to the open boundaries. In the case 
+\end{tabular}    \end{center}
+\caption{     \label{Tab_obc_param}
+Names of different indices relating to the open boundaries. In the case 
 of a completely open ocean domain with four ocean boundaries, the parameters 
 take exactly the values indicated.}
 \end{table}
+%------------------------------------------------------------------------------------------------------------
 
 The open boundaries must be along coordinate lines. On the C-grid, the boundary 
@@ -496,9 +500,9 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_obc_north}  \begin{center}
+\begin{figure}[!t]     \begin{center}
 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf}
-\caption {Localization of the North open boundary points.}
-\end{center}  
-\end{figure}
+\caption{    \label{Fig_obc_north}
+Localization of the North open boundary points.}
+\end{center}     \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -559,8 +563,5 @@
 
 %--------------------------------------------------TABLE--------------------------------------------------
-
-\begin{table}[htbp]  \label{Tab_obc_ind}
-\begin{center}
-\begin{tabular}{|l|c|c|c|c|c|}
+\begin{table}[htbp]     \begin{center}      \begin{tabular}{|l|c|c|c|c|c|}
 \hline
 OBC  & Variable   & file name      & Index  & Start  & end  \\
@@ -581,7 +582,7 @@
       &    V      &   obcnorth\_V.nc  &  $je$-2     &  $ib$+1 &  $ie-1$  \\  
 \hline
-\end{tabular}
-\end{center}
-\caption{Requirements for creating open boundary files from a global configuration, 
+\end{tabular}     \end{center}
+\caption{    \label{Tab_obc_ind}
+Requirements for creating open boundary files from a global configuration, 
 appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the 
 $i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global 
@@ -590,4 +591,5 @@
 $-F$ $-d\;y,je-2$  $-d\;x,ib+1,ie-1$ } 
 \end{table}
+%-----------------------------------------------------------------------------------------------------------
 
 It is assumed that the open boundary files contain the variables for the period of 
@@ -878,7 +880,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_LBC_bdy_geom}  \begin{center}
+\begin{figure}[!t]      \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_bdy_geom.pdf}
-\caption {Example of geometry of unstructured open boundary}
+\caption {      \label{Fig_LBC_bdy_geom}
+Example of geometry of unstructured open boundary}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -923,7 +926,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_LBC_nc_header}  \begin{center}
+\begin{figure}[!t]     \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_nc_header.pdf}
-\caption {Example of header of netcdf input data file for BDY}
+\caption {     \label{Fig_LBC_nc_header} 
+Example of header of netcdf input data file for BDY}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_LDF.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 2376)
@@ -353,8 +353,9 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \label{Fig_LDF_ZDF1}  \begin{center}
+\begin{figure}[!ht]      \begin{center}
 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_LDF_ZDF1.pdf}
-\caption {averaging procedure for isopycnal slope computation.}
-\end{center}   \end{figure}
+\caption {    \label{Fig_LDF_ZDF1}
+averaging procedure for isopycnal slope computation.}
+\end{center}    \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -380,7 +381,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \label{Fig_eiv_slp}  \begin{center}
+\begin{figure}[!ht]     \begin{center}
 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_eiv_slp.pdf}
-\caption {Vertical profile of the slope used for lateral mixing in the mixed layer : 
+\caption {     \label{Fig_eiv_slp}
+Vertical profile of the slope used for lateral mixing in the mixed layer : 
 \textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior, 
 which has to be adjusted at the surface boundary (i.e. it must tend to zero at the 
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_MISC.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 2376)
@@ -58,8 +58,9 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tbp] \label{Fig_MISC_strait_hand}  \begin{center}
+\begin{figure}[!tbp] 	 \begin{center}
 \includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_Gibraltar.pdf}
 \includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_Gibraltar2.pdf}
-\caption {Example of the Gibraltar strait defined in a $1\deg \times 1\deg$ mesh. 
+\caption{	\label{Fig_MISC_strait_hand} 
+Example of the Gibraltar strait defined in a $1\deg \times 1\deg$ mesh. 
 \textit{Top}: using partially open cells. The meridional scale factor at $v$-point 
 is reduced on both sides of the strait to account for the real width of the strait 
@@ -134,34 +135,11 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \label{Fig_LBC_zoom}  \begin{center}
+\begin{figure}[!ht] 	  \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_zoom.pdf}
-\caption {Position of a model domain compared to the data input domain when the zoom functionality is used.}
+\caption{	\label{Fig_LBC_zoom}
+Position of a model domain compared to the data input domain when the zoom functionality is used.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-
-% ================================================================
-% 1D model functionality
-% ================================================================
-\section{Water column model: 1D model (\key{c1d})}
-\label{MISC_1D}
-
-The 1D model option simulates a stand alone water column within the 3D \NEMO system. 
-It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers 
-or a biogeochemical model. It is set up by defining the \key{c1d} CPP key. 
-The 1D model is a very useful tool  
-\textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes ; 
-\textit{(b)} to investigate suitable parameterisations of unresolved turbulence (wind steering, 
-langmuir circulation, skin layers) ; 
-\textit{(c)} to compare the behaviour of different vertical mixing schemes  ; 
-\textit{(d)} to perform sensitivity studies on the vertical diffusion at a particular point of an ocean domain ; 
-\textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model.
-
-The methodology is based on the use of the zoom functionality over the smallest possible 
-domain : a 3 x 3 domain centred on the grid point of interest (see \S\ref{MISC_zoom}), 
-with some extra routines. There is no need to define a new mesh, bathymetry, 
-initial state or forcing, since the 1D model will use those of the configuration it is a zoom of. 
-The chosen grid point is set in par\_oce.F90 module by setting the jpizoom and jpjzoom 
-parameters to the indices of the location of the chosen grid point.
 
 % ================================================================
@@ -260,10 +238,12 @@
 The "bit comparison" option has been introduced in order to be able to check that 
 mono-processor and multi-processor runs give exactly the same results. 
+%THIS is to be updated with the mpp_sum_glo  introduced in v3.3
+% nn_bit_cmp  today only check that the nn_cla = 0 (no cross land advection)
 
 $\bullet$  Benchmark (\np{nn\_bench}). This option defines a benchmark run based on 
-a GYRE configuration in which the resolution remains the same whatever the domain 
-size. This allows a very large model domain to be used, just by changing the domain 
-size (\jp{jpiglo}, \jp{jpjglo}) and without adjusting either the time-step or the physical 
-parameterisations. 
+a GYRE configuration (see \S\ref{CFG_gyre}) in which the resolution remains the same 
+whatever the domain size. This allows a very large model domain to be used, just by 
+changing the domain size (\jp{jpiglo}, \jp{jpjglo}) and without adjusting either the time-step 
+or the physical parameterisations. 
 
 
@@ -607,6 +587,6 @@
 volume ratio of each processing region.
 
-\begin{table} 
-\begin{tabular}{lrrr}
+%------------------------------------------TABLE----------------------------------------------------
+\begin{table} 	\begin{tabular}{lrrr}
 Filename & NetCDF3 & NetCDF4 & Reduction\\
          &filesize & filesize & \% \\
@@ -638,7 +618,8 @@
 ORCA2\_2d\_grid\_W\_0007.nc & 4416 & 1368 & 70\%\\
 \end{tabular}
-\caption{\label{Tab_NC4} Filesize comparison between NetCDF3 and NetCDF4 
-with chunking and compression}
+\caption{ 	\label{Tab_NC4} 
+Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression}
 \end{table}
+%----------------------------------------------------------------------------------------------------
 
 Since version 3.2, an I/O server has been added which provides more
@@ -758,7 +739,8 @@
 %-------------------------------------------------------------------------------------------------------------
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_mask_subasins}  \begin{center}
+\begin{figure}[!t] 	  \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf}
-\caption {Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute
+\caption{	\label{Fig_mask_subasins}
+Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute
 the heat and salt transports as well as the meridional stream-function: Atlantic basin (red), 
 Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green). 
@@ -931,203 +913,2 @@
 
 % ================================================================
-% predefined configurations
-% ================================================================
-\section{predefined configurations}
-\label{MISC_config}
-
-There is several predefined ocean configuration which use is controlled by a specific CPP key. 
-
-The key set the domain sizes (\jp{jpiglo}, \jp{jpjglo}, \jp{jpk}), the mesh and the bathymetry, 
-and, in some cases, add to the model physics some specific treatments.
-
-% -------------------------------------------------------------------------------------------------------------
-%       ORCA family configurations
-% -------------------------------------------------------------------------------------------------------------
-\subsection{ORCA family: global ocean with tripolar grid}
-\label{MISC_config_orca}
-
-The NEMO system is provided with four built-in ORCA configurations which differ in the 
-horizontal resolution used:
-\begin{description}
-\item[\key{orca\_r4}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~4
-\item[\key{orca\_r2}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~2
-\item[\key{orca\_r1}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~1
-\item[\key{orca\_r05}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~05
-\item[\key{orca\_r025}]  \jp{cp\_cfg}~=~orca ; \jp{jp\_cfg}~=~025
-\end{description}
-
-\subsubsection{ORCA mesh}
-
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_MISC_ORCA_msh}  \begin{center}
-\includegraphics[width=0.98\textwidth]{./TexFiles/Figures/Fig_ORCA_NH_mesh.pdf}
-\caption {ORCA mesh conception. The departure from an isotropic Mercator grid start poleward of 20\deg N.
-The two "north pole" are the foci of a series of embedded ellipses (blue curves) 
-which are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes). 
-Then, following \citet{Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed 
-which provide the j-lines of the mesh (pseudo longitudes).
- }
-\end{center}   \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tbp] \label{Fig_MISC_ORCA_e1e2}  \begin{center}
-\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_ORCA_NH_msh05_e1_e2.pdf}
-\includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_ORCA_aniso.pdf}
-\caption {\textit{Top}: Horizontal scale factors ($e_1$, $e_2$) and 
-\textit{Bottom}: ratio of anisotropy ($e_1 / e_2$)
-for ORCA 0.5\deg ~mesh. South of 20\deg N a Mercator grid is used ($e_1 = e_2$) 
-so that the anisotropy ratio is 1. Poleward of 20\deg N, the two "north pole" 
-introduce a weak anisotropy over the ocean areas ($< 1.2$) except in vicinity of Victoria Island 
-(Canadian Arctic Archipelago). }
-\end{center}   \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-%--------------------------------------------------TABLE--------------------------------------------------
-\begin{table}[htbp]  \label{Tab_ORCA}
-\begin{center}
-\begin{tabular}{ccccc}
-key          	             & \jp{jp\_cfg} &  \jp{jpiglo} & \jp{jpiglo} &       \\  
-\hline  \hline
-\key{orca\_r4} 	    &        4         &         92     &      76      &       \\
-\key{orca\_r2}       &        2         &       182     &    149      &        \\
-%\key{orca\_r1}       &        1         &       362     &     511     &        \\
-\key{orca\_r05}     &        05       &       722     &     261     &        \\
-\key{orca\_r025}   &        025     &      1442    &   1021     &        \\
-%\key{orca\_r8}       &        8         &      2882    &   2042     &        \\
-%\key{orca\_r12}     &      12         &      4322    &   3062      &       \\
-\hline
-\hline
-\end{tabular}
-\caption {Set of predefined ORCA parameters. }
-\end{center}
-\end{table}
-%--------------------------------------------------------------------------------------------------------------
-
-The tripolar grid used in ORCA configuration ....
-
-NB: the two north poles position has been chosen to minimise the anisotropy ratio in 
-the Gulf Stream and kuroshio areas, two highly turbulent regions.
-
-ORCA~2 : a $2\deg$ zonal resolution, and a meridional resolution varying from $0.5\deg$ at the 
-equator to $2\deg cos\phi$ south of $20\deg$S (Fig. 1). The grid features two points of convergence in the 
-Northern Hemisphere, both situated on continents. Minimum resolution in high latitudes is about 
-65~km in the Arctic and 50~km in the Antarctic. Local mesh refinements are applied to the 
-Mediterranean, Red, Black and Caspian Seas. None of them appears to be of particular 
-importance for the study of high latitude climate, but the fine resolution is needed in order to have 
-their local circulation and their role in the World Ocean's circulation considered correctly. 
-
-
- 
-ORCA2-LIM (global ocean sea-ice configuration \citep{Timmermann_al_OM05}. 
-The horizontal mesh is based on a $2\deg \times 2\deg$ Mercator grid ($i.e.$ same zonal and 
-meridional grid spacing) which has been modified poleward  
-of $20\deg$N in order to include two numerical inland poles \citep{Murray_JCP96}. 
-This modification is semi-analytical \citep{Madec_Imbard_CD96} 
-and based on a series of embedded ellipses. It insures that the mesh remains 
-close to isotropy and that the smallest grid cell is along Antarctica. 
-In order to refine the meridional resolution up to $0.5\deg$ at the equator, 
-additional local transformations were applied with in the Tropics. 
-Local mesh refinements are also applied to the Mediterranean, Red, Black 
-and Caspian Seas so that the resolution is $1\deg \time 1\deg$ there. 
-There are 31 levels in the vertical, with the highest resolution (10m) 
-in the upper 150m. The bottom topography and the coastlines are derived 
-from the global atlas of Smith and Sandwell (1997).
-
-\key{orca\_lev10} 10 time more vertical levels
-
-\key{agrif}  : ORCA2-LIM plus an AGRIF zoom over the Agulhas current area
-
-\key{arctic}, \key{antarctic}  (not used in ORCA\_R4)
-
-
-We thus only provide a brief introduction in this chapter. 
-The global coupled ocean-ice configuration is very similar to that used as part of the climate 
-model developed at GFDL for the 4th IPCC assessment of climate change (Griffies et al., 2005; 
-Gnanadesikan et al., 2006). 
-The ORCA2-LIM configuration is also the basis for the \NEMO contribution to the 
-Coordinate Ocean-ice Reference Experiments (COREs) documented in \citet{Griffies_al_OM09}. 
-These experiments employ the boundary forcing from \citet{Large_Yeager_Rep04} (see \S\ref{SBC_blk_core}), 
-which was developed for the purpose of running global coupled ocean-ice simulations without an 
-interactive atmosphere. This \citet{Large_Yeager_Rep04} dataset is available through the GFDL web 
-site \footnote{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}. 
-The "normal year" of \citet{Large_Yeager_Rep04} has been chosen of the \NEMO distribution 
-since release v3.3. 
-
-% -------------------------------------------------------------------------------------------------------------
-%       GYRE family configuration
-% -------------------------------------------------------------------------------------------------------------
-\subsection{GYRE family: double gyre basin (\key{gyre})}
-\label{MISC_config_gyre}
-
-The GYRE configuration \citep{Levy_al_OM10} have been built to simulated 
-the seasonal cycle of a double-gyre box model. It consist in an idealized domain 
-similar to that used in the studies of \citet{Drijfhout_JPO94} and \citet{Hazeleger_Drijfhout_JPO98, 
-Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00}, 
-over which an analytical seasonal forcing is applied. This allows to investigate the 
-spontaneous generation of a large number of interacting, transient mesoscale eddies 
-and their contribution to the large scale circulation. 
-
-The domain geometry is a closed rectangular basin on the $\beta$-plane centred 
-at $\sim 30\deg$N and rotated by 45\deg, 3180~km long, 2120~km wide 
-and 4~km deep (Fig.~\ref{Fig_MISC_strait_hand}). 
-The domain is bounded by vertical walls and by a ßat bottom. The configuration is 
-meant to represent an idealized North Atlantic or North Pacific basin. 
-The circulation is forced by analytical profiles of wind and buoyancy ßuxes. 
-The applied forcings vary seasonally in a sinusoidal manner between winter 
-and summer extrema \citep{Levy_al_OM10}. 
-The wind stress is zonal and its curl changes sign at 22\deg N and 36\deg N. 
-It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain 
-and a small recirculation gyre in the southern corner. 
-The net heat ßux takes the form of a restoring toward a zonal apparent air 
-temperature profile. A portion of the net heat ßux which comes from the solar radiation
-is allowed to penetrate within the water column. 
-The fresh water ßux is also prescribed and varies zonally. 
-It is determined such as, at each time step, the basin-integrated ßux is zero. 
-The basin is initialised at rest with vertical profiles of temperature and salinity 
-uniformly applied to the whole domain.
-
-The GYRE configuration is set through the \key{gyre} CPP key. Its horizontal resolution 
-(and thus the size of the domain) is determined by setting \jp{jp\_cfg} in \hf{par\_GYRE} file: \\
-\jp{jpiglo} $= 30 \times$ \jp{jp\_cfg} + 2   \\
-\jp{jpjglo} $= 20 \times$ \jp{jp\_cfg} + 2   \\
-Obviously, the namelist parameters have to be adjusted to the chosen resolution.
-In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}=31) (Fig.~\ref{Fig_zgr}).
-
-The GYRE configuration is also used in benchmark test as it is very simple to increase 
-its resolution and as it does not requires any input file. For example, keeping a same model size 
-on each processor while increasing the number of processor used is very easy, even though the 
-physical integrity of the solution can be compromised.
-
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_GYRE}  \begin{center}
-\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_GYRE.pdf}
-\caption {Snapshot of relative vorticity at the surface of the model domain 
-in GYRE R9, R27 and R54. From \citet{Levy_al_OM10}.}
-\end{center}   \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-% -------------------------------------------------------------------------------------------------------------
-%       EEL family configuration
-% -------------------------------------------------------------------------------------------------------------
-\subsection{EEL family: periodic channel}
-\label{MISC_config_EEL}
-
-\begin{description}
-\item[\key{eel\_r2}]  
-\item[\key{eel\_r5}]  
-\item[\key{eel\_r6}]  
-\end{description}
-
-% -------------------------------------------------------------------------------------------------------------
-%       POMME configuration
-% -------------------------------------------------------------------------------------------------------------
-\subsection{POMME: mid-latitude sub-domain}
-\label{MISC_config_POMME}
-
-
-\key{pomme\_r025} 
-
-
-
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_Model_Basics.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_Model_Basics.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_Model_Basics.tex	(revision 2376)
@@ -6,5 +6,4 @@
 \label{PE}
 \minitoc
-
 
 \newpage
@@ -114,9 +113,10 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!ht] \label{Fig_ocean_bc}  \begin{center}
+\begin{figure}[!ht]	 \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_I_ocean_bc.pdf}
-\caption{The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, where $H$ 
-is the depth of the sea floor and $\eta$ the height of the sea surface. Both $H$ and $\eta $ 
-are referenced to $z=0$.}
+\caption{	 \label{Fig_ocean_bc} 
+The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,t)$, where $H$ 
+is the depth of the sea floor and $\eta$ the height of the sea surface. 
+Both $H$ and $\eta$ are referenced to $z=0$.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -167,6 +167,6 @@
 
 
-\newpage
-$\ $\newline    % force a new ligne
+%\newpage
+%$\ $\newline    % force a new ligne
 
 % ================================================================
@@ -371,7 +371,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tb] \label{Fig_referential}  \begin{center}
+\begin{figure}[!tb] 	 \begin{center}
 \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_I_earth_referential.pdf}
-\caption{the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 
+\caption{	\label{Fig_referential} 
+the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 
 coordinate system (\textbf{i},\textbf{j},\textbf{k}). }
 \end{center}   \end{figure}
@@ -703,5 +704,4 @@
 \label{PE_gco}
 
-%\gmcomment{
 The ocean domain presents a huge diversity of situation in the vertical. First the ocean surface is a time dependent surface (moving surface). Second the ocean floor depends on the geographical position, varying from more than 6,000 meters in abyssal trenches to zero at the coast. Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 
 Therefore, in order to represent the ocean with respect to the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height $e.g.$ an $z$*-coordinate; for the second point, a space variation to fit the change of bottom topography $e.g.$ a terrain-following or $\sigma$-coordinate; and for the third point, one will be tempted to use a space and time dependent coordinate that follows the isopycnal surfaces, $e.g.$ an isopycnic coordinate.
@@ -717,15 +717,35 @@
 The coordinate is also sometime referenced as an adaptive coordinate \citep{Hofmeister_al_OM09}, since the coordinate system is adapted in the course of the simulation. Its most often used implementation is via an ALE algorithm, in which a pure lagrangian step is followed by regridding and remapping steps, the later step implicitly embedding the vertical advection \citep{Hirt_al_JCP74, Chassignet_al_JPO03, White_al_JCP09}. Here we follow the \citep{Kasahara_MWR74} strategy : a regridding step (an update of the vertical coordinate) followed by an eulerian step with an explicit computation of vertical advection relative to the moving s-surfaces.
 
-A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient...
-
-the generalized vertical coordinates used in ocean modelling are not orthogonal, which contrasts with many other applications in mathematical physics. Hence, it is useful to keep in mind the following properties that may seem odd on initial encounter.
-
-the horizontal velocity in ocean models measures motions in the horizontal plane, perpendicular to the local gravitational field. That is, horizontal velocity is mathematically the same regardless the vertical coordinate, be it geopotential, isopycnal, pressure, or terrain following. The key motivation for maintaining the same horizontal velocity component is that the hydrostatic and geostrophic balances are dominant in the large-scale ocean. Use of an alternative quasi-horizontal velocity, for example one oriented parallel to the generalized surface, would lead to unacceptable numerical errors. Correspondingly, the vertical direction is anti-parallel to the gravitational force in all of the coordinate systems. We do not choose the alternative of a quasi-vertical direction oriented normal to the surface of a constant generalized vertical coordinate. 
-
-It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between the vertical coordinate choices. That is, computation of the dia-surface velocity component represents the fundamental distinction between the various coordinates. In some models, such as geopotential, pressure, 
-and terrain following, this transport is typically diagnosed from volume or mass conservation. In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about the physical processes producing a flux across the layer interfaces. 
-
-
-In this section we first establish the PE in the generalised vertical $s$-coordinate, then we discuss the particular cases available in \NEMO, namely $z$, $z$*, $s$, and $\tilde z$.  
+%\gmcomment{
+
+%A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient...
+
+the generalized vertical coordinates used in ocean modelling are not orthogonal, 
+which contrasts with many other applications in mathematical physics. 
+Hence, it is useful to keep in mind the following properties that may seem 
+odd on initial encounter.
+
+The horizontal velocity in ocean models measures motions in the horizontal plane, 
+perpendicular to the local gravitational field. That is, horizontal velocity is mathematically 
+the same regardless the vertical coordinate, be it geopotential, isopycnal, pressure, 
+or terrain following. The key motivation for maintaining the same horizontal velocity 
+component is that the hydrostatic and geostrophic balances are dominant in the large-scale ocean. 
+Use of an alternative quasi-horizontal velocity, for example one oriented parallel 
+to the generalized surface, would lead to unacceptable numerical errors. 
+Correspondingly, the vertical direction is anti-parallel to the gravitational force in all 
+of the coordinate systems. We do not choose the alternative of a quasi-vertical 
+direction oriented normal to the surface of a constant generalized vertical coordinate. 
+
+It is the method used to measure transport across the generalized vertical coordinate 
+surfaces which differs between the vertical coordinate choices. That is, computation 
+of the dia-surface velocity component represents the fundamental distinction between 
+the various coordinates. In some models, such as geopotential, pressure, and 
+terrain following, this transport is typically diagnosed from volume or mass conservation. 
+In other models, such as isopycnal layered models, this transport is prescribed based 
+on assumptions about the physical processes producing a flux across the layer interfaces. 
+
+
+In this section we first establish the PE in the generalised vertical $s$-coordinate, 
+then we discuss the particular cases available in \NEMO, namely $z$, $z$*, $s$, and $\tilde z$.  
 %}
 
@@ -821,8 +841,10 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!b] \label{Fig_z_zstar}  \begin{center}
+\begin{figure}[!b] 	 \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zstar.pdf}
-\caption{(a) $z$-coordinate in linear free-surface case ; (b) $z-$coordinate in non-linear 
-free surface case (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate 
+\caption{	\label{Fig_z_zstar} 
+(a) $z$-coordinate in linear free-surface case ; 
+(b) $z-$coordinate in non-linear free surface case ; 
+(c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate 
 \citep{Adcroft_Campin_OM04} ).}
 \end{center}   \end{figure}
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex	(revision 2376)
@@ -117,8 +117,19 @@
 
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
-\begin{figure}[!t] \label{Fig_DYN_dynspg_ts}
-\begin{center}
+\begin{figure}[!t] 	\begin{center}
 \includegraphics[width=0.90\textwidth]{./Figures/Fig_DYN_dynspg_ts.pdf}
-\caption{Schematic of the split-explicit time stepping scheme for the barotropic and baroclinic modes, after \citet{Griffies2004}. Time increases to the right. Baroclinic time steps are denoted by $t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$. The curved line represents a leap-frog time step, and the smaller barotropic time steps $N \Delta t=2\Delta t$ are denoted by the zig-zag line. The vertically integrated forcing \textbf{M}(t) computed at baroclinic time step t represents the interaction between the barotropic and baroclinic motions. While keeping the total depth, tracer, and freshwater forcing fields fixed, a leap-frog integration carries the surface height and vertically integrated velocity from t to $t+2 \Delta t$ using N barotropic time steps of length $\Delta t$. Time averaging the barotropic fields over the N+1 time steps (endpoints included) centers the vertically integrated velocity at the baroclinic timestep $t+\Delta t$. A baroclinic leap-frog time step carries the surface height to $t+\Delta t$ using the convergence of the time averaged vertically integrated velocity taken from baroclinic time step t. }
+\caption{	 \label{Fig_DYN_dynspg_ts}
+Schematic of the split-explicit time stepping scheme for the barotropic and baroclinic modes, 
+after \citet{Griffies2004}. Time increases to the right. Baroclinic time steps are denoted by 
+$t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$. The curved line represents a leap-frog time step, 
+and the smaller barotropic time steps $N \Delta t=2\Delta t$ are denoted by the zig-zag line. 
+The vertically integrated forcing \textbf{M}(t) computed at baroclinic time step t represents 
+the interaction between the barotropic and baroclinic motions. While keeping the total depth, 
+tracer, and freshwater forcing fields fixed, a leap-frog integration carries the surface height 
+and vertically integrated velocity from t to $t+2 \Delta t$ using N barotropic time steps of length 
+$\Delta t$. Time averaging the barotropic fields over the N+1 time steps (endpoints included) 
+centers the vertically integrated velocity at the baroclinic timestep $t+\Delta t$. 
+A baroclinic leap-frog time step carries the surface height to $t+\Delta t$ using the convergence 
+of the time averaged vertically integrated velocity taken from baroclinic time step t. }
 \end{center}
 \end{figure}
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_OBS.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_OBS.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_OBS.tex	(revision 2376)
@@ -678,12 +678,11 @@
 \subsubsection{Geographical distribution of observations among processors}
 
-\begin{figure}
-\begin{center}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure} 	\begin{center}
 \includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_ASM_obsdist_local}
-\end{center}
-\caption{Example of the distribution of observations with the geographical
-distribution of observational data.} 
-\label{fig:obslocal}
-\end{figure}
+\caption{ 	\label{fig:obslocal}
+Example of the distribution of observations with the geographical distribution of observational data.} 
+\end{center} 	\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 This is the simplest option in which the observations are distributed according 
@@ -706,12 +705,11 @@
 \subsubsection{Round-robin distribution of observations among processors}
 
-\begin{figure}
-\begin{center}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}	\begin{center}
 \includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_ASM_obsdist_global}
-\end{center}
-\caption{Example of the distribution of observations with the round-robin
-  distribution of observational data.}
-\label{fig:obsglobal}
-\end{figure}
+\caption{ 	\label{fig:obsglobal}
+Example of the distribution of observations with the round-robin distribution of observational data.}
+\end{center}	\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 An alternative approach is to distribute the observations equally
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2376)
@@ -29,9 +29,8 @@
 (\np{ln\_core}~=~true) or CLIO (\np{ln\_clio}~=~true) bulk formulae) and a coupled 
 formulation (exchanges with a atmospheric model via the OASIS coupler) 
-(\np{ln\_cpl}~=~true). The optional atmospheric pressure can be used either 
-to force ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true), or in the bulk 
-formulae computation (\np{ln\_apr\_dyn}~=~true)
-\footnote{None of the two current bulk formulea (CLIO and CORE) uses the 
-atmospheric pressure field.}. 
+(\np{ln\_cpl}~=~true). When used, the atmospheric pressure forces both 
+ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true)
+\footnote{The surface pressure field could be use in bulk formulae, nevertheless 
+none of the current bulk formulea (CLIO and CORE) uses the it.}. 
 The frequency at which the six or seven fields have to be updated is the \np{nn\_fsbc} 
 namelist parameter. 
@@ -39,5 +38,5 @@
 need not be supplied on the model grid.  Instead a file of coordinates and weights can 
 be supplied which maps the data from the supplied grid to the model points 
-(so called "Interpolation on the Fly").
+(so called "Interpolation on the Fly", see \S\ref{SBC_iof}).
 In addition, the resulting fields can be further modified using several namelist options. 
 These options control  the rotation of vector components supplied relative to an east-north 
@@ -52,5 +51,6 @@
 
 In this chapter, we first discuss where the surface boundary condition appears in the
-model equations. Then we present the four ways of providing the surface boundary condition. 
+model equations. Then we present the four ways of providing the surface boundary condition, 
+followed by the description of the atmospheric pressure and the river runoff. 
 Next the scheme for interpolation on the fly is described.
 Finally, the different options that further modify the fluxes applied to the ocean are discussed.
@@ -157,7 +157,5 @@
 
 %-------------------------------------------------TABLE---------------------------------------------------
-\begin{table}[tb]  \label{Tab_ssm}
-\begin{center}
-\begin{tabular}{|l|l|l|l|}
+\begin{table}[tb]   \begin{center}   \begin{tabular}{|l|l|l|l|}
 \hline
 Variable description					& Model variable	& Units	& point \\	\hline
@@ -167,14 +165,12 @@
 Sea surface salinty					& sss\_m	& $psu$			& T \\	\hline
 \end{tabular}
-\caption{Ocean variables provided by the ocean to the surface module (SBC). 
+\caption{  \label{Tab_ssm}   
+Ocean variables provided by the ocean to the surface module (SBC). 
 The variable are averaged over nf{\_}sbc time step, $i.e.$ the frequency of 
 computation of surface fluxes.}
-\end{center}
-\end{table}
+\end{center}   \end{table}
 %--------------------------------------------------------------------------------------------------------------
 
-
-
-%\colorbox{yellow}{Penser a} mettre dans le restant l'info nf{\_}sbc ET nf{\_}sbc*rdt de sorte de reinitialiser la moyenne si on change la frequence ou le pdt
+%\colorbox{yellow}{Penser a} mettre dans le restant l'info nn{\_}fsbc ET nn{\_}fsbc*rdt de sorte de reinitialiser la moyenne si on change la frequence ou le pdt
 
 
@@ -385,9 +381,11 @@
 %-------------------------------------------------------------------------------------------------------------
 
-The optional atmospheric pressure can be used either to force ocean and ice dynamics 
-(\np{ln\_apr\_dyn}~=~true), or in the bulk formulae computation (\np{ln\_apr\_dyn}~=~true).
-The input atmospheric forcing is interpolated in time to the model time step, and optionally 
-in space when interpolation on-the-fly is used. When used to force the dynamics, it is further 
-transformed into an equivalent inverse barometer sea surface height, $\eta_{ib}$, using:
+The optional atmospheric pressure can be used to force ocean and ice dynamics 
+(\np{ln\_apr\_dyn}~=~true, \textit{namsbc} namelist ).
+The input atmospheric forcing defined via \np{sn\_apr} structure (\textit{namsbc\_apr} namelist) 
+can be interpolated in time to the model time step, and even in space when the 
+interpolation on-the-fly is used. When used to force the dynamics, the atmospheric 
+pressure is further transformed into an equivalent inverse barometer sea surface height, 
+$\eta_{ib}$, using:
 \begin{equation} \label{SBC_ssh_ib}
 	\eta_{ib} = -  \frac{1}{g\,\rho_o}  \left( P_{atm} - P_o \right) 
@@ -398,10 +396,9 @@
 $\eta_{ib}$ is kept to zero at all time step.
 
-A gradient of $\eta_{ib}$ is added to the RHS of the ocean momentum equation 
+The gradient of $\eta_{ib}$ is added to the RHS of the ocean momentum equation 
 (see \mdl{dynspg} for the ocean). For sea-ice, the sea surface height, $\eta_m$, 
 which is provided to the sea ice model is set to $\eta - \eta_{ib}$ (see \mdl{sbcssr} module).
-Furthermore, $\eta_{ib}$ can be set in the output. This simplifies the altirmetry data 
-and model comparison as inverse barometer sea surface height is usually removed 
-from thise date prior to their distribution.
+$\eta_{ib}$ can be set in the output. This can simplify the altirmetry data and model comparison 
+as inverse barometer sea surface height is usually removed from these date prior to their distribution.
 
 % ================================================================
@@ -522,8 +519,141 @@
 
 
+% ================================================================
+% Interpolation on the Fly
+% ================================================================
+
+\section [Interpolation on the Fly] {Interpolation on the Fly}
+\label{SBC_iof}
+
+Interpolation on the Fly allows the user to supply input files required
+for the surface forcing on grids other than the model grid.
+To do this he or she must supply, in addition to the source data file,
+a file of weights to be used to interpolate from the data grid to the model
+grid.
+The original development of this code used the SCRIP package (freely available 
+under a copyright agreement from http://climate.lanl.gov/Software/SCRIP).
+In principle, any package can be used to generate the weights, but the
+variables in the input weights file must have the same names and meanings as
+assumed by the model.
+Two methods are currently available: bilinear and bicubic interpolation.
+
+\subsection{Bilinear Interpolation}
+\label{SBC_iof_bilinear}
+
+The input weights file in this case has two sets of variables: src01, src02,
+src03, src04 and wgt01, wgt02, wgt03, wgt04.
+The "src" variables correspond to the point in the input grid to which the weight
+"wgt" is to be applied. Each src value is an integer corresponding to the index of a 
+point in the input grid when written as a one dimensional array.  For example, for an input grid
+of size 5x10, point (3,2) is referenced as point 8, since (2-1)*5+3=8.
+There are four of each variable because bilinear interpolation uses the four points defining
+the grid box containing the point to be interpolated.
+All of these arrays are on the model grid, so that values src01(i,j) and
+wgt01(i,j) are used to generate a value for point (i,j) in the model.
+
+Symbolically, the algorithm used is:
+
+\begin{equation}
+f_{m}(i,j) = f_{m}(i,j) + \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}
+\end{equation}
+where function idx() transforms a one dimensional index src(k) into a two dimensional index,
+and wgt(1) corresponds to variable "wgt01" for example.
+
+\subsection{Bicubic Interpolation}
+\label{SBC_iof_bicubic}
+
+Again there are two sets of variables: "src" and "wgt".
+But in this case there are 16 of each.
+The symbolic algorithm used to calculate values on the model grid is now:
+
+\begin{equation*} \begin{split}
+f_{m}(i,j) =  f_{m}(i,j) +& \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}     
+              +   \sum_{k=5}^{8} {wgt(k)\left.\frac{\partial f}{\partial i}\right| _{idx(src(k))} }    \\
+              +& \sum_{k=9}^{12} {wgt(k)\left.\frac{\partial f}{\partial j}\right| _{idx(src(k))} }   
+              +   \sum_{k=13}^{16} {wgt(k)\left.\frac{\partial ^2 f}{\partial i \partial j}\right| _{idx(src(k))} }
+\end{split}
+\end{equation*}
+The gradients here are taken with respect to the horizontal indices and not distances since the spatial dependency has been absorbed into the weights.
+
+\subsection{Implementation}
+\label{SBC_iof_imp}
+
+To activate this option, a non-empty string should be supplied in the weights filename column 
+of the relevant namelist; if this is left as an empty string no action is taken.
+In the model, weights files are read in and stored in a structured type (WGT) in the fldread 
+module, as and when they are first required.
+This initialisation procedure determines whether the input data grid should be treated 
+as cyclical or not by inspecting a global attribute stored in the weights input file.
+This attribute must be called "ew\_wrap" and be of integer type.
+If it is negative, the input non-model grid is assumed not to be cyclic.
+If zero or greater, then the value represents the number of columns that overlap.
+$E.g.$ if the input grid has columns at longitudes 0, 1, 2, .... , 359, then ew\_wrap should be set to 0;
+if longitudes are 0.5, 2.5, .... , 358.5, 360.5, 362.5, ew\_wrap should be 2.
+If the model does not find attribute ew\_wrap, then a value of -999 is assumed.
+In this case the \rou{fld\_read} routine defaults ew\_wrap to value 0 and therefore the grid 
+is assumed to be cyclic with no overlapping columns.
+(In fact this only matters when bicubic interpolation is required.)
+Note that no testing is done to check the validity in the model, since there is no way 
+of knowing the name used for the longitude variable,
+so it is up to the user to make sure his or her data is correctly represented.
+
+Next the routine reads in the weights.
+Bicubic interpolation is assumed if it finds a variable with name "src05", otherwise 
+bilinear interpolation is used. The WGT structure includes dynamic arrays both for 
+the storage of the weights (on the model grid), and when required, for reading in 
+the variable to be interpolated (on the input data grid).
+The size of the input data array is determined by examining the values in the "src" 
+arrays to find the minimum and maximum i and j values required.
+Since bicubic interpolation requires the calculation of gradients at each point on the grid, 
+the corresponding arrays are dimensioned with a halo of width one grid point all the way around.
+When the array of points from the data file is adjacent to an edge of the data grid, 
+the halo is either a copy of the row/column next to it (non-cyclical case), or is a copy 
+of one from the first few columns on the opposite side of the grid (cyclical case).
+
+\subsection{Limitations}
+\label{SBC_iof_lim}
+
+\begin{description}
+\item
+The case where input data grids are not logically rectangular has not been tested.
+\item
+This code is not guaranteed to produce positive definite answers from positive definite inputs.
+\item
+The cyclic condition is only applied on left and right columns, and not to top and bottom rows.
+\item
+The gradients across the ends of a cyclical grid assume that the grid spacing between the two columns involved are consistent with the weights used.
+\item
+Neither interpolation scheme is conservative.
+(There is a conservative scheme available in SCRIP, but this has not been implemented.)
+\end{description}
+
+\subsection{Utilities}
+\label{SBC_iof_util}
+
+% to be completed
+A set of utilities to create a weights file for a rectilinear input grid is available 
+(see the directory NEMOGCM/TOOLS/WEIGHTS).
+
+% ================================================================
+% Miscellanea options
+% ================================================================
+\section{Miscellaneous options}
+\label{SBC_misc}
+
+% -------------------------------------------------------------------------------------------------------------
+%        Diurnal cycle
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Diurnal  cycle (\textit{sbcdcy})]
+			{Diurnal cycle (\mdl{sbcdcy})}
+\label{SBC_dcy}
+%------------------------------------------namsbc_rnf----------------------------------------------------
+%\namdisplay{namsbc} 
+%-------------------------------------------------------------------------------------------------------------
+
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_SBC_diurnal}  \begin{center}
+\begin{figure}[!t]    \begin{center}
 \includegraphics[width=0.8\textwidth]{./TexFiles/Figures/Fig_SBC_diurnal.pdf}
-\caption{Example of recontruction of the diurnal cycle variation of short wave flux  
+\caption{ \label{Fig_SBC_diurnal}    
+Example of recontruction of the diurnal cycle variation of short wave flux  
 from daily mean values. The reconstructed diurnal cycle (black line) is chosen 
 as the mean value of the analytical cycle (blue line) over a time step, not 
@@ -531,14 +661,4 @@
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-% ================================================================
-%        Diurnal cycle
-% ================================================================
-\section   [Diurnal  cycle (\textit{sbcdcy})]
-			{Diurnal cycle (\mdl{sbcdcy})}
-\label{SBC_dcy}
-%------------------------------------------namsbc_rnf----------------------------------------------------
-%\namdisplay{namsbc} 
-%-------------------------------------------------------------------------------------------------------------
 
 \cite{Bernie_al_JC05} have shown that to capture 90$\%$ of the diurnal variability of 
@@ -565,7 +685,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_SBC_dcy}  \begin{center}
+\begin{figure}[!t]  \begin{center}
 \includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_SBC_dcy.pdf}
-\caption{Example of recontruction of the diurnal cycle variation of short wave flux  
+\caption{ \label{Fig_SBC_dcy}   
+Example of recontruction of the diurnal cycle variation of short wave flux  
 from daily mean values on an ORCA2 grid with a time sampling of 2~hours (from 1am to 11pm). 
 The display is on (i,j) plane. }
@@ -577,125 +698,4 @@
 appear due to an inconsistency between the scale of the vertical resolution 
 and the forcing acting on that scale.
-
-% ================================================================
-% Interpolation on the Fly
-% ================================================================
-
-\section [Interpolation on the Fly] {Interpolation on the Fly}
-\label{SBC_iof}
-
-Interpolation on the Fly allows the user to supply input files required
-for the surface forcing on grids other than the model grid.
-To do this he or she must supply, in addition to the source data file,
-a file of weights to be used to interpolate from the data grid to the model
-grid.
-The original development of this code used the SCRIP package (freely available 
-under a copyright agreement from http://climate.lanl.gov/Software/SCRIP).
-In principle, any package can be used to generate the weights, but the
-variables in the input weights file must have the same names and meanings as
-assumed by the model.
-Two methods are currently available: bilinear and bicubic interpolation.
-
-\subsection{Bilinear Interpolation}
-\label{SBC_iof_bilinear}
-
-The input weights file in this case has two sets of variables: src01, src02,
-src03, src04 and wgt01, wgt02, wgt03, wgt04.
-The "src" variables correspond to the point in the input grid to which the weight
-"wgt" is to be applied. Each src value is an integer corresponding to the index of a 
-point in the input grid when written as a one dimensional array.  For example, for an input grid
-of size 5x10, point (3,2) is referenced as point 8, since (2-1)*5+3=8.
-There are four of each variable because bilinear interpolation uses the four points defining
-the grid box containing the point to be interpolated.
-All of these arrays are on the model grid, so that values src01(i,j) and
-wgt01(i,j) are used to generate a value for point (i,j) in the model.
-
-Symbolically, the algorithm used is:
-
-\begin{equation}
-f_{m}(i,j) = f_{m}(i,j) + \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}
-\end{equation}
-where function idx() transforms a one dimensional index src(k) into a two dimensional index,
-and wgt(1) corresponds to variable "wgt01" for example.
-
-\subsection{Bicubic Interpolation}
-\label{SBC_iof_bicubic}
-
-Again there are two sets of variables: "src" and "wgt".
-But in this case there are 16 of each.
-The symbolic algorithm used to calculate values on the model grid is now:
-
-\begin{equation*} \begin{split}
-f_{m}(i,j) =  f_{m}(i,j) +& \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}     
-              +   \sum_{k=5}^{8} {wgt(k)\left.\frac{\partial f}{\partial i}\right| _{idx(src(k))} }    \\
-              +& \sum_{k=9}^{12} {wgt(k)\left.\frac{\partial f}{\partial j}\right| _{idx(src(k))} }   
-              +   \sum_{k=13}^{16} {wgt(k)\left.\frac{\partial ^2 f}{\partial i \partial j}\right| _{idx(src(k))} }
-\end{split}
-\end{equation*}
-The gradients here are taken with respect to the horizontal indices and not distances since the spatial dependency has been absorbed into the weights.
-
-\subsection{Implementation}
-\label{SBC_iof_imp}
-
-To activate this option, a non-empty string should be supplied in the weights filename column 
-of the relevant namelist; if this is left as an empty string no action is taken.
-In the model, weights files are read in and stored in a structured type (WGT) in the fldread 
-module, as and when they are first required.
-This initialisation procedure determines whether the input data grid should be treated 
-as cyclical or not by inspecting a global attribute stored in the weights input file.
-This attribute must be called "ew\_wrap" and be of integer type.
-If it is negative, the input non-model grid is assumed not to be cyclic.
-If zero or greater, then the value represents the number of columns that overlap.
-$E.g.$ if the input grid has columns at longitudes 0, 1, 2, .... , 359, then ew\_wrap should be set to 0;
-if longitudes are 0.5, 2.5, .... , 358.5, 360.5, 362.5, ew\_wrap should be 2.
-If the model does not find attribute ew\_wrap, then a value of -999 is assumed.
-In this case the \rou{fld\_read} routine defaults ew\_wrap to value 0 and therefore the grid 
-is assumed to be cyclic with no overlapping columns.
-(In fact this only matters when bicubic interpolation is required.)
-Note that no testing is done to check the validity in the model, since there is no way 
-of knowing the name used for the longitude variable,
-so it is up to the user to make sure his or her data is correctly represented.
-
-Next the routine reads in the weights.
-Bicubic interpolation is assumed if it finds a variable with name "src05", otherwise 
-bilinear interpolation is used. The WGT structure includes dynamic arrays both for 
-the storage of the weights (on the model grid), and when required, for reading in 
-the variable to be interpolated (on the input data grid).
-The size of the input data array is determined by examining the values in the "src" 
-arrays to find the minimum and maximum i and j values required.
-Since bicubic interpolation requires the calculation of gradients at each point on the grid, 
-the corresponding arrays are dimensioned with a halo of width one grid point all the way around.
-When the array of points from the data file is adjacent to an edge of the data grid, 
-the halo is either a copy of the row/column next to it (non-cyclical case), or is a copy 
-of one from the first few columns on the opposite side of the grid (cyclical case).
-
-\subsection{Limitations}
-\label{SBC_iof_lim}
-
-\begin{description}
-\item
-The case where input data grids are not logically rectangular has not been tested.
-\item
-This code is not guaranteed to produce positive definite answers from positive definite inputs.
-\item
-The cyclic condition is only applied on left and right columns, and not to top and bottom rows.
-\item
-The gradients across the ends of a cyclical grid assume that the grid spacing between the two columns involved are consistent with the weights used.
-\item
-Neither interpolation scheme is conservative.
-(There is a conservative scheme available in SCRIP, but this has not been implemented.)
-\end{description}
-
-\subsection{Utilities}
-\label{SBC_iof_util}
-
-% to be completed
-A set of utilities to create a weights file for a rectilinear input grid is available.
-
-% ================================================================
-% Miscellanea options
-% ================================================================
-\section{Miscellaneous options}
-\label{SBC_misc}
 
 % -------------------------------------------------------------------------------------------------------------
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_STP.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_STP.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_STP.tex	(revision 2376)
@@ -195,7 +195,8 @@
 %\gmcomment{ 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_TimeStep_flowchart}  \begin{center}
+\begin{figure}[!t] 	  \begin{center}
 \includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_TimeStepping_flowchart.pdf}
-\caption{Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}. 
+\caption{ 	\label{Fig_TimeStep_flowchart}
+Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}. 
 The use of a semi-implicit computation of the hydrostatic pressure gradient requires
 the tracer equation to be stepped forward prior to the momentum equation. 
@@ -286,7 +287,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_MLF_forcing}  \begin{center}
+\begin{figure}[!t] 	  \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_MLF_forcing.pdf}
-\caption{Illustration of forcing integration methods. 
+\caption{ 	\label{Fig_MLF_forcing}
+Illustration of forcing integration methods. 
 (top) ''Traditional'' formulation : the forcing is defined at the same time as the variable 
 to which it is applied (integer value of the time step index) and it is applied over a $2\rdt$ period. 
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2376)
@@ -90,7 +90,8 @@
 the continuity equation which is used to calculate the vertical velocity.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_adv_scheme}  \begin{center}
+\begin{figure}[!t] 	 \begin{center}
 \includegraphics[width=0.9\textwidth]{./TexFiles/Figures/Fig_adv_scheme.pdf}
-\caption{Schematic representation of some ways used to evaluate the tracer value 
+\caption{	\label{Fig_adv_scheme} 
+Schematic representation of some ways used to evaluate the tracer value 
 at $u$-point and the amount of tracer exchanged between two neighbouring grid 
 points. Upsteam biased scheme (ups): the upstream value is used and the black 
@@ -836,9 +837,10 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_traqsr_irradiance}  \begin{center}
+\begin{figure}[!t] 	  \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_Irradiance.pdf}
-\caption{Penetration profile of the downward solar irradiance 
-calculated by four models. Two waveband chlorophyll-independent formulation (blue), 
-a chlorophyll-dependent monochromatic formulation (green), 4 waveband RGB formulation (red), 
+\caption{ 	 \label{Fig_traqsr_irradiance}
+Penetration profile of the downward solar irradiance calculated by four models. 
+Two waveband chlorophyll-independent formulation (blue), a chlorophyll-dependent 
+monochromatic formulation (green), 4 waveband RGB formulation (red), 
 61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of 
 (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. From \citet{Lengaigne_al_CD07}.}
@@ -856,7 +858,8 @@
 %--------------------------------------------------------------------------------------------------------------
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_geothermal}  \begin{center}
+\begin{figure}[!t] 	  \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_geoth.pdf}
-\caption{Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
+\caption{ 	\label{Fig_geothermal}
+Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
 It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}.}
 \end{center}   \end{figure}
@@ -963,7 +966,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_bbl}  \begin{center}
+\begin{figure}[!t] 	\begin{center}
 \includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf}
-\caption{Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is 
+\caption{ 	\label{Fig_bbl}  
+Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is 
 activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$. 
 Red arrows indicate the additional overturning circulation due to the advective BBL. 
@@ -1316,7 +1320,13 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!p] \label{Fig_Partial_step_scheme}  \begin{center}
+\begin{figure}[!p] 	 \begin{center}
 \includegraphics[width=0.9\textwidth]{./TexFiles/Figures/Partial_step_scheme.pdf}
-\caption{ Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate (\np{ln\_zps}=true) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. The horizontal difference is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$.  }
+\caption{ 	\label{Fig_Partial_step_scheme} 
+Discretisation of the horizontal difference and average of tracers in the $z$-partial 
+step coordinate (\np{ln\_zps}=true) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. 
+A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value 
+at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 
+The horizontal difference is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ 
+and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$.  }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2376)
@@ -154,4 +154,23 @@
 The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist parameter.
 
+At the sea surface, the value of $\bar{e}$ is prescribed from the wind 
+stress field as $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn\_ebb} 
+namelist parameter. The default value of $e_{bb}$ is 3.75. \citep{Gaspar1990}), 
+however a much larger value can be used when taking into account the 
+surface wave breaking (see below Eq. \eqref{ZDF_Esbc}). 
+The bottom value of TKE is assumed to be equal to the value of the level just above. 
+The time integration of the $\bar{e}$ equation may formally lead to negative values 
+because the numerical scheme does not ensure its positivity. To overcome this 
+problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} 
+namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set 
+to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations 
+to match that of \citet{Gargett1984} for the diffusion in the thermocline and 
+deep ocean :  $K_\rho = 10^{-3} / N$. 
+In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical 
+instabilities associated with too weak vertical diffusion. They must be 
+specified at least larger than the molecular values, and are set through 
+\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
+
+\subsubsection{Turbulent length scale}
 For computational efficiency, the original formulation of the turbulent length 
 scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 
@@ -187,7 +206,8 @@
 mixing length scales as (and note that here we use numerical indexing):
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_mixing_length}  \begin{center}
+\begin{figure}[!t] \begin{center}
 \includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_mixing_length.pdf}
-\caption {Illustration of the mixing length computation. }
+\caption{ \label{Fig_mixing_length} 
+Illustration of the mixing length computation. }
 \end{center}  
 \end{figure}
@@ -204,7 +224,7 @@
 $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
 
-In the \np{nn\_mxl}=2 case, the dissipation and mixing length scales take the same 
+In the \np{nn\_mxl}~=~2 case, the dissipation and mixing length scales take the same 
 value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the 
-\np{nn\_mxl}=2 case, the dissipation and mixing turbulent length scales are give 
+\np{nn\_mxl}~=~3 case, the dissipation and mixing turbulent length scales are give 
 as in \citet{Gaspar1990}:
 \begin{equation} \label{Eq_tke_mxl_gaspar}
@@ -215,23 +235,102 @@
 \end{equation}
 
-At the sea surface the value of $\bar{e}$ is prescribed from the wind 
-stress field: $\bar{e}=rn\_ebb\;\left| \tau \right|$ (\np{rn\_ebb}=60 by default) 
-with a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist
-parameters). Its value at the bottom of the ocean is assumed to be 
-equal to the value of the level just above. The time integration of the 
-$\bar{e}$ equation may formally lead to negative values because the 
-numerical scheme does not ensure its positivity. To overcome this 
-problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} 
-namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set 
-to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations 
-to match that of \citet{Gargett1984} for the diffusion in the thermocline and 
-deep ocean :  $K_\rho = 10^{-3} / N$. 
-In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical 
-instabilities associated with too weak vertical diffusion. They must be 
-specified at least larger than the molecular values, and are set through 
-\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
-
-% -------------------------------------------------------------------------------------------------------------
-%        TKE Turbulent Closure Scheme : new organization to energetic considerations
+At the ocean surface, a non zero length scale is set through the  \np{rn\_lmin0} namelist 
+parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$ 
+where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness 
+parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} 
+leads to a 0.04~m, the default value of \np{rn\_lsurf}. In the ocean interior 
+a minimum length scale is set to recover the molecular viscosity when $\bar{e}$ 
+reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
+
+
+\subsubsection{Surface wave breaking parameterization}
+%-----------------------------------------------------------------------%
+
+Following \citet{Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified 
+to include the effect of surface wave breaking energetics. This results in a reduction of summertime 
+surface temperature when the mixed layer is relatively shallow. The \citet{Mellor_Blumberg_JPO04} 
+modifications acts on surface length scale and TKE values and air-sea drag coefficient. 
+The latter concerns the bulk formulea and is not discussed here. 
+
+Following \citet{Craig_Banner_JPO94}, the boundary condition on surface TKE value is :
+\begin{equation}  \label{ZDF_Esbc}
+\bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
+\end{equation}
+where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of proportionality 
+which depends on the ''wave age'', ranging from 57 for mature waves to 146 for 
+younger waves \citep{Mellor_Blumberg_JPO04}. 
+The boundary condition on the turbulent length scale follows the Charnock's relation:
+\begin{equation} \label{ZDF_Lsbc}
+l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
+\end{equation}
+where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant.
+\citet{Mellor_Blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by \citet{Stacey_JPO99}
+citing observation evidence, and $\alpha_{CB} = 100$ the Craig and Banner's value.
+As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 
+with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds 
+to $\alpha_{CB} = 100$. further setting  \np{ln\_lsurf} to true applies \eqref{ZDF_Lsbc} 
+as surface boundary condition on length scale, with $\beta$ hard coded to the Stacet's value.
+Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) 
+is applied on surface $\bar{e}$ value.
+
+
+\subsubsection{Langmuir cells}
+%--------------------------------------%
+Langmuir circulations (LC) can be described as ordered large-scale vertical motions 
+in the surface layer of the oceans. Although LC have nothing to do with convection, 
+the circulation pattern is rather similar to so-called convective rolls in the atmospheric 
+boundary layer. The detailed physics behind LC is described in, for example, 
+\citet{Craik_Leibovich_JFM76}. The prevailing explanation is that LC arise from 
+a nonlinear interaction between the Stokes drift and wind drift currents. 
+
+Here we introduced in the TKE turbulent closure the simple parameterization of 
+Langmuir circulations proposed by \citep{Axell_JGR02} for a $k-\epsilon$ turbulent closure. 
+The parameterization, tuned against large-eddy simulation, includes the whole effect
+of LC in an extra source terms of TKE, $P_{LC}$.
+The presence of $P_{LC}$ in \eqref{Eq_zdftke_e}, the TKE equation, is controlled 
+by setting \np{ln\_lc} to \textit{true} in the namtke namelist.
+ 
+By making an analogy with the characteristic convective velocity scale 
+($e.g.$, \citet{D'Alessio_al_JPO98}), $P_{LC}$ is assumed to be : 
+\begin{equation}
+P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}}
+\end{equation}
+where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth.
+With no information about the wave field, $w_{LC}$ is assumed to be proportional to 
+the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module 
+\footnote{Following \citet{Li_Garrett_JMR93}, the surface Stoke drift velocity
+may be expressed as $u_s =  0.016 \,|U_{10m}|$. Assuming an air density of 
+$\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of $1.5~10^{-3}$ give the expression 
+used of $u_s$ as a function of the module of surface stress}. 
+For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as 
+at a finite depth $H_{LC}$ (which is often close to the mixed layer depth), and simply 
+varies as a sine function in between (a first-order profile for the Langmuir cell structures). 
+The resulting expression for $w_{LC}$ is :
+\begin{equation}
+w_{LC}  = \begin{cases}
+                   c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$} 	\\
+                   0                 				 &      \text{otherwise} 
+                 \end{cases}
+\end{equation}
+where $c_{LC} = 0.15$ has been chosen by \citep{Axell_JGR02} as a good compromise 
+to fit LES data. The chosen value yields maximum vertical velocities $w_{LC}$ of the order 
+of a few centimeters per second. The value of $c_{LC}$ is set through the \np{rn\_lc} 
+namelist parameter, having in mind that it should stay between 0.15 and 0.54 \citep{Axell_JGR02}. 
+
+The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations:
+$H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift
+can reach on its own by converting its kinetic energy to potential energy, according to 
+\begin{equation}
+- \int_{-H_{LC}}^0 { N^2\;z  \;dz} = \frac{1}{2} u_s^2
+\end{equation}
+
+
+%\subsubsection{Mixing just below the mixed layer}
+%---------------------------------------------------------------%
+
+% add here a description of "penetration of TKE" and the associated namelist parameters
+
+% -------------------------------------------------------------------------------------------------------------
+%        TKE discretization considerations
 % -------------------------------------------------------------------------------------------------------------
 \subsection{TKE discretization considerations (\key{zdftke})}
@@ -239,7 +338,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_TKE_time_scheme}  \begin{center}
+\begin{figure}[!t]   \begin{center}
 \includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_ZDF_TKE_time_scheme.pdf}
-\caption {Illustration of the TKE time integration and its links to the momentum and tracer time integration. }
+\caption{ \label{Fig_TKE_time_scheme} 
+Illustration of the TKE time integration and its links to the momentum and tracer time integration. }
 \end{center}  
 \end{figure}
@@ -389,6 +489,5 @@
 
 %--------------------------------------------------TABLE--------------------------------------------------
-\begin{table}[htbp]  \label{Tab_GLS}
-\begin{center}
+\begin{table}[htbp]	\begin{center}
 %\begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
 \begin{tabular}{ccccc}
@@ -408,7 +507,8 @@
 \hline
 \end{tabular}
-\caption {Set of predefined GLS parameters, or equivalently predefined turbulence models available with \key{zdfgls} and controlled by the \np{nn\_clos} namelist parameter.}
-\end{center}
-\end{table}
+\caption{   \label{Tab_GLS} 
+Set of predefined GLS parameters, or equivalently predefined turbulence models available 
+with \key{zdfgls} and controlled by the \np{nn\_clos} namelist parameter.}
+\end{center}	\end{table}
 %--------------------------------------------------------------------------------------------------------------
 
@@ -417,9 +517,12 @@
 value near physical boundaries (logarithmic boundary layer law). $C_{\mu}$ and $C_{\mu'}$ 
 are calculated from stability function proposed by \citet{Galperin_al_JAS88}, or by \citet{Kantha_Clayson_1994} 
-or one of the two functions suggested by \citet{Canuto_2001}  (\np{nn\_stab\_func} = 0, 1, 2 or 3, resp.}). The value of $C_{0\mu}$ depends of the choice of the stability function.
+or one of the two functions suggested by \citet{Canuto_2001}  (\np{nn\_stab\_func} = 0, 1, 2 or 3, resp.}). 
+The value of $C_{0\mu}$ depends of the choice of the stability function.
 
 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated 
 thanks to Dirichlet or Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp. 
-The wave effect on the mixing could be also being considered \citep{Craig_Banner_1994}.
+As for TKE closure , the wave effect on the mixing is considered when \np{ln\_crban}~=~true
+\citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}. The \np{rn\_crban} namelist parameter 
+is $\alpha_{CB}$ in \eqref{ZDF_Esbc} and \np{rn\_charn} provides the value of $\beta$ in \eqref{ZDF_Lsbc}. 
 
 The $\psi$ equation is known to fail in stably stratified flows, and for this reason 
@@ -433,4 +536,8 @@
 if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value.
 
+The time and space discretization of the GLS equations follows the same energetic 
+consideration as for the TKE case described in \S\ref{ZDF_tke_ene}  \citep{Burchard_OM02}. 
+Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}.
+
 % -------------------------------------------------------------------------------------------------------------
 %        K Profile Parametrisation (KPP) 
@@ -479,7 +586,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!htb] \label{Fig_npc}  	\begin{center}
+\begin{figure}[!htb]  	\begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_npc.pdf}
-\caption {Example of an unstable density profile treated by the non penetrative 
+\caption{  \label{Fig_npc} 
+Example of an unstable density profile treated by the non penetrative 
 convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from 
 the surface to the bottom. It is found to be unstable between levels 3 and 4. 
@@ -641,7 +749,8 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_zdfddm}  \begin{center}
+\begin{figure}[!t]   \begin{center}
 \includegraphics[width=0.99\textwidth]{./TexFiles/Figures/Fig_zdfddm.pdf}
-\caption {From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ 
+\caption{  \label{Fig_zdfddm}
+From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ 
 and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy 
 curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves 
@@ -986,9 +1095,9 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \label{Fig_ZDF_M2_K1_tmx}  \begin{center}
+\begin{figure}[!t] 	\begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_ZDF_M2_K1_tmx.pdf}
-\caption {(a) M2 and (b) K2 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$). }
-\end{center}  
-\end{figure}
+\caption{  \label{Fig_ZDF_M2_K1_tmx} 
+(a) M2 and (b) K2 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$). }
+\end{center}  	\end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Introduction.tex
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Introduction.tex	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Introduction.tex	(revision 2376)
@@ -103,15 +103,16 @@
 is a module related to the TRAcers equation, computing the Lateral DiFfussion. 
 The complete list of module names is presented in Appendix~\ref{Apdx_D}. 
-Furthermore, modules are organized in a few directories
- that correspond to their category, as indicated by the first three letters of their name. 
+Furthermore, modules are organized in a few directories that correspond to their category, 
+as indicated by the first three letters of their name (Tab.~\ref{Tab_chap}). 
 
 The manual mirrors the organization of the model. 
 After the presentation of the continuous equations (Chapter \ref{PE}), the following chapters 
-refer to specific terms of the equations each associated with a group of modules .
+refer to specific terms of the equations each associated with a group of modules (Tab.~\ref{Tab_chap}).
 
 
-\begin{table}[htbp] \label{tab1}
+%--------------------------------------------------TABLE--------------------------------------------------
+\begin{table}[!t]	
 %\begin{center} \begin{tabular}{|p{143pt}|l|l|}	\hline
-\begin{center} \begin{tabular}{|l|l|l|}	\hline
+\begin{center} 	\begin{tabular}{|l|l|l|}	\hline
 Chapter \ref{STP}	& -   		        & model time STePping environment \\		\hline
 Chapter \ref{DOM}	& DOM		& model DOMain \\		\hline
@@ -121,19 +122,24 @@
 Chapter \ref{LBC}	& LBC		& Lateral Boundary Conditions (also OBC and BDY)  \\ 		\hline
 Chapter \ref{LDF}	& LDF		& Lateral DiFfusion (parameterisations) \\	\hline
-Chapter \ref{ZDF}	& ZDF		& vertical (Z) DiFfusion  \\		\hline
+Chapter \ref{ZDF}	& ZDF		& vertical (Z) DiFfusion (parameterisations)  \\		\hline
 Chapter \ref{OBS}	& OBS		& OBServation and model comparison  \\		\hline
 Chapter \ref{ASM}	& ASM		& ASsimilation increment  \\		\hline
-Chapter \ref{MISC}	& ...		& Miscellaneous  topics (DIA, DTA, IOM, SOL, TRD, FLO...)	 \\			\hline
-\end{tabular} 	\end{center}
-\end{table}
+Chapter \ref{MISC}	& ...			& Miscellaneous  topics (DIA, DTA, IOM,   \\
+						&				   & SOL, TRD, FLO...)   \\       \hline
+Chapter \ref{CFG}	&  - 			& predefined configurations  \\		\hline
+\end{tabular} 	
+\caption{	 \label{Tab_chap}
+Organization of Chapters which miminc the one of the model directories. }
+\end{center}	\end{table}
+%--------------------------------------------------------------------------------------------------------------
 
- \vspace{1cm}   Nota Bene : \vspace{0.25cm}
 
 \subsubsection{Changes between releases}
 NEMO/OPA, like all research tools, is in perpetual evolution. The present document describes 
 the OPA version include in the release 3.3 of NEMO.  This release differs significantly
-from version 8, documented in \citet{Madec1998}.
+from version 8, documented in \citet{Madec1998}.\\
 
 $\bullet$ The main modifications from OPA v8 and NEMO/OPA v3.2 are :\\
+\\
 (1) transition to full native \textsc{Fortran} 90, deep code restructuring and drastic 
 reduction of CPP keys; \\
@@ -161,4 +167,5 @@
  \vspace{1cm}
 $\bullet$ The main modifications from NEMO/OPA v3.2 and  v3.2 are :\\
+\\
 (1) introduction of a modified leapfrog-Asselin filter time stepping scheme \citep{Leclair_Madec_OM09}; \\
 (2) additional scheme for  iso-neutral mixing \citep{Griffies_al_JPO98}, although it is still a "work in progress"; \\
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Namelist/namsbc
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Namelist/namsbc	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Namelist/namsbc	(revision 2376)
@@ -14,7 +14,4 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                         ("key_lim3" or "key_lim2)
-   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
-                           !                       =1 stress recomputed each ocean time step ("key_lim3" only)
-                           !                       =2 combination of 0 and 1 cases           ("key_lim3" only)
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
    ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf)
Index: branches/nemo_v3_3_beta/DOC/TexFiles/Namelist/namzdf_tke
===================================================================
--- branches/nemo_v3_3_beta/DOC/TexFiles/Namelist/namzdf_tke	(revision 2375)
+++ branches/nemo_v3_3_beta/DOC/TexFiles/Namelist/namzdf_tke	(revision 2376)
@@ -4,25 +4,24 @@
    rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
    rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
-   rn_ebb      =  60.      !  coef. of the surface input of tke
+   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
    rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
    rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
-   rn_bshear   =   1.e-20  !  background shear (>0)
    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                            !                 = 1 bounded by the local vertical scale factor
                            !                 = 2 first vertical derivative of mixing length bounded by 1
-                           !                 = 3 same criteria as case 2 but applied in a different way
+                           !                 = 3 as =2 with distinct disspipative an mixing length scale
    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
-   ln_mxl0     = .false.   !  mixing length scale surface value as function of wind stress (T) or not (F)
-   rn_lmin     =   0.001   !  interior buoyancy lenght scale minimum value
-   rn_lmin0    =   0.01    !  surface  buoyancy lenght scale minimum value
-   nn_etau     =   0       !  exponentially deceasing penetration of tke due to internal & intertial waves
-                           !        = 0 no penetration ( O(2 km) resolution)
-                           !        = 1 additional tke source
-                           !        = 2 additional tke source applied only at the base of the mixed layer
-   nn_htau     =   1       !  type of exponential decrease of tke penetration
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
+   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
+   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+   nn_etau     =   0       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+                           !        = 0 no penetration
+                           !        = 1 add a tke source below the ML
+                           !        = 2 add a tke source just at the base of the ML
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                            !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m at high latitudes
-   rn_efr      =   0.05    !  fraction of surface tke value which penetrates inside the ocean
-   ln_lc       = .false.   !  Langmuir cell parameterisation
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
