Index: branches/2017/dev_merge_2017/DOC/tex_sub/annex_C.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/annex_C.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/annex_C.tex	(revision 9392)
@@ -363,5 +363,5 @@
 %       Vorticity Term with ENE scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity Term with ENE scheme (\protect\np{ln\_dynvor\_ene}=.true.)}
+\subsubsection{Vorticity Term with ENE scheme (\protect\np{ln_dynvor_ene}=.true.)}
 \label{Apdx_C_vorENE} 
 
@@ -400,5 +400,5 @@
 %       Vorticity Term with EEN scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity Term with EEN scheme (\protect\np{ln\_dynvor\_een}=.true.)}
+\subsubsection{Vorticity Term with EEN scheme (\protect\np{ln_dynvor_een}=.true.)}
 \label{Apdx_C_vorEEN} 
 
@@ -883,5 +883,5 @@
 %       Vorticity Term with ENS scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity Term with ENS scheme  (\protect\np{ln\_dynvor\_ens}=.true.)}
+\subsubsection{Vorticity Term with ENS scheme  (\protect\np{ln_dynvor_ens}=.true.)}
 \label{Apdx_C_vorENS} 
 
@@ -943,5 +943,5 @@
 %       Vorticity Term with EEN scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity Term with EEN scheme (\protect\np{ln\_dynvor\_een}=.true.)}
+\subsubsection{Vorticity Term with EEN scheme (\protect\np{ln_dynvor_een}=.true.)}
 \label{Apdx_C_vorEEN} 
 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/annex_E.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/annex_E.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/annex_E.tex	(revision 9392)
@@ -19,5 +19,5 @@
 %        UBS scheme  
 % -------------------------------------------------------------------------------------------------------------
-\section{Upstream Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs}=T)}
+\section{Upstream Biased Scheme (UBS) (\protect\forcode{ln_traadv_ubs = .true.})}
 \label{TRA_adv_ubs}
 
@@ -59,5 +59,5 @@
 where the control of artificial diapycnal fluxes is of paramount importance. 
 It has therefore been preferred to evaluate the vertical flux using the TVD 
-scheme when \np{ln\_traadv\_ubs}=T.
+scheme when \forcode{ln_traadv_ubs = .true.}.
 
 For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex	(revision 9392)
@@ -15,11 +15,11 @@
 
 Two scheme are available to perform the iso-neutral diffusion. 
-If the namelist logical \np{ln\_traldf\_triad} is set true, 
+If the namelist logical \np{ln_traldf_triad} is set true, 
 \NEMO updates both active and passive tracers using the Griffies triad representation 
 of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes. 
-If the namelist logical \np{ln\_traldf\_iso} is set true, 
+If the namelist logical \np{ln_traldf_iso} is set true, 
 the filtered version of Cox's original scheme (the Standard scheme) is employed (\S\ref{LDF_slp}). 
 In the present implementation of the Griffies scheme, 
-the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false.
+the advective skew fluxes are implemented even if \np{ln_traldf_eiv} is false.
 
 Values of iso-neutral diffusivity and GM coefficient are set as
@@ -31,24 +31,24 @@
 The options specific to the Griffies scheme include:
 \begin{description}[font=\normalfont]
-\item[\np{ln\_triad\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then
+\item[\np{ln_triad_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then
   `iso-neutral' mixing is accomplished within the surface mixed-layer
   along slopes linearly decreasing with depth from the value immediately below
   the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}). 
   This is the same treatment as used in the default implementation \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}.  
-  Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced 
+  Where \np{ln_triad_iso} is set true, the vertical skew flux is further reduced 
   to ensure no vertical buoyancy flux, giving an almost pure
   horizontal diffusive tracer flux within the mixed layer. This is similar to
   the tapering suggested by \citet{Gerdes1991}. See \S\ref{sec:triad:Gerdes-taper}
-\item[\np{ln\_botmix\_triad}] See \S\ref{sec:triad:iso_bdry}. 
+\item[\np{ln_botmix_triad}] See \S\ref{sec:triad:iso_bdry}. 
   If this is set false (the default) then the lateral diffusive fluxes
   associated with triads partly masked by topography are neglected. 
   If it is set true, however, then these lateral diffusive fluxes are applied, 
   giving smoother bottom tracer fields at the cost of introducing diapycnal mixing.
-\item[\np{rn\_sw\_triad}]  blah blah to be added....
+\item[\np{rn_sw_triad}]  blah blah to be added....
 \end{description}
 The options shared with the Standard scheme include:
 \begin{description}[font=\normalfont]
-\item[\np{ln\_traldf\_msc}]   blah blah to be added
-\item[\np{rn\_slpmax}]  blah blah to be added
+\item[\np{ln_traldf_msc}]   blah blah to be added
+\item[\np{rn_slpmax}]  blah blah to be added
 \end{description}
 \section{Triad formulation of iso-neutral diffusion}
@@ -651,11 +651,11 @@
 or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is
 masked. The associated lateral fluxes (grey-black dashed line) are
-masked if \np{ln\_botmix\_triad}=false, but left unmasked,
-giving bottom mixing, if \np{ln\_botmix\_triad}=true.
-
-The default option \np{ln\_botmix\_triad}=false is suitable when the
-bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}=1),
+masked if \forcode{ln_botmix_triad = .false.}, but left unmasked,
+giving bottom mixing, if \forcode{ln_botmix_triad = .true.}.
+
+The default option \forcode{ln_botmix_triad = .false.} is suitable when the
+bbl mixing option is enabled (\key{trabbl}, with \forcode{nn_bbl_ldf = 1}),
 or  for simple idealized  problems. For setups with topography without
-bbl mixing, \np{ln\_botmix\_triad}=true may be necessary.
+bbl mixing, \forcode{ln_botmix_triad = .true.} may be necessary.
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[h] \begin{center}
@@ -674,6 +674,6 @@
       or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point
       is masked. The associated lateral fluxes (grey-black dashed
-      line) are masked if \protect\np{botmix\_triad}=.false., but left
-      unmasked, giving bottom mixing, if \protect\np{botmix\_triad}=.true.}
+      line) are masked if \protect\np{botmix_triad}=.false., but left
+      unmasked, giving bottom mixing, if \protect\np{botmix_triad}=.true.}
  \end{center} \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -710,5 +710,5 @@
 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper}
 This is the option activated by the default choice
-\np{ln\_triad\_iso}=false. Slopes $\tilde{r}_i$ relative to
+\forcode{ln_triad_iso = .false.}. Slopes $\tilde{r}_i$ relative to
 geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the
 surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values
@@ -842,5 +842,5 @@
   components}
 \label{sec:triad:Gerdes-taper}
-The alternative option is activated by setting \np{ln\_triad\_iso} =
+The alternative option is activated by setting \np{ln_triad_iso} =
   true. This retains the same tapered slope $\rML$  described above for the
 calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the
@@ -917,5 +917,5 @@
 it to the Eulerian velocity prior to computing the tracer
 advection. This is implemented if \key{traldf\_eiv} is set in the
-default implementation, where \np{ln\_traldf\_triad} is set
+default implementation, where \np{ln_traldf_triad} is set
 false. This allows us to take advantage of all the advection schemes
 offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$
@@ -924,5 +924,5 @@
 paramount importance.
 
-However, when \np{ln\_traldf\_triad} is set true, \NEMO instead
+However, when \np{ln_traldf_triad} is set true, \NEMO instead
 implements eddy induced advection according to the so-called skew form
 \citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes
@@ -1123,5 +1123,5 @@
 and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the
 $i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$
-$u$-point is masked. The namelist parameter \np{ln\_botmix\_triad} has
+$u$-point is masked. The namelist parameter \np{ln_botmix_triad} has
 no effect on the eddy-induced skew-fluxes.
 
@@ -1138,5 +1138,5 @@
 option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the
 slopes used to calculate the eddy-induced fluxes is
-unaffected by the value of \np{ln\_triad\_iso}.
+unaffected by the value of \np{ln_triad_iso}.
 
 The justification for this linear slope tapering is that, for $A_e$
@@ -1153,5 +1153,5 @@
 
 \subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag}
-Where the namelist parameter \np{ln\_traldf\_gdia}=true, diagnosed
+Where the namelist parameter \forcode{ln_traldf_gdia = .true.}, diagnosed
 mean eddy-induced velocities are output. Each time step,
 streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_ASM.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_ASM.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_ASM.tex	(revision 9392)
@@ -30,5 +30,5 @@
 Direct initialization (DI) refers to the instantaneous correction
 of the model background state using the analysis increment.
-DI is used when \np{ln\_asmdin} is set to true.
+DI is used when \np{ln_asmdin} is set to true.
 
 \section{Incremental Analysis Updates}
@@ -40,5 +40,5 @@
 is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}.
 IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI.
-IAU is used when \np{ln\_asmiau} is set to true.
+IAU is used when \np{ln_asmiau} is set to true.
 
 With IAU, the model state trajectory ${\bf x}$ in the assimilation window 
@@ -117,5 +117,5 @@
 integration \citep{Talagrand_JAS72, Dobricic_al_OS07}. Diffusion coefficients are defined as 
 $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$. The divergence damping is activated by
-assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist a value greater than zero. 
+assigning to \np{nn_divdmp} in the \textit{nam\_asminc} namelist a value greater than zero. 
 By choosing this value to be of the order of 100 the increments in the vertical velocity will 
 be significantly reduced.
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_CONFIG.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_CONFIG.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_CONFIG.tex	(revision 9392)
@@ -86,7 +86,7 @@
 the LIM sea-ice model (ORCA-LIM) and possibly with PISCES biogeochemical model 
 (ORCA-LIM-PISCES), using various resolutions.
-An appropriate namelist is available in \textit{CONFIG/ORCA2\_LIM3\_PISCES/EXP00/namelist\_cfg} 
+An appropriate namelist is available in \path{CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg} 
 for ORCA2.
-The domain of ORCA2 configuration is defined in ORCA\_R2\_zps\_domcfg.nc file, this file is available in tar file in the wiki of NEMO : \\
+The domain of ORCA2 configuration is defined in \ifile{ORCA\_R2\_zps\_domcfg} file, this file is available in tar file in the wiki of NEMO : \\
 https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ReferenceConfigurations/ORCA2\_LIM3\_PISCES \\
 In this namelist\_cfg the name of domain input file is set in \ngn{namcfg} block of namelist. 
@@ -156,5 +156,5 @@
 horizontal resolution. The value of the resolution is given by the resolution at the Equator 
 expressed in degrees. Each of configuration is set through the \textit{domain\_cfg} domain configuration file, 
-which sets the grid size and configuration name parameters. The NEMO System Team provides only ORCA2 domain input file "ORCA\_R2\_zps\_domcfg.nc" file  (Tab. \ref{Tab_ORCA}).
+which sets the grid size and configuration name parameters. The NEMO System Team provides only ORCA2 domain input file "\ifile{ORCA\_R2\_zps\_domcfg}" file  (Tab. \ref{Tab_ORCA}).
 
 
@@ -164,5 +164,5 @@
 \begin{table}[!t]     \begin{center}
 \begin{tabular}{p{4cm} c c c c}
-Horizontal Grid         	             & \np{ORCA\_index} &  \np{jpiglo} & \np{jpjglo} &       \\  
+Horizontal Grid         	             & \np{ORCA_index} &  \np{jpiglo} & \np{jpjglo} &       \\  
 \hline  \hline
 \~4\deg	   &        4         &         92     &      76      &       \\
@@ -214,5 +214,5 @@
 
 ORCA\_R2 pre-defined configuration can also be run with an AGRIF zoom over the Agulhas 
-current area ( \key{agrif}  defined) and, by setting the appropriate variables, see \textit{CONFIG/SHARED/namelist\_ref}
+current area ( \key{agrif}  defined) and, by setting the appropriate variables, see \path{CONFIG/SHARED/namelist_ref}
 a regional Arctic or peri-Antarctic configuration is extracted from an ORCA\_R2 or R05 configurations
 using sponge layers at open boundaries. 
@@ -252,7 +252,7 @@
 
 The GYRE configuration is set like an analytical configuration. Through \np{ln\_read\_cfg\textit{=false}} in \textit{namcfg} namelist defined in the reference configuration \textit{CONFIG/GYRE/EXP00/namelist\_cfg} anaylitical definition of grid in GYRE is done in usrdef\_hrg, usrdef\_zgr routines. Its horizontal resolution 
-(and thus the size of the domain) is determined by setting \np{nn\_GYRE} in  \ngn{namusr\_def}: \\
-\np{jpiglo} $= 30 \times$ \np{nn\_GYRE} + 2   \\
-\np{jpjglo} $= 20 \times$ \np{nn\_GYRE} + 2   \\
+(and thus the size of the domain) is determined by setting \np{nn_GYRE} in  \ngn{namusr\_def}: \\
+\np{jpiglo} $= 30 \times$ \np{nn_GYRE} + 2   \\
+\np{jpjglo} $= 20 \times$ \np{nn_GYRE} + 2   \\
 Obviously, the namelist parameters have to be adjusted to the chosen resolution, see the Configurations 
 pages on the NEMO web site (Using NEMO\/Configurations) .
@@ -296,5 +296,5 @@
 In addition to the tidal boundary condition the model may also take
 open boundary conditions from a North Atlantic model. Boundaries may be
-completely omitted by setting \np{ln\_bdy} to false.
+completely omitted by setting \np{ln_bdy} to false.
 Sample surface fluxes, river forcing and a sample initial restart file
 are included to test a realistic model run. The Baltic boundary is
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIA.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIA.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIA.tex	(revision 9392)
@@ -217,5 +217,5 @@
 For example:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <field_definition>
       <!-- T grid -->
@@ -226,5 +226,5 @@
       ...
    </field_definition> 
-\end{xmlcode}
+\end{xmllines}
 Note your definition must be added to the field\_group whose reference grid is consistent 
 with the size of the array passed to iomput. 
@@ -233,7 +233,7 @@
 or defined in the domain\_def.xml file. $e.g.$:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
      <grid id="grid_T_3D" domain_ref="grid_T" axis_ref="deptht"/>
-\end{xmlcode}
+\end{xmllines}
 Note, if your array is computed within the surface module each nn\_fsbc time\_step, 
 add the field definition within the field\_group defined with the id ''SBC'': $<$field\_group id=''SBC''...$>$ 
@@ -242,5 +242,5 @@
 \item[4.] add your field in one of the output files defined in iodef.xml (again see subsequent sections for syntax and rules)   \\
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <file id="file1" .../>   
       ...
@@ -248,5 +248,5 @@
       ...
    </file>   
-\end{xmlcode}
+\end{xmllines}
 
 \end{description}
@@ -398,10 +398,10 @@
 example 1: Direct inheritance.
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <field_definition operation="average" >
      <field id="sst"                    />   <!-- averaged      sst --> 
      <field id="sss" operation="instant"/>   <!-- instantaneous sss --> 
    </field_definition> 
-\end{xmlcode}
+\end{xmllines}
 The field ''sst'' which is part (or a child) of the field\_definition will inherit the value ''average'' 
 of the attribute ''operation'' from its parent. Note that a child can overwrite 
@@ -411,5 +411,5 @@
 example 2: Inheritance by reference.
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <field_definition>
      <field id="sst" long_name="sea surface temperature" />   
@@ -423,5 +423,5 @@
      </file>   
    </file_definition> 
-\end{xmlcode}
+\end{xmllines}
 Inherit (and overwrite, if needed) the attributes of a tag you are refering to.
 
@@ -433,5 +433,5 @@
 Note that for the field ''toce'', we overwrite the grid definition inherited from the group by ''grid\_T\_3D''.
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <field_group id="grid_T" grid_ref="grid_T_2D">
     <field id="toce" long_name="temperature"             unit="degC" grid_ref="grid_T_3D"/>
@@ -440,5 +440,5 @@
     <field id="ssh"  long_name="sea surface height"      unit="m"                        />
          ...
-\end{xmlcode}
+\end{xmllines}
 
 Secondly, the group can be used to replace a list of elements. 
@@ -446,5 +446,5 @@
 For example, a short list of the usual variables related to the U grid:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <field_group id="groupU" >
     <field field_ref="uoce"  />
@@ -452,13 +452,13 @@
     <field field_ref="utau"  />
    </field_group>
-\end{xmlcode}
+\end{xmllines}
 that can be directly included in a file through the following syntax:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <file id="myfile_U" output_freq="1d" />   
     <field_group group_ref="groupU"/>  
     <field field_ref="uocetr_eff"  />  <!-- add another field -->
    </file>   
-\end{xmlcode}
+\end{xmllines}
 
 \subsection{Detailed functionalities }
@@ -473,15 +473,15 @@
 of a 5 by 5 box with the bottom left corner at point (10,10).
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <domain_group id="grid_T">
     <domain id="myzoom" zoom_ibegin="10" zoom_jbegin="10" zoom_ni="5" zoom_nj="5" />
-\end{xmlcode}
+\end{xmllines}
 The use of this subdomain is done through the redefinition of the attribute domain\_ref of the tag family field. For example:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <file id="myfile_vzoom" output_freq="1d" >
       <field field_ref="toce" domain_ref="myzoom"/>
    </file>
-\end{xmlcode}
+\end{xmllines}
 Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain. 
 The Equatorial section, the TAO, RAMA and PIRATA moorings are alredy registered in the code 
@@ -491,9 +491,9 @@
 by ''T'' (for example: ''8s137eT'', ''1.5s80.5eT'' ...)
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <file id="myfile_vzoom" output_freq="1d" >
       <field field_ref="toce" domain_ref="0n180wT"/>
    </file>
-\end{xmlcode}
+\end{xmllines}
 Note that if the domain decomposition used in XIOS cuts the subdomain in several parts and if you use the ''multiple\_file'' type for your output files, you will endup with several files you will need to rebuild using unprovided tools (like ncpdq and ncrcat, \href{http://nco.sourceforge.net/nco.html#Concatenation}{see nco manual}). We are therefore advising to use the ''one\_file'' type in this case.
 
@@ -501,16 +501,16 @@
 Vertical zooms are defined through the attributs zoom\_begin and zoom\_end of the tag family axis. It must therefore be done in the axis part of the XML file. For example, in NEMOGCM/CONFIG/ORCA2\_LIM/iodef\_demo.xml, we provide the following example:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <axis_group id="deptht" long_name="Vertical T levels" unit="m" positive="down" >
       <axis id="deptht" />
       <axis id="deptht_myzoom" zoom_begin="1" zoom_end="10" />
-\end{xmlcode}
+\end{xmllines}
 The use of this vertical zoom is done through the redefinition of the attribute axis\_ref of the tag family field. For example:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <file id="myfile_hzoom" output_freq="1d" >
       <field field_ref="toce" axis_ref="deptht_myzoom"/>
    </file>
-\end{xmlcode}
+\end{xmllines}
 
 \subsubsection{Control of the output file names}
@@ -518,5 +518,5 @@
 The output file names are defined by the attributs ''name'' and ''name\_suffix'' of the tag family file. for example:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
    <file_group id="1d" output_freq="1d" name="myfile_1d" > 
       <file id="myfileA" name_suffix="_AAA" > <!-- will create file "myfile_1d_AAA"  -->
@@ -527,5 +527,5 @@
       </file>
    </file_group>
-\end{xmlcode}
+\end{xmllines}
 However it is often very convienent to define the file name with the name of the experiment, the output file frequency and the date of the beginning and the end of the simulation (which are informations stored either in the namelist or in the XML file). To do so, we added the following rule: if the id of the tag file is ''fileN''(where N = 1 to 999 on 1 to 3 digits) or one of the predefined sections or moorings (see next subsection), the following part of the name and the name\_suffix (that can be inherited) will be automatically replaced by:\\
 \\
@@ -589,11 +589,11 @@
    \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
-    \multicolumn{2}{|c|}{ptrc\_T}           & freq\_op & \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\
-   \hline
-    \multicolumn{2}{|c|}{diad\_T}           & freq\_op & \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\
+    \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
+    \multicolumn{2}{|c|}{ptrc\_T}           & freq\_op & \np{rn_rdt} $\times$ \np{nn_dttrc} \\
+   \hline
+    \multicolumn{2}{|c|}{diad\_T}           & freq\_op & \np{rn_rdt} $\times$ \np{nn_dttrc} \\
    \hline
     \multicolumn{2}{|c|}{EqT, EqU, EqW} & jbegin, ni,      & according to the grid    \\
@@ -613,9 +613,9 @@
 
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
  <field field\_ref="sst"  name="tosK"  unit="degK" > sst + 273.15 </field>
  <field field\_ref="taum" name="taum2" unit="N2/m4" long\_name="square of wind stress module" > taum * taum </field>
  <field field\_ref="qt"   name="stupid\_check" > qt - qsr - qns </field>
-\end{xmlcode}
+\end{xmllines}
 
 (2) Simple computation: define a new variable and use it in the file definition.
@@ -623,12 +623,12 @@
 in field\_definition:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
  <field id="sst2" long\_name="square of sea surface temperature" unit="degC2" >  sst * sst </field >
-\end{xmlcode}
+\end{xmllines}
 in file\_definition:
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
  <field field\_ref="sst2" > sst2 </field>
-\end{xmlcode}
+\end{xmllines}
 Note that in this case, the following syntaxe $<$field field\_ref="sst2" /$>$ is not working as sst2 won't be evaluated.
 
@@ -636,10 +636,10 @@
 
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
      <!-- force to keep real 8 -->
  <field field\_ref="sst" name="tos\_r8" prec="8" />
       <!-- integer 2  with add\_offset and scale\_factor attributes -->
  <field field\_ref="sss" name="sos\_i2" prec="2" add\_offset="20." scale\_factor="1.e-3" />
-\end{xmlcode}
+\end{xmllines}
 Note that, then the code is crashing, writting real4 variables forces a numerical convection from real8 to real4 which will create an internal error in NetCDF and will avoid the creation of the output files. Forcing double precision outputs with prec="8" (for example in the field\_definition) will avoid this problem.
 
@@ -647,5 +647,5 @@
 
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
       <file\_group id="1d" output\_freq="1d" output\_level="10" enabled=".TRUE."> <!-- 1d files --> 
 	<file id="file1" name\_suffix="\_grid\_T" description="ocean T grid variables" >
@@ -658,5 +658,5 @@
        </file>
      </file\_group> 
-\end{xmlcode}
+\end{xmllines}
 
 (5) use of the ``@'' function: example 1, weighted temporal average
@@ -664,10 +664,10 @@
  - define a new variable in field\_definition
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
  <field id="toce\_e3t" long\_name="temperature * e3t" unit="degC*m" grid\_ref="grid\_T\_3D" > toce * e3t </field >
-\end{xmlcode}
+\end{xmllines}
  - use it when defining your file.  
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
 <file\_group id="5d" output\_freq="5d"  output\_level="10" enabled=".TRUE." >  <!-- 5d files -->  
  <file id="file1" name\_suffix="\_grid\_T" description="ocean T grid variables" >
@@ -675,5 +675,5 @@
  </file>
 </file\_group> 
-\end{xmlcode}
+\end{xmllines}
 The freq\_op="5d" attribute is used to define the operation frequency of the ``@'' function: here 5 day. The temporal operation done by the ``@'' is the one defined in the field definition: here we use the default, average. So, in the above case, @toce\_e3t will do the 5-day mean of toce*e3t. Operation="instant" refers to the temporal operation to be performed on the field''@toce\_e3t / @e3t'': here the temporal average is alreday done by the ``@'' function so we just use instant to do the ratio of the 2 mean values. field\_ref="toce" means that attributes not explicitely defined, are inherited from toce field. Note that in this case, freq\_op must be equal to the file output\_freq.
 
@@ -682,10 +682,10 @@
  - define a new variable in field\_definition
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
  <field id="ssh2" long\_name="square of sea surface temperature" unit="degC2" >  ssh * ssh </field >
-\end{xmlcode}
+\end{xmllines}
  - use it when defining your file.  
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
 <file\_group id="1m" output\_freq="1m"  output\_level="10" enabled=".TRUE." >  <!-- 1m files -->  
  <file id="file1" name\_suffix="\_grid\_T" description="ocean T grid variables" >
@@ -693,5 +693,5 @@
  </file>
 </file\_group> 
-\end{xmlcode}
+\end{xmllines}
 The freq\_op="1m" attribute is used to define the operation frequency of the ``@'' function: here 1 month. The temporal operation done by the ``@'' is the one defined in the field definition: here we use the default, average. So, in the above case, @ssh2 will do the monthly mean of ssh*ssh. Operation="instant" refers to the temporal operation to be performed on the field ''sqrt( @ssh2 - @ssh * @ssh )'': here the temporal average is alreday done by the ``@'' function so we just use instant. field\_ref="ssh" means that attributes not explicitely defined, are inherited from ssh field. Note that in this case, freq\_op must be equal to the file output\_freq.
 
@@ -700,11 +700,11 @@
  - define 2 new variables in field\_definition
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
  <field id="sstmax" field\_ref="sst" long\_name="max of sea surface temperature" operation="maximum" />
  <field id="sstmin" field\_ref="sst" long\_name="min of sea surface temperature" operation="minimum" />
-\end{xmlcode}
+\end{xmllines}
  - use these 2 new variables when defining your file.  
 \vspace{-20pt}
-\begin{xmlcode}
+\begin{xmllines}
 <file\_group id="1m" output\_freq="1m"  output\_level="10" enabled=".TRUE." >  <!-- 1m files -->  
  <file id="file1" name\_suffix="\_grid\_T" description="ocean T grid variables" >
@@ -712,5 +712,5 @@
  </file>
 </file\_group> 
-\end{xmlcode}
+\end{xmllines}
 The freq\_op="1d" attribute is used to define the operation frequency of the ``@'' function: here 1 day. The temporal operation done by the ``@'' is the one defined in the field definition: here maximum for sstmax and minimum for sstmin. So, in the above case, @sstmax will do the daily max and @sstmin the daily min. Operation="average" refers to the temporal operation to be performed on the field ``@sstmax - @sstmin'': here monthly mean (of daily max - daily min of the sst). field\_ref="sst" means that attributes not explicitely defined, are inherited from sst field.
 
@@ -1024,5 +1024,5 @@
 Output from the XIOS-1.0 IO server is compliant with \href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of the CF metadata standard. Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of section \ref{IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard.
 
-Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by the namelist parameter \np{ln\_cfmeta} in the \ngn{namrun} namelist. This must be set to true if these metadata are to be included in the output files.
+Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by the namelist parameter \np{ln_cfmeta} in the \ngn{namrun} namelist. This must be set to true if these metadata are to be included in the output files.
 
 
@@ -1048,5 +1048,5 @@
 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} 
+files by setting the \np{ln_nc4zip} logical to false in the \textit{namnc4} 
 namelist:
 
@@ -1056,5 +1056,5 @@
 
 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
+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.
@@ -1106,29 +1106,29 @@
          &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\%\\
+\ifile{ORCA2\_restart\_0000} & 16420 & 8860 & 47\%\\
+\ifile{ORCA2\_restart\_0001} & 16064 & 11456 & 29\%\\
+\ifile{ORCA2\_restart\_0002} & 16064 & 9744 & 40\%\\
+\ifile{ORCA2\_restart\_0003} & 16420 & 9404 & 43\%\\
+\ifile{ORCA2\_restart\_0004} & 16200 & 5844 & 64\%\\
+\ifile{ORCA2\_restart\_0005} & 15848 & 8172 & 49\%\\
+\ifile{ORCA2\_restart\_0006} & 15848 & 8012 & 50\%\\
+\ifile{ORCA2\_restart\_0007} & 16200 & 5148 & 69\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0000} & 2200 & 1504 & 32\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0001} & 2200 & 1748 & 21\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0002} & 2200 & 1592 & 28\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0003} & 2200 & 1540 & 30\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0004} & 2200 & 1204 & 46\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0005} & 2200 & 1444 & 35\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0006} & 2200 & 1428 & 36\%\\
+\ifile{ORCA2\_2d\_grid\_T\_0007} & 2200 & 1148 & 48\%\\
+            ...         & ...  &  ... & ...  \\
+\ifile{ORCA2\_2d\_grid\_W\_0000} & 4416 & 2240 & 50\%\\
+\ifile{ORCA2\_2d\_grid\_W\_0001} & 4416 & 2924 & 34\%\\
+\ifile{ORCA2\_2d\_grid\_W\_0002} & 4416 & 2512 & 44\%\\
+\ifile{ORCA2\_2d\_grid\_W\_0003} & 4416 & 2368 & 47\%\\
+\ifile{ORCA2\_2d\_grid\_W\_0004} & 4416 & 1432 & 68\%\\
+\ifile{ORCA2\_2d\_grid\_W\_0005} & 4416 & 1972 & 56\%\\
+\ifile{ORCA2\_2d\_grid\_W\_0006} & 4416 & 2028 & 55\%\\
+\ifile{ORCA2\_2d\_grid\_W\_0007} & 4416 & 1368 & 70\%\\
 \end{tabular}
 \caption{ 	\protect\label{Tab_NC4} 
@@ -1138,5 +1138,5 @@
 
 When \key{iomput} is activated with \key{netcdf4} chunking and
-compression parameters for fields produced via \np{iom\_put} calls are
+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
@@ -1167,15 +1167,15 @@
 What is done depends on the \ngn{namtrd} logical set to \textit{true}:
 \begin{description}
-\item[\np{ln\_glo\_trd}] : at each \np{nn\_trd} time-step a check of the basin averaged properties 
+\item[\np{ln_glo_trd}] : at each \np{nn_trd} time-step a check of the basin averaged properties 
 of the momentum and tracer equations is performed. This also includes a check of $T^2$, $S^2$, 
 $\tfrac{1}{2} (u^2+v2)$, and potential energy time evolution equations properties ; 
-\item[\np{ln\_dyn\_trd}] : each 3D trend of the evolution of the two momentum components is output ; 
-\item[\np{ln\_dyn\_mxl}] : each 3D trend of the evolution of the two momentum components averaged 
+\item[\np{ln_dyn_trd}] : each 3D trend of the evolution of the two momentum components is output ; 
+\item[\np{ln_dyn_mxl}] : each 3D trend of the evolution of the two momentum components averaged 
                            over the mixed layer is output  ; 
-\item[\np{ln\_vor\_trd}] : a vertical summation of the moment tendencies is performed, 
+\item[\np{ln_vor_trd}] : a vertical summation of the moment tendencies is performed, 
                            then the curl is computed to obtain the barotropic vorticity tendencies which are output ;
-\item[\np{ln\_KE\_trd}]  : each 3D trend of the Kinetic Energy equation is output ;
-\item[\np{ln\_tra\_trd}] : each 3D trend of the evolution of temperature and salinity is output ;
-\item[\np{ln\_tra\_mxl}] : each 2D trend of the evolution of temperature and salinity averaged 
+\item[\np{ln_KE_trd}]  : each 3D trend of the Kinetic Energy equation is output ;
+\item[\np{ln_tra_trd}] : each 3D trend of the evolution of temperature and salinity is output ;
+\item[\np{ln_tra_mxl}] : each 2D trend of the evolution of temperature and salinity averaged 
                            over the mixed layer is output ;
 \end{description}
@@ -1185,5 +1185,5 @@
 
 \textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested. 
-In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} 
+In particular, options associated with \np{ln_dyn_mxl}, \np{ln_vor_trd}, and \np{ln_tra_mxl} 
 are not working, and none of the option have been tested with variable volume ($i.e.$ \key{vvl} defined).
 
@@ -1203,5 +1203,5 @@
 namelis variables. 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} 
+Runge-Hutta algorithm (\forcode{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. 
@@ -1209,9 +1209,9 @@
 \subsubsection{ Input data: initial coordinates }
 
-Initial coordinates can be given with Ariane Tools convention ( IJK coordinates ,(\np{ln\_ariane}=true) )
+Initial coordinates can be given with Ariane Tools convention ( IJK coordinates ,(\forcode{ln_ariane = .true.}) )
 or with longitude and latitude.
 
 
-In case of Ariane convention, input filename is \np{init\_float\_ariane}. Its format is:
+In case of Ariane convention, input filename is \np{init_float_ariane}. Its format is:
 
 \texttt{ I J K nisobfl itrash itrash }
@@ -1258,16 +1258,16 @@
 
 \np{jpnfl} is the total number of floats during the run.
-When initial positions are read in a restart file ( \np{ln\_rstflo}= .TRUE. ),  \np{jpnflnewflo}
+When initial positions are read in a restart file ( \np{ln_rstflo}= .TRUE. ),  \np{jpnflnewflo}
 can be added in the initialization file. 
 
 \subsubsection{ Output data }
 
-\np{nn\_writefl} is the frequency of writing in float output file and \np{nn\_stockfl} 
+\np{nn_writefl} is the frequency of writing in float output file and \np{nn_stockfl} 
 is the frequency of creation of the float restart file.
 
-Output data can be written in ascii files (\np{ln\_flo\_ascii} = .TRUE. ). In that case, 
+Output data can be written in ascii files (\np{ln_flo_ascii} = .TRUE. ). In that case, 
 output filename is trajec\_float.
 
-Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii} = .FALSE. ). There are 2 possibilities:
+Another possiblity of writing format is Netcdf (\np{ln_flo_ascii} = .FALSE. ). There are 2 possibilities:
 
  - if (\key{iomput}) is used, outputs are selected in  iodef.xml. Here it is an example of specification 
@@ -1275,5 +1275,5 @@
 
 \vspace{-30pt}
-\begin{xmlcode}
+\begin{xmllines}
      <group id="1d\_grid\_T" name="auto" description="ocean T grid variables" >   }
        <file id="floats"  description="floats variables"> }\\
@@ -1287,8 +1287,8 @@
        </file>}
   </group>}
-\end{xmlcode}
-
-
- -  if (\key{iomput}) is not used, a file called trajec\_float.nc will be created by IOIPSL library.
+\end{xmllines}
+
+
+ -  if (\key{iomput}) is not used, a file called \ifile{trajec\_float} will be created by IOIPSL library.
 
 
@@ -1312,15 +1312,15 @@
 Some parameters are available in namelist \ngn{namdia\_harm} :
 
-- \np{nit000\_han} is the first time step used for harmonic analysis
-
-- \np{nitend\_han} is the last time step used for harmonic analysis
-
-- \np{nstep\_han} is the time step frequency for harmonic analysis
-
-- \np{nb\_ana} is the number of harmonics to analyse
+- \np{nit000_han} is the first time step used for harmonic analysis
+
+- \np{nitend_han} is the last time step used for harmonic analysis
+
+- \np{nstep_han} is the time step frequency for harmonic analysis
+
+- \np{nb_ana} is the number of harmonics to analyse
 
 - \np{tname} is an array with names of tidal constituents to analyse
 
-\np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation.
+\np{nit000_han} and \np{nitend_han} must be between \np{nit000} and \np{nitend} of the simulation.
 The restart capability is not implemented.
 
@@ -1369,9 +1369,9 @@
 and the time scales over which they are averaged, as well as the level of output for debugging:
 
-\np{nn\_dct}: frequency of instantaneous transports computing
-
-\np{nn\_dctwri}: frequency of writing ( mean of instantaneous transports )
-
-\np{nn\_debug}: debugging of the section
+\np{nn_dct}: frequency of instantaneous transports computing
+
+\np{nn_dctwri}: frequency of writing ( mean of instantaneous transports )
+
+\np{nn_debug}: debugging of the section
 
 \subsubsection{ Creating a binary file containing the pathway of each section }
@@ -1681,6 +1681,6 @@
 The poleward heat and salt transports, their advective and diffusive component, and 
 the meriodional stream function can be computed on-line in \mdl{diaptr} 
-\np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below).  
-When \np{ln\_subbas}~=~true, transports and stream function are computed 
+\np{ln_diaptr} to true (see the \textit{\ngn{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 
@@ -1756,5 +1756,5 @@
 in the zonal, meridional and vertical directions respectively. The vertical component is included although it is not strictly valid as the vertical velocity is calculated from the continuity equation rather than as a prognostic variable. Physically this represents the rate at which information is propogated across a grid cell. Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step.
 
-The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of where the maximum value occurs. At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. 
+The variables can be activated by setting the \np{nn_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of where the maximum value occurs. At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. 
 
 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIU.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIU.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIU.tex	(revision 9392)
@@ -37,8 +37,8 @@
 This namelist contains only two variables:
 \begin{description}
-\item[\np{ln\_diurnal}] A logical switch for turning on/off both the cool skin and warm layer.
-\item[\np{ln\_diurnal\_only}] A logical switch which if .TRUE. will run the diurnal model
+\item[\np{ln_diurnal}] A logical switch for turning on/off both the cool skin and warm layer.
+\item[\np{ln_diurnal_only}] A logical switch which if .TRUE. will run the diurnal model
 without the other dynamical parts of NEMO.  
-\np{ln\_diurnal\_only} must be .FALSE. if \np{ln\_diurnal} is .FALSE.
+\np{ln_diurnal_only} must be .FALSE. if \np{ln_diurnal} is .FALSE.
 \end{description}
 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex	(revision 9392)
@@ -454,5 +454,5 @@
 (d) hybrid $s-z$ coordinate, 
 (e) hybrid $s-z$ coordinate with partial step, and 
-(f) same as (e) but in the non-linear free surface (\protect\np{ln_linssh}=false). 
+(f) same as (e) but in the non-linear free surface (\protect\forcode{ln_linssh = .false.}). 
 Note that the non-linear free surface can be used with any of the 
 5 coordinates (a) to (e).}
@@ -469,5 +469,5 @@
 (Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). By default a non-linear free surface is used:
 the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 
-$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln_linssh}=true), 
+$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). When a linear free surface is assumed (\forcode{ln_linssh = .true.}), 
 the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface 
 (in other words, the top of the ocean in not a rigid-lid). 
@@ -503,10 +503,10 @@
 %%%
 
-Unless a linear free surface is used (\np{ln_linssh}=false), the arrays describing 
+Unless a linear free surface is used (\forcode{ln_linssh = .false.}), the arrays describing 
 the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ 
 defined at \textit{before}, \textit{now} and \textit{after} time step. The time at which they are
 defined is indicated by a suffix:$\_b$, $\_n$, or $\_a$, respectively. They are updated at each model time step
 using a fixed reference coordinate system which computer names have a $\_0$ suffix. 
-When the linear free surface option is used (\np{ln_linssh}=true), \textit{before}, \textit{now} 
+When the linear free surface option is used (\forcode{ln_linssh = .true.}), \textit{before}, \textit{now} 
 and \textit{after} arrays are simply set one for all to their reference counterpart. 
 
@@ -551,5 +551,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection[$z$-coordinate (\protect\np{ln_zco}]
-	  	  {$z$-coordinate (\protect\np{ln_zco}=true) and reference coordinate}
+	  	  {$z$-coordinate (\protect\forcode{ln_zco = .true.}) and reference coordinate}
 \label{DOM_zco}
 
@@ -629,5 +629,5 @@
 Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is 
 possible to recalculate them. In that case the user sets 
-\np{ppsur}=\np{ppa0}=\np{ppa1}=999999., in \ngn{namcfg} namelist, 
+\np{ppsur}=\np{ppa0}=\forcode{ppa1 = 999999}., in \ngn{namcfg} namelist, 
 and specifies instead the four following parameters:
 \begin{itemize}
@@ -640,6 +640,6 @@
 \end{itemize}
 As an example, for the $45$ layers used in the DRAKKAR configuration those 
-parameters are: \jp{jpk}=46, \np{ppacr}=9, \np{ppkth}=23.563, \np{ppdzmin}=6m, 
-\np{pphmax}=5750m.
+parameters are: \jp{jpk}=46, \forcode{ppacr = 9}, \forcode{ppkth = 23}.563, \forcode{ppdzmin = 6}m, 
+\forcode{pphmax = 5750}m.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -722,5 +722,5 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [$s$-coordinate (\protect\np{ln_sco})]
-		     {$s$-coordinate (\protect\np{ln_sco}=true)}
+		     {$s$-coordinate (\protect\forcode{ln_sco = .true.})}
 \label{DOM_sco}
 %------------------------------------------nam_zgr_sco---------------------------------------------------
@@ -850,5 +850,5 @@
 %        z*- or s*-coordinate
 % -------------------------------------------------------------------------------------------------------------
-\subsection{$z^*$- or $s^*$-coordinate (\protect\np{ln_linssh}=false) }
+\subsection{$z^*$- or $s^*$-coordinate (\protect\forcode{ln_linssh = .false.}) }
 \label{DOM_zgr_star}
 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex	(revision 9392)
@@ -200,5 +200,5 @@
 %                 enstrophy conserving scheme
 %-------------------------------------------------------------
-\subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}=true)}
+\subsubsection{Enstrophy conserving scheme (\protect\forcode{ln_dynvor_ens = .true.})}
 \label{DYN_vor_ens}
 
@@ -221,5 +221,5 @@
 %                 energy conserving scheme
 %-------------------------------------------------------------
-\subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}=true)}
+\subsubsection{Energy conserving scheme (\protect\forcode{ln_dynvor_ene = .true.})}
 \label{DYN_vor_ene}
 
@@ -238,5 +238,5 @@
 %                 mix energy/enstrophy conserving scheme
 %-------------------------------------------------------------
-\subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}=true) }
+\subsubsection{Mixed energy/enstrophy conserving scheme (\protect\forcode{ln_dynvor_mix = .true.}) }
 \label{DYN_vor_mix}
 
@@ -261,5 +261,5 @@
 %                 energy and enstrophy conserving scheme
 %-------------------------------------------------------------
-\subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}=true) }
+\subsubsection{Energy and enstrophy conserving scheme (\protect\forcode{ln_dynvor_een = .true.}) }
 \label{DYN_vor_een}
 
@@ -305,6 +305,6 @@
 A key point in \eqref{Eq_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. 
 It uses the sum of masked t-point vertical scale factor divided either 
-by the sum of the four t-point masks (\np{nn\_een\_e3f}~=~1), 
-or  just by $4$ (\np{nn\_een\_e3f}~=~true).
+by the sum of the four t-point masks (\np{nn_een_e3f}~=~1), 
+or  just by $4$ (\np{nn_een_e3f}~=~true).
 The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ 
 tends to zero and extends by continuity the value of $e_{3f}$ into the land areas. 
@@ -377,10 +377,10 @@
 \end{aligned}         \right.
 \end{equation} 
-When \np{ln\_dynzad\_zts}~=~\textit{true}, a split-explicit time stepping with 5 sub-timesteps is used 
+When \np{ln_dynzad_zts}~=~\textit{true}, a split-explicit time stepping with 5 sub-timesteps is used 
 on the vertical advection term.
 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}. 
 Note that in this case, a similar split-explicit time stepping should be used on 
 vertical advection of tracer to ensure a better stability, 
-an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \S\ref{TRA_adv_tvd}).
+an option which is only available with a TVD scheme (see \np{ln_traadv_tvd_zts} in \S\ref{TRA_adv_tvd}).
 
 
@@ -451,5 +451,5 @@
 difference scheme, CEN2, or a $3^{rd}$ order upstream biased scheme, UBS. 
 The latter is described in \citet{Shchepetkin_McWilliams_OM05}. The schemes are 
-selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. 
+selected using the namelist logicals \np{ln_dynadv_cen2} and \np{ln_dynadv_ubs}. 
 In flux form, the schemes differ by the choice of a space and time interpolation to 
 define the value of $u$ and $v$ at the centre of each face of $u$- and $v$-cells, 
@@ -460,5 +460,5 @@
 %                 2nd order centred scheme
 %-------------------------------------------------------------
-\subsubsection{$2^{nd}$ order centred scheme (cen2) (\protect\np{ln\_dynadv\_cen2}=true)}
+\subsubsection{$2^{nd}$ order centred scheme (cen2) (\protect\forcode{ln_dynadv_cen2 = .true.})}
 \label{DYN_adv_cen2}
 
@@ -481,5 +481,5 @@
 %                 UBS scheme
 %-------------------------------------------------------------
-\subsubsection{Upstream Biased Scheme (UBS) (\protect\np{ln\_dynadv\_ubs}=true)}
+\subsubsection{Upstream Biased Scheme (UBS) (\protect\forcode{ln_dynadv_ubs = .true.})}
 \label{DYN_adv_ubs}
 
@@ -501,5 +501,5 @@
 those in the centred second order method. As the scheme already includes 
 a diffusion component, it can be used without explicit  lateral diffusion on momentum 
-($i.e.$ \np{ln\_dynldf\_lap}=\np{ln\_dynldf\_bilap}=false), and it is recommended to do so.
+($i.e.$ \np{ln_dynldf_lap}=\forcode{ln_dynldf_bilap = .false.}), and it is recommended to do so.
 
 The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$ 
@@ -554,6 +554,6 @@
 %           z-coordinate with full step
 %--------------------------------------------------------------------------------------------------------------
-\subsection   [$z$-coordinate with full step (\protect\np{ln\_dynhpg\_zco}) ]
-			{$z$-coordinate with full step (\protect\np{ln\_dynhpg\_zco}=true)}
+\subsection   [$z$-coordinate with full step (\protect\np{ln_dynhpg_zco}) ]
+			{$z$-coordinate with full step (\protect\forcode{ln_dynhpg_zco = .true.})}
 \label{DYN_hpg_zco}
 
@@ -595,6 +595,6 @@
 %           z-coordinate with partial step
 %--------------------------------------------------------------------------------------------------------------
-\subsection   [$z$-coordinate with partial step (\protect\np{ln\_dynhpg\_zps})]
-			{$z$-coordinate with partial step (\protect\np{ln\_dynhpg\_zps}=true)}
+\subsection   [$z$-coordinate with partial step (\protect\np{ln_dynhpg_zps})]
+			{$z$-coordinate with partial step (\protect\forcode{ln_dynhpg_zps = .true.})}
 \label{DYN_hpg_zps}
 
@@ -624,5 +624,5 @@
 cubic polynomial method is currently disabled whilst known bugs are under investigation.
 
-$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=true)
+$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\forcode{ln_dynhpg_sco = .true.})
 \begin{equation} \label{Eq_dynhpg_sco}
 \left\{ \begin{aligned}
@@ -639,16 +639,16 @@
 ($e_{3w}$).
  
-$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}=true).
-This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}=true).
-
-$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true)
+$\bullet$ Traditional coding with adaptation for ice shelf cavities (\forcode{ln_dynhpg_isf = .true.}).
+This scheme need the activation of ice shelf cavities (\forcode{ln_isfcav = .true.}).
+
+$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\forcode{ln_dynhpg_prj = .true.})
 
 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05} 
-(\np{ln\_dynhpg\_djc}=true) (currently disabled; under development)
+(\forcode{ln_dynhpg_djc = .true.}) (currently disabled; under development)
 
 Note that expression \eqref{Eq_dynhpg_sco} is commonly used when the variable volume formulation is
 activated (\key{vvl}) because in that case, even with a flat bottom, the coordinate surfaces are not
 horizontal but follow the free surface \citep{Levier2007}. The pressure jacobian scheme
-(\np{ln\_dynhpg\_prj}=true) is available as an improved option to \np{ln\_dynhpg\_sco}=true when
+(\forcode{ln_dynhpg_prj = .true.}) is available as an improved option to \forcode{ln_dynhpg_sco = .true.} when
 \key{vvl} is active.  The pressure Jacobian scheme uses a constrained cubic spline to reconstruct
 the density profile across the water column. This method maintains the monotonicity between the
@@ -660,6 +660,6 @@
 \label{DYN_hpg_isf}
 Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and
- the pressure gradient due to the ocean load. If cavity opened (\np{ln\_isfcav}~=~true) these 2 terms can be
- calculated by setting \np{ln\_dynhpg\_isf}~=~true. No other scheme are working with the ice shelf.\\
+ the pressure gradient due to the ocean load. If cavity opened (\np{ln_isfcav}~=~true) these 2 terms can be
+ calculated by setting \np{ln_dynhpg_isf}~=~true. No other scheme are working with the ice shelf.\\
 
 $\bullet$ The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium.
@@ -673,6 +673,6 @@
 %           Time-scheme
 %--------------------------------------------------------------------------------------------------------------
-\subsection   [Time-scheme (\protect\np{ln\_dynhpg\_imp}) ]
-			{Time-scheme (\protect\np{ln\_dynhpg\_imp}= true/false)}
+\subsection   [Time-scheme (\protect\np{ln_dynhpg_imp}) ]
+			{Time-scheme (\protect\np{ln_dynhpg_imp}= true/false)}
 \label{DYN_hpg_imp}
 
@@ -689,5 +689,5 @@
 time level $t$ only, as in the standard leapfrog scheme. 
 
-$\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}=true):
+$\bullet$ leapfrog scheme (\forcode{ln_dynhpg_imp = .true.}):
 
 \begin{equation} \label{Eq_dynhpg_lf}
@@ -696,5 +696,5 @@
 \end{equation}
 
-$\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}=true):
+$\bullet$ semi-implicit scheme (\forcode{ln_dynhpg_imp = .true.}):
 \begin{equation} \label{Eq_dynhpg_imp}
 \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
@@ -713,5 +713,5 @@
 the stability limits associated with advection or diffusion.
 
-In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}=true. 
+In practice, the semi-implicit scheme is used when \forcode{ln_dynhpg_imp = .true.}. 
 In this case, we choose to apply the time filter to temperature and salinity used in 
 the equation of state, instead of applying it to the hydrostatic pressure or to the 
@@ -727,5 +727,5 @@
 Note that in the semi-implicit case, it is necessary to save the filtered density, an 
 extra three-dimensional field, in the restart file to restart the model with exact 
-reproducibility. This option is controlled by  \np{nn\_dynhpg\_rst}, a namelist parameter.
+reproducibility. This option is controlled by  \np{nn_dynhpg_rst}, a namelist parameter.
 
 % ================================================================
@@ -806,9 +806,9 @@
 variables (Fig.~\ref{Fig_DYN_dynspg_ts}). 
 The size of the small time step, $\rdt_e$ (the external mode or barotropic time step)
- is provided through the \np{nn\_baro} namelist parameter as: 
-$\rdt_e = \rdt / nn\_baro$. This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}=true) 
+ is provided through the \np{nn_baro} namelist parameter as: 
+$\rdt_e = \rdt / nn\_baro$. This parameter can be optionally defined automatically (\forcode{ln_bt_nn_auto = .true.}) 
 considering that the stability of the barotropic system is essentially controled by external waves propagation. 
 Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry.
-Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn\_bt\_cmax}.
+Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn_bt_cmax}.
 
 %%%
@@ -839,18 +839,18 @@
 The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged 
 transports to advect tracers.
-a) Forward time integration: \protect\np{ln\_bt\_fw}=true,  \protect\np{ln\_bt\_av}=true. 
-b) Centred time integration: \protect\np{ln\_bt\_fw}=false, \protect\np{ln\_bt\_av}=true. 
-c) Forward time integration with no time filtering (POM-like scheme): \protect\np{ln\_bt\_fw}=true, \protect\np{ln\_bt\_av}=false. }
+a) Forward time integration: \protect\forcode{ln_bt_fw = .true.},  \protect\forcode{ln_bt_av = .true.}. 
+b) Centred time integration: \protect\forcode{ln_bt_fw = .false.}, \protect\forcode{ln_bt_av = .true.}. 
+c) Forward time integration with no time filtering (POM-like scheme): \protect\forcode{ln_bt_fw = .true.}, \protect\forcode{ln_bt_av = .false.}. }
 \end{center}    \end{figure}
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
 
-In the default case (\np{ln\_bt\_fw}=true), the external mode is integrated 
+In the default case (\forcode{ln_bt_fw = .true.}), the external mode is integrated 
 between \textit{now} and  \textit{after} baroclinic time-steps (Fig.~\ref{Fig_DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic 
-quantities (\np{ln\_bt\_av}=true). In that case, the integration is extended slightly beyond  \textit{after} time step to provide time filtered quantities. 
+quantities (\forcode{ln_bt_av = .true.}). In that case, the integration is extended slightly beyond  \textit{after} time step to provide time filtered quantities. 
 These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. 
 Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, 
 asselin filtering is not applied to barotropic quantities. \\
 Alternatively, one can choose to integrate barotropic equations starting 
-from \textit{before} time step (\np{ln\_bt\_fw}=false). Although more computationaly expensive ( \np{nn\_baro} additional iterations are indeed necessary), the baroclinic to barotropic forcing term given at \textit{now} time step 
+from \textit{before} time step (\forcode{ln_bt_fw = .false.}). Although more computationaly expensive ( \np{nn_baro} additional iterations are indeed necessary), the baroclinic to barotropic forcing term given at \textit{now} time step 
 become centred in the middle of the integration window. It can easily be shown that this property 
 removes part of splitting errors between modes, which increases the overall numerical robustness.
@@ -868,5 +868,5 @@
 %%%
 
-One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_av}=false). 
+One can eventually choose to feedback instantaneous values by not using any time filter (\forcode{ln_bt_av = .false.}). 
 In that case, external mode equations are continuous in time, ie they are not re-initialized when starting a new 
 sub-stepping sequence. This is the method used so far in the POM model, the stability being maintained by refreshing at (almost) 
@@ -1036,6 +1036,6 @@
 
 % ================================================================
-\subsection   [Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}) ]
-			{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}=true)}
+\subsection   [Iso-level laplacian operator (\protect\np{ln_dynldf_lap}) ]
+			{Iso-level laplacian operator (\protect\forcode{ln_dynldf_lap = .true.})}
 \label{DYN_ldf_lap}
 
@@ -1060,14 +1060,14 @@
 %           Rotated laplacian operator
 %--------------------------------------------------------------------------------------------------------------
-\subsection   [Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}) ]
-			{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}=true)}
+\subsection   [Rotated laplacian operator (\protect\np{ln_dynldf_iso}) ]
+			{Rotated laplacian operator (\protect\forcode{ln_dynldf_iso = .true.})}
 \label{DYN_ldf_iso}
 
 A rotation of the lateral momentum diffusion operator is needed in several cases: 
-for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}=true) and for 
-either iso-neutral (\np{ln\_dynldf\_iso}=true) or geopotential 
-(\np{ln\_dynldf\_hor}=true) diffusion in the $s$-coordinate. In the partial step 
+for iso-neutral diffusion in the $z$-coordinate (\forcode{ln_dynldf_iso = .true.}) and for 
+either iso-neutral (\forcode{ln_dynldf_iso = .true.}) or geopotential 
+(\forcode{ln_dynldf_hor = .true.}) diffusion in the $s$-coordinate. In the partial step 
 case, coordinates are horizontal except at the deepest level and no 
-rotation is performed when \np{ln\_dynldf\_hor}=true. The diffusion operator 
+rotation is performed when \forcode{ln_dynldf_hor = .true.}. The diffusion operator 
 is defined simply as the divergence of down gradient momentum fluxes on each 
 momentum component. It must be emphasized that this formulation ignores 
@@ -1129,6 +1129,6 @@
 %           Iso-level bilaplacian operator
 %--------------------------------------------------------------------------------------------------------------
-\subsection   [Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap})]
-			{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}=true)}
+\subsection   [Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap})]
+			{Iso-level bilaplacian operator (\protect\forcode{ln_dynldf_bilap = .true.})}
 \label{DYN_ldf_bilap}
 
@@ -1157,8 +1157,8 @@
 would be too restrictive a constraint on the time step. Two time stepping schemes 
 can be used for the vertical diffusion term : $(a)$ a forward time differencing 
-scheme (\np{ln\_zdfexp}=true) using a time splitting technique 
-(\np{nn\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme 
-(\np{ln\_zdfexp}=false) (see \S\ref{STP}). Note that namelist variables 
-\np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 
+scheme (\forcode{ln_zdfexp = .true.}) using a time splitting technique 
+(\np{nn_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme 
+(\forcode{ln_zdfexp = .false.}) (see \S\ref{STP}). Note that namelist variables 
+\np{ln_zdfexp} and \np{nn_zdfexp} apply to both tracers and dynamics. 
 
 The formulation of the vertical subgrid scale physics is the same whatever 
@@ -1206,11 +1206,11 @@
 may enter the dynamical equations by affecting the surface pressure gradient. 
 
-(1) When \np{ln\_apr\_dyn}~=~true (see \S\ref{SBC_apr}), the atmospheric pressure is taken 
+(1) When \np{ln_apr_dyn}~=~true (see \S\ref{SBC_apr}), the atmospheric pressure is taken 
 into account when computing the surface pressure gradient.
 
-(2) When \np{ln\_tide\_pot}~=~true and \np{ln\_tide}~=~true (see \S\ref{SBC_tide}), 
+(2) When \np{ln_tide_pot}~=~true and \np{ln_tide}~=~true (see \S\ref{SBC_tide}), 
 the tidal potential is taken into account when computing the surface pressure gradient.
 
-(3) When \np{nn\_ice\_embd}~=~2 and LIM or CICE is used ($i.e.$ when the sea-ice is embedded in the ocean), 
+(3) When \np{nn_ice_embd}~=~2 and LIM or CICE is used ($i.e.$ when the sea-ice is embedded in the ocean), 
 the snow-ice mass is taken into account when computing the surface pressure gradient.
 
@@ -1238,5 +1238,5 @@
 weighted velocity (see \S\ref{Apdx_A_momentum})  
 
-$\bullet$ vector invariant form or linear free surface (\np{ln\_dynhpg\_vec}=true ; \key{vvl} not defined):
+$\bullet$ vector invariant form or linear free surface (\forcode{ln_dynhpg_vec = .true.} ; \key{vvl} not defined):
 \begin{equation} \label{Eq_dynnxt_vec}
 \left\{   \begin{aligned}
@@ -1246,5 +1246,5 @@
 \end{equation} 
 
-$\bullet$ flux form and nonlinear free surface (\np{ln\_dynhpg\_vec}=false ; \key{vvl} defined):
+$\bullet$ flux form and nonlinear free surface (\forcode{ln_dynhpg_vec = .false.} ; \key{vvl} defined):
 \begin{equation} \label{Eq_dynnxt_flux}
 \left\{   \begin{aligned}
@@ -1256,5 +1256,5 @@
 where RHS is the right hand side of the momentum equation, the subscript $f$ 
 denotes filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is 
-initialized as \np{nn\_atfp} (namelist parameter). Its default value is \np{nn\_atfp} = $10^{-3}$.
+initialized as \np{nn_atfp} (namelist parameter). Its default value is \np{nn_atfp} = $10^{-3}$.
 In both cases, the modified Asselin filter is not applied since perfect conservation 
 is not an issue for the momentum equations.
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_LBC.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_LBC.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_LBC.tex	(revision 9392)
@@ -17,5 +17,5 @@
 % Boundary Condition at the Coast
 % ================================================================
-\section{Boundary Condition at the Coast (\protect\np{rn\_shlat})}
+\section{Boundary Condition at the Coast (\protect\np{rn_shlat})}
 \label{LBC_coast}
 %--------------------------------------------nam_lbc-------------------------------------------------------
@@ -72,5 +72,5 @@
 condition influences the relative vorticity and momentum diffusive trends, and is 
 required in order to compute the vorticity at the coast. Four different types of 
-lateral boundary condition are available, controlled by the value of the \np{rn\_shlat} 
+lateral boundary condition are available, controlled by the value of the \np{rn_shlat} 
 namelist parameter. (The value of the mask$_{f}$ array along the coastline is set 
 equal to this parameter.) These are:
@@ -88,10 +88,10 @@
 \begin{description}
 
-\item[free-slip boundary condition (\np{rn\_shlat}=0): ]  the tangential velocity at the 
+\item[free-slip boundary condition (\forcode{rn_shlat = 0}): ]  the tangential velocity at the 
 coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the 
 tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set 
 to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a).
 
-\item[no-slip boundary condition (\np{rn\_shlat}=2): ] the tangential velocity vanishes 
+\item[no-slip boundary condition (\forcode{rn_shlat = 2}): ] the tangential velocity vanishes 
 at the coastline. Assuming that the tangential velocity decreases linearly from 
 the closest ocean velocity grid point to the coastline, the normal derivative is 
@@ -112,5 +112,5 @@
 \end{equation}
 
-\item["partial" free-slip boundary condition (0$<$\np{rn\_shlat}$<$2): ] the tangential 
+\item["partial" free-slip boundary condition (0$<$\np{rn_shlat}$<$2): ] the tangential 
 velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral 
 friction but not strong enough to make the tangential velocity at the coast vanish 
@@ -118,5 +118,5 @@
 strictly inbetween $0$ and $2$.
 
-\item["strong" no-slip boundary condition (2$<$\np{rn\_shlat}): ] the viscous boundary 
+\item["strong" no-slip boundary condition (2$<$\np{rn_shlat}): ] the viscous boundary 
 layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d). 
 The friction is thus larger than in the no-slip case.
@@ -331,5 +331,5 @@
 the model output files is undefined. Note that this is a problem for the meshmask file 
 which requires to be defined over the whole domain. Therefore, user should not eliminate 
-land processors when creating a meshmask file ($i.e.$ when setting a non-zero value to \np{nn\_msh}).
+land processors when creating a meshmask file ($i.e.$ when setting a non-zero value to \np{nn_msh}).
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -387,15 +387,15 @@
 \label{BDY_namelist}
 
-The BDY module is activated by setting \np{ln\_bdy} to true.
+The BDY module is activated by setting \np{ln_bdy} to true.
 It is possible to define more than one boundary ``set'' and apply
 different boundary conditions to each set. The number of boundary
-sets is defined by \np{nb\_bdy}.  Each boundary set may be defined
+sets is defined by \np{nb_bdy}.  Each boundary set may be defined
 as a set of straight line segments in a namelist
-(\np{ln\_coords\_file}=.false.) or read in from a file
-(\np{ln\_coords\_file}=.true.). If the set is defined in a namelist,
+(\np{ln_coords_file}=.false.) or read in from a file
+(\np{ln_coords_file}=.true.). If the set is defined in a namelist,
 then the namelists nambdy\_index must be included separately, one for
 each set. If the set is defined by a file, then a
-``coordinates.bdy.nc'' file must be provided. The coordinates.bdy file
-is analagous to the usual NEMO ``coordinates.nc'' file. In the example
+``\ifile{coordinates.bdy}'' file must be provided. The coordinates.bdy file
+is analagous to the usual NEMO ``\ifile{coordinates}'' file. In the example
 above, there are two boundary sets, the first of which is defined via
 a file and the second is defined in a namelist. For more details of
@@ -410,6 +410,6 @@
 (``tra''). For each set of variables there is a choice of algorithm
 and a choice for the data, eg. for the active tracers the algorithm is
-set by \np{nn\_tra} and the choice of data is set by
-\np{nn\_tra\_dta}. 
+set by \np{nn_tra} and the choice of data is set by
+\np{nn_tra_dta}. 
 
 The choice of algorithm is currently as follows:
@@ -429,6 +429,6 @@
 
 The main choice for the boundary data is
-to use initial conditions as boundary data (\np{nn\_tra\_dta}=0) or to
-use external data from a file (\np{nn\_tra\_dta}=1). For the
+to use initial conditions as boundary data (\forcode{nn_tra_dta = 0}) or to
+use external data from a file (\forcode{nn_tra_dta = 1}). For the
 barotropic solution there is also the option to use tidal
 harmonic forcing either by itself or in addition to other external
@@ -492,5 +492,5 @@
 \end{equation}
 The width of the FRS zone is specified in the namelist as 
-\np{nn\_rimwidth}. This is typically set to a value between 8 and 10. 
+\np{nn_rimwidth}. This is typically set to a value between 8 and 10. 
 
 %----------------------------------------------
@@ -534,5 +534,5 @@
 
 The boundary geometry for each set may be defined in a namelist
-nambdy\_index or by reading in a ``coordinates.bdy.nc'' file. The
+nambdy\_index or by reading in a ``\ifile{coordinates.bdy}'' file. The
 nambdy\_index namelist defines a series of straight-line segments for
 north, east, south and west boundaries. For the northern boundary,
@@ -546,5 +546,5 @@
 
 The boundary geometry may also be defined from a
-``coordinates.bdy.nc'' file. Figure \ref{Fig_LBC_nc_header}
+``\ifile{coordinates.bdy}'' file. Figure \ref{Fig_LBC_nc_header}
 gives an example of the header information from such a file. The file
 should contain the index arrays for each of the $T$, $U$ and $V$
@@ -561,5 +561,5 @@
 shelf break, then the areas of ocean outside of this boundary will
 need to be masked out. This can be done by reading a mask file defined
-as \np{cn\_mask\_file} in the nam\_bdy namelist. Only one mask file is
+as \np{cn_mask_file} in the nam\_bdy namelist. Only one mask file is
 used even if multiple boundary sets are defined.
 
@@ -609,5 +609,5 @@
 \includegraphics[width=1.0\textwidth]{Fig_LBC_nc_header}
 \caption {     \protect\label{Fig_LBC_nc_header} 
-Example of the header for a coordinates.bdy.nc file}
+Example of the header for a \ifile{coordinates.bdy} file}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -618,9 +618,9 @@
 
 There is an option to force the total volume in the regional model to be constant, 
-similar to the option in the OBC module. This is controlled  by the \np{nn\_volctl} 
-parameter in the namelist. A value of \np{nn\_volctl}~=~0 indicates that this option is not used. 
-If  \np{nn\_volctl}~=~1 then a correction is applied to the normal velocities 
+similar to the option in the OBC module. This is controlled  by the \np{nn_volctl} 
+parameter in the namelist. A value of \np{nn_volctl}~=~0 indicates that this option is not used. 
+If  \np{nn_volctl}~=~1 then a correction is applied to the normal velocities 
 around the boundary at each timestep to ensure that the integrated volume flow 
-through the boundary is zero. If \np{nn\_volctl}~=~2 then the calculation of 
+through the boundary is zero. If \np{nn_volctl}~=~2 then the calculation of 
 the volume change on the timestep includes the change due to the freshwater 
 flux across the surface and the correction velocity corrects for this as well.
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_LDF.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_LDF.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_LDF.tex	(revision 9392)
@@ -25,5 +25,5 @@
 Note that this chapter describes the standard implementation of iso-neutral
 tracer mixing, and Griffies's implementation, which is used if
-\np{traldf\_grif}=true, is described in Appdx\ref{sec:triad}
+\forcode{traldf_grif = .true.}, is described in Appdx\ref{sec:triad}
 
 %-----------------------------------nam_traldf - nam_dynldf--------------------------------------------
@@ -46,5 +46,5 @@
 A direction for lateral mixing has to be defined when the desired operator does 
 not act along the model levels. This occurs when $(a)$ horizontal mixing is 
-required on tracer or momentum (\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor}) 
+required on tracer or momentum (\np{ln_traldf_hor} or \np{ln_dynldf_hor}) 
 in $s$- or mixed $s$-$z$- coordinates, and $(b)$ isoneutral mixing is required 
 whatever the vertical coordinate is. This direction of mixing is defined by its 
@@ -88,6 +88,6 @@
 %gm%  caution I'm not sure the simplification was a good idea! 
 
-These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}=True, 
-and either \np{ln\_traldf\_hor}=True or \np{ln\_dynldf\_hor}=True. 
+These slopes are computed once in \rou{ldfslp\_init} when \forcode{ln_sco = .true.}rue, 
+and either \forcode{ln_traldf_hor = .true.}rue or \forcode{ln_dynldf_hor = .true.}rue. 
 
 \subsection{Slopes for tracer iso-neutral mixing}\label{LDF_slp_iso}
@@ -147,7 +147,7 @@
 \item[$s$- or hybrid $s$-$z$- coordinate : ] in the current release of \NEMO, 
 iso-neutral mixing is only employed for $s$-coordinates if the
-Griffies scheme is used (\np{traldf\_grif}=true; see Appdx \ref{sec:triad}). 
+Griffies scheme is used (\forcode{traldf_grif = .true.}; see Appdx \ref{sec:triad}). 
 In other words, iso-neutral mixing will only be accurately represented with a 
-linear equation of state (\np{nn\_eos}=1 or 2). In the case of a "true" equation 
+linear equation of state (\forcode{nn_eos = 1} or 2). In the case of a "true" equation 
 of state, the evaluation of $i$ and $j$ derivatives in \eqref{Eq_ldfslp_iso} 
 will include a pressure dependent part, leading to the wrong evaluation of 
@@ -212,5 +212,5 @@
 ocean model are modified \citep{Weaver_Eby_JPO97,
   Griffies_al_JPO98}. Griffies's scheme is now available in \NEMO if
-\np{traldf\_grif\_iso} is set true; see Appdx \ref{sec:triad}. Here,
+\np{traldf_grif_iso} is set true; see Appdx \ref{sec:triad}. Here,
 another strategy is presented \citep{Lazar_PhD97}: a local
 filtering of the iso-neutral slopes (made on 9 grid-points) prevents
@@ -347,5 +347,5 @@
 When none of the \textbf{key\_dynldf\_...} and \textbf{key\_traldf\_...} keys are 
 defined, a constant value is used over the whole ocean for momentum and 
-tracers, which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist 
+tracers, which is specified through the \np{rn_ahm0} and \np{rn_aht0} namelist 
 parameters.
 
@@ -356,5 +356,5 @@
 mixing coefficients will require 3D arrays. In the 1D option, a hyperbolic variation 
 of the lateral mixing coefficient is introduced in which the surface value is 
-\np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value, 
+\np{rn_aht0} (\np{rn_ahm0}), the bottom value is 1/4 of the surface value, 
 and the transition takes place around z=300~m with a width of 300~m 
 ($i.e.$ both the depth and the width of the inflection point are set to 300~m). 
@@ -372,5 +372,5 @@
 \end{equation}
 where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked 
-ocean domain, and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer) 
+ocean domain, and $A_o^l$ is the \np{rn_ahm0} (momentum) or \np{rn_aht0} (tracer) 
 namelist parameter. This variation is intended to reflect the lesser need for subgrid 
 scale eddy mixing where the grid size is smaller in the domain. It was introduced in 
@@ -384,6 +384,6 @@
 Other formulations can be introduced by the user for a given configuration. 
 For example, in the ORCA2 global ocean model (see Configurations), the laplacian 
-viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ 
-north and south and decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s 
+viscosity operator uses \np{rn_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ 
+north and south and decreases linearly to \np{rn_aht0}~= 2.10$^3$ m$^2$/s 
 at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}. This modification 
 can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}. 
@@ -423,5 +423,5 @@
 (3) for isopycnal diffusion on momentum or tracers, an additional purely 
 horizontal background diffusion with uniform coefficient can be added by 
-setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0}, a background horizontal 
+setting a non zero value of \np{rn_ahmb0} or \np{rn_ahtb0}, a background horizontal 
 eddy viscosity or diffusivity coefficient (namelist parameters whose default 
 values are $0$). However, the technique used to compute the isopycnal 
@@ -438,7 +438,7 @@
 (6) it is possible to use both the laplacian and biharmonic operators concurrently.
 
-(7) it is possible to run without explicit lateral diffusion on momentum (\np{ln\_dynldf\_lap} = 
-\np{ln\_dynldf\_bilap} = false). This is recommended when using the UBS advection 
-scheme on momentum (\np{ln\_dynadv\_ubs} = true, see \ref{DYN_adv_ubs}) 
+(7) it is possible to run without explicit lateral diffusion on momentum (\np{ln_dynldf_lap} = 
+\np{ln_dynldf_bilap} = false). This is recommended when using the UBS advection 
+scheme on momentum (\np{ln_dynadv_ubs} = true, see \ref{DYN_adv_ubs}) 
 and can be useful for testing purposes.
 
@@ -455,6 +455,6 @@
 described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd},
 N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and
-GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and
-\np{rn\_aeiv\_0}. If 2D-varying coefficients are set with
+GM diffusivity $A_e$ are directly set by \np{rn_aeih_0} and
+\np{rn_aeiv_0}. If 2D-varying coefficients are set with
 \key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal
 scale factor according to \eqref{Eq_title} \footnote{Except in global ORCA
@@ -467,5 +467,5 @@
   case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further
   reduced by a factor $|f/f_{20}|$, where $f_{20}$ is the value of $f$
-  at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored
+  at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn_aeiv_0} is ignored
 unless it is zero.
 }
@@ -485,5 +485,5 @@
 \end{equation}
 where $A^{eiv}$ is the eddy induced velocity coefficient whose value is set 
-through \np{rn\_aeiv}, a \textit{nam\_traldf} namelist parameter. 
+through \np{rn_aeiv}, a \textit{nam\_traldf} namelist parameter. 
 The three components of the eddy induced velocity are computed and add 
 to the eulerian velocity in \mdl{traadv\_eiv}. This has been preferred to a 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_OBS.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_OBS.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_OBS.tex	(revision 9392)
@@ -24,5 +24,5 @@
 The OBS code is called from \mdl{nemogcm} for model initialisation and to calculate the model
 equivalent values for observations on the 0th timestep. The code is then called again after
-each timestep from \mdl{step}. The code is only activated if the namelist logical \np{ln\_diaobs}
+each timestep from \mdl{step}. The code is only activated if the namelist logical \np{ln_diaobs}
 is set to true.
 
@@ -34,7 +34,7 @@
 Some profile observation types (e.g. tropical moored buoys) are made available as daily averaged quantities.
 The observation operator code can be set-up to calculate the equivalent daily average model temperature fields
-using the \np{nn\_profdavtypes} namelist array. Some SST observations are equivalent to a night-time
+using the \np{nn_profdavtypes} namelist array. Some SST observations are equivalent to a night-time
 average value and the observation operator code can calculate equivalent night-time average model SST fields by
-setting the namelist value \np{ln\_sstnight} to true. Otherwise the model value from the nearest timestep to the
+setting the namelist value \np{ln_sstnight} to true. Otherwise the model value from the nearest timestep to the
 observation time is used.
 
@@ -88,14 +88,14 @@
 
 Options are defined through the  \ngn{namobs} namelist variables.
-The options \np{ln\_t3d} and \np{ln\_s3d} switch on the temperature and salinity
+The options \np{ln_t3d} and \np{ln_s3d} switch on the temperature and salinity
 profile observation operator code. The filename or array of filenames are
-specified using the \np{cn\_profbfiles} variable. The model grid points for a
+specified using the \np{cn_profbfiles} variable. The model grid points for a
 particular  observation latitude and longitude are found using the grid
 searching part of the code. This can be expensive, particularly for large
-numbers of observations, setting \np{ln\_grid\_search\_lookup} allows the use of
+numbers of observations, setting \np{ln_grid_search_lookup} allows the use of
 a lookup table which is saved into an ``xypos`` file (or files). This will need
 to be generated the first time if it does not exist in the run directory.
 However, once produced it will significantly speed up future grid searches.
-Setting \np{ln\_grid\_global} means that the code distributes the observations
+Setting \np{ln_grid_global} means that the code distributes the observations
 evenly between processors. Alternatively each processor will work with
 observations located within the model subdomain (see section~\ref{OBS_parallel}).
@@ -406,5 +406,5 @@
 The mean dynamic
 topography (MDT) must be provided in a separate file defined on the model grid
- called {\it slaReferenceLevel.nc}. The MDT is required in
+ called \ifile{slaReferenceLevel}. The MDT is required in
 order to produce the model equivalent sea level anomaly from the model sea
 surface height. Below is an example header for this file (on the ORCA025 grid).
@@ -556,17 +556,17 @@
 NEMO therefore has the capability to specify either an interpolation or an averaging (for surface observation types only). 
 
-The main namelist option associated with the interpolation/averaging is \np{nn\_2dint}. This default option can be set to values from 0 to 6. 
+The main namelist option associated with the interpolation/averaging is \np{nn_2dint}. This default option can be set to values from 0 to 6. 
 Values between 0 to 4 are associated with interpolation while values 5 or 6 are associated with averaging.
 \begin{itemize}
-\item \np{nn\_2dint}=0: Distance-weighted interpolation
-\item \np{nn\_2dint}=1: Distance-weighted interpolation (small angle)
-\item \np{nn\_2dint}=2: Bilinear interpolation (geographical grid)
-\item \np{nn\_2dint}=3: Bilinear remapping interpolation (general grid)
-\item \np{nn\_2dint}=4: Polynomial interpolation
-\item \np{nn\_2dint}=5: Radial footprint averaging with diameter specified in the namelist as \np{rn\_???\_avglamscl} in degrees or metres (set using \np{ln\_???\_fp\_indegs})
-\item \np{nn\_2dint}=6: Rectangular footprint averaging with E/W and N/S size specified in the namelist as \np{rn\_???\_avglamscl} and \np{rn\_???\_avgphiscl} in degrees or metres (set using \np{ln\_???\_fp\_indegs})
+\item \forcode{nn_2dint = 0}: Distance-weighted interpolation
+\item \forcode{nn_2dint = 1}: Distance-weighted interpolation (small angle)
+\item \forcode{nn_2dint = 2}: Bilinear interpolation (geographical grid)
+\item \forcode{nn_2dint = 3}: Bilinear remapping interpolation (general grid)
+\item \forcode{nn_2dint = 4}: Polynomial interpolation
+\item \forcode{nn_2dint = 5}: Radial footprint averaging with diameter specified in the namelist as \np{rn\_???\_avglamscl} in degrees or metres (set using \np{ln\_???\_fp\_indegs})
+\item \forcode{nn_2dint = 6}: Rectangular footprint averaging with E/W and N/S size specified in the namelist as \np{rn\_???\_avglamscl} and \np{rn\_???\_avgphiscl} in degrees or metres (set using \np{ln\_???\_fp\_indegs})
 \end{itemize}
 The ??? in the last two options indicate these options should be specified for each observation type for which the averaging is to be performed (see namelist example above).
-The \np{nn\_2dint} default option can be overridden for surface observation types using namelist values \np{nn\_2dint\_???} where ??? is one of sla,sst,sss,sic.
+The \np{nn_2dint} default option can be overridden for surface observation types using namelist values \np{nn\_2dint\_???} where ??? is one of sla,sst,sss,sic.
 
 Below is some more detail on the various options for interpolation and averaging available in NEMO.
@@ -956,5 +956,5 @@
 
 The above namelist will result in feedback files whose first 12 hours contain
-the first field of foo.nc and the second 12 hours contain the second field.
+the first field of \ifile{foo} and the second 12 hours contain the second field.
 
 %\begin{framed}
@@ -998,5 +998,5 @@
 \noindent
 \linebreak
-\textbf{\$\{prefix\}\_\$\{yyyymmdd\}\_\$\{sys\}\_\$\{cfg\}\_\$\{vn\}\_\$\{kind\}\_\$\{nproc\}.nc}
+\ifile{\textbf{\$\{prefix\}\_\$\{yyyymmdd\}\_\$\{sys\}\_\$\{cfg\}\_\$\{vn\}\_\$\{kind\}\_\$\{nproc\}}}
 
 \noindent
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_SBC.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_SBC.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_SBC.tex	(revision 9392)
@@ -28,13 +28,13 @@
 
 Five different ways to provide the first six fields to the ocean are available which 
-are controlled by namelist \ngn{namsbc} variables: an analytical formulation (\np{ln\_ana}~=~true), 
-a flux formulation (\np{ln\_flx}~=~true), a bulk formulae formulation (CORE 
-(\np{ln\_blk\_core}~=~true), CLIO (\np{ln\_blk\_clio}~=~true) or MFS
+are controlled by namelist \ngn{namsbc} variables: an analytical formulation (\np{ln_ana}~=~true), 
+a flux formulation (\np{ln_flx}~=~true), a bulk formulae formulation (CORE 
+(\np{ln_blk_core}~=~true), CLIO (\np{ln_blk_clio}~=~true) or MFS
 \footnote { Note that MFS bulk formulae compute fluxes only for the ocean component}
-(\np{ln\_blk\_mfs}~=~true) bulk formulae) and a coupled or mixed forced/coupled formulation 
-(exchanges with a atmospheric model via the OASIS coupler) (\np{ln\_cpl} or \np{ln\_mixcpl}~=~true). 
-When used ($i.e.$ \np{ln\_apr\_dyn}~=~true), the atmospheric pressure forces both ocean and ice dynamics.
-
-The frequency at which the forcing fields have to be updated is given by the \np{nn\_fsbc} namelist parameter. 
+(\np{ln_blk_mfs}~=~true) bulk formulae) and a coupled or mixed forced/coupled formulation 
+(exchanges with a atmospheric model via the OASIS coupler) (\np{ln_cpl} or \np{ln_mixcpl}~=~true). 
+When used ($i.e.$ \np{ln_apr_dyn}~=~true), the atmospheric pressure forces both ocean and ice dynamics.
+
+The frequency at which the forcing fields have to be updated is given by the \np{nn_fsbc} namelist parameter. 
 When the fields are supplied from data files (flux and bulk formulations), the input fields 
 need not be supplied on the model grid. Instead a file of coordinates and weights can 
@@ -50,11 +50,11 @@
 \item the rotation of vector components supplied relative to an east-north 
 coordinate system onto the local grid directions in the model ; 
-\item the addition of a surface restoring term to observed SST and/or SSS (\np{ln\_ssr}~=~true) ; 
-\item the modification of fluxes below ice-covered areas (using observed ice-cover or a sea-ice model) (\np{nn\_ice}~=~0,1, 2 or 3) ; 
-\item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln\_rnf}~=~true) ; 
-\item the addition of isf melting as lateral inflow (parameterisation) or as fluxes applied at the land-ice ocean interface (\np{ln\_isf}) ; 
-\item the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2) ; 
-\item the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle (\np{ln\_dm2dc}~=~true) ; 
-and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}~=~true). 
+\item the addition of a surface restoring term to observed SST and/or SSS (\np{ln_ssr}~=~true) ; 
+\item the modification of fluxes below ice-covered areas (using observed ice-cover or a sea-ice model) (\np{nn_ice}~=~0,1, 2 or 3) ; 
+\item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln_rnf}~=~true) ; 
+\item the addition of isf melting as lateral inflow (parameterisation) or as fluxes applied at the land-ice ocean interface (\np{ln_isf}) ; 
+\item the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn_fwb}~=~0,~1~or~2) ; 
+\item the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle (\np{ln_dm2dc}~=~true) ; 
+and a neutral drag coefficient can be read from an external wave model (\np{ln_cdgw}~=~true). 
 \end{itemize}
 The latter option is possible only in case core or mfs bulk formulas are selected.
@@ -91,5 +91,5 @@
 and \eqref{Eq_tra_sbc_lin} in \S\ref{TRA_sbc}). 
 The latter is the penetrative part of the heat flux. It is applied as a 3D 
-trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=\textit{true}.
+trends of the temperature equation (\mdl{traqsr} module) when \np{ln_traqsr}=\textit{true}.
 The way the light penetrates inside the water column is generally a sum of decreasing 
 exponentials (see \S\ref{TRA_qsr}). 
@@ -110,5 +110,5 @@
 %created!)
 %
-%Especially the \np{nn\_fsbc}, the \mdl{sbc\_oce} module (fluxes + mean sst sss ssu 
+%Especially the \np{nn_fsbc}, the \mdl{sbc\_oce} module (fluxes + mean sst sss ssu 
 %ssv) i.e. information required by flux computation or sea-ice
 %
@@ -130,7 +130,7 @@
 The ocean model provides, at each time step, to the surface module (\mdl{sbcmod}) 
 the surface currents, temperature and salinity.  
-These variables are averaged over \np{nn\_fsbc} time-step (\ref{Tab_ssm}), 
+These variables are averaged over \np{nn_fsbc} time-step (\ref{Tab_ssm}), 
 and it is these averaged fields which are used to computes the surface fluxes 
-at a frequency of \np{nn\_fsbc} time-step.
+at a frequency of \np{nn_fsbc} time-step.
 
 
@@ -185,5 +185,5 @@
 
 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} 
+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='./'.
@@ -214,6 +214,6 @@
 \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
+clim = false	& \ifile{fn\_yYYYYmMMdDD}  &   \ifile{fn\_yYYYYmMM}   &   \ifile{fn\_yYYYY}  \\   \hline
+clim = true	 	   & not possible 	              &  \ifile{fn\_m??}             &   fn                \\   \hline
 \end{tabular}
 \end{center}
@@ -271,7 +271,7 @@
 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, 
+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, 
+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 
@@ -438,5 +438,5 @@
 \item  Development of sea-ice algorithms or parameterizations.
 \item  spinup of the iceberg floats
-\item  ocean/sea-ice simulation with both media running in parallel (\np{ln\_mixcpl}~=~\textit{true})
+\item  ocean/sea-ice simulation with both media running in parallel (\np{ln_mixcpl}~=~\textit{true})
 \end{itemize}
 
@@ -492,8 +492,8 @@
 In this case, all the six fluxes needed by the ocean are assumed to 
 be uniform in space. They take constant values given in the namelist 
-\ngn{namsbc{\_}ana} by the variables \np{rn\_utau0}, \np{rn\_vtau0}, \np{rn\_qns0}, 
-\np{rn\_qsr0}, and \np{rn\_emp0} ($\textit{emp}=\textit{emp}_S$). The runoff is set to zero. 
+\ngn{namsbc{\_}ana} by the variables \np{rn_utau0}, \np{rn_vtau0}, \np{rn_qns0}, 
+\np{rn_qsr0}, and \np{rn_emp0} ($\textit{emp}=\textit{emp}_S$). The runoff is set to zero. 
 In addition, the wind is allowed to reach its nominal value within a given number 
-of time steps (\np{nn\_tau000}).
+of time steps (\np{nn_tau000}).
 
 If a user wants to apply a different analytical forcing, the \mdl{sbcana} 
@@ -513,5 +513,5 @@
 %-------------------------------------------------------------------------------------------------------------
 
-In the flux formulation (\np{ln\_flx}=true), the surface boundary 
+In the flux formulation (\forcode{ln_flx = .true.}), the surface boundary 
 condition fields are directly read from input files. The user has to define 
 in the namelist \ngn{namsbc{\_}flx} the name of the file, the name of the variable 
@@ -537,5 +537,5 @@
 The atmospheric fields used depend on the bulk formulae used. Three bulk formulations 
 are available : the CORE, the CLIO and the MFS bulk formulea. The choice is made by setting to true
-one of the following namelist variable : \np{ln\_core} ; \np{ln\_clio} or  \np{ln\_mfs}.
+one of the following namelist variable : \np{ln_core} ; \np{ln_clio} or  \np{ln_mfs}.
 
 Note : in forced mode, when a sea-ice model is used, a bulk formulation (CLIO or CORE) have to be used. 
@@ -546,6 +546,6 @@
 %        CORE Bulk formulea
 % -------------------------------------------------------------------------------------------------------------
-\subsection    [CORE Bulk formulea (\protect\np{ln\_core}=true)]
-		      {CORE Bulk formulea (\protect\np{ln\_core}=true, \protect\mdl{sbcblk\_core})}
+\subsection    [CORE Bulk formulea (\protect\forcode{ln_core = .true.})]
+		      {CORE Bulk formulea (\protect\forcode{ln_core = .true.}, \protect\mdl{sbcblk\_core})}
 \label{SBC_blk_core}
 %------------------------------------------namsbc_core----------------------------------------------------
@@ -592,17 +592,17 @@
 or larger than the one of the input atmospheric fields.
 
-The \np{sn\_wndi}, \np{sn\_wndj}, \np{sn\_qsr}, \np{sn\_qlw}, \np{sn\_tair}, \np{sn\_humi},
-\np{sn\_prec}, \np{sn\_snow}, \np{sn\_tdif} parameters describe the fields 
+The \np{sn_wndi}, \np{sn_wndj}, \np{sn_qsr}, \np{sn_qlw}, \np{sn_tair}, \np{sn_humi},
+\np{sn_prec}, \np{sn_snow}, \np{sn_tdif} parameters describe the fields 
 and the way they have to be used (spatial and temporal interpolations). 
 
-\np{cn\_dir} is the directory of location of bulk files
-\np{ln\_taudif} is the flag to specify if we use Hight Frequency (HF) tau information (.true.) or not (.false.)
-\np{rn\_zqt}: is the height of humidity and temperature measurements (m)
-\np{rn\_zu}: is the height of wind measurements (m)
+\np{cn_dir} is the directory of location of bulk files
+\np{ln_taudif} is the flag to specify if we use Hight Frequency (HF) tau information (.true.) or not (.false.)
+\np{rn_zqt}: is the height of humidity and temperature measurements (m)
+\np{rn_zu}: is the height of wind measurements (m)
 
 Three multiplicative factors are availables : 
-\np{rn\_pfac} and \np{rn\_efac} allows to adjust (if necessary) the global freshwater budget 
+\np{rn_pfac} and \np{rn_efac} allows to adjust (if necessary) the global freshwater budget 
 by increasing/reducing the precipitations (total and snow) and or evaporation, respectively.
-The third one,\np{rn\_vfac}, control to which extend the ice/ocean velocities are taken into account 
+The third one,\np{rn_vfac}, control to which extend the ice/ocean velocities are taken into account 
 in the calculation of surface wind stress. Its range should be between zero and one, 
 and it is recommended to set it to 0.
@@ -611,6 +611,6 @@
 %        CLIO Bulk formulea
 % -------------------------------------------------------------------------------------------------------------
-\subsection    [CLIO Bulk formulea (\protect\np{ln\_clio}=true)]
-		      {CLIO Bulk formulea (\protect\np{ln\_clio}=true, \protect\mdl{sbcblk\_clio})}
+\subsection    [CLIO Bulk formulea (\protect\forcode{ln_clio = .true.})]
+		      {CLIO Bulk formulea (\protect\forcode{ln_clio = .true.}, \protect\mdl{sbcblk\_clio})}
 \label{SBC_blk_clio}
 %------------------------------------------namsbc_clio----------------------------------------------------
@@ -652,6 +652,6 @@
 %        MFS Bulk formulae
 % -------------------------------------------------------------------------------------------------------------
-\subsection    [MFS Bulk formulea (\protect\np{ln\_mfs}=true)]
-		      {MFS Bulk formulea (\protect\np{ln\_mfs}=true, \protect\mdl{sbcblk\_mfs})}
+\subsection    [MFS Bulk formulea (\protect\forcode{ln_mfs = .true.})]
+		      {MFS Bulk formulea (\protect\forcode{ln_mfs = .true.}, \protect\mdl{sbcblk\_mfs})}
 \label{SBC_blk_mfs}
 %------------------------------------------namsbc_mfs----------------------------------------------------
@@ -679,11 +679,11 @@
 The required 7 input fields must be provided on the model Grid-T and  are:
 \begin{itemize}
-\item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windi})
-\item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windj})
-\item          Total Claud Cover (\%)  (\np{sn\_clc})
-\item          2m Air Temperature ($K$) (\np{sn\_tair})
-\item          2m Dew Point Temperature ($K$)  (\np{sn\_rhm})
-\item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn\_prec})
-\item          Mean Sea Level Pressure (${Pa}$) (\np{sn\_msl})
+\item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn_windi})
+\item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn_windj})
+\item          Total Claud Cover (\%)  (\np{sn_clc})
+\item          2m Air Temperature ($K$) (\np{sn_tair})
+\item          2m Dew Point Temperature ($K$)  (\np{sn_rhm})
+\item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn_prec})
+\item          Mean Sea Level Pressure (${Pa}$) (\np{sn_msl})
 \end{itemize}
 % -------------------------------------------------------------------------------------------------------------
@@ -709,5 +709,5 @@
 as well as to \href{http://wrf-model.org/}{WRF} (Weather Research and Forecasting Model).
 
-Note that in addition to the setting of \np{ln\_cpl} to true, the \key{coupled} have to be defined. 
+Note that in addition to the setting of \np{ln_cpl} to true, the \key{coupled} have to be defined. 
 The CPP key is mainly used in sea-ice to ensure that the atmospheric fluxes are 
 actually recieved by the ice-ocean system (no calculation of ice sublimation in coupled mode).
@@ -738,6 +738,6 @@
 
 The optional atmospheric pressure can be used to force ocean and ice dynamics 
-(\np{ln\_apr\_dyn}~=~true, \textit{\ngn{namsbc}} namelist ).
-The input atmospheric forcing defined via \np{sn\_apr} structure (\textit{namsbc\_apr} namelist) 
+(\np{ln_apr_dyn}~=~true, \textit{\ngn{namsbc}} namelist ).
+The input atmospheric forcing defined via \np{sn_apr} structure (\textit{namsbc\_apr} namelist) 
 can be interpolated in time to the model time step, and even in space when the 
 interpolation on-the-fly is used. When used to force the dynamics, the atmospheric 
@@ -748,5 +748,5 @@
 \end{equation}
 where $P_{atm}$ is the atmospheric pressure and $P_o$ a reference atmospheric pressure.
-A value of $101,000~N/m^2$ is used unless \np{ln\_ref\_apr} is set to true. In this case $P_o$ 
+A value of $101,000~N/m^2$ is used unless \np{ln_ref_apr} is set to true. In this case $P_o$ 
 is set to the value of $P_{atm}$ averaged over the ocean domain, $i.e.$ the mean value of 
 $\eta_{ib}$ is kept to zero at all time step.
@@ -760,5 +760,5 @@
 When using time-splitting and BDY package for open boundaries conditions, the equivalent 
 inverse barometer sea surface height $\eta_{ib}$ can be added to BDY ssh data: 
-\np{ln\_apr\_obc}  might be set to true.
+\np{ln_apr_obc}  might be set to true.
 
 % ================================================================
@@ -774,13 +774,13 @@
 
 A module is available to compute the tidal potential and use it in the momentum equation.
-This option is activated when \np{ln\_tide} is set to true in \ngn{nam\_tide}.
+This option is activated when \np{ln_tide} is set to true in \ngn{nam\_tide}.
 
 Some parameters are available in namelist \ngn{nam\_tide}:
 
-- \np{ln\_tide\_load} activate the load potential forcing and \np{filetide\_load} is  the associated file 
-
-- \np{ln\_tide\_pot} activate the tidal potential forcing
-
-- \np{nb\_harmo} is the number of constituent used
+- \np{ln_tide_load} activate the load potential forcing and \np{filetide_load} is  the associated file 
+
+- \np{ln_tide_pot} activate the tidal potential forcing
+
+- \np{nb_harmo} is the number of constituent used
 
 - \np{clname} is the name of constituent
@@ -863,5 +863,5 @@
 depth (in metres) which the river should be added to.
 
-Namelist variables in \ngn{namsbc\_rnf}, \np{ln\_rnf\_depth}, \np{ln\_rnf\_sal} and \np{ln\_rnf\_temp} control whether 
+Namelist variables in \ngn{namsbc\_rnf}, \np{ln_rnf_depth}, \np{ln_rnf_sal} and \np{ln_rnf_temp} control whether 
 the river attributes (depth, salinity and temperature) are read in and used.  If these are set 
 as false the river is added to the surface box only, assumed to be fresh (0~psu), and/or 
@@ -876,5 +876,5 @@
 to give the heat and salt content of the river runoff.
 After the user specified depth is read ini, the number of grid boxes this corresponds to is 
-calculated and stored in the variable \np{nz\_rnf}.
+calculated and stored in the variable \np{nz_rnf}.
 The variable \textit{h\_dep} is then calculated to be the depth (in metres) of the bottom of the 
 lowest box the river water is being added to (i.e. the total depth that river water is being added to in the model).
@@ -943,15 +943,15 @@
 \forfile{../namelists/namsbc_isf}
 %--------------------------------------------------------------------------------------------------------
-Namelist variable in \ngn{namsbc}, \np{nn\_isf}, controls the ice shelf representation used. 
+Namelist variable in \ngn{namsbc}, \np{nn_isf}, controls the ice shelf representation used. 
 \begin{description}
-\item[\np{nn\_isf}~=~1]
-The ice shelf cavity is represented (\np{ln\_isfcav}~=~true needed). The fwf and heat flux are computed. 
+\item[\np{nn_isf}~=~1]
+The ice shelf cavity is represented (\np{ln_isfcav}~=~true needed). The fwf and heat flux are computed. 
 Two different bulk formula are available:
 	\begin{description}
-	\item[\np{nn\_isfblk}~=~1]
+	\item[\np{nn_isfblk}~=~1]
 	The bulk formula used to compute the melt is based the one described in \citet{Hunter2006}.
         This formulation is based on a balance between the upward ocean heat flux and the latent heat flux at the ice shelf base.
 
-	\item[\np{nn\_isfblk}~=~2] 
+	\item[\np{nn_isfblk}~=~2] 
 	The bulk formula used to compute the melt is based the one described in \citet{Jenkins1991}.
         This formulation is based on a 3 equations formulation (a heat flux budget, a salt flux budget
@@ -962,9 +962,9 @@
 	\begin{description}
         \item[\np{nn\_gammablk~=~0~}]
-	The salt and heat exchange coefficients are constant and defined by \np{rn\_gammas0} and \np{rn\_gammat0}
+	The salt and heat exchange coefficients are constant and defined by \np{rn_gammas0} and \np{rn_gammat0}
 
 	\item[\np{nn\_gammablk~=~1~}]
 	The salt and heat exchange coefficients are velocity dependent and defined as $rn\_gammas0 \times u_{*}$ and $rn\_gammat0 \times u_{*}$
-        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn\_hisf\_tbl} meters).
+        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl} meters).
         See \citet{Jenkins2010} for all the details on this formulation.
 	
@@ -972,5 +972,5 @@
 	The salt and heat exchange coefficients are velocity and stability dependent and defined as 
         $\gamma_{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}}$
-        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn\_hisf\_tbl} meters), 
+        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl} meters), 
         $\Gamma_{Turb}$ the contribution of the ocean stability and 
         $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion.
@@ -978,38 +978,38 @@
         \end{description}
 
-\item[\np{nn\_isf}~=~2]
+\item[\np{nn_isf}~=~2]
 A parameterisation of isf is used. The ice shelf cavity is not represented. 
 The fwf is distributed along the ice shelf edge between the depth of the average grounding line (GL)
-(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}) as in (\np{nn\_isf}~=~3). 
+(\np{sn_depmax_isf}) and the base of the ice shelf along the calving front (\np{sn_depmin_isf}) as in (\np{nn_isf}~=~3). 
 Furthermore the fwf and heat flux are computed using the \citet{Beckmann2003} parameterisation of isf melting. 
-The effective melting length (\np{sn\_Leff\_isf}) is read from a file.
-
-\item[\np{nn\_isf}~=~3]
+The effective melting length (\np{sn_Leff_isf}) is read from a file.
+
+\item[\np{nn_isf}~=~3]
 A simple parameterisation of isf is used. The ice shelf cavity is not represented. 
-The fwf (\np{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between the depth of the average grounding line (GL)
-(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}). 
+The fwf (\np{sn_rnfisf}) is prescribed and distributed along the ice shelf edge between the depth of the average grounding line (GL)
+(\np{sn_depmax_isf}) and the base of the ice shelf along the calving front (\np{sn_depmin_isf}). 
 The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.
 
-\item[\np{nn\_isf}~=~4]
-The ice shelf cavity is opened (\np{ln\_isfcav}~=~true needed). However, the fwf is not computed but specified from file \np{sn\_fwfisf}). 
+\item[\np{nn_isf}~=~4]
+The ice shelf cavity is opened (\np{ln_isfcav}~=~true needed). However, the fwf is not computed but specified from file \np{sn_fwfisf}). 
 The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.\\
 \end{description}
 
 
-$\bullet$ \np{nn\_isf}~=~1 and \np{nn\_isf}~=~2 compute a melt rate based on the water mass properties, ocean velocities and depth.
+$\bullet$ \np{nn_isf}~=~1 and \np{nn_isf}~=~2 compute a melt rate based on the water mass properties, ocean velocities and depth.
  This flux is thus highly dependent of the model resolution (horizontal and vertical), realism of the water masses onto the shelf ...\\
 
 
-$\bullet$ \np{nn\_isf}~=~3 and \np{nn\_isf}~=~4 read the melt rate from a file. You have total control of the fwf forcing.
+$\bullet$ \np{nn_isf}~=~3 and \np{nn_isf}~=~4 read the melt rate from a file. You have total control of the fwf forcing.
 This can be usefull if the water masses on the shelf are not realistic or the resolution (horizontal/vertical) are too 
 coarse to have realistic melting or for studies where you need to control your heat and fw input.\\ 
 
 A namelist parameters control over how many meters the heat and fw fluxes are spread. 
-\np{rn\_hisf\_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}. 
-This parameter is only used if \np{nn\_isf}~=~1 or \np{nn\_isf}~=~4
-
-If \np{rn\_hisf\_tbl} = 0., the fluxes are put in the top level whatever is its tickness. 
-
-If \np{rn\_hisf\_tbl} $>$ 0., the fluxes are spread over the first \np{rn\_hisf\_tbl} m (ie over one or several cells).\\
+\np{rn_hisf_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}. 
+This parameter is only used if \np{nn_isf}~=~1 or \np{nn_isf}~=~4
+
+If \np{rn_hisf_tbl} = 0., the fluxes are put in the top level whatever is its tickness. 
+
+If \np{rn_hisf_tbl} $>$ 0., the fluxes are spread over the first \np{rn_hisf_tbl} m (ie over one or several cells).\\
 
 The ice shelf melt is implemented as a volume flux with in the same way as for the runoff.
@@ -1043,5 +1043,5 @@
 	set mask to 1, T/S is extrapolated from neighbours, ssh is extrapolated from neighbours and U/V set to 0. If no neighbour, T/S/U/V and mask set to 0.
 \end{description}
-The extrapolation is call \np{nn\_drown} times. It means that if the grounding line retreat by more than \np{nn\_drown} cells between 2 coupling steps,
+The extrapolation is call \np{nn_drown} times. It means that if the grounding line retreat by more than \np{nn_drown} cells between 2 coupling steps,
  the code will be unable to fill all the new wet cells properly. The default number is set up for the MISOMIP idealised experiments.\\
 This coupling procedure is able to take into account grounding line and calving front migration. However, it is a non-conservative processe. 
@@ -1049,6 +1049,6 @@
  a simple conservation scheme is available with \np{ln\_hsb = ~true}. The heat/salt/vol. gain/loss is diagnose, as well as the location. 
 Based on what is done on sbcrnf to prescribed a source of heat/salt/vol., the heat/salt/vol. gain/loss is removed/added,
- over a period of \np{rn\_fiscpl} time step, into the system. 
-So after \np{rn\_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\
+ over a period of \np{rn_fiscpl} time step, into the system. 
+So after \np{rn_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\
 
 As the before and now fields are not compatible (modification of the geometry), the restart time step is prescribed to be an euler time step instead of a leap frog and $fields_b = fields_n$.
@@ -1068,20 +1068,20 @@
 Icebergs are initially spawned into one of ten classes which have specific mass and thickness as described 
 in the \ngn{namberg} namelist: 
-\np{rn\_initial\_mass} and \np{rn\_initial\_thickness}.
-Each class has an associated scaling (\np{rn\_mass\_scaling}), which is an integer representing how many icebergs 
+\np{rn_initial_mass} and \np{rn_initial_thickness}.
+Each class has an associated scaling (\np{rn_mass_scaling}), which is an integer representing how many icebergs 
 of this class are being described as one lagrangian point (this reduces the numerical problem of tracking every single iceberg).
-They are enabled by setting \np{ln\_icebergs}~=~true.
+They are enabled by setting \np{ln_icebergs}~=~true.
 
 Two initialisation schemes are possible.
 \begin{description}
-\item[\np{nn\_test\_icebergs}~$>$~0]
-In this scheme, the value of \np{nn\_test\_icebergs} represents the class of iceberg to generate 
-(so between 1 and 10), and \np{nn\_test\_icebergs} provides a lon/lat box in the domain at each 
+\item[\np{nn_test_icebergs}~$>$~0]
+In this scheme, the value of \np{nn_test_icebergs} represents the class of iceberg to generate 
+(so between 1 and 10), and \np{nn_test_icebergs} provides a lon/lat box in the domain at each 
 grid point of which an iceberg is generated at the beginning of the run. 
-(Note that this happens each time the timestep equals \np{nn\_nit000}.)
-\np{nn\_test\_icebergs} is defined by four numbers in \np{nn\_test\_box} representing the corners 
+(Note that this happens each time the timestep equals \np{nn_nit000}.)
+\np{nn_test_icebergs} is defined by four numbers in \np{nn_test_box} representing the corners 
 of the geographical box: lonmin,lonmax,latmin,latmax
-\item[\np{nn\_test\_icebergs}~=~-1]
-In this scheme the model reads a calving file supplied in the \np{sn\_icb} parameter.
+\item[\np{nn_test_icebergs}~=~-1]
+In this scheme the model reads a calving file supplied in the \np{sn_icb} parameter.
 This should be a file with a field on the configuration grid (typically ORCA) representing ice accumulation rate at each model point. 
 These should be ocean points adjacent to land where icebergs are known to calve.
@@ -1095,5 +1095,5 @@
 Icebergs are influenced by wind, waves and currents, bottom melt and erosion.
 The latter act to disintegrate the iceberg. This is either all melted freshwater, or 
-(if \np{rn\_bits\_erosion\_fraction}~$>$~0) into melt and additionally small ice bits
+(if \np{rn_bits_erosion_fraction}~$>$~0) into melt and additionally small ice bits
 which are assumed to propagate with their larger parent and thus delay fluxing into the ocean.
 Melt water (and other variables on the configuration grid) are written into the main NEMO model output files.
@@ -1101,14 +1101,14 @@
 Extensive diagnostics can be produced.
 Separate output files are maintained for human-readable iceberg information.
-A separate file is produced for each processor (independent of \np{ln\_ctl}).
+A separate file is produced for each processor (independent of \np{ln_ctl}).
 The amount of information is controlled by two integer parameters:
 \begin{description}
-\item[\np{nn\_verbose\_level}]  takes a value between one and four and represents 
+\item[\np{nn_verbose_level}]  takes a value between one and four and represents 
 an increasing number of points in the code at which variables are written, and an 
 increasing level of obscurity.
-\item[\np{nn\_verbose\_write}] is the number of timesteps between writes
+\item[\np{nn_verbose_write}] is the number of timesteps between writes
 \end{description}
 
-Iceberg trajectories can also be written out and this is enabled by setting \np{nn\_sample\_rate}~$>$~0.
+Iceberg trajectories can also be written out and this is enabled by setting \np{nn_sample_rate}~$>$~0.
 A non-zero value represents how many timesteps between writes of information into the output file.
 These output files are in NETCDF format.
@@ -1155,6 +1155,6 @@
 the diurnal cycle of SWF is a scaling of the top of the atmosphere diurnal cycle 
 of incident SWF. The \cite{Bernie_al_CD07} reconstruction algorithm is available
-in \NEMO by setting \np{ln\_dm2dc}~=~true (a \textit{\ngn{namsbc}} namelist variable) when using 
-CORE bulk formulea (\np{ln\_blk\_core}~=~true) or the flux formulation (\np{ln\_flx}~=~true). 
+in \NEMO by setting \np{ln_dm2dc}~=~true (a \textit{\ngn{namsbc}} namelist variable) when using 
+CORE bulk formulea (\np{ln_blk_core}~=~true) or the flux formulation (\np{ln_flx}~=~true). 
 The reconstruction is performed in the \mdl{sbcdcy} module. The detail of the algoritm used 
 can be found in the appendix~A of \cite{Bernie_al_CD07}. The algorithm preserve the daily 
@@ -1162,5 +1162,5 @@
 of the analytical cycle over this time step (Fig.\ref{Fig_SBC_diurnal}). 
 The use of diurnal cycle reconstruction requires the input SWF to be daily 
-($i.e.$ a frequency of 24 and a time interpolation set to true in \np{sn\_qsr} namelist parameter).
+($i.e.$ a frequency of 24 and a time interpolation set to true in \np{sn_qsr} namelist parameter).
 Furthermore, it is recommended to have a least 8 surface module time step per day,
 that is  $\rdt \ nn\_fsbc < 10,800~s = 3~h$. An example of recontructed SWF 
@@ -1189,5 +1189,5 @@
 \label{SBC_rotation}
 
-When using a flux (\np{ln\_flx}=true) or bulk (\np{ln\_clio}=true or \np{ln\_core}=true) formulation, 
+When using a flux (\forcode{ln_flx = .true.}) or bulk (\forcode{ln_clio = .true.} or \forcode{ln_core = .true.}) formulation, 
 pairs of vector components can be rotated from east-north directions onto the local grid directions.  
 This is particularly useful when interpolation on the fly is used since here any vectors are likely to be defined 
@@ -1213,5 +1213,5 @@
 
 IOptions are defined through the  \ngn{namsbc\_ssr} namelist variables.
-n forced mode using a flux formulation (\np{ln\_flx}~=~true), a 
+n forced mode using a flux formulation (\np{ln_flx}~=~true), a 
 feedback term \emph{must} be added to the surface heat flux $Q_{ns}^o$:
 \begin{equation} \label{Eq_sbc_dmp_q}
@@ -1251,5 +1251,5 @@
 The presence at the sea surface of an ice covered area modifies all the fluxes 
 transmitted to the ocean. There are several way to handle sea-ice in the system 
-depending on the value of the \np{nn\_ice} namelist parameter found in \ngn{namsbc} namelist.  
+depending on the value of the \np{nn_ice} namelist parameter found in \ngn{namsbc} namelist.  
 \begin{description}
 \item[nn{\_}ice = 0]  there will never be sea-ice in the computational domain. 
@@ -1287,5 +1287,5 @@
 \textit{calc\_strair~=~true} and \textit{calc\_Tsfc~=~true} in the CICE name-list), or alternatively when NEMO 
 is coupled to the HadGAM3 atmosphere model (with \textit{calc\_strair~=~false} and \textit{calc\_Tsfc~=~false}).
-The code is intended to be used with \np{nn\_fsbc} set to 1 (although coupling ocean and ice less frequently 
+The code is intended to be used with \np{nn_fsbc} set to 1 (although coupling ocean and ice less frequently 
 should work, it is possible the calculation of some of the ocean-ice fluxes needs to be modified slightly - the
 user should check that results are not significantly different to the standard case).
@@ -1311,9 +1311,9 @@
 in the freshwater fluxes. In \NEMO, two way of controlling the the freshwater budget. 
 \begin{description}
-\item[\np{nn\_fwb}=0]  no control at all. The mean sea level is free to drift, and will 
+\item[\forcode{nn_fwb = 0}]  no control at all. The mean sea level is free to drift, and will 
 certainly do so.
-\item[\np{nn\_fwb}=1]  global mean \textit{emp} set to zero at each model time step. 
+\item[\forcode{nn_fwb = 1}]  global mean \textit{emp} set to zero at each model time step. 
 %Note that with a sea-ice model, this technique only control the mean sea level with linear free surface (\key{vvl} not defined) and no mass flux between ocean and ice (as it is implemented in the current ice-ocean coupling). 
-\item[\np{nn\_fwb}=2]  freshwater budget is adjusted from the previous year annual 
+\item[\forcode{nn_fwb = 2}]  freshwater budget is adjusted from the previous year annual 
 mean budget which is read in the \textit{EMPave\_old.dat} file. As the model uses the 
 Boussinesq approximation, the annual mean fresh water budget is simply evaluated 
@@ -1333,6 +1333,6 @@
 
 In order to read a neutral drag coeff, from an external data source ($i.e.$ a wave model), the 
-logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \textit{true}. 
-The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the
+logical variable \np{ln_cdgw} in \ngn{namsbc} namelist must be set to \textit{true}. 
+The \mdl{sbcwave} module containing the routine \np{sbc_wave} reads the
 namelist \ngn{namsbc\_wave} (for external data names, locations, frequency, interpolation and all 
 the miscellanous options allowed by Input Data generic Interface see \S\ref{SBC_input}) 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_STO.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_STO.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_STO.tex	(revision 9392)
@@ -155,22 +155,22 @@
 Parameters for the processes can be specified through the following \ngn{namsto} namelist parameters:
 \begin{description}
-   \item[\np{nn\_sto\_eos}]   : number of independent random walks 
-   \item[\np{rn\_eos\_stdxy}] : random walk horz. standard deviation (in grid points)
-   \item[\np{rn\_eos\_stdz}]  : random walk vert. standard deviation (in grid points)
-   \item[\np{rn\_eos\_tcor}]  : random walk time correlation (in timesteps)
-   \item[\np{nn\_eos\_ord}]   : order of autoregressive processes
-   \item[\np{nn\_eos\_flt}]   : passes of Laplacian filter
-   \item[\np{rn\_eos\_lim}]   : limitation factor (default = 3.0)
+   \item[\np{nn_sto_eos}]   : number of independent random walks 
+   \item[\np{rn_eos_stdxy}] : random walk horz. standard deviation (in grid points)
+   \item[\np{rn_eos_stdz}]  : random walk vert. standard deviation (in grid points)
+   \item[\np{rn_eos_tcor}]  : random walk time correlation (in timesteps)
+   \item[\np{nn_eos_ord}]   : order of autoregressive processes
+   \item[\np{nn_eos_flt}]   : passes of Laplacian filter
+   \item[\np{rn_eos_lim}]   : limitation factor (default = 3.0)
 \end{description}
 This routine also includes the initialization (seeding) of the random number generator.
 
 The third routine (\rou{sto\_rst\_write}) writes a restart file (which suffix name is 
-given by \np{cn\_storst\_out} namelist parameter) containing the current value of 
+given by \np{cn_storst_out} namelist parameter) containing the current value of 
 all autoregressive processes to allow restarting a simulation from where it has been interrupted.
 This file also contains the current state of the random number generator.
-When \np{ln\_rststo} is set to \textit{true}), the restart file (which suffix name is 
-given by \np{cn\_storst\_in} namelist parameter) is read by the initialization routine 
+When \np{ln_rststo} is set to \textit{true}), the restart file (which suffix name is 
+given by \np{cn_storst_in} namelist parameter) is read by the initialization routine 
 (\rou{sto\_par\_init}). The simulation will continue exactly as if it was not interrupted
-only  when \np{ln\_rstseed} is set to \textit{true}, $i.e.$ when the state of 
+only  when \np{ln_rstseed} is set to \textit{true}, $i.e.$ when the state of 
 the random number generator is read in the restart file.
 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_TRA.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_TRA.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_TRA.tex	(revision 9392)
@@ -57,5 +57,5 @@
 
 The user has the option of extracting each tendency term on the RHS of the tracer 
-equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{DIA}.
+equation for output (\np{ln_tra_trd} or \np{ln_tra_mxl}~=~true), as described in Chap.~\ref{DIA}.
 
 $\ $\newline    % force a new ligne
@@ -70,5 +70,5 @@
 %-------------------------------------------------------------------------------------------------------------
 
-When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \textit{true}), 
+When considered ($i.e.$ when \np{ln_traadv_NONE} is not set to \textit{true}), 
 the advection tendency of a tracer is expressed in flux form, 
 $i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by :
@@ -84,5 +84,5 @@
 by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
 which results from the use of the continuity equation,  $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 
-(which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}=true). 
+(which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \forcode{ln_linssh = .true.}). 
 Therefore it is of paramount importance to design the discrete analogue of the 
 advection tendency so that it is consistent with the continuity equation in order to 
@@ -114,5 +114,5 @@
 boundary condition depends on the type of sea surface chosen: 
 \begin{description}
-\item [linear free surface:] (\np{ln\_linssh}=true) the first level thickness is constant in time: 
+\item [linear free surface:] (\forcode{ln_linssh = .true.}) the first level thickness is constant in time: 
 the vertical boundary condition is applied at the fixed surface $z=0$ 
 rather than on the moving surface $z=\eta$. There is a non-zero advective 
@@ -120,5 +120,5 @@
 $\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$ 
 the product of surface velocity (at $z=0$) by the first level tracer value.
-\item [non-linear free surface:] (\np{ln\_linssh}=false) 
+\item [non-linear free surface:] (\forcode{ln_linssh = .false.}) 
 convergence/divergence in the first ocean level moves the free surface 
 up/down. There is no tracer advection through it so that the advective 
@@ -174,13 +174,13 @@
 %        2nd and 4th order centred schemes
 % -------------------------------------------------------------------------------------------------------------
-\subsection [Centred schemes (CEN) (\protect\np{ln\_traadv\_cen})]
-			   {Centred schemes (CEN) (\protect\np{ln\_traadv\_cen}=true)}
+\subsection [Centred schemes (CEN) (\protect\np{ln_traadv_cen})]
+			   {Centred schemes (CEN) (\protect\forcode{ln_traadv_cen = .true.})}
 \label{TRA_adv_cen}
 
 %        2nd order centred scheme  
 
-The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~=~\textit{true}. 
+The centred advection scheme (CEN) is used when \np{ln_traadv_cen}~=~\textit{true}. 
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) 
-and vertical direction by setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$. 
+and vertical direction by setting \np{nn_cen_h} and \np{nn_cen_v} to $2$ or $4$. 
 CEN implementation can be found in the \mdl{traadv\_cen} module.
 
@@ -212,5 +212,5 @@
 =\overline{   T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2}
 \end{equation}
-In the vertical direction (\np{nn\_cen\_v}=$4$), a $4^{th}$ COMPACT interpolation 
+In the vertical direction (\np{nn_cen_v}=$4$), a $4^{th}$ COMPACT interpolation 
 has been prefered \citep{Demange_PhD2014}.
 In the COMPACT scheme, both the field and its derivative are interpolated, 
@@ -246,11 +246,11 @@
 %        FCT scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Flux Corrected Transport schemes (FCT) (\protect\np{ln\_traadv\_fct})]
-			{Flux Corrected Transport schemes (FCT) (\protect\np{ln\_traadv\_fct}=true)}
+\subsection   [Flux Corrected Transport schemes (FCT) (\protect\np{ln_traadv_fct})]
+			{Flux Corrected Transport schemes (FCT) (\protect\forcode{ln_traadv_fct = .true.})}
 \label{TRA_adv_tvd}
 
-The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~=~\textit{true}. 
+The Flux Corrected Transport schemes (FCT) is used when \np{ln_traadv_fct}~=~\textit{true}. 
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) 
-and vertical direction by setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$.
+and vertical direction by setting \np{nn_fct_h} and \np{nn_fct_v} to $2$ or $4$.
 FCT implementation can be found in the \mdl{traadv\_fct} module.
 
@@ -269,5 +269,5 @@
 where $c_u$ is a flux limiter function taking values between 0 and 1. 
 The FCT order is the one of the centred scheme used ($i.e.$ it depends on the setting of
-\np{nn\_fct\_h} and \np{nn\_fct\_v}.
+\np{nn_fct_h} and \np{nn_fct_v}.
 There exist many ways to define $c_u$, each corresponding to a different 
 FCT scheme. The one chosen in \NEMO is described in \citet{Zalesak_JCP79}. 
@@ -277,8 +277,8 @@
 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}. 
 
-An additional option has been added controlled by \np{nn\_fct\_zts}. By setting this integer to 
+An additional option has been added controlled by \np{nn_fct_zts}. By setting this integer to 
 a value larger than zero, a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction, 
 but on the latter, a split-explicit time stepping is used, with a number of sub-timestep equals
-to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited 
+to \np{nn_fct_zts}. This option can be useful when the size of the timestep is limited 
 by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar split-explicit 
 time stepping should be used on vertical advection of momentum to insure a better stability
@@ -293,9 +293,9 @@
 %        MUSCL scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection[MUSCL scheme  (\protect\np{ln\_traadv\_mus})]
-	{Monotone Upstream Scheme for Conservative Laws (MUSCL) (\protect\np{ln\_traadv\_mus}=T)}
+\subsection[MUSCL scheme  (\protect\np{ln_traadv_mus})]
+	{Monotone Upstream Scheme for Conservative Laws (MUSCL) (\protect\forcode{ln_traadv_mus = .true.})}
 \label{TRA_adv_mus}
 
-The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~=~\textit{true}. 
+The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln_traadv_mus}~=~\textit{true}. 
 MUSCL implementation can be found in the \mdl{traadv\_mus} module.
 
@@ -321,14 +321,14 @@
 the \textit{positive} character of the scheme. 
 In addition, fluxes round a grid-point where a runoff is applied can optionally be 
-computed using upstream fluxes (\np{ln\_mus\_ups}~=~\textit{true}).
+computed using upstream fluxes (\np{ln_mus_ups}~=~\textit{true}).
 
 % -------------------------------------------------------------------------------------------------------------
 %        UBS scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Upstream-Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs})]
-			{Upstream-Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs}=true)}
+\subsection   [Upstream-Biased Scheme (UBS) (\protect\np{ln_traadv_ubs})]
+			{Upstream-Biased Scheme (UBS) (\protect\forcode{ln_traadv_ubs = .true.})}
 \label{TRA_adv_ubs}
 
-The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~=~\textit{true}. 
+The Upstream-Biased Scheme (UBS) is used when \np{ln_traadv_ubs}~=~\textit{true}. 
 UBS implementation can be found in the \mdl{traadv\_mus} module.
 
@@ -358,5 +358,5 @@
 where the control of artificial diapycnal fluxes is of paramount importance \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}. 
 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme 
-or a $4^th$ order COMPACT scheme (\np{nn\_cen\_v}=2 or 4).
+or a $4^th$ order COMPACT scheme (\forcode{nn_cen_v = 2} or 4).
 
 For stability reasons  (see \S\ref{STP}),
@@ -401,11 +401,11 @@
 %        QCK scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [QUICKEST scheme (QCK) (\protect\np{ln\_traadv\_qck})]
-			{QUICKEST scheme (QCK) (\protect\np{ln\_traadv\_qck}=true)}
+\subsection   [QUICKEST scheme (QCK) (\protect\np{ln_traadv_qck})]
+			{QUICKEST scheme (QCK) (\protect\forcode{ln_traadv_qck = .true.})}
 \label{TRA_adv_qck}
 
 The Quadratic Upstream Interpolation for Convective Kinematics with 
 Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} 
-is used when \np{ln\_traadv\_qck}~=~\textit{true}. 
+is used when \np{ln_traadv_qck}~=~\textit{true}. 
 QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
 
@@ -449,5 +449,5 @@
 except for the pure vertical component that appears when a rotation tensor is used. 
 This latter component is solved implicitly together with the vertical diffusion term (see \S\ref{STP}). 
-When \np{ln\_traldf\_msc}~=~\textit{true}, a Method of Stabilizing Correction is used in which 
+When \np{ln_traldf_msc}~=~\textit{true}, a Method of Stabilizing Correction is used in which 
 the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
 
@@ -456,17 +456,17 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Type of operator (\protect\np{ln\_traldf\{\_NONE, \_lap, \_blp\}})]
-			     {Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp} = true) } 
+			     {Type of operator (\protect\np{ln_traldf_NONE}, \protect\np{ln_traldf_lap}, or \protect\np{ln_traldf_blp} = true) } 
 \label{TRA_ldf_op}
 
 Three operator options are proposed and, one and only one of them must be selected:
 \begin{description}
-\item [\np{ln\_traldf\_NONE}] = true : no operator selected, the lateral diffusive tendency will not be 
+\item [\np{ln_traldf_NONE}] = true : no operator selected, the lateral diffusive tendency will not be 
 applied to the tracer equation. This option can be used when the selected advection scheme 
 is diffusive enough (MUSCL scheme for example).
-\item [ \np{ln\_traldf\_lap}] = true : a laplacian operator is selected. This harmonic operator 
+\item [ \np{ln_traldf_lap}] = true : a laplacian operator is selected. This harmonic operator 
 takes the following expression:  $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $, 
 where the gradient operates along the selected direction (see \S\ref{TRA_ldf_dir}),
 and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see Chap.~\ref{LDF}).
-\item [\np{ln\_traldf\_blp}] = true : a bilaplacian operator is selected. This biharmonic operator 
+\item [\np{ln_traldf_blp}] = true : a bilaplacian operator is selected. This biharmonic operator 
 takes the following expression:  
 $\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 
@@ -489,19 +489,19 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Direction of action (\protect\np{ln\_traldf\{\_lev, \_hor, \_iso, \_triad\}})]
-			     {Direction of action (\protect\np{ln\_traldf\_lev}, \textit{...\_hor}, \textit{...\_iso}, or \textit{...\_triad} = true) } 
+			     {Direction of action (\protect\np{ln_traldf_lev}, \textit{...\_hor}, \textit{...\_iso}, or \textit{...\_triad} = true) } 
 \label{TRA_ldf_dir}
 
 The choice of a direction of action determines the form of operator used. 
 The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane 
-when iso-level option is used (\np{ln\_traldf\_lev}~=~\textit{true})
+when iso-level option is used (\np{ln_traldf_lev}~=~\textit{true})
 or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}-coordinate 
-(\np{ln\_traldf\_hor} and \np{ln\_zco} equal \textit{true}). 
+(\np{ln_traldf_hor} and \np{ln_zco} equal \textit{true}). 
 The associated code can be found in the \mdl{traldf\_lap\_blp} module.
 The operator is a rotated (re-entrant) laplacian when the direction along which it acts 
 does not coincide with the iso-level surfaces, 
-that is when standard or triad iso-neutral option is used (\np{ln\_traldf\_iso} or 
- \np{ln\_traldf\_triad} equals \textit{true}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.), 
+that is when standard or triad iso-neutral option is used (\np{ln_traldf_iso} or 
+ \np{ln_traldf_triad} equals \textit{true}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.), 
 or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}-coordinate 
-(\np{ln\_traldf\_hor} and \np{ln\_sco} equal \textit{true})
+(\np{ln_traldf_hor} and \np{ln_sco} equal \textit{true})
 \footnote{In this case, the standard iso-neutral operator will be automatically selected}. 
 In that case, a rotation is applied to the gradient(s) that appears in the operator 
@@ -515,6 +515,6 @@
 %       iso-level operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Iso-level (bi-)laplacian operator ( \protect\np{ln\_traldf\_iso})]
-			{Iso-level (bi-)laplacian operator ( \protect\np{ln\_traldf\_iso}) }
+\subsection   [Iso-level (bi-)laplacian operator ( \protect\np{ln_traldf_iso})]
+			{Iso-level (bi-)laplacian operator ( \protect\np{ln_traldf_iso}) }
 \label{TRA_ldf_lev}
 
@@ -534,10 +534,10 @@
 It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate 
 with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 
-It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~=~\textit{true}, 
-we have \np{ln\_traldf\_lev}~=~\textit{true} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~=~\textit{true}. 
+It is thus used when, in addition to \np{ln_traldf_lap} or \np{ln_traldf_blp}~=~\textit{true}, 
+we have \np{ln_traldf_lev}~=~\textit{true} or \np{ln_traldf_hor}~=~\np{ln_zco}~=~\textit{true}. 
 In both cases, it significantly contributes to diapycnal mixing. 
 It is therefore never recommended, even when using it in the bilaplacian case.
 
-Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), tracers in horizontally 
+Note that in the partial step $z$-coordinate (\forcode{ln_zps = .true.}), tracers in horizontally 
 adjacent cells are located at different depths in the vicinity of the bottom. 
 In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level 
@@ -584,6 +584,6 @@
 ($z$- or $s$-surfaces) and the surface along which the diffusion operator 
 acts ($i.e.$ horizontal or iso-neutral surfaces).  It is thus used when, 
-in addition to \np{ln\_traldf\_lap}= true, we have \np{ln\_traldf\_iso}=true, 
-or both \np{ln\_traldf\_hor}=true and \np{ln\_zco}=true. The way these 
+in addition to \np{ln_traldf_lap}= true, we have \forcode{ln_traldf_iso = .true.}, 
+or both \forcode{ln_traldf_hor = .true.} and \forcode{ln_zco = .true.}. The way these 
 slopes are evaluated is given in \S\ref{LDF_slp}. At the surface, bottom 
 and lateral boundaries, the turbulent fluxes of heat and salt are set to zero 
@@ -603,5 +603,5 @@
 background horizontal diffusion \citep{Guilyardi_al_CD01}. 
 
-Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal derivatives 
+Note that in the partial step $z$-coordinate (\forcode{ln_zps = .true.}), the horizontal derivatives 
 at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific treatment. 
 They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
@@ -609,12 +609,12 @@
 %&&     Triad rotated (bi-)laplacian operator
 %&&  -------------------------------------------
-\subsubsection   [Triad rotated (bi-)laplacian operator (\protect\np{ln\_traldf\_triad})]
-		        	  {Triad rotated (bi-)laplacian operator (\protect\np{ln\_traldf\_triad})}
+\subsubsection   [Triad rotated (bi-)laplacian operator (\protect\np{ln_traldf_triad})]
+		        	  {Triad rotated (bi-)laplacian operator (\protect\np{ln_traldf_triad})}
 \label{TRA_ldf_triad}
 
-If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}=true ; see App.\ref{sec:triad}) 
+If the Griffies triad scheme is employed (\forcode{ln_traldf_triad = .true.} ; see App.\ref{sec:triad}) 
 
 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 
-is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of 
+is also available in \NEMO (\forcode{ln_traldf_grif = .true.}). A complete description of 
 the algorithm is given in App.\ref{sec:triad}.
 
@@ -635,13 +635,13 @@
 \label{TRA_ldf_options}
 
-\np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
-
-\np{rn\_slpmax} = slope limit (both operators)
-
-\np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only)
-
-\np{rn\_sw\_triad} =1 switching triad ; =0 all 4 triads used (triad only) 
-
-\np{ln\_botmix\_triad} = lateral mixing on bottom (triad only)
+\np{ln_traldf_msc} = Method of Stabilizing Correction (both operators)
+
+\np{rn_slpmax} = slope limit (both operators)
+
+\np{ln_triad_iso} = pure horizontal mixing in ML (triad only)
+
+\np{rn_sw_triad} =1 switching triad ; =0 all 4 triads used (triad only) 
+
+\np{ln_botmix_triad} = lateral mixing on bottom (triad only)
 
 % ================================================================
@@ -685,10 +685,10 @@
 The large eddy coefficient found in the mixed layer together with high 
 vertical resolution implies that in the case of explicit time stepping 
-(\np{ln\_zdfexp}=true) there would be too restrictive a constraint on 
+(\forcode{ln_zdfexp = .true.}) there would be too restrictive a constraint on 
 the time step. Therefore, the default implicit time stepping is preferred 
 for the vertical diffusion since it overcomes the stability constraint. 
-A forward time differencing scheme (\np{ln\_zdfexp}=true) using a time 
-splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative. 
-Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both 
+A forward time differencing scheme (\forcode{ln_zdfexp = .true.}) using a time 
+splitting technique (\np{nn_zdfexp} $> 1$) is provided as an alternative. 
+Namelist variables \np{ln_zdfexp} and \np{nn_zdfexp} apply to both 
 tracers and dynamics. 
 
@@ -750,5 +750,5 @@
 divergence of odd and even time step (see \S\ref{STP}).
 
-In the linear free surface case (\np{ln\_linssh}~=~\textit{true}), 
+In the linear free surface case (\np{ln_linssh}~=~\textit{true}), 
 an additional term has to be added on both temperature and salinity. 
 On temperature, this term remove the heat content associated with mass exchange
@@ -781,7 +781,7 @@
 
 Options are defined through the  \ngn{namtra\_qsr} namelist variables.
-When the penetrative solar radiation option is used (\np{ln\_flxqsr}=true), 
+When the penetrative solar radiation option is used (\forcode{ln_flxqsr = .true.}), 
 the solar radiation penetrates the top few tens of meters of the ocean. If it is not used 
-(\np{ln\_flxqsr}=false) all the heat flux is absorbed in the first ocean level. 
+(\forcode{ln_flxqsr = .false.}) all the heat flux is absorbed in the first ocean level. 
 Thus, in the former case a term is added to the time evolution equation of 
 temperature \eqref{Eq_PE_tra_T} and the surface boundary condition is 
@@ -805,12 +805,12 @@
 wavelengths contribute to heating the upper few tens of centimetres. The fraction of $Q_{sr}$ 
 that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ (specified 
-through namelist parameter \np{rn\_abs}).  It is assumed to penetrate the ocean 
+through namelist parameter \np{rn_abs}).  It is assumed to penetrate the ocean 
 with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 
-of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namelist).
+of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn_si0} in the namtra\_qsr namelist).
 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy 
 propagates to larger depths where it contributes to 
 local heating. 
 The way this second part of the solar energy penetrates into the ocean depends on 
-which formulation is chosen. In the simple 2-waveband light penetration scheme  (\np{ln\_qsr\_2bd}=true) 
+which formulation is chosen. In the simple 2-waveband light penetration scheme  (\forcode{ln_qsr_2bd = .true.}) 
 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 
 leading to the following expression  \citep{Paulson1977}:
@@ -819,5 +819,5 @@
 \end{equation}
 where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths.  
-It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter. 
+It is usually chosen to be 23~m by setting the \np{rn_si0} namelist parameter. 
 The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in 
 Jerlov's (1968) classification (oligotrophic waters).
@@ -839,18 +839,18 @@
 computational efficiency. The 2-bands formulation does not reproduce the full model very well. 
 
-The RGB formulation is used when \np{ln\_qsr\_rgb}=true. The RGB attenuation coefficients
+The RGB formulation is used when \forcode{ln_qsr_rgb = .true.}. The RGB attenuation coefficients
 ($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform 
 chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb} 
 in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation:
 \begin{description} 
-\item[\np{nn\_chdta}=0] 
+\item[\forcode{nn_chdta = 0}] 
 a constant 0.05 g.Chl/L value everywhere ; 
-\item[\np{nn\_chdta}=1]  
+\item[\forcode{nn_chdta = 1}]  
 an observed time varying chlorophyll deduced from satellite surface ocean color measurement 
 spread uniformly in the vertical direction ; 
-\item[\np{nn\_chdta}=2]  
+\item[\forcode{nn_chdta = 2}]  
 same as previous case except that a vertical profile of chlorophyl is used. 
 Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ;
-\item[\np{ln\_qsr\_bio}=true]  
+\item[\forcode{ln_qsr_bio = .true.}]  
 simulated time varying chlorophyll by TOP biogeochemical model. 
 In this case, the RGB formulation is used to calculate both the phytoplankton 
@@ -913,8 +913,8 @@
 Options are defined through the  \ngn{namtra\_bbc} namelist variables.
 The presence of geothermal heating is controlled by setting the namelist 
-parameter  \np{ln\_trabbc} to true. Then, when \np{nn\_geoflx} is set to 1, 
+parameter  \np{ln_trabbc} to true. Then, when \np{nn_geoflx} is set to 1, 
 a constant geothermal heating is introduced whose value is given by the 
-\np{nn\_geoflx\_cst}, which is also a namelist parameter. 
-When  \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is 
+\np{nn_geoflx_cst}, which is also a namelist parameter. 
+When  \np{nn_geoflx} is set to 2, a spatially varying geothermal heat flux is 
 introduced which is provided in the \ifile{geothermal\_heating} NetCDF file 
 (Fig.\ref{Fig_geothermal}) \citep{Emile-Geay_Madec_OS09}.
@@ -959,8 +959,8 @@
 %        Diffusive BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Diffusive Bottom Boundary layer (\protect\np{nn\_bbl\_ldf}=1)}
+\subsection{Diffusive Bottom Boundary layer (\protect\forcode{nn_bbl_ldf = 1})}
 \label{TRA_bbl_diff}
 
-When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1), 
+When applying sigma-diffusion (\key{trabbl} defined and \np{nn_bbl_ldf} set to 1), 
 the diffusive flux between two adjacent cells at the ocean floor is given by 
 \begin{equation} \label{Eq_tra_bbl_diff}
@@ -978,5 +978,5 @@
 \end{equation} 
 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist 
-parameter \np{rn\_ahtbbl} and usually set to a value much larger 
+parameter \np{rn_ahtbbl} and usually set to a value much larger 
 than the one used for lateral mixing in the open ocean. The constraint in \eqref{Eq_tra_bbl_coef} 
 implies that sigma-like diffusion only occurs when the density above the sea floor, at the top of 
@@ -994,5 +994,5 @@
 %        Advective BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection   {Advective Bottom Boundary Layer  (\protect\np{nn\_bbl\_adv}= 1 or 2)}
+\subsection   {Advective Bottom Boundary Layer  (\protect\np{nn_bbl_adv}= 1 or 2)}
 \label{TRA_bbl_adv}
 
@@ -1022,10 +1022,10 @@
 %%%gmcomment   :  this section has to be really written
 
-When applying an advective BBL (\np{nn\_bbl\_adv} = 1 or 2), an overturning 
+When applying an advective BBL (\np{nn_bbl_adv} = 1 or 2), an overturning 
 circulation is added which connects two adjacent bottom grid-points only if dense 
 water overlies less dense water on the slope. The density difference causes dense 
 water to move down the slope. 
 
-\np{nn\_bbl\_adv} = 1 : the downslope velocity is chosen to be the Eulerian
+\np{nn_bbl_adv} = 1 : the downslope velocity is chosen to be the Eulerian
 ocean velocity just above the topographic step (see black arrow in Fig.\ref{Fig_bbl}) 
 \citep{Beckmann_Doscher1997}. It is a \textit{conditional advection}, that is, advection
@@ -1034,5 +1034,5 @@
 greater depth ($i.e.$ $\vect{U}  \cdot  \nabla H>0$).
 
-\np{nn\_bbl\_adv} = 2 : the downslope velocity is chosen to be proportional to $\Delta \rho$,
+\np{nn_bbl_adv} = 2 : the downslope velocity is chosen to be proportional to $\Delta \rho$,
 the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}.
 The advection is allowed only  if dense water overlies less dense water on the slope ($i.e.$ 
@@ -1044,5 +1044,5 @@
 \end{equation}
 where $\gamma$, expressed in seconds, is the coefficient of proportionality 
-provided as \np{rn\_gambbl}, a namelist parameter, and \textit{kup} and \textit{kdwn} 
+provided as \np{rn_gambbl}, a namelist parameter, and \textit{kup} and \textit{kdwn} 
 are the vertical index of the higher and lower cells, respectively.
 The parameter $\gamma$ should take a different value for each bathymetric 
@@ -1101,10 +1101,10 @@
 are given temperature and salinity fields (usually a climatology). 
 Options are defined through the  \ngn{namtra\_dmp} namelist variables.
-The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true. 
-It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true
-in \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are 
+The restoring term is added when the namelist parameter \np{ln_tradmp} is set to true. 
+It also requires that both \np{ln_tsd_init} and \np{ln_tsd_tradmp} are set to true
+in \textit{namtsd} namelist as well as \np{sn_tem} and \np{sn_sal} structures are 
 correctly set  ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read 
 using \mdl{fldread}, see \S\ref{SBC_fldread}). 
-The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. The file name is specified by the namelist variable \np{cn\_resto}. The DMP\_TOOLS tool is provided to allow users to generate the netcdf file.
+The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. The file name is specified by the namelist variable \np{cn_resto}. The DMP\_TOOLS tool is provided to allow users to generate the netcdf file.
 
 The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} 
@@ -1128,13 +1128,13 @@
 by stabilising the water column too much.
 
-The namelist parameter \np{nn\_zdmp} sets whether the damping should be applied in the whole water column or only below the mixed layer (defined either on a density or $S_o$ criterion). It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here \citep{Madec_al_JPO96}.
-
-\subsection[DMP\_TOOLS]{Generating resto.nc using DMP\_TOOLS}
+The namelist parameter \np{nn_zdmp} sets whether the damping should be applied in the whole water column or only below the mixed layer (defined either on a density or $S_o$ criterion). It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here \citep{Madec_al_JPO96}.
+
+\subsection[DMP\_TOOLS]{Generating \ifile{resto} using DMP\_TOOLS}
 
 DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. 
 Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled 
-and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input. 
-This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. 
-The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. 
+and run on the same machine as the NEMO model. A \ifile{mesh\_mask} file for the model configuration is required as an input. 
+This can be generated by carrying out a short model run with the namelist parameter \np{nn_msh} set to 1. 
+The namelist parameter \np{ln_tradmp} will also need to be set to .false. for this to work. 
 The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
 
@@ -1143,25 +1143,25 @@
 %-------------------------------------------------------------------------------------------------------
 
-\np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in the \nl{nam\_zoom\_dmp} name list.
+\np{cp_cfg}, \np{cp_cpz}, \np{jp_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom_n}, \np{lzoom_e}, \np{lzoom_s}, \np{lzoom_w} in the \nl{nam\_zoom\_dmp} name list.
 
 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations. 
-\np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. 
-\np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea 
+\np{ln_full_field} specifies that newtonian damping should be applied to the whole model domain. 
+\np{ln_med_red_seas} specifies grid specific restoration coefficients in the Mediterranean Sea 
 for the ORCA4, ORCA2 and ORCA05 configurations. 
-If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as 
+If \np{ln_old_31_lev_code} is set then the depth variation of the coeffients will be specified as 
 a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference 
 configurations with previous model versions. 
-\np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. 
-This option only has an effect if \np{ln\_full\_field} is true. 
-\np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. 
-Finally \np{ln\_custom} specifies that the custom module will be called. 
+\np{ln_coast} specifies that the restoration coefficient should be reduced near to coastlines. 
+This option only has an effect if \np{ln_full_field} is true. 
+\np{ln_zero_top_layer} specifies that the restoration coefficient should be zero in the surface layer. 
+Finally \np{ln_custom} specifies that the custom module will be called. 
 This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.
 
-The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}. 
+The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn_hdmp}. 
 Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to 
 the full values of a 10\deg latitud band. 
 This is often used because of the short adjustment time scale in the equatorial region 
 \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a 
-hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.  
+hyperbolic tangent, with \np{rn_surf} as surface value, \np{rn_bot} as bottom value and a transition depth of \np{rn_dep}.  
 
 % ================================================================
@@ -1191,6 +1191,6 @@
 the subscript $f$ denotes filtered values, $\gamma$ is the Asselin coefficient,
 and $S$ is the total forcing applied on $T$ ($i.e.$ fluxes plus content in mass exchanges). 
-$\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 
-Its default value is \np{rn\_atfp}=$10^{-3}$. Note that the forcing correction term in the filter
+$\gamma$ is initialized as \np{rn_atfp} (\textbf{namelist} parameter). 
+Its default value is \np{rn_atfp}=$10^{-3}$. Note that the forcing correction term in the filter
 is not applied in linear free surface (\jp{lk\_vvl}=false) (see \S\ref{TRA_sbc}.
 Not also that in constant volume case, the time stepping is performed on $T$, 
@@ -1217,5 +1217,5 @@
 %        Equation of State
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Equation Of Seawater (\protect\np{nn\_eos} = -1, 0, or 1)}
+\subsection{Equation Of Seawater (\protect\np{nn_eos} = -1, 0, or 1)}
 \label{TRA_eos}
 
@@ -1248,9 +1248,9 @@
 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}.
 
-Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos} 
+Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn_eos} 
 which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS).
 \begin{description}
 
-\item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.  
+\item[\np{nn_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.  
 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 
 but it is optimized for a boussinesq fluid and the polynomial expressions have simpler 
@@ -1268,9 +1268,9 @@
 $\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as 
 \textit{Conservative} Temperature and \textit{Absolute} Salinity. 
-In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST 
+In addition, setting \np{ln_useCT} to \textit{true} convert the Conservative SST to potential SST 
 prior to either computing the air-sea and ice-sea fluxes (forced mode) 
 or sending the SST field to the atmosphere (coupled mode).
 
-\item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used.
+\item[\np{nn_eos}$=0$] the polyEOS80-bsq equation of seawater is used.
 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized 
 to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80 
@@ -1283,5 +1283,5 @@
 value, the TEOS10 value. 
  
-\item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 
+\item[\np{nn_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 
 the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.) 
 (see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both 
@@ -1315,11 +1315,11 @@
 \hline
 coeff.	& computer name   & S-EOS		&  description		                  \\ \hline
-$a_0$       & \np{rn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff. 	\\ \hline
-$b_0$	      & \np{rn\_b0} 	   & 7.6554 $10^{-1}$ &  linear haline  expansion coeff. 	\\ \hline
-$\lambda_1$	& \np{rn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$ 	      \\ \hline
-$\lambda_2$	& \np{rn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$	 	   \\ \hline
-$\nu$       & \np{rn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$ 	   \\ \hline
-$\mu_1$     & \np{rn\_mu1} 	& 1.4970 $10^{-4}$ &  thermobaric coeff. in T    	   \\ \hline
-$\mu_2$     & \np{rn\_mu2} 	& 1.1090 $10^{-5}$ &  thermobaric coeff. in S   	      \\ \hline
+$a_0$       & \np{rn_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff. 	\\ \hline
+$b_0$	      & \np{rn_b0} 	   & 7.6554 $10^{-1}$ &  linear haline  expansion coeff. 	\\ \hline
+$\lambda_1$	& \np{rn_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$ 	      \\ \hline
+$\lambda_2$	& \np{rn_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$	 	   \\ \hline
+$\nu$       & \np{rn_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$ 	   \\ \hline
+$\mu_1$     & \np{rn_mu1} 	& 1.4970 $10^{-4}$ &  thermobaric coeff. in T    	   \\ \hline
+$\mu_2$     & \np{rn_mu2} 	& 1.1090 $10^{-5}$ &  thermobaric coeff. in S   	      \\ \hline
 \end{tabular}
 \caption{ \protect\label{Tab_SEOS}
@@ -1333,5 +1333,5 @@
 %        Brunt-V\"{a}is\"{a}l\"{a} Frequency
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Brunt-V\"{a}is\"{a}l\"{a} Frequency (\protect\np{nn\_eos} = 0, 1 or 2)}
+\subsection{Brunt-V\"{a}is\"{a}l\"{a} Frequency (\protect\np{nn_eos} = 0, 1 or 2)}
 \label{TRA_bn2}
 
@@ -1395,9 +1395,9 @@
                    I've changed "derivative" to "difference" and "mean" to "average"}
 
-With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general, 
+With partial cells (\forcode{ln_zps = .true.}) at bottom and top (\forcode{ln_isfcav = .true.}), in general, 
 tracers in horizontally adjacent cells live at different depths. 
 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) 
 and the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 
-The partial cell properties at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom. 
+The partial cell properties at the top (\forcode{ln_isfcav = .true.}) are computed in the same way as for the bottom. 
 So, only the bottom interpolation is explained below.
 
@@ -1413,5 +1413,5 @@
 \caption{ 	\protect\label{Fig_Partial_step_scheme} 
 Discretisation of the horizontal difference and average of tracers in the $z$-partial 
-step coordinate (\protect\np{ln\_zps}=true) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. 
+step coordinate (\protect\forcode{ln_zps = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. 
 A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value 
 at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_ZDF.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_ZDF.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_ZDF.tex	(revision 9392)
@@ -42,6 +42,6 @@
 general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 
 These trends can be computed using either a forward time stepping scheme 
-(namelist parameter \np{ln\_zdfexp}=true) or a backward time stepping 
-scheme (\np{ln\_zdfexp}=false) depending on the magnitude of the mixing 
+(namelist parameter \forcode{ln_zdfexp = .true.}) or a backward time stepping 
+scheme (\forcode{ln_zdfexp = .false.}) depending on the magnitude of the mixing 
 coefficients, and thus of the formulation used (see \S\ref{STP}).
 
@@ -65,5 +65,5 @@
 \end{align*}
 
-These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. 
+These values are set through the \np{rn_avm0} and \np{rn_avt0} namelist parameters. 
 In all cases, do not use values smaller that those associated with the molecular 
 viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, 
@@ -103,9 +103,9 @@
 is the maximum value that can be reached by the coefficient when $Ri\leq 0$, 
 $a=5$ and $n=2$. The last three values can be modified by setting the 
-\np{rn\_avmri}, \np{rn\_alp} and \np{nn\_ric} namelist parameters, respectively.
+\np{rn_avmri}, \np{rn_alp} and \np{nn_ric} namelist parameters, respectively.
 
 A simple mixing-layer model to transfer and dissipate the atmospheric
  forcings (wind-stress and buoyancy fluxes) can be activated setting 
-the \np{ln\_mldw} =.true. in the namelist.
+the \np{ln_mldw} =.true. in the namelist.
 
 In this case, the local depth of turbulent wind-mixing or "Ekman depth"
@@ -125,7 +125,7 @@
 
 is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$.
-The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}.
+The final $h_{e}$ is further constrained by the adjustable bounds \np{rn_mldmin} and \np{rn_mldmax}.
 Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to 
-the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{Lermusiaux2001}.
+the empirical values \np{rn_wtmix} and \np{rn_wvmix} \citep{Lermusiaux2001}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -170,5 +170,5 @@
 and diffusivity coefficients. The constants $C_k =  0.1$ and $C_\epsilon = \sqrt {2} /2$  
 $\approx 0.7$ are designed to deal with vertical mixing at any depth \citep{Gaspar1990}. 
-They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}. 
+They are set through namelist parameters \np{nn_ediff} and \np{nn_ediss}. 
 $P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function 
 of the local Richardson number, $R_i$:
@@ -181,8 +181,8 @@
 \end{align*}
 Options are defined through the  \ngn{namzdfy\_tke} namelist variables.
-The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable.
+The choice of $P_{rt}$ is controlled by the \np{nn_pdl} namelist variable.
 
 At the sea surface, the value of $\bar{e}$ is prescribed from the wind 
-stress field as $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn\_ebb} 
+stress field as $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn_ebb} 
 namelist parameter. The default value of $e_{bb}$ is 3.75. \citep{Gaspar1990}), 
 however a much larger value can be used when taking into account the 
@@ -191,5 +191,5 @@
 The time integration of the $\bar{e}$ equation may formally lead to negative values 
 because the numerical scheme does not ensure its positivity. To overcome this 
-problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} 
+problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn_emin} 
 namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set 
 to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations 
@@ -199,10 +199,10 @@
 instabilities associated with too weak vertical diffusion. They must be 
 specified at least larger than the molecular values, and are set through 
-\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
+\np{rn_avm0} and \np{rn_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
 
 \subsubsection{Turbulent length scale}
 For computational efficiency, the original formulation of the turbulent length 
 scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 
-are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist 
+are proposed, the choice of which is controlled by the \np{nn_mxl} namelist 
 parameter. The first two are based on the following first order approximation 
 \citep{Blanke1993}:
@@ -212,10 +212,10 @@
 which is valid in a stable stratified region with constant values of the Brunt-
 Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance 
-to the surface or to the bottom (\np{nn\_mxl} = 0) or by the local vertical scale factor 
-(\np{nn\_mxl} = 1). \citet{Blanke1993} notice that this simplification has two major 
+to the surface or to the bottom (\np{nn_mxl} = 0) or by the local vertical scale factor 
+(\np{nn_mxl} = 1). \citet{Blanke1993} notice that this simplification has two major 
 drawbacks: it makes no sense for locally unstable stratification and the 
 computation no longer uses all the information contained in the vertical density 
 profile. To overcome these drawbacks, \citet{Madec1998} introduces the 
-\np{nn\_mxl} = 2 or 3 cases, which add an extra assumption concerning the vertical 
+\np{nn_mxl} = 2 or 3 cases, which add an extra assumption concerning the vertical 
 gradient of the computed length scale. So, the length scales are first evaluated 
 as in \eqref{Eq_tke_mxl0_1} and then bounded such that:
@@ -253,7 +253,7 @@
 $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
 
-In the \np{nn\_mxl}~=~2 case, the dissipation and mixing length scales take the same 
+In the \np{nn_mxl}~=~2 case, the dissipation and mixing length scales take the same 
 value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the 
-\np{nn\_mxl}~=~3 case, the dissipation and mixing turbulent length scales are give 
+\np{nn_mxl}~=~3 case, the dissipation and mixing turbulent length scales are give 
 as in \citet{Gaspar1990}:
 \begin{equation} \label{Eq_tke_mxl_gaspar}
@@ -264,9 +264,9 @@
 \end{equation}
 
-At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} namelist 
+At the ocean surface, a non zero length scale is set through the  \np{rn_mxl0} namelist 
 parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$ 
 where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness 
 parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} 
-leads to a 0.04~m, the default value of \np{rn\_mxl0}. In the ocean interior 
+leads to a 0.04~m, the default value of \np{rn_mxl0}. In the ocean interior 
 a minimum length scale is set to recover the molecular viscosity when $\bar{e}$ 
 reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
@@ -296,8 +296,8 @@
 citing observation evidence, and $\alpha_{CB} = 100$ the Craig and Banner's value.
 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 
-with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds 
-to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc} 
+with $e_{bb}$ the \np{rn_ebb} namelist parameter, setting \np{rn_ebb}~=~67.83 corresponds 
+to $\alpha_{CB} = 100$. Further setting  \np{ln_mxl0} to true applies \eqref{ZDF_Lsbc} 
 as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value.
-Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) 
+Note that a minimal threshold of \np{rn_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) 
 is applied on surface $\bar{e}$ value.
 
@@ -317,5 +317,5 @@
 of LC in an extra source terms of TKE, $P_{LC}$.
 The presence of $P_{LC}$ in \eqref{Eq_zdftke_e}, the TKE equation, is controlled 
-by setting \np{ln\_lc} to \textit{true} in the namtke namelist.
+by setting \np{ln_lc} to \textit{true} in the namtke namelist.
  
 By making an analogy with the characteristic convective velocity scale 
@@ -343,5 +343,5 @@
 where $c_{LC} = 0.15$ has been chosen by \citep{Axell_JGR02} as a good compromise 
 to fit LES data. The chosen value yields maximum vertical velocities $w_{LC}$ of the order 
-of a few centimeters per second. The value of $c_{LC}$ is set through the \np{rn\_lc} 
+of a few centimeters per second. The value of $c_{LC}$ is set through the \np{rn_lc} 
 namelist parameter, having in mind that it should stay between 0.15 and 0.54 \citep{Axell_JGR02}. 
 
@@ -366,5 +366,5 @@
 ($i.e.$ near-inertial oscillations and ocean swells and waves).
 
-When using this parameterization ($i.e.$ when \np{nn\_etau}~=~1), the TKE input to the ocean ($S$) 
+When using this parameterization ($i.e.$ when \np{nn_etau}~=~1), the TKE input to the ocean ($S$) 
 imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized 
 by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by:
@@ -379,10 +379,10 @@
 and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely 
 covered by sea-ice).
-The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. 
-The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}~=~0) 
+The value of $f_r$, usually a few percents, is specified through \np{rn_efr} namelist parameter. 
+The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn_etau}~=~0) 
 or a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m 
-at high latitudes (\np{nn\_etau}~=~1). 
-
-Note that two other option existe, \np{nn\_etau}~=~2, or 3. They correspond to applying 
+at high latitudes (\np{nn_etau}~=~1). 
+
+Note that two other option existe, \np{nn_etau}~=~2, or 3. They correspond to applying 
 \eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part 
 of the stress to evaluate the fraction of TKE that penetrate the ocean. 
@@ -558,5 +558,5 @@
 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) 
 depends of the choice of the turbulence model. Four different turbulent models are pre-defined 
-(Tab.\ref{Tab_GLS}). They are made available through the \np{nn\_clo} namelist parameter. 
+(Tab.\ref{Tab_GLS}). They are made available through the \np{nn_clo} namelist parameter. 
 
 %--------------------------------------------------TABLE--------------------------------------------------
@@ -567,5 +567,5 @@
 %                        & \citep{Mellor_Yamada_1982} &  \citep{Rodi_1987}       & \citep{Wilcox_1988} &                 \\  
 \hline  \hline 
-\np{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\  
+\np{nn_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\  
 \hline 
 $( p , n , m )$	       &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
@@ -581,5 +581,5 @@
 \caption{   \protect\label{Tab_GLS} 
 Set of predefined GLS parameters, or equivalently predefined turbulence models available 
-with \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls} .}
+with \protect\key{zdfgls} and controlled by the \protect\np{nn_clos} namelist variable in \protect\ngn{namzdf\_gls} .}
 \end{center}	\end{table}
 %--------------------------------------------------------------------------------------------------------------
@@ -589,12 +589,12 @@
 value near physical boundaries (logarithmic boundary layer law). $C_{\mu}$ and $C_{\mu'}$ 
 are calculated from stability function proposed by \citet{Galperin_al_JAS88}, or by \citet{Kantha_Clayson_1994} 
-or one of the two functions suggested by \citet{Canuto_2001}  (\np{nn\_stab\_func} = 0, 1, 2 or 3, resp.). 
+or one of the two functions suggested by \citet{Canuto_2001}  (\np{nn_stab_func} = 0, 1, 2 or 3, resp.). 
 The value of $C_{0\mu}$ depends of the choice of the stability function.
 
 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated 
-thanks to Dirichlet or Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp. 
-As for TKE closure , the wave effect on the mixing is considered when \np{ln\_crban}~=~true
-\citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}. The \np{rn\_crban} namelist parameter 
-is $\alpha_{CB}$ in \eqref{ZDF_Esbc} and \np{rn\_charn} provides the value of $\beta$ in \eqref{ZDF_Lsbc}. 
+thanks to Dirichlet or Neumann condition through \np{nn_tkebc_surf} and \np{nn_tkebc_bot}, resp. 
+As for TKE closure , the wave effect on the mixing is considered when \np{ln_crban}~=~true
+\citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}. The \np{rn_crban} namelist parameter 
+is $\alpha_{CB}$ in \eqref{ZDF_Esbc} and \np{rn_charn} provides the value of $\beta$ in \eqref{ZDF_Lsbc}. 
 
 The $\psi$ equation is known to fail in stably stratified flows, and for this reason 
@@ -606,5 +606,5 @@
 stably stratified situations, and that its value has to be chosen in accordance 
 with the algebraic model for the turbulent fluxes. The clipping is only activated 
-if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value.
+if \forcode{ln_length_lim = .true.}, and the $c_{lim}$ is set to the \np{rn_clim_galp} value.
 
 The time and space discretization of the GLS equations follows the same energetic 
@@ -646,6 +646,6 @@
 %       Non-Penetrative Convective Adjustment 
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Non-Penetrative Convective Adjustment (\protect\np{ln\_tranpc}) ]
-			{Non-Penetrative Convective Adjustment (\protect\np{ln\_tranpc}=.true.) }
+\subsection   [Non-Penetrative Convective Adjustment (\protect\np{ln_tranpc}) ]
+			{Non-Penetrative Convective Adjustment (\protect\np{ln_tranpc}=.true.) }
 \label{ZDF_npc}
 
@@ -671,6 +671,6 @@
 
 Options are defined through the  \ngn{namzdf} namelist variables.
-The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}~=~\textit{true}. 
-It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously 
+The non-penetrative convective adjustment is used when \np{ln_zdfnpc}~=~\textit{true}. 
+It is applied at each \np{nn_npc} time step and mixes downwards instantaneously 
 the statically unstable portion of the water column, but only until the density 
 structure becomes neutrally stable ($i.e.$ until the mixed portion of the water 
@@ -713,6 +713,6 @@
 %       Enhanced Vertical Diffusion 
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Enhanced Vertical Diffusion (\protect\np{ln\_zdfevd})]
-		     	  {Enhanced Vertical Diffusion (\protect\np{ln\_zdfevd}=true)}
+\subsection   [Enhanced Vertical Diffusion (\protect\np{ln_zdfevd})]
+		     	  {Enhanced Vertical Diffusion (\protect\forcode{ln_zdfevd = .true.})}
 \label{ZDF_evd}
 
@@ -722,19 +722,19 @@
 
 Options are defined through the  \ngn{namzdf} namelist variables.
-The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}=true. 
+The enhanced vertical diffusion parameterisation is used when \forcode{ln_zdfevd = .true.}. 
 In this case, the vertical eddy mixing coefficients are assigned very large values 
 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 
 ($i.e.$ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) 
 \citep{Lazar_PhD97, Lazar_al_JPO99}. This is done either on tracers only 
-(\np{nn\_evdm}=0) or on both momentum and tracers (\np{nn\_evdm}=1).
+(\forcode{nn_evdm = 0}) or on both momentum and tracers (\forcode{nn_evdm = 1}).
 
 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and 
-if \np{nn\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
-values also, are set equal to the namelist parameter \np{rn\_avevd}. A typical value 
+if \forcode{nn_evdm = 1}, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
+values also, are set equal to the namelist parameter \np{rn_avevd}. A typical value 
 for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of 
 convective processes is less time consuming than the convective adjustment 
 algorithm presented above when mixing both tracers and momentum in the 
 case of static instabilities. It requires the use of an implicit time stepping on 
-vertical diffusion terms (i.e. \np{ln\_zdfexp}=false). 
+vertical diffusion terms (i.e. \forcode{ln_zdfexp = .false.}). 
 
 Note that the stability test is performed on both \textit{before} and \textit{now} 
@@ -761,10 +761,10 @@
 because the mixing length scale is bounded by the distance to the sea surface. 
 It can thus be useful to combine the enhanced vertical 
-diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc} 
+diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln_zdfnpc} 
 namelist parameter to true and defining the turbulent closure CPP key all together.
 
 The KPP turbulent closure scheme already includes enhanced vertical diffusion 
 in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ 
-found in \mdl{zdfkpp}, therefore \np{ln\_zdfevd}=false should be used with the KPP 
+found in \mdl{zdfkpp}, therefore \forcode{ln_zdfevd = .false.} should be used with the KPP 
 scheme. %gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
 
@@ -918,5 +918,5 @@
 %       Linear Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Linear Bottom Friction (\protect\np{nn\_botfr} = 0 or 1) }
+\subsection{Linear Bottom Friction (\protect\np{nn_botfr} = 0 or 1) }
 \label{ZDF_bfr_linear}
 
@@ -940,5 +940,5 @@
 $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 
 This is the default value used in \NEMO. It corresponds to a decay time scale 
-of 115~days. It can be changed by specifying \np{rn\_bfri1} (namelist parameter).
+of 115~days. It can be changed by specifying \np{rn_bfri1} (namelist parameter).
 
 For the linear friction case the coefficients defined in the general 
@@ -950,17 +950,17 @@
 \end{split}
 \end{equation}
-When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfri1}. 
-Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip 
+When \forcode{nn_botfr = 1}, the value of $r$ used is \np{rn_bfri1}. 
+Setting \forcode{nn_botfr = 0} is equivalent to setting $r=0$ and leads to a free-slip 
 bottom boundary condition. These values are assigned in \mdl{zdfbfr}. 
 From v3.2 onwards there is support for local enhancement of these values 
-via an externally defined 2D mask array (\np{ln\_bfr2d}=true) given
+via an externally defined 2D mask array (\forcode{ln_bfr2d = .true.}) given
 in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1. 
 Locations with a non-zero mask value will have the friction coefficient increased 
-by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}.
+by $mask\_value$*\np{rn_bfrien}*\np{rn_bfri1}.
 
 % -------------------------------------------------------------------------------------------------------------
 %       Non-Linear Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Non-Linear Bottom Friction (\protect\np{nn\_botfr} = 2)}
+\subsection{Non-Linear Bottom Friction (\protect\np{nn_botfr} = 2)}
 \label{ZDF_bfr_nonlinear}
 
@@ -977,5 +977,5 @@
 $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} 
 uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. 
-The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} 
+The CME choices have been set as default values (\np{rn_bfri2} and \np{rn_bfeb2} 
 namelist parameters).
 
@@ -993,19 +993,19 @@
 
 The coefficients that control the strength of the non-linear bottom friction are
-initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}.
+initialised as namelist parameters: $C_D$= \np{rn_bfri2}, and $e_b$ =\np{rn_bfeb2}.
 Note for applications which treat tides explicitly a low or even zero value of
-\np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ is possible
-via an externally defined 2D mask array (\np{ln\_bfr2d}=true).  This works in the same way
+\np{rn_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ is possible
+via an externally defined 2D mask array (\forcode{ln_bfr2d = .true.}).  This works in the same way
 as for the linear bottom friction case with non-zero masked locations increased by
-$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}.
+$mask\_value$*\np{rn_bfrien}*\np{rn_bfri2}.
 
 % -------------------------------------------------------------------------------------------------------------
 %       Bottom Friction Log-layer
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Log-layer Bottom Friction enhancement (\protect\np{ln\_loglayer} = .true.)]{Log-layer Bottom Friction enhancement (\protect\np{nn\_botfr} = 2, \protect\np{ln\_loglayer} = .true.)}
+\subsection[Log-layer Bottom Friction enhancement (\protect\np{ln_loglayer} = .true.)]{Log-layer Bottom Friction enhancement (\protect\np{nn_botfr} = 2, \protect\np{ln_loglayer} = .true.)}
 \label{ZDF_bfr_loglayer}
 
 In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally
-enhanced using a "law of the wall" scaling. If  \np{ln\_loglayer} = .true., $C_D$ is no
+enhanced using a "law of the wall" scaling. If  \np{ln_loglayer} = .true., $C_D$ is no
 longer constant but is related to the thickness of the last wet layer in each column by:
 
@@ -1014,10 +1014,10 @@
 \end{equation}
 
-\noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness
+\noindent where $\kappa$ is the von-Karman constant and \np{rn_bfrz0} is a roughness
 length provided via the namelist.
 
 For stability, the drag coefficient is bounded such that it is kept greater or equal to
-the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional
-namelist parameter: \np{rn\_bfri2\_max}, i.e.:
+the base \np{rn_bfri2} value and it is not allowed to exceed the value of an additional
+namelist parameter: \np{rn_bfri2_max}, i.e.:
 
 \begin{equation}
@@ -1026,7 +1026,7 @@
 
 \noindent Note also that a log-layer enhancement can also be applied to the top boundary
-friction if under ice-shelf cavities are in use (\np{ln\_isfcav}=.true.).  In this case, the
-relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2}
-and \np{rn\_tfri2\_max}.
+friction if under ice-shelf cavities are in use (\np{ln_isfcav}=.true.).  In this case, the
+relevant namelist parameters are \np{rn_tfrz0}, \np{rn_tfri2}
+and \np{rn_tfri2_max}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1082,11 +1082,11 @@
 %       Implicit Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Implicit Bottom Friction (\protect\np{ln\_bfrimp})]{Implicit Bottom Friction (\protect\np{ln\_bfrimp}$=$\textit{T})}
+\subsection[Implicit Bottom Friction (\protect\np{ln_bfrimp})]{Implicit Bottom Friction (\protect\np{ln_bfrimp}$=$\textit{T})}
 \label{ZDF_bfr_imp}
 
 An optional implicit form of bottom friction has been implemented to improve
 model stability. We recommend this option for shelf sea and coastal ocean applications, especially 
-for split-explicit time splitting. This option can be invoked by setting \np{ln\_bfrimp} 
-to \textit{true} in the \textit{nambfr} namelist. This option requires \np{ln\_zdfexp} to be \textit{false} 
+for split-explicit time splitting. This option can be invoked by setting \np{ln_bfrimp} 
+to \textit{true} in the \textit{nambfr} namelist. This option requires \np{ln_zdfexp} to be \textit{false} 
 in the \textit{namzdf} namelist. 
 
@@ -1135,5 +1135,5 @@
 %       Bottom Friction with split-explicit time splitting
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Bottom Friction with split-explicit time splitting]{Bottom Friction with split-explicit time splitting (\protect\np{ln\_bfrimp})}
+\subsection[Bottom Friction with split-explicit time splitting]{Bottom Friction with split-explicit time splitting (\protect\np{ln_bfrimp})}
 \label{ZDF_bfr_ts}
 
@@ -1144,7 +1144,7 @@
 \key{dynspg\_flt}). Extra attention is required, however, when using 
 split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface 
-equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, while the three 
+equation is solved with a small time step \np{rn_rdt}/\np{nn_baro}, while the three 
 dimensional prognostic variables are solved with the longer time step 
-of \np{rn\_rdt} seconds. The trend in the barotropic momentum due to bottom 
+of \np{rn_rdt} seconds. The trend in the barotropic momentum due to bottom 
 friction appropriate to this method is that given by the selected parameterisation 
 ($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities 
@@ -1176,5 +1176,5 @@
 limiting is thought to be having a major effect (a more likely prospect in coastal and shelf seas
 applications) then the fully implicit form of the bottom friction should be used (see \S\ref{ZDF_bfr_imp} ) 
-which can be selected by setting \np{ln\_bfrimp} $=$ \textit{true}.
+which can be selected by setting \np{ln_bfrimp} $=$ \textit{true}.
 
 Otherwise, the implicit formulation takes the form:
@@ -1220,7 +1220,7 @@
 and $F(z)$ the vertical structure function. 
 
-The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter)
+The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn_me} namelist parameter)
 and is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980). 
-The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter) 
+The tidal dissipation efficiency is given by the parameter $q$ (\np{rn_tfe} namelist parameter) 
 represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally, 
 with the remaining $1-q$ radiating away as low mode internal waves and 
@@ -1229,5 +1229,5 @@
 The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical. 
 It is implemented as a simple exponential decaying upward away from the bottom, 
-with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter, with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04}, 
+with a vertical scale of $h_o$ (\np{rn_htmx} namelist parameter, with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04}, 
 \begin{equation} \label{Eq_Fz}
 F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }
@@ -1238,5 +1238,5 @@
 diffusivity assuming a Prandtl number of 1, $i.e.$ $A^{vm}_{tides}=A^{vT}_{tides}$. 
 In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity 
-is capped at $300\,cm^2/s$ and impose a lower limit on $N^2$ of \np{rn\_n2min} 
+is capped at $300\,cm^2/s$ and impose a lower limit on $N^2$ of \np{rn_n2min} 
 usually set to $10^{-8} s^{-2}$. These bounds are usually rarely encountered.
 
@@ -1266,14 +1266,14 @@
 %        Indonesian area specific treatment 
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})}
+\subsection{Indonesian area specific treatment (\protect\np{ln_zdftmx_itf})}
 \label{ZDF_tmx_itf}
 
 When the Indonesian Through Flow (ITF) area is included in the model domain,
 a specific treatment of tidal induced mixing in this area can be used. 
-It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide
+It is activated through the namelist logical \np{ln_tmx_itf}, and the user must provide
 an input NetCDF file, \ifile{mask\_itf}, which contains a mask array defining the ITF area
 where the specific treatment is applied.
 
-When \np{ln\_tmx\_itf}=true, the two key parameters $q$ and $F(z)$ are adjusted following 
+When \forcode{ln_tmx_itf = .true.}, the two key parameters $q$ and $F(z)$ are adjusted following 
 the parameterisation developed by \citet{Koch-Larrouy_al_GRL07}:
 
@@ -1285,5 +1285,5 @@
 So it is assumed that $q = 1$, $i.e.$ all the energy generated is available for mixing.
 Note that for test purposed, the ITF tidal dissipation efficiency is a 
-namelist parameter (\np{rn\_tfe\_itf}). A value of $1$ or close to is
+namelist parameter (\np{rn_tfe_itf}). A value of $1$ or close to is
 this recommended for this parameter.
 
@@ -1329,5 +1329,5 @@
 \end{equation}
 where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution 
-of the energy available for mixing. If the \np{ln\_mevar} namelist parameter is set to false, 
+of the energy available for mixing. If the \np{ln_mevar} namelist parameter is set to false, 
 the mixing efficiency is taken as constant and equal to 1/6 \citep{Osborn_JPO80}. 
 In the opposite (recommended) case, $R_f$ is instead a function of the turbulence intensity parameter 
@@ -1338,5 +1338,5 @@
 
 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 
-as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice). 
+as a function of $Re_b$ by setting the \np{ln_tsdiff} parameter to true, a recommended choice). 
 This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014}, 
 is implemented as in \cite{de_lavergne_JPO2016_efficiency}.
@@ -1356,5 +1356,5 @@
 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
 \end{equation*}
-The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)  controls the stratification-dependence of the pycnocline-intensified dissipation. 
+The $n_p$ parameter (given by \np{nn_zpyc} in \ngn{namzdf\_tmx\_new} namelist)  controls the stratification-dependence of the pycnocline-intensified dissipation. 
 It can take values of 1 (recommended) or 2.
 Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_misc.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_misc.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_misc.tex	(revision 9392)
@@ -100,16 +100,14 @@
 existing 'zoom' options are overly complex for this task and marked for deletion anyway.
 This alternative subsetting operates for the j-direction only and works by optionally
-looking for and using a global file attribute (named: \np{open\_ocean\_jstart}) to
+looking for and using a global file attribute (named: \np{open_ocean_jstart}) to
 determine the starting j-row for input. The use of this option is best explained with an
 example: Consider an ORCA1 configuration using the extended grid bathymetry and coordinate
 files:
 \vspace{-10pt}
-\begin{verbatim}
-eORCA1_bathymetry_v2.nc
-eORCA1_coordinates.nc
-\end{verbatim}
+\ifile{eORCA1\_bathymetry\_v2}
+\ifile{eORCA1\_coordinates}
 \noindent These files define a horizontal domain of 362x332. Assuming the first row with
 open ocean wet points in the non-isf bathymetry for this set is row 42 (Fortran indexing)
-then the formally correct setting for \np{open\_ocean\_jstart} is 41. Using this value as the
+then the formally correct setting for \np{open_ocean_jstart} is 41. Using this value as the
 first row to be read will result in a 362x292 domain which is the same size as the original
 ORCA1 domain. Thus the extended coordinates and bathymetry files can be used with all the
@@ -130,5 +128,5 @@
 
 \noindent Note the j-size of the global domain is the (extended j-size minus
-\np{open\_ocean\_jstart} + 1 ) and this must match the size of all datasets other than
+\np{open_ocean_jstart} + 1 ) and this must match the size of all datasets other than
 bathymetry and coordinates currently. However the option can be extended to any global, 2D
 and 3D, netcdf, input field by adding the:
@@ -137,5 +135,5 @@
 lrowattr=ln_use_jattr
 \end{forlines}
-optional argument to the appropriate \np{iom\_get} call and the \np{open\_ocean\_jstart} attribute to the corresponding input files. It remains the users responsibility to set \np{jpjdta} and \np{jpjglo} values in the \np{namelist\_cfg} file according to their needs.
+optional argument to the appropriate \np{iom_get} call and the \np{open_ocean_jstart} attribute to the corresponding input files. It remains the users responsibility to set \np{jpjdta} and \np{jpjglo} values in the \np{namelist_cfg} file according to their needs.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -225,5 +223,5 @@
 required by each individual for the fold operation. This alternative method should give identical
 results to the default \textsc{ALLGATHER} method and is recommended for large values of \np{jpni}.
-The new method is activated by setting \np{ln\_nnogather} to be true ({\bf nammpp}). The
+The new method is activated by setting \np{ln_nnogather} to be true ({\bf nammpp}). The
 reproducibility of results using the two methods should be confirmed for each new, non-reference
 configuration.
@@ -253,13 +251,13 @@
 $\bullet$ Control print %: describe here 4 things:
 
-1- \np{ln\_ctl} : compute and print the trends averaged over the interior domain 
+1- \np{ln_ctl} : compute and print the trends averaged over the interior domain 
 in all TRA, DYN, LDF and ZDF modules. This option is very helpful when 
 diagnosing the origin of an undesired change in model results. 
 
-2- also \np{ln\_ctl} but using the nictl and njctl namelist parameters to check 
+2- also \np{ln_ctl} but using the nictl and njctl namelist parameters to check 
 the source of differences between mono and multi processor runs.
 
 %%gm   to be removed both here and in the code
-3- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of 
+3- last digit comparison (\np{nn_bit_cmp}). In an MPP simulation, the computation of 
 a sum over the whole domain is performed as the summation over all processors of 
 each of their sums over their interior domains. This double sum never gives exactly 
@@ -271,5 +269,5 @@
 %%gm end
 
-$\bullet$  Benchmark (\np{nn\_bench}). This option defines a benchmark run based on 
+$\bullet$  Benchmark (\np{nn_bench}). This option defines a benchmark run based on 
 a GYRE configuration (see \S\ref{CFG_gyre}) in which the resolution remains the same 
 whatever the domain size. This allows a very large model domain to be used, just by 
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_model_basics_zstar.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_model_basics_zstar.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_model_basics_zstar.tex	(revision 9392)
@@ -247,5 +247,5 @@
 documented in \S\ref{MISC}. The amplitude of the extra term is given by the namelist variable \np{rnu}. The default value is 1, as recommended by \citet{Roullet2000}
 
-\colorbox{red}{\np{rnu}=1 to be suppressed from namelist !}
+\colorbox{red}{\forcode{rnu = 1} to be suppressed from namelist !}
 
 %-------------------------------------------------------------
Index: branches/2017/dev_merge_2017/DOC/tex_sub/chap_time_domain.tex
===================================================================
--- branches/2017/dev_merge_2017/DOC/tex_sub/chap_time_domain.tex	(revision 9389)
+++ branches/2017/dev_merge_2017/DOC/tex_sub/chap_time_domain.tex	(revision 9392)
@@ -91,6 +91,6 @@
 \end{equation} 
 where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin 
-coefficient. $\gamma$ is initialized as \np{rn\_atfp} (namelist parameter). 
-Its default value is \np{rn\_atfp}=$10^{-3}$ (see \S~\ref{STP_mLF}), 
+coefficient. $\gamma$ is initialized as \np{rn_atfp} (namelist parameter). 
+Its default value is \np{rn_atfp}=$10^{-3}$ (see \S~\ref{STP_mLF}), 
 causing only a weak dissipation of high frequency motions (\citep{Farge1987}). 
 The addition of a time filter degrades the accuracy of the 
@@ -141,6 +141,6 @@
 constraint on the time step. Two solutions are available in \NEMO to overcome 
 the stability constraint: $(a)$ a forward time differencing scheme using a 
-time splitting technique (\np{ln\_zdfexp} = true) or $(b)$ a backward (or implicit) 
-time differencing scheme (\np{ln\_zdfexp} = false). In $(a)$, the master 
+time splitting technique (\np{ln_zdfexp} = true) or $(b)$ a backward (or implicit) 
+time differencing scheme (\np{ln_zdfexp} = false). In $(a)$, the master 
 time step $\Delta $t is cut into $N$ fractional time steps so that the 
 stability criterion is reduced by a factor of $N$. The computation is performed as 
@@ -156,5 +156,5 @@
 \end{equation}
 with DF a vertical diffusion term. The number of fractional time steps, $N$, is given 
-by setting \np{nn\_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally 
+by setting \np{nn_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally 
 stable but diffusive. It can be written as follows:
 \begin{equation} \label{Eq_STP_imp}
@@ -309,5 +309,5 @@
 gradient (see \S\ref{DYN_hpg_imp}), an extra three-dimensional field has to be 
 added to the restart file to ensure an exact restartability. This is done optionally 
-via the  \np{nn\_dynhpg\_rst} namelist parameter, so that the size of the
+via the  \np{nn_dynhpg_rst} namelist parameter, so that the size of the
 restart file can be reduced when restartability is not a key issue (operational 
 oceanography or in ensemble simulations for seasonal forecasting).
