Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -55,9 +57,9 @@
 
 \begin{description}
- \item[\np{jphgr_mesh}{jphgr\_mesh}=0]  The most general curvilinear orthogonal grids.
+ \item[{\np{jphgr_mesh}{jphgr\_mesh}=0}]  The most general curvilinear orthogonal grids.
   The coordinates and their first derivatives with respect to $i$ and $j$ are provided
   in a input file (\ifile{coordinates}), read in \rou{hgr\_read} subroutine of the domhgr module.
   This is now the only option available within \NEMO\ itself from v4.0 onwards.
-\item[\np{jphgr_mesh}{jphgr\_mesh}=1 to 5] A few simple analytical grids are provided (see below).
+\item[{\np{jphgr_mesh}{jphgr\_mesh}=1 to 5}] A few simple analytical grids are provided (see below).
   For other analytical grids, the \mdl{domhgr} module (\texttt{DOMAINcfg} variant) must be
   modified by the user. In most cases, modifying the \mdl{usrdef\_hgr} module of \NEMO\ is
@@ -136,5 +138,5 @@
 
 It is possible to define a simple regular vertical grid by giving zero stretching
-(\np{ppacr}{ppacr}\forcode{ = 0}).  In that case, the parameters \jp{jpk} (number of $w$-levels)
+(\np[=0]{ppacr}{ppacr}).  In that case, the parameters \jp{jpk} (number of $w$-levels)
 and \np{pphmax}{pphmax} (total ocean depth in meters) fully define the grid.
 
@@ -152,5 +154,5 @@
 top and bottom with a smooth hyperbolic tangent transition in between (\autoref{fig:DOMCFG_zgr}).
 
-A double hyperbolic tangent version (\np{ldbletanh}{ldbletanh}\forcode{ = .true.}) is also available
+A double hyperbolic tangent version (\np[=.true.]{ldbletanh}{ldbletanh}) is also available
 which permits finer control and is used, typically, to obtain a well resolved upper ocean
 without compromising on resolution at depth using a moderate number of levels.
@@ -170,5 +172,5 @@
 \end{gather}
 
-If the ice shelf cavities are opened (\np{ln_isfcav}{ln\_isfcav}\forcode{ = .true.}), the definition
+If the ice shelf cavities are opened (\np[=.true.]{ln_isfcav}{ln\_isfcav}), the definition
 of $z_0$ is the same.  However, definition of $e_3^0$ at $t$- and $w$-points is
 respectively changed to:
@@ -314,12 +316,12 @@
 \np{nn_bathy}{nn\_bathy} (found in \nam{dom}{dom} namelist (\texttt{DOMAINCFG} variant) ):
 \begin{description}
-\item[\np{nn_bathy}{nn\_bathy}\forcode{ = 0}]:
+\item[{\np[=0]{nn_bathy}{nn\_bathy}}]:
   a flat-bottom domain is defined.
   The total depth $z_w (jpk)$ is given by the coordinate transformation.
   The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}{jperio}.
-\item[\np{nn_bathy}{nn\_bathy}\forcode{ = -1}]:
+\item[{\np[=-1]{nn_bathy}{nn\_bathy}}]:
   a domain with a bump of topography one third of the domain width at the central latitude.
   This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount.
-\item[\np{nn_bathy}{nn\_bathy}\forcode{ = 1}]:
+\item[{\np[=1]{nn_bathy}{nn\_bathy}}]:
   read a bathymetry and ice shelf draft (if needed).
   The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at
@@ -332,5 +334,5 @@
   The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at
   each grid point of the model grid.
-  This file is only needed if \np{ln_isfcav}{ln\_isfcav}\forcode{ = .true.}.
+  This file is only needed if \np[=.true.]{ln_isfcav}{ln\_isfcav}.
   Defining the ice shelf draft will also define the ice shelf edge and the grounding line position.
 \end{description}
@@ -412,5 +414,5 @@
 %--------------------------------------------------------------------------------------------------------------
 Options are defined in \nam{zgr_sco}{zgr\_sco} (\texttt{DOMAINcfg} only).
-In $s$-coordinate (\np{ln_sco}{ln\_sco}\forcode{ = .true.}), the depth and thickness of the model levels are defined from
+In $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}), the depth and thickness of the model levels are defined from
 the product of a depth field and either a stretching function or its derivative, respectively:
 
@@ -435,5 +437,5 @@
 
 The original default \NEMO\ s-coordinate stretching is available if neither of the other options are specified as true
-(\np{ln_s_SH94}{ln\_s\_SH94}\forcode{ = .false.} and \np{ln_s_SF12}{ln\_s\_SF12}\forcode{ = .false.}).
+(\np[=.false.]{ln_s_SH94}{ln\_s\_SH94} and \np[=.false.]{ln_s_SF12}{ln\_s\_SF12}).
 This uses a depth independent $\tanh$ function for the stretching \citep{madec.delecluse.ea_JPO96}:
 
@@ -455,5 +457,5 @@
 
 A stretching function,
-modified from the commonly used \citet{song.haidvogel_JCP94} stretching (\np{ln_s_SH94}{ln\_s\_SH94}\forcode{ = .true.}),
+modified from the commonly used \citet{song.haidvogel_JCP94} stretching (\np[=.true.]{ln_s_SH94}{ln\_s\_SH94}),
 is also available and is more commonly used for shelf seas modelling:
 
@@ -558,7 +560,7 @@
 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO\ web site.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -17,5 +19,5 @@
 %        UBS scheme
 % -------------------------------------------------------------------------------------------------------------
-\section{Upstream Biased Scheme (UBS) (\protect\np{ln_traadv_ubs}{ln\_traadv\_ubs}\forcode{ = .true.})}
+\section{Upstream Biased Scheme (UBS) (\protect\np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs})}
 \label{sec:ALGOS_tra_adv_ubs}
 
@@ -59,5 +61,5 @@
 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}{ln\_traadv\_ubs}\forcode{ = .true.}.
+\np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs}.
 
 For stability reasons, in \autoref{eq:TRA_adv_ubs}, the first term which corresponds to
@@ -841,7 +843,7 @@
 \ie\ the variance of the tracer is preserved by the discretisation of the skew fluxes.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -411,7 +413,7 @@
 that is a Laplacian diffusion is applied on momentum along the coordinate directions.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -366,5 +368,5 @@
 %       Vorticity Term with ENE scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity term with ENE scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene}\forcode{ = .true.})}
+\subsubsection{Vorticity term with ENE scheme (\protect\np[=.true.]{ln_dynvor_ene}{ln\_dynvor\_ene})}
 \label{subsec:INVARIANTS_vorENE}
 
@@ -406,5 +408,5 @@
 %       Vorticity Term with EEN scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity term with EEN scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een}\forcode{ = .true.})}
+\subsubsection{Vorticity term with EEN scheme (\protect\np[=.true.]{ln_dynvor_een}{ln\_dynvor\_een})}
 \label{subsec:INVARIANTS_vorEEN_vect}
 
@@ -878,5 +880,5 @@
 %       Vorticity Term with ENS scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity term with ENS scheme  (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens}\forcode{ = .true.})}
+\subsubsection{Vorticity term with ENS scheme  (\protect\np[=.true.]{ln_dynvor_ens}{ln\_dynvor\_ens})}
 \label{subsec:INVARIANTS_vorENS}
 
@@ -947,5 +949,5 @@
 %       Vorticity Term with EEN scheme
 % -------------------------------------------------------------------------------------------------------------
-\subsubsection{Vorticity Term with EEN scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een}\forcode{ = .true.})}
+\subsubsection{Vorticity Term with EEN scheme (\protect\np[=.true.]{ln_dynvor_een}{ln\_dynvor\_een})}
 \label{subsec:INVARIANTS_vorEEN}
 
@@ -1528,7 +1530,7 @@
 %%%%  end of appendix in gm comment
 %}
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -594,7 +596,7 @@
 the expression of the 3D divergence in the $s-$coordinates established above.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex	(revision 11582)
@@ -11,4 +11,6 @@
 \newcommand{\rtriadt}[1]{\ensuremath{\triadt{i}{k}{#1}{i_p}{k_p}}}
 
+\onlyinsubfile{\makeindex}
+
 \begin{document}
 % ================================================================
@@ -42,5 +44,5 @@
 The options specific to the Griffies scheme include:
 \begin{description}
-\item[\np{ln_triad_iso}{ln\_triad\_iso}]
+\item[{\np{ln_triad_iso}{ln\_triad\_iso}}]
   See \autoref{sec:TRIADS_taper}.
   If this is set false (the default),
@@ -53,5 +55,5 @@
   giving an almost pure horizontal diffusive tracer flux within the mixed layer.
   This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:TRIADS_Gerdes-taper}
-\item[\np{ln_botmix_triad}{ln\_botmix\_triad}]
+\item[{\np{ln_botmix_triad}{ln\_botmix\_triad}}]
   See \autoref{sec:TRIADS_iso_bdry}.
   If this is set false (the default) then the lateral diffusive fluxes
@@ -59,11 +61,11 @@
   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}{rn\_sw\_triad}]
+\item[{\np{rn_sw_triad}{rn\_sw\_triad}}]
   blah blah to be added....
 \end{description}
 The options shared with the Standard scheme include:
 \begin{description}
-\item[\np{ln_traldf_msc}{ln\_traldf\_msc}]   blah blah to be added
-\item[\np{rn_slpmax}{rn\_slpmax}]  blah blah to be added
+\item[{\np{ln_traldf_msc}{ln\_traldf\_msc}}]   blah blah to be added
+\item[{\np{rn_slpmax}{rn\_slpmax}}]  blah blah to be added
 \end{description}
 
@@ -646,10 +648,10 @@
 Note that both near bottom triad slopes \triad{i}{k}{R}{1/2}{1/2} and \triad{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, \ie\ the $i,k+1$ $u$-point is masked.
-The associated lateral fluxes (grey-black dashed line) are masked if \np{ln_botmix_triad}{ln\_botmix\_triad}\forcode{ = .false.},
-but left unmasked, giving bottom mixing, if \np{ln_botmix_triad}{ln\_botmix\_triad}\forcode{ = .true.}.
-
-The default option \np{ln_botmix_triad}{ln\_botmix\_triad}\forcode{ = .false.} is suitable when the bbl mixing option is enabled
-(\np{ln_trabbl}{ln\_trabbl}\forcode{ = .true.}, with \np{nn_bbl_ldf}{nn\_bbl\_ldf}\forcode{ = 1}), or for simple idealized problems.
-For setups with topography without bbl mixing, \np{ln_botmix_triad}{ln\_botmix\_triad}\forcode{ = .true.} may be necessary.
+The associated lateral fluxes (grey-black dashed line) are masked if \np[=.false.]{ln_botmix_triad}{ln\_botmix\_triad},
+but left unmasked, giving bottom mixing, if \np[=.true.]{ln_botmix_triad}{ln\_botmix\_triad}.
+
+The default option \np[=.false.]{ln_botmix_triad}{ln\_botmix\_triad} is suitable when the bbl mixing option is enabled
+(\np[=.true.]{ln_trabbl}{ln\_trabbl}, with \np[=1]{nn_bbl_ldf}{nn\_bbl\_ldf}), or for simple idealized problems.
+For setups with topography without bbl mixing, \np[=.true.]{ln_botmix_triad}{ln\_botmix\_triad} may be necessary.
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[h]
@@ -672,6 +674,6 @@
     \ie\ the $i,k+1$ $u$-point is masked.
     The associated lateral fluxes (grey-black dashed line) are masked if
-    \protect\np{ln_botmix_triad}{ln\_botmix\_triad}\forcode{ = .false.}, but left unmasked,
-    giving bottom mixing, if \protect\np{ln_botmix_triad}{ln\_botmix\_triad}\forcode{ = .true.}}
+    \protect\np[=.false.]{ln_botmix_triad}{ln\_botmix\_triad}, but left unmasked,
+    giving bottom mixing, if \protect\np[=.true.]{ln_botmix_triad}{ln\_botmix\_triad}}
   \label{fig:TRIADS_bdry_triads}
 \end{figure}
@@ -715,5 +717,5 @@
 \label{sec:TRIADS_lintaper}
 
-This is the option activated by the default choice \np{ln_triad_iso}{ln\_triad\_iso}\forcode{ = .false.}.
+This is the option activated by the default choice \np[=.false.]{ln_triad_iso}{ln\_triad\_iso}.
 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 \autoref{fig:LDF_eiv_slp}, to values
@@ -1145,5 +1147,5 @@
 \label{sec:TRIADS_sfdiag}
 
-Where the namelist parameter \np{ln_traldf_gdia}{ln\_traldf\_gdia}\forcode{ = .true.},
+Where the namelist parameter \np[=.true.]{ln_traldf_gdia}{ln\_traldf\_gdia},
 diagnosed mean eddy-induced velocities are output.
 Each time step, streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
@@ -1171,7 +1173,7 @@
 \]
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -194,7 +196,7 @@
 \end{clines}
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -125,5 +127,5 @@
 
 XIOS may be used to read single file restart produced by \NEMO. Currently only the variables written to
-file \forcode{numror} can be handled by XIOS. To activate restart reading using XIOS, set \np{ln_xios_read}{ln\_xios\_read}\forcode{=.true. }
+file \forcode{numror} can be handled by XIOS. To activate restart reading using XIOS, set \np[=.true. ]{ln_xios_read}{ln\_xios\_read}
 in \textit{namelist\_cfg}. This setting will be ignored when multiple restart files are present, and default \NEMO
 functionality will be used for reading. There is no need to change iodef.xml file to use XIOS to read
@@ -143,6 +145,6 @@
 type of restart \NEMO\ will write. If it is set to 0, default \NEMO\ functionality will be used - each
 processor writes its own restart file; if it is set to 1 XIOS will write restart into a single file;
-for \np{nn_wxios}{nn\_wxios}\forcode{=2} the restart will be written by XIOS into multiple files, one for each XIOS server.
-Note, however, that \textbf{\NEMO\ will not read restart generated by XIOS when \np{nn_wxios}{nn\_wxios}\forcode{=2}}. The restart will
+for \np[=2]{nn_wxios}{nn\_wxios} the restart will be written by XIOS into multiple files, one for each XIOS server.
+Note, however, that \textbf{\NEMO\ will not read restart generated by XIOS when \np[=2]{nn_wxios}{nn\_wxios}}. The restart will
 have to be rebuild before continuing the run. This option aims to reduce number of restart files generated by \NEMO\ only,
 and may be useful when there is a need to change number of processors used to run simulation.
@@ -1470,21 +1472,21 @@
 
 \begin{description}
-\item[\np{ln_glo_trd}{ln\_glo\_trd}]:
+\item[{\np{ln_glo_trd}{ln\_glo\_trd}}]:
   at each \np{nn_trd}{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}{ln\_dyn\_trd}]:
+\item[{\np{ln_dyn_trd}{ln\_dyn\_trd}}]:
   each 3D trend of the evolution of the two momentum components is output;
-\item[\np{ln_dyn_mxl}{ln\_dyn\_mxl}]:
+\item[{\np{ln_dyn_mxl}{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}{ln\_vor\_trd}]:
+\item[{\np{ln_vor_trd}{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}{ln\_KE\_trd}] :
+\item[{\np{ln_KE_trd}{ln\_KE\_trd}}] :
   each 3D trend of the Kinetic Energy equation is output;
-\item[\np{ln_tra_trd}{ln\_tra\_trd}]:
+\item[{\np{ln_tra_trd}{ln\_tra\_trd}}]:
   each 3D trend of the evolution of temperature and salinity is output;
-\item[\np{ln_tra_mxl}{ln\_tra\_mxl}]:
+\item[{\np{ln_tra_mxl}{ln\_tra\_mxl}}]:
   each 2D trend of the evolution of temperature and salinity averaged over the mixed layer is output;
 \end{description}
@@ -1495,5 +1497,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}{ln\_dyn\_mxl}, \np{ln_vor_trd}{ln\_vor\_trd}, and \np{ln_tra_mxl}{ln\_tra\_mxl} are not working,
-and none of the options have been tested with variable volume (\ie\ \np{ln_linssh}{ln\_linssh}\forcode{=.true.}).
+and none of the options have been tested with variable volume (\ie\ \np[=.true.]{ln_linssh}{ln\_linssh}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1515,5 +1517,5 @@
 Options are defined by \nam{flo}{flo} namelist 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}{ln\_flork4}\forcode{=.true.}).
+or on a $4^th$ Runge-Hutta algorithm (\np[=.true.]{ln_flork4}{ln\_flork4}).
 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.
@@ -1522,5 +1524,5 @@
 
 Initial coordinates can be given with Ariane Tools convention
-(IJK coordinates, (\np{ln_ariane}{ln\_ariane}\forcode{=.true.}) ) or with longitude and latitude.
+(IJK coordinates, (\np[=.true.]{ln_ariane}{ln\_ariane}) ) or with longitude and latitude.
 
 In case of Ariane convention, input filename is \textit{init\_float\_ariane}.
@@ -1573,5 +1575,5 @@
 
 \np{jpnfl}{jpnfl} is the total number of floats during the run.
-When initial positions are read in a restart file (\np{ln_rstflo}{ln\_rstflo}\forcode{=.true.} ),
+When initial positions are read in a restart file (\np[=.true.]{ln_rstflo}{ln\_rstflo} ),
 \np{jpnflnewflo}{jpnflnewflo} can be added in the initialization file.
 
@@ -1581,8 +1583,8 @@
 creation of the float restart file.
 
-Output data can be written in ascii files (\np{ln_flo_ascii}{ln\_flo\_ascii}\forcode{=.true.}).
+Output data can be written in ascii files (\np[=.true.]{ln_flo_ascii}{ln\_flo\_ascii}).
 In that case, output filename is trajec\_float.
 
-Another possiblity of writing format is Netcdf (\np{ln_flo_ascii}{ln\_flo\_ascii}\forcode{=.false.}) with
+Another possiblity of writing format is Netcdf (\np[=.false.]{ln_flo_ascii}{ln\_flo\_ascii}) with
 \key{iomput} and outputs selected in iodef.xml.
 Here it is an example of specification to put in files description section:
@@ -1939,5 +1941,5 @@
 
 Third, the discretisation of \autoref{eq:DIA_steric_Bq} depends on the type of free surface which is considered.
-In the non linear free surface case, \ie\ \np{ln_linssh}{ln\_linssh}\forcode{=.true.}, it is given by
+In the non linear free surface case, \ie\ \np[=.true.]{ln_linssh}{ln\_linssh}, it is given by
 
 \[
@@ -2034,8 +2036,8 @@
 sea water pressure at sea floor (botpres), dynamic sea surface height (sshdyn).
 
-In \mdl{diaptr} when \np{ln_diaptr}{ln\_diaptr}\forcode{=.true.}
+In \mdl{diaptr} when \np[=.true.]{ln_diaptr}{ln\_diaptr}
 (see the \nam{ptr}{ptr} namelist below) can be computed on-line the poleward heat and salt transports,
 their advective and diffusive component, and the meriodional stream function .
-When \np{ln_subbas}{ln\_subbas}\forcode{=.true.}, transports and stream function are computed for the Atlantic, Indian,
+When \np[=.true.]{ln_subbas}{ln\_subbas}, transports and stream function are computed for the Atlantic, Indian,
 Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean.
 The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays,
@@ -2119,7 +2121,7 @@
 % ================================================================
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -46,7 +48,7 @@
 This namelist contains only two variables:
 \begin{description}
-\item[\np{ln_diurnal}{ln\_diurnal}]
+\item[{\np{ln_diurnal}{ln\_diurnal}}]
   A logical switch for turning on/off both the cool skin and warm layer.
-\item[\np{ln_diurnal_only}{ln\_diurnal\_only}]
+\item[{\np{ln_diurnal_only}{ln\_diurnal\_only}}]
   A logical switch which if \forcode{.true.} will run the diurnal model without the other dynamical parts of \NEMO.
   \np{ln_diurnal_only}{ln\_diurnal\_only} must be \forcode{.false.} if \np{ln_diurnal}{ln\_diurnal} is \forcode{.false.}.
@@ -159,7 +161,7 @@
 \]
 
-\biblio
+\onlyinsubfile{\bibliography{../main/bibliography}}
 
-\pindex
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -487,5 +489,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}{ln\_linssh}\forcode{=.false.}).
+    (f) same as (e) but in the non-linear free surface (\protect\np[=.false.]{ln_linssh}{ln\_linssh}).
     Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).}
   \label{fig:DOM_z_zps_s_sps}
@@ -524,7 +526,7 @@
 
 \begin{itemize}
-\item $z$-coordinate with full step bathymetry (\np{ln_zco}{ln\_zco}\forcode{=.true.}),
-\item $z$-coordinate with partial step ($zps$) bathymetry (\np{ln_zps}{ln\_zps}\forcode{=.true.}),
-\item Generalized, $s$-coordinate (\np{ln_sco}{ln\_sco}\forcode{=.true.}).
+\item $z$-coordinate with full step bathymetry (\np[=.true.]{ln_zco}{ln\_zco}),
+\item $z$-coordinate with partial step ($zps$) bathymetry (\np[=.true.]{ln_zps}{ln\_zps}),
+\item Generalized, $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}).
 \end{itemize}
 
@@ -544,5 +546,5 @@
 They are updated at each model time step.
 The initial fixed reference coordinate system is held in variable names with a $\_0$ suffix.
-When the linear free surface option is used (\np{ln_linssh}{ln\_linssh}\forcode{=.true.}),
+When the linear free surface option is used (\np[=.true.]{ln_linssh}{ln\_linssh}),
 \textit{before}, \textit{now} and \textit{after} arrays are initially set to
 their reference counterpart and remain fixed.
@@ -682,5 +684,5 @@
 
 \begin{description}
-\item[\np{ln_tsd_init}{ln\_tsd\_init}\forcode{= .true.}]
+\item[{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}]
   Use T and S input files that can be given on the model grid itself or on their native input data grids.
   In the latter case, the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid
@@ -688,5 +690,5 @@
   The information relating to the input files are specified in the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures.
   The computation is done in the \mdl{dtatsd} module.
-\item[\np{ln_tsd_init}{ln\_tsd\_init}\forcode{= .false.}]
+\item[{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}]
   Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}.
   The default version sets horizontally uniform T and profiles as used in the GYRE configuration
@@ -694,7 +696,7 @@
 \end{description}
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -205,5 +207,5 @@
 (EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}).
 In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of
-vorticity term with analytical equations (\np{ln_dynvor_con}{ln\_dynvor\_con}\forcode{=.true.}).
+vorticity term with analytical equations (\np[=.true.]{ln_dynvor_con}{ln\_dynvor\_con}).
 The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module.
 
@@ -327,5 +329,5 @@
 A key point in \autoref{eq:DYN_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}{nn\_een\_e3f}\forcode{=1}), or just by $4$ (\np{nn_een_e3f}{nn\_een\_e3f}\forcode{=.true.}).
+(\np[=1]{nn_een_e3f}{nn\_een\_e3f}), or just by $4$ (\np[=.true.]{nn_een_e3f}{nn\_een\_e3f}).
 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.
@@ -407,5 +409,5 @@
   \right.
 \]
-When \np{ln_dynzad_zts}{ln\_dynzad\_zts}\forcode{=.true.},
+When \np[=.true.]{ln_dynzad_zts}{ln\_dynzad\_zts},
 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.debreu.ea_OM15}.
@@ -534,5 +536,5 @@
 But the amplitudes of the false extrema are significantly reduced over 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
-(\ie\ \np{ln_dynldf_lap}{ln\_dynldf\_lap}\forcode{=}\np{ln_dynldf_bilap}{ln\_dynldf\_bilap}\forcode{=.false.}),
+(\ie\ \np[=]{ln_dynldf_lap}{ln\_dynldf\_lap}\np[=.false.]{ln_dynldf_bilap}{ln\_dynldf\_bilap}),
 and it is recommended to do so.
 
@@ -665,5 +667,5 @@
 density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation.
 
-$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln_dynhpg_sco}{ln\_dynhpg\_sco}\forcode{=.true.})
+$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco})
 \begin{equation}
   \label{eq:DYN_hpg_sco}
@@ -683,17 +685,17 @@
 ($e_{3w}$).
 
-$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln_dynhpg_isf}{ln\_dynhpg\_isf}\forcode{=.true.}).
-This scheme need the activation of ice shelf cavities (\np{ln_isfcav}{ln\_isfcav}\forcode{=.true.}).
-
-$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln_dynhpg_prj}{ln\_dynhpg\_prj}\forcode{=.true.})
+$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}).
+This scheme need the activation of ice shelf cavities (\np[=.true.]{ln_isfcav}{ln\_isfcav}).
+
+$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj})
 
 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05}
-(\np{ln_dynhpg_djc}{ln\_dynhpg\_djc}\forcode{=.true.}) (currently disabled; under development)
+(\np[=.true.]{ln_dynhpg_djc}{ln\_dynhpg\_djc}) (currently disabled; under development)
 
 Note that expression \autoref{eq:DYN_hpg_sco} is commonly used when the variable volume formulation is activated
 (\texttt{vvl?}) because in that case, even with a flat bottom,
 the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}.
-The pressure jacobian scheme (\np{ln_dynhpg_prj}{ln\_dynhpg\_prj}\forcode{=.true.}) is available as
-an improved option to \np{ln_dynhpg_sco}{ln\_dynhpg\_sco}\forcode{=.true.} when \texttt{vvl?} is active.
+The pressure jacobian scheme (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj}) is available as
+an improved option to \np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco} when \texttt{vvl?} is active.
 The pressure Jacobian scheme uses a constrained cubic spline to
 reconstruct the density profile across the water column.
@@ -707,5 +709,5 @@
 
 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 (\np{ln_dynhpg_isf}{ln\_dynhpg\_isf}\forcode{=.true.}).\\
+the pressure gradient due to the ocean load (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}).\\
 
 The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium.
@@ -738,5 +740,5 @@
 rather than at the central time level $t$ only, as in the standard leapfrog scheme.
 
-$\bullet$ leapfrog scheme (\np{ln_dynhpg_imp}{ln\_dynhpg\_imp}\forcode{=.true.}):
+$\bullet$ leapfrog scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
 
 \begin{equation}
@@ -746,5 +748,5 @@
 \end{equation}
 
-$\bullet$ semi-implicit scheme (\np{ln_dynhpg_imp}{ln\_dynhpg\_imp}\forcode{=.true.}):
+$\bullet$ semi-implicit scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
 \begin{equation}
   \label{eq:DYN_hpg_imp}
@@ -764,5 +766,5 @@
 such as the stability limits associated with advection or diffusion.
 
-In practice, the semi-implicit scheme is used when \np{ln_dynhpg_imp}{ln\_dynhpg\_imp}\forcode{=.true.}.
+In practice, the semi-implicit scheme is used when \np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}.
 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 density,
@@ -865,5 +867,5 @@
 The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through
 the \np{nn_baro}{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$.
-This parameter can be optionally defined automatically (\np{ln_bt_nn_auto}{ln\_bt\_nn\_auto}\forcode{=.true.}) considering that
+This parameter can be optionally defined automatically (\np[=.true.]{ln_bt_nn_auto}{ln\_bt\_nn\_auto}) 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.
@@ -904,5 +906,5 @@
     In this particular exemple,
     a boxcar averaging window over \np{nn_baro}{nn\_baro} barotropic time steps is used
-    (\np{nn_bt_flt}{nn\_bt\_flt}\forcode{=1}) and \np{nn_baro}{nn\_baro}\forcode{=5}.
+    (\np[=1]{nn_bt_flt}{nn\_bt\_flt}) and \np[=5]{nn_baro}{nn\_baro}.
     Internal mode time steps (which are also the model time steps) are denoted by
     $t-\rdt$, $t$ and $t+\rdt$.
@@ -913,18 +915,18 @@
     the latter are used to obtain time averaged transports to advect tracers.
     a) Forward time integration:
-    \protect\np{ln_bt_fw}{ln\_bt\_fw}\forcode{=.true.},  \protect\np{ln_bt_av}{ln\_bt\_av}\forcode{=.true.}.
+    \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw},  \protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}.
     b) Centred time integration:
-    \protect\np{ln_bt_fw}{ln\_bt\_fw}\forcode{=.false.}, \protect\np{ln_bt_av}{ln\_bt\_av}\forcode{=.true.}.
+    \protect\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}, \protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}.
     c) Forward time integration with no time filtering (POM-like scheme):
-    \protect\np{ln_bt_fw}{ln\_bt\_fw}\forcode{=.true.},  \protect\np{ln_bt_av}{ln\_bt\_av}\forcode{=.false.}.}
+    \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw},  \protect\np[=.false.]{ln_bt_av}{ln\_bt\_av}.}
   \label{fig:DYN_spg_ts}
 \end{figure}
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
 
-In the default case (\np{ln_bt_fw}{ln\_bt\_fw}\forcode{=.true.}),
+In the default case (\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}),
 the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps
 (\autoref{fig:DYN_spg_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}{ln\_bt\_av}\forcode{=.true.}).
+time filtering is eventually applied on barotropic quantities (\np[=.true.]{ln_bt_av}{ln\_bt\_av}).
 In that case, the integration is extended slightly beyond \textit{after} time step to
 provide time filtered quantities.
@@ -933,5 +935,5 @@
 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}{ln\_bt\_fw}\forcode{=.false.}).
+(\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}).
 Although more computationaly expensive ( \np{nn_baro}{nn\_baro} additional iterations are indeed necessary),
 the baroclinic to barotropic forcing term given at \textit{now} time step become centred in
@@ -958,5 +960,5 @@
 
 One can eventually choose to feedback instantaneous values by not using any time filter
-(\np{ln_bt_av}{ln\_bt\_av}\forcode{=.false.}).
+(\np[=.false.]{ln_bt_av}{ln\_bt\_av}).
 In that case, external mode equations are continuous in time,
 \ie\ they are not re-initialized when starting a new sub-stepping sequence.
@@ -1191,9 +1193,9 @@
 
 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}{ln\_dynldf\_iso}\forcode{=.true.}) and
-for either iso-neutral (\np{ln_dynldf_iso}{ln\_dynldf\_iso}\forcode{=.true.}) or
-geopotential (\np{ln_dynldf_hor}{ln\_dynldf\_hor}\forcode{=.true.}) diffusion in the $s$-coordinate.
+for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and
+for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or
+geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) 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}{ln\_dynldf\_hor}\forcode{=.true.}.
+no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.
 The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on
 each momentum component.
@@ -1269,6 +1271,6 @@
 Two time stepping schemes can be used for the vertical diffusion term:
 $(a)$ a forward time differencing scheme
-(\np{ln_zdfexp}{ln\_zdfexp}\forcode{=.true.}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or
-$(b)$ a backward (or implicit) time differencing scheme (\np{ln_zdfexp}{ln\_zdfexp}\forcode{=.false.})
+(\np[=.true.]{ln_zdfexp}{ln\_zdfexp}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or
+$(b)$ a backward (or implicit) time differencing scheme (\np[=.false.]{ln_zdfexp}{ln\_zdfexp})
 (see \autoref{chap:TD}).
 Note that namelist variables \np{ln_zdfexp}{ln\_zdfexp} and \np{nn_zdfexp}{nn\_zdfexp} apply to both tracers and dynamics.
@@ -1320,11 +1322,11 @@
 three other forcings may enter the dynamical equations by affecting the surface pressure gradient.
 
-(1) When \np{ln_apr_dyn}{ln\_apr\_dyn}\forcode{=.true.} (see \autoref{sec:SBC_apr}),
+(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
 the atmospheric pressure is taken into account when computing the surface pressure gradient.
 
-(2) When \np{ln_tide_pot}{ln\_tide\_pot}\forcode{=.true.} and \np{ln_tide}{ln\_tide}\forcode{=.true.} (see \autoref{sec:SBC_tide}),
+(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_tide}),
 the tidal potential is taken into account when computing the surface pressure gradient.
 
-(3) When \np{nn_ice_embd}{nn\_ice\_embd}\forcode{=2} and LIM or CICE is used
+(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and LIM or CICE is used
 (\ie\ when the sea-ice is embedded in the ocean),
 the snow-ice mass is taken into account when computing the surface pressure gradient.
@@ -1406,7 +1408,7 @@
 
 The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).
-If the user sets \np{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp}\forcode{=.false.} then zuwdmask is 1 when the
+If the user sets \np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then zuwdmask is 1 when the
 flux is from a cell with water depth greater than \np{rn_wdmin1}{rn\_wdmin1} and 0 otherwise. If the user sets
-\np{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp}\forcode{=.true.} the flux across the face is ramped down as the water depth decreases
+\np[=.true.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} the flux across the face is ramped down as the water depth decreases
 from 2 * \np{rn_wdmin1}{rn\_wdmin1} to \np{rn_wdmin1}{rn\_wdmin1}. The use of this ramp reduced grid-scale noise in idealised test cases.
 
@@ -1425,5 +1427,5 @@
 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because
 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts
-to equal their mean value during the barotropic steps. If the user sets \np{ln_wd_dl_bc}{ln\_wd\_dl\_bc}\forcode{=.true.}, the
+to equal their mean value during the barotropic steps. If the user sets \np[=.true.]{ln_wd_dl_bc}{ln\_wd\_dl\_bc}, the
 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.
 
@@ -1658,5 +1660,5 @@
 
 $\bullet$ vector invariant form or linear free surface
-(\np{ln_dynhpg_vec}{ln\_dynhpg\_vec}\forcode{=.true.} ; \texttt{vvl?} not defined):
+(\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} not defined):
 \[
   % \label{eq:DYN_nxt_vec}
@@ -1670,5 +1672,5 @@
 
 $\bullet$ flux form and nonlinear free surface
-(\np{ln_dynhpg_vec}{ln\_dynhpg\_vec}\forcode{=.false.} ; \texttt{vvl?} defined):
+(\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} defined):
 \[
   % \label{eq:DYN_nxt_flux}
@@ -1684,5 +1686,5 @@
 the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient.
 $\gamma$ is initialized as \np{nn_atfp}{nn\_atfp} (namelist parameter).
-Its default value is \np{nn_atfp}{nn\_atfp}\forcode{=10.e-3}.
+Its default value is \np[=10.e-3]{nn_atfp}{nn\_atfp}.
 In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for
 the momentum equations.
@@ -1693,7 +1695,7 @@
 
 % ================================================================
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -101,7 +103,7 @@
   \caption[Lateral boundary conditions]{
     Lateral boundary conditions
-    (a) free-slip                       (\protect\np{rn_shlat}{rn\_shlat}\forcode{=0});
-    (b) no-slip                         (\protect\np{rn_shlat}{rn\_shlat}\forcode{=2});
-    (c) "partial" free-slip (\forcode{0<}\protect\np{rn_shlat}{rn\_shlat}\forcode{<2}) and
+    (a) free-slip                       (\protect\np[=0]{rn_shlat}{rn\_shlat});
+    (b) no-slip                         (\protect\np[=2]{rn_shlat}{rn\_shlat});
+    (c) "partial" free-slip (\forcode{0<}\protect\np[<2]{rn_shlat}{rn\_shlat}) and
     (d) "strong" no-slip    (\forcode{2<}\protect\np{rn_shlat}{rn\_shlat}).
     Implied "ghost" velocity inside land area is display in grey.}
@@ -112,5 +114,5 @@
 \begin{description}
 
-\item[free-slip boundary condition (\np{rn_shlat}{rn\_shlat}\forcode{=0}):] the tangential velocity at
+\item[free-slip boundary condition ({\np[=0]{rn_shlat}{rn\_shlat}})] the tangential velocity at
   the coastline is equal to the offshore velocity,
   \ie\ the normal derivative of the tangential velocity is zero at the coast,
@@ -118,5 +120,5 @@
   (\autoref{fig:LBC_shlat}-a).
 
-\item[no-slip boundary condition (\np{rn_shlat}{rn\_shlat}\forcode{=2}):] the tangential velocity vanishes at the coastline.
+\item[no-slip boundary condition ({\np[=2]{rn_shlat}{rn\_shlat}})] the tangential velocity vanishes at the coastline.
   Assuming that the tangential velocity decreases linearly from
   the closest ocean velocity grid point to the coastline,
@@ -139,10 +141,10 @@
   \]
 
-\item["partial" free-slip boundary condition (0$<$\np{rn_shlat}{rn\_shlat}$<$2):] the tangential velocity at
+\item["partial" free-slip boundary condition (0$<$\np{rn_shlat}{rn\_shlat}$<$2)] the tangential velocity at
   the coastline is smaller than the offshore velocity, \ie\ there is a lateral friction but
   not strong enough to make the tangential velocity at the coast vanish (\autoref{fig:LBC_shlat}-c).
   This can be selected by providing a value of mask$_{f}$ strictly inbetween $0$ and $2$.
 
-\item["strong" no-slip boundary condition (2$<$\np{rn_shlat}{rn\_shlat}):] the viscous boundary layer is assumed to
+\item["strong" no-slip boundary condition (2$<$\np{rn_shlat}{rn\_shlat})] the viscous boundary layer is assumed to
   be smaller than half the grid size (\autoref{fig:LBC_shlat}-d).
   The friction is thus larger than in the no-slip case.
@@ -393,9 +395,9 @@
 \label{subsec:LBC_bdy_namelist}
 
-The BDY module is activated by setting \np{ln_bdy}{ln\_bdy}\forcode{=.true.} .
+The BDY module is activated by setting \np[=.true.]{ln_bdy}{ln\_bdy} .
 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}{nb\_bdy}.
 Each boundary set can be either defined as a series of straight line segments directly in the namelist
-(\np{ln_coords_file}{ln\_coords\_file}\forcode{=.false.}, and a namelist block \nam{bdy_index}{bdy\_index} must be included for each set) or read in from a file (\np{ln_coords_file}{ln\_coords\_file}\forcode{=.true.}, and a ``\ifile{coordinates.bdy}'' file must be provided).
+(\np[=.false.]{ln_coords_file}{ln\_coords\_file}, and a namelist block \nam{bdy_index}{bdy\_index} must be included for each set) or read in from a file (\np[=.true.]{ln_coords_file}{ln\_coords\_file}, and a ``\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
@@ -422,9 +424,9 @@
 
 The boundary data is either set to initial conditions
-(\np{nn_tra_dta}{nn\_tra\_dta}\forcode{=0}) or forced with external data from a file (\np{nn_tra_dta}{nn\_tra\_dta}\forcode{=1}).
+(\np[=0]{nn_tra_dta}{nn\_tra\_dta}) or forced with external data from a file (\np[=1]{nn_tra_dta}{nn\_tra\_dta}).
 In case the 3d velocity data contain the total velocity (ie, baroclinic and barotropic velocity),
-the bdy code can derived baroclinic and barotropic velocities by setting \np{ln_full_vel}{ln\_full\_vel}\forcode{=.true.}
+the bdy code can derived baroclinic and barotropic velocities by setting \np[=.true.]{ln_full_vel}{ln\_full\_vel}
 For the barotropic solution there is also the option to use tidal harmonic forcing either by
-itself (\np{nn_dyn2d_dta}{nn\_dyn2d\_dta}\forcode{=2}) or in addition to other external data (\np{nn_dyn2d_dta}{nn\_dyn2d\_dta}\forcode{=3}).\\
+itself (\np[=2]{nn_dyn2d_dta}{nn\_dyn2d\_dta}) or in addition to other external data (\np[=3]{nn_dyn2d_dta}{nn\_dyn2d\_dta}).\\
 If not set to initial conditions, sea-ice salinity, temperatures and melt ponds data at the boundary can either be read in a file or defined as constant (by \np{rn_ice_sal}{rn\_ice\_sal}, \np{rn_ice_tem}{rn\_ice\_tem}, \np{rn_ice_apnd}{rn\_ice\_apnd}, \np{rn_ice_hpnd}{rn\_ice\_hpnd}). Ice age is constant and defined by \np{rn_ice_age}{rn\_ice\_age}.
 
@@ -602,5 +604,5 @@
 \jp{jpinft} give the start and end $i$ indices for each segment with similar for the other boundaries.
 These segments define a list of $T$ grid points along the outermost row of the boundary ($nbr\,=\, 1$).
-The code deduces the $U$ and $V$ points and also the points for $nbr\,>\, 1$ if \np{nn_rimwidth}{nn\_rimwidth}\forcode{>1}.
+The code deduces the $U$ and $V$ points and also the points for $nbr\,>\, 1$ if \np[>1]{nn_rimwidth}{nn\_rimwidth}.
 
 The boundary geometry may also be defined from a ``\ifile{coordinates.bdy}'' file.
@@ -673,9 +675,9 @@
 There is an option to force the total volume in the regional model to be constant.
 This is controlled  by the \np{ln_vol}{ln\_vol} parameter in the namelist.
-A value of \np{ln_vol}{ln\_vol}\forcode{=.false.} indicates that this option is not used.
+A value of \np[=.false.]{ln_vol}{ln\_vol} indicates that this option is not used.
 Two options to control the volume are available (\np{nn_volctl}{nn\_volctl}).
-If \np{nn_volctl}{nn\_volctl}\forcode{=0} then a correction is applied to the normal barotropic velocities around the boundary at
+If \np[=0]{nn_volctl}{nn\_volctl} then a correction is applied to the normal barotropic velocities around the boundary at
 each timestep to ensure that the integrated volume flow through the boundary is zero.
-If \np{nn_volctl}{nn\_volctl}\forcode{=1} then the calculation of the volume change on
+If \np[=1]{nn_volctl}{nn\_volctl} 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.
@@ -741,7 +743,7 @@
 direction of rotation). %, e.g. anticlockwise or clockwise.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -24,5 +26,5 @@
 (see the \nam{tra_ldf}{tra\_ldf} and \nam{dyn_ldf}{dyn\_ldf} below).
 Note that this chapter describes the standard implementation of iso-neutral tracer mixing.
-Griffies's implementation, which is used if \np{ln_traldf_triad}{ln\_traldf\_triad}\forcode{=.true.},
+Griffies's implementation, which is used if \np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad},
 is described in \autoref{apdx:TRIADS}
 
@@ -40,16 +42,16 @@
 \subsection[No lateral mixing (\forcode{ln_traldf_OFF} \& \forcode{ln_dynldf_OFF})]{No lateral mixing (\protect\np{ln_traldf_OFF}{ln\_traldf\_OFF} \& \protect\np{ln_dynldf_OFF}{ln\_dynldf\_OFF})}
 
-It is possible to run without explicit lateral diffusion on tracers (\protect\np{ln_traldf_OFF}{ln\_traldf\_OFF}\forcode{=.true.}) and/or
-momentum (\protect\np{ln_dynldf_OFF}{ln\_dynldf\_OFF}\forcode{=.true.}). The latter option is even recommended if using the
-UBS advection scheme on momentum (\np{ln_dynadv_ubs}{ln\_dynadv\_ubs}\forcode{=.true.},
+It is possible to run without explicit lateral diffusion on tracers (\protect\np[=.true.]{ln_traldf_OFF}{ln\_traldf\_OFF}) and/or
+momentum (\protect\np[=.true.]{ln_dynldf_OFF}{ln\_dynldf\_OFF}). The latter option is even recommended if using the
+UBS advection scheme on momentum (\np[=.true.]{ln_dynadv_ubs}{ln\_dynadv\_ubs},
 see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.
 
 \subsection[Laplacian mixing (\forcode{ln_traldf_lap} \& \forcode{ln_dynldf_lap})]{Laplacian mixing (\protect\np{ln_traldf_lap}{ln\_traldf\_lap} \& \protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})}
-Setting \protect\np{ln_traldf_lap}{ln\_traldf\_lap}\forcode{=.true.} and/or \protect\np{ln_dynldf_lap}{ln\_dynldf\_lap}\forcode{=.true.} enables
+Setting \protect\np[=.true.]{ln_traldf_lap}{ln\_traldf\_lap} and/or \protect\np[=.true.]{ln_dynldf_lap}{ln\_dynldf\_lap} enables
 a second order diffusion on tracers and momentum respectively. Note that in \NEMO\ 4, one can not combine
 Laplacian and Bilaplacian operators for the same variable.
 
 \subsection[Bilaplacian mixing (\forcode{ln_traldf_blp} \& \forcode{ln_dynldf_blp})]{Bilaplacian mixing (\protect\np{ln_traldf_blp}{ln\_traldf\_blp} \& \protect\np{ln_dynldf_blp}{ln\_dynldf\_blp})}
-Setting \protect\np{ln_traldf_blp}{ln\_traldf\_blp}\forcode{=.true.} and/or \protect\np{ln_dynldf_blp}{ln\_dynldf\_blp}\forcode{=.true.} enables
+Setting \protect\np[=.true.]{ln_traldf_blp}{ln\_traldf\_blp} and/or \protect\np[=.true.]{ln_dynldf_blp}{ln\_dynldf\_blp} enables
 a fourth order diffusion on tracers and momentum respectively. It is implemented by calling the above Laplacian operator twice.
 We stress again that from \NEMO\ 4, the simultaneous use Laplacian and Bilaplacian operators is not allowed.
@@ -107,6 +109,6 @@
 %gm%  caution I'm not sure the simplification was a good idea!
 
-These slopes are computed once in \rou{ldf\_slp\_init} when \np{ln_sco}{ln\_sco}\forcode{=.true.},
-and either \np{ln_traldf_hor}{ln\_traldf\_hor}\forcode{=.true.} or \np{ln_dynldf_hor}{ln\_dynldf\_hor}\forcode{=.true.}.
+These slopes are computed once in \rou{ldf\_slp\_init} when \np[=.true.]{ln_sco}{ln\_sco},
+and either \np[=.true.]{ln_traldf_hor}{ln\_traldf\_hor} or \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.
 
 \subsection{Slopes for tracer iso-neutral mixing}
@@ -164,8 +166,8 @@
 \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{ln_traldf_triad}{ln\_traldf\_triad}\forcode{=.true.};
+  the Griffies scheme is used (\np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad};
   see \autoref{apdx:TRIADS}).
   In other words, iso-neutral mixing will only be accurately represented with a linear equation of state
-  (\np{ln_seos}{ln\_seos}\forcode{=.true.}).
+  (\np[=.true.]{ln_seos}{ln\_seos}).
   In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:LDF_slp_iso}
   will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes.
@@ -222,5 +224,5 @@
 To overcome this problem, several techniques have been proposed in which the numerical schemes of
 the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}.
-Griffies's scheme is now available in \NEMO\ if \np{ln_traldf_triad}{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:TRIADS}.
+Griffies's scheme is now available in \NEMO\ if \np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad}; see \autoref{apdx:TRIADS}.
 Here, another strategy is presented \citep{lazar_phd97}:
 a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of
@@ -326,5 +328,5 @@
 The way the mixing coefficients are set in the reference version can be described as follows:
 
-\subsection[Mixing coefficients read from file (\forcode{=-20, -30})]{ Mixing coefficients read from file (\protect\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=-20, -30} \& \protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=-20, -30})}
+\subsection[Mixing coefficients read from file (\forcode{=-20, -30})]{ Mixing coefficients read from file (\protect\np[=-20, -30]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=-20, -30]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
 Mixing coefficients can be read from file if a particular geographical variation is needed. For example, in the ORCA2 global ocean model,
@@ -332,5 +334,5 @@
 decreases linearly to $A^l$~= 2.10$^3$ m$^2$/s at the equator \citep{madec.delecluse.ea_JPO96, delecluse.madec_icol99}.
 Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of ORCA2 and ORCA05.
-The provided fields can either be 2d (\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=-20}, \np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=-20}) or 3d (\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=-30},  \np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=-30}). They must be given at U, V points for tracers and T, F points for momentum (see \autoref{tab:LDF_files}).
+The provided fields can either be 2d (\np[=-20]{nn_aht_ijk_t}{nn\_aht\_ijk\_t}, \np[=-20]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}) or 3d (\np[=-30]{nn_aht_ijk_t}{nn\_aht\_ijk\_t},  \np[=-30]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}). They must be given at U, V points for tracers and T, F points for momentum (see \autoref{tab:LDF_files}).
 
 %-------------------------------------------------TABLE---------------------------------------------------
@@ -340,8 +342,8 @@
     \hline
     Namelist parameter				            & Input filename	                             & dimensions	& variable names                      \\	\hline
-    \np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=-20}	    & \forcode{eddy_viscosity_2D.nc }            &     $(i,j)$	        & \forcode{ahmt_2d, ahmf_2d}  \\	\hline
-    \np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=-20}           & \forcode{eddy_diffusivity_2D.nc }           &     $(i,j)$	        & \forcode{ahtu_2d, ahtv_2d}    \\	\hline
-    \np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=-30} 	    & \forcode{eddy_viscosity_3D.nc }            &     $(i,j,k)$	        & \forcode{ahmt_3d, ahmf_3d}  \\	\hline
-    \np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=-30}	    & \forcode{eddy_diffusivity_3D.nc }           &     $(i,j,k)$	        & \forcode{ahtu_3d, ahtv_3d}    \\	\hline
+    \np[=-20]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}	    & \forcode{eddy_viscosity_2D.nc }            &     $(i,j)$	        & \forcode{ahmt_2d, ahmf_2d}  \\	\hline
+    \np[=-20]{nn_aht_ijk_t}{nn\_aht\_ijk\_t}           & \forcode{eddy_diffusivity_2D.nc }           &     $(i,j)$	        & \forcode{ahtu_2d, ahtv_2d}    \\	\hline
+    \np[=-30]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t} 	    & \forcode{eddy_viscosity_3D.nc }            &     $(i,j,k)$	        & \forcode{ahmt_3d, ahmf_3d}  \\	\hline
+    \np[=-30]{nn_aht_ijk_t}{nn\_aht\_ijk\_t}	    & \forcode{eddy_diffusivity_3D.nc }           &     $(i,j,k)$	        & \forcode{ahtu_3d, ahtv_3d}    \\	\hline
   \end{tabular}
   \caption{Description of expected input files if mixing coefficients are read from NetCDF files}
@@ -350,5 +352,5 @@
 %--------------------------------------------------------------------------------------------------------------
 
-\subsection[Constant mixing coefficients (\forcode{=0})]{ Constant mixing coefficients (\protect\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=0} \& \protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=0})}
+\subsection[Constant mixing coefficients (\forcode{=0})]{ Constant mixing coefficients (\protect\np[=0]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=0]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
 If constant, mixing coefficients are set thanks to a velocity and a length scales ($U_{scl}$, $L_{scl}$) such that:
@@ -366,5 +368,5 @@
  $U_{scl}$ and $L_{scl}$ are given by the namelist parameters \np{rn_Ud}{rn\_Ud}, \np{rn_Uv}{rn\_Uv}, \np{rn_Ld}{rn\_Ld} and \np{rn_Lv}{rn\_Lv}.
 
-\subsection[Vertically varying mixing coefficients (\forcode{=10})]{Vertically varying mixing coefficients (\protect\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=10} \& \protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=10})}
+\subsection[Vertically varying mixing coefficients (\forcode{=10})]{Vertically varying mixing coefficients (\protect\np[=10]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=10]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
 In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which
@@ -373,5 +375,5 @@
 This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users.
 
-\subsection[Mesh size dependent mixing coefficients (\forcode{=20})]{Mesh size dependent mixing coefficients (\protect\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=20} \& \protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=20})}
+\subsection[Mesh size dependent mixing coefficients (\forcode{=20})]{Mesh size dependent mixing coefficients (\protect\np[=20]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=20]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
 In that case, the horizontal variation of the eddy coefficient depends on the local mesh size and
@@ -398,5 +400,5 @@
 \colorbox{yellow}{CASE \np{nn_aht_ijk_t}{nn\_aht\_ijk\_t} = 21 to be added}
 
-\subsection[Mesh size and depth dependent mixing coefficients (\forcode{=30})]{Mesh size and depth dependent mixing coefficients (\protect\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=30} \& \protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=30})}
+\subsection[Mesh size and depth dependent mixing coefficients (\forcode{=30})]{Mesh size and depth dependent mixing coefficients (\protect\np[=30]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=30]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
 The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases above,
@@ -404,5 +406,5 @@
 the magnitude of the coefficient.
 
-\subsection[Velocity dependent mixing coefficients (\forcode{=31})]{Flow dependent mixing coefficients (\protect\np{nn_aht_ijk_t}{nn\_aht\_ijk\_t}\forcode{=31} \& \protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=31})}
+\subsection[Velocity dependent mixing coefficients (\forcode{=31})]{Flow dependent mixing coefficients (\protect\np[=31]{nn_aht_ijk_t}{nn\_aht\_ijk\_t} \& \protect\np[=31]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 In that case, the eddy coefficient is proportional to the local velocity magnitude so that the Reynolds number $Re =  \lvert U \rvert  e / A_l$ is constant (and here hardcoded to $12$):
 \colorbox{yellow}{JC comment: The Reynolds is effectively set to 12 in the code for both operators but shouldn't it be 2 for Laplacian ?}
@@ -418,5 +420,5 @@
 \end{equation}
 
-\subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t=32})]{Deformation rate dependent viscosities (\protect\np{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t}\forcode{=32})}
+\subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t=32})]{Deformation rate dependent viscosities (\protect\np[=32]{nn_ahm_ijk_t}{nn\_ahm\_ijk\_t})}
 
 This option refers to the \citep{smagorinsky_MW63} scheme which is here implemented for momentum only. Smagorinsky chose as a
@@ -505,5 +507,5 @@
 }
 
-When  \citet{gent.mcwilliams_JPO90} diffusion is used (\np{ln_ldfeiv}{ln\_ldfeiv}\forcode{=.true.}),
+When  \citet{gent.mcwilliams_JPO90} diffusion is used (\np[=.true.]{ln_ldfeiv}{ln\_ldfeiv}),
 an eddy induced tracer advection term is added,
 the formulation of which depends on the slopes of iso-neutral surfaces.
@@ -512,5 +514,5 @@
 and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $s$-coordinates.
 
-If isopycnal mixing is used in the standard way, \ie\ \np{ln_traldf_triad}{ln\_traldf\_triad}\forcode{=.false.}, the eddy induced velocity is given by:
+If isopycnal mixing is used in the standard way, \ie\ \np[=.false.]{ln_traldf_triad}{ln\_traldf\_triad}, the eddy induced velocity is given by:
 \begin{equation}
   \label{eq:LDF_eiv}
@@ -536,5 +538,5 @@
 \colorbox{yellow}{CASE \np{nn_aei_ijk_t}{nn\_aei\_ijk\_t} = 21 to be added}
 
-In case of setting \np{ln_traldf_triad}{ln\_traldf\_triad}\forcode{ = .true.}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:TRIADS}.
+In case of setting \np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:TRIADS}.
 
 % ================================================================
@@ -554,11 +556,11 @@
 %--------------------------------------------------------------------------------------------------------------
 
-If  \np{ln_mle}{ln\_mle}\forcode{=.true.} in \nam{tra_mle}{tra\_mle} namelist, a parameterization of the mixing due to unresolved mixed layer instabilities is activated (\citet{foxkemper.ferrari_JPO08}). Additional transport is computed in \rou{ldf\_mle\_trp} and added to the eulerian transport in \rou{tra\_adv} as done for eddy induced advection.
+If  \np[=.true.]{ln_mle}{ln\_mle} in \nam{tra_mle}{tra\_mle} namelist, a parameterization of the mixing due to unresolved mixed layer instabilities is activated (\citet{foxkemper.ferrari_JPO08}). Additional transport is computed in \rou{ldf\_mle\_trp} and added to the eulerian transport in \rou{tra\_adv} as done for eddy induced advection.
 
 \colorbox{yellow}{TBC}
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -569,12 +571,12 @@
 Values between 0 to 4 are associated with interpolation while values 5 or 6 are associated with averaging.
 \begin{itemize}
-\item \np{nn_2dint}{nn\_2dint}\forcode{ = 0}: Distance-weighted interpolation
-\item \np{nn_2dint}{nn\_2dint}\forcode{ = 1}: Distance-weighted interpolation (small angle)
-\item \np{nn_2dint}{nn\_2dint}\forcode{ = 2}: Bilinear interpolation (geographical grid)
-\item \np{nn_2dint}{nn\_2dint}\forcode{ = 3}: Bilinear remapping interpolation (general grid)
-\item \np{nn_2dint}{nn\_2dint}\forcode{ = 4}: Polynomial interpolation
-\item \np{nn_2dint}{nn\_2dint}\forcode{ = 5}: Radial footprint averaging with diameter specified in the namelist as
+\item \np[=0]{nn_2dint}{nn\_2dint}: Distance-weighted interpolation
+\item \np[=1]{nn_2dint}{nn\_2dint}: Distance-weighted interpolation (small angle)
+\item \np[=2]{nn_2dint}{nn\_2dint}: Bilinear interpolation (geographical grid)
+\item \np[=3]{nn_2dint}{nn\_2dint}: Bilinear remapping interpolation (general grid)
+\item \np[=4]{nn_2dint}{nn\_2dint}: Polynomial interpolation
+\item \np[=5]{nn_2dint}{nn\_2dint}: Radial footprint averaging with diameter specified in the namelist as
   \texttt{rn\_[var]\_avglamscl} in degrees or metres (set using \texttt{ln\_[var]\_fp\_indegs})
-\item \np{nn_2dint}{nn\_2dint}\forcode{ = 6}: Rectangular footprint averaging with E/W and N/S size specified in
+\item \np[=6]{nn_2dint}{nn\_2dint}: Rectangular footprint averaging with E/W and N/S size specified in
   the namelist as \texttt{rn\_[var]\_avglamscl} and \texttt{rn\_[var]\_avgphiscl} in degrees or metres
   (set using \texttt{ln\_[var]\_fp\_indegs})
@@ -1168,7 +1170,7 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -41,12 +43,12 @@
 \begin{itemize}
 \item
-  a bulk formulation (\np{ln_blk}{ln\_blk}\forcode{=.true.} with four possible bulk algorithms),
-\item
-  a flux formulation (\np{ln_flx}{ln\_flx}\forcode{=.true.}),
+  a bulk formulation (\np[=.true.]{ln_blk}{ln\_blk} with four possible bulk algorithms),
+\item
+  a flux formulation (\np[=.true.]{ln_flx}{ln\_flx}),
 \item
   a coupled or mixed forced/coupled formulation (exchanges with a atmospheric model via the OASIS coupler),
-(\np{ln_cpl}{ln\_cpl} or \np{ln_mixcpl}{ln\_mixcpl}\forcode{=.true.}),
-\item
-  a user defined formulation (\np{ln_usr}{ln\_usr}\forcode{=.true.}).
+(\np{ln_cpl}{ln\_cpl} or \np[=.true.]{ln_mixcpl}{ln\_mixcpl}),
+\item
+  a user defined formulation (\np[=.true.]{ln_usr}{ln\_usr}).
 \end{itemize}
 
@@ -69,41 +71,41 @@
   the local grid directions in the model,
 \item
-  the use of a land/sea mask for input fields (\np{nn_lsm}{nn\_lsm}\forcode{=.true.}),
-\item
-  the addition of a surface restoring term to observed SST and/or SSS (\np{ln_ssr}{ln\_ssr}\forcode{=.true.}),
+  the use of a land/sea mask for input fields (\np[=.true.]{nn_lsm}{nn\_lsm}),
+\item
+  the addition of a surface restoring term to observed SST and/or SSS (\np[=.true.]{ln_ssr}{ln\_ssr}),
 \item
   the modification of fluxes below ice-covered areas (using climatological ice-cover or a sea-ice model)
-  (\np{nn_ice}{nn\_ice}\forcode{=0..3}),
-\item
-  the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln_rnf}{ln\_rnf}\forcode{=.true.}),
+  (\np[=0..3]{nn_ice}{nn\_ice}),
+\item
+  the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np[=.true.]{ln_rnf}{ln\_rnf}),
 \item
   the addition of ice-shelf melting as lateral inflow (parameterisation) or
-  as fluxes applied at the land-ice ocean interface (\np{ln_isf}{ln\_isf}\forcode{=.true.}),
+  as fluxes applied at the land-ice ocean interface (\np[=.true.]{ln_isf}{ln\_isf}),
 \item
   the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift
-  (\np{nn_fwb}{nn\_fwb}\forcode{=0..2}),
+  (\np[=0..2]{nn_fwb}{nn\_fwb}),
 \item
   the transformation of the solar radiation (if provided as daily mean) into an analytical diurnal cycle
-  (\np{ln_dm2dc}{ln\_dm2dc}\forcode{=.true.}),
-\item
-  the activation of wave effects from an external wave model  (\np{ln_wave}{ln\_wave}\forcode{=.true.}),
-\item
-  a neutral drag coefficient is read from an external wave model (\np{ln_cdgw}{ln\_cdgw}\forcode{=.true.}),
-\item
-  the Stokes drift from an external wave model is accounted for (\np{ln_sdw}{ln\_sdw}\forcode{=.true.}),
-\item
-  the choice of the Stokes drift profile parameterization (\np{nn_sdrift}{nn\_sdrift}\forcode{=0..2}),
-\item
-  the surface stress given to the ocean is modified by surface waves (\np{ln_tauwoc}{ln\_tauwoc}\forcode{=.true.}),
-\item
-  the surface stress given to the ocean is read from an external wave model (\np{ln_tauw}{ln\_tauw}\forcode{=.true.}),
-\item
-  the Stokes-Coriolis term is included (\np{ln_stcor}{ln\_stcor}\forcode{=.true.}),
-\item
-  the light penetration in the ocean (\np{ln_traqsr}{ln\_traqsr}\forcode{=.true.} with namelist \nam{tra_qsr}{tra\_qsr}),
-\item
-  the atmospheric surface pressure gradient effect on ocean and ice dynamics (\np{ln_apr_dyn}{ln\_apr\_dyn}\forcode{=.true.} with namelist \nam{sbc_apr}{sbc\_apr}),
-\item
-  the effect of sea-ice pressure on the ocean (\np{ln_ice_embd}{ln\_ice\_embd}\forcode{=.true.}).
+  (\np[=.true.]{ln_dm2dc}{ln\_dm2dc}),
+\item
+  the activation of wave effects from an external wave model  (\np[=.true.]{ln_wave}{ln\_wave}),
+\item
+  a neutral drag coefficient is read from an external wave model (\np[=.true.]{ln_cdgw}{ln\_cdgw}),
+\item
+  the Stokes drift from an external wave model is accounted for (\np[=.true.]{ln_sdw}{ln\_sdw}),
+\item
+  the choice of the Stokes drift profile parameterization (\np[=0..2]{nn_sdrift}{nn\_sdrift}),
+\item
+  the surface stress given to the ocean is modified by surface waves (\np[=.true.]{ln_tauwoc}{ln\_tauwoc}),
+\item
+  the surface stress given to the ocean is read from an external wave model (\np[=.true.]{ln_tauw}{ln\_tauw}),
+\item
+  the Stokes-Coriolis term is included (\np[=.true.]{ln_stcor}{ln\_stcor}),
+\item
+  the light penetration in the ocean (\np[=.true.]{ln_traqsr}{ln\_traqsr} with namelist \nam{tra_qsr}{tra\_qsr}),
+\item
+  the atmospheric surface pressure gradient effect on ocean and ice dynamics (\np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} with namelist \nam{sbc_apr}{sbc\_apr}),
+\item
+  the effect of sea-ice pressure on the ocean (\np[=.true.]{ln_ice_embd}{ln\_ice\_embd}).
 \end{itemize}
 
@@ -142,5 +144,5 @@
 The latter is the penetrative part of the heat flux.
 It is applied as a 3D trend of the temperature equation (\mdl{traqsr} module) when
-\np{ln_traqsr}{ln\_traqsr}\forcode{=.true.}.
+\np[=.true.]{ln_traqsr}{ln\_traqsr}.
 The way the light penetrates inside the water column is generally a sum of decreasing exponentials
 (see \autoref{subsec:TRA_qsr}).
@@ -278,7 +280,7 @@
                                   &  daily or weekLL     &  monthly           &  yearly        \\
       \hline
-      \np{clim}{clim}\forcode{=.false.} &  fn\_yYYYYmMMdDD.nc  &  fn\_yYYYYmMM.nc   &  fn\_yYYYY.nc  \\
+      \np[=.false.]{clim}{clim} &  fn\_yYYYYmMMdDD.nc  &  fn\_yYYYYmMM.nc   &  fn\_yYYYY.nc  \\
       \hline
-      \np{clim}{clim}\forcode{=.true.}  &  not possible        &  fn\_m??.nc        &  fn            \\
+      \np[=.true.]{clim}{clim}  &  not possible        &  fn\_m??.nc        &  fn            \\
       \hline
     \end{tabular}
@@ -351,5 +353,5 @@
 However, for forcing data related to the surface module,
 values are not needed at every time-step but at every \np{nn_fsbc}{nn\_fsbc} time-step.
-For example with \np{nn_fsbc}{nn\_fsbc}\forcode{=3}, the surface module will be called at time-steps 1, 4, 7, etc.
+For example with \np[=3]{nn_fsbc}{nn\_fsbc}, the surface module will be called at time-steps 1, 4, 7, etc.
 The date used for the time interpolation is thus redefined to the middle of \np{nn_fsbc}{nn\_fsbc} time-step period.
 In the previous example, this leads to: 1h30'00", 4h30'00", 7h30'00", etc. \\
@@ -550,5 +552,5 @@
   Spinup of the iceberg floats
 \item
-  Ocean/sea-ice simulation with both models running in parallel (\np{ln_mixcpl}{ln\_mixcpl}\forcode{=.true.})
+  Ocean/sea-ice simulation with both models running in parallel (\np[=.true.]{ln_mixcpl}{ln\_mixcpl})
 \end{itemize}
 
@@ -605,5 +607,5 @@
 
 The user can also choose in the \nam{sbc_sas}{sbc\_sas} namelist to read the mean (nn\_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level using
- (\np{ln_flx}{ln\_flx}\forcode{=.true.}) and to provide 3D oceanic velocities instead of 2D ones (\np{ln_flx}{ln\_flx}\forcode{=.true.}). In that last case, only the 1st level will be read in.
+ (\np[=.true.]{ln_flx}{ln\_flx}) and to provide 3D oceanic velocities instead of 2D ones (\np{ln_flx}{ln\_flx}\forcode{=.true.}). In that last case, only the 1st level will be read in.
 
 
@@ -623,5 +625,5 @@
 %-------------------------------------------------------------------------------------------------------------
 
-In the flux formulation (\np{ln_flx}{ln\_flx}\forcode{=.true.}),
+In the flux formulation (\np[=.true.]{ln_flx}{ln\_flx}),
 the surface boundary condition fields are directly read from input files.
 The user has to define in the namelist \nam{sbc_flx}{sbc\_flx} the name of the file,
@@ -731,5 +733,5 @@
 \begin{itemize}
 \item
-  NCAR (\np{ln_NCAR}{ln\_NCAR}\forcode{=.true.}):
+  NCAR (\np[=.true.]{ln_NCAR}{ln\_NCAR}):
   The NCAR bulk formulae have been developed by \citet{large.yeager_rpt04}.
   They have been designed to handle the NCAR forcing, a mixture of NCEP reanalysis and satellite data.
@@ -741,11 +743,11 @@
   This is the so-called DRAKKAR Forcing Set (DFS) \citep{brodeau.barnier.ea_OM10}.
 \item
-  COARE 3.0 (\np{ln_COARE_3p0}{ln\_COARE\_3p0}\forcode{=.true.}):
+  COARE 3.0 (\np[=.true.]{ln_COARE_3p0}{ln\_COARE\_3p0}):
   See \citet{fairall.bradley.ea_JC03} for more details
 \item
-  COARE 3.5 (\np{ln_COARE_3p5}{ln\_COARE\_3p5}\forcode{=.true.}):
+  COARE 3.5 (\np[=.true.]{ln_COARE_3p5}{ln\_COARE\_3p5}):
   See \citet{edson.jampana.ea_JPO13} for more details
 \item
-  ECMWF (\np{ln_ECMWF}{ln\_ECMWF}\forcode{=.true.}):
+  ECMWF (\np[=.true.]{ln_ECMWF}{ln\_ECMWF}):
   Based on \href{https://www.ecmwf.int/node/9221}{IFS (Cy31)} implementation and documentation.
   Surface roughness lengths needed for the Obukhov length are computed following \citet{beljaars_QJRMS95}.
@@ -762,8 +764,8 @@
 \begin{itemize}
 \item
-  Constant value (\np{constant value}{constant\ value}\forcode{ Cd_ice = 1.4e-3 }):
+  Constant value (\np[ Cd_ice=1.4e-3 ]{constant value}{constant\ value}):
   default constant value used for momentum and heat neutral transfer coefficients
 \item
-  \citet{lupkes.gryanik.ea_JGR12} (\np{ln_Cd_L12}{ln\_Cd\_L12}\forcode{=.true.}):
+  \citet{lupkes.gryanik.ea_JGR12} (\np[=.true.]{ln_Cd_L12}{ln\_Cd\_L12}):
   This scheme adds a dependency on edges at leads, melt ponds and flows
   of the constant neutral air-ice drag. After some approximations,
@@ -773,5 +775,5 @@
   It is theoretically applicable to all ice conditions (not only MIZ).
 \item
-  \citet{lupkes.gryanik_JGR15} (\np{ln_Cd_L15}{ln\_Cd\_L15}\forcode{=.true.}):
+  \citet{lupkes.gryanik_JGR15} (\np[=.true.]{ln_Cd_L15}{ln\_Cd\_L15}):
   Alternative turbulent transfer coefficients formulation between sea-ice
   and atmosphere with distinct momentum and heat coefficients depending
@@ -842,5 +844,5 @@
 
 The optional atmospheric pressure can be used to force ocean and ice dynamics
-(\np{ln_apr_dyn}{ln\_apr\_dyn}\forcode{=.true.}, \nam{sbc}{sbc} namelist).
+(\np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn}, \nam{sbc}{sbc} namelist).
 The input atmospheric forcing defined via \np{sn_apr}{sn\_apr} structure (\nam{sbc_apr}{sbc\_apr} namelist)
 can be interpolated in time to the model time step, and even in space when the interpolation on-the-fly is used.
@@ -914,6 +916,6 @@
 computationally too expensive. Here, two options are available:
 $\Pi_{sal}$ generated by an external model can be read in
-(\np{ln_read_load}{ln\_read\_load}\forcode{ =.true.}), or a ``scalar approximation'' can be
-used (\np{ln_scal_load}{ln\_scal\_load}\forcode{ =.true.}). In the latter case
+(\np[=.true.]{ln_read_load}{ln\_read\_load}), or a ``scalar approximation'' can be
+used (\np[=.true.]{ln_scal_load}{ln\_scal\_load}). In the latter case
 \[
   \Pi_{sal} = \beta \eta,
@@ -1076,6 +1078,6 @@
 \begin{description}
 
-  \item[\np{nn_isf}{nn\_isf}\forcode{=1}]:
-  The ice shelf cavity is represented (\np{ln_isfcav}{ln\_isfcav}\forcode{=.true.} needed).
+  \item[{\np[=1]{nn_isf}{nn\_isf}}]:
+  The ice shelf cavity is represented (\np[=.true.]{ln_isfcav}{ln\_isfcav} needed).
   The fwf and heat flux are depending of the local water properties.
 
@@ -1083,8 +1085,8 @@
 
 	\begin{description}
-	\item[\np{nn_isfblk}{nn\_isfblk}\forcode{=1}]:
+	\item[{\np[=1]{nn_isfblk}{nn\_isfblk}}]:
      The melt rate is based on a balance between the upward ocean heat flux and
      the latent heat flux at the ice shelf base. A complete description is available in \citet{hunter_rpt06}.
-	\item[\np{nn_isfblk}{nn\_isfblk}\forcode{=2}]:
+	\item[{\np[=2]{nn_isfblk}{nn\_isfblk}}]:
      The melt rate and the heat flux are based on a 3 equations formulation
      (a heat flux budget at the ice base, a salt flux budget at the ice base and a linearised freezing point temperature equation).
@@ -1103,5 +1105,5 @@
      There are 3 different ways to compute the exchange coeficient:
 	\begin{description}
-        \item[\np{nn_gammablk}{nn\_gammablk}\forcode{=0}]:
+        \item[{\np[=0]{nn_gammablk}{nn\_gammablk}}]:
      The salt and heat exchange coefficients are constant and defined by \np{rn_gammas0}{rn\_gammas0} and \np{rn_gammat0}{rn\_gammat0}.
      \begin{gather*}
@@ -1111,5 +1113,5 @@
      \end{gather*}
      This is the recommended formulation for ISOMIP.
-	\item[\np{nn_gammablk}{nn\_gammablk}\forcode{=1}]:
+	\item[{\np[=1]{nn_gammablk}{nn\_gammablk}}]:
      The salt and heat exchange coefficients are velocity dependent and defined as
      \begin{gather*}
@@ -1119,5 +1121,5 @@
      where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters).
      See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. It is the recommended formulation for realistic application.
-	\item[\np{nn_gammablk}{nn\_gammablk}\forcode{=2}]:
+	\item[{\np[=2]{nn_gammablk}{nn\_gammablk}}]:
      The salt and heat exchange coefficients are velocity and stability dependent and defined as:
 \[
@@ -1130,12 +1132,12 @@
      This formulation has not been extensively tested in \NEMO\ (not recommended).
    \end{description}
-  \item[\np{nn_isf}{nn\_isf}\forcode{=2}]:
+  \item[{\np[=2]{nn_isf}{nn\_isf}}]:
    The ice shelf cavity is not represented.
    The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting.
    The fluxes are distributed along the ice shelf edge between the depth of the average grounding line (GL)
    (\np{sn_depmax_isf}{sn\_depmax\_isf}) and the base of the ice shelf along the calving front
-   (\np{sn_depmin_isf}{sn\_depmin\_isf}) as in (\np{nn_isf}{nn\_isf}\forcode{=3}).
+   (\np{sn_depmin_isf}{sn\_depmin\_isf}) as in (\np[=3]{nn_isf}{nn\_isf}).
    The effective melting length (\np{sn_Leff_isf}{sn\_Leff\_isf}) is read from a file.
-  \item[\np{nn_isf}{nn\_isf}\forcode{=3}]:
+  \item[{\np[=3]{nn_isf}{nn\_isf}}]:
    The ice shelf cavity is not represented.
    The fwf (\np{sn_rnfisf}{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between
@@ -1143,17 +1145,17 @@
    the base of the ice shelf along the calving front (\np{sn_depmin_isf}{sn\_depmin\_isf}).
    The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.
-  \item[\np{nn_isf}{nn\_isf}\forcode{=4}]:
-   The ice shelf cavity is opened (\np{ln_isfcav}{ln\_isfcav}\forcode{=.true.} needed).
+  \item[{\np[=4]{nn_isf}{nn\_isf}}]:
+   The ice shelf cavity is opened (\np[=.true.]{ln_isfcav}{ln\_isfcav} needed).
    However, the fwf is not computed but specified from file \np{sn_fwfisf}{sn\_fwfisf}).
    The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.
-   As in \np{nn_isf}{nn\_isf}\forcode{=1}, the fluxes are spread over the top boundary layer thickness (\np{rn_hisf_tbl}{rn\_hisf\_tbl})\\
+   As in \np[=1]{nn_isf}{nn\_isf}, the fluxes are spread over the top boundary layer thickness (\np{rn_hisf_tbl}{rn\_hisf\_tbl})\\
 \end{description}
 
-$\bullet$ \np{nn_isf}{nn\_isf}\forcode{=1} and \np{nn_isf}{nn\_isf}\forcode{=2} compute a melt rate based on
+$\bullet$ \np[=1]{nn_isf}{nn\_isf} and \np[=2]{nn_isf}{nn\_isf} 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}{nn\_isf}\forcode{=3} and \np{nn_isf}{nn\_isf}\forcode{=4} read the melt rate from a file.
+$\bullet$ \np[=3]{nn_isf}{nn\_isf} and \np[=4]{nn_isf}{nn\_isf} read the melt rate from a file.
 You have total control of the fwf forcing.
 This can be useful if the water masses on the shelf are not realistic or
@@ -1204,5 +1206,5 @@
 \end{description}
 
-If \np{ln_iscpl}{ln\_iscpl}\forcode{=.true.}, the isf draft is assume to be different at each restart step with
+If \np[=.true.]{ln_iscpl}{ln\_iscpl}, the isf draft is assume to be different at each restart step with
 potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics.
 The wetting and drying scheme applied on the restart is very simple and described below for the 6 different possible cases:
@@ -1240,5 +1242,5 @@
 
 In order to remove the trend and keep the conservation level as close to 0 as possible,
-a simple conservation scheme is available with \np{ln_hsb}{ln\_hsb}\forcode{=.true.}.
+a simple conservation scheme is available with \np[=.true.]{ln_hsb}{ln\_hsb}.
 The heat/salt/vol. gain/loss is diagnosed, as well as the location.
 A correction increment is computed and apply each time step during the next \np{rn_fiscpl}{rn\_fiscpl} time steps.
@@ -1270,9 +1272,9 @@
 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}{ln\_icebergs}\forcode{=.true.}.
+They are enabled by setting \np[=.true.]{ln_icebergs}{ln\_icebergs}.
 
 Two initialisation schemes are possible.
 \begin{description}
-\item[\np{nn_test_icebergs}{nn\_test\_icebergs}~$>$~0]
+\item[{\np{nn_test_icebergs}{nn\_test\_icebergs}~$>$~0}]
   In this scheme, the value of \np{nn_test_icebergs}{nn\_test\_icebergs} represents the class of iceberg to generate
   (so between 1 and 10), and \np{nn_test_icebergs}{nn\_test\_icebergs} provides a lon/lat box in the domain at each grid point of
@@ -1281,5 +1283,5 @@
   \np{nn_test_icebergs}{nn\_test\_icebergs} is defined by four numbers in \np{nn_test_box}{nn\_test\_box} representing the corners of
   the geographical box: lonmin,lonmax,latmin,latmax
-\item[\np{nn_test_icebergs}{nn\_test\_icebergs}\forcode{=-1}]
+\item[{\np[=-1]{nn_test_icebergs}{nn\_test\_icebergs}}]
   In this scheme, the model reads a calving file supplied in the \np{sn_icb}{sn\_icb} parameter.
   This should be a file with a field on the configuration grid (typically ORCA)
@@ -1306,8 +1308,8 @@
 The amount of information is controlled by two integer parameters:
 \begin{description}
-\item[\np{nn_verbose_level}{nn\_verbose\_level}] takes a value between one and four and
+\item[{\np{nn_verbose_level}{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}{nn\_verbose\_write}] is the number of timesteps between writes
+\item[{\np{nn_verbose_write}{nn\_verbose\_write}}] is the number of timesteps between writes
 \end{description}
 
@@ -1343,5 +1345,5 @@
 
 Physical processes related to ocean surface waves can be accounted by setting the logical variable
-\np{ln_wave}{ln\_wave}\forcode{=.true.} in \nam{sbc}{sbc} namelist. In addition, specific flags accounting for
+\np[=.true.]{ln_wave}{ln\_wave} in \nam{sbc}{sbc} namelist. In addition, specific flags accounting for
 different processes should be activated as explained in the following sections.
 
@@ -1351,5 +1353,5 @@
 for external data names, locations, frequency, interpolation and all the miscellanous options allowed by
 Input Data generic Interface (see \autoref{sec:SBC_input}).
-\item[coupled mode]: \NEMO\ and an external wave model can be coupled by setting \np{ln_cpl}{ln\_cpl} \forcode{= .true.}
+\item[coupled mode]: \NEMO\ and an external wave model can be coupled by setting \np[=.true.]{ln_cpl}{ln\_cpl}
 in \nam{sbc}{sbc} namelist and filling the \nam{sbc_cpl}{sbc\_cpl} namelist.
 \end{description}
@@ -1364,5 +1366,5 @@
 
 The neutral surface drag coefficient provided from an external data source (\ie\ a wave model),
-can be used by setting the logical variable \np{ln_cdgw}{ln\_cdgw} \forcode{= .true.} in \nam{sbc}{sbc} namelist.
+can be used by setting the logical variable \np[=.true.]{ln_cdgw}{ln\_cdgw} in \nam{sbc}{sbc} namelist.
 Then using the routine \rou{sbcblk\_algo\_ncar} and starting from the neutral drag coefficent provided,
 the drag coefficient is computed according to the stable/unstable conditions of the
@@ -1408,5 +1410,5 @@
 
 \begin{description}
-\item[\np{nn_sdrift}{nn\_sdrift} = 0]: exponential integral profile parameterization proposed by
+\item[{\np{nn_sdrift}{nn\_sdrift} = 0}]: exponential integral profile parameterization proposed by
 \citet{breivik.janssen.ea_JPO14}:
 
@@ -1427,5 +1429,5 @@
 where $H_s$ is the significant wave height and $\omega$ is the wave frequency.
 
-\item[\np{nn_sdrift}{nn\_sdrift} = 1]: velocity profile based on the Phillips spectrum which is considered to be a
+\item[{\np{nn_sdrift}{nn\_sdrift} = 1}]: velocity profile based on the Phillips spectrum which is considered to be a
 reasonable estimate of the part of the spectrum mostly contributing to the Stokes drift velocity near the surface
 \citep{breivik.bidlot.ea_OM16}:
@@ -1439,5 +1441,5 @@
 where $erf$ is the complementary error function and $k_p$ is the peak wavenumber.
 
-\item[\np{nn_sdrift}{nn\_sdrift} = 2]: velocity profile based on the Phillips spectrum as for \np{nn_sdrift}{nn\_sdrift} = 1
+\item[{\np{nn_sdrift}{nn\_sdrift} = 2}]: velocity profile based on the Phillips spectrum as for \np{nn_sdrift}{nn\_sdrift} = 1
 but using the wave frequency from a wave model.
 
@@ -1477,5 +1479,5 @@
 In order to include this term, once evaluated the Stokes drift (using one of the 3 possible
 approximations described in \autoref{subsec:SBC_wave_sdw}),
-\np{ln_stcor}{ln\_stcor}\forcode{=.true.} has to be set.
+\np[=.true.]{ln_stcor}{ln\_stcor} has to be set.
 
 
@@ -1517,6 +1519,6 @@
 
 The wave stress derived from an external wave model can be provided either through the normalized
-wave stress into the ocean by setting \np{ln_tauwoc}{ln\_tauwoc}\forcode{=.true.}, or through the zonal and
-meridional stress components by setting \np{ln_tauw}{ln\_tauw}\forcode{=.true.}.
+wave stress into the ocean by setting \np[=.true.]{ln_tauwoc}{ln\_tauwoc}, or through the zonal and
+meridional stress components by setting \np[=.true.]{ln_tauw}{ln\_tauw}.
 
 
@@ -1561,7 +1563,7 @@
 assuming that the diurnal cycle of SWF is a scaling of the top of the atmosphere diurnal cycle of incident SWF.
 The \cite{bernie.guilyardi.ea_CD07} reconstruction algorithm is available in \NEMO\ by
-setting \np{ln_dm2dc}{ln\_dm2dc}\forcode{=.true.} (a \textit{\nam{sbc}{sbc}} namelist variable) when
-using a bulk formulation (\np{ln_blk}{ln\_blk}\forcode{=.true.}) or
-the flux formulation (\np{ln_flx}{ln\_flx}\forcode{=.true.}).
+setting \np[=.true.]{ln_dm2dc}{ln\_dm2dc} (a \textit{\nam{sbc}{sbc}} namelist variable) when
+using a bulk formulation (\np[=.true.]{ln_blk}{ln\_blk}) or
+the flux formulation (\np[=.true.]{ln_flx}{ln\_flx}).
 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.guilyardi.ea_CD07}.
@@ -1598,5 +1600,5 @@
 \label{subsec:SBC_rotation}
 
-When using a flux (\np{ln_flx}{ln\_flx}\forcode{=.true.}) or bulk (\np{ln_blk}{ln\_blk}\forcode{=.true.}) formulation,
+When using a flux (\np[=.true.]{ln_flx}{ln\_flx}) or bulk (\np[=.true.]{ln_blk}{ln\_blk}) 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
@@ -1627,5 +1629,5 @@
 
 Options are defined through the \nam{sbc_ssr}{sbc\_ssr} namelist variables.
-On forced mode using a flux formulation (\np{ln_flx}{ln\_flx}\forcode{=.true.}),
+On forced mode using a flux formulation (\np[=.true.]{ln_flx}{ln\_flx}),
 a feedback term \emph{must} be added to the surface heat flux $Q_{ns}^o$:
 \[
@@ -1746,12 +1748,12 @@
 
 \begin{description}
-\item[\np{nn_fwb}{nn\_fwb}\forcode{=0}]
+\item[{\np[=0]{nn_fwb}{nn\_fwb}}]
   no control at all.
   The mean sea level is free to drift, and will certainly do so.
-\item[\np{nn_fwb}{nn\_fwb}\forcode{=1}]
+\item[{\np[=1]{nn_fwb}{nn\_fwb}}]
   global mean \textit{emp} set to zero at each model time step.
   %GS: comment below still relevant ?
   %Note that with a sea-ice model, this technique only controls the mean sea level with linear free surface and no mass flux between ocean and ice (as it is implemented in the current ice-ocean coupling).
-\item[\np{nn_fwb}{nn\_fwb}\forcode{=2}]
+\item[{\np[=2]{nn_fwb}{nn\_fwb}}]
   freshwater budget is adjusted from the previous year annual mean budget which
   is read in the \textit{EMPave\_old.dat} file.
@@ -1783,7 +1785,7 @@
 
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -186,17 +188,17 @@
 
 \begin{description}
-\item[\np{nn_sto_eos}{nn\_sto\_eos}:]   number of independent random walks
-\item[\np{rn_eos_stdxy}{rn\_eos\_stdxy}:] random walk horizontal standard deviation (in grid points)
-\item[\np{rn_eos_stdz}{rn\_eos\_stdz}:]  random walk vertical standard deviation (in grid points)
-\item[\np{rn_eos_tcor}{rn\_eos\_tcor}:]  random walk time correlation (in timesteps)
-\item[\np{nn_eos_ord}{nn\_eos\_ord}:]   order of autoregressive processes
-\item[\np{nn_eos_flt}{nn\_eos\_flt}:]   passes of Laplacian filter
-\item[\np{rn_eos_lim}{rn\_eos\_lim}:]   limitation factor (default = 3.0)
+\item[{\np{nn_sto_eos}{nn\_sto\_eos}:}]   number of independent random walks
+\item[{\np{rn_eos_stdxy}{rn\_eos\_stdxy}:}] random walk horizontal standard deviation (in grid points)
+\item[{\np{rn_eos_stdz}{rn\_eos\_stdz}:}]  random walk vertical standard deviation (in grid points)
+\item[{\np{rn_eos_tcor}{rn\_eos\_tcor}:}]  random walk time correlation (in timesteps)
+\item[{\np{nn_eos_ord}{nn\_eos\_ord}:}]   order of autoregressive processes
+\item[{\np{nn_eos_flt}{nn\_eos\_flt}:}]   passes of Laplacian filter
+\item[{\np{rn_eos_lim}{rn\_eos\_lim}:}]   limitation factor (default = 3.0)
 \end{description}
 
 The first four parameters define the stochastic part of equation of state.
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -55,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}{ln\_tra\_trd} or \np{ln_tra_mxl}{ln\_tra\_mxl}\forcode{=.true.}), as described in \autoref{chap:DIA}.
+(\np{ln_tra_trd}{ln\_tra\_trd} or \np[=.true.]{ln_tra_mxl}{ln\_tra\_mxl}), as described in \autoref{chap:DIA}.
 
 % ================================================================
@@ -85,5 +87,5 @@
 Indeed, it is obtained by using the following equality: $\nabla \cdot (\vect U \, T) = \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, \ie\ \np{ln_linssh}{ln\_linssh}\forcode{=.true.}).
+(which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, \ie\ \np[=.true.]{ln_linssh}{ln\_linssh}).
 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 enforce the conservation properties of
@@ -121,5 +123,5 @@
 \begin{description}
 \item[linear free surface:]
-  (\np{ln_linssh}{ln\_linssh}\forcode{=.true.})
+  (\np[=.true.]{ln_linssh}{ln\_linssh})
   the first level thickness is constant in time:
   the vertical boundary condition is applied at the fixed surface $z = 0$ rather than on
@@ -129,5 +131,5 @@
   the first level tracer value.
 \item[non-linear free surface:]
-  (\np{ln_linssh}{ln\_linssh}\forcode{=.false.})
+  (\np[=.false.]{ln_linssh}{ln\_linssh})
   convergence/divergence in the first ocean level moves the free surface up/down.
   There is no tracer advection through it so that the advective fluxes through the surface are also zero.
@@ -190,5 +192,5 @@
 %        2nd order centred scheme
 
-The centred advection scheme (CEN) is used when \np{ln_traadv_cen}{ln\_traadv\_cen}\forcode{=.true.}.
+The centred advection scheme (CEN) is used when \np[=.true.]{ln_traadv_cen}{ln\_traadv\_cen}.
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by
 setting \np{nn_cen_h}{nn\_cen\_h} and \np{nn_cen_v}{nn\_cen\_v} to $2$ or $4$.
@@ -222,5 +224,5 @@
   \tau_u^{cen4} = \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \, \Big]}^{\,i + 1/2}
 \end{equation}
-In the vertical direction (\np{nn_cen_v}{nn\_cen\_v}\forcode{=4}),
+In the vertical direction (\np[=4]{nn_cen_v}{nn\_cen\_v}),
 a $4^{th}$ COMPACT interpolation has been prefered \citep{demange_phd14}.
 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion,
@@ -255,5 +257,5 @@
 \label{subsec:TRA_adv_tvd}
 
-The Flux Corrected Transport schemes (FCT) is used when \np{ln_traadv_fct}{ln\_traadv\_fct}\forcode{=.true.}.
+The Flux Corrected Transport schemes (FCT) is used when \np[=.true.]{ln_traadv_fct}{ln\_traadv\_fct}.
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by
 setting \np{nn_fct_h}{nn\_fct\_h} and \np{nn_fct_v}{nn\_fct\_v} to $2$ or $4$.
@@ -298,5 +300,5 @@
 \label{subsec:TRA_adv_mus}
 
-The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln_traadv_mus}{ln\_traadv\_mus}\forcode{=.true.}.
+The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np[=.true.]{ln_traadv_mus}{ln\_traadv\_mus}.
 MUSCL implementation can be found in the \mdl{traadv\_mus} module.
 
@@ -326,5 +328,5 @@
 This choice ensure 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}{ln\_mus\_ups}\forcode{=.true.}).
+(\np[=.true.]{ln_mus_ups}{ln\_mus\_ups}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -334,5 +336,5 @@
 \label{subsec:TRA_adv_ubs}
 
-The Upstream-Biased Scheme (UBS) is used when \np{ln_traadv_ubs}{ln\_traadv\_ubs}\forcode{=.true.}.
+The Upstream-Biased Scheme (UBS) is used when \np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs}.
 UBS implementation can be found in the \mdl{traadv\_mus} module.
 
@@ -364,5 +366,5 @@
 \citep{shchepetkin.mcwilliams_OM05, demange_phd14}.
 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme
-(\np{nn_ubs_v}{nn\_ubs\_v}\forcode{=2 or 4}).
+(\np[=2 or 4]{nn_ubs_v}{nn\_ubs\_v}).
 
 For stability reasons (see \autoref{chap:TD}), the first term  in \autoref{eq:TRA_adv_ubs}
@@ -407,5 +409,5 @@
 
 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme
-proposed by \citet{leonard_CMAME79} is used when \np{ln_traadv_qck}{ln\_traadv\_qck}\forcode{=.true.}.
+proposed by \citet{leonard_CMAME79} is used when \np[=.true.]{ln_traadv_qck}{ln\_traadv\_qck}.
 QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
 
@@ -452,5 +454,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 \autoref{chap:TD}).
-When \np{ln_traldf_msc}{ln\_traldf\_msc}\forcode{=.true.}, a Method of Stabilizing Correction is used in which
+When \np[=.true.]{ln_traldf_msc}{ln\_traldf\_msc}, a Method of Stabilizing Correction is used in which
 the pure vertical component is split into an explicit and an implicit part \citep{lemarie.debreu.ea_OM12}.
 
@@ -464,13 +466,13 @@
 
 \begin{description}
-\item[\np{ln_traldf_OFF}{ln\_traldf\_OFF}\forcode{=.true.}:]
+\item[{\np[=.true.]{ln_traldf_OFF}{ln\_traldf\_OFF}}]
   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}{ln\_traldf\_lap}\forcode{=.true.}:]
+\item[{\np[=.true.]{ln_traldf_lap}{ln\_traldf\_lap}}]
   a laplacian operator is selected.
   This harmonic operator takes the following expression:  $\mathcal{L}(T) = \nabla \cdot A_{ht} \; \nabla T $,
   where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}),
   and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}).
-\item[\np{ln_traldf_blp}{ln\_traldf\_blp}\forcode{=.true.}]:
+\item[{\np[=.true.]{ln_traldf_blp}{ln\_traldf\_blp}}]:
   a bilaplacian operator is selected.
   This biharmonic operator takes the following expression:
@@ -497,7 +499,7 @@
 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}{ln\_traldf\_lev}\forcode{=.true.}) or
+iso-level option is used (\np[=.true.]{ln_traldf_lev}{ln\_traldf\_lev}) or
 when a horizontal (\ie\ geopotential) operator is demanded in \textit{z}-coordinate
-(\np{ln_traldf_hor}{ln\_traldf\_hor} and \np{ln_zco}{ln\_zco}\forcode{=.true.}).
+(\np{ln_traldf_hor}{ln\_traldf\_hor} and \np[=.true.]{ln_zco}{ln\_zco}).
 The associated code can be found in the \mdl{traldf\_lap\_blp} module.
 The operator is a rotated (re-entrant) laplacian when
@@ -536,10 +538,10 @@
 It is a \textit{horizontal} operator (\ie 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}{ln\_traldf\_lap} or \np{ln_traldf_blp}{ln\_traldf\_blp}\forcode{=.true.},
-we have \np{ln_traldf_lev}{ln\_traldf\_lev}\forcode{=.true.} or \np{ln_traldf_hor}{ln\_traldf\_hor}~=~\np{ln_zco}{ln\_zco}\forcode{=.true.}.
+It is thus used when, in addition to \np{ln_traldf_lap}{ln\_traldf\_lap} or \np[=.true.]{ln_traldf_blp}{ln\_traldf\_blp},
+we have \np[=.true.]{ln_traldf_lev}{ln\_traldf\_lev} or \np{ln_traldf_hor}{ln\_traldf\_hor}~=~\np[=.true.]{ln_zco}{ln\_zco}.
 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}{ln\_zps}\forcode{=.true.}),
+Note that in the partial step $z$-coordinate (\np[=.true.]{ln_zps}{ln\_zps}),
 tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom.
 In this case, horizontal derivatives in (\autoref{eq:TRA_ldf_lap}) at the bottom level require a specific treatment.
@@ -573,7 +575,7 @@
 $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and
 the surface along which the diffusion operator acts (\ie\ horizontal or iso-neutral surfaces).
-It is thus used when, in addition to \np{ln_traldf_lap}{ln\_traldf\_lap}\forcode{=.true.},
-we have \np{ln_traldf_iso}{ln\_traldf\_iso}\forcode{=.true.},
-or both \np{ln_traldf_hor}{ln\_traldf\_hor}\forcode{=.true.} and \np{ln_zco}{ln\_zco}\forcode{=.true.}.
+It is thus used when, in addition to \np[=.true.]{ln_traldf_lap}{ln\_traldf\_lap},
+we have \np[=.true.]{ln_traldf_iso}{ln\_traldf\_iso},
+or both \np[=.true.]{ln_traldf_hor}{ln\_traldf\_hor} and \np[=.true.]{ln_zco}{ln\_zco}.
 The way these slopes are evaluated is given in \autoref{sec:LDF_slp}.
 At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using
@@ -591,5 +593,5 @@
 any additional background horizontal diffusion \citep{guilyardi.madec.ea_CD01}.
 
-Note that in the partial step $z$-coordinate (\np{ln_zps}{ln\_zps}\forcode{=.true.}),
+Note that in the partial step $z$-coordinate (\np[=.true.]{ln_zps}{ln\_zps}),
 the horizontal derivatives at the bottom level in \autoref{eq:TRA_ldf_iso} require a specific treatment.
 They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}.
@@ -601,5 +603,5 @@
 
 An alternative scheme developed by \cite{griffies.gnanadesikan.ea_JPO98} which ensures tracer variance decreases
-is also available in \NEMO\ (\np{ln_traldf_triad}{ln\_traldf\_triad}\forcode{=.true.}).
+is also available in \NEMO\ (\np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad}).
 A complete description of the algorithm is given in \autoref{apdx:TRIADS}.
 
@@ -647,5 +649,5 @@
 respectively.
 Generally, $A_w^{vT} = A_w^{vS}$ except when double diffusive mixing is parameterised
-(\ie\ \np{ln_zdfddm}{ln\_zdfddm}\forcode{=.true.},).
+(\ie\ \np[=.true.]{ln_zdfddm}{ln\_zdfddm},).
 The way these coefficients are evaluated is given in \autoref{chap:ZDF} (ZDF).
 Furthermore, when iso-neutral mixing is used, both mixing coefficients are increased by
@@ -722,5 +724,5 @@
 Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:TD}).
 
-In the linear free surface case (\np{ln_linssh}{ln\_linssh}\forcode{=.true.}), an additional term has to be added on
+In the linear free surface case (\np[=.true.]{ln_linssh}{ln\_linssh}), an additional term has to be added on
 both temperature and salinity.
 On temperature, this term remove the heat content associated with mass exchange that has been added to $Q_{ns}$.
@@ -757,7 +759,7 @@
 
 Options are defined through the \nam{tra_qsr}{tra\_qsr} namelist variables.
-When the penetrative solar radiation option is used (\np{ln_traqsr}{ln\_traqsr}\forcode{=.true.}),
+When the penetrative solar radiation option is used (\np[=.true.]{ln_traqsr}{ln\_traqsr}),
 the solar radiation penetrates the top few tens of meters of the ocean.
-If it is not used (\np{ln_traqsr}{ln\_traqsr}\forcode{=.false.}) all the heat flux is absorbed in the first ocean level.
+If it is not used (\np[=.false.]{ln_traqsr}{ln\_traqsr}) 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 \autoref{eq:MB_PE_tra_T} and
 the surface boundary condition is modified to take into account only the non-penetrative part of the surface
@@ -788,5 +790,5 @@
 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}{ln\_qsr\_2bd}\forcode{=.true.})
+In the simple 2-waveband light penetration scheme (\np[=.true.]{ln_qsr_2bd}{ln\_qsr\_2bd})
 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths,
 leading to the following expression \citep{paulson.simpson_JPO77}:
@@ -816,5 +818,5 @@
 The 2-bands formulation does not reproduce the full model very well.
 
-The RGB formulation is used when \np{ln_qsr_rgb}{ln\_qsr\_rgb}\forcode{=.true.}.
+The RGB formulation is used when \np[=.true.]{ln_qsr_rgb}{ln\_qsr\_rgb}.
 The RGB attenuation coefficients (\ie\ the inverses of the extinction length scales) are tabulated over
 61 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L
@@ -823,13 +825,13 @@
 
 \begin{description}
-\item[\np{nn_chldta}{nn\_chldta}\forcode{=0}]
+\item[{\np[=0]{nn_chldta}{nn\_chldta}}]
   a constant 0.05 g.Chl/L value everywhere ;
-\item[\np{nn_chldta}{nn\_chldta}\forcode{=1}]
+\item[{\np[=1]{nn_chldta}{nn\_chldta}}]
   an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in
   the vertical direction;
-\item[\np{nn_chldta}{nn\_chldta}\forcode{=2}]
+\item[{\np[=2]{nn_chldta}{nn\_chldta}}]
   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}{ln\_qsr\_bio}\forcode{=.true.}]
+\item[{\np[=.true.]{ln_qsr_bio}{ln\_qsr\_bio}}]
   simulated time varying chlorophyll by TOP biogeochemical model.
   In this case, the RGB formulation is used to calculate both the phytoplankton light limitation in
@@ -944,8 +946,8 @@
 %        Diffusive BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Diffusive bottom boundary layer (\forcode{nn_bbl_ldf=1})]{Diffusive bottom boundary layer (\protect\np{nn_bbl_ldf}{nn\_bbl\_ldf}\forcode{=1})}
+\subsection[Diffusive bottom boundary layer (\forcode{nn_bbl_ldf=1})]{Diffusive bottom boundary layer (\protect\np[=1]{nn_bbl_ldf}{nn\_bbl\_ldf})}
 \label{subsec:TRA_bbl_diff}
 
-When applying sigma-diffusion (\np{ln_trabbl}{ln\_trabbl}\forcode{=.true.} and \np{nn_bbl_ldf}{nn\_bbl\_ldf} set to 1),
+When applying sigma-diffusion (\np[=.true.]{ln_trabbl}{ln\_trabbl} and \np{nn_bbl_ldf}{nn\_bbl\_ldf} set to 1),
 the diffusive flux between two adjacent cells at the ocean floor is given by
 \[
@@ -983,5 +985,5 @@
 %        Advective BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Advective bottom boundary layer (\forcode{nn_bbl_adv=1,2})]{Advective bottom boundary layer (\protect\np{nn_bbl_adv}{nn\_bbl\_adv}\forcode{=1,2})}
+\subsection[Advective bottom boundary layer (\forcode{nn_bbl_adv=1,2})]{Advective bottom boundary layer (\protect\np[=1,2]{nn_bbl_adv}{nn\_bbl\_adv})}
 \label{subsec:TRA_bbl_adv}
 
@@ -1014,9 +1016,9 @@
 %%%gmcomment   :  this section has to be really written
 
-When applying an advective BBL (\np{nn_bbl_adv}{nn\_bbl\_adv}\forcode{=1..2}), an overturning circulation is added which
+When applying an advective BBL (\np[=1..2]{nn_bbl_adv}{nn\_bbl\_adv}), 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}{nn\_bbl\_adv}\forcode{=1}:
+\np[=1]{nn_bbl_adv}{nn\_bbl\_adv}:
 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step
 (see black arrow in \autoref{fig:TRA_bbl}) \citep{beckmann.doscher_JPO97}.
@@ -1025,5 +1027,5 @@
 if the velocity is directed towards greater depth (\ie\ $\vect U \cdot \nabla H > 0$).
 
-\np{nn_bbl_adv}{nn\_bbl\_adv}\forcode{=2}:
+\np[=2]{nn_bbl_adv}{nn\_bbl\_adv}:
 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_T99}.
@@ -1153,5 +1155,5 @@
 (\ie\ fluxes plus content in mass exchanges).
 $\gamma$ is initialized as \np{rn_atfp}{rn\_atfp} (\textbf{namelist} parameter).
-Its default value is \np{rn_atfp}{rn\_atfp}\forcode{=10.e-3}.
+Its default value is \np[=10.e-3]{rn_atfp}{rn\_atfp}.
 Note that the forcing correction term in the filter is not applied in linear free surface
 (\jp{ln\_linssh}\forcode{=.true.}) (see \autoref{subsec:TRA_sbc}).
@@ -1216,5 +1218,5 @@
 
 \begin{description}
-\item[\np{ln_teos10}{ln\_teos10}\forcode{=.true.}]
+\item[{\np[=.true.]{ln_teos10}{ln\_teos10}}]
   the polyTEOS10-bsq equation of seawater \citep{roquet.madec.ea_OM15} is used.
   The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,
@@ -1235,5 +1237,5 @@
   either computing the air-sea and ice-sea fluxes (forced mode) or
   sending the SST field to the atmosphere (coupled mode).
-\item[\np{ln_eos80}{ln\_eos80}\forcode{=.true.}]
+\item[{\np[=.true.]{ln_eos80}{ln\_eos80}}]
   the polyEOS80-bsq equation of seawater is used.
   It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized to
@@ -1247,5 +1249,5 @@
   Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which
   is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value.
-\item[\np{ln_seos}{ln\_seos}\forcode{=.true.}]
+\item[{\np[=.true.]{ln_seos}{ln\_seos}}]
   a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is chosen,
   the coefficients of which has been optimized to fit the behavior of TEOS10
@@ -1367,9 +1369,9 @@
 I've changed "derivative" to "difference" and "mean" to "average"}
 
-With partial cells (\np{ln_zps}{ln\_zps}\forcode{=.true.}) at bottom and top (\np{ln_isfcav}{ln\_isfcav}\forcode{=.true.}),
+With partial cells (\np[=.true.]{ln_zps}{ln\_zps}) at bottom and top (\np[=.true.]{ln_isfcav}{ln\_isfcav}),
 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}{ln\_isfcav}\forcode{=.true.}) are computed in the same way as
+The partial cell properties at the top (\np[=.true.]{ln_isfcav}{ln\_isfcav}) are computed in the same way as
 for the bottom.
 So, only the bottom interpolation is explained below.
@@ -1387,5 +1389,5 @@
   the $z$-partial step coordinate]{
     Discretisation of the horizontal difference and average of tracers in
-    the $z$-partial step coordinate (\protect\np{ln_zps}{ln\_zps}\forcode{=.true.}) in
+    the $z$-partial step coordinate (\protect\np[=.true.]{ln_zps}{ln\_zps}) in
     the case $(e3w_k^{i + 1} - e3w_k^i) > 0$.
     A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$,
@@ -1459,7 +1461,7 @@
 %%%
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 11582)
@@ -3,4 +3,6 @@
 %% Custom aliases
 \newcommand{\cf}{\ensuremath{C\kern-0.14em f}}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -42,6 +44,6 @@
 are computed and added to the 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}{ln\_zdfexp}\forcode{=.true.}) or a backward time stepping scheme
-%(\np{ln_zdfexp}{ln\_zdfexp}\forcode{=.false.}) depending on the magnitude of the mixing coefficients,
+%(namelist parameter \np[=.true.]{ln_zdfexp}{ln\_zdfexp}) or a backward time stepping scheme
+%(\np[=.false.]{ln_zdfexp}{ln\_zdfexp}) depending on the magnitude of the mixing coefficients,
 %and thus of the formulation used (see \autoref{chap:TD}).
 
@@ -92,5 +94,5 @@
 %--------------------------------------------------------------------------------------------------------------
 
-When \np{ln_zdfric}{ln\_zdfric}\forcode{=.true.}, a local Richardson number dependent formulation for the vertical momentum and
+When \np[=.true.]{ln_zdfric}{ln\_zdfric}, a local Richardson number dependent formulation for the vertical momentum and
 tracer eddy coefficients is set through the \nam{zdf_ric}{zdf\_ric} namelist variables.
 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model.
@@ -118,5 +120,5 @@
 
 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}{ln\_mldw}\forcode{=.true.} in the namelist.
+(wind-stress and buoyancy fluxes) can be activated setting the \np[=.true.]{ln_mldw}{ln\_mldw} in the namelist.
 
 In this case, the local depth of turbulent wind-mixing or "Ekman depth" $h_{e}(x,y,t)$ is evaluated and
@@ -225,9 +227,9 @@
 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}{nn\_mxl}\forcode{=0}) or by the local vertical scale factor (\np{nn_mxl}{nn\_mxl}\forcode{=1}).
+(\np[=0]{nn_mxl}{nn\_mxl}) or by the local vertical scale factor (\np[=1]{nn_mxl}{nn\_mxl}).
 \citet{blanke.delecluse_JPO93} 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{madec.delecluse.ea_NPM98} introduces the \np{nn_mxl}{nn\_mxl}\forcode{=2, 3} cases,
+To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np[=2, 3]{nn_mxl}{nn\_mxl} cases,
 which add an extra assumption concerning the vertical gradient of the computed length scale.
 So, the length scales are first evaluated as in \autoref{eq:ZDF_tke_mxl0_1} and then bounded such that:
@@ -267,6 +269,6 @@
 where $l^{(k)}$ is computed using \autoref{eq:ZDF_tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
 
-In the \np{nn_mxl}{nn\_mxl}\forcode{=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}{nn\_mxl}\forcode{=3} case,
+In the \np[=2]{nn_mxl}{nn\_mxl} 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[=3]{nn_mxl}{nn\_mxl} case,
 the dissipation and mixing turbulent length scales are give as in \citet{gaspar.gregoris.ea_JGR90}:
 \[
@@ -312,7 +314,7 @@
 $\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}{rn\_ebb} namelist parameter, setting \np{rn_ebb}{rn\_ebb}\forcode{ = 67.83} corresponds
+with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter, setting \np[=67.83]{rn_ebb}{rn\_ebb} corresponds
 to $\alpha_{CB} = 100$.
-Further setting  \np{ln_mxl0}{ln\_mxl0}\forcode{ =.true.},  applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,
+Further setting  \np[=.true.]{ln_mxl0}{ln\_mxl0},  applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,
 with $\beta$ hard coded to the Stacey's value.
 Note that a minimal threshold of \np{rn_emin0}{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the
@@ -385,5 +387,5 @@
 (\ie\ near-inertial oscillations and ocean swells and waves).
 
-When using this parameterization (\ie\ when \np{nn_etau}{nn\_etau}\forcode{=1}),
+When using this parameterization (\ie\ when \np[=1]{nn_etau}{nn\_etau}),
 the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations,
 swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition,
@@ -398,9 +400,9 @@
 (no penetration if $f_i=1$, \ie\ if the ocean is entirely covered by sea-ice).
 The value of $f_r$, usually a few percents, is specified through \np{rn_efr}{rn\_efr} namelist parameter.
-The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn_etau}{nn\_etau}\forcode{=0}) or
+The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np[=0]{nn_etau}{nn\_etau}) 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}{nn\_etau}\forcode{=1}).
-
-Note that two other option exist, \np{nn_etau}{nn\_etau}\forcode{=2, 3}.
+(\np[=1]{nn_etau}{nn\_etau}).
+
+Note that two other option exist, \np[=2, 3]{nn_etau}{nn\_etau}.
 They correspond to applying \autoref{eq: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 penetrates the ocean.
@@ -508,5 +510,5 @@
   \caption[Set of predefined GLS parameters or equivalently predefined turbulence models available]{
     Set of predefined GLS parameters, or equivalently predefined turbulence models available with
-    \protect\np{ln_zdfgls}{ln\_zdfgls}\forcode{=.true.} and controlled by
+    \protect\np[=.true.]{ln_zdfgls}{ln\_zdfgls} and controlled by
     the \protect\np{nn_clos}{nn\_clos} namelist variable in \protect\nam{zdf_gls}{zdf\_gls}.}
   \label{tab:ZDF_GLS}
@@ -519,5 +521,5 @@
 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88},
 or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01}
-(\np{nn_stab_func}{nn\_stab\_func}\forcode{=0, 3}, resp.).
+(\np[=0, 3]{nn_stab_func}{nn\_stab\_func}, resp.).
 The value of $C_{0\mu}$ depends on the choice of the stability function.
 
@@ -525,5 +527,5 @@
 Neumann condition through \np{nn_bc_surf}{nn\_bc\_surf} and \np{nn_bc_bot}{nn\_bc\_bot}, resp.
 As for TKE closure, the wave effect on the mixing is considered when
-\np{rn_crban}{rn\_crban}\forcode{ > 0.} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}.
+\np[ > 0.]{rn_crban}{rn\_crban} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}.
 The \np{rn_crban}{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and
 \np{rn_charn}{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}.
@@ -536,5 +538,5 @@
 the entrainment depth predicted in 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}{ln\_length\_lim}\forcode{=.true.},
+The clipping is only activated if \np[=.true.]{ln_length_lim}{ln\_length\_lim},
 and the $c_{lim}$ is set to the \np{rn_clim_galp}{rn\_clim\_galp} value.
 
@@ -707,5 +709,5 @@
 
 Options are defined through the \nam{zdf}{zdf} namelist variables.
-The non-penetrative convective adjustment is used when \np{ln_zdfnpc}{ln\_zdfnpc}\forcode{=.true.}.
+The non-penetrative convective adjustment is used when \np[=.true.]{ln_zdfnpc}{ln\_zdfnpc}.
 It is applied at each \np{nn_npc}{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
@@ -751,12 +753,12 @@
 
 Options are defined through the  \nam{zdf}{zdf} namelist variables.
-The enhanced vertical diffusion parameterisation is used when \np{ln_zdfevd}{ln\_zdfevd}\forcode{=.true.}.
+The enhanced vertical diffusion parameterisation is used when \np[=.true.]{ln_zdfevd}{ln\_zdfevd}.
 In this case, the vertical eddy mixing coefficients are assigned very large values
 in regions where the stratification is unstable
 (\ie\ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}.
-This is done either on tracers only (\np{nn_evdm}{nn\_evdm}\forcode{=0}) or
-on both momentum and tracers (\np{nn_evdm}{nn\_evdm}\forcode{=1}).
-
-In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np{nn_evdm}{nn\_evdm}\forcode{=1},
+This is done either on tracers only (\np[=0]{nn_evdm}{nn\_evdm}) or
+on both momentum and tracers (\np[=1]{nn_evdm}{nn\_evdm}).
+
+In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np[=1]{nn_evdm}{nn\_evdm},
 the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ values also, are set equal to
 the namelist parameter \np{rn_avevd}{rn\_avevd}.
@@ -795,5 +797,5 @@
 The OSMOSIS 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}{ln\_zdfevd}\forcode{=.false.} should be used with the OSMOSIS scheme.
+therefore \np[=.false.]{ln_zdfevd}{ln\_zdfevd} should be used with the OSMOSIS scheme.
 % gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
 
@@ -1002,10 +1004,10 @@
     c_b^T = - r
 \]
-When \np{ln_lin}{ln\_lin} \forcode{= .true.}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}.
-Setting \np{ln_OFF}{ln\_OFF} \forcode{= .true.} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.
+When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}.
+Setting \np[=.true.]{ln_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.
 
 These values are assigned in \mdl{zdfdrg}.
 Note that there is support for local enhancement of these values via an externally defined 2D mask array
-(\np{ln_boost}{ln\_boost}\forcode{=.true.}) given in the \ifile{bfr\_coef} input NetCDF file.
+(\np[=.true.]{ln_boost}{ln\_boost}) 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
@@ -1043,5 +1045,5 @@
 $C_D$= \np{rn_Cd0}{rn\_Cd0}, and $e_b$ =\np{rn_bfeb2}{rn\_bfeb2}.
 Note that for applications which consider tides explicitly, a low or even zero value of \np{rn_bfeb2}{rn\_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array
-(\np{ln_boost}{ln\_boost}\forcode{=.true.}).
+(\np[=.true.]{ln_boost}{ln\_boost}).
 This works in the same way as for the linear friction case with non-zero masked locations increased by
 $mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}.
@@ -1055,5 +1057,5 @@
 In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using
 a "law of the wall" scaling. This assumes that the model vertical resolution can capture the logarithmic layer which typically occur for layers thinner than 1 m or so.
-If  \np{ln_loglayer}{ln\_loglayer} \forcode{= .true.}, $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness):
+If  \np[=.true.]{ln_loglayer}{ln\_loglayer}, $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness):
 \[
   C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2
@@ -1070,5 +1072,5 @@
 
 \noindent The log-layer enhancement can also be applied to the top boundary friction if
-under ice-shelf cavities are activated (\np{ln_isfcav}{ln\_isfcav}\forcode{=.true.}).
+under ice-shelf cavities are activated (\np[=.true.]{ln_isfcav}{ln\_isfcav}).
 %In this case, the relevant namelist parameters are \np{rn_tfrz0}{rn\_tfrz0}, \np{rn_tfri2}{rn\_tfri2} and \np{rn_tfri2_max}{rn\_tfri2\_max}.
 
@@ -1076,8 +1078,8 @@
 %       Explicit bottom Friction
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np{ln_drgimp}{ln\_drgimp}\forcode{=.false.})}
+\subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np[=.false.]{ln_drgimp}{ln\_drgimp})}
 \label{subsec:ZDF_drg_stability}
 
-Setting \np{ln_drgimp}{ln\_drgimp} \forcode{= .false.} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads:
+Setting \np[=.false.]{ln_drgimp}{ln\_drgimp} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads:
 
 At the top (below an ice shelf cavity):
@@ -1137,5 +1139,5 @@
 %       Implicit Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np{ln_drgimp}{ln\_drgimp}\forcode{=.true.})}
+\subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np[=.true.]{ln_drgimp}{ln\_drgimp})}
 \label{subsec:ZDF_drg_imp}
 
@@ -1170,5 +1172,5 @@
 \label{subsec:ZDF_drg_ts}
 
-With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \np{ln_drgimp}{ln\_drgimp}\forcode{= .false.} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie\ \forcode{ln_bt_fw=.false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \np{ln_drgimp}{ln\_drgimp}\forcode{= .true.},  stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions.
+With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \np[=.false.]{ln_drgimp}{ln\_drgimp} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie\ \forcode{ln_bt_fw=.false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \np[=.true.]{ln_drgimp}{ln\_drgimp},  stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions.
 
 The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO\ is as follows:
@@ -1508,7 +1510,7 @@
 % ================================================================
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_cfgs.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_cfgs.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_cfgs.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -253,5 +255,5 @@
 
 The GYRE configuration is set like an analytical configuration.
-Through \np{ln_read_cfg}{ln\_read\_cfg}\forcode{ = .false.} in \nam{cfg}{cfg} namelist defined in
+Through \np[=.false.]{ln_read_cfg}{ln\_read\_cfg} in \nam{cfg}{cfg} namelist defined in
 the reference configuration \path{./cfgs/GYRE_PISCES/EXPREF/namelist_cfg}
 analytical definition of grid in GYRE is done in usrdef\_hrg, usrdef\_zgr routines.
@@ -271,5 +273,5 @@
 For example, keeping a same model size on each processor while increasing the number of processor used is very easy,
 even though the physical integrity of the solution can be compromised.
-Benchmark is activate via \np{ln_bench}{ln\_bench}\forcode{ = .true.} in \nam{usr_def}{usr\_def} in
+Benchmark is activate via \np[=.true.]{ln_bench}{ln\_bench} in \nam{usr_def}{usr\_def} in
 namelist \path{./cfgs/GYRE_PISCES/EXPREF/namelist_cfg}.
 
@@ -299,5 +301,5 @@
 In particular, the AMM uses $s$-coordinates in the vertical rather than $z$-coordinates and
 is forced with tidal lateral boundary conditions using a Flather boundary condition from the BDY module.
-Also specific to the AMM configuration is the use of the GLS turbulence scheme (\np{ln_zdfgls}{ln\_zdfgls} \forcode{= .true.}).
+Also specific to the AMM configuration is the use of the GLS turbulence scheme (\np[=.true.]{ln_zdfgls}{ln\_zdfgls}).
 
 In addition to the tidal boundary condition the model may also take open boundary conditions from
@@ -308,7 +310,7 @@
 Unlike ordinary river points the Baltic inputs also include salinity and temperature data.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -329,7 +331,7 @@
 It has not been implemented.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -416,7 +418,7 @@
 
 % ================================================================
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex	(revision 11582)
@@ -1173,7 +1173,7 @@
 Nevertheless it is currently not available in the iso-neutral case.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\input{../../global/printindex}}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -299,5 +301,5 @@
 The default value is 1, as recommended by \citet{Roullet2000?}
 
-\colorbox{red}{\np{rnu}{rnu}\forcode{=1} to be suppressed from namelist !}
+\colorbox{red}{\np[=1]{rnu}{rnu} to be suppressed from namelist !}
 
 %-------------------------------------------------------------
@@ -313,7 +315,7 @@
 In particular, this means that in filtered case, the matrix to be inverted has to be recomputed at each time-step.
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex	(revision 11581)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex	(revision 11582)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+
+\onlyinsubfile{\makeindex}
 
 \begin{document}
@@ -86,5 +88,5 @@
 where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin coefficient.
 $\gamma$ is initialized as \np{rn_atfp}{rn\_atfp} (namelist parameter).
-Its default value is \np{rn_atfp}{rn\_atfp}\forcode{ = 10.e-3} (see \autoref{sec:TD_mLF}),
+Its default value is \np[=10.e-3]{rn_atfp}{rn\_atfp} (see \autoref{sec:TD_mLF}),
 causing only a weak dissipation of high frequency motions (\citep{farge-coulombier_phd87}).
 The addition of a time filter degrades the accuracy of the calculation from second to first order.
@@ -172,8 +174,8 @@
 
 The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated
-at \textit{now} time step. This refers to as the explicit free surface case in the code (\np{ln_dynspg_exp}{ln\_dynspg\_exp}\forcode{=.true.}).
+at \textit{now} time step. This refers to as the explicit free surface case in the code (\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}).
 This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation
 of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface
-(\np{ln_dynspg_ts}{ln\_dynspg\_ts}\forcode{=.true.}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc
+(\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc
 time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of
 the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}).
@@ -297,5 +299,5 @@
 When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step
 is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting
-the namelist variable \np{nn_euler}{nn\_euler}\forcode{=0}. Other options to control the time integration of the model
+the namelist variable \np[=0]{nn_euler}{nn\_euler}. Other options to control the time integration of the model
 are defined through the  \nam{run}{run} namelist variables.
 %%%
@@ -386,7 +388,7 @@
 }
 
-\biblio
-
-\pindex
+\onlyinsubfile{\bibliography{../main/bibliography}}
+
+\onlyinsubfile{\printindex}
 
 \end{document}
