Index: /trunk/DOC/NEMO_book.tex
===================================================================
--- /trunk/DOC/NEMO_book.tex	(revision 2540)
+++ /trunk/DOC/NEMO_book.tex	(revision 2541)
@@ -207,9 +207,10 @@
  {\Huge NEMO ocean engine} 
 \rule{345pt}{1.5pt} \\ 
-{\small  -- version 3.3 --}   }
+ }
 %{ -- Draft --}   }
 %\date{\today}
 \date{
-October 2010  \\
+January 2011  \\
+{\small  -- version 3.3 --} \\
 ~  \\
 \textit{\small Note du P\^ole de mod\'{e}lisation de l'Institut Pierre-Simon Laplace No 27 }\\
@@ -241,5 +242,5 @@
 % ================================================================
 
-%\includeonly{./TexFiles/Chapters/Chap_MISC, ./TexFiles/Chapters/Chap_CFG}
+%\includeonly{./TexFiles/Chapters/Chap_MISC}
 %\includeonly{./TexFiles/Chapters/Chap_ZDF}
 %\includeonly{./TexFiles/Chapters/Chap_STP,./TexFiles/Chapters/Chap_SBC,./TexFiles/Chapters/Chap_TRA}
@@ -295,4 +296,6 @@
 
 \include{./TexFiles/Chapters/Chap_ZDF}			% Vertical diffusion
+
+\include{./TexFiles/Chapters/Chap_DIA}			% Miscellaneous topics
 
 \include{./TexFiles/Chapters/Chap_OBS}                  % Observation operator
@@ -317,4 +320,6 @@
 \include{./TexFiles/Chapters/Annex_ISO}	                  % Isoneutral diffusion using triads 
 %\include{./TexFiles/Chapters/Annex_E}	                  % Notes on some on going staff (no included in the DOC) 
+%\include{./TexFiles/Chapters/Annex_Fox-Kemper}   % Notes on Fox-Kemper (no included in the DOC) 
+%\include{./TexFiles/Chapters/Annex_EVP}        	% Notes on EVP (no included in the DOC) 
 
 % ================================================================
Index: /trunk/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- /trunk/DOC/TexFiles/Biblio/Biblio.bib	(revision 2540)
+++ /trunk/DOC/TexFiles/Biblio/Biblio.bib	(revision 2541)
@@ -55,4 +55,6 @@
 
 @STRING{JPO = {J. Phys. Oceanogr.}}
+
+@STRING{JSC = {J. Supercomputing.}}
 
 @STRING{JWS = {John Wiley and Sons}}
@@ -103,6 +105,9 @@
   title = {Re-scaled height coordinates for accurate representation of free-surface flows in ocean circulation models},
   journal = OM,
+  year = {2004},
   volume = {7},
-  year = {2004}
+  pages = {269--284},
+  doi = {10.1016/j.ocemod.2003.09.003},
+  url = {http://dx.doi.org/10.1016/j.ocemod.2003.09.003}
 }
 
@@ -445,4 +450,15 @@
 }
 
+@ARTICLE{Bouillon_al_OM09,
+  author = {S. Bouillon and M.A. Morales Maqueda and V. Legat and T. Fichefet},
+  title = {An Elastic-Viscous-Plastic Sea Ice Model formulated on Arakawa B and C Grids},
+  journal = OM,
+  year = {2009},
+  volume = {27},
+  pages = {174--184},
+  doi = {10.1016/j.ocemod.2009.01.004},
+  url = {http://dx.doi.org/10.1016/j.ocemod.2009.01.004}
+}
+
 @ARTICLE{Boulanger_al_GRL01,
   author = {J.-P. Boulanger and E. Durand and J.-P. Duvel and C. Menkes and P.
@@ -548,5 +564,5 @@
   author = {J.-M. Campin and J. Marshall and D. Ferreira},
   title = {Sea ice-ocean coupling using a rescaled vertical coordinate z*},
-  journal = {Ocean Modelling},
+  journal = OM,
   year = {2008},
   volume = {24},  number = {1-2},
@@ -588,5 +604,5 @@
   title = {North Atlantic Simulations with the Hybrid Coordinate Ocean Model
 	(HYCOM): Impact of the Vertical Coordinate Choice, Reference Pressure, and Thermobaricity},
-  journal = {Journal of Physical Oceanography},
+  journal = JPO,
   year = {2003},
   volume = {33},
@@ -664,5 +680,5 @@
   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},
+  journal = JPO,
   year = {1998},
   volume = {28},
@@ -684,5 +700,5 @@
  title = {A lateral boundary formulation for multi-level prediction models},
  year = {1976},
- journal = {QJRMS},
+ journal = QJRMS,
  volume = {102},
  pages = {405--418}
@@ -808,5 +824,5 @@
   title = {Helium Isotopic Constraints on Simulated Ocean Circulations - Implications
 	for abyssal theories},
-  journal = {Environmental Fluid Mechanics},
+  journal = EFM,
   year = {2009},
   volume = {10},  number = {1-2},
@@ -944,5 +960,5 @@
  year = {1994},
  title = {A storm surge prediction model for the northern Bay of Bengal with application to the cyclone disaster in April 1991},
- journal = {JPO},
+ journal = JPO,
  volume = {24},
  pages = {172--190}
@@ -1231,4 +1247,13 @@
 }
 
+@ARTICLE{He_Ding_JSC01,
+  author = {Y. He and C. H. Q. Ding},
+  title = {Using Accurate Arithmetics to Improve Numerical Reproducibility and Stability in Parallel Applications},
+  journal = JSC,
+  year = {2001},
+  volume = {18},
+  pages = {259--277},
+}
+
 @ARTICLE{Hirt_al_JCP74,
   author = {C. W. Hirt and A. A. Amsden and J. L. Cook},
@@ -1243,5 +1268,5 @@
   author = {R. Hofmeister and H. Burchard and J.-M. Beckers},
   title = {Non-uniform adaptive vertical grids for 3D numerical ocean models},
-  journal = {Ocean Modelling},
+  journal = OM,
   year = {2009},
   volume = {33},
@@ -1392,5 +1417,5 @@
   author = {L. Kantha and S. Carniel},
   title = {Comment on ''Generic length-scale equation for geophysical turbulence models'' by L. Umlauf and H. Burchard},
-  journal = {Journal of Marine Systems}, 
+  journal = JMS, 
   year = {2005},
   volume = {61},
@@ -2011,5 +2036,5 @@
 author = {G. Mellor and A. Blumberg},
 title = {Wave Breaking and Ocean Surface Layer Thermal Response},
-journal = {JPO},
+journal = JPO,
 volume = {34}, number = {3},
 pages = {693--698},
@@ -2307,5 +2332,5 @@
   title = {The regional oceanic modeling system (ROMS) - a split-explicit, free-surface,
 	topography-following-coordinate oceanic modelr},
-  journal = {Ocean Modelling},
+  journal = OM,
   year = {2005},
   volume = {9},  number = {4},
@@ -2549,5 +2574,5 @@
   author = {L. Umlauf and H. Burchard},
   title = {A generic length-scale equation for geophysical turbulence models},
-  journal = {JMS}, 
+  journal = JMS, 
   year = {2003},
   volume = {61},  number = {2},
@@ -2558,5 +2583,5 @@
   author = {L. Umlauf and H. Burchard},
   title = {Second-order turbulence closure models for geophysical boundary layers. A review of recent work},
-  journal = {Journal of Marine Systems}, 
+  journal = JMS, 
   year = {2005},
   volume = {25},
@@ -2591,5 +2616,5 @@
   title = {Simulating the mass balance and salinity of Arctic and Antarctic sea ice. 2. 
      Importance of sea ice salinity variations},
-  journal = {0M},
+  journal = OM,
   year = {2009},
   volume = {27},
Index: /trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex	(revision 2541)
@@ -12,13 +12,14 @@
 equation model adapted to regional and global ocean circulation problems. It is intended to 
 be a flexible tool for studying the ocean and its interactions with the others components of 
-the earth climate system (atmosphere, sea-ice, biogeochemical tracers, ...) over a wide range 
-of space and time scales. Prognostic variables are the three-dimensional velocity field, a linear 
+the earth climate system over a wide range of space and time scales. 
+Prognostic variables are the three-dimensional velocity field, a linear 
 or non-linear sea surface height, the temperature and the salinity. In the horizontal direction, 
 the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step 
 $z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a 
 three-dimensional Arakawa C-type grid. Various physical choices are available to describe 
-ocean physics, including TKE and KPP vertical physics. Within NEMO, the ocean is interfaced 
-with a sea-ice model (LIM v2 and v3), passive tracer and biogeochemical models (TOP) 
-and, via the OASIS coupler, with several atmospheric general circulation models. 
+ocean physics, including TKE, GLS and KPP vertical physics. Within NEMO, the ocean is 
+interfaced with a sea-ice model (LIM v2 and v3), passive tracer and biogeochemical models (TOP) 
+and, via the OASIS coupler, with several atmospheric general circulation models. It also 
+support two-way grid embedding via the AGRIF software.
 
 % ================================================================
@@ -28,14 +29,15 @@
 mod\`{e}le aux \'{e}quations primitives de la circulation oc\'{e}anique r\'{e}gionale et globale. 
 Il se veut un outil flexible pour \'{e}tudier sur un vaste spectre spatiotemporel l'oc\'{e}an et ses 
-interactions avec les autres composantes du syst\`{e}me climatique terrestre (atmosph\`{e}re, 
-glace de mer, traceurs biog\'{e}ochimiques...). Les variables pronostiques sont le champ 
-tridimensionnel de vitesse, une hauteur de la mer lin\'{e}aire ou non, la temperature et la salinit\'{e}. 
+interactions avec les autres composantes du syst\`{e}me climatique terrestre. 
+Les variables pronostiques sont le champ tridimensionnel de vitesse, une hauteur de la mer 
+lin\'{e}aire ou non, la temperature et la salinit\'{e}. 
 La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une 
 coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore 
 une combinaison des deux. Diff\'{e}rents choix sont propos\'{e}s pour d\'{e}crire la physique 
-oc\'{e}anique, incluant notamment des physiques verticales TKE et KPP. A travers l'infrastructure 
-NEMO, l'oc\'{e}an est interfac\'{e} avec un mod\`{e}le de glace de mer, des mod\`{e}les 
-biog\'{e}ochimiques et de traceur passif, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les 
-de circulation g\'{e}n\'{e}rale atmosph\'{e}rique. 
+oc\'{e}anique, incluant notamment des physiques verticales TKE, GLS et KPP. A travers l'infrastructure 
+NEMO, l'oc\'{e}an est interfac\'{e} avec des mod\`{e}les de glace de mer, de biog\'{e}ochimie 
+et de traceurs passifs, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les de circulation 
+g\'{e}n\'{e}rale atmosph\'{e}rique. Il supporte \'{e}galement l'embo\^{i}tement interactif de 
+maillages via le logiciel AGRIF.
 } 
 
Index: /trunk/DOC/TexFiles/Chapters/Annex_D.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Annex_D.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Annex_D.tex	(revision 2541)
@@ -88,6 +88,6 @@
 
 - In the declaration of a PUBLIC variable, the comment part at the end of the line 
-should start with the two characters "\verb?!:?". the following UNIX command, 
-\verb?grep var_name *90 \ grep \!: ? 
+should start with the two characters "\verb?!:?". the following UNIX command, \\
+\verb?grep var_name *90 \ grep \!: ? \\
 will display the module name and the line where the var\_name declaration is.
 
Index: /trunk/DOC/TexFiles/Chapters/Chap_CFG.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_CFG.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Chap_CFG.tex	(revision 2541)
@@ -3,5 +3,5 @@
 % ================================================================
 \chapter{Configurations}
-\label{MISC}
+\label{CFG}
 \minitoc
 
@@ -69,5 +69,5 @@
 % ORCA family configurations
 % ================================================================
-\section{ORCA family: global ocean with tripolar grid}
+\section{ORCA family: global ocean with tripolar grid (\key{orca\_rX})}
 \label{CFG_orca}
 
@@ -140,6 +140,7 @@
 The NEMO system is provided with five built-in ORCA configurations which differ in the 
 horizontal resolution. The value of the resolution is given by the resolution at the Equator 
-expressed in degrees. Each of configuration is set through a CPP key with set the grid size 
-and configuration name parameters  (Tab.~\ref{Tab_ORCA}).
+expressed in degrees. Each of configuration is set through a CPP key, \key{orca\_rX} 
+(with X being an indicator of the resolution), which set the grid size and configuration 
+name parameters  (Tab.~\ref{Tab_ORCA}).
 .
 
Index: /trunk/DOC/TexFiles/Chapters/Chap_DIA.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_DIA.tex	(revision 2541)
+++ /trunk/DOC/TexFiles/Chapters/Chap_DIA.tex	(revision 2541)
@@ -0,0 +1,903 @@
+% ================================================================
+% Chapter Ñ I/O & Diagnostics
+% ================================================================
+\chapter{Ouput and Diagnostics (IOM, DIA, TRD, FLO)}
+\label{DIA}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new ligne
+
+% ================================================================
+%       Old Model Output 
+% ================================================================
+\section{Old Model Output (default or \key{dimgout})}
+\label{DIA_io_old}
+
+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.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 presented in next section. 
+
+
+%\gmcomment{                    % start of gmcomment
+
+% ================================================================
+% Diagnostics
+% ================================================================
+\section{Standard model Output (IOM)}
+\label{DIA_iom}
+
+
+Since version 3.2, iom\_put is the NEMO output interface. It was designed to be simple to use, 
+flexible and efficient. Two main functionalities are covered by iom\_put: 
+(1) the control of the output files through an external xml file defined by the user ; 
+(2) the distribution (or not) of all task related to output files on dedicated processors.
+The first functionality allows the user to specify, without touching anything into the code, 
+the way he want to output data: \\
+- choice of output frequencies that can be different for each file (including real months and years) \\
+- choice of file contents: decide which data will be written in which file (the same data can be 
+outputted in different files)  \\
+- possibility to extract a subdomain (for example all TAO-PIRATA-RAMA moorings are already defined)  \\
+- choice of the temporal operation to perform: mean, instantaneous, min, max  \\
+- extremely large choice of data available   \\
+- redefine variables name and long\_name  \\
+In addition, iom\_put allows the user to output any variable (scalar, 2D or 3D) in the code 
+in a very easy way. All details of iom\_put functionalities are listed in the following subsections.
+An example of the iodef.xml file that control the outputs can be found here: 
+NEMOGCM/CONFIG/ORCA2\_LIM/EXP00/iodef.xml
+
+The second functionality targets outputs performances when running on a very large number of processes. 
+The idea is to dedicate N specific processes to write the outputs, where N is defined by the user. 
+In the current version, this functionality is technically working however, its performance are usually poor 
+(for known reasons). Users can therefore test this functionality but they must be aware that expected 
+performance improvement will not be achieved before the release 3.4.
+An example of xmlio\_server.def NEMOGCM/CONFIG/ORCA2\_LIM/EXP00/xmlio\_server.def
+ 
+
+\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:
+\begin{description}
+\item[field definition]: define all variables that can be output \\
+all lines in between the following two tags\\
+\verb?   <field\_definition ...> ?  \\
+\verb?   </field\_definition ...> ? 
+\item[file definition]: define the netcdf files to be created and the variables they will contain \\
+all lines in between the following two tags \\
+\verb?   <field\_definition> ?  \\
+\verb?   </field\_definition> ?
+\item[axis and grid definitions]: define the horizontal and vertical grids \\
+all lines in between the following two set of two tags\\
+\verb?   <axis\_definition ...> ?  \\
+\verb?   </axis\_definition ...> ?
+and \\
+\verb?   <grid\_definition ...> ?  \\
+\verb?   </grid\_definition ...> ?
+\end{description}
+
+\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: \\
+\\
+\begin{tabular}{|l|c|c|c|}
+   \hline
+ \multicolumn{2}{|c|}{tag ids affected by automatic           }  & name      & attribute value \\
+  \multicolumn{2}{|c|}{definition of some of their attributes }  & attribute  &        \\
+   \hline
+   \hline
+    \multicolumn{2}{|c|}{field\_definition} & freq\_op & \np{rn\_rdt} \\
+   \hline
+    \multicolumn{2}{|c|}{SBC}                  & freq\_op & \np{rn\_rdt} $\times$ \np{nn\_fsbc} \\
+   \hline
+    1h, 2h, 3h, 4h, 6h, 12h & \_grid\_T, \_grid\_U,  &  name & filename defined by   \\
+    1d, 3d, 5d                     & \_grid\_V, \_grid\_W, &            & a call to rou{dia\_nam}  \\
+    1m, 2m, 3m, 4m, 6m    & \_icemod, \_ptrc\_T,  &            & following NEMO \\
+    1y, 2y, 5y, 10y               &  \_diad\_T, \_scalar   &            & nomenclature \\
+   \hline
+    \multicolumn{2}{|c|}{EqT, EqU, EqW} & jbegin, ni,      & according to the grid    \\
+    \multicolumn{2}{|c|}{                         } & name\_suffix &                                      \\
+   \hline
+    \multicolumn{2}{|c|}{TAO, RAMA and PIRATA moorings} & ibegin, jbegin,      & according to the grid    \\
+    \multicolumn{2}{|c|}{                                                       } & name\_suffix &                                      \\
+   \hline
+\end{tabular}
+
+
+\subsection{ Detailed functionalities }
+
+\subsubsection{Tag list}
+
+\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
+
+
+% ================================================================
+%       NetCDF4 support
+% ================================================================
+\section{NetCDF4 Support (\key{netcdf4})}
+\label{DIA_iom}
+
+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 \textit{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}
+%----------------------------------------------------------------------------------------------------
+
+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 \textit{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
+% -------------------------------------------------------------------------------------------------------------
+\section[Tracer/Dynamics Trends (TRD)]
+                  {Tracer/Dynamics Trends  (\key{trdtra}, \key{trddyn},    \\ 
+                                                             \key{trddvor}, \key{trdmld})}
+\label{DIA_trd}
+
+%------------------------------------------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
+% -------------------------------------------------------------------------------------------------------------
+\section{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
+% -------------------------------------------------------------------------------------------------------------
+\section{Other Diagnostics (\key{diahth}, \key{diaar5})}
+\label{DIA_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-basins. 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{Diagnosing the Steric effect in sea surface height}
+\label{DIA_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.
+
+% ================================================================
+
+
+
+
+
+
+
+
+
+
+
+
Index: /trunk/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2541)
@@ -900,8 +900,19 @@
 
 The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation. 
-The extra term introduced in the equations (see {\S}I.2.2) is solved implicitly. 
+The extra term introduced in the equations (see \S\ref{PE_free_surface}) is solved implicitly. 
 The elliptic solvers available in the code are documented in \S\ref{MISC}.
 
 %% gm %%======>>>>   given here the discrete eqs provided to the solver
+\gmcomment{               %%% copy from chap-model basics 
+\begin{equation} \label{Eq_spg_flt}
+\frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}
+- g \nabla \left( \tilde{\rho} \ \eta \right) 
+- g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) 
+\end{equation}
+where $T_c$, is a parameter with dimensions of time which characterizes the force, 
+$\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ 
+represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 
+non-linear and viscous terms in \eqref{Eq_PE_dyn}.
+}   %end gmcomment
 
 Note that in the linear free surface formulation (\key{vvl} not defined), the ocean depth 
Index: /trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 2540)
+++ /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
-% ================================================================
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
Index: /trunk/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2541)
@@ -63,5 +63,4 @@
 \label{SBC_general}
 
-
 The surface ocean stress is the stress exerted by the wind and the sea-ice 
 on the ocean. The two components of stress are assumed to be interpolated 
@@ -176,4 +175,260 @@
 
 % ================================================================
+%       Input Data 
+% ================================================================
+\section{Input Data generic interface}
+\label{SBC_input}
+
+A generic interface has been introduced to manage the way input data (2D or 3D fields, 
+like surface forcing or ocean T and S) are specify in \NEMO. This task is archieved by fldread.F90. 
+The module was design with four main objectives in mind: 
+\begin{enumerate}  
+\item optionally provide a time interpolation of the input data at model time-step, 
+whatever their input frequency is, and according to the different calendars available in the model.
+\item optionally provide an on-the-fly space interpolation from the native input data grid to the model grid.
+\item make the run duration independent from the period cover by the input files.
+\item provide a simple user interface and a rather simple developer interface by limiting the 
+ number of prerequisite information. 
+\end{enumerate}  
+
+As a results the user have only to fill in for each variable a structure in the namelist file 
+to defined the input data file and variable names, the frequency of the data (in hours or months), 
+whether its is climatological data or not, the period covered by the input file (one year, month, week or day), 
+and two additional parameters for on-the-fly interpolation. When adding a new input variable, 
+the developer has to add the associated structure in the namelist, read this information 
+by mirroring the namelist read in \rou{sbc\_blk\_init} for example, and simply call \rou{fld\_read} 
+to obtain the desired input field at the model time-step and grid points.
+
+The only constraints are that the input file is a NetCDF file, the file name follows a nomenclature 
+(see \S\ref{SBC_fldread}), the period it cover is one year, month, week or day, and, if on-the-fly 
+interpolation is used, a file of weights must be supplied (see \S\ref{SBC_iof}).
+
+Note that when an input data is archived on a disc which is accessible directly 
+from the workspace where the code is executed, then the use can set the \np{cn\_dir} 
+to the pathway leading to the data. By default, the data are assumed to have been 
+copied so that cn\_dir='./'.
+
+% -------------------------------------------------------------------------------------------------------------
+% Input Data specification (\mdl{fldread})
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Input Data specification (\mdl{fldread})}
+\label{SBC_fldread}
+
+The structure associated with an input variable contains the following information:
+\begin{alltt}  {{\tiny    
+\begin{verbatim}
+!  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! 
+!             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+\end{verbatim}
+}}\end{alltt} 
+where 
+\begin{description}  
+\item[File name]: the stem name of the NetCDF file to be open. 
+This stem will be completed automatically by the model, with the addition of a '.nc' at its end 
+and by date information and possibly a prefix (when using AGRIF). 
+Tab.\ref{Tab_fldread} provides the resulting file name in all possible cases according to whether 
+it is a climatological file or not, and to the open/close frequency (see below for definition). 
+
+%--------------------------------------------------TABLE--------------------------------------------------
+\begin{table}[htbp] 
+\begin{center}
+\begin{tabular}{|l|c|c|c|}
+\hline
+                         & daily or weekLLL	        & monthly                   &   yearly          \\   \hline
+clim = false	& fn\_yYYYYmMMdDD  &   fn\_yYYYYmMM   &   fn\_yYYYY  \\   \hline
+clim = true	 	   & not possible 	              &  fn\_m??.nc             &   fn                \\   \hline
+\end{tabular}
+\end{center}
+\caption{ \label{Tab_fldread}   naming nomenclature for climatological or interannual input file, 
+as a function of the Open/close frequency. The stem name is assumed to be 'fn'. 
+For weekly files, the 'LLL' corresponds to the first three letters of the first day of the week ($i.e.$ 'sun','sat','fri','thu','wed','tue','mon'). The 'YYYY', 'MM' and 'DD' should be replaced by the 
+actual year/month/day, always coded with 4 or 2 digits. Note that (1) in mpp, if the file is split 
+over each subdomain, the suffix '.nc' is replaced by '\_PPPP.nc', where 'PPPP' is the 
+process number coded with 4 digits; (2) when using AGRIF, the prefix ÔN\_Õ is added to files, 
+where 'N'  is the child grid number.}
+\end{table}
+%--------------------------------------------------------------------------------------------------------------
+  
+
+\item[Record frequency]: the frequency of the records contained in the input file. 
+Its unit is in hours if it is positive (for example 24 for daily forcing) or in months if negative 
+(for example -1 for monthly forcing or -12 for annual forcing). 
+Note that this frequency must really be an integer and not a real. 
+On some computers, seting it to '24.' can be interpreted as 240!
+
+\item[Variable name]: the name of the variable to be read in the input NetCDF file.
+
+\item[Time interpolation]: a logical to activate, or not, the time interpolation. If set to 'false', 
+the forcing will have a steplike shape remaining constant during each forcing period. 
+For example, when using a daily forcing without time interpolation, the forcing remaining 
+constant from 00h00'00'' to 23h59'59". If set to 'true', the forcing will have a broken line shape. 
+Records are assumed to be dated the middle of the forcing period. 
+For example, when using a daily forcing with time interpolation, linear interpolation will 
+be performed between mid-day of two consecutive days. 
+
+\item[Climatological forcing]: a logical to specify if a input file contains climatological forcing 
+which can be cycle in time, or an interannual forcing which will requires additional files 
+if the period covered by the simulation exceed the one of the file. See the above the file 
+naming strategy which impacts the expected name of the file to be opened. 
+
+\item[Open/close frequency]: the frequency at which forcing files must be opened/closed. 
+Four cases are coded: 'daily', 'weekLLL' (with 'LLL' the first 3 letters of the first day of the week), 
+'monthly' and 'yearly' which means the forcing files will contain data for one day, one week, 
+one month or one year. Files are assumed to contain data from the beginning of the open/close period. 
+For example, the first record of a yearly file containing daily data is Jan 1st even if the experiment 
+is not starting at the beginning of the year. 
+
+\item[Others]: 'weights filename' and 'pairing rotation' are associted with on-the-fly interpolation 
+which is described in \S\ref{SBC_iof}.
+
+\end{description}
+
+Additional remarks:\\
+(1) The time interpolation is a simple linear interpolation between two consecutive records of 
+the input data. The only tricky point is therefore to specify the date at which we need to do 
+the interpolation and the date of the records read in the input files. 
+Following \citet{Leclair_Madec_OM09}, the date of a time step is set at the middle of the 
+time step. For example, for an experiment starting at 0h00'00" with a one hour time-step, 
+a time interpolation will be performed at the following time: 0h30'00", 1h30'00", 2h30'00", etc.
+However, for forcing data related to the surface module, values are not needed at every 
+time-step but at every \np{nn\_fsbc} time-step. For example with \np{nn\_fsbc}~=~3, 
+the surface module will be called at time-steps 1, 4, 7, etc. The date used for the time interpolation 
+is thus redefined to be at the middle of \np{nn\_fsbc} time-step period. In the previous example, 
+this leads to: 1h30'00", 4h30'00", 7h30'00", etc. \\ 
+(2) For code readablility and maintenance issues, we don't take into account the NetCDF input file 
+calendar. The calendar associated with the forcing field is build according to the information 
+provided by user in the record frequency, the open/close frequency and the type of temporal interpolation. 
+For example, the first record of a yearly file containing daily data that will be interpolated in time 
+is assumed to be start Jan 1st at 12h00'00" and end Dec 31st at 12h00'00". \\
+(3) If a time interpolation is requested, the code will pick up the needed data in the previous (next) file 
+when interpolating data with the first (last) record of the open/close period. 
+For example, if the input file specifications are ''yearly, containing daily data to be interpolated in time'', 
+the values given by the code between 00h00'00" and 11h59'59" on Jan 1st will be interpolated values 
+between Dec 31st 12h00'00" and Jan 1st 12h00'00". If the forcing is climatological, Dec and Jan will 
+be keep-up from the same year. However, if the forcing is not climatological, at the end of the 
+open/close period the code will automatically close the current file and open the next one. 
+Note that, if the experiment is starting (ending) at the beginning (end) of an open/close period 
+we do accept that the previous (next) file is not existing. In this case, the time interpolation 
+will be performed between two identical values. For example, when starting an experiment on 
+Jan 1st of year Y with yearly files and daily data to be interpolated, we do accept that the file 
+related to year Y-1 is not existing. The value of Jan 1st will be used as the missing one for 
+Dec 31st of year Y-1. If the file of year Y-1 exists, the code will read its last record. 
+Therefore, this file can contain only one record corresponding to Dec 31st, a useful feature for 
+user considering that it is too heavy to manipulate the complete file for year Y-1.
+
+
+% -------------------------------------------------------------------------------------------------------------
+% Interpolation on the Fly
+% -------------------------------------------------------------------------------------------------------------
+\subsection [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 
+\href{http://climate.lanl.gov/Software/SCRIP}{here} under a copyright agreement).
+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.
+
+\subsubsection{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.
+
+\subsubsection{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.
+
+\subsubsection{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).
+
+\subsubsection{Limitations}
+\label{SBC_iof_lim}
+
+\begin{enumerate}  
+\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
+          when a bicubic interpolation method is used.
+\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{enumerate}
+
+\subsubsection{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).
+
+
+% ================================================================
 % Analytical formulation (sbcana module) 
 % ================================================================
@@ -215,19 +470,5 @@
 read in the file, the time frequency at which it is given (in hours), and a logical 
 setting whether a time interpolation to the model time step is required 
-for this field). (fld\_i namelist structure).
-
-\textbf{Caution}: when the frequency is set to --12, the data are monthly 
-values. These are assumed to be climatological values, so time interpolation 
-between December the 15$^{th}$ and January the 15$^{th}$ is done using 
-records 12 and 1
-
-When higher frequency is set and time interpolation is demanded, the model 
-will try to read the last (first) record of previous (next) year in a file 
-having the same name but a suffix {\_}prev{\_}year ({\_}next{\_}year) being 
-added (e.g. "{\_}1989"). These files must only contain a single record. If they don't exist, 
-the model assumes that the last record of the previous year is equal to the first 
-record of the current year, and similarly, that the first record of the 
-next year is equal to the last record of the current year. This will cause 
-the forcing to remain constant over the first and last half fld\_frequ hours.
+for this field. See \S\ref{SBC_fldread} for a more detailed description of the parameters.
 
 Note that in general, a flux formulation is used in associated with a 
@@ -281,5 +522,5 @@
 \begin{table}[htbp]   \label{Tab_CORE}
 \begin{center}
-\begin{tabular}{|l|l|l|l|}
+\begin{tabular}{|l|c|c|c|}
 \hline
 Variable desciption					& Model variable	& Units	 & point \\		\hline
@@ -297,6 +538,8 @@
 %--------------------------------------------------------------------------------------------------------------
 
-Note that the air velocity is provided at a tracer ocean point, not at a velocity ocean point ($u$- and $v$-points). It is simpler and faster (less fields to be read), but it is not the recommended method when the ocean grid
-size is the same or larger than the one of the input atmospheric fields.
+Note that the air velocity is provided at a tracer ocean point, not at a velocity ocean 
+point ($u$- and $v$-points). It is simpler and faster (less fields to be read), 
+but it is not the recommended method when the ocean grid size is the same 
+or larger than the one of the input atmospheric fields.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -338,6 +581,5 @@
 As for the flux formulation, information about the input data required by the 
 model is provided in the namsbc\_blk\_core or namsbc\_blk\_clio 
-namelist (via the structure fld\_i). The first and last record assumption is also made 
-(see \S\ref{SBC_flx})
+namelist (see \S\ref{SBC_fldread}). 
 
 % ================================================================
@@ -399,5 +641,5 @@
 (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).
-$\eta_{ib}$ can be set in the output. This can simplify the altirmetry data and model comparison 
+$\eta_{ib}$ can be set in the output. This can simplify altimetry data and model comparison 
 as inverse barometer sea surface height is usually removed from these date prior to their distribution.
 
@@ -433,5 +675,6 @@
 River runoff generally enters the ocean at a nonzero depth rather than through the surface.
 Many models, however, have traditionally inserted river runoff to the top model cell.
-This was the case in \NEMO prior to the version 3.3, and was combined with an option to increase vertical mixing near the river mouth.
+This was the case in \NEMO prior to the version 3.3, and was combined with an option 
+to increase vertical mixing near the river mouth.
 
 However, with this method numerical and physical problems arise when the top grid cells are 
@@ -517,121 +760,4 @@
 
 }
-
-
-% ================================================================
-% 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 
-\href{http://climate.lanl.gov/Software/SCRIP}{here} under a copyright agreement).
-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).
 
 % ================================================================
Index: /trunk/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2541)
@@ -1053,4 +1053,6 @@
 %--------------------------------------------namtra_dmp-------------------------------------------------
 \namdisplay{namtra_dmp}
+\namdisplay{namdta_tem}
+\namdisplay{namdta_sal}
 %--------------------------------------------------------------------------------------------------------------
 
@@ -1067,7 +1069,9 @@
 The restoring term is added when \key{tradmp} is defined. 
 It also requires that both \key{dtatem} and \key{dtasal} are defined 
-($i.e.$ that $T_o$ and $S_o$ are read). The restoring coefficient 
-$\gamma$ is a three-dimensional array initialized by the user in routine 
-\rou{dtacof} also located in module \mdl{tradmp}. 
+and fill in \textit{namdta\_tem} and \textit{namdta\_sal namelists} 
+($i.e.$ that $T_o$ and $S_o$ are read using \mdl{fldread}, 
+see \S\ref{SBC_fldread}). The restoring coefficient $\gamma$ is 
+a three-dimensional array initialized by the user in routine \rou{dtacof} 
+also located in module \mdl{tradmp}. 
 
 The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} 
Index: /trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2541)
@@ -246,5 +246,4 @@
 \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 
@@ -329,9 +328,9 @@
 %--------------------------------------------------------------%
 
-% add here a description of "penetration of TKE" and the associated namelist parameters
-% \np{nn\_etau}, \np{rn\_efr}, \np{nn\_htau}
+To be add here a description of "penetration of TKE" and the associated namelist parameters
+ \np{nn\_etau}, \np{rn\_efr} and \np{nn\_htau}.
 
 % from Burchard et al OM 2008 : 
-% the most critical process not reproduced by statistical turbulence models is the activity of internal waves and their interaction with turbulence. After the Reynolds decomposition, internal waves are in principle included in the RANS equations, but later partially excluded by the hydrostatic assumption and the model resolution. Thus far, the representation of internal wave mixing in ocean models has been relatively crude (e.g. Mellor, 1989; Largeetal., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
+% the most critical process not reproduced by statistical turbulence models is the activity of internal waves and their interaction with turbulence. After the Reynolds decomposition, internal waves are in principle included in the RANS equations, but later partially excluded by the hydrostatic assumption and the model resolution. Thus far, the representation of internal wave mixing in ocean models has been relatively crude (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
 
 
@@ -355,5 +354,6 @@
 of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with
 the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}). To do so a special care 
-have to be taken for both the time and space discretization of the TKE equation \citep{Burchard_OM02}.
+have to be taken for both the time and space discretization of the TKE equation 
+\citep{Burchard_OM02,Marsaleix_al_OM08}.
 
 Let us first address the time stepping issue. Fig.~\ref{Fig_TKE_time_scheme} shows 
Index: /trunk/DOC/TexFiles/Chapters/Introduction.tex
===================================================================
--- /trunk/DOC/TexFiles/Chapters/Introduction.tex	(revision 2540)
+++ /trunk/DOC/TexFiles/Chapters/Introduction.tex	(revision 2541)
@@ -41,14 +41,14 @@
 the tendency terms of the equations are evaluated either centered  in time, or forward, 
 or backward depending of the nature of the term.
-Chapter~\ref{DOM} presents the space domain. The model is discretised on a staggered grid 
-(Arakawa C grid) with masking of land areas and uses a Leap-frog environment for time-stepping. 
-Vertical discretisation used depends on both how the bottom topography is represented and 
-whether the free surface is linear or not. Full step or partial step $z$-coordinate or 
-$s$- (terrain-following) coordinate is used with linear free surface (level position are then 
-fixed in time). In non-linear free surface, the corresponding rescaled height coordinate 
-formulation (\textit{z*} or \textit{s*}) is used (the level position then vary in time as a 
-function of the sea surface heigh). The following two chapters (\ref{TRA} and \ref{DYN}) 
-describe the discretisation of the prognostic equations for the active tracers and the 
-momentum. Explicit, split-explicit and filtered free surface formulations are implemented. 
+Chapter~\ref{DOM} presents the space domain. The model is discretised on a staggered 
+grid (Arakawa C grid) with masking of land areas. Vertical discretisation used depends 
+on both how the bottom topography is represented and whether the free surface is linear or not. 
+Full step or partial step $z$-coordinate or $s$- (terrain-following) coordinate is used 
+with linear free surface (level position are then fixed in time). In non-linear free surface, 
+the corresponding rescaled height coordinate formulation (\textit{z*} or \textit{s*}) is used 
+(the level position then vary in time as a function of the sea surface heigh). 
+The following two chapters (\ref{TRA} and \ref{DYN}) describe the discretisation of the 
+prognostic equations for the active tracers and the momentum. Explicit, split-explicit 
+and filtered free surface formulations are implemented. 
 A number of numerical schemes are available for momentum advection, for the computation 
 of the pressure gradients, as well as for the advection of tracers (second or higher 
@@ -79,12 +79,16 @@
 or \citet{Umlauf_Burchard_JMS03} mixing schemes.
 
-Chapter~\ref{OBS} describes a tool which reads in observation files (profile temperature and salinity, 
-sea surface temperature, sea level anomaly and sea ice concentration) and calculates an interpolated 
-model equivalent value at the observation location and nearest model timestep. Originally 
-developed of data assimilation, it is a fantastic tool for model and data comparison.
-Other Specific online diagnostics (not documented yet) are available in the model: output of all 
-the tendencies of the momentum and tracers equations, output of tracers tendencies 
-averaged over the time evolving mixed layer, output of the tendencies of the barotropic 
-vorticity equation, on-line floats trajectories...
+Model outputs management and specific online diagnostics are described in chapters~\ref{DIA}.
+The diagnostics includes the output of all the tendencies of the momentum and tracers equations, 
+the output of tracers tendencies averaged over the time evolving mixed layer, the output of 
+the tendencies of the barotropic vorticity equation, the computation of on-line floats trajectories... 
+Chapter~\ref{OBS} describes a tool which reads in observation files (profile temperature 
+and salinity, sea surface temperature, sea level anomaly and sea ice concentration) 
+and calculates an interpolated model equivalent value at the observation location 
+and nearest model timestep. Originally developed of data assimilation, it is a fantastic 
+tool for model and data comparison. Chapter~\ref{ASM} describes how increments 
+produced by data assimilation may be applied to the model equations.
+Finally, Chapter~\ref{CFG} provides a brief introduction to the pre-defined model 
+configurations (water column model, ORCA and GYRE families of configurations).
 
 The model is implemented in \textsc{Fortran 90}, with preprocessing (C-pre-processor). 
@@ -102,5 +106,5 @@
 around the code, the module names follow a three-letter rule. For example, \mdl{traldf} 
 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}. 
+%The complete list of module names is presented in Appendix~\ref{Apdx_D}.      %====>>>> to be done !
 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}). 
@@ -114,4 +118,5 @@
 \begin{table}[!t]	
 %\begin{center} \begin{tabular}{|p{143pt}|l|l|}	\hline
+\caption{ \label{Tab_chap}   Organization of Chapters mimicking the one of the model directories. }
 \begin{center} 	\begin{tabular}{|l|l|l|}	\hline
 Chapter \ref{STP}	& -   		        & model time STePping environment \\		\hline
@@ -123,12 +128,10 @@
 Chapter \ref{LDF}	& LDF		& Lateral DiFfusion (parameterisations) \\	\hline
 Chapter \ref{ZDF}	& ZDF		& vertical (Z) DiFfusion (parameterisations)  \\		\hline
+Chapter \ref{DIA}	& DIA		& I/O and DIAgnostics (also IOM, FLO and TRD) \\		\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
-Chapter \ref{CFG}	&  - 			& predefined configurations  \\		\hline
+Chapter \ref{ASM}	& ASM		& ASsiMilation increment  \\		\hline
+Chapter \ref{MISC}	& SOL		& Miscellaneous  topics (including solvers)  \\       \hline
+Chapter \ref{CFG}	&  - 			& predefined configurations (including C1D) \\		\hline
 \end{tabular} 	
-\caption{	 \label{Tab_chap}
-Organization of Chapters which miminc the one of the model directories. }
 \end{center}	\end{table}
 %--------------------------------------------------------------------------------------------------------------
@@ -141,5 +144,4 @@
 
 $\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; \\
@@ -150,5 +152,5 @@
 coordinate and for the new options for horizontal pressure gradient computation with 
 a non-linear equation of state.}; \\ 
-(4) more choices for the treatment of the free surface: full explicit, split-explicit and filtered. \\
+(4) more choices for the treatment of the free surface: full explicit, split-explicit or filtered schemes. \\
 (5) suppression of the rigid-lid option;\\
 (6) non linear free surface option (associated with the rescaled height coordinate  
@@ -162,18 +164,31 @@
 (12) surface module (SBC) that simplify the way the ocean is forced and include two
 bulk formulea (CLIO and CORE) and which includes an on-the-fly interpolation of input forcing fields\\
-(13) introduction of LIM 3, the new Louvain-la-Neuve sea-ice model (C-grid rheology and
+(13) RGB light penetration and optional use of ocean color 
+(14) major changes in the TKE schemes: it now includes a Langmuir cell parameterization  \citep{Axell_JGR02}, 
+the \citet{Mellor_Blumberg_JPO04} surface wave breaking parameterization, and has a time discretization 
+which is energetically consistent with the ocean model equations \citep{Burchard_OM02, Marsaleix_al_OM08}; \\
+(15) tidal mixing parametrisation (bottom intensification) + Indonesian specific tidal mixing \citep{Koch-Larrouy_al_GRL07}; \\
+(16) introduction of LIM-3, the new Louvain-la-Neuve sea-ice model (C-grid rheology and
 new thermodynamics including bulk ice salinity) \citep{Vancoppenolle_al_OM09a, Vancoppenolle_al_OM09b}
 
  \vspace{1cm}
-$\bullet$ The main modifications from NEMO/OPA v3.2 and  v3.2 are :\\
-\\
+$\bullet$ The main modifications from NEMO/OPA v3.2 and  v3.3 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"; \\
-(3) a rewriting of the bottom boundary scheme, following \citet{Campin_Goosse_Tel99}; \\
-(4) addition of the atmospheric pressure as an external forcing on both ocean and sea-ice dynamics; \\
-(5) addition of a diurnal cycle on solar radiation \citep{Bernie_al_CD07}; \\
-(6) addition of an on-line observation and model comparison (thanks to NEMOVAR project); \\
-(7) optional application of an assimilation increment (thanks to NEMOVAR project); \\
-(8) introduction of .....    
+(2) additional scheme for iso-neutral mixing \citep{Griffies_al_JPO98}, although it is still a "work in progress"; \\
+(3) a rewriting of the bottom boundary layer scheme, following \citet{Campin_Goosse_Tel99}; \\
+(4) addition of a Generic Length Scale vertical mixing scheme, following \citet{Umlauf_Burchard_JMS03};
+(5) addition of the atmospheric pressure as an external forcing on both ocean and sea-ice dynamics; \\
+(6) addition of a diurnal cycle on solar radiation \citep{Bernie_al_CD07}; \\
+(7) river runoffs added through a non-zero depth, and having its own temperature and salinity; \\
+(8) CORE II normal year forcing set as the default forcing of ORCA2-LIM configuration ; \\
+(9) generalisation of the use of \mdl{fldread} for all input fields (ocean, climatology, sea-ice damping...) 
+(10) addition of an on-line observation and model comparison (thanks to NEMOVAR project); \\
+(11) optional application of an assimilation increment (thanks to NEMOVAR project); \\
+(12) coupling interface adjusted for WRF atmospheric model
+(13) C-grid ice rheology now available fro both LIM-2 and LIM-3 \citep{Bouillon_al_OM09}; \\
+(14) a deep re-writting and simplification of the off-line tracer component (OFF\_SRC) ;  \\
+(15) the merge of passive and active advection and diffusion modules \\
+(16)  Use of the Flexible Configuration Manager (FCM) to build configurations, generate the Makefile and produce the executable ; \\
+(17) Linear-tangent and Adjoint component (TAM) added, phased with v3.0
 
  \vspace{1cm}
Index: unk/DOC/TexFiles/namlist_opa_New
===================================================================
--- /trunk/DOC/TexFiles/namlist_opa_New	(revision 2540)
+++ 	(revision )
@@ -1,513 +1,0 @@
-!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-!! OPA namelists :   1 - run manager      (namrun, nam_ctl, nam_mpp)
-!!                   2 - Domain           (nam_zgr, nam_zgr_sco, namdom)
-!!                   3 - Surface boundary (namtau, namflx, namqsr, namfwb, namalb, namcpl)
-!!                   4 - lateral boundary (namlbc, namobc, namcla)
-!!                   5 - bottom  boundary (nambfr, nambbc, nambbl)
-!!                   6 - Tracer           (nam_traadv, nam_traldf, namtdp)
-!!                   7 - dynamics         (nam_dynvor, nam_dynhpg, namcfg, nam_dynldf)
-!!                   8 - Verical physics  (namzdf, namnpc, namgap, namznl, namspr, namflo)
-!!                   9 - diagnostics      (namtrd, namgap, namznl, namptr, namspr, namflo)
-!!                   9 - miscellaneous    (solver)
-!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-!  CAUTION: some scripts does not support CAPITALs for logical definition 
-!  *******       use .true. or .false. and NOT .TRUE. or .FALSE.
-
-!!======================================================================
-!!                   ***  Run management namelists  ***
-!!======================================================================
-!!   namrun       parameters of the run
-!!   nam_ctl      Control prints & Benchmark
-!!   nam_mpp      Massively Parallel Processing
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&namrun   !   parameters of the run
-!-----------------------------------------------------------------------
-   no         =       0    !  job number
-   cexper     =  "ORCA2"   !  experience name for vairmer format
-   ln_rstart  = .false.    !  boolean term for restart (true or false)
-   nrstdt     =       0    !  restart control = 0 restart, do not control nit000 in the restart file.
-                           !                  = 1 restart, control nit000 in the restart file. Do not
-                           !                      use the date in the restart file (use ndate0 in namelist)
-                           !                  = 2 restart, control nit000 in the restart file, use the date
-                           !                      in the restart file. ndate0 in the namelist is ignored.
-   nit000     =       1    !  number of the first time step
-   nitend     =    5475    !  number of the first time step
-   ndate0     =  010101    !  initial calendar date aammjj
-   nleapy     =       0    !  Leap year calendar (0/1)
-   ninist     =       0    !  initial state output flag (0/1)
-   nstock     =    5475    !  number of the first time step
-   nwrite     =    5475    !  frequency of OUTPUT file
-   nrunoff    =       2    !  = 0 no, 1 runoff, 2 runoff+river mouth ups adv
-/
-!-----------------------------------------------------------------------
-&nam_ctl      !   Control prints & Benchmark
-!-----------------------------------------------------------------------
-   ln_ctl =  .false.       !  trends control print (expensive!)
-   nprint =       0        !  level of print (0 no print)
-   nictls =       0        !  start i indice to make the control SUM (very usefull to compare mono-
-   nictle =       0        !  end   i indice to make the control SUM (-versus multi processor runs)
-   njctls =       0        !  start j indice to make the control SUM (very usefull to compare mono-
-   njctle =       0        !  end   j indice to make the control SUM (-versus multi processor runs)
-   isplt  =       1        !  number of processors following i
-   jsplt  =       1        !  number of processors following j
-   nbench =       0        !  Bench parameter (0/1): CAUTION it must be zero except for bench
-                           !     for which we don't care about physical meaning of the results
-   nbit_cmp =     0        !  bit comparison mode parameter (0/1): enables bit comparison between 
-                           !             single and multiple processor runs
-/
-!-----------------------------------------------------------------------
-&nam_mpp      !   Massively Parallel Processing
-!-----------------------------------------------------------------------
-   c_mpi_send =  'S'       !  mpi send/recieve type   ='S', 'B', or 'I' for standard send,
-                           !  buffer blocking send or immediate non-blocking sends, resp.
-/
-
-!!======================================================================
-!!                      ***  Domain namelists  ***
-!!======================================================================
-!!   nam_zgr       vertical coordinate
-!!   nam_zgr_sco   s-coordinate or hybrid z-s-coordinate
-!!   namdom        space and time domain (bathymetry, mesh, timestep)
-!!======================================================================
-!-----------------------------------------------------------------------
-&nam_zgr       !   vertical coordinate
-!-----------------------------------------------------------------------
-   ln_zco   =  .false.     !  z-coordinate - full    steps   (T/F)      ("key_zco" may also be defined)
-   ln_zps   =  .true.      !  z-coordinate - partial steps   (T/F)
-   ln_sco   =  .false.     !  s- or hybrid z-s-coordinate    (T/F)
-/
-!-----------------------------------------------------------------------
-&nam_zgr_sco   !   s-coordinate or hybrid z-s-coordinate
-!-----------------------------------------------------------------------
-   sbot_min =  300.        !  minimum depth of s-bottom surface (>0) (m)
-   sbot_max = 5250.        !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   theta    =    6.0       !  surface control parameter (0<=theta<=20)
-   thetb    =    0.75      !  bottom control parameter  (0<=thetb<= 1)
-   r_max    =    0.15      !  maximum cut-off r-value allowed (0<r_max<1)
-/
-!-----------------------------------------------------------------------
-&namdom   !   space and time domain (bathymetry, mesh, timestep)
-!-----------------------------------------------------------------------
-   ntopo     =     1       !  = 0/1 ,compute/read the bathymetry file (mbathy, nb of T-ocean levels)
-   e3zps_min =     5.      !  the thickness of the partial step is set larger than the minimum
-   e3zps_rat =     0.1     !  of e3zps_min and e3zps_rat * e3t   (N.B. 0<e3zps_rat<1)
-   nmsh      =     0       !  =1 create a mesh file (coordinates, scale factors, masks)
-   nacc      =     0       !  the acceleration of convergence method   = 0, no acceleration, rdt = rdttra
-                           !                                           = 1, acceleration used, rdt < rdttra(k)
-   atfp      =     0.1     !  asselin time filter parameter
-   rdt       =  5760.      !  time step for the dynamics (and tracer if nacc=0)
-   rdtmin    =  5760.      !  minimum time step on tracers
-   rdtmax    =  5760.      !  maximum time step on tracers
-   rdth      =   800.      !  depth variation of tracer time step
-   rdtbt     =    90.      !  barotropic time step (for the time splitting algorithm)
-   nfice     =     5       !  frequency of ice model call
-   nfbulk    =     5       !  frequency of bulk formulea call (not used if ice used)
-   nclosea   =     0       !  = 0 no closed sea in the model domain
-                           !  = 1 closed sea (Black Sea, Caspian Sea, Great US Lakes...) 
-/
-
-!!======================================================================
-!!            ***  Surface Boundary Condition namelists  ***
-!!======================================================================
-!!   namtau      surface wind stress
-!!   namflx      surface fluxes
-!!   namqsr      penetrative solar radiation
-!!   namfwb      freshwater budget correction
-!!   namalb      albedo parameters
-!!   namcpl      coupled ocean/atmosphere model (#ifdef "key_coupled")
-!!======================================================================
-!-----------------------------------------------------------------------
-&namtau   !   surface wind stress
-!-----------------------------------------------------------------------
-   ntau000 =      0        !  gently increase the stress over the first ntau_rst time-steps
-   tau0x   =      0.e0     !  uniform value used as default surface heat flux
-   tau0y   =      0.e0     !  uniform value used as default solar radiation flux
-/
-!-----------------------------------------------------------------------
-&namflx   !   surface fluxes
-!-----------------------------------------------------------------------
-   q0      =      0.e0     !  uniform value used as default surface heat flux
-   qsr0    =      0.e0     !  uniform value used as default solar radiation flux
-   emp0    =      0.e0     !  uniform value used as default surface freswater budget (E-P)
-   dqdt0   =      -40.     !  feedback coefficient for SST damping (W/m2/K)
-   deds0   =        0.     !  feedback coefficient for SSS damping (mm/day)
-/
-!-----------------------------------------------------------------------
-&namqsr   !   penetrative solar radiation
-!-----------------------------------------------------------------------
-   ln_traqsr = .true.      !  penetrative solar radiation (T) or not (F)     (Default=T)
-   rabs     =   0.58       !  fraction of qsr associated with xsi1
-   xsi1     =   0.35       !  first depth of extinction
-   xsi2     =   23.0       !  second depth of extinction
-/
-!-----------------------------------------------------------------------
-&namfwb   !   freshwater budget correction
-!-----------------------------------------------------------------------
-   ln_fwb    = .true.      !  flag for freshwater budget correction (0 annual mean)
-/
-!-----------------------------------------------------------------------
-&namalb   !   albedo parameters
-!-----------------------------------------------------------------------
-   cgren    =      0.06    !  correction of the snow or ice albedo to take into account
-   albice   =      0.5     !  albedo of melting ice in the arctic and antarctic
-   alphd    =      0.80    !  coefficients for linear interpolation used to
-   alphc    =      0.65    !  compute albedo between two extremes values 
-   alphdi   =      0.72    !  (Pyane, 1972)
-/
-!-----------------------------------------------------------------------
-&namcpl   !   coupled ocean/atmosphere model                            ("key_coupled")
-!-----------------------------------------------------------------------
-/
-
-!!======================================================================
-!!               ***  Lateral boundary condition  ***
-!!======================================================================
-!!   namlbc        lateral momentum boundary condition
-!!   namcro        cross land advection
-!!   namobc        open boundaries parameters                           ("key_obc")
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&namlbc   !   lateral momentum boundary condition
-!-----------------------------------------------------------------------
-   shlat  =      2.        !      shlat = 0 : free slip
-                           !  0 < shlat < 2 : partial slip
-                           !      shlat = 2 : no slip
-                           !  2 < shlat     : strong slip
-/
-!-----------------------------------------------------------------------
-&namcro   !   cross land advection
-!-----------------------------------------------------------------------
-!  n_cla   advection between 2 ocean pts separates by land 
-   n_cla   = 1
-/
-!-----------------------------------------------------------------------
-&namobc   !   open boundaries parameters                            ("key_obc")
-!-----------------------------------------------------------------------
-    nobc_dta    =    0     !  = 0 the obc data are equal to the initial state
-                           !  = 1 the obc data are read in 'obc   .dta' files
-    rdpein      =    1.    !  ???
-    rdpwin      =    1.    !  ???
-    rdpnin      =   30.    !  ???
-    rdpsin      =    1.    !  ???
-    rdpeob      = 1500.    !  time relaxation (days) for the east  open boundary
-    rdpwob      =   15.    !    "        "           for the west  open boundary
-    rdpnob      =  150.    !    "        "           for the north open boundary
-    rdpsob      =   15.    !    "        "           for the south open boundary
-    zbsic1      =  140.e+6 !  barotropic stream function on first  isolated coastline
-    zbsic2      =    1.e+6 !    "                   "    on second isolated coastline
-    zbsic3      =    0.    !    "                   "    on thrid  isolated coastline
-    ln_obc_clim = .true.   !  climatological obc data files (default T)
-    ln_vol_cst  = .true.   !  total volume conserved
-/
-
-!!======================================================================
-!!                 ***  Bottom boundary condition  ***
-!!======================================================================
-!!   nambfr        bottom friction
-!!   nambbc        bottom temperature boundary condition                ("key_trabbc")
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&nambfr   !   bottom friction
-!-----------------------------------------------------------------------
-   nbotfr =       1        !  type of bottom friction :   = 0 : no slip        ,  = 2 : nonlinear friction
-                           !                              = 1 : linear friction,  = 3 : free slip
-   bfri1  =   4.e-4        !  bottom drag coefficient (linear case)
-   bfri2  =   1.e-3        !  bottom drag coefficient (non linear case)
-   bfeb2  =  2.5e-3        !  bottom turbulent kinetic energy background  (m^2/s^2)
-/
-!-----------------------------------------------------------------------
-&nambbc   !   bottom temperature boundary condition
-!-----------------------------------------------------------------------
-   ngeo_flux =  2             !  = 0 no geothermal heat flux
-                              !  = 1 constant geothermal heat flux
-                              !  = 2 variable geothermal heat flux (read in geothermal_heating.nc in mW/m2) 
-   ngeo_flux_const = 86.4e-3  !  Constant value of geothermal heat flux (W/m2)
-/
-
-!!======================================================================
-!!                        Tracer (T & S ) namelists
-!!======================================================================
-!!   nameos        equation of state
-!!   nam_traadv    advection scheme
-!!   nam_traldf    lateral diffusion scheme
-!!   nambbl        bottom boundary layer scheme                         ("key_trabbl_dif","key_trabbl_adv")
-!!   namdmp        tracer newtonian damping                             ("key_tradmp")
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&nameos   !   ocean physical parameters
-!-----------------------------------------------------------------------
-   neos   =      0         !  type of equation of state and Brunt-Vaisala frequency
-                           !     = 0, UNESCO (formulation of Jackett and McDougall (1994) and of McDougall (1987) )
-                           !     = 1, linear: rho(T)   = rau0 * ( 1.028 - ralpha * T )
-                           !     = 2, linear: rho(T,S) = rau0 * ( rbeta * S - ralpha * T )
-   ralpha =  2.e-4         !  thermal expension coefficient (neos= 1 or 2)
-   rbeta  =  0.001         !  saline  expension coefficient (neos= 2)
-/
-!-----------------------------------------------------------------------
-&nam_traadv   !   advection scheme for tracer 
-!-----------------------------------------------------------------------
-   ln_traadv_cen2   =  .true.    !  2nd order centered scheme    (default T)
-   ln_traadv_tvd    =  .false.   !  TVD scheme                   (default F)
-   ln_traadv_muscl  =  .false.   !  MUSCL scheme                 (default F)
-   ln_traadv_muscl2 =  .false.   !  MUSCL2 scheme                (default F)
-   ln_traadv_ubs    =  .false.   !  UBS scheme                   (default F)
-/
-!-----------------------------------------------------------------------
-&nam_traldf   !   lateral diffusion scheme for tracer 
-!-----------------------------------------------------------------------
-!                                !  Type of the operator : 
-   ln_traldf_lap    =  .true.    !     laplacian operator          (default T)
-   ln_traldf_bilap  =  .false.   !     bilaplacian operator        (default F)
-                                 !  Direction of action  :
-   ln_traldf_level  =  .false.   !     iso-level                   (default F)
-   ln_traldf_hor    =  .false.   !     horizontal (geopotential)   (default F)   (require "key_ldfslp" when ln_sco=T)
-   ln_traldf_iso    =  .true.    !     iso-neutral                 (default T)   (require "key_ldfslp")
-!                                !  Coefficient
-   aht0    =  2000.              !     horizontal eddy diffusivity for tracers (m2/s)
-   ahtb0   =     0.              !     background eddy diffusivity for isopycnal diffusion (m2/s)
-   aeiv0   =  2000.              !     eddy induced velocity coefficient (m2/s)  (require "key_traldf_eiv")
-/
-!-----------------------------------------------------------------------
-&nam_trabbl   !   bottom boundary layer scheme                          ("key_trabbl_dif","key_trabbl_adv")
-!-----------------------------------------------------------------------
-!  "key_trabbl	                  ! diffusive bbl 
-!  "key_trabbl_adv               ! advective bbl 
-   atrbbl = 10000.               !  lateral tracer coeff. for bottom boundary layer scheme (m2/s)
-/
-!-----------------------------------------------------------------------
-&namdmp   tracer newtonian damping ('key_tradmp')
-!-----------------------------------------------------------------------
-   ndmp   =   -1                 !  type of damping in temperature and salinity 
-                                 !     ='latitude', damping poleward of 'ndmp' degrees and function 
-                                 !                  of the distance-to-coast. Red and Med Seas as ndmp=-1
-                                 !     =-1 damping only in Med and Red Seas
-   ndmpf  =    1                 !  =1 create a damping.coeff NetCDF file (the 3D damping array)
-   nmldmp =    1                 !  type of damping in the mixed layer
-                                 !     =0 damping throughout the water column
-                                 !     =1 no damping in the mixed layer defined by avt >5cm2/s )
-                                 !     =2 no damping in the mixed layer defined rho<rho(surf)+.01 )
-   sdmp   =  50.                 !  surface time scale for internal damping (days)
-   bdmp   = 360.                 !  bottom  time scale for internal damping (days)
-   hdmp   = 800.                 !  depth of transition between sdmp and bdmp (meters)
-/
-
-!!======================================================================
-!!                      ***  Dynamics namelists  ***
-!!======================================================================
-!!   nam_dynvor    advection scheme
-!!   nam_dynhpg    hydrostatic pressure gradient
-!!   namflg        hydrostatic pressure gradient time stepping
-!!   nam_dynldf    lateral diffusion scheme
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&nam_dynvor   !   option of physics/algorithm 
-!-----------------------------------------------------------------------
-   ln_dynvor_ene = .false.    !  enstrophy conserving scheme (default T)
-   ln_dynvor_ens = .true.     !  energy conserving scheme    (default F)
-   ln_dynvor_mix = .false.    !  mixed scheme                (default F)
-   ln_dynvor_een = .false.    !  energy & enstrophy scheme   (default F)
-/
-!-----------------------------------------------------------------------
-&nam_dynhpg   !   Hydrostatic pressure gradient option
-!-----------------------------------------------------------------------
-   ln_hpg_zco = .false.       !  z-coordinate - full steps                   
-   ln_hpg_zps = .true.        !  z-coordinate - partial steps (interpolation)
-   ln_hpg_sco = .false.       !  s-coordinate (standard jacobian formulation)
-   ln_hpg_hel = .false.       !  s-coordinate (helsinki modification)
-   ln_hpg_wdj = .false.       !  s-coordinate (weighted density jacobian)
-   ln_hpg_djc = .false.       !  s-coordinate (Density Jacobian with Cubic polynomial)
-   ln_hpg_rot = .false.       !  s-coordinate (ROTated axes scheme)
-   gamm       = 0.e0          !  weighting coefficient (wdj scheme)
-/
-!-----------------------------------------------------------------------
-&namflg   !   algorithm flags (algorithm not control by CPP keys)
-!-----------------------------------------------------------------------
-   ln_dynhpg_imp   =  .false. !  hydrostatic pressure gradient: semi-implicit time scheme  (T)
-                              !                                 centered      time scheme  (F)
-   nn_dynhpg_rst   =  0       !  add dynhpg implicit variables in restart ot not (1/0)
-/
-!-----------------------------------------------------------------------
-!nam_dynspg   !   surface pressure gradient   (CPP key only)
-!-----------------------------------------------------------------------
-!                                !  explicit free surface               ("key_dynspg_exp")
-!                                !  filtered free surface               ("key_dynspg_flt")
-!                                !  split-explicit free surface         ("key_dynspg_ts")
-!                                !  rigid-lid                           ("key_dynspg_rl")
-
-!-----------------------------------------------------------------------
-&nam_dynldf   !   lateral diffusion on momentum
-!-----------------------------------------------------------------------
-                                 !  Type of the operator : 
-   ln_dynldf_lap    =  .true.    !     laplacian operator          (default T)
-   ln_dynldf_bilap  =  .false.   !     bilaplacian operator        (default F)
-                                 !  Direction of action  : 
-   ln_dynldf_level  =  .false.   !     iso-level                   (default F)
-   ln_dynldf_hor    =  .true.    !     horizontal (geopotential)   (default F)   (require "key_ldfslp" in s-coord.)
-   ln_dynldf_iso    =  .false.   !     iso-neutral                 (default T)   (require "key_ldfslp")
-                                 !  Coefficient
-   ahm0    = 40000.              !     horizontal eddy viscosity for the dynamics (m2/s)
-   ahmb0   =     0.              !     background eddy viscosity for isopycnal diffusion (m2/s)
-/
-
-!!======================================================================
-!!             Tracers & Dynamics vertical physics namelists
-!!======================================================================
-!!       namzdf        vertical physics
-!!       namnpc        non penetrative convection                       
-!!       nambbl        bottom boundary layer scheme                     ("key_trabbl_dif", "key_trabbl_adv")
-!!       namric        richardson number dependent vertical mixing      ("key_zdfric"      )
-!!       namtke        TKE dependent vertical mixing                    ("key_zdftke"      )
-!!       namkpp        KPP dependent vertical mixing                    ("key_zdfkpp"      )
-!!       namddm        double diffusive mixing parameterization         ("key_zdfddm"      )
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&namzdf   !   vertical physics
-!-----------------------------------------------------------------------
-   ln_zdfnpc = .false.        !  Non-Penetrative Convection          (default T)
-   avm0      = 1.2e-4         !  vertical eddy viscosity   for the dynamic (m2/s)   (background Kz if not "key_zdfcst")
-   avt0      = 1.2e-5         !  vertical eddy diffusivity for tracers     (m2/s)   (background Kz if not "key_zdfcst")
-   ln_zdfevd = .true.         !  enhanced vertical diffusion         (default T)
-   avevd     =   100.         !  vertical coefficient for enhanced diffusion scheme (m2/s)
-   n_evdm    =     0          !  enhanced mixing apply on tracer (=0) or on tracer and momentum (=1)
-   ln_zdfexp =  .false.       !  =T/F  split explicit / implicit      (Default=F)
-   n_zdfexp  =     3          !  number of sub-timestep for ln_zdfexp=T
-/
-!-----------------------------------------------------------------------
-&namnpc   !   non penetrative convection
-!-----------------------------------------------------------------------
-   nnpc1  =      1            !  non penetrative convective scheme computation frequency
-   nnpc2  =    365            !  non penetrative convective scheme print frequency
-/
-!-----------------------------------------------------------------------
-&namric   richardson number dependent vertical diffusion                ("key_zdfric" )
-!-----------------------------------------------------------------------
-   avmri = 100.e-4            !  maximum value of the vertical viscosity
-   alp   =      5.            !  coefficient of the parameterization
-   nric  =       2            !  coefficient of the parameterization
-/
-!-----------------------------------------------------------------------
-&namtke   !   turbulent eddy kinetic dependent vertical diffusion       ("key_zdftke" )
-!-----------------------------------------------------------------------
-   ln_rstke = .false.         !  flag to restart with tke from a run without tke (default F)
-   ediff =       0.1          !  coef. to compute vertical eddy coef. (avt=ediff*mxl*sqrt(e) )
-   ediss =       0.7          !  coef. of the Kolmogoroff dissipation
-   ebb   =      3.75          !  coef. of the surface input of tke
-   efave =        1.          !  coef. to applied to the tke diffusion ( avtke=efave*avm )
-   emin  =     1.e-6          !  minimum value of tke (m^2/s^2)
-   emin0 =     1.e-4          !  surface minimum value of tke (m^2/s^2)
-   nitke =        50          !  number of restart iterative loops
-!! ri_c  =         0.22222222 !  critic richardson number (default = 2/9)
-   nmxl  =         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
-   npdl  =         1          !  prandtl number   = 0 no vertical prandtl number (avt=avm)
-                              !                   = 1 prandtl number function of richarson number (avt=pdl*avm)
-                              !                   = 2 same as = 1 but a shapiro filter is applied on pdl
-   nave  =         1          !  horizontal averaged (=1) or not (=0) on avt and amv   (default =1)
-   navb  =         0          !  constant (=0) or profile (=1) background on avt
-/
-!-----------------------------------------------------------------------
-&namkpp   !   K-Profile Parameterization dependent vertical mixing      ("key_zdfkpp", and optionnally:
-!                                                                        "key_kppcustom" or "key_kpplktb")
-!-----------------------------------------------------------------------
-   ln_kpprimix  = .true.      !  shear instability mixing  (default T)
-   difmiw       =  1.0e-04    !  constant internal wave viscosity (m2/s)
-   difsiw       =  0.1e-04    !  constant internal wave diffusivity (m2/s)
-   Riinfty      =  0.8        !  local Richardson Number limit for shear instability
-   difri        =  0.0050     !  maximum shear mixing at Rig = 0    (m2/s)
-   bvsqcon      = -0.01e-07   !  Brunt-Vaisala squared (1/s**2) for maximum convection
-   difcon       =  1.         !  maximum mixing in interior convection (m2/s)
-   navb         =  0          !  horizontal averaged (=1) or not (=0) on avt and amv   (default =1)
-   nave         =  1          !  constant (=0) or profile (=1) background on avt
-/
-!-----------------------------------------------------------------------
-&namddm   !   double diffusive mixing parameterization                  ("key_zdfddm")
-!-----------------------------------------------------------------------
-      avts  = 1.e-4           !  maximum avs (vertical mixing on salinity)
-      hsbfr = 1.6             !  heat/salt buoyancy flux ratio
-/
-
-!!======================================================================
-!!                  ***  Miscelaneous namelists  ***
-!!======================================================================
-!!   namsol      elliptic solver / island / free surface 
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&namsol   !   elliptic solver / island / free surface 
-!-----------------------------------------------------------------------
-   nsolv     =      1         !  type of elliptic solver  =1 preconditioned conjugate gradient (pcg)
-                              !                           =2 successive-over-relaxation (sor)
-                              !                           =3 FETI (fet)                               ("key_feti")
-                              !                           =4 sor with extra outer halo
-   nsol_arp  =      0         !  absolute/relative (0/1) precision convergence test
-   nmin      =    300         !  minimum of iterations for the SOR solver
-   nmax      =    800         !  maximum of iterations for the SOR solver
-   nmod      =     10         !  frequency of test for the SOR solver
-   eps       =  1.E-6         !  absolute precision of the solver
-   resmax    = 1.E-10         !  absolute precision for the SOR solver
-   sor       =   1.92         !  optimal coefficient for SOR solver
-   epsisl    = 1.e-10         !  absolute precision on stream function solver
-   nmisl     =   4000         !  maximum pcg iterations for island                                    ("key_islands")
-   rnu       =     1.         !  strength of the additional force used in filtered free surface
-/
-
-!!======================================================================
-!!                       ***  Diagnostics namelists  ***
-!!======================================================================
-!!   namtrd       dynamics and/or tracer trends                         ("key_trddyn","key_trdtra","key_trdmld")
-!!   namgap       level mean model-data gap                             ("key_diagap")
-!!   namptr       Poleward Transport Diagnostics                        ("key_diaspr")
-!!   namspr       surface pressure diagnosed in rigid-lid               ("key_diaspr")
-!!   namflo       float parameters                                      ("key_float")
-!!======================================================================
-
-!-----------------------------------------------------------------------
-&namtrd    !   diagnostics on dynamics and/or tracer trends             ("key_trddyn" and/or "key_trdtra")
-!                          or mixed-layer trends                        ('key_trdmld')
-!                          or barotropic vorticity                      ("key_trdvor")
-!-----------------------------------------------------------------------
-   ntrd  = 365                   !  time step frequency dynamics and tracers trends
-   nctls =   0                   !  control surface type in mixed-layer trends (0,1 or n<jpk)
-   ln_trdmld_restart = .false.   !  restart for ML diagnostics
-   ucf   =  1.                   !  unit conversion factor (=1 -> /seconds | =86400. -> /day)
-   ln_trdmld_instant = .false.   !  flag to diagnose trends of instantantaneous or mean ML T/S
-/
-!-----------------------------------------------------------------------
-&namgap    !   level mean model-data gap                                ('key_diagap')
-!-----------------------------------------------------------------------
-   ngap =  15                    !  time-step frequency of model-data gap computation
-   nprg =  10                    !  time-step frequency of gap print in model output
-/
-!-----------------------------------------------------------------------
-&namptr   !   Poleward Transport Diagnostic
-!-----------------------------------------------------------------------
-   ln_diaptr = .false.           !  logical flag for Poleward transport computation
-   ln_subbas = .false.           !  logical flag for Atlantic/Pacific/Indian basins computation
-                                 !  need input basins mask file named "subbasins.nc"
-   nf_ptr    =  15               !  Frequency of computation
-/
-!-----------------------------------------------------------------------
-&namspr  !   surface pressure diagnostic
-!-----------------------------------------------------------------------
-   nmaxp   =   1000              !  maximum of iterations for the solver
-   epsp    =  1.e-3              !  absolute precision of the solver
-   niterp  =    400              !  number of iteration done by the solver
-/
-!-----------------------------------------------------------------------
-&namflo    !   float parameters                                         ("key_float")
-!-----------------------------------------------------------------------
-    ln_rstflo = .false.          !  boolean term for float restart (true or false)
-    nwritefl  =      75          !  frequency of float output file 
-    nstockfl  =    5475          !  frequency of float restart file 
-    ln_argo   = .false.          !  Argo type floats (stay at the surface each 10 days)
-    ln_flork4 = .false.          !  = T trajectories computed with a 4th order Runge-Kutta
-                                 !  = F  (default)   computed with Blanke' scheme
-/
