Index: NEMO/branches/2021/dev_r14312_MPI_Interface/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/cfgs/SHARED/namelist_ref	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/cfgs/SHARED/namelist_ref	(revision 14403)
@@ -644,5 +644,5 @@
    rn_test_box             = 108.0,  116.0, -66.0, -58.0
    ln_use_calving          = .false. ! Use calving data even when nn_test_icebergs > 0
-   rn_speed_limit          = 0.      ! CFL speed limit for a berg
+   rn_speed_limit          = 0.      ! CFL speed limit for a berg (safe value is 0.4, see #2581)
    !
    ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data)
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/.svnignore
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/.svnignore	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/.svnignore	(revision 14403)
@@ -9,4 +9,7 @@
 *.lo*
 *.out
+*.pdf
+*.pyg
+*.tdo
 *.toc
 *.xdv
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/main/settings.tex
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/main/settings.tex	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/main/settings.tex	(revision 14403)
@@ -5,12 +5,10 @@
 \def\spcup{\vspace*{2.15cm}}
 \def\hdg{NEMO ocean engine}
-%\def\shdg{}
-\def\spcdn{\vspace*{1cm}}
-\def\autwd{0.25\linewidth}
-\def\lnlg{270pt}
-\def\abswd{0.65\linewidth}
+%\def\shdg{} %% No subheading
+\def\spcdn{\vspace*{1.00cm}}
+\def\autwd{0.25\linewidth}\def\lnlg{270pt}\def\abswd{0.65\linewidth}
 
 %% Color in cmyk model for manual theme (frontpage banner, links and chapter boxes)
-\def\clr{1,.60,0,.4}
+\def\clr{1.0,0.6,0.0,0.4}
 
 %% IPSL publication number
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex	(revision 14403)
@@ -6,6 +6,4 @@
 \label{apdx:DOMCFG}
 
-%    {\em 4.0} & {\em Andrew Coward} & {\em Created at v4.0 from materials removed from chap\_DOM that are still relevant to the \forcode{DOMAINcfg} tool and which illustrate and explain the choices to be made by the user when setting up new domains }  \\
-
 \chaptertoc
 
@@ -14,11 +12,10 @@
 {\footnotesize
   \begin{tabularx}{\textwidth}{l||X|X}
-    Release & Author(s) & Modifications \\
-    \hline
-    {\em   next}& {\em Pierre Mathiot} & {\em add ice shelf and closed sea option description } \\
-    {\em   4.0} & {\em Andrew Coward}  & {\em Created at v4.0 from materials removed from chap\_DOM that are still relevant to the \forcode{DOMAINcfg} tool and which illustrate and explain the choices to be made by the user when setting up new domains }  \\
-    {\em   3.6} & {\em ...} & {\em ...} \\
-    {\em   3.4} & {\em ...} & {\em ...} \\
-    {\em <=3.4} & {\em ...} & {\em ...}
+    Release     & Author(s)            & Modifications                                                 \\
+    \hline
+    {\em  next} & {\em Pierre Mathiot} & {\em Add ice shelf and closed sea option description        } \\
+    {\em   4.0} & {\em  Andrew Coward} & {\em Creation from materials removed from \autoref{chap:DOM}
+                                              that are still relevant to the DOMAINcfg tool
+                                              when setting up new domains                            }
   \end{tabularx}
 }
@@ -592,5 +589,5 @@
   the \textit{isfdraft\_meter} file (Netcdf format). This file need to include the \textit{isf\_draft} variable. 
   A positive value will mean ice shelf/ocean or ice shelf bedrock interface below the reference 0m ssh. 
-  The exact shape of the ice shelf cavity (grounding line position and minimum thickness of the water column under an ice shelf, ...) can be specify in \nam{zgr_isf}{zgr_isf}.
+  The exact shape of the ice shelf cavity (grounding line position and minimum thickness of the water column under an ice shelf, ...) can be specify in \nam{zgr_isf}{zgr\_isf}.
 
 \begin{listing}
@@ -616,64 +613,64 @@
 \end{listing}
 
-   The options available to define the shape of the under ice shelf cavities are listed in \nam{zgr_isf}{zgr_isf} (\texttt{DOMAINcfg} only, \autoref{lst:namzgr_isf}).
-
-   \subsection{Model ice shelf draft definition}
-   \label{subsec:zgrisf_isfd}
-
-   First of all, the tool make sure, the ice shelf draft ($h_{isf}$) is sensible and compatible with the bathymetry.
-   There are 3 compulsory steps to achieve this:
-
-   \begin{description}
-   \item{\np{rn_isfdep_min}{rn\_isfdep\_min}:} this is the minimum ice shelf draft. This is to make sure there is no ridiculous thin ice shelf. If \np{rn_isfdep_min}{rn\_isfdep\_min} is smaller than the surface level, \np{rn_isfdep_min}{rn\_isfdep\_min} is set to $e3t\_1d(1)$. 
-   Where $h_{isf} < MAX(e3t\_1d(1),\np{rn_isfdep_min}{rn\_isfdep\_min}$), $h_{isf}$ is set to \np{rn_isfdep_min}{rn\_isfdep\_min}.
-
-   \item{\np{rn_glhw_min}{rn\_glhw\_min}:} This parameter is used to define the grounding line position.
-   Where the difference between the bathymetry and the ice shelf draft is smaller than \np{rn_glhw_min}{rn\_glhw\_min}, the cell are grounded (ie masked). 
-   This step is needed to take into account possible small mismatch between ice shelf draft value and bathymetry value (sources are coming from different grid, different data processes, rounding error, ...).
-
-   \item{\np{rn_isfhw_min}{rn\_isfhw\_min}:} This parameter is the minimum water column thickness in the cavity. 
-   Where the water column thickness is lower than \np{rn_isfhw_min}{rn\_isfhw\_min}, the ice shelf draft is adjusted to match this criterion. 
-   If for any reason, this adjustement break the minimum ice shelf draft allowed (\np{rn_isfdep_min}{rn\_isfdep\_min}), the cell is masked.
-   \end{description}
-
-   Once all these adjustements are made, if the water column thickness contains one cell wide channels, these channels can be closed using \np{ln_isfchannel}{ln\_isfchannel}.  
+   The options available to define the shape of the under ice shelf cavities are listed in \nam{zgr_isf}{zgr\_isf} (\texttt{DOMAINcfg} only, \autoref{lst:namzgr_isf}).
+
+\subsection{Model ice shelf draft definition}
+\label{subsec:zgrisf_isfd}
+
+First of all, the tool make sure, the ice shelf draft ($h_{isf}$) is sensible and compatible with the bathymetry.
+There are 3 compulsory steps to achieve this:
+
+\begin{description}
+\item{\np{rn_isfdep_min}{rn\_isfdep\_min}:} this is the minimum ice shelf draft. This is to make sure there is no ridiculous thin ice shelf. If \np{rn_isfdep_min}{rn\_isfdep\_min} is smaller than the surface level, \np{rn_isfdep_min}{rn\_isfdep\_min} is set to $e3t\_1d(1)$. 
+  Where $h_{isf} < MAX(e3t\_1d(1),rn\_isfdep\_min)$, $h_{isf}$ is set to \np{rn_isfdep_min}{rn\_isfdep\_min}.
+
+\item{\np{rn_glhw_min}{rn\_glhw\_min}:} This parameter is used to define the grounding line position.
+  Where the difference between the bathymetry and the ice shelf draft is smaller than \np{rn_glhw_min}{rn\_glhw\_min}, the cell are grounded (ie masked). 
+  This step is needed to take into account possible small mismatch between ice shelf draft value and bathymetry value (sources are coming from different grid, different data processes, rounding error, ...).
+
+\item{\np{rn_isfhw_min}{rn\_isfhw\_min}:} This parameter is the minimum water column thickness in the cavity. 
+  Where the water column thickness is lower than \np{rn_isfhw_min}{rn\_isfhw\_min}, the ice shelf draft is adjusted to match this criterion. 
+  If for any reason, this adjustement break the minimum ice shelf draft allowed (\np{rn_isfdep_min}{rn\_isfdep\_min}), the cell is masked.
+\end{description}
+
+Once all these adjustements are made, if the water column thickness contains one cell wide channels, these channels can be closed using \np{ln_isfchannel}{ln\_isfchannel}.  
  
-   \subsection{Model top level definition}
-   After the definition of the ice shelf draft, the tool defines the top level. 
-   The compulsory criterion is that the water column needs at least 2 wet cells in the water column at U- and V-points.
-   To do so, if there one cell wide water column, the tools adjust the ice shelf draft to fillful the requierement.\\
-
-   The process is the following:
-   \begin{description}
-   \item{step 1:} The top level is defined in the same way as the bottom level is defined.
-   \item{step 2:} The isolated grid point in the bathymetry are filled (as it is done in a domain without ice shelf)
-   \item{step 3:} The tools make sure, the top level is above or equal to the bottom level
-   \item{step 4:} If the water column at a U- or V- point is one wet cell wide, the ice shelf draft is adjusted. So the actual top cell become fully open and the new
-   top cell thickness is set to the minimum cell thickness allowed (following the same logic as for the bottom partial cell). This step is iterated 4 times to ensure the condition is fullfill along the 4 sides of the cell.
-   \end{description}
-
-   In case of steep slope and shallow water column, it likely that 2 cells are disconnected (bathymetry above its neigbourging ice shelf draft). 
-   The option \np{ln_isfconnect}{ln\_isfconnect} allow the tool to force the connection between these 2 cells.
-   Some limiters in meter or levels on the digging allowed by the tool are available (respectively, \np{rn_zisfmax}{rn\_zisfmax} or \np{rn_kisfmax}{rn\_kisfmax}).
-   This will prevent the formation of subglacial lakes at the expense of long vertical pipe to connect cells at very different levels.
-
-   \subsection{Subglacial lakes}
-   Despite careful setting of your ice shelf draft and bathymetry input file as well as setting described in \autoref{subsec:zgrisf_isfd}, some situation are unavoidable.
-   For exemple if you setup your ice shelf draft and bathymetry to do ocean/ice sheet coupling, 
-   you may decide to fill the whole antarctic with a bathymetry and an ice shelf draft value (ice/bedrock interface depth when grounded). 
-   If you do so, the subglacial lakes will show up (Vostock for example). An other possibility is with coarse vertical resolution, some ice shelves could be cut in 2 parts: 
-   one connected to the main ocean and an other one closed which can be considered as a subglacial sea be the model.\\
-
-   The namelist option \np{ln_isfsubgl}{ln\_isfsubgl} allow you to remove theses subglacial lakes.
-   This may be useful for esthetical reason or for stability reasons:
-
-   \begin{description}
-   \item $\bullet$ In a subglacial lakes, in case of very weak circulation (often the case), the only heat flux is the conductive heat flux through the ice sheet. 
-         This will lead to constant freezing until water reaches -20C. 
-         This is one of the defitiency of the 3 equation melt formulation (for details on this formulation, see: \autoref{sec:isf}).
-   \item $\bullet$ In case of coupling with an ice sheet model, 
-         the ssh in the subglacial lakes and the main ocean could be very different (ssh initial adjustement for example), 
-         and so if for any reason both a connected at some point, the model is likely to fall over.\\
-   \end{description}
+\subsection{Model top level definition}
+After the definition of the ice shelf draft, the tool defines the top level. 
+The compulsory criterion is that the water column needs at least 2 wet cells in the water column at U- and V-points.
+To do so, if there one cell wide water column, the tools adjust the ice shelf draft to fillful the requierement.\\
+
+The process is the following:
+\begin{description}
+\item{step 1:} The top level is defined in the same way as the bottom level is defined.
+\item{step 2:} The isolated grid point in the bathymetry are filled (as it is done in a domain without ice shelf)
+\item{step 3:} The tools make sure, the top level is above or equal to the bottom level
+\item{step 4:} If the water column at a U- or V- point is one wet cell wide, the ice shelf draft is adjusted. So the actual top cell become fully open and the new
+  top cell thickness is set to the minimum cell thickness allowed (following the same logic as for the bottom partial cell). This step is iterated 4 times to ensure the condition is fullfill along the 4 sides of the cell.
+\end{description}
+
+In case of steep slope and shallow water column, it likely that 2 cells are disconnected (bathymetry above its neigbourging ice shelf draft). 
+The option \np{ln_isfconnect}{ln\_isfconnect} allow the tool to force the connection between these 2 cells.
+Some limiters in meter or levels on the digging allowed by the tool are available (respectively, \np{rn_zisfmax}{rn\_zisfmax} or \np{rn_kisfmax}{rn\_kisfmax}).
+This will prevent the formation of subglacial lakes at the expense of long vertical pipe to connect cells at very different levels.
+
+\subsection{Subglacial lakes}
+Despite careful setting of your ice shelf draft and bathymetry input file as well as setting described in \autoref{subsec:zgrisf_isfd}, some situation are unavoidable.
+For exemple if you setup your ice shelf draft and bathymetry to do ocean/ice sheet coupling, 
+you may decide to fill the whole antarctic with a bathymetry and an ice shelf draft value (ice/bedrock interface depth when grounded). 
+If you do so, the subglacial lakes will show up (Vostock for example). An other possibility is with coarse vertical resolution, some ice shelves could be cut in 2 parts: 
+one connected to the main ocean and an other one closed which can be considered as a subglacial sea be the model.\\
+
+The namelist option \np{ln_isfsubgl}{ln\_isfsubgl} allow you to remove theses subglacial lakes.
+This may be useful for esthetical reason or for stability reasons:
+
+\begin{description}
+\item $\bullet$ In a subglacial lakes, in case of very weak circulation (often the case), the only heat flux is the conductive heat flux through the ice sheet. 
+  This will lead to constant freezing until water reaches -20C. 
+  This is one of the defitiency of the 3 equation melt formulation (for details on this formulation, see: \autoref{sec:isf}).
+\item $\bullet$ In case of coupling with an ice sheet model, 
+  the ssh in the subglacial lakes and the main ocean could be very different (ssh initial adjustement for example), 
+  and so if for any reason both a connected at some point, the model is likely to fall over.\\
+\end{description}
 
 \section{Closed sea definition}
@@ -707,5 +704,5 @@
 \end{listing}
 
-The options available to define the closed seas and how closed sea net fresh water input will be redistributed by NEMO are listed in \nam{clo}{dom_clo} (\texttt{DOMAINcfg} only).
+The options available to define the closed seas and how closed sea net fresh water input will be redistributed by NEMO are listed in \nam{dom_clo}{dom\_clo} (\texttt{DOMAINcfg} only).
 The individual definition of each closed sea is managed by \np{sn_lake}{sn\_lake}. In this fields the user needs to define:\\
    \begin{description}
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 14403)
@@ -1179,5 +1179,5 @@
 %% =================================================================================================
 \section[Ice Shelf (ISF)]{Interaction with ice shelves (ISF)}
-\label{sec:isf}
+\label{sec:SBC_isf}
 
 \begin{listing}
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/chap_TRA.tex
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 14403)
@@ -733,5 +733,5 @@
   (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
 \item [\textit{fwfisf}] The mass flux associated with ice shelf melt,
-  (see \autoref{sec:isf} for further details on how the ice shelf melt is computed and applied).
+  (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied).
 \end{labeling}
 
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/SI3/main/settings.tex
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/SI3/main/settings.tex	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/SI3/main/settings.tex	(revision 14403)
@@ -1,20 +1,18 @@
-%% Engine (folder name)
-\def \engine{SI3}
+%% Engine
+\def\eng{SI3}
 
-%% Title and cover page settings
-\def \spacetop{  \vspace*{1.2cm} }
-\def \heading{Sea Ice modelling Integrated Initiative (SI$^3$)}
-\def \subheading{The NEMO sea ice engine}
-\def \spacedown{ \vspace*{  1cm} }
-\def \authorswidth{0.2\linewidth}
-\def \rulelenght{230pt}
-\def \abstractwidth{0.65\linewidth}
+%% Cover page
+\def\spcup{\vspace*{1.20cm}}
+\def \hdg{Sea Ice modelling Integrated Initiative (SI$^3$)}
+\def\shdg{The NEMO sea ice engine                         }
+\def\spcdn{\vspace*{1.00cm}}
+\def\autwd{0.20\linewidth}\def\lnlg{230pt}\def\abswd{0.65\linewidth}
 
-%% Color for document (frontpage banner, links and chapter boxes)
-\def \setmanualcolor{ \definecolor{manualcolor}{cmyk}{0, 0, 0, 0.4} }
+%% Color in cmyk model for manual theme (frontpage banner, links and chapter boxes)
+\def\clr{0.0,0.0,0.0,0.4}
 
 %% IPSL publication number
-\def \ipslnum{31}
+\def\ipsl{31}
 
-%% Zenodo ID, i.e. doi:10.5281/zenodo.\([0-9]*\)
-\def \zid{1471689}
+%% Zenodo ID, i.e. doi:10.5281/zenodo.\zid
+\def\zid{1471689}
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/main/bibliography.bib
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/main/bibliography.bib	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/main/bibliography.bib	(revision 14403)
@@ -185,4 +185,16 @@
   journal	= {Global Biogeochemical Cycles},
   publisher	= {American Geophysical Union (AGU)}
+}
+
+@techreport{      gibson_trpt86,
+  title         = "Standards for software development and maintenance",
+  pages         = "21",
+  series        = "ECMWF Technical Memoranda",
+  number        = "120",
+  author        = "J. K. Gibson",
+  institution   = "ECMWF Operations Department; Reading, United Kingdom",
+  year          = "1986",
+  month         = "aug",
+  doi           = "10.21957/gi113q4gn"
 }
 
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/main/settings.tex
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/main/settings.tex	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/main/settings.tex	(revision 14403)
@@ -1,20 +1,18 @@
-%% Engine (folder name)
-\def \engine{TOP}
+%% Engine
+\def\eng{TOP}
 
-%% Title and cover page settings
-\def \spacetop{  \vspace*{1.3cm} }
-\def \heading{Tracers in Ocean Paradigm (TOP)}
-\def \subheading{The NEMO passive tracers engine}
-\def \spacedown{ \vspace*{  1cm} }
-\def \authorswidth{0.15\linewidth}
-\def \rulelenght{110pt}
-\def \abstractwidth{0.7\linewidth}
+%% Cover page
+\def\spcup{\vspace*{1.30cm}}
+\def \hdg{Tracers in Ocean Paradigm (TOP)}
+\def\shdg{The NEMO passive tracers engine}
+\def\spcdn{\vspace*{1.00cm}}
+\def\autwd{0.15\linewidth}\def\lnlg{110pt}\def\abswd{0.70\linewidth}
 
-%% Color for document (frontpage banner, links and chapter boxes)
-\def \setmanualcolor{ \definecolor{manualcolor}{cmyk}{1, 0, 1, .4} }
+%% Color in cmyk model for manual theme (frontpage banner, links and chapter boxes)
+\def\clr{1.0,0.0,1.0,0.4}
 
 %% IPSL publication number
-\def \ipslnum{28}
+\def\ipsl{28}
 
-%% Zenodo ID, i.e. doi:10.5281/zenodo.\([0-9]*\)
-\def \zid{1471700}
+%% Zenodo ID, i.e. doi:10.5281/zenodo.\zid
+\def\zid{1471700}
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/subfiles/model_description.tex
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/subfiles/model_description.tex	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/TOP/subfiles/model_description.tex	(revision 14403)
@@ -9,6 +9,4 @@
 \newcommand{\Dcq}{\Delta ^{14}\mathrm{C}}
 \newcommand{\Rq}{\mathrm{^{14}{R}}}
-\newcommand{\CODE}[1]{\textsc{#1}}
-%\newcommand{\CODE}[1]{\textcolor{black}{\textsc{#1}}\xspace}
 
 \chapter{Model Description}
@@ -28,7 +26,11 @@
 where expressions of $D^{lC}$ and $D^{vC}$ depend on the choice for the lateral and vertical subgrid scale parameterizations, see equations 5.10 and 5.11 in \citep{nemo_manual}
 
-{S(C)} , the first term on the right hand side of \ref{Eq_tracer}; is the SMS - Source Minus Sink - inherent to the tracer.  In the case of biological tracer such as phytoplankton, {S(C)} is the balance between phytoplankton growth and its decay through mortality and grazing. In the case of a tracer comprising carbon,  {S(C)} accounts for gas exchange, river discharge, flux to the sediments, gravitational sinking and other biological processes. In the case of a radioactive tracer, {S(C)} is simply loss due to radioactive decay.
-
-The second term (within brackets) represents the advection of the tracer in the three directions. It can be interpreted as the budget between the incoming and outgoing tracer fluxes in a volume $T$-cells $b_t= e_{1t}\,e_{2t}\,e_{3t}$
+{S(C)} , the first term on the right hand side of \autoref{Eq_tracer}; is the SMS - Source Minus Sink - inherent to the tracer.
+In the case of biological tracer such as phytoplankton, {S(C)} is the balance between phytoplankton growth and its decay through mortality and grazing.
+In the case of a tracer comprising carbon,  {S(C)} accounts for gas exchange, river discharge, flux to the sediments, gravitational sinking and other biological processes.
+In the case of a radioactive tracer, {S(C)} is simply loss due to radioactive decay.
+
+The second term (within brackets) represents the advection of the tracer in the three directions.
+It can be interpreted as the budget between the incoming and outgoing tracer fluxes in a volume $T$-cells $b_t= e_{1t}\,e_{2t}\,e_{3t}$
 
 The third term  represents the change due to lateral diffusion.
@@ -46,5 +48,6 @@
 \label{sec:TopInt}
 
-TOP is the NEMO hardwired interface toward biogeochemical models and provide the physical constraints/boundaries for oceanic tracers. It consists of a modular framework to handle multiple ocean tracers, including also a variety of built-in modules.
+TOP is the NEMO hardwired interface toward biogeochemical models and provide the physical constraints/boundaries for oceanic tracers.
+It consists of a modular framework to handle multiple ocean tracers, including also a variety of built-in modules.
 
 This component of the NEMO framework allows one to exploit available modules  and further develop a range of applications, spanning from the implementation of a dye passive tracer to evaluate dispersion processes (by means of MY\_TRC), track water masses age (AGE module), assess the ocean interior penetration of persistent chemical compounds (e.g., gases like CFC or even PCBs), up to the full set of equations involving marine biogeochemical cycles.
@@ -61,5 +64,6 @@
         \item \textbf{AGE}	    :    Water age tracking
         \item \textbf{MY\_TRC}  :	Template for creation of new modules and external BGC models coupling
-        \item \textbf{PISCES}    :	 Built in BGC model. See \citep{aumont_2015} for a throughout description.
+        \item \textbf{PISCES}    :	 Built in BGC model.
+See \citep{aumont_2015} for a throughout description.
 \end{itemize}
 %  ----------------------------------------------------------
@@ -69,17 +73,21 @@
 The passive tracer transport component  shares the same advection/diffusion routines with the dynamics, with specific treatment of some features like the surface boundary conditions, or the positivity of passive tracers concentrations.
 
- \subsection{ Advection}
+\subsection{Advection}
 %------------------------------------------namtrc_adv----------------------------------------------------
 \nlst{namtrc_adv}
 %-------------------------------------------------------------------------------------------------------------
-The advection schemes used for the passive tracers are the same than the ones for $T$ and $S$ and described in section 5.1 of \citep{nemo_manual}. The choice of an advection scheme  can be selected independently and  can differ from the ones used for active tracers. This choice is made in the \textit{namtrc\_adv} namelist, by  setting to \textit{true} one and only one of the logicals \textit{ln\_trcadv\_xxx}, the same way of what is done for dynamics.
-cen2, MUSCL2, and UBS are not \textit{positive} schemes meaning that negative values can appear in an initially strictly positive tracer field which is advected, implying that false extrema are permitted. Their use is not recommended on passive tracers
-
- \subsection{ Lateral diffusion}
+The advection schemes used for the passive tracers are the same than the ones for $T$ and $S$ and described in section 5.1 of \citep{nemo_manual}.
+The choice of an advection scheme  can be selected independently and  can differ from the ones used for active tracers.
+This choice is made in the \textit{namtrc\_adv} namelist, by  setting to \textit{true} one and only one of the logicals \textit{ln\_trcadv\_xxx}, the same way of what is done for dynamics.
+cen2, MUSCL2, and UBS are not \textit{positive} schemes meaning that negative values can appear in an initially strictly positive tracer field which is advected, implying that false extrema are permitted.
+Their use is not recommended on passive tracers
+
+\subsection{Lateral diffusion}
 %------------------------------------------namtrc_ldf----------------------------------------------------
 \nlst{namtrc_ldf}
 %-------------------------------------------------------------------------------------------------------------
-In NEMO v4.0, the passive tracer diffusion has necessarily the same form as the active tracer diffusion, meaning that the numerical scheme must be the same. However the passive tracer mixing coefficient can be chosen as a multiple of the active ones by changing the value of \textit{rn\_ldf\_multi} in namelist \textit{namtrc\_ldf}. The choice of numerical scheme is then set  in the \nam{namtra_ldf}{namtra\_ldf} namelist for the dynamic described in section 5.2 of \citep{nemo_manual}.
-
+In NEMO v4.0, the passive tracer diffusion has necessarily the same form as the active tracer diffusion, meaning that the numerical scheme must be the same.
+However the passive tracer mixing coefficient can be chosen as a multiple of the active ones by changing the value of \textit{rn\_ldf\_multi} in namelist \textit{namtrc\_ldf}.
+The choice of numerical scheme is then set in the \forcode{&namtra_ldf} namelist for the dynamic described in section 5.2 of \citep{nemo_manual}.
 
 %-----------------We also offers the possibility to increase zonal equatorial diffusion for passive tracers by introducing an enhanced zonal diffusivity coefficent in the equatorial domain which can be defined by the equation below :
@@ -88,5 +96,5 @@
 %-----------------\end{equation}
 
- \subsection{ Tracer damping}
+\subsection{Tracer damping}
 
 %------------------------------------------namtrc_dmp----------------------------------------------------
@@ -94,8 +102,12 @@
 %-------------------------------------------------------------------------------------------------------------
 
-The use of newtonian damping  to climatological fields or observations is also coded, sharing the same routine dans active tracers. Boolean variables are defined in the namelist\_top\_ref to select the tracers on which restoring is applied
-Options are defined through the \nam{namtrc_dmp}{namtrc\_dmp} namelist variables. The restoring term is added when the namelist parameter \np{ln\_trcdmp} is set to true. The restoring coefficient is a three-dimensional array read in a file, which name is specified by the namelist variable \np{cn\_resto\_tr}. This netcdf file can be generated using the DMP\_TOOLS tool.
-
- \subsection{ Tracer positivity}
+The use of newtonian damping  to climatological fields or observations is also coded, sharing the same routine dans active tracers.
+Boolean variables are defined in the namelist\_top\_ref to select the tracers on which restoring is applied
+Options are defined through the \nam{trc_dmp}{trc\_dmp} namelist variables.
+The restoring term is added when the namelist parameter \np{ln\_trcdmp} is set to true.
+The restoring coefficient is a three-dimensional array read in a file, which name is specified by the namelist variable \np{cn\_resto\_tr}.
+This netcdf file can be generated using the DMP\_TOOLS tool.
+
+\subsection{Tracer positivity}
 
 %------------------------------------------namtrc_rad----------------------------------------------------
@@ -103,6 +115,8 @@
 %-------------------------------------------------------------------------------------------------------------
 
-Sometimes, numerical scheme can generates negative values of passive tracers concentration that must be positive. For exemple,  isopycnal diffusion can created extrema. The trcrad routine artificially corrects negative concentrations with a very crude solution that either sets negative concentration to zero without adjusting the tracer budget, or by removing negative concentration and keeping mass conservation.
-The treatment of negative concentrations is an option and can be selected in the namelist \nam{namtrc_rad}{namtrc\_rad} by setting the parameter \np{ln\_trcrad}  to true.
+Sometimes, numerical scheme can generates negative values of passive tracers concentration that must be positive.
+For exemple,  isopycnal diffusion can created extrema.
+The trcrad routine artificially corrects negative concentrations with a very crude solution that either sets negative concentration to zero without adjusting the tracer budget, or by removing negative concentration and keeping mass conservation.
+The treatment of negative concentrations is an option and can be selected in the namelist \nam{trc_rad}{trc\_rad} by setting the parameter \np{ln\_trcrad}  to true.
 
 \section{The SMS modules}
@@ -119,31 +133,42 @@
 %----------------------------------------------------------------------------------------------------------
 
-
- An `ideal age' tracer is integrated online in TOP when \textit{ln\_age} = \texttt{.true.} in namelist \textit{namtrc}. This tracer marks the length of time in units of years that fluid has spent in the interior of the ocean, insulated from exposure to the atmosphere. Thus, away from the surface for $z<-H_{\mathrm{Age}}$ where $H_{\mathrm{Age}}$ is specified by the \textit{namage} namelist variable \textit{rn\_age\_depth}, whose default value is 10~m, there is a source $\mathrm{SMS_{\mathrm{Age}}}$ of the age tracer $A$:
+An `ideal age' tracer is integrated online in TOP when \textit{ln\_age} = \texttt{.true.} in namelist \textit{namtrc}.
+This tracer marks the length of time in units of years that fluid has spent in the interior of the ocean, insulated from exposure to the atmosphere.
+Thus, away from the surface for $z<-H_{\mathrm{Age}}$ where $H_{\mathrm{Age}}$ is specified by the \textit{namage} namelist variable \textit{rn\_age\_depth}, whose default value is 10~m, there is a source $\mathrm{SMS_{\mathrm{Age}}}$ of the age tracer $A$:
+
 \begin{equation}
   \label{eq:TOP-age-interior}
   \mathrm{SMS_{\mathrm{Age}}} = 1 \mathrm{yr}\;^{-1} = 1/T_{\mathrm{year}},
- \end{equation}
- where the length of the current year $T_{\mathrm{year}} = 86400*N_{\mathrm{days\;in\;current\; year}}\;\mathrm{s}$, where $N_{\mathrm{days\;in\;current\; year}}$ may be 366 or 365 depending on whether the current year is a leap year or not.
- Near the surface, for $z>-H_{\mathrm{Age}}$, ideal age is relaxed back to zero:
+\end{equation}
+
+where the length of the current year $T_{\mathrm{year}} = 86400*N_{\mathrm{days\;in\;current\; year}}\;\mathrm{s}$, where $N_{\mathrm{days\;in\;current\; year}}$ may be 366 or 365 depending on whether the current year is a leap year or not.
+Near the surface, for $z>-H_{\mathrm{Age}}$, ideal age is relaxed back to zero:
+
 \begin{equation}
   \label{eq:TOP-age-surface}
    \mathrm{SMS_{\mathrm{Age}}} = -\lambda_{\mathrm{Age}}A,
- \end{equation}
- where the relaxation rate $\lambda_{\mathrm{Age}}$  (units $\mathrm{s}\;^{-1}$) is specified by the \textit{namage} namelist variable \textit{rn\_age\_kill\_rate} and has a default value of 1/7200~s. Since this relaxation is applied explicitly, this relaxation rate in principle should not exceed $1/\Delta t$, where $\Delta t$ is the time step used to step forward passive tracers (2 * \textit{nn\_dttrc * rn\_rdt} when the default  leapfrog time-stepping scheme is employed).
-
- Currently the 1-dimensional reference depth of the grid boxes is used rather than the dynamically evolving depth to determine whether the age tracer is incremented or relaxed to zero. This means that the tracer only works correctly in z-coordinates. To ensure that the forcing is independent of the level thicknesses, where the tracer cell at level $k$ has its upper face $z=-depw(k)$ above the depth $-H_{\mathrm{Age}}$, but its lower face $z=-depw(k+1)$ below that depth, then the age source
- \begin{equation}
+\end{equation}
+
+where the relaxation rate $\lambda_{\mathrm{Age}}$  (units $\mathrm{s}\;^{-1}$) is specified by the \textit{namage} namelist variable \textit{rn\_age\_kill\_rate} and has a default value of 1/7200~s.
+Since this relaxation is applied explicitly, this relaxation rate in principle should not exceed $1/\Delta t$, where $\Delta t$ is the time step used to step forward passive tracers (2 * \textit{nn\_dttrc * rn\_rdt} when the default  leapfrog time-stepping scheme is employed).
+
+Currently the 1-dimensional reference depth of the grid boxes is used rather than the dynamically evolving depth to determine whether the age tracer is incremented or relaxed to zero.
+This means that the tracer only works correctly in z-coordinates.
+To ensure that the forcing is independent of the level thicknesses, where the tracer cell at level $k$ has its upper face $z=-depw(k)$ above the depth $-H_{\mathrm{Age}}$, but its lower face $z=-depw(k+1)$ below that depth, then the age source
+
+\begin{equation}
   \label{eq:TOP-age-mixed}
    \mathrm{SMS_{\mathrm{Age}}} = -f_{\mathrm{kill}}\lambda_{\mathrm{Age}}A +f_{\mathrm{add}}/T_{\mathrm{year}} ,
- \end{equation}
- where
- \begin{align}
+\end{equation}
+
+where
+
+\begin{align}
     f_{\mathrm{kill}} &= e3t_k^{-1}(H_{\mathrm{Age}} - depw(k)) , \\
     f_{\mathrm{add}} &= 1 - f_{\mathrm{kill}}.
- \end{align}
-
-
- This implementation was first used in the CORE-II intercomparison runs described e.g.\ in \citet{danabasoglu_2014}.
+\end{align}
+
+
+This implementation was first used in the CORE-II intercomparison runs described e.g.\ in \citet{danabasoglu_2014}.
 
 \subsection{Inert carbons tracer}
@@ -155,10 +180,14 @@
 
 Chlorofluorocarbons 11 and 12 (CFC-11 and CFC-12), and sulfur hexafluoride (SF6), are synthetic chemicals manufactured for industrial and domestic applications from the early 20th century onwards.
-CFC-11 (CCl$_{3}$F) is a volatile liquid at room temperature, and was widely used in refrigeration. CFC-12 (CCl$_{2}$F$_{2}$) is a gas at room temperature, and, like CFC-11, was widely used as a refrigerant,
-and additionally as an aerosol propellant. SF6 (SF$_{6}$) is also a gas at room temperature, with a range of applications based around its property as an excellent electrical insulator (often replacing more toxic alternatives).
-All three are relatively inert chemicals that are both non-toxic and non-flammable, and their wide use has led to their accumulation within the Earth's atmosphere. Large-scale production of CFC-11 and CFC-12 began in the 1930s, while production of SF6 began in the 1950s, and their atmospheric concentration time-histories are shown in Figure \ref{img_cfcatm}.
+CFC-11 (CCl$_{3}$F) is a volatile liquid at room temperature, and was widely used in refrigeration.
+CFC-12 (CCl$_{2}$F$_{2}$) is a gas at room temperature, and, like CFC-11, was widely used as a refrigerant,
+and additionally as an aerosol propellant.
+SF6 (SF$_{6}$) is also a gas at room temperature, with a range of applications based around its property as an excellent electrical insulator (often replacing more toxic alternatives).
+All three are relatively inert chemicals that are both non-toxic and non-flammable, and their wide use has led to their accumulation within the Earth's atmosphere.
+Large-scale production of CFC-11 and CFC-12 began in the 1930s, while production of SF6 began in the 1950s, and their atmospheric concentration time-histories are shown in Figure \autoref{img_cfcatm}.
 As can be seen in the figure, while the concentration of SF6 continues to rise to the present  day, the concentrations of both CFC-11 and CFC-12 have levelled off and declined since around the 1990s.
 These declines have been driven by the Montreal Protocol (effective since August 1989), which has banned the production of CFC-11 and CFC-12 (as well as other CFCs) because of their role in the depletion of
-stratospheric ozone (O$_{3}$), critical in decreasing the flux of ultraviolet radiation to the Earth's surface. Separate to this role in ozone-depletion, all three chemicals are significantly more potent greenhouse gases
+stratospheric ozone (O$_{3}$), critical in decreasing the flux of ultraviolet radiation to the Earth's surface.
+Separate to this role in ozone-depletion, all three chemicals are significantly more potent greenhouse gases
 than CO$_{2}$ (especially SF6), although their relatively low atmospheric concentrations limit their role in climate change. \\
 
@@ -171,5 +200,5 @@
 % This release began in the 1930s for CFC-11 and CFC-12, and the 1950s for SF6, and
 % regularly increasing their atmospheric concentration until the 1090s, 2000s for respectively CFC11, CFC12,
-% and is still increasing, and SF6 (see Figure \ref{img_cfcatm}).  \\
+% and is still increasing, and SF6 (see Figure \autoref{img_cfcatm}).  \\
 
 The ocean is a notable sink for all three gases, and their relatively recent occurrence in the atmosphere, coupled to the ease of making high precision measurements of their dissolved concentrations, has made them
@@ -177,8 +206,10 @@
 Because they only enter the ocean via surface air-sea exchange, and are almost completely chemically and biologically inert, their distribution within the ocean interior reveals its ventilation via transport and mixing.
 Measuring the dissolved concentrations of the gases -- as well as the mixing ratios between them -- shows circulation pathways within the ocean as well as water mass ages (i.e. the time since last contact with the
-atmosphere). This feature of the gases has made them valuable across a wide range of oceanographic problems. One use lies in ocean modelling, where they can be used to evaluate the realism of the circulation and
+atmosphere).
+This feature of the gases has made them valuable across a wide range of oceanographic problems.
+One use lies in ocean modelling, where they can be used to evaluate the realism of the circulation and
 ventilation of models, key for understanding the behaviour of wider modelled marine biogeochemistry (e.g. \citep{dutay_2002,palmieri_2015}). \\
 
-Modelling these gases (henceforth CFCs) in NEMO is done within the passive tracer transport module, TOP, using the conservation state equation \ref{Eq_tracer}
+Modelling these gases (henceforth CFCs) in NEMO is done within the passive tracer transport module, TOP, using the conservation state equation \autoref{Eq_tracer}
 
 Advection and diffusion of the CFCs in NEMO are calculated by the physical module, OPA,
@@ -198,6 +229,6 @@
 \end{eqnarray}
 
-Where $K_{w}$ is the piston velocity (in m~s$^{-1}$), as defined in Equation \ref{equ_Kw};
-$C_{sat}$ is the saturation concentration of the CFC tracer, as defined in Equation \ref{equ_C_sat};
+Where $K_{w}$ is the piston velocity (in m~s$^{-1}$), as defined in Equation \autoref{equ_Kw};
+$C_{sat}$ is the saturation concentration of the CFC tracer, as defined in Equation \autoref{equ_C_sat};
 $C_{surf}$ is the local surface concentration of the CFC tracer within the model (in mol~m$^{-3}$);
 and $f_{i}$ is the fractional sea-ice cover of the local ocean (ranging between 0.0 for ice-free ocean,
@@ -211,5 +242,5 @@
 \end{eqnarray}
 
-Where $Sol$ is the gas solubility in mol~m$^{-3}$~pptv$^{-1}$, as defined in Equation \ref{equ_Sol_CFC};
+Where $Sol$ is the gas solubility in mol~m$^{-3}$~pptv$^{-1}$, as defined in Equation \autoref{equ_Sol_CFC};
 and $P_{cfc}$ is the atmosphere concentration of the CFC (in parts per trillion by volume, pptv).
 This latter concentration is provided to the model by the historical time-series of \citet{bullister_2017}.
@@ -231,5 +262,5 @@
 $a$ is a constant re-estimated by \citet{wanninkhof_2014} to 0.251 (in $\frac{cm~h^{-1}}{(m~s^{-1})^{2}}$);
 and $u$ is the 10~m wind speed in m~s$^{-1}$ from either an atmosphere model or reanalysis atmospheric forcing.
-$Sc$ is the Schmidt number, and is calculated as follow, using coefficients from \citet{wanninkhof_2014} (see Table \ref{tab_Sc}).
+$Sc$ is the Schmidt number, and is calculated as follow, using coefficients from \citet{wanninkhof_2014} (see Table \autoref{tab_Sc}).
 
 \begin{eqnarray}
@@ -238,5 +269,5 @@
 \end{eqnarray}
 
-The solubility, $Sol$, used in Equation \ref{equ_C_sat} is calculated in mol~l$^{-1}$~atm$^{-1}$,
+The solubility, $Sol$, used in Equation \autoref{equ_C_sat} is calculated in mol~l$^{-1}$~atm$^{-1}$,
 and is specific for each gas.
 It has been experimentally estimated by \citet{warner_1985} as a function of temperature
@@ -260,13 +291,13 @@
 
 Where $T_{X}$ is $\frac{T + 273.16}{100}$, a function of temperature;
-and the $a_{x}$ and $b_{x}$ coefficients are specific for each gas (see Table \ref{tab_ref_CFC}).
+and the $a_{x}$ and $b_{x}$ coefficients are specific for each gas (see Table \autoref{tab_ref_CFC}).
 This is then converted to mol~m$^{-3}$~pptv$^{-1}$ assuming a constant atmospheric surface pressure of 1~atm.
 The solubility of CFCs thus decreases with rising $T$ while being relatively insensitive to salinity changes.
-Consequently, this translates to a pattern of solubility where it is greatest in cold, polar regions (see Figure \ref{img_cfcsol}).
+Consequently, this translates to a pattern of solubility where it is greatest in cold, polar regions (see Figure \autoref{img_cfcsol}).
 
 % AXY: not 100% sure about the units below; they might be in nanomol, or in seconds or years
 
 The standard outputs of the CFC module are seawater CFC concentrations (in mol~m$^{-3}$), the net air-sea flux (in mol~m$^{-2}$~d$^{-1}$) and the cumulative net air-sea flux (in mol~m$^{-2}$).
-Using XIOS, it is possible to obtain outputs such as the vertical integral of CFC concentrations (in mol~m$^{-2}$; see Figure \ref{img_cfcinv}).
+Using XIOS, it is possible to obtain outputs such as the vertical integral of CFC concentrations (in mol~m$^{-2}$; see Figure \autoref{img_cfcinv}).
 This property, when divided by the surface CFC concentration, estimates the local penetration depth (in m) of the CFC.
 
@@ -285,5 +316,5 @@
 
 \begin{table}[!t]
-\caption{Coefficients for fit of the CFCs solubility (Eq. \ref{equ_Sol_CFC}).}
+\caption{Coefficients for fit of the CFCs solubility (Eq. \autoref{equ_Sol_CFC}).}
 \vskip4mm
 \centering
@@ -302,5 +333,5 @@
 
 \begin{table}[!t]
-\caption{Coefficients for fit of the CFCs Schmidt number (Eq. \ref{equ_Sc}). }
+\caption{Coefficients for fit of the CFCs Schmidt number (Eq. \autoref{equ_Sc}). }
 \vskip4mm
 \centering
@@ -353,19 +384,28 @@
 %----------------------------------------------------------------------------------------------------------
 
-The C14 package implemented in NEMO by Anne Mouchet models ocean $\Dcq$. It offers several possibilities: $\Dcq$ as a physical tracer of the ocean ventilation (natural $\cq$), assessment of bomb radiocarbon uptake, as well as transient studies of paleo-historical ocean radiocarbon distributions.
+The C14 package implemented in NEMO by Anne Mouchet models ocean $\Dcq$.
+It offers several possibilities: $\Dcq$ as a physical tracer of the ocean ventilation (natural $\cq$), assessment of bomb radiocarbon uptake, as well as transient studies of paleo-historical ocean radiocarbon distributions.
 
 \subsubsection{Method}
 
- Let  $\Rq$ represent the ratio of $\cq$ atoms to the total number of carbon atoms in the sample, i.e. $\cq/\mathrm{C}$. Then, radiocarbon anomalies are reported as
+Let  $\Rq$ represent the ratio of $\cq$ atoms to the total number of carbon atoms in the sample, i.e. $\cq/\mathrm{C}$.
+Then, radiocarbon anomalies are reported as
+
 \begin{equation}
 \Dcq = \left( \frac{\Rq}{\Rq_\mathrm{ref}} - 1 \right) 10^3, \label{eq:c14dcq}
 \end{equation}
-where $\Rq_{\textrm{ref}}$ is a reference ratio. For the purpose of ocean ventilation studies $\Rq_{\textrm{ref}}$ is set to one.
+
+where $\Rq_{\textrm{ref}}$ is a reference ratio.
+For the purpose of ocean ventilation studies $\Rq_{\textrm{ref}}$ is set to one.
 
 Here we adopt the approach of \cite{fiadeiro_1982} and \cite{toggweiler_1989a,toggweiler_1989b} in which  the ratio $\Rq$ is transported rather than the individual concentrations C and $\cq$.
-This approach calls for a strong assumption, i.e., that of a homogeneous and constant dissolved inorganic carbon (DIC) field \citep{toggweiler_1989a,mouchet_2013}. While in terms of
-oceanic $\Dcq$, it yields similar results to approaches involving carbonate chemistry, it underestimates the bomb radiocarbon inventory because it assumes a constant air-sea $\cd$ disequilibrium (Mouchet, 2013). Yet, field reconstructions of the ocean bomb $\cq$ inventory are also biased low \citep{naegler_2009} since they assume that the anthropogenic perturbation did not affect ocean DIC since the pre-bomb epoch. For these reasons, bomb $\cq$ inventories obtained with the present method are directly comparable to reconstructions based on field measurements.
-
-This simplified approach also neglects the effects of fractionation (e.g.,  air-sea exchange) and of biological processes. Previous studies by \cite{bacastow_1990} and \cite{joos_1997} resulted in nearly identical $\Dcq$ distributions among experiments considering biology or not.
+This approach calls for a strong assumption, i.e., that of a homogeneous and constant dissolved inorganic carbon (DIC) field \citep{toggweiler_1989a,mouchet_2013}.
+While in terms of
+oceanic $\Dcq$, it yields similar results to approaches involving carbonate chemistry, it underestimates the bomb radiocarbon inventory because it assumes a constant air-sea $\cd$ disequilibrium (Mouchet, 2013).
+Yet, field reconstructions of the ocean bomb $\cq$ inventory are also biased low \citep{naegler_2009} since they assume that the anthropogenic perturbation did not affect ocean DIC since the pre-bomb epoch.
+For these reasons, bomb $\cq$ inventories obtained with the present method are directly comparable to reconstructions based on field measurements.
+
+This simplified approach also neglects the effects of fractionation (e.g.,  air-sea exchange) and of biological processes.
+Previous studies by \cite{bacastow_1990} and \cite{joos_1997} resulted in nearly identical $\Dcq$ distributions among experiments considering biology or not.
 Since observed $\Rq$ ratios are corrected for the isotopic fractionation when converted to the standard $\Dcq$ notation \citep{stuiver_1977} the model results are directly comparable to observations.
 
@@ -373,28 +413,37 @@
 
 The equation governing the transport of $\Rq$  in the ocean is
+
 \begin{equation}
 \frac{\partial}{\partial t} {\Rq} =  - \bigtriangledown \cdot ( \mathbf{u} \Rq - \mathbf{K} \cdot \bigtriangledown \Rq )  - \lambda \Rq, \label{eq:quick}
 \end{equation}
+
 where $\lambda$ is the radiocarbon decay rate, ${\mathbf{u}}$ the 3-D velocity field, and $\mathbf{K}$ the diffusivity tensor.
 
-At the air-sea interface a Robin boundary condition \citep{haine_2006} is applied to \eqref{eq:quick}, i.e., the flux
+At the air-sea interface a Robin boundary condition \citep{haine_2006} is applied to \autoref{eq:quick}, i.e., the flux
 through the interface is proportional to the difference in the ratios between
 the ocean and the atmosphere
+
 \begin{equation}
 \mathcal{\!F} =  \kappa_{R}  (\Rq  - \Rq_{a} ), \label{eq:BCR}
 \end{equation}
-where $\mathcal{\!F}$ is the flux out of the ocean, and $\Rq_{a}$ is the atmospheric $\cq/\mathrm{C}$ ratio. The transfer velocity $ \kappa_{R} $ for the radiocarbon ratio in \eqref{eq:BCR} is computed as
+
+where $\mathcal{\!F}$ is the flux out of the ocean, and $\Rq_{a}$ is the atmospheric $\cq/\mathrm{C}$ ratio.
+The transfer velocity $ \kappa_{R} $ for the radiocarbon ratio in \autoref{eq:BCR} is computed as
+
 \begin{equation}
  \kappa_{R} =  \frac{\kappa_{\cd} K_0}{\overline{\Ct}} \, \pacd   \label{eq:Rspeed}
 \end{equation}
+
 with $\kappa_{\cd}$ the carbon dioxide transfer or piston velocity, $K_0$ the $\cd$ solubility in seawater, $\pacd$ the atmospheric $\cd$ pressure at sea level, and $\overline{\Ct}$ the average sea-surface dissolved inorganic carbon concentration.
 
-
-The $\cd$ transfer velocity is based on the empirical formulation of \cite{wanninkhof_1992} with chemical enhancement \citep{wanninkhof_1996,wanninkhof_2014}. The original formulation is modified to account for the reduction of the  air-sea exchange rate in the presence of sea ice. Hence
+The $\cd$ transfer velocity is based on the empirical formulation of \cite{wanninkhof_1992} with chemical enhancement \citep{wanninkhof_1996,wanninkhof_2014}.
+The original formulation is modified to account for the reduction of the  air-sea exchange rate in the presence of sea ice.
+Hence
+
 \begin{equation}
 \kappa_\cd=\left( K_W\,\mathrm{w}^2 + b  \right)\, (1-f_\mathrm{ice})\,\sqrt{660/Sc}, \label{eq:wanc14}
 \end{equation}
 with $\mathrm{w}$ the wind magnitude, $f_\mathrm{ice}$ the fractional ice cover, and $Sc$ the Schmidt number.
-$K_W$ in \eqref{eq:wanc14} is an empirical coefficient with dimension of an inverse velocity.
+$K_W$ in \autoref{eq:wanc14} is an empirical coefficient with dimension of an inverse velocity.
 The chemical enhancement term $b$ is represented as a function of temperature $T$ \citep{wanninkhof_1992}
 \begin{equation}
@@ -402,5 +451,5 @@
 \end{equation}
 
-%We compare the results of equilibrium and transient experiments obtained with both methods in section \ref{sec:UNDEU}.
+%We compare the results of equilibrium and transient experiments obtained with both methods in section \autoref{sec:UNDEU}.
 
 %
@@ -413,35 +462,41 @@
 \label{sec:param}
  %
-The radiocarbon decay rate (\CODE{rlam14}; in \texttt{trcnam\_c14} module) is set to $\lambda=(1/8267)$ yr$^{-1}$ \citep{stuiver_1977}, which corresponds to a half-life of 5730 yr.\\[1pt]
-%
-The Schmidt number $Sc$, Eq. \eqref{eq:wanc14}, is calculated with the help of the formulation of \cite{wanninkhof_2014}. The $\cd$ solubility $K_0$ in \eqref{eq:Rspeed} is taken from \cite{weiss_1974}. $K_0$ and $Sc$ are computed with the OGCM temperature and salinity fields (\texttt{trcsms\_c14} module).\\[1pt]
+The radiocarbon decay rate (\forcode{rlam14}; in \texttt{trcnam\_c14} module) is set to $\lambda=(1/8267)$ yr$^{-1}$ \citep{stuiver_1977}, which corresponds to a half-life of 5730 yr.\\[1pt]
+%
+The Schmidt number $Sc$, Eq. \autoref{eq:wanc14}, is calculated with the help of the formulation of \cite{wanninkhof_2014}.
+The $\cd$ solubility $K_0$ in \autoref{eq:Rspeed} is taken from \cite{weiss_1974}. $K_0$ and $Sc$ are computed with the OGCM temperature and salinity fields (\texttt{trcsms\_c14} module).\\[1pt]
 %
 The following parameters intervening in the air-sea exchange rate are set in \texttt{namelist\_c14}:
+
 \begin{itemize}
-\item The reference DIC concentration $\overline{\Ct}$ (\CODE{xdicsur}) intervening in \eqref{eq:Rspeed} is classically set to 2 mol m$^{-3}$ \citep{toggweiler_1989a,orr_2001,butzin_2005}.
-%
-\item The value of the empirical coefficient $K_W$ (\CODE{xkwind}) in \eqref{eq:wanc14} depends on the wind field and on the model upper ocean mixing rate \citep{toggweiler_1989a,wanninkhof_1992,naegler_2009,wanninkhof_2014}.
+\item The reference DIC concentration $\overline{\Ct}$ (\forcode{xdicsur}) intervening in \autoref{eq:Rspeed} is classically set to 2 mol m$^{-3}$ \citep{toggweiler_1989a,orr_2001,butzin_2005}.
+%
+\item The value of the empirical coefficient $K_W$ (\forcode{xkwind}) in \autoref{eq:wanc14} depends on the wind field and on the model upper ocean mixing rate \citep{toggweiler_1989a,wanninkhof_1992,naegler_2009,wanninkhof_2014}.
 It should be adjusted so that the globally averaged $\cd$ piston velocity is $\kappa_\cd = 16.5\pm 3.2$ cm/h \citep{naegler_2009}.
-%The sensitivity to this parametrization is discussed in section \ref{sec:result}.
-%
-\item Chemical enhancement (term $b$  in Eq. \ref{eq:wanchem}) may be set on/off by means of the logical variable \CODE{ln\_chemh}.
+%The sensitivity to this parametrization is discussed in section \autoref{sec:result}.
+%
+\item Chemical enhancement (term $b$  in Eq. \autoref{eq:wanchem}) may be set on/off by means of the logical variable \forcode{ln\_chemh}.
 \end{itemize}
 
 %
 \paragraph{Experiment type}
-The type of experiment is determined by the value given to \CODE{kc14typ} in \texttt{namelist\_c14}. There are three possibilities:
+The type of experiment is determined by the value given to \forcode{kc14typ} in \texttt{namelist\_c14}.
+There are three possibilities:
+
 \begin{enumerate}
-\item natural $\Dcq$: \CODE{kc14typ}=0
-\item bomb $\Dcq$: \CODE{kc14typ}=1
-\item transient paleo-historical $\Dcq$: \CODE{kc14typ}=2
+\item natural                    $\Dcq$: \forcode{kc14typ}=0
+\item bomb                       $\Dcq$: \forcode{kc14typ}=1
+\item transient paleo-historical $\Dcq$: \forcode{kc14typ}=2
 \end{enumerate}
-%
+
+% 
 \textbf{Natural or Equilibrium radiocarbon}
-\CODE{kc14typ}=0
-
-Unless otherwise specified in \texttt{namelist\_c14}, the atmospheric $\Rq_a$ (\CODE{rc14at}) is set to one, the atmospheric $\cd$ (\CODE{pco2at}) to 280 ppm, and the ocean $\Rq$ is initialized with \CODE{rc14init=0.85}, i.e., $\Dcq=$-150\textperthousand  \cite[typical for deep-ocean, Fig 6 in][]{key_2004}.
-
-Equilibrium experiment should last until 98\% of the ocean volume exhibit a drift of less than 0.001\textperthousand/year \citep{orr_2000}; this is usually achieved after few kyr (Fig. \ref{fig:drift}).
-%
+\forcode{kc14typ}=0
+
+Unless otherwise specified in \texttt{namelist\_c14}, the atmospheric $\Rq_a$ (\forcode{rc14at}) is set to one, the atmospheric $\cd$ (\forcode{pco2at}) to 280 ppm, and the ocean $\Rq$ is initialized with \forcode{rc14init=0.85}, i.e., $\Dcq=$-150\textperthousand \cite[typical for deep-ocean, Fig 6 in][]{key_2004}.
+
+Equilibrium experiment should last until 98\% of the ocean volume exhibit a drift of less than 0.001\textperthousand/year \citep{orr_2000}; this is usually achieved after few kyr (Fig. \autoref{fig:drift}).
+%
+
 \begin{figure}[!h]
 \begin{center}
@@ -449,10 +504,12 @@
 \end{center}
 \vspace{-4ex}
-\caption{ Time evolution of $\Rq$ inventory anomaly for equilibrium run with homogeneous ocean initial state. The anomaly (or drift) is given in \%  change in total ocean inventory per 50 years. Time on x-axis is in simulation year.\label{fig:drift} }
+\caption{ Time evolution of $\Rq$ inventory anomaly for equilibrium run with homogeneous ocean initial state.
+The anomaly (or drift) is given in \%  change in total ocean inventory per 50 years.
+Time on x-axis is in simulation year.\label{fig:drift} }
 \end{figure}
 
 \textbf{Transient: Bomb}
 \label{sec:bomb}
-\CODE{kc14typ}=1
+\forcode{kc14typ}=1
 
 \begin{figure}[!h]
@@ -461,27 +518,35 @@
 \end{center}
 \vspace{-4ex}
-\caption{Atmospheric $\Dcq$ (solid; left axis) and $\cd$ (dashed; right axis)  forcing for the $\cq$-bomb experiments. The $\Dcq$ is illustrated for the three zonal bands (upper, middle, and lower curves correspond to latitudes $> 20$N, $\in [20\mathrm{S},20\mathrm{N}]$, and $< 20$S, respectively.} \label{fig:bomb}
+\caption{Atmospheric $\Dcq$ (solid; left axis) and $\cd$ (dashed; right axis)  forcing for the $\cq$-bomb experiments.
+The $\Dcq$ is illustrated for the three zonal bands (upper, middle, and lower curves correspond to latitudes $> 20$N, $\in [20\mathrm{S},20\mathrm{N}]$, and $< 20$S, respectively.} \label{fig:bomb}
 \end{figure}
 
-Performing this type of experiment requires that a pre-industrial equilibrium run be performed beforehand (\CODE{ln\_rsttr} should be set to \texttt{.TRUE.}).
-
-An exception to this rule is when wishing to perform a perturbation bomb experiment as was possible with the package \texttt{C14b}. It is still possible to easily set-up that type of transient experiment for which no previous run is needed.  In addition to the instructions as given in this section it is however necessary to adapt the \texttt{atmc14.dat} file so that it does no longer contain any negative $\Dcq$ values (Suess effect in the pre-bomb period).
-
-The model  is integrated from a given initial date following the observed records provided from 1765 AD on ( Fig. \ref{fig:bomb}).
-The file \texttt{atmc14.dat}  \cite[][\& I. Levin, personal comm.]{enting_1994} provides atmospheric $\Dcq$ for three latitudinal bands: 90S-20S,    20S-20N \&    20N-90N.
-Atmospheric $\cd$ in the file \texttt{splco2.dat} is obtained from a spline fit through ice core data and direct atmospheric measurements \cite[][\& J. Orr, personal comm.]{orr_2000}.
+Performing this type of experiment requires that a pre-industrial equilibrium run be performed beforehand (\forcode{ln\_rsttr} should be set to \texttt{.TRUE.}).
+
+An exception to this rule is when wishing to perform a perturbation bomb experiment as was possible with the package \texttt{C14b}.
+It is still possible to easily set-up that type of transient experiment for which no previous run is needed.
+In addition to the instructions as given in this section it is however necessary to adapt the \texttt{atmc14.dat} file so that it does no longer contain any negative $\Dcq$ values (Suess effect in the pre-bomb period).
+
+The model  is integrated from a given initial date following the observed records provided from 1765 AD on ( Fig. \autoref{fig:bomb}).
+The file \texttt{atmc14.dat}  \cite[][\& I.
+Levin, personal comm.]{enting_1994} provides atmospheric $\Dcq$ for three latitudinal bands: 90S-20S,    20S-20N \&    20N-90N.
+Atmospheric $\cd$ in the file \texttt{splco2.dat} is obtained from a spline fit through ice core data and direct atmospheric measurements \cite[][\& J.
+Orr, personal comm.]{orr_2000}.
 Dates in these forcing files are expressed as yr AD.
 
 To ensure that the atmospheric forcing is applied properly as well as that output files contain consistent dates and inventories the experiment should be set up carefully:
+
 \begin{itemize}
-\item Specify the starting date of the experiment: \CODE{nn\_date0} in \texttt{namelist}.  \CODE{nn\_date0} is written as Year0101 where Year may take any positive value (AD).
-\item Then the parameters \CODE{nn\_rstctl} in  \texttt{namelist} (on-line) and \CODE{nn\_rsttr} in \texttt{namelist\_top} (off-line)  must be \textbf{set to 0} at the start of the experiment (force the date to \CODE{nn\_date0} for the \textbf{first} experiment year).
-\item These two parameters (\CODE{nn\_rstctl} and \CODE{nn\_rsttr}) have then to be \textbf{set to 2} for the following years (the date must be read in the restart file).
+\item Specify the starting date of the experiment: \forcode{nn\_date0} in \texttt{namelist}.  \forcode{nn\_date0} is written as Year0101 where Year may take any positive value (AD).
+\item Then the parameters \forcode{nn\_rstctl} in  \texttt{namelist} (on-line) and \forcode{nn\_rsttr} in \texttt{namelist\_top} (off-line)  must be \textbf{set to 0} at the start of the experiment (force the date to \forcode{nn\_date0} for the \textbf{first} experiment year).
+\item These two parameters (\forcode{nn\_rstctl} and \forcode{nn\_rsttr}) have then to be \textbf{set to 2} for the following years (the date must be read in the restart file).
 \end{itemize}
- If the experiment date is outside the data time span then the first or last atmospheric concentrations are prescribed depending on whether the date is earlier or later. Note that \CODE{tyrc14\_beg} (\texttt{namelist\_c14}) is not used in this context.
+
+If the experiment date is outside the data time span then the first or last atmospheric concentrations are prescribed depending on whether the date is earlier or later.
+Note that \forcode{tyrc14\_beg} (\texttt{namelist\_c14}) is not used in this context.
 
 %
 \textbf{Transient: Past}
-\CODE{kc14typ}=2
+\forcode{kc14typ}=2
 %
 \begin{figure}[!h]
@@ -490,43 +555,52 @@
 \end{center}
 \vspace{-4ex}
-\caption{Atmospheric $\Dcq$ (solid) and $\cd$ (dashed)  forcing for the Paleo experiments. The $\cd$ scale is given on the right axis.} \label{fig:paleo}
+\caption{Atmospheric $\Dcq$ (solid) and $\cd$ (dashed)  forcing for the Paleo experiments.
+The $\cd$ scale is given on the right axis.} \label{fig:paleo}
 \end{figure}
 
-This experiment type does not need a previous equilibrium run. It should start 5--6 kyr earlier than the period to be analyzed.
-Atmospheric $\Rq_a$ and $\cd$ are prescribed from forcing files. The ocean $\Rq$ is initialized with the value attributed to \CODE{rc14init} in \texttt{namelist\_c14}.
+This experiment type does not need a previous equilibrium run.
+It should start 5--6 kyr earlier than the period to be analyzed.
+Atmospheric $\Rq_a$ and $\cd$ are prescribed from forcing files.
+The ocean $\Rq$ is initialized with the value attributed to \forcode{rc14init} in \texttt{namelist\_c14}.
 
 The file \texttt{intcal13.14c} \citep{reimer_2013} contains atmospheric $\Dcq$ from 0 to 50 kyr cal BP\footnote{cal BP: number of years before 1950 AD}.
-The $\cd$ forcing is provided in file \texttt{ByrdEdcCO2.txt}. The content of this file is based on  the high resolution record from EPICA Dome C \citep{monnin_2004} for the Holocene and the Transition, and on Byrd Ice Core CO2 Data for 20--90 kyr BP  \citep{ahn_2008}. These atmospheric values are reproduced in Fig. \ref{fig:paleo}. Dates in these files are expressed as yr BP.
+The $\cd$ forcing is provided in file \texttt{ByrdEdcCO2.txt}.
+The content of this file is based on  the high resolution record from EPICA Dome C \citep{monnin_2004} for the Holocene and the Transition, and on Byrd Ice Core CO2 Data for 20--90 kyr BP  \citep{ahn_2008}.
+These atmospheric values are reproduced in Fig. \autoref{fig:paleo}.
+Dates in these files are expressed as yr BP.
 
 To ensure that the atmospheric forcing is applied properly as well as that output files contain consistent dates and inventories the experiment should be set up carefully.
-The true experiment starting date is given by \CODE{tyrc14\_beg} (in yr BP) in \texttt{namelist\_c14}. In consequence, \CODE{nn\_date0} in \texttt{namelist} MUST be set to 00010101.\\
-Then the parameters \CODE{nn\_rstctl} in  \texttt{namelist} (on-line) and \CODE{nn\_rsttr} in \texttt{namelist\_top} (off-line)  must be set to 0 at the start of the experiment (force the date to \CODE{nn\_date0} for the first experiment year). These two parameters have then to be set to 2 for the following years (read the date in the restart file). \\
- If the experiment date is outside the data time span then the first or last atmospheric concentrations are prescribed depending on whether the date is earlier or later.
+The true experiment starting date is given by \forcode{tyrc14\_beg} (in yr BP) in \texttt{namelist\_c14}.
+In consequence, \forcode{nn\_date0} in \texttt{namelist} MUST be set to 00010101.\\
+Then the parameters \forcode{nn\_rstctl} in  \texttt{namelist} (on-line) and \forcode{nn\_rsttr} in \texttt{namelist\_top} (off-line)  must be set to 0 at the start of the experiment (force the date to \forcode{nn\_date0} for the first experiment year).
+These two parameters have then to be set to 2 for the following years (read the date in the restart file). \\
+If the experiment date is outside the data time span then the first or last atmospheric concentrations are prescribed depending on whether the date is earlier or later.
 
 %
 \paragraph{Model output}
 \label{sec:output}
-All output fields in Table \ref{tab:out} are routinely computed. It depends on the actual settings in \texttt{iodef.xml} whether they are stored or not.
+
+All output fields in Table \autoref{tab:out} are routinely computed.
+It depends on the actual settings in \texttt{iodef.xml} whether they are stored or not.
 %
 \begin{table}[!h]
 \begin{center}
-\caption{Standard output fields for the C14 package \label{tab:out}.
-}
+\caption{Standard output fields for the C14 package \label{tab:out}.}
 %\begin{small}
 \renewcommand{\arraystretch}{1.3}%
 \begin{tabular}[h]{|l*{3}{|c}|l|}
 \hline
-Field & Type & Dim & Units & Description \\ \hline
-RC14 & ptrc & 3-D & -        & Radiocarbon ratio \\
-DeltaC14 & diad & 3-D & \textperthousand & $\Dcq$\\
-C14Age & diad & 3-D & yr &   Radiocarbon age \\
-RAge & diad & 2-D & yr & Reservoir age\\
-qtr\_c14 &  diad & 2-D & m$^{-2}$ yr$^{-1}$ & Air-to-sea net $\Rq$ flux\\
-qint\_c14 & diad & 2-D &   m$^{-2}$ &  Cumulative air-to-sea $\Rq$ flux \\
-AtmCO2 & scalar & 0-D & ppm & Global atmospheric $\cd$ \\
-AtmC14 & scalar & 0-D & \textperthousand  & Global atmospheric $\Dcq$\\
-K\_CO2 & scalar & 0-D & cm h$^{-1}$  & Global $\cd$ piston velocity ($ \overline{\kappa_{\cd}}$) \\
-K\_C14 & scalar & 0-D &m yr$^{-1}$ & Global $\Rq$ transfer velocity  ($ \overline{\kappa_R}$)\\
-C14Inv & scalar & 0-D & $10^{26}$ atoms & Ocean radiocarbon inventory \\ \hline
+Field     & Type   & Dim & Units              & Description                                               \\ \hline
+RC14      & ptrc   & 3-D & -                  & Radiocarbon ratio                                         \\
+DeltaC14  & diad   & 3-D & \textperthousand   & $\Dcq$                                                    \\
+C14Age    & diad   & 3-D & yr                 & Radiocarbon age                                           \\
+RAge      & diad   & 2-D & yr                 & Reservoir age                                             \\
+qtr\_c14  & diad   & 2-D & m$^{-2}$ yr$^{-1}$ & Air-to-sea net $\Rq$ flux                                 \\
+qint\_c14 & diad   & 2-D & m$^{-2}$           & Cumulative air-to-sea $\Rq$ flux                          \\
+AtmCO2    & scalar & 0-D & ppm                & Global atmospheric $\cd$                                  \\
+AtmC14    & scalar & 0-D & \textperthousand   & Global atmospheric $\Dcq$                                 \\
+K\_CO2    & scalar & 0-D & cm h$^{-1}$        & Global $\cd$ piston velocity ($ \overline{\kappa_{\cd}}$) \\
+K\_C14    & scalar & 0-D & m yr$^{-1}$        & Global $\Rq$ transfer velocity  ($ \overline{\kappa_R}$)  \\
+C14Inv    & scalar & 0-D & $10^{26}$ atoms    & Ocean radiocarbon inventory                               \\ \hline
 \end{tabular}
 %\end{small}
@@ -539,20 +613,30 @@
 The radiocarbon age is computed as  $(-1/\lambda) \ln{ \left( \Rq \right)}$, with zero age corresponding to $\Rq=1$.
 
-The reservoir age is the age difference between the ocean uppermost layer and the atmosphere. It is usually reported as conventional radiocarbon age; i.e., computed by means of the Libby radiocarbon mean life \cite[8033 yr;][]{stuiver_1977}
+The reservoir age is the age difference between the ocean uppermost layer and the atmosphere.
+It is usually reported as conventional radiocarbon age; i.e., computed by means of the Libby radiocarbon mean life \cite[8033 yr;][]{stuiver_1977}
+
 \begin{align}
 {^{14}\tau_\mathrm{c}}= -8033 \; \ln \left(1 + \frac{\Dcq}{10^3}\right), \label{eq:convage}
 \end{align}
-where ${^{14}\tau_\mathrm{c}}$ is expressed in years B.P. Here we do not use that convention and compute reservoir ages with the mean lifetime $1/\lambda$. Conversion from one scale to the other is readily performed. The conventional radiocarbon age is lower than the radiocarbon age by $\simeq3\%$.
+
+where ${^{14}\tau_\mathrm{c}}$ is expressed in years B.P.
+Here we do not use that convention and compute reservoir ages with the mean lifetime $1/\lambda$.
+Conversion from one scale to the other is readily performed.
+The conventional radiocarbon age is lower than the radiocarbon age by $\simeq3\%$.
 
 The ocean radiocarbon  inventory is computed as
+
 \begin{equation}
 N_A \Rq_\mathrm{oxa} \overline{\Ct} \left( \int_\Omega \Rq d\Omega \right) /10^{26}, \label{eq:inv}
 \end{equation}
-where $N_A$ is the Avogadro's number ($N_A=6.022\times10^{23}$ at/mol), $\Rq_\mathrm{oxa}$ is the oxalic acid radiocarbon standard \cite[$\Rq_\mathrm{oxa}=1.176\times10^{-12}$;][]{stuiver_1977}, and $\Omega$ is the ocean volume.  Bomb $\cq$ inventories are traditionally reported in units of $10^{26}$ atoms, hence the denominator in \eqref{eq:inv}.
-
-All transformations from second to year, and inversely, are performed with the help of the physical constant \CODE{rsiyea} the sideral year length expressed in seconds\footnote{The variable (\CODE{nyear\_len}) which reports the length in days of the previous/current/future year (see \textrm{oce\_trc.F90}) is not a constant. }.
+
+where $N_A$ is the Avogadro's number ($N_A=6.022\times10^{23}$ at/mol), $\Rq_\mathrm{oxa}$ is the oxalic acid radiocarbon standard \cite[$\Rq_\mathrm{oxa}=1.176\times10^{-12}$;][]{stuiver_1977}, and $\Omega$ is the ocean volume.
+Bomb $\cq$ inventories are traditionally reported in units of $10^{26}$ atoms, hence the denominator in \autoref{eq:inv}.
+
+All transformations from second to year, and inversely, are performed with the help of the physical constant \forcode{rsiyea} the sideral year length expressed in seconds\footnote{The variable (\forcode{nyear\_len}) which reports the length in days of the previous/current/future year (see \textrm{oce\_trc.F90}) is not a constant. }.
 
 The global transfer velocities represent time-averaged\footnote{the actual duration is set in \texttt{iodef.xml}} global integrals of the transfer rates:
- \begin{equation}
+
+\begin{equation}
  \overline{\kappa_{\cd}}= \int_\mathcal{S} \kappa_{\cd} d\mathcal{S}  \text{ and } \overline{\kappa_R}= \int_\mathcal{S} \kappa_R d\mathcal{S}
 \end{equation}
@@ -561,19 +645,41 @@
 \subsection{PISCES biogeochemical model}
 
-PISCES is a biogeochemical model which simulates the lower trophic levels of marine ecosystem (phytoplankton, microzooplankton and mesozooplankton) and the biogeochemical cycles of carbonand of the main nutrients (P, N, Fe, and Si). The  model is intended to be used for both regional and global configurations at high or low spatial resolutions as well as for  short-term (seasonal, interannual) and long-term (climate change, paleoceanography) analyses.
+PISCES is a biogeochemical model which simulates the lower trophic levels of marine ecosystem (phytoplankton, microzooplankton and mesozooplankton) and the biogeochemical cycles of carbonand of the main nutrients (P, N, Fe, and Si).
+The  model is intended to be used for both regional and global configurations at high or low spatial resolutions as well as for  short-term (seasonal, interannual) and long-term (climate change, paleoceanography) analyses.
 Two versions of PISCES are available in NEMO v4.0 :
 
-PISCES-v2, by setting in namelist\_pisces\_ref  \np{ln\_p4z} to true,  can be seen as one of the many Monod models \citep{monod_1958}. It assumes a constant Redfield ratio and phytoplankton growth depends on the external concentration in nutrients. There are twenty-four prognostic variables (tracers) including two phytoplankton compartments  (diatoms and nanophytoplankton), two zooplankton size-classes (microzooplankton and  mesozooplankton) and a description of the carbonate chemistry. Formulations in PISCES-v2 are based on a mixed Monod/Quota formalism: On one hand, stoichiometry of C/N/P is fixed and growth rate of phytoplankton is limited by the external availability in N, P and Si. On the other hand, the iron and silicium quotas are variable and growth rate of phytoplankton is limited by the internal availability in Fe. Various parameterizations can be activated in PISCES-v2, setting for instance the complexity of iron chemistry or the description of particulate organic materials.
-
-PISCES-QUOTA has been built on the PISCES-v2 model described in \citet{aumont_2015}. PISCES-QUOTA has thirty-nine prognostic compartments. Phytoplankton growth can be controlled by five modeled limiting nutrients: Nitrate and Ammonium, Phosphate, Silicate and Iron. Five living compartments are represented: Three phytoplankton size classes/groups corresponding to picophytoplankton, nanophytoplankton and diatoms, and two zooplankton size classes which are microzooplankton and mesozooplankton. For phytoplankton, the prognostic variables are the carbon, nitrogen, phosphorus,  iron, chlorophyll and silicon biomasses (the latter only for diatoms). This means that the N/C, P/C, Fe/C and Chl/C ratios of both phytoplankton groups as well as the Si/C ratio of diatoms are prognostically predicted  by the model. Zooplankton are assumed to be strictly homeostatic \citep[e.g.,][]{sterner_2003,woods_2013,meunier_2014}. As a consequence, the C/N/P/Fe ratios of these groups are maintained constant and are not allowed to vary. In PISCES, the Redfield ratios C/N/P are set to 122/16/1 \citep{takahashi_1985} and the -O/C ratio is set to 1.34 \citep{kortzinger_2001}. No silicified zooplankton is assumed. The bacterial pool is not yet explicitly modeled.
-
-There are three non-living compartments: Semi-labile dissolved organic matter, small sinking particles, and large sinking particles. As a consequence of the variable stoichiometric ratios of phytoplankton and of the stoichiometric regulation of zooplankton, elemental ratios in organic matter cannot be supposed constant anymore as that was the case in PISCES-v2. Indeed, the nitrogen, phosphorus, iron, silicon and calcite pools of the particles are now all explicitly modeled. The sinking speed of the particles is not altered by their content in calcite and biogenic silicate (''The ballast effect'', \citep{honjo_1996,armstrong_2001}). The latter particles are assumed to sink at the same speed as the large organic matter particles. All the non-living compartments experience aggregation due to turbulence and differential settling as well as Brownian coagulation for DOM.
-
+PISCES-v2, by setting in namelist\_pisces\_ref  \np{ln\_p4z} to true,  can be seen as one of the many Monod models \citep{monod_1958}.
+It assumes a constant Redfield ratio and phytoplankton growth depends on the external concentration in nutrients.
+There are twenty-four prognostic variables (tracers) including two phytoplankton compartments  (diatoms and nanophytoplankton), two zooplankton size-classes (microzooplankton and  mesozooplankton) and a description of the carbonate chemistry.
+Formulations in PISCES-v2 are based on a mixed Monod/Quota formalism: On one hand, stoichiometry of C/N/P is fixed and growth rate of phytoplankton is limited by the external availability in N, P and Si.
+On the other hand, the iron and silicium quotas are variable and growth rate of phytoplankton is limited by the internal availability in Fe.
+Various parameterizations can be activated in PISCES-v2, setting for instance the complexity of iron chemistry or the description of particulate organic materials.
+
+PISCES-QUOTA has been built on the PISCES-v2 model described in \citet{aumont_2015}.
+PISCES-QUOTA has thirty-nine prognostic compartments.
+Phytoplankton growth can be controlled by five modeled limiting nutrients: Nitrate and Ammonium, Phosphate, Silicate and Iron.
+Five living compartments are represented: Three phytoplankton size classes/groups corresponding to picophytoplankton, nanophytoplankton and diatoms, and two zooplankton size classes which are microzooplankton and mesozooplankton.
+For phytoplankton, the prognostic variables are the carbon, nitrogen, phosphorus,  iron, chlorophyll and silicon biomasses (the latter only for diatoms).
+This means that the N/C, P/C, Fe/C and Chl/C ratios of both phytoplankton groups as well as the Si/C ratio of diatoms are prognostically predicted  by the model.
+Zooplankton are assumed to be strictly homeostatic \citep[e.g.,][]{sterner_2003,woods_2013,meunier_2014}.
+As a consequence, the C/N/P/Fe ratios of these groups are maintained constant and are not allowed to vary.
+In PISCES, the Redfield ratios C/N/P are set to 122/16/1 \citep{takahashi_1985} and the -O/C ratio is set to 1.34 \citep{kortzinger_2001}.
+No silicified zooplankton is assumed.
+The bacterial pool is not yet explicitly modeled.
+
+There are three non-living compartments: Semi-labile dissolved organic matter, small sinking particles, and large sinking particles.
+As a consequence of the variable stoichiometric ratios of phytoplankton and of the stoichiometric regulation of zooplankton, elemental ratios in organic matter cannot be supposed constant anymore as that was the case in PISCES-v2.
+Indeed, the nitrogen, phosphorus, iron, silicon and calcite pools of the particles are now all explicitly modeled.
+The sinking speed of the particles is not altered by their content in calcite and biogenic silicate (''The ballast effect'', \citep{honjo_1996,armstrong_2001}).
+The latter particles are assumed to sink at the same speed as the large organic matter particles.
+All the non-living compartments experience aggregation due to turbulence and differential settling as well as Brownian coagulation for DOM.
 
 \subsection{MY\_TRC interface for coupling external BGC models}
 \label{Mytrc}
 
-The NEMO-TOP has only one built-in biogeochemical model - PISCES - but there are several BGC models - MEDUSA, ERSEM, BFM or ECO3M - which are meant to be coupled with the NEMO dynamics. Therefore it was necessary to provide to the users a framework for easily add their own BGCM model, that can be a single passive tracer.
-The generalized interface is pivoted on MY\_TRC module that contains template files to build the coupling between NEMO and any external BGC model. The call to MY\_TRC is activated by setting  \textit{ln\_my\_trc} = \texttt{.true.} in namelist \textit{namtrc}
+The NEMO-TOP has only one built-in biogeochemical model - PISCES - but there are several BGC models - MEDUSA, ERSEM, BFM or ECO3M - which are meant to be coupled with the NEMO dynamics.
+Therefore it was necessary to provide to the users a framework for easily add their own BGCM model, that can be a single passive tracer.
+The generalized interface is pivoted on MY\_TRC module that contains template files to build the coupling between NEMO and any external BGC model.
+The call to MY\_TRC is activated by setting  \textit{ln\_my\_trc} = \texttt{.true.} in namelist \textit{namtrc}
 
 The following 6 fortran files are available in MY\_TRC with the specific purposes here described.
@@ -581,11 +687,15 @@
 \begin{itemize}
    \item \textit{par\_my\_trc.F90} :  This module allows to define additional arrays and public variables to be used within the MY\_TRC interface
-   \item \textit{trcini\_my\_trc.F90} :  Here are initialized user defined namelists and the call to the external BGC model initialization procedures to populate general tracer array (trn and trb). Here are also likely to be defined suport arrays related to system metrics that could be needed by the BGC model.
+   \item \textit{trcini\_my\_trc.F90} :  Here are initialized user defined namelists and the call to the external BGC model initialization procedures to populate general tracer array (trn and trb).
+Here are also likely to be defined suport arrays related to system metrics that could be needed by the BGC model.
   \item \textit{trcnam\_my\_trc.F90} :  This routine is called at the beginning of trcini\_my\_trc and should contain the initialization of additional namelists for the BGC model or user-defined code.
-  \item \textit{trcsms\_my\_trc.F90} :  The routine performs the call to Boundary Conditions and its main purpose is to contain the Source-Minus-Sinks terms due to the biogeochemical processes of the external model. Be aware that lateral boundary conditions are applied in trcnxt routine. IMPORTANT: the routines to compute the light penetration along the water column and the tracer vertical sinking should be defined/called in here, as generalized modules are still missing in the code.
- \item \textit{trcice\_my\_trc.F90} : Here it is possible to prescribe the tracers concentrations in the seaice that will be used as boundary conditions when ice melting occurs (nn\_ice\_tr =1 in namtrc\_ice). See e.g. the correspondent PISCES subroutine.
- \item \textit{trcwri\_my\_trc.F90} : This routine performs the output of the model tracers using IOM module (see Manual Chapter Output and Diagnostics). It is possible to place here the output of additional variables produced by the model, if not done elsewhere in the code, using the call to \textit{iom\_put}.
+  \item \textit{trcsms\_my\_trc.F90} :  The routine performs the call to Boundary Conditions and its main purpose is to contain the Source-Minus-Sinks terms due to the biogeochemical processes of the external model.
+Be aware that lateral boundary conditions are applied in trcnxt routine.
+IMPORTANT: the routines to compute the light penetration along the water column and the tracer vertical sinking should be defined/called in here, as generalized modules are still missing in the code.
+ \item \textit{trcice\_my\_trc.F90} : Here it is possible to prescribe the tracers concentrations in the seaice that will be used as boundary conditions when ice melting occurs (nn\_ice\_tr =1 in namtrc\_ice).
+See e.g. the correspondent PISCES subroutine.
+ \item \textit{trcwri\_my\_trc.F90} : This routine performs the output of the model tracers using IOM module (see Manual Chapter Output and Diagnostics).
+It is possible to place here the output of additional variables produced by the model, if not done elsewhere in the code, using the call to \textit{iom\_put}.
 \end{itemize}
-
 
 \section{The Offline Option}
@@ -596,9 +706,12 @@
 %-------------------------------------------------------------------------------------------------------------
 
-Coupling passive tracers offline with NEMO requires precomputed  physical fields from OGCM. Those fields are read from files and interpolated on-the-fly at each model time step
+Coupling passive tracers offline with NEMO requires precomputed  physical fields from OGCM.
+Those fields are read from files and interpolated on-the-fly at each model time step
 At least the following dynamical parameters should be absolutely passed to the transport : ocean velocities, temperature, salinity, mixed layer depth and for ecosystem models like PISCES, sea ice concentration, short wave radiation at the ocean surface, wind speed (or at least, wind stress).
 The so-called offline mode is useful since it has lower computational costs for example to perform very longer simulations - about 3000 years - to reach equilibrium of CO2 sinks for climate-carbon studies.
 
-The offline interface is located in the code repository : \path{<repository>/src/OFF/}. It is activated by adding the CPP key  \textit{key\_offline} to the CPP keys list. There are two specifics routines for the Offline code :
+The offline interface is located in the code repository : \path{<repository>/src/OFF/}.
+It is activated by adding the CPP key  \textit{key\_offline} to the CPP keys list.
+There are two specifics routines for the Offline code :
 
 \begin{itemize}
@@ -606,5 +719,4 @@
    \item \textit{nemogcm.F90} :  a degraded version of the main nemogcm.F90 code of NEMO to manage the time-stepping
 \end{itemize}
-
 
 %-
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/global/latexmk.pl
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/global/latexmk.pl	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/latex/global/latexmk.pl	(revision 14403)
@@ -9,5 +9,9 @@
 
 ## Custom cmds
-set_tex_cmds('-shell-escape -interaction=batchmode');
-#set_tex_cmds('-shell-escape');
-$makeindex = 'makeindex -s %R.ist %O -o %D %S';
+set_tex_cmds('-shell-escape -file-line-error -interaction=batchmode');
+#set_tex_cmds('-shell-escape -file-line-error');
+$makeindex = 'makeindex %O -s %R.ist -o %D %S';
+## %D: Destination file (.ind for index)
+## %O: Options
+## %R: Root filename (${model}_manual)
+## %S: Source file      (.idx for index)
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/doc/rst/source/conf.py
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/doc/rst/source/conf.py	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/doc/rst/source/conf.py	(revision 14403)
@@ -279,2 +279,5 @@
 # Default language to highlight set to fortran
 highlight_language = 'fortran'
+
+# Extra setting for sphinxcontrib.bibtex upgrade to 2.X.X branche
+bibtex_bibfiles = ['cfgs.bib', 'ref.bib', 'tests.bib', 'zooms.bib']
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbdia.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbdia.F90	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbdia.F90	(revision 14403)
@@ -86,5 +86,5 @@
    INTEGER                       ::  nbergs_start, nbergs_end, nbergs_calved
    INTEGER                       ::  nbergs_melted
-   INTEGER                       ::  nspeeding_tickets
+   INTEGER                       ::  nspeeding_tickets, nspeeding_tickets_all
    INTEGER , DIMENSION(nclasses) ::  nbergs_calved_by_class
 
@@ -125,4 +125,5 @@
       nbergs_calved_by_class(:) = 0
       nspeeding_tickets         = 0
+      nspeeding_tickets_all     = 0
       stored_heat_end           = 0._wp
       floating_heat_end         = 0._wp
@@ -271,8 +272,8 @@
             CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+4), nclasses+4 )
             !
-            nbergs_end        = nsumbuf(1)
-            nbergs_calved     = nsumbuf(2)
-            nbergs_melted     = nsumbuf(3)
-            nspeeding_tickets = nsumbuf(4)
+            nbergs_end            = nsumbuf(1)
+            nbergs_calved         = nsumbuf(2)
+            nbergs_melted         = nsumbuf(3)
+            nspeeding_tickets_all = nsumbuf(4)
             DO ik = 1,nclasses
                nbergs_calved_by_class(ik)= nsumbuf(4+ik)
@@ -329,5 +330,8 @@
          IF (nn_verbose_level > 0) THEN
             WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses)
-            IF( nspeeding_tickets > 0 )   WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets
+            IF( nspeeding_tickets_all > 0 ) THEN
+                WRITE( numicb, '("speeding tickets issued (this domain)  = ",i6)') nspeeding_tickets
+                WRITE( numicb, '("speeding tickets issued (all domains)  = ",i6)') nspeeding_tickets_all
+            END IF
          ENDIF
          !
@@ -338,4 +342,5 @@
          nbergs_calved_by_class(:) = 0
          nspeeding_tickets         = 0
+         nspeeding_tickets_all     = 0
          stored_heat_start         = stored_heat_end
          floating_heat_start       = floating_heat_end
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbdyn.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbdyn.F90	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbdyn.F90	(revision 14403)
@@ -85,6 +85,6 @@
 
          !                                         !**   A1 = A(X1,V1)
-         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
-            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
+         CALL icb_accel( kt, berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
+            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 0.5_wp )
          !
          zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
@@ -102,6 +102,6 @@
 
          !                                         !**   A2 = A(X2,V2)
-         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
-            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
+         CALL icb_accel( kt, berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
+            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 0.5_wp )
          !
          zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
@@ -114,10 +114,10 @@
          zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
          !
-         CALL icb_ground( berg, zxi3, zxi1, zu3,   &
-            &                   zyj3, zyj1, zv3, ll_bounced )
+         CALL icb_ground( berg, zxi3, zxi1, zu2,   &
+            &                   zyj3, zyj1, zv2, ll_bounced )
 
          !                                         !**   A3 = A(X3,V3)
-         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
-            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
+         CALL icb_accel( kt, berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
+            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 1._wp )
          !
          zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
@@ -130,10 +130,10 @@
          zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
 
-         CALL icb_ground( berg, zxi4, zxi1, zu4,   &
-            &                   zyj4, zyj1, zv4, ll_bounced )
+         CALL icb_ground( berg, zxi4, zxi1, zu3,   &
+            &                   zyj4, zyj1, zv3, ll_bounced )
 
          !                                         !**   A4 = A(X4,V4)
-         CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
-            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
+         CALL icb_accel( kt, berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
+            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 1._wp )
 
          zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
@@ -255,6 +255,6 @@
 
 
-   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
-      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
+   SUBROUTINE icb_accel( kt, berg , pxi, pe1, puvel, puvel0, pax,                 &
+      &                             pyj, pe2, pvvel, pvvel0, pay, pdt, pcfl_scale )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_accel  ***
@@ -265,4 +265,6 @@
       !!----------------------------------------------------------------------
       TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
+      INTEGER                , INTENT(in   ) ::   kt               ! time step
+      REAL(wp)               , INTENT(in   ) ::   pcfl_scale
       REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
       REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
@@ -404,9 +406,20 @@
          zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
          IF( zspeed > 0._wp ) THEN
-            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
-            zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
+            zloc_dx = MIN( pe1, pe2 )                                ! minimum grid spacing
+            ! cfl scale is function of the RK4 step
+            zspeed_new = zloc_dx / pdt * rn_speed_limit * pcfl_scale ! Speed limit as a factor of dx / dt
             IF( zspeed_new < zspeed ) THEN
-               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
-               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
+               zuveln = zuveln * ( zspeed_new / zspeed )             ! Scale velocity to reduce speed
+               zvveln = zvveln * ( zspeed_new / zspeed )             ! without changing the direction
+               pax = (zuveln - puvel0)/pdt
+               pay = (zvveln - pvvel0)/pdt
+               !
+               ! print speeding ticket
+               IF (nn_verbose_level > 0) THEN
+                  WRITE(numicb, 9200) 'icb speeding : ',kt, nknberg, zspeed, &
+                       &                pxi, pyj, zuo, zvo, zua, zva, zui, zvi
+                  9200 FORMAT(a,i9,i6,f9.2,1x,4(1x,2f9.2))
+               END IF
+               !
                CALL icb_dia_speed()
             ENDIF
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbutl.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbutl.F90	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/ICB/icbutl.F90	(revision 14403)
@@ -320,9 +320,12 @@
          !
          IF ( ierr > 0 ) THEN
-            WRITE(numout,*) 'bottom left corner T point out of bound'
-            WRITE(numout,*) pi, kii, mig( 1 ), mig(jpi)
-            WRITE(numout,*) pj, kij, mjg( 1 ), mjg(jpj)
-            WRITE(numout,*) pmsk
-            CALL ctl_stop('STOP','icb_utl_bilin_h: an icebergs coordinates is out of valid range (out of bound error)')
+            WRITE(numicb,*) 'bottom left corner T point out of bound'
+            WRITE(numicb,*) pi, kii, mig( 1 ), mig(jpi)
+            WRITE(numicb,*) pj, kij, mjg( 1 ), mjg(jpj)
+            WRITE(numicb,*) pmsk
+            CALL FLUSH(numicb)
+            CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error).'       , &
+                 &                                'This can be fixed using rn_speed_limit=0.4 in &namberg.'                   , &
+                 &                                'More details in the corresponding iceberg.stat file (nn_verbose_level > 0).' )
          END IF
       END IF
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/LBC/lib_mpp.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/LBC/lib_mpp.F90	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/LBC/lib_mpp.F90	(revision 14403)
@@ -1325,5 +1325,5 @@
          END DO
          IF ( crname_lbc(n_sequence_lbc) /= 'already counted' ) THEN
-            WRITE(numcom,'(A, I4, A, A)') ' - ', 1,' times by subroutine ', TRIM(crname_lbc(ncom_rec_max))
+            WRITE(numcom,'(A, I4, A, A)') ' - ', 1,' times by subroutine ', TRIM(crname_lbc(n_sequence_lbc))
          END IF
          WRITE(numcom,*) ' '
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcblk.F90	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcblk.F90	(revision 14403)
@@ -40,5 +40,5 @@
    USE sbcdcy         ! surface boundary condition: diurnal cycle
    USE sbcwave , ONLY :   cdn_wave ! wave module
-   USE lib_fortran    ! to use key_nosignedzero
+   USE lib_fortran    ! to use key_nosignedzero and glob_sum
    !
 #if defined key_si3
@@ -348,4 +348,7 @@
       !                                      !- fill the bulk structure with namelist informations
       CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
+      sf(jp_wndi )%zsgn = -1._wp   ;   sf(jp_wndj )%zsgn = -1._wp   ! vector field at T point: overwrite default definition of zsgn
+      sf(jp_uoatm)%zsgn = -1._wp   ;   sf(jp_voatm)%zsgn = -1._wp   ! vector field at T point: overwrite default definition of zsgn
+      sf(jp_hpgi )%zsgn = -1._wp   ;   sf(jp_hpgj )%zsgn = -1._wp   ! vector field at T point: overwrite default definition of zsgn
       !
       DO jfpr= 1, jpfld
@@ -509,32 +512,29 @@
       ! Sanity/consistence test on humidity at first time step to detect potential screw-up:
       IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*) ''
-#if defined key_agrif
-         IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ==='
-#else
-         ztst = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain
-         IF( ztst > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points!
-            ztst = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztst ! mean humidity over ocean on proc
-            llerr = .FALSE.
-            SELECT CASE( nhumi )
-            CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!)
-               IF( (ztst <   0._wp) .OR. (ztst > 0.065_wp) )   llerr = .TRUE.
-            CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K]
-               IF( (ztst < 110._wp) .OR. (ztst >  320._wp) )   llerr = .TRUE.
-            CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%]
-               IF( (ztst <   0._wp) .OR. (ztst >  100._wp) )   llerr = .TRUE.
-            END SELECT
-            IF(llerr) THEN
-               WRITE(ctmp1,'("   Error on mean humidity value: ",f10.5)') ztst
-               CALL ctl_stop( 'STOP', ctmp1, 'Something is wrong with air humidity!!!', &
-                  &   ' ==> check the unit in your input files'       , &
-                  &   ' ==> check consistence of namelist choice: specific? relative? dew-point?', &
-                  &   ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' )
-            END IF
-         END IF
-         IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ==='
-#endif
-         IF(lwp) WRITE(numout,*) ''
-      END IF !IF( kt == nit000 )
+         ! mean humidity over ocean on proc
+         ztst = glob_sum( 'sbcblk', sf(jp_humi)%fnow(:,:,1) * e1e2t(:,:) * tmask(:,:,1) ) / glob_sum( 'sbcblk', e1e2t(:,:) * tmask(:,:,1) )
+         llerr = .FALSE.
+         SELECT CASE( nhumi )
+         CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!)
+            IF( (ztst <   0._wp) .OR. (ztst > 0.065_wp) )   llerr = .TRUE.
+         CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K]
+            IF( (ztst < 110._wp) .OR. (ztst >  320._wp) )   llerr = .TRUE.
+         CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%]
+            IF( (ztst <   0._wp) .OR. (ztst >  100._wp) )   llerr = .TRUE.
+         END SELECT
+         IF(llerr) THEN
+            WRITE(ctmp1,'("   Error on mean humidity value: ",f10.5)') ztst
+            CALL ctl_stop( 'STOP', ctmp1, 'Something is wrong with air humidity!!!', &
+               &   ' ==> check the unit in your input files'       , &
+               &   ' ==> check consistence of namelist choice: specific? relative? dew-point?', &
+               &   ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' )
+         ENDIF
+         IF(lwp) THEN
+            WRITE(numout,*) ''
+            WRITE(numout,*) ' Global mean humidity at kt = nit000: ', ztst
+            WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ==='
+            WRITE(numout,*) ''
+         ENDIF
+      ENDIF   !IF( kt == nit000 )
       !                                            ! compute the surface ocean fluxes using bulk formulea
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
@@ -622,6 +622,6 @@
       !!---------------------------------------------------------------------
       INTEGER , INTENT(in   )                 ::   kt     ! time step index
-      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndi  ! atmospheric wind at U-point              [m/s]
-      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndj  ! atmospheric wind at V-point              [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndi  ! atmospheric wind at T-point              [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndj  ! atmospheric wind at T-point              [m/s]
       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqair  ! specific humidity at T-points            [kg/kg]
       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   ptair  ! potential temperature at T-points        [Kelvin]
@@ -832,7 +832,7 @@
 
          IF( ln_crt_fbk ) THEN
-            CALL lbc_lnk( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. )
+            CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp, taum, 'T', 1._wp )
          ELSE
-            CALL lbc_lnk( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
+            CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp )
          ENDIF
 
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcflx.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcflx.F90	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcflx.F90	(revision 14403)
@@ -119,4 +119,6 @@
          !                                         ! fill sf with slf_i and control print
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
+         sf(jp_utau)%cltype = 'U'   ;   sf(jp_utau)%zsgn = -1._wp   ! vector field at U point: overwrite default definition of cltype and zsgn
+         sf(jp_vtau)%cltype = 'V'   ;   sf(jp_vtau)%zsgn = -1._wp   ! vector field at V point: overwrite default definition of cltype and zsgn
          !
       ENDIF
@@ -129,9 +131,9 @@
             qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
          ELSE
-            DO_2D( 0, 0, 0, 0 )
-               qsr(ji,jj) =          sf(jp_qsr)%fnow(ji,jj,1)   * tmask(ji,jj,1)
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+               qsr(ji,jj) =     sf(jp_qsr)%fnow(ji,jj,1) * tmask(ji,jj,1)
             END_2D
          ENDIF
-         DO_2D( 0, 0, 0, 0 )                                      ! set the ocean fluxes from read fields
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                  ! set the ocean fluxes from read fields
             utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
             vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
@@ -143,8 +145,4 @@
          !!clem: I do not think it is needed
          !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
-         !
-         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x)
-         CALL lbc_lnk( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
-            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  )
          !
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
Index: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcwave.F90
===================================================================
--- NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcwave.F90	(revision 14402)
+++ NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/SBC/sbcwave.F90	(revision 14403)
@@ -503,4 +503,5 @@
                !
                CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
+               sf_sd(jp_usd)%zsgn = -1._wp   ;  sf_sd(jp_vsd)%zsgn = -1._wp   ! vector field at T point: overwrite default definition of zsgn
             ENDIF
             !
