Index: trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 2531)
+++ trunk/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 2541)
@@ -193,4 +193,65 @@
 a more clever choice.
 
+
+% ================================================================
+% Accuracy and Reproducibility
+% ================================================================
+\section{Accuracy and Reproducibility (\mdl{lib\_fortran})}
+\label{MISC_fortran}
+
+\subsection{Issues with intrinsinc SIGN function (\key{nosignedzero})}
+\label{MISC_sign}
+
+The SIGN(A, B) is the \textsc {Fortran} intrinsic function delivers the magnitude 
+of A with the sign of B. For example, SIGN(-3.0,2.0) has the value 3.0.
+The problematic case is when the second argument is zero, because, on platforms 
+that support IEEE arithmetic, zero is actually a signed number. 
+There is a positive zero and a negative zero.
+
+In \textsc{Fortran}~90, the processor was required always to deliver a positive result for SIGN(A, B) 
+if B was zero. Nevertheless, in \textsc{Fortran}~95, the processor is allowed to do the correct thing 
+and deliver ABS(A) when B is a positive zero and -ABS(A) when B is a negative zero. 
+This change in the specification becomes apparent only when B is of type real, and is zero, 
+and the processor is capable of distinguishing between positive and negative zero, 
+and B is negative real zero. Then SIGN delivers a negative result where, under \textsc{Fortran}~90 
+rules,  it used to return a positive result. 
+This change may be especially sensitive for the ice model, so we overwrite the intrinsinc 
+function with our own function simply performing :   \\
+\verb?   IF( B >= 0.e0 ) THEN   ;   SIGN(A,B) = ABS(A)  ?    \\
+\verb?   ELSE                   ;   SIGN(A,B) =-ABS(A)     ?  \\
+\verb?   ENDIF    ? \\
+This feature can be found in \mdl{lib\_fortran} module and is effective when \key{nosignedzero}
+is defined. We use a CPP key as the overwritting of a intrinsic function can present 
+performance issues with some computers/compilers.
+
+
+\subsection{MPP reproducibility}
+\label{MISC_glosum}
+
+The numerical reproducibility of simulations on distributed memory parallel computers 
+is a critical issue. In particular, within NEMO global summation of distributed arrays 
+is most susceptible to rounding errors, and their propagation and accumulation cause 
+uncertainty in final simulation reproducibility on different numbers of processors.
+To avoid so, based on \citet{He_Ding_JSC01} review of different technics, 
+we use a so called self-compensated summation method. The idea is to estimate 
+the roundoff error, store it in a buffer, and then add it back in the next addition. 
+
+Suppose we need to calculate $b = a_1 + a_2 + a_3$. The following algorithm 
+will allow to split the sum in two ($sum_1 = a_{1} + a_{2}$ and $b = sum_2 = sum_1 + a_3$) 
+with exactly the same rounding errors as the sum performed all at once.
+\begin{align*}
+	sum_1 \ \  &= a_1 + a_2 \\
+	error_1     &= a_2 + ( a_1 - sum_1 ) \\
+	sum_2 \ \  &= sum_1 + a_3 + error_1 \\
+	error_2     &= a_3 + error_1 + ( sum_1 - sum_2 ) \\
+	b \qquad \ &= sum_2 \\
+\end{align*}
+This feature can be found in \mdl{lib\_fortran} module and is effective when \key{mpp\_rep}.
+In that case, all calls to glob\_sum function (summation over the entire basin excluding 
+duplicated rows and columns due to cyclic or north fold boundary condition as well as 
+overlap MPP areas). 
+Note this implementation may be sensitive to the optimization level. 
+
+
 % ================================================================
 % Model optimisation, Control Print and Benchmark
@@ -202,5 +263,5 @@
 %--------------------------------------------------------------------------------------------------------------
 
-% \gmcomment{why not make these bullets into subsections?}
+ \gmcomment{why not make these bullets into subsections?}
 
 
@@ -267,5 +328,5 @@
 It is a fast and rather easy method to use; which are attractive features for a large 
 number of ocean situations (variable bottom topography, complex coastal geometry, 
-variable grid spacing, islands, open or cyclic boundaries, etc ...). It does not require 
+variable grid spacing, open or cyclic boundaries, etc ...). It does not require 
 a search for an optimal parameter as in the SOR method. However, the SOR has 
 been retained because it is a linear solver, which is a very useful property when 
@@ -284,5 +345,5 @@
 \end{aligned}
 \end{equation}
-the five-point finite difference equation \eqref{Eq_psi_total} can be rewritten as:
+the resulting five-point finite difference equation is given by:
 \begin{equation}  \label{Eq_solmat}
 \begin{split}
@@ -480,889 +541,14 @@
 
 % ================================================================
-% Diagnostics
-% ================================================================
-\section{Diagnostics (DIA, IOM, TRD, FLO)}
-\label{MISC_diag}
-
-% -------------------------------------------------------------------------------------------------------------
-%       Standard Model Output 
-% -------------------------------------------------------------------------------------------------------------
-\subsection{Model Output (default or \key{iomput} or \key{dimgout} or \key{netcdf4})}
-\label{MISC_iom}
-
-%to be updated with Seb documentation on the IO
-
-The model outputs are of three types: the restart file, the output listing, 
-and the output file(s). The restart file is used internally by the code when 
-the user wants to start the model with initial conditions defined by a 
-previous simulation. It contains all the information that is necessary in 
-order for there to be no changes in the model results (even at the computer 
-precision) between a run performed with several restarts and the same run 
-performed in one step. It should be noted that this requires that the restart file 
-contain two consecutive time steps for all the prognostic variables, and 
-that it is saved in the same binary format as the one used by the computer 
-that is to read it (in particular, 32 bits binary IEEE format must not be used for 
-this file). The output listing and file(s) are predefined but should be checked 
-and eventually adapted to the user's needs. The output listing is stored in 
-the $ocean.output$ file. The information is printed from within the code on the 
-logical unit $numout$. To locate these prints, use the UNIX command 
-"\textit{grep -i numout}" in the source code directory.
-
-In the standard configuration, the user will find the model results in 
-NetCDF files containing mean values (or instantaneous values if 
-\key{diainstant} is defined) for every time-step where output is demanded.
-These outputs are defined in the \mdl{diawri} module. 
-When defining \key{dimgout}, the output are written in DIMG format,
-an IEEE output format.
-
-Since version 3.3, support for NetCDF4 chunking and (loss-less) compression has
-been included.  These options build on the standard NetCDF output and allow
-the user control over the size of the chunks via namelist settings. Chunking
-and compression can lead to significant reductions in file sizes for a small
-runtime overhead. For a fuller discussion on chunking and other performance
-issues the reader is referred to the NetCDF4 documentation found 
-\href{http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}.
-
-The new features are only available when the code has been linked with a
-NetCDF4 library (version 4.1 onwards, recommended) which has been built
-with HDF5 support (version 1.8.4 onwards, recommended). Datasets created
-with chunking and compression are not backwards compatible with NetCDF3
-"classic" format but most analysis codes can be relinked simply with the
-new libraries and will then read both NetCDF3 and NetCDF4 files. NEMO
-executables linked with NetCDF4 libraries can be made to produce NetCDF3
-files by setting the \np{ln\_nc4zip} logical to false in the \np{namnc4} 
-namelist:
-
-%------------------------------------------namnc4----------------------------------------------------
-\namdisplay{namnc4}
-%-------------------------------------------------------------------------------------------------------------
-
-If \key{netcdf4} has not been defined, these namelist parameters are not read. 
-In this case, \np{ln\_nc4zip} is set false and dummy routines for a few
-NetCDF4-specific functions are defined. These functions will not be used but
-need to be included so that compilation is possible with NetCDF3 libraries.
-
-When using NetCDF4 libraries, \key{netcdf4} should be defined even if the
-intention is to create only NetCDF3-compatible files. This is necessary to
-avoid duplication between the dummy routines and the actual routines present
-in the library. Most compilers will fail at compile time when faced with
-such duplication. Thus when linking with NetCDF4 libraries the user must
-define \key{netcdf4} and control the type of NetCDF file produced via the
-namelist parameter.
-
-Chunking and compression is applied only to 4D fields and there is no
-advantage in chunking across more than one time dimension since previously
-written chunks would have to be read back and decompressed before being
-added to. Therefore, user control over chunk sizes is provided only for the
-three space dimensions. The user sets an approximate number of chunks along
-each spatial axis. The actual size of the chunks will depend on global domain
-size for mono-processors or, more likely, the local processor domain size for
-distributed processing. The derived values are subject to practical minimum
-values (to avoid wastefully small chunk sizes) and cannot be greater than the
-domain size in any dimension. The algorithm used is:
-
-\begin{alltt}  {{\scriptsize 
-\begin{verbatim}
-     ichunksz(1) = MIN( idomain_size,MAX( (idomain_size-1)/nn_nchunks_i + 1 ,16 ) )
-     ichunksz(2) = MIN( jdomain_size,MAX( (jdomain_size-1)/nn_nchunks_j + 1 ,16 ) )
-     ichunksz(3) = MIN( kdomain_size,MAX( (kdomain_size-1)/nn_nchunks_k + 1 , 1 ) )
-     ichunksz(4) = 1
-\end{verbatim}
-}}\end{alltt} 
-
-\noindent As an example, setting:
-\vspace{-20pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-     nn_nchunks_i=4, nn_nchunks_j=4 and nn_nchunks_k=31
-\end{verbatim}
-}}\end{alltt} \vspace{-10pt}
-
-\noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1}
-respectively in the mono-processor case (i.e. global domain of {\small\tt 182x149x31}).
-An illustration of the potential space savings that NetCDF4 chunking and compression
-provides is given in table \ref{Tab_NC4} which compares the results of two short
-runs of the ORCA2\_LIM reference configuration with a 4x2 mpi partitioning. Note
-the variation in the compression ratio achieved which reflects chiefly the dry to wet 
-volume ratio of each processing region.
-
-%------------------------------------------TABLE----------------------------------------------------
-\begin{table} 	\begin{tabular}{lrrr}
-Filename & NetCDF3 & NetCDF4 & Reduction\\
-         &filesize & filesize & \% \\
-         &(KB)     & (KB)     & \\
-ORCA2\_restart\_0000.nc & 16420 & 8860 & 47\%\\
-ORCA2\_restart\_0001.nc & 16064 & 11456 & 29\%\\
-ORCA2\_restart\_0002.nc & 16064 & 9744 & 40\%\\
-ORCA2\_restart\_0003.nc & 16420 & 9404 & 43\%\\
-ORCA2\_restart\_0004.nc & 16200 & 5844 & 64\%\\
-ORCA2\_restart\_0005.nc & 15848 & 8172 & 49\%\\
-ORCA2\_restart\_0006.nc & 15848 & 8012 & 50\%\\
-ORCA2\_restart\_0007.nc & 16200 & 5148 & 69\%\\
-ORCA2\_2d\_grid\_T\_0000.nc & 2200 & 1504 & 32\%\\
-ORCA2\_2d\_grid\_T\_0001.nc & 2200 & 1748 & 21\%\\
-ORCA2\_2d\_grid\_T\_0002.nc & 2200 & 1592 & 28\%\\
-ORCA2\_2d\_grid\_T\_0003.nc & 2200 & 1540 & 30\%\\
-ORCA2\_2d\_grid\_T\_0004.nc & 2200 & 1204 & 46\%\\
-ORCA2\_2d\_grid\_T\_0005.nc & 2200 & 1444 & 35\%\\
-ORCA2\_2d\_grid\_T\_0006.nc & 2200 & 1428 & 36\%\\
-ORCA2\_2d\_grid\_T\_0007.nc & 2200 & 1148 & 48\%\\
-             ...            &  ... &  ... & ..  \\
-ORCA2\_2d\_grid\_W\_0000.nc & 4416 & 2240 & 50\%\\
-ORCA2\_2d\_grid\_W\_0001.nc & 4416 & 2924 & 34\%\\
-ORCA2\_2d\_grid\_W\_0002.nc & 4416 & 2512 & 44\%\\
-ORCA2\_2d\_grid\_W\_0003.nc & 4416 & 2368 & 47\%\\
-ORCA2\_2d\_grid\_W\_0004.nc & 4416 & 1432 & 68\%\\
-ORCA2\_2d\_grid\_W\_0005.nc & 4416 & 1972 & 56\%\\
-ORCA2\_2d\_grid\_W\_0006.nc & 4416 & 2028 & 55\%\\
-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}
-\end{table}
-%----------------------------------------------------------------------------------------------------
-
-Since version 3.2, an I/O server has been added which provides more
-flexibility in the choice of the fields to be output as well as how the 
-writing work is distributed over the processors in massively parallel
-computing. It is activated when \key{iomput} is defined. 
-
-When \key{iomput} is activated with \key{netcdf4} chunking and
-compression parameters for fields produced via \np{iom\_put} calls are
-set via an equivalent and identically named namelist to \np{namnc4} in
-\np{xmlio\_server.def}. Typically this namelist serves the mean files
-whilst the \np{ namnc4} in the main namelist file continues to serve the
-restart files. This duplication is unfortunate but appropriate since, if
-using io\_servers, the domain sizes of the individual files produced by the
-io\_server processes may be different to those produced by the invidual
-processing regions and different chunking choices may be desired.
-{ 
-
-% -------------------------------------------------------------------------------------------------------------
-%       Tracer/Dynamics Trends
-% -------------------------------------------------------------------------------------------------------------
-\subsection[Tracer/Dynamics Trends (TRD)]
-                  {Tracer/Dynamics Trends (\key{trdmld}, \key{trdtra}, \key{trddyn}, \key{trdmld\_trc})}
-\label{MISC_tratrd}
-
-%------------------------------------------namtrd----------------------------------------------------
-\namdisplay{namtrd} 
-%-------------------------------------------------------------------------------------------------------------
-
-When \key{trddyn} and/or \key{trddyn} CPP variables are defined, each 
-trend of the dynamics and/or temperature and salinity time evolution equations 
-is stored in three-dimensional arrays just after their computation ($i.e.$ at the end 
-of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). These trends are then 
-used in \mdl{trdmod} (see TRD directory) every \textit{nn\_trd } time-steps.
-
-What is done depends on the CPP keys defined:
-\begin{description}
-\item[\key{trddyn}, \key{trdtra}] : a check of the basin averaged properties of the momentum 
-and/or tracer equations is performed ; 
-\item[\key{trdvor}] : a vertical summation of the moment tendencies is performed, 
-then the curl is computed to obtain the barotropic vorticity tendencies which are output ;
-\item[\key{trdmld}] : output of the tracer tendencies averaged vertically  
-either over the mixed layer (\np{nn\_ctls}=0), 
-or       over a fixed number of model levels (\np{nn\_ctls}$>$1 provides the number of level), 
-or       over a spatially varying but temporally fixed number of levels (typically the base 
-of the winter mixed layer) read in \ifile{ctlsurf\_idx} (\np{nn\_ctls}=1). 
-\end{description}
-
-The units in the output file can be changed using the \np{nn\_ucf} namelist parameter. 
-For example, in case of salinity tendency the units are given by PSU/s/\np{nn\_ucf}.
-Setting \np{nn\_ucf}=86400 ($i.e.$ the number of second in a day) provides the tendencies in PSU/d.
-
-When \key{trdmld} is defined, two time averaging procedure are proposed.
-Setting \np{ln\_trdmld\_instant} to \textit{true}, a simple time averaging is performed, 
-so that the resulting tendency is the contribution to the change of a quantity between 
-the two instantaneous values taken at the extremities of the time averaging period.
-Setting \np{ln\_trdmld\_instant} to \textit{false}, a double time averaging is performed, 
-so that the resulting tendency is the contribution to the change of a quantity between 
-two \textit{time mean} values. The later option requires the use of an extra file, \ifile{restart\_mld}  
-(\np{ln\_trdmld\_restart}=true), to restart a run.
-
-
-Note that the mixed layer tendency diagnostic can also be used on biogeochemical models 
-via the \key{trdtrc} and \key{trdmld\_trc} CPP keys.
-
-% -------------------------------------------------------------------------------------------------------------
-%       On-line Floats trajectories
-% -------------------------------------------------------------------------------------------------------------
-\subsection{On-line Floats trajectories (FLO) (\key{floats})}
-\label{FLO}
-%--------------------------------------------namflo-------------------------------------------------------
-\namdisplay{namflo} 
-%--------------------------------------------------------------------------------------------------------------
-
-The on-line computation of floats advected either by the three dimensional velocity 
-field or constraint to remain at a given depth ($w = 0$ in the computation) have been 
-introduced in the system during the CLIPPER project. The algorithm used is based 
-either on the work of \cite{Blanke_Raynaud_JPO97} (default option), or on a $4^th$
-Runge-Hutta algorithm (\np{ln\_flork4}=true). Note that the \cite{Blanke_Raynaud_JPO97} 
-algorithm have the advantage of providing trajectories which are consistent with the 
-numeric of the code, so that the trajectories never intercept the bathymetry. 
-
-See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing 
-the off-line use of this marvellous diagnostic tool.
-
-% -------------------------------------------------------------------------------------------------------------
-%       Other Diagnostics
-% -------------------------------------------------------------------------------------------------------------
-\subsection{Other Diagnostics (\key{diahth}, \key{diaar5})}
-\label{MISC_diag_others}
-
-
-Aside from the standard model variables, other diagnostics can be computed 
-on-line. The available ready-to-add diagnostics routines can be found in directory DIA. 
-Among the available diagnostics the following ones are obtained when defining 
-the \key{diahth} CPP key: 
-
-- the mixed layer depth (based on a density criterion, \citet{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
-
-- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})
-
-- the depth of the 20\deg C isotherm (\mdl{diahth})
-
-- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth})
-
-The poleward heat and salt transports, their advective and diffusive component, and 
-the meriodional stream function can be computed on-line in \mdl{diaptr} by setting 
-\np{ln\_diaptr} to true (see the \textit{namptr} namelist below).  
-When \np{ln\_subbas}~=~true, transports and stream function are computed 
-for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S) 
-as well as for the World Ocean. The sub-basin decomposition requires an input file 
-(\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask 
-been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}). 
-
-%------------------------------------------namptr----------------------------------------------------
-\namdisplay{namptr} 
-%-------------------------------------------------------------------------------------------------------------
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] 	  \begin{center}
-\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf}
-\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). 
-Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay 
-are removed from the sub-basin. Note also that the Arctic Ocean has been split 
-into Atlantic and Pacific basins along the North fold line.  }
-\end{center}   \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-In addition, a series of diagnostics has been added in the \mdl{diaar5}. 
-They corresponds to outputs that are required for AR5 simulations 
-(see Section \ref{MISC_steric} below for one of them). 
-Activating those outputs requires to define the \key{diaar5} CPP key.
-
-
-
-% ================================================================
-% Steric effect in sea surface height
-% ================================================================
-\section{Steric effect in sea surface height}
-\label{MISC_steric}
-
-
-Changes in steric sea level are caused when changes in the density of the water 
-column imply an expansion or contraction of the column. It is essentially produced 
-through surface heating/cooling and to a lesser extent through non-linear effects of 
-the equation of state (cabbeling, thermobaricity...).
-Non-Boussinesq models contain all ocean effects within the ocean acting 
-on the sea level. In particular, they include the steric effect. In contrast, 
-Boussinesq models, such as \NEMO, conserve volume, rather than mass, 
-and so do not properly represent expansion or contraction. The steric effect is 
-therefore not explicitely represented.
-This approximation does not represent a serious error with respect to the flow field 
-calculated by the model \citep{Greatbatch_JGR94}, but extra attention is required
-when investigating sea level, as steric changes are an important 
-contribution to local changes in sea level on seasonal and climatic time scales.
-This is especially true for investigation into sea level rise due to global warming. 
-
-Fortunately, the steric contribution to the sea level consists of a spatially uniform 
-component that can be diagnosed by considering the mass budget of the world 
-ocean \citep{Greatbatch_JGR94}. 
-In order to better understand how global mean sea level evolves and thus how
-the steric sea level can be diagnosed, we compare, in the following, the 
-non-Boussinesq and Boussinesq cases.
-
-Let denote 
-$\mathcal{M}$ the total mass of liquid seawater ($\mathcal{M}=\int_D \rho dv$), 
-$\mathcal{V}$ the total volume of seawater ($\mathcal{V}=\int_D dv$), 
-$\mathcal{A}$ the total surface of the ocean ($\mathcal{A}=\int_S ds$), 
-$\bar{\rho}$ the global mean seawater (\textit{in situ}) density ($\bar{\rho}= 1/\mathcal{V} \int_D \rho \,dv$), and
-$\bar{\eta}$ the global mean sea level ($\bar{\eta}=1/\mathcal{A}\int_S \eta \,ds$).
-
-A non-Boussinesq fluid conserves mass. It satisfies the following relations:
-\begin{equation} \label{Eq_MV_nBq} 
-\begin{split} 
-\mathcal{M} &=  \mathcal{V}  \;\bar{\rho}      	\\
-\mathcal{V} &=  \mathcal{A}  \;\bar{\eta}	
-\end{split}
-\end{equation}
-Temporal changes in total mass is obtained from the density conservation equation :
-\begin{equation}  \label{Eq_Co_nBq}
-\frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} ) = \left. \frac{\textit{emp}}{e_3}\right|_\textit{surface}
-\end{equation}
-where $\rho$ is the \textit{in situ} density, and \textit{emp} the surface mass 
-exchanges with the other media of the Earth system (atmosphere, sea-ice, land). 
-Its global averaged leads to the total mass change 
-\begin{equation}  \label{Eq_Mass_nBq}
-\partial_t \mathcal{M} = \mathcal{A} \;\overline{\textit{emp}}
-\end{equation}
-where $\overline{\textit{emp}}=\int_S \textit{emp}\,ds$ is the net mass flux 
-through the ocean surface.
-Bringing \eqref{Eq_Mass_nBq} and the time derivative of \eqref{Eq_MV_nBq} 
-together leads to the evolution equation of the mean sea level
-\begin{equation} \label{Eq_ssh_nBq}
-  \partial_t \bar{\eta} =  \frac{\overline{\textit{emp}}}{ \bar{\rho}} 
-  					- \frac{\mathcal{V}}{\mathcal{A}}  \;\frac{\partial_t \bar{\rho} }{\bar{\rho}}
-\end{equation}
-The first term in equation \eqref{Eq_ssh_nBq} alters sea level by adding or 
-subtracting mass from the ocean. 
-The second term arises from temporal changes in the global mean 
-density; $i.e.$ from steric effects. 
-
-In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$ 
-appears multiplied by the gravity ($i.e.$ in the hydrostatic balance of the primitive Equations). 
-In particular, the mass conservation equation, \eqref{Eq_Co_nBq}, degenerates into 
-the incompressibility equation:
-\begin{equation}  \label{Eq_Co_Bq}
-\frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} ) =  \left. \frac{\textit{emp}}{\rho_o \,e_3}\right|_ \textit{surface}
-\end{equation}
-and the global average of this equation now gives the temporal change of the total volume,
-\begin{equation}  \label{Eq_V_Bq}
-  \partial_t \mathcal{V} =   \mathcal{A} \;\frac{\overline{\textit{emp}}}{\rho_o} 
-\end{equation}
-Only the volume is conserved, not mass, or, more precisely, the mass which is conserved is the 
-Boussinesq mass, $\mathcal{M}_o = \rho_o \mathcal{V}$. The total volume (or equivalently  
-the global mean sea level) is altered only by net volume fluxes across the ocean surface,  
-not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid.
- 
-Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be 
-diagnosed by considering the mass budget of the ocean. 
-The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface 
-mass flux must be compensated by a spatially uniform change in the mean sea level due to 
-expansion/contraction of the ocean \citep{Greatbatch_JGR94}. In others words, the Boussinesq 
-mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$, the  total mass of the ocean seen 
-by the Boussinesq model, via the steric contribution to the sea level, $\eta_s$, a spatially 
-uniform variable, as follows:
-\begin{equation}  \label{Eq_M_Bq}
-   \mathcal{M}_o  =  \mathcal{M} + \rho_o \,\eta_s \,\mathcal{A} 
-\end{equation}
-Any change in $\mathcal{M}$ which cannot be explained by the net mass flux through 
-the ocean surface is converted into a mean change in sea level. Introducing the total density 
-anomaly, $\mathcal{D}= \int_D d_a \,dv$, where $d_a= (\rho -\rho_o ) / \rho_o$  
-is the density anomaly used in \NEMO (cf. \S\ref{TRA_eos}) in \eqref{Eq_M_Bq}
-leads to a very simple form for the steric height:
-\begin{equation}  \label{Eq_steric_Bq}
-   \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D} 
-\end{equation}
-
-The above formulation of the steric height of a Boussinesq ocean requires four remarks.
-First, one can be tempted to define $\rho_o$ as the initial value of $\mathcal{M}/\mathcal{V}$,
-$i.e.$ set $\mathcal{D}_{t=0}=0$, so that the initial steric height is zero. We do not
-recommend that. Indeed, in this case $\rho_o$ depends on the initial state of the ocean. 
-Since $\rho_o$ has a direct effect on the dynamics of the ocean (it appears in the pressure 
-gradient term of the momentum equation) it is definitively not a good idea when 
-inter-comparing experiments. 
-We better recommend to fixe once for all $\rho_o$ to $1035\;Kg\,m^{-3}$. This value is a 
-sensible choice for the reference density used in a Boussinesq ocean climate model since, 
-with the exception of only a small percentage of the ocean, density in the World Ocean 
-varies by no more than 2$\%$ from this value (\cite{Gill1982}, page 47).
-
-Second, we have assumed here that the total ocean surface, $\mathcal{A}$, does not
-change when the sea level is changing as it is the case in all global ocean GCMs 
-(wetting and drying of grid point is not allowed). 
-  
-Third, the discretisation of \eqref{Eq_steric_Bq} depends on the type of free surface
-which is considered. In the non linear free surface case, $i.e.$ \key{vvl} defined, it is
-given by
-\begin{equation}  \label{Eq_discrete_steric_Bq}
-   \eta_s =  - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t} e_{2t} e_{3t} }
-   					{ \sum_{i,\,j,\,k}         e_{1t} e_{2t} e_{3t} } 
-\end{equation}
-whereas in the linear free surface, the volume above the \textit{z=0} surface must be explicitly taken 
-into account to better approximate the total ocean mass and thus the steric sea level:
-\begin{equation}  \label{Eq_discrete_steric_Bq}
-   \eta_s =  - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} d_a\; e_{1t}e_{2t} \eta }
-   						{\sum_{i,\,j,\,k} e_{1t}e_{2t}e_{3t} + \sum_{i,\,j}           e_{1t}e_{2t} \eta } 
-\end{equation}
-
-The fourth and last remark concerns the effective sea level and the presence of sea-ice.
-In the real ocean, sea ice (and snow above it)  depresses the liquid seawater through 
-its mass loading. This depression is a result of the mass of sea ice/snow system acting 
-on the liquid ocean. There is, however, no dynamical effect associated with these depressions 
-in the liquid ocean sea level, so that there are no associated ocean currents. Hence, the 
-dynamically relevant sea level is the effective sea level, $i.e.$ the sea level as if sea ice 
-(and snow) were converted to liquid seawater \citep{Campin_al_OM08}. However,
-in the current version of \NEMO the sea-ice is levitating above the ocean without 
-mass exchanges between ice and ocean. Therefore the model effective sea level
-is always given by $\eta + \eta_s$, whether or not there is sea ice present.
-
-In AR5 outputs, the thermosteric sea level is demanded. It is steric sea level due to 
-changes in ocean density arising just from changes in temperature. It is given by:
-\begin{equation}  \label{Eq_thermosteric_Bq}
-   \eta_s = - \frac{1}{\mathcal{A}} \int_D d_a(T,S_o,p_o) \,dv
-\end{equation}
-where $S_o$ and $p_o$ are the initial salinity and pressure, respectively.
-
-Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs
-the \key{diaar5} defined to be called.
-
-
-
-\gmcomment{                    % start of gmcomment
-
-
-% ================================================================
-% Diagnostics
-% ================================================================
-\section{Standard model Output (IOM)}
-\label{MISC_iom}
-
-% -------------------------------------------------------------------------------------------------------------
-%       Standard Model Output 
-% -------------------------------------------------------------------------------------------------------------
-%\subsection{Model Output (default or \key{iomput} }
-%\label{MISC_iom}
-
-
-
-\subsection{Basic knowledge}
-
-
-\subsubsection{ XML basic rules}
-
-XML tags begin with the less-than character ("$<$") and end with the greater-than character (''$>$''). 
-You use tags to mark the start and end of elements, which are the logical units of information 
-in an XML document. In addition to marking the beginning of an element, XML start tags also 
-provide a place to specify attributes. An attribute specifies a single property for an element, 
-using a name/value pair, for example: $<$a b="x" c="y" b="z"$>$ ... $<$/a$>$.
-See \href{http://www.xmlnews.org/docs/xml-basics.html}{here} for more details.
-
-\subsubsection{Structure of the xml file used in NEMO}
-
-The xml file is split into 3 parts:
-
-\textbf{field definition}: define all variables that can be output (all lines between   
-\texttt{$<$field\_definition$>$} and \texttt{$<$/field\_definition$>$})
-
-\textbf{file definition}: define the netcdf files to be created and the variables they will contain 
-(all lines between \texttt{ $<$file\_definition$>$} and \texttt{$<$/file\_definition$>$})
-
-\textbf{axis and grid definitions}: define the horizontal and vertical grids (all lines between 
-\texttt{$<$axis\_definition$>$} and \texttt{$<$/axis\_definition$>$} and all lines between 
-\texttt{$<$grid\_definition$>$} and \texttt{$<$/grid\_definition$>$})
-
-\subsubsection{Inheritance and group }
-
- Xml extensively uses the concept of inheritance. \\
-\\
-example 1: \\
-\vspace{-30pt}
-\begin{alltt}  {{\scriptsize    
-\begin{verbatim}
-   <field_definition operation="ave(X)" >
-      <field id="sst"                    />   <!-- averaged      sst --> 
-      <field id="sss" operation="inst(X)"/>   <!-- instantaneous sss --> 
-   </field_definition> 
-\end{verbatim}
-}}\end{alltt} 
-
-The field ''sst'' which is part (or a child) of the field\_definition will inherit the value ''ave(X)'' 
-of the attribute ''operation'' from its parent ''field definition''. Note that a child can overwrite 
-the attribute definition inherited from its parents. In the example above, the field ''sss'' will 
-therefore output instantaneous values instead of average values.
-
-example 2: Use (or overwrite) attributes value of a field when listing the variables included in a file
-\vspace{-20pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-   <field_definition>
-      <field id="sst" description="sea surface temperature" />   
-      <field id="sss" description="sea surface salinity"    />  
-   </field_definition>      
-
-   <file_definition>
-      <file id="file_1" />   
-            <field ref="sst"                              />  <!-- default def -->
-            <field ref="sss" description="my description" />  <!-- overwrite   -->
-      </file>   
-   </file_definition> 
-\end{verbatim}
-}}\end{alltt} 
-
-With the help of the inheritance, the concept of group allow to define a set of attributes 
-for several fields or files. 
-
-example 3, group of fields: define a group ''T\_grid\_variables'' identified with the name 
-''grid\_T''. By default variables of this group have no vertical axis but, following inheritance 
-rules, ''axis\_ref'' can be redefined for the field ''toce'' that is a 3D variable.
-\vspace{-30pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-   <field_definition>
-      <group id="grid_T" axis_ref="none" grid_ref="T_grid_variables">
-            <field id="sst"/>  
-            <field id="sss"/>  
-            <field id="toce" axis_ref="deptht"/>  <!-- overwrite axis def -->
-      </group>
-   </field_definition> 
-\end{verbatim}
-}}\end{alltt} 
-
-example 4, group of files: define a group of file with the attribute output\_freq equal to 432000 (5 days)
-\vspace{-30pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-   <file_definition>
-      <group id="5d" output_freq="432000">    <!-- 5d files -->
-         <file id="5d_grid_T" name="auto">   <!-- T grid file -->
-         ...
-         </file>
-         <file id="5d_grid_U" name="auto">   <!-- U grid file -->
-         ...
-         </file>
-      </group>
-   </file_definition> 
-\end{verbatim}
-}}\end{alltt} 
-
-\subsubsection{Control of the xml attributes from NEMO}
-
-The values of some attributes are automatically defined by NEMO (and any definition 
-given in the xml file is overwritten). By convention, these attributes are defined to ''auto'' 
-(for string) or ''0000'' (for integer) in the xml file (but this is not necessary). 
-
-Here is the list of these attributes: \\
-
-%table to be created here....
-
-tag ids affected by automatic definition of some of their attributes
-
-name attribute
-attribute value
-field\_definition
-freq\_op
-\np{rn\_rdt} (namelist)
-SBC
-freq\_op
-\np{rn\_rdt} $\times$ \np{nn\_fsbc} (namelist)
-1h, 2h, 3h, 4h, 6h, 12h
-1d, 3d, 5d
-1m, 2m, 3m, 4m, 6m
-1y, 2y, 5y, 10y
-\_grid\_T, \_grid\_U, \_grid\_V, \_grid\_W, \_icemod, \_ptrc\_T, \_diad\_T, \_scalar
-name
-filename defined by a call to rou{dia\_nam} following NEMO nomenclature
-EqT, EqU, EqW
-jbegin, ni, name\_suffix
-according to the grid
-TAO, RAMA and PIRATA moorings
-ibegin, jbegin, name\_suffix
-according to the grid
-
-
-\subsection{ Detailed functionalities }
-
-\subsubsection{Tag list}
-
-
-Table might be easier to read:        % table to create
-
-Tag
-Description
-Accepted attribute
-Accepted attribute value(s)
-Parent tag
-context
-define the model using the xml file
-id
-"nemo" or "n\_nemo" for the nth AGRIF zoom.
-simulation
-
-What do you think, Seb?
-
-
-\begin{description}
-
-\item[context]: define the model using the xml file. Id is the only attribute accepted. 
-Its value must be ''nemo'' or ''n\_nemo'' for the nth AGRIF zoom. Child of simulation tag.
-
-\item[field]: define the field to be output. Accepted attributes are axis\_ref, description, enable, 
-freq\_op, grid\_ref, id (if child of field\_definition), level, operation, name, ref (if child of file), 
-unit, zoom\_ref. Child of field\_definition, file or group of fields tag.
-
-\item[field\_definition]: definition of the part of the xml file corresponding to the field definition. 
-Accept the same attributes as field tag. Child of context tag.
-
-\item[group]: define a group of file or field. Accept the same attributes as file or field.
-
-\item[file]: define the output fileÕs characteristics. Accepted attributes are description, enable, 
-output\_freq, output\_level, id, name, name\_suffix. Child of file\_definition or group of files tag.
-
-\item[file\_definition]: definition of the part of the xml file corresponding to the file definition. 
-Accept the same attributes as file tag. Child of context tag.
-
-\item[axis]: definition of the vertical axis. Accepted attributes are description, id, positive, size, unit. 
-Child of axis\_definition tag.
-
-\item[axis\_definition]: definition of the part of the xml file corresponding to the vertical axis definition. 
-Accept the same attributes as axis tag. Child of context tag
-
-\item[grid]: definition of the horizontal grid. Accepted attributes are description and id. 
-Child of axis\_definition tag.
-
-\item[grid\_definition]: definition of the part of the xml file corresponding to the horizontal grid definition. 
-Accept the same attributes as grid tag. Child of context tag
-
-\item[zoom]: definition of a subdomain of an horizontal grid. Accepted attributes are description, id, 
-i/jbegin, ni/j. Child of grid tag.
-
-\end{description}
-
-
-\subsubsection{Attributes list}
-
-Applied to a tag or a group of tags.
-
-% table to be added ?
-Another table, perhaps?
-
-%%%%
-
-Attribute
-Applied to?
-Definition
-Comment
-axis\_ref
-field
-String defining the vertical axis of the variable. It refers to the id of the vertical axis defined in the axis tag.
-Use ''non'' if the variable has no vertical axis
-
-%%%%%%
-
-\begin{description}
-
-\item[axis\_ref]: field attribute. String defining the vertical axis of the variable. 
-It refers to the id of the vertical axis defined in the axis tag. 
-Use ''none'' if the variable has no vertical axis
-
-\item[description]: this attribute can be applied to all tags but it is used only with the field tag. 
-In this case, the value of description will be used to define, in the output netcdf file, 
-the attributes long\_name and standard\_name of the variable.
-
-\item[enabled]: field and file attribute. Logical to switch on/off the output of a field or a file. 
-
-\item[freq\_op]: field attribute (automatically defined, see part 1.4 (''control of the xml attributes'')). 
-An integer defining the frequency in seconds at which NEMO is calling iom\_put for this variable. 
-It corresponds to the model time step (rn\_rdt in the namelist) except for the variables computed 
-at the frequency of the surface boundary condition (rn\_rdt ? nn\_fsbc in the namelist).   
-
-\item[grid\_ref]: field attribute. String defining the horizontal grid of the variable. 
-It refers to the id of the grid tag.
-
-\item[ibegin]: zoom attribute. Integer defining the zoom starting point along x direction. 
-Automatically defined for TAO/RAMA/PIRATA moorings (see part 1.4).  
-
-\item[id]: exists for all tag. This is a string defining the name to a specific tag that will be used 
-later to refer to this tag. Tags of the same category must have different ids.
-
-\item[jbegin]: zoom attribute. Integer defining the zoom starting point along y direction. 
-Automatically defined for TAO/RAMA/PIRATA moorings and equatorial section (see part 1.4). 
-
-\item[level]: field attribute. Integer from 0 to 10 defining the output priority of a field. 
-See output\_level attribute definition
-
-\item[operation]: field attribute. String defining the type of temporal operation to perform on a variable. 
-Possible choices are ''ave(X)'' for temporal mean, ''inst(X)'' for instantaneous, ''t\_min(X)'' for temporal min 
-and ''t\_max(X)'' for temporal max.
-
-\item[output\_freq]: file attribute. Integer defining the operation frequency in seconds. 
-For example 86400 for daily mean.
-
-\item[output\_level]: file attribute. Integer from 0 to 10 defining the output priority of variables in a file: 
-all variables listed in the file with a level smaller or equal to output\_level will be output. 
-Other variables wonÕt be output even if they are listed in the file.  
-
-\item[positive]: axis attribute (always .FALSE.). Logical defining the vertical axis convention used 
-in \NEMO (positive downward). Define the attribute positive of the variable in the netcdf output file.
-
-\item[prec]: field attribute. Integer defining the output precision. 
-Not implemented, we always output real4 arrays. 
-
-\item[name]: field or file attribute. String defining the name of a variable or a file. 
-If the name of a file is undefined, its id is used as a name. 2 files must have different names. 
-Files with specific ids will have their name automatically defined (see part 1.4). 
-Note that is name will be automatically completed by the cpu number (if needed) and ''.nc''
-
-\item[name\_suffix]: file attribute. String defining a suffix to be inserted after the name 
-and before the cpu number and the ''.nc'' termination. Files with specific ids have an 
-automatic definition of their suffix (see part 1.4).
-
-\item[ni]: zoom attribute. Integer defining the zoom extent along x direction. 
-Automatically defined for equatorial sections (see part 1.4).  
-
-\item[nj]: zoom attribute. Integer defining the zoom extent along x direction. 
-
-\item[ref]: field attribute. String referring to the id of the field we want to add in a file.
-
-\item[size]: axis attribute. use unknown...
-
-\item[unit]: field attribute. String defining the unit of a variable and the associated 
-attribute in the netcdf output file.
-
-\item[zoom\_ref]: field attribute. String defining the subdomain of data on which 
-the file should be written (to ouput data only in a limited area). 
-It refers to the id of a zoom defined in the zoom tag.
-\end{description}
-
-
-\subsection{IO\_SERVER}
-
-\subsubsection{Attached or detached mode?}
-
-Iom\_put is based on the io\_server developed by Yann Meurdesoif from IPSL 
-(see \href{http://forge.ipsl.jussieu.fr/ioserver/browser}{here} for the source code or 
-see its copy in NEMOGCM/EXTERNAL directory). 
-This server can be used in ''attached mode'' (as a library) or in ''detached mode'' 
-(as an external executable on n cpus). In attached mode, each cpu of NEMO will output 
-its own subdomain. In detached mode, the io\_server will gather data from NEMO 
-and output them split over n files with n the number of cpu dedicated to the io\_server.
-
-\subsubsection{Control the io\_server: the namelist file xmlio\_server.def}
-
-%
-%Again, a small table might be more readable?
-%Name
-%Type
-%Description
-%Comment
-%Using_server
-%Logical 
-%Switch to use the server in attached or detached mode
-%(.TRUE. corresponding to detached mode).
-
-The control of the use of the io\_server is done through the namelist file of the io\_server 
-called xmlio\_server.def.
-
-\textbf{using\_server}: logical, switch to use the server in attached or detached mode 
-(.TRUE. corresponding to detached mode).
-
-\textbf{using\_oasis}: logical, set to .TRUE. if NEMO is used in coupled mode.
-
-\textbf{client\_id} = ''oceanx'' : character, used only in coupled mode. 
-Specify the id used in OASIS to refer to NEMO. The same id must be used to refer to NEMO 
-in the \$NBMODEL part of OASIS namcouple in the call of prim\_init\_comp\_proto in cpl\_oasis3f90
-
-\textbf{server\_id} = ''ionemo'' : character, used only in coupled mode. 
-Specify the id used in OASIS to refer to the IO\_SERVER when used in detached mode. 
-Use the same id to refer to the io\_server in the \$NBMODEL part of OASIS namcouple.
-
-\textbf{global\_mpi\_buffer\_size}: integer; define the size in Mb of the MPI buffer used by the io\_server.
-
-\subsubsection{Number of cpu used by the io\_server in detached mode}
-
-The number of cpu used by the io\_server is specified only when launching the model.
-Here is an example of 2 cpus for the io\_server and 6 cpu for opa using mpirun:
-
-\texttt{ -p 2 -e ./ioserver}
-
-\texttt{ -p 6 -e ./opa }
-
-
-\subsection{Practical issues}
-
-\subsubsection{Add your own outputs}
-
-It is very easy to add you own outputs with iom\_put. 4 points must be followed.
-\begin{description}
-\item[1-] in NEMO code, add a \\
-\texttt{      CALL iom\_put( 'identifier', array ) } \\
-where you want to output a 2D or 3D array.
-
-\item[2-] don't forget to add \\
-\texttt{   USE iom            ! I/O manager library }  \\
-in the list of used modules in the upper part of your module. 
-
-\item[3-] in the file\_definition part of the xml file, add the definition of your variable using the same identifier you used in the f90 code.
-\vspace{-20pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-   <field_definition>
-      ...
-      <field id="identifier" description="blabla" />   
-      ...
-   </field_definition> 
-\end{verbatim}
-}}\end{alltt} 
-attributes axis\_ref and grid\_ref must be consistent with the size of the array to pass to iom\_put.
-if your array is computed within the surface module each nn\_fsbc time\_step, 
-add the field definition within the group defined with the id ''SBC'': $<$group id=''SBC''...$>$
-
-\item[4-] add your field in one of the output files   \\
-\vspace{-20pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-   <file id="file_1" .../>   
-      ...
-      <field ref="identifier" />   
-      ...
-   </file>   
-\end{verbatim}
-}}\end{alltt} 
-
-\end{description}
-
-\subsubsection{Several time axes in the output file}
-
-If your output file contains variables with different operations (see operation definition), 
-IOIPSL will create one specific time axis for each operation. Note that inst(X) will have 
-a time axis corresponding to the end each output period whereas all other operators 
-will have a time axis centred in the middle of the output periods.
-
-\subsubsection{Error/bug messages from IOIPSL}
-
-If you get the following error in the standard output file:
-\vspace{-20pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-FATAL ERROR FROM ROUTINE flio_dom_set
- --> too many domains simultaneously defined 
- --> please unset useless domains 
- --> by calling flio_dom_unset 
-\end{verbatim}
-}}\end{alltt} 
-
-You must increase the value of dom\_max\_nb in fliocom.f90 (multiply it by 10 for example). 
-
-If you mix, in the same file, variables with different freq\_op (see definition above), 
-like for example variables from the surface module with other variables, 
-IOIPSL will print in the standard output file warning messages saying there may be a bug. 
-\vspace{-20pt}
-\begin{alltt}  {{\scriptsize
-\begin{verbatim}
-WARNING FROM ROUTINE histvar_seq   
- --> There were 10 errors in the learned sequence of variables  
- --> for file   4 
- --> This looks like a bug, please report it.
-\end{verbatim}
-}}\end{alltt} 
-
-Don't worry, there is no bug, everything is properly working!
-
-     }      %end  \gmcomment
-% ================================================================
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
