Changeset 11123 for NEMO/trunk

NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.sty

-                      r11043
+                      r11123
 %% LaTeX packages
 %% ==============================================================================
 \usepackage{natbib}           %% bib
 \usepackage{caption}          %% caption
 \usepackage{xcolor}           %% color
 \usepackage{times}            %% font
 \usepackage{hyperref}         %% hyper
 \usepackage{idxlayout}        %% index
 \usepackage{enumitem}         %% list
 \usepackage{minted}           %% listing
 \usepackage{amsmath}          %% maths
 \usepackage{fancyhdr}         %% page
 \usepackage{minitoc}          %% toc
 \usepackage{subfiles}         %% subdocs
 \usepackage[utf8]{inputenc}   %% input encoding
 \usepackage{draftwatermark}   %% watermark
 \usepackage{textcomp}         %% Companion fonts
+\usepackage{natbib}                      %% bib
+\usepackage{caption}                     %% caption
+\usepackage{xcolor}                      %% color
+\usepackage{times}                       %% font
+\usepackage{hyperref}                    %% hyper
+\usepackage{idxlayout}                   %% index
+\usepackage{enumitem}                    %% list
+\usepackage[outputdir=../build]{minted}   %% listing
+\usepackage{amsmath}                     %% maths
+\usepackage{fancyhdr}                    %% page
+\usepackage{minitoc}                     %% toc
+\usepackage{subfiles}                    %% subdocs
+\usepackage[utf8]{inputenc}              %% input encoding
+\usepackage{draftwatermark}              %% watermark
+\usepackage{textcomp}                    %% Companion fonts
 %% Extensions in bundle package
 …
 \newcommand{\pd}[2][]{\ensuremath{\frac{\partial #1}{\partial #2}}}
-%% Shortened DOI in bibliography
-\newcommand{\doi}[1]{\href{http://dx.doi.org/#1}{doi:#1}}
 %% Namelists inclusion
 \newcommand{\nlst}[1]{\forfile{../../../namelists/#1}}

NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.tex

-                      r11043
+                      r11123
 %% Chapters
 \subfile{../subfiles/chap_model_basics}
 \subfile{../subfiles/chap_time_domain}   % Time discretisation (time stepping strategy)
 \subfile{../subfiles/chap_DOM}           % Space discretisation
 \subfile{../subfiles/chap_TRA}           % Tracer advection/diffusion equation
 \subfile{../subfiles/chap_DYN}           % Dynamics : momentum equation
 \subfile{../subfiles/chap_SBC}           % Surface Boundary Conditions
 \subfile{../subfiles/chap_LBC}           % Lateral Boundary Conditions
 \subfile{../subfiles/chap_LDF}           % Lateral diffusion
 \subfile{../subfiles/chap_ZDF}           % Vertical diffusion
 \subfile{../subfiles/chap_DIA}           % Outputs and Diagnostics
 \subfile{../subfiles/chap_OBS}           % Observation operator
 \subfile{../subfiles/chap_ASM}           % Assimilation increments
 \subfile{../subfiles/chap_STO}           % Stochastic param.
 \subfile{../subfiles/chap_misc}          % Miscellaneous topics
 \subfile{../subfiles/chap_CONFIG}        % Predefined configurations
+\subfile{../subfiles/chap_time_domain}    % Time discretisation (time stepping strategy)
+\subfile{../subfiles/chap_DOM}            % Space discretisation
+\subfile{../subfiles/chap_TRA}            % Tracer advection/diffusion equation
+\subfile{../subfiles/chap_DYN}            % Dynamics : momentum equation
+\subfile{../subfiles/chap_SBC}            % Surface Boundary Conditions
+\subfile{../subfiles/chap_LBC}            % Lateral Boundary Conditions
+\subfile{../subfiles/chap_LDF}            % Lateral diffusion
+\subfile{../subfiles/chap_ZDF}            % Vertical diffusion
+\subfile{../subfiles/chap_DIA}            % Outputs and Diagnostics
+\subfile{../subfiles/chap_OBS}            % Observation operator
+\subfile{../subfiles/chap_ASM}            % Assimilation increments
+\subfile{../subfiles/chap_STO}            % Stochastic param.
+\subfile{../subfiles/chap_misc}           % Miscellaneous topics
+\subfile{../subfiles/chap_CONFIG}         % Predefined configurations
 %% Appendix
 \appendix
 \subfile{../subfiles/annex_A}             % Generalised vertical coordinate
 \subfile{../subfiles/annex_B}             % Diffusive operator
 \subfile{../subfiles/annex_C}             % Discrete invariants of the eqs.
 \subfile{../subfiles/annex_iso}           % Isoneutral diffusion using triads
 \subfile{../subfiles/annex_D}             % Coding rules
+\subfile{../subfiles/annex_A}             % Generalised vertical coordinate
+\subfile{../subfiles/annex_B}             % Diffusive operator
+\subfile{../subfiles/annex_C}             % Discrete invariants of the eqs.
+\subfile{../subfiles/annex_iso}            % Isoneutral diffusion using triads
+\subfile{../subfiles/annex_D}             % Coding rules
 %% Not included
+%\subfile{../subfiles/chap_conservation}  %
+%\subfile{../subfiles/chap_model_basics_zstar}
+%\subfile{../subfiles/chap_DIU}
+%\subfile{../subfiles/chap_conservation}
 %\subfile{../subfiles/annex_E}            % Notes on some on going staff
 %% Backmatter

NEMO/trunk/doc/latex/NEMO/subfiles

Property svn:ignore deleted

NEMO/trunk/doc/latex/NEMO/subfiles/annex_A.tex

r10442	r11123
399	399
400	400	As in $z$-coordinate,
401		the horizontal pressure gradient can be split in two parts following \citet{~~Marsaleix_al~~_OM08}.
	401	the horizontal pressure gradient can be split in two parts following \citet{marsaleix.auclair.ea_OM08}.
402	402	Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$,
403	403	and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$.

NEMO/trunk/doc/latex/NEMO/subfiles/annex_B.tex

-                      r10442
+                      r11123
 the ($i$,$j$,$k$) curvilinear coordinate system in which
 the equations of the ocean circulation model are formulated,
 takes the following form \citep{Redi_JPO82}:
+takes the following form \citep{redi_JPO82}:
 \begin{equation}
 …
 In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean,
 so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}:
+so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}:
 \begin{subequations}
   \label{apdx:B4}

NEMO/trunk/doc/latex/NEMO/subfiles/annex_D.tex

r10442	r11123
32	32
33	33	To satisfy part of these aims, \NEMO is written with a coding standard which is close to the ECMWF rules,
34		named DOCTOR \citep{~~Gibson_TR~~86}.
	34	named DOCTOR \citep{gibson_rpt86}.
35	35	These rules present some advantages like:
36	36

NEMO/trunk/doc/latex/NEMO/subfiles/annex_E.tex

-                      r10442
+                      r11123
 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error
 \citep{Shchepetkin_McWilliams_OM05}.
 The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}.
+\citep{shchepetkin.mcwilliams_OM05}.
+The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}.
 It is a relatively good compromise between accuracy and smoothness.
 It is not a \emph{positive} scheme meaning false extrema are permitted but
 …
 the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity
 (forward in time).
 This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.
+This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
 UBS and QUICK schemes only differ by one coefficient.
 Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
+Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme.
 …
 $\tau_w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme,
 or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following
 \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS.
+\citet{shchepetkin.mcwilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS.
 The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme.
 …
 \subsection{Griffies iso-neutral diffusion operator}
 Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} scheme,
+Let try to define a scheme that get its inspiration from the \citet{griffies.gnanadesikan.ea_JPO98} scheme,
 but is formulated within the \NEMO framework
 (\ie using scale factors rather than grid-size and having a position of $T$-points that
 …
 Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers
 (\ie they enter the computation of density), but it does not work for a passive tracer.
 \citep{Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that
+\citep{griffies.gnanadesikan.ea_JPO98} introduce a different way to discretise the off-diagonal terms that
 nicely solve the problem.
 The idea is to get rid of combinations of an averaged in one direction combined with
 …
 \]
 \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
+\citep{griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
 It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows:
 …
 The horizontal component reduces to the one use for an horizontal laplacian operator and
 the vertical one keeps the same complexity, but not more.
 This property has been used to reduce the computational time \citep{Griffies_JPO98},
+This property has been used to reduce the computational time \citep{griffies_JPO98},
 but it is not of practical use as usually $A \neq A_e$.
 Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which

NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex

-                      r10442
+                      r11123
   the vertical skew flux is further reduced to ensure no vertical buoyancy flux,
   giving an almost pure horizontal diffusive tracer flux within the mixed layer.
   This is similar to the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper}
+  This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:Gerdes-taper}
 \item[\np{ln\_botmix\_triad}]
   See \autoref{sec:iso_bdry}.
 …
 \label{sec:iso}
 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98},
+We have implemented into \NEMO a scheme inspired by \citet{griffies.gnanadesikan.ea_JPO98},
 but formulated within the \NEMO framework, using scale factors rather than grid-sizes.
 …
 \subsection{Expression of the skew-flux in terms of triad slopes}
 \citep{Griffies_al_JPO98} introduce a different discretization of the off-diagonal terms that
+\citep{griffies.gnanadesikan.ea_JPO98} introduce a different discretization of the off-diagonal terms that
 nicely solves the problem.
 % Instead of multiplying the mean slope calculated at the $u$-point by
 …
 To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$.
 \citet{Griffies_al_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells,
+\citet{griffies.gnanadesikan.ea_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells,
 defined in terms of the distances between $T$, $u$,$f$ and $w$-points.
 This is the natural discretization of \autoref{eq:cts-var}.
 …
 As discussed in \autoref{subsec:LDF_slp_iso},
 iso-neutral slopes relative to geopotentials must be bounded everywhere,
 both for consistency with the small-slope approximation and for numerical stability \citep{Cox1987, Griffies_Bk04}.
+both for consistency with the small-slope approximation and for numerical stability \citep{cox_OM87, griffies_bk04}.
 The bound chosen in \NEMO is applied to each component of the slope separately and
 has a value of $1/100$ in the ocean interior.
 …
 \footnote{
   To ensure good behaviour where horizontal density gradients are weak,
   we in fact follow \citet{Gerdes1991} and
+  we in fact follow \citet{gerdes.koberle.ea_CD91} and
   set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.
+}
 …
 This approach is similar to multiplying the iso-neutral diffusion coefficient by
 $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep slopes,
 as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}).
+as suggested by \citet{gerdes.koberle.ea_CD91} (see also \citet{griffies_bk04}).
 Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$
 …
 However, when \np{ln\_traldf\_triad} is set true,
 \NEMO instead implements eddy induced advection according to the so-called skew form \citep{Griffies_JPO98}.
+\NEMO instead implements eddy induced advection according to the so-called skew form \citep{griffies_JPO98}.
 It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
 For example in the (\textbf{i},\textbf{k}) plane,
 …
 it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer
 \autoref{eq:eiv_v}.
 This ensures that the eiv velocities do not restratify the mixed layer \citep{Treguier1997,Danabasoglu_al_2008}.
+This ensures that the eiv velocities do not restratify the mixed layer \citep{treguier.held.ea_JPO97,danabasoglu.ferrari.ea_JC08}.
 Equivantly, in terms of the skew-flux formulation we use here,
 the linear slope tapering within the mixed-layer gives a linearly varying vertical flux,
 …
 $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$)
 points (see Table \autoref{tab:cell}) respectively.
 We follow \citep{Griffies_Bk04} and calculate the streamfunction at a given $uw$-point from
+We follow \citep{griffies_bk04} and calculate the streamfunction at a given $uw$-point from
 the surrounding four triads according to:
 \[

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex

-                      r10442
+                      r11123
 it may be preferable to introduce the increment gradually into the ocean model in order to
 minimize spurious adjustment processes.
 This technique is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}.
+This technique is referred to as Incremental Analysis Updates (IAU) \citep{bloom.takacs.ea_MWR96}.
 IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI.
 IAU is used when \np{ln\_asmiau} is set to true.
 …
 This type of the initialisation reduces the vertical velocity magnitude and
 alleviates the problem of the excessive unphysical vertical mixing in the first steps of the model integration
 \citep{Talagrand_JAS72, Dobricic_al_OS07}.
+\citep{talagrand_JAS72, dobricic.pinardi.ea_OS07}.
 Diffusion coefficients are defined as $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$.
 The divergence damping is activated by assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist

NEMO/trunk/doc/latex/NEMO/subfiles/chap_CONFIG.tex

-                      r11112
+                      r11123
       The two "north pole" are the foci of a series of embedded ellipses (blue curves) which
       are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes).
       Then, following \citet{Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed which
+      Then, following \citet{madec.imbard_CD96}, the normal to the series of ellipses (red curves) is computed which
       provides the j-lines of the mesh (pseudo longitudes).
+    }
 …
 \label{subsec:CFG_orca_grid}
 The ORCA grid is a tripolar grid based on the semi-analytical method of \citet{Madec_Imbard_CD96}.
+The ORCA grid is a tripolar grid based on the semi-analytical method of \citet{madec.imbard_CD96}.
 It allows to construct a global orthogonal curvilinear ocean mesh which has no singularity point inside
 the computational domain since two north mesh poles are introduced and placed on lands.
 …
 %climate change (Marti et al., 2009).
 %It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs)
 %documented in \citet{Griffies_al_OM09}.
+%documented in \citet{griffies.biastoch.ea_OM09}.
 This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m) in the upper 150m
 (see \autoref{tab:orca_zgr} and \autoref{fig:zgr}).
 The bottom topography and the coastlines are derived from the global atlas of Smith and Sandwell (1997).
 The default forcing uses the boundary forcing from \citet{Large_Yeager_Rep04} (see \autoref{subsec:SBC_blk_core}),
+The default forcing uses the boundary forcing from \citet{large.yeager_rpt04} (see \autoref{subsec:SBC_blk_core}),
 which was developed for the purpose of running global coupled ocean-ice simulations without
 an interactive atmosphere.
 This \citet{Large_Yeager_Rep04} dataset is available through
+This \citet{large.yeager_rpt04} dataset is available through
 the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}.
 The "normal year" of \citet{Large_Yeager_Rep04} has been chosen of the NEMO distribution since release v3.3.
 …
 \label{sec:CFG_gyre}
 The GYRE configuration \citep{Levy_al_OM10} has been built to
+The GYRE configuration \citep{levy.klein.ea_OM10} has been built to
 simulate the seasonal cycle of a double-gyre box model.
 It consists in an idealized domain similar to that used in the studies of \citet{Drijfhout_JPO94} and
 \citet{Hazeleger_Drijfhout_JPO98, Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00},
+It consists in an idealized domain similar to that used in the studies of \citet{drijfhout_JPO94} and
+\citet{hazeleger.drijfhout_JPO98, hazeleger.drijfhout_JPO99, hazeleger.drijfhout_JGR00, hazeleger.drijfhout_JPO00},
 over which an analytical seasonal forcing is applied.
 This allows to investigate the spontaneous generation of a large number of interacting, transient mesoscale eddies
 …
 The configuration is meant to represent an idealized North Atlantic or North Pacific basin.
 The circulation is forced by analytical profiles of wind and buoyancy fluxes.
 The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{Levy_al_OM10}.
+The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{levy.klein.ea_OM10}.
 The wind stress is zonal and its curl changes sign at 22\deg{N} and 36\deg{N}.
 It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain and
 …
       \protect\label{fig:GYRE}
       Snapshot of relative vorticity at the surface of the model domain in GYRE R9, R27 and R54.
       From \citet{Levy_al_OM10}.
+      From \citet{levy.klein.ea_OM10}.
+    }
   \end{center}

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex

-                      r10509
+                      r11123
 remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project.
 Options are defined by \ngn{namflo} namelis variables.
 The algorithm used is based either on the work of \cite{Blanke_Raynaud_JPO97} (default option),
+The algorithm used is based either on the work of \cite{blanke.raynaud_JPO97} (default option),
 or on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}).
 Note that the \cite{Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which
+Note that the \cite{blanke.raynaud_JPO97} algorithm have the advantage of providing trajectories which
 are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry.
 …
 The steric effect is therefore not explicitely represented.
 This approximation does not represent a serious error with respect to the flow field calculated by the model
 \citep{Greatbatch_JGR94}, but extra attention is required when investigating sea level,
+\citep{greatbatch_JGR94}, but extra attention is required when investigating sea level,
 as steric changes are an important contribution to local changes in sea level on seasonal and climatic time scales.
 This is especially true for investigation into sea level rise due to global warming.
 Fortunately, the steric contribution to the sea level consists of a spatially uniform component that
 can be diagnosed by considering the mass budget of the world ocean \citep{Greatbatch_JGR94}.
+can be diagnosed by considering the mass budget of the world ocean \citep{greatbatch_JGR94}.
 In order to better understand how global mean sea level evolves and thus how the steric sea level can be diagnosed,
 we compare, in the following, the non-Boussinesq and Boussinesq cases.
 …
 the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid.
 Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by
+Nevertheless, following \citep{greatbatch_JGR94}, the steric effect on the volume can be diagnosed by
 considering the mass budget of the ocean.
 The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface mass flux
 must be compensated by a spatially uniform change in the mean sea level due to expansion/contraction of the ocean
 \citep{Greatbatch_JGR94}.
+\citep{greatbatch_JGR94}.
 In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$,
 the total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level,
 …
 This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since,
 with the exception of only a small percentage of the ocean, density in the World Ocean varies by no more than
 $\%$ from this value (\cite{Gill1982}, page 47).
+$\%$ from this value (\cite{gill_bk82}, page 47).
 Second, we have assumed here that the total ocean surface, $\mathcal{A}$,
 …
 so that there are no associated ocean currents.
 Hence, the dynamically relevant sea level is the effective sea level,
 \ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{Campin_al_OM08}.
+\ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{campin.marshall.ea_OM08}.
 However, in the current version of \NEMO the sea-ice is levitating above the ocean without mass exchanges between
 ice and ocean.
 …
 Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key:
 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
+- the mixed layer depth (based on a density criterion \citep{de-boyer-montegut.madec.ea_JGR04}) (\mdl{diahth})
 - the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex

-                      r10442
+                      r11123
 %===============================================================
 The warm layer is calculated using the model of \citet{Takaya_al_JGR10} (TAKAYA10 model hereafter).
+The warm layer is calculated using the model of \citet{takaya.bidlot.ea_JGR10} (TAKAYA10 model hereafter).
 This is a simple flux based model that is defined by the equations
 \begin{align}
 …
 where $Q_{\rm{h}}$ is the sensible and latent heat flux, $Q_{\rm{lw}}$ is the long wave flux,
 and $Q_{\rm{sol}}$ is the solar flux absorbed within the diurnal warm layer.
 For $Q_{\rm{sol}}$ the 9 term representation of \citet{Gentemann_al_JGR09} is used.
+For $Q_{\rm{sol}}$ the 9 term representation of \citet{gentemann.minnett.ea_JGR09} is used.
 In equation \autoref{eq:ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$,
 where $L_a=0.3$\footnote{
 …
 %===============================================================
 The cool skin is modelled using the framework of \citet{Saunders_JAS82} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$.
+The cool skin is modelled using the framework of \citet{saunders_JAS67} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$.
 As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\rm{cs}}$ becomes
 \[
 …
 \end{equation}
 where $\mu$ is the kinematic viscosity of sea water and $\lambda$ is a constant of proportionality which
 \citet{Saunders_JAS82} suggested varied between 5 and 10.
+\citet{saunders_JAS67} suggested varied between 5 and 10.
 The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{Artale_al_JGR02},
 which is shown in \citet{Tu_Tsuang_GRL05} to outperform a number of other parametrisations at
+The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02},
+which is shown in \citet{tu.tsuang_GRL05} to outperform a number of other parametrisations at
 both low and high wind speeds.
 Specifically,

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex

-                      r10502
+                      r11123
 the centre of each face of the cells (\autoref{fig:cell}).
 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification
 \citep{Mesinger_Arakawa_Bk76}.
+\citep{mesinger.arakawa_bk76}.
 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and
 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points.
 …
 (\ie as the analytical first derivative of the transformation that
 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$)
 is specific to the \NEMO model \citep{Marti_al_JGR92}.
+is specific to the \NEMO model \citep{marti.madec.ea_JGR92}.
 As an example, $e_{1t}$ is defined locally at a $t$-point,
 whereas many other models on a C grid choose to define such a scale factor as
 …
 since they are first introduced in the continuous equations;
 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids
 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{Treguier1996}.
+(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}.
 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true
 (\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}).
 This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
+This uses a depth independent $\tanh$ function for the stretching \citep{madec.delecluse.ea_JPO96}:
 \[
 …
 A stretching function,
 modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}),
+modified from the commonly used \citet{song.haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}),
 is also available and is more commonly used for shelf seas modelling:
 …
 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in
 an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}.
+an analytical terrain-following stretching \citet{siddorn.furner_OM13}.
 In this case the a stretching function $\gamma$ is defined such that:
 …
   \includegraphics[]{Fig_DOM_compare_coordinates_surface}
   \caption{
     A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines),
+    A comparison of the \citet{song.haidvogel_JCP94} $S$-coordinate (solid lines),
     a 50 level $Z$-coordinate (contoured surfaces) and
     the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for
+    the \citet{siddorn.furner_OM13} $S$-coordinate (dashed lines) in the surface $100~m$ for
     a idealised bathymetry that goes from $50~m$ to $5500~m$ depth.
     For clarity every third coordinate surface is shown.
 …
 creating a non-analytical vertical coordinate that
 therefore may suffer from large gradients in the vertical resolutions.
 This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching,
+This stretching is less straightforward to implement than the \citet{song.haidvogel_JCP94} stretching,
 but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes.
 As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than
+As with the \citet{song.haidvogel_JCP94} stretching the stretch is only applied at depths greater than
 the critical depth $h_c$.
 In this example two options are available in depths shallower than $h_c$,
 …
 Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as
 large slopes lead to hydrostatic consistency.
 A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented,
+A hydrostatic consistency parameter diagnostic following \citet{haney_JPO91} has been implemented,
 and is output as part of the model mesh file at the start of the run.
 …
 Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with
 steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}.
+steps that follow the face of the model cells (step like topography) \citep{madec.delecluse.ea_JPO96}.
 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which
 gives the number of ocean levels (\ie those that are not masked) at each $t$-point.

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex

-                      r10499
+                      r11123
 Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to
 the sea surface height equation otherwise tracer content will not be conserved
 \citep{Griffies_al_MWR01, Leclair_Madec_OM09}.
+\citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}.
 The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom,
 …
 Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....}
 A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}.
+A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}.
 The idea is to get rid of the double averaging by considering triad combinations of vorticity.
 It is noteworthy that this solution is conceptually quite similar to the one proposed by
 \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).
 The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified
 for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.
+\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).
+The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified
+for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme.
 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
 \[
 …
 (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry)
 that tends to reinforce the topostrophy of the flow
 (\ie the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.
+(\ie the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}.
 Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as
 …
 (\ie $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).
 Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of
 the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}.
+the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}.
 Furthermore, used in combination with a partial steps representation of bottom topography,
 it improves the interaction between current and topography,
 leading to a larger topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}.
+leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}.
 %--------------------------------------------------------------------------------------------------------------
 …
 When \np{ln\_dynzad\_zts}\forcode{ = .true.},
 a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term.
 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.
+This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.
 Note that in this case,
 a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability,
 …
 a $2^{nd}$ order centered finite difference scheme, CEN2,
 or a $3^{rd}$ order upstream biased scheme, UBS.
 The latter is described in \citet{Shchepetkin_McWilliams_OM05}.
+The latter is described in \citet{shchepetkin.mcwilliams_OM05}.
 The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}.
 In flux form, the schemes differ by the choice of a space and time interpolation to define the value of
 …
 where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$.
 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error
 \citep{Shchepetkin_McWilliams_OM05}.
 The overall performance of the advection scheme is similar to that reported in \citet{Farrow1995}.
+\citep{shchepetkin.mcwilliams_OM05}.
+The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}.
 It is a relatively good compromise between accuracy and smoothness.
 It is not a \emph{positive} scheme, meaning that false extrema are permitted.
 …
 while the second term, which is the diffusion part of the scheme,
 is evaluated using the \textit{before} velocity (forward in time).
 This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.
+This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
 Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by
 one coefficient.
 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
+Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme.
 …
 Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers
 (\eg, \citet{Song1998, Shchepetkin_McWilliams_OM05}).
+(\eg, \citet{song_MWR98, shchepetkin.mcwilliams_OM05}).
 A number of different pressure gradient options are coded but the ROMS-like,
 density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation.
 $\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.})
+$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.})
 \begin{equation}
   \label{eq:dynhpg_sco}
 …
 $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ = .true.})
 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05}
+$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05}
 (\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development)
 Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated
 (\key{vvl}) because in that case, even with a flat bottom,
 the coordinate surfaces are not horizontal but follow the free surface \citep{Levier2007}.
+the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}.
 The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) is available as
 an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \key{vvl} is active.
 …
 corresponds to the water replaced by the ice shelf.
 This top pressure is constant over time.
 A detailed description of this method is described in \citet{Losch2008}.\\
+A detailed description of this method is described in \citet{losch_JGR08}.\\
 The pressure gradient due to ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in
 …
 the physical phenomenon that controls the time-step is internal gravity waves (IGWs).
 A semi-implicit scheme for doubling the stability limit associated with IGWs can be used
 \citep{Brown_Campana_MWR78, Maltrud1998}.
+\citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}.
 It involves the evaluation of the hydrostatic pressure gradient as
 an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$
 …
 which imposes a very small time step when an explicit time stepping is used.
 Two methods are proposed to allow a longer time step for the three-dimensional equations:
 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt}),
+the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt?}),
 and the split-explicit free surface described below.
 The extra term introduced in the filtered method is calculated implicitly,
 …
 The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined),
 also called the time-splitting formulation, follows the one proposed by \citet{Shchepetkin_McWilliams_OM05}.
+also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}.
 The general idea is to solve the free surface equation and the associated barotropic velocity equations with
 a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables
 …
 (see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration.
 Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm
 detailed in \citet{Shchepetkin_McWilliams_OM05}.
+detailed in \citet{shchepetkin.mcwilliams_OM05}.
 AB3-AM4 coefficients used in \NEMO follow the second-order accurate,
 "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08}
+"multi-purpose" stability compromise as defined in \citet{shchepetkin.mcwilliams_ibk09}
 (see their figure 12, lower left).
 …
 and time splitting not compatible.
 Advective barotropic velocities are obtained by using a secondary set of filtering weights,
 uniquely defined from the filter coefficients used for the time averaging (\citet{Shchepetkin_McWilliams_OM05}).
+uniquely defined from the filter coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}).
 Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to
 obtain exact conservation.
 …
 external gravity waves in idealized or weakly non-linear cases.
 Although the damping is lower than for the filtered free surface,
 it is still significant as shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave.
+it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave.
 %>>>>>===============
 …
 the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}.
 We have tried various forms of such filtering,
 with the following method discussed in \cite{Griffies_al_MWR01} chosen due to
+with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to
 its stability and reasonably good maintenance of tracer conservation properties (see ??).
 …
 \label{subsec:DYN_spg_fltp}
 The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation.
+The filtered formulation follows the \citet{roullet.madec_JGR00} implementation.
 The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly.
 The elliptic solvers available in the code are documented in \autoref{chap:MISC}.
 …
 There are two main options for wetting and drying code (wd):
 (a) an iterative limiter (il) and (b) a directional limiter (dl).
 The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for RO
+The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO
 MS
 which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative
+which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative
 limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$
 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated
 …
 \cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic
+\cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic
 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than
 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex

-                      r10614
+                      r11123
 The BDY module was modelled on the OBC module (see NEMO 3.4) and shares many features and
 a similar coding structure \citep{Chanut2005}.
+a similar coding structure \citep{chanut_rpt05}.
 The specification of the location of the open boundary is completely flexible and
 allows for example the open boundary to follow an isobath or other irregular contour.
 …
 \label{subsec:BDY_FRS_scheme}
 The Flow Relaxation Scheme (FRS) \citep{Davies_QJRMS76,Engerdahl_Tel95},
+The Flow Relaxation Scheme (FRS) \citep{davies_QJRMS76,engedahl_T95},
 applies a simple relaxation of the model fields to externally-specified values over
 a zone next to the edge of the model domain.
 …
 \label{subsec:BDY_flather_scheme}
 The \citet{Flather_JPO94} scheme is a radiation condition on the normal,
+The \citet{flather_JPO94} scheme is a radiation condition on the normal,
 depth-mean transport across the open boundary.
 It takes the form
 …
 \label{subsec:BDY_orlanski_scheme}
 The Orlanski scheme is based on the algorithm described by \citep{Marchesiello2001}, hereafter MMS.
+The Orlanski scheme is based on the algorithm described by \citep{marchesiello.mcwilliams.ea_OM01}, hereafter MMS.
 The adaptive Orlanski condition solves a wave plus relaxation equation at the boundary:

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex

-                      r10442
+                      r11123
 \gmcomment{
   we should emphasize here that the implementation is a rather old one.
   Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme.
+  Better work can be achieved by using \citet{griffies.gnanadesikan.ea_JPO98, griffies_bk04} iso-neutral scheme.
+}
 …
 %In practice, \autoref{eq:ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.
 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
+%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
 %In the $z$-coordinate, the derivative of the  \autoref{eq:ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.
 …
   thus the $in situ$ density can be used.
   This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$,
   where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{McDougall1987}
+  where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{mcdougall_JPO87}
   (see \autoref{subsec:TRA_bn2}).
 …
   Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for
   the constraint on iso-neutral fluxes.
   Following \citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral diffusive flux of
+  Following \citet{griffies_bk04}, instead of specifying directly that there is a zero neutral diffusive flux of
   locally referenced potential density, we stay in the $T$-$S$ plane and consider the balance between
   the neutral direction diffusive fluxes of potential temperature and salinity:
 …
 a minimum background horizontal diffusion for numerical stability reasons.
 To overcome this problem, several techniques have been proposed in which the numerical schemes of
 the ocean model are modified \citep{Weaver_Eby_JPO97, Griffies_al_JPO98}.
+the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}.
 Griffies's scheme is now available in \NEMO if \np{traldf\_grif\_iso} is set true; see Appdx \autoref{apdx:triad}.
 Here, another strategy is presented \citep{Lazar_PhD97}:
+Here, another strategy is presented \citep{lazar_phd97}:
 a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of
 grid point noise generated by the iso-neutral diffusion operator (\autoref{fig:LDF_ZDF1}).
 …
 Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,
 contrary to the \citet{Griffies_al_JPO98} operator which has that property.
+contrary to the \citet{griffies.gnanadesikan.ea_JPO98} operator which has that property.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 % In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04},
+% In addition and also for numerical stability reasons \citep{cox_OM87, griffies_bk04},
 % the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly
 % to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the
 % surface motivates this flattening of isopycnals near the surface).
 For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also be bounded by
+For numerical stability reasons \citep{cox_OM87, griffies_bk04}, the slopes must also be bounded by
 $1/100$ everywhere.
 This constraint is applied in a piecewise linear fashion, increasing from zero at the surface to
 …
 This variation is intended to reflect the lesser need for subgrid scale eddy mixing where
 the grid size is smaller in the domain.
 It was introduced in the context of the DYNAMO modelling project \citep{Willebrand_al_PO01}.
+It was introduced in the context of the DYNAMO modelling project \citep{willebrand.barnier.ea_PO01}.
 Note that such a grid scale dependance of mixing coefficients significantly increase the range of stability of
 model configurations presenting large changes in grid pacing such as global ocean models.
 …
 For example, in the ORCA2 global ocean model (see Configurations),
 the laplacian viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and
 decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}.
+decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s at the equator \citep{madec.delecluse.ea_JPO96, delecluse.madec_icol99}.
 This modification can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}.
 Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of
 …
 since it allows us to take advantage of all the advection schemes offered for the tracers
 (see \autoref{sec:TRA_adv}) and not just the $2^{nd}$ order advection scheme as in
 previous releases of OPA \citep{Madec1998}.
+previous releases of OPA \citep{madec.delecluse.ea_NPM98}.
 This is particularly useful for passive tracers where \emph{positivity} of the advection scheme is of
 paramount importance.

NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex

-                      r10442
+                      r11123
    and $M$ corresponds to $B$, $C$ or $D$.
    A more stable form of the great-circle distance formula for small distances ($x$ near 1)
    involves the arcsine function (\eg see p.~101 of \citet{Daley_Barker_Bk01}:
+   involves the arcsine function (\eg see p.~101 of \citet{daley.barker_bk01}:
    \begin{align*}
      s\left( {\rm P}, {\rm M} \right) & \hspace{-2mm} = \hspace{-2mm} & \sin^{-1} \! \left\{ \sqrt{ 1 - x^2 } \right\}
 …
   An iterative scheme that involves first mapping a quadrilateral cell into
   a cell with coordinates (0,0), (1,0), (0,1) and (1,1).
   This method is based on the SCRIP interpolation package \citep{Jones_1998}.
+  This method is based on the \href{https://github.com/SCRIP-Project/SCRIP}{SCRIP interpolation package}.
 \end{enumerate}
 …
 where ${{\bf r}_{}}_{\rm PA}$, ${{\bf r}_{}}_{\rm PB}$, etc. correspond to
 the vectors between points P and A, P and B, etc..
 The method used is similar to the method used in the SCRIP interpolation package \citep{Jones_1998}.
+The method used is similar to the method used in the \href{https://github.com/SCRIP-Project/SCRIP}{SCRIP interpolation package}.
 In order to speed up the grid search, there is the possibility to construct a lookup table for a user specified resolution.

NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex

-                      r10614
+                      r11123
 The only tricky point is therefore to specify the date at which we need to do the interpolation and
 the date of the records read in the input files.
 Following \citet{Leclair_Madec_OM09}, the date of a time step is set at the middle of the time step.
+Following \citet{leclair.madec_OM09}, the date of a time step is set at the middle of the time step.
 For example, for an experiment starting at 0h00'00" with a one hour time-step,
 a time interpolation will be performed at the following time: 0h30'00", 1h30'00", 2h30'00", etc.
 …
 %-------------------------------------------------------------------------------------------------------------
 The CORE bulk formulae have been developed by \citet{Large_Yeager_Rep04}.
+The CORE bulk formulae have been developed by \citet{large.yeager_rpt04}.
 They have been designed to handle the CORE forcing, a mixture of NCEP reanalysis and satellite data.
 They use an inertial dissipative method to compute the turbulent transfer coefficients
 (momentum, sensible heat and evaporation) from the 10 metre wind speed, air temperature and specific humidity.
 This \citet{Large_Yeager_Rep04} dataset is available through
+This \citet{large.yeager_rpt04} dataset is available through
 the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}.
 Note that substituting ERA40 to NCEP reanalysis fields does not require changes in the bulk formulea themself.
 This is the so-called DRAKKAR Forcing Set (DFS) \citep{Brodeau_al_OM09}.
+This is the so-called DRAKKAR Forcing Set (DFS) \citep{brodeau.barnier.ea_OM10}.
 Options are defined through the  \ngn{namsbc\_core} namelist variables.
 …
 The CLIO bulk formulae were developed several years ago for the Louvain-la-neuve coupled ice-ocean model
 (CLIO, \cite{Goosse_al_JGR99}).
+(CLIO, \cite{goosse.deleersnijder.ea_JGR99}).
 They are simpler bulk formulae.
 They assume the stress to be known and compute the radiative fluxes from a climatological cloud cover.
 …
 The SAL term should in principle be computed online as it depends on
 the model tidal prediction itself (see \citet{Arbic2004} for a
+the model tidal prediction itself (see \citet{arbic.garner.ea_DSR04} for a
 discussion about the practical implementation of this term).
 Nevertheless, the complex calculations involved would make this
 …
 %coastal modelling and becomes more and more often open ocean and climate modelling
 %\footnote{At least a top cells thickness of 1~meter and a 3 hours forcing frequency are
 %required to properly represent the diurnal cycle \citep{Bernie_al_JC05}. see also \autoref{fig:SBC_dcy}.}.
+%required to properly represent the diurnal cycle \citep{bernie.woolnough.ea_JC05}. see also \autoref{fig:SBC_dcy}.}.
 …
 \footnote{
   At least a top cells thickness of 1~meter and a 3 hours forcing frequency are required to
   properly represent the diurnal cycle \citep{Bernie_al_JC05}.
+  properly represent the diurnal cycle \citep{bernie.woolnough.ea_JC05}.
   see also \autoref{fig:SBC_dcy}.}.
 …
 %--------------------------------------------------------------------------------------------------------
 The namelist variable in \ngn{namsbc}, \np{nn\_isf}, controls the ice shelf representation.
 Description and result of sensitivity test to \np{nn\_isf} are presented in \citet{Mathiot2017}.
+Description and result of sensitivity test to \np{nn\_isf} are presented in \citet{mathiot.jenkins.ea_GMD17}.
 The different options are illustrated in \autoref{fig:SBC_isf}.
 …
    \item[\np{nn\_isfblk}\forcode{ = 1}]:
      The melt rate is based on a balance between the upward ocean heat flux and
      the latent heat flux at the ice shelf base. A complete description is available in \citet{Hunter2006}.
+     the latent heat flux at the ice shelf base. A complete description is available in \citet{hunter_rpt06}.
    \item[\np{nn\_isfblk}\forcode{ = 2}]:
      The melt rate and the heat flux are based on a 3 equations formulation
      (a heat flux budget at the ice base, a salt flux budget at the ice base and a linearised freezing point temperature equation).
      A complete description is available in \citet{Jenkins1991}.
+     A complete description is available in \citet{jenkins_JGR91}.
    \end{description}
      Temperature and salinity used to compute the melt are the average temperature in the top boundary layer \citet{Losch2008}.
+     Temperature and salinity used to compute the melt are the average temperature in the top boundary layer \citet{losch_JGR08}.
      Its thickness is defined by \np{rn\_hisf\_tbl}.
      The fluxes and friction velocity are computed using the mean temperature, salinity and velocity in the the first \np{rn\_hisf\_tbl} m.
 …
 \]
      where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn\_hisf\_tbl} meters).
      See \citet{Jenkins2010} for all the details on this formulation. It is the recommended formulation for realistic application.
+     See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. It is the recommended formulation for realistic application.
    \item[\np{nn\_gammablk}\forcode{ = 2}]:
      The salt and heat exchange coefficients are velocity and stability dependent and defined as:
 …
      $\Gamma_{Turb}$ the contribution of the ocean stability and
      $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion.
      See \citet{Holland1999} for all the details on this formulation.
+     See \citet{holland.jenkins_JPO99} for all the details on this formulation.
      This formulation has not been extensively tested in NEMO (not recommended).
    \end{description}
  \item[\np{nn\_isf}\forcode{ = 2}]:
    The ice shelf cavity is not represented.
    The fwf and heat flux are computed using the \citet{Beckmann2003} parameterisation of isf melting.
+   The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting.
    The fluxes are distributed along the ice shelf edge between the depth of the average grounding line (GL)
    (\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front
 …
 %-------------------------------------------------------------------------------------------------------------
 Icebergs are modelled as lagrangian particles in NEMO \citep{Marsh_GMD2015}.
 Their physical behaviour is controlled by equations as described in \citet{Martin_Adcroft_OM10} ).
+Icebergs are modelled as lagrangian particles in NEMO \citep{marsh.ivchenko.ea_GMD15}.
+Their physical behaviour is controlled by equations as described in \citet{martin.adcroft_OM10} ).
 (Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO).
 Icebergs are initially spawned into one of ten classes which have specific mass and thickness as
 …
 Then using the routine \rou{turb\_ncar} and starting from the neutral drag coefficent provided,
 the drag coefficient is computed according to the stable/unstable conditions of the
 air-sea interface following \citet{Large_Yeager_Rep04}.
+air-sea interface following \citet{large.yeager_rpt04}.
 …
 \label{subsec:SBC_wave_sdw}
 The Stokes drift is a wave driven mechanism of mass and momentum transport \citep{Stokes_1847}.
+The Stokes drift is a wave driven mechanism of mass and momentum transport \citep{stokes_ibk09}.
 It is defined as the difference between the average velocity of a fluid parcel (Lagrangian velocity)
 and the current measured at a fixed point (Eulerian velocity).
 …
 \begin{description}
 \item[\np{nn\_sdrift} = 0]: exponential integral profile parameterization proposed by
 \citet{Breivik_al_JPO2014}:
+\citet{breivik.janssen.ea_JPO14}:
 \[
 …
 \item[\np{nn\_sdrift} = 1]: velocity profile based on the Phillips spectrum which is considered to be a
 reasonable estimate of the part of the spectrum most contributing to the Stokes drift velocity near the surface
 \citep{Breivik_al_OM2016}:
+\citep{breivik.bidlot.ea_OM16}:
 \[
 …
 The surface stress felt by the ocean is the atmospheric stress minus the net stress going
 into the waves \citep{Janssen_al_TM13}. Therefore, when waves are growing, momentum and energy is spent and is not
+into the waves \citep{janssen.breivik.ea_rpt13}. Therefore, when waves are growing, momentum and energy is spent and is not
 available for forcing the mean circulation, while in the opposite case of a decaying sea
 state more momentum is available for forcing the ocean.
 …
       the mean value of the analytical cycle (blue line) over a time step,
       not as the mid time step value of the analytically cycle (red square).
       From \citet{Bernie_al_CD07}.
+      From \citet{bernie.guilyardi.ea_CD07}.
+    }
   \end{center}
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \cite{Bernie_al_JC05} have shown that to capture 90$\%$ of the diurnal variability of SST requires a vertical resolution in upper ocean of 1~m or better and a temporal resolution of the surface fluxes of 3~h or less.
+\cite{bernie.woolnough.ea_JC05} have shown that to capture 90$\%$ of the diurnal variability of SST requires a vertical resolution in upper ocean of 1~m or better and a temporal resolution of the surface fluxes of 3~h or less.
 Unfortunately high frequency forcing fields are rare, not to say inexistent.
 Nevertheless, it is possible to obtain a reasonable diurnal cycle of the SST knowning only short wave flux (SWF) at
 high frequency \citep{Bernie_al_CD07}.
+high frequency \citep{bernie.guilyardi.ea_CD07}.
 Furthermore, only the knowledge of daily mean value of SWF is needed,
 as higher frequency variations can be reconstructed from them,
 assuming that the diurnal cycle of SWF is a scaling of the top of the atmosphere diurnal cycle of incident SWF.
 The \cite{Bernie_al_CD07} reconstruction algorithm is available in \NEMO by
+The \cite{bernie.guilyardi.ea_CD07} reconstruction algorithm is available in \NEMO by
 setting \np{ln\_dm2dc}\forcode{ = .true.} (a \textit{\ngn{namsbc}} namelist variable) when
 using CORE bulk formulea (\np{ln\_blk\_core}\forcode{ = .true.}) or
 the flux formulation (\np{ln\_flx}\forcode{ = .true.}).
 The reconstruction is performed in the \mdl{sbcdcy} module.
 The detail of the algoritm used can be found in the appendix~A of \cite{Bernie_al_CD07}.
+The detail of the algoritm used can be found in the appendix~A of \cite{bernie.guilyardi.ea_CD07}.
 The algorithm preserve the daily mean incoming SWF as the reconstructed SWF at
 a given time step is the mean value of the analytical cycle over this time step (\autoref{fig:SBC_diurnal}).
 …
 (observed, climatological or an atmospheric model product),
 \textit{SSS}$_{Obs}$ is a sea surface salinity
 (usually a time interpolation of the monthly mean Polar Hydrographic Climatology \citep{Steele2001}),
+(usually a time interpolation of the monthly mean Polar Hydrographic Climatology \citep{steele.morley.ea_JC01}),
 $\left.S\right|_{k=1}$ is the model surface layer salinity and
 $\gamma_s$ is a negative feedback coefficient which is provided as a namelist parameter.
 Unlike heat flux, there is no physical justification for the feedback term in \autoref{eq:sbc_dmp_emp} as
 the atmosphere does not care about ocean surface salinity \citep{Madec1997}.
+the atmosphere does not care about ocean surface salinity \citep{madec.delecluse_IWN97}.
 The SSS restoring term should be viewed as a flux correction on freshwater fluxes to
 reduce the uncertainties we have on the observed freshwater budget.

NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex

-                      r10442
+                      r11123
 The stochastic parametrization module aims to explicitly simulate uncertainties in the model.
 More particularly, \cite{Brankart_OM2013} has shown that,
+More particularly, \cite{brankart_OM13} has shown that,
 because of the nonlinearity of the seawater equation of state, unresolved scales represent a major source of
 uncertainties in the computation of the large scale horizontal density gradient (from T/S large scale fields),
 …
 A generic approach is thus to add one single new module in NEMO,
 generating processes with appropriate statistics to simulate each kind of uncertainty in the model
 (see \cite{Brankart_al_GMD2015} for more details).
+(see \cite{brankart.candille.ea_GMD15} for more details).
 In practice, at every model grid point,

NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex

-                      r10544
+                      r11123
 Nevertheless, in the latter case, it is achieved to a good approximation since
 the non-conservative term is the product of the time derivative of the tracer and the free surface height,
 two quantities that are not correlated \citep{Roullet_Madec_JGR00, Griffies_al_MWR01, Campin2004}.
 The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco}) is
+two quantities that are not correlated \citep{roullet.madec_JGR00, griffies.pacanowski.ea_MWR01, campin.adcroft.ea_OM04}.
+The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco?}) is
 the centred (\textit{now}) \textit{effective} ocean velocity, \ie the \textit{eulerian} velocity
 (see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) and/or
 …
 \end{equation}
 In the vertical direction (\np{nn\_cen\_v}~\forcode{= 4}),
 a $4^{th}$ COMPACT interpolation has been prefered \citep{Demange_PhD2014}.
+a $4^{th}$ COMPACT interpolation has been prefered \citep{demange_phd14}.
 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion,
 spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}.
+spectral characteristics similar to schemes of higher order \citep{lele_JCP92}.
 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but
 …
 (\ie it depends on the setting of \np{nn\_fct\_h} and \np{nn\_fct\_v}).
 There exist many ways to define $c_u$, each corresponding to a different FCT scheme.
 The one chosen in \NEMO is described in \citet{Zalesak_JCP79}.
+The one chosen in \NEMO is described in \citet{zalesak_JCP79}.
 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field.
 The resulting scheme is quite expensive but \textit{positive}.
 It can be used on both active and passive tracers.
 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}.
+A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{levy.estublier.ea_GRL01}.
 An additional option has been added controlled by \np{nn\_fct\_zts}.
 …
 a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction, but on the latter,
 a split-explicit time stepping is used, with a number of sub-timestep equals to \np{nn\_fct\_zts}.
 This option can be useful when the size of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.
+This option can be useful when the size of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.
 Note that in this case, a similar split-explicit time stepping should be used on vertical advection of momentum to
 insure a better stability (see \autoref{subsec:DYN_zad}).
 …
 MUSCL implementation can be found in the \mdl{traadv\_mus} module.
 MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}.
+MUSCL has been first implemented in \NEMO by \citet{levy.estublier.ea_GRL01}.
 In its formulation, the tracer at velocity points is evaluated assuming a linear tracer variation between
 two $T$-points (\autoref{fig:adv_scheme}).
 …
 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error
 \citep{Shchepetkin_McWilliams_OM05}.
 The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}.
+\citep{shchepetkin.mcwilliams_OM05}.
+The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}.
 It is a relatively good compromise between accuracy and smoothness.
 Nevertheless the scheme is not \textit{positive}, meaning that false extrema are permitted,
 …
 The intrinsic diffusion of UBS makes its use risky in the vertical direction where
 the control of artificial diapycnal fluxes is of paramount importance
 \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}.
+\citep{shchepetkin.mcwilliams_OM05, demange_phd14}.
 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme
 (\np{nn\_cen\_v}~\forcode{= 2 or 4}).
 …
 (which is the diffusive part of the scheme),
 is evaluated using the \textit{before} tracer (forward in time).
 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the QUICK advection scheme.
+This choice is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the QUICK advection scheme.
 UBS and QUICK schemes only differ by one coefficient.
 Replacing 1/6 with 1/8 in \autoref{eq:tra_adv_ubs} leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
+Replacing 1/6 with 1/8 in \autoref{eq:tra_adv_ubs} leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
 …
 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme
 proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}~\forcode{= .true.}.
+proposed by \citet{leonard_CMAME79} is used when \np{ln\_traadv\_qck}~\forcode{= .true.}.
 QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
 QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST limiter
 \citep{Leonard1991}.
+\citep{leonard_CMAME91}.
 It has been implemented in NEMO by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.
 The resulting scheme is quite expensive but \textit{positive}.
 …
 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}).
 When \np{ln\_traldf\_msc}~\forcode{= .true.}, a Method of Stabilizing Correction is used in which
 the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
+the pure vertical component is split into an explicit and an implicit part \citep{lemarie.debreu.ea_OM12}.
 % -------------------------------------------------------------------------------------------------------------
 …
 This formulation conserves the tracer but does not ensure the decrease of the tracer variance.
 Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without
 any additional background horizontal diffusion \citep{Guilyardi_al_CD01}.
+any additional background horizontal diffusion \citep{guilyardi.madec.ea_CD01}.
 Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}),
 …
 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}~\forcode{= .true.}; see \autoref{apdx:triad})
 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
+An alternative scheme developed by \cite{griffies.gnanadesikan.ea_JPO98} which ensures tracer variance decreases
 is also available in \NEMO (\np{ln\_traldf\_grif}~\forcode{= .true.}).
 A complete description of the algorithm is given in \autoref{apdx:triad}.
 …
 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface.
 In the linear free surface case, there is a small imbalance.
 The imbalance is larger than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.
+The imbalance is larger than the imbalance associated with the Asselin time filter \citep{leclair.madec_OM09}.
 This is the reason why the modified filter is not applied in the linear free surface case (see \autoref{chap:STP}).
 …
 In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}~\forcode{= .true.})
 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths,
 leading to the following expression \citep{Paulson1977}:
+leading to the following expression \citep{paulson.simpson_JPO77}:
 \[
   % \label{eq:traqsr_iradiance}
 …
 Such assumptions have been shown to provide a very crude and simplistic representation of
 observed light penetration profiles (\cite{Morel_JGR88}, see also \autoref{fig:traqsr_irradiance}).
+observed light penetration profiles (\cite{morel_JGR88}, see also \autoref{fig:traqsr_irradiance}).
 Light absorption in the ocean depends on particle concentration and is spectrally selective.
 \cite{Morel_JGR88} has shown that an accurate representation of light penetration can be provided by
+\cite{morel_JGR88} has shown that an accurate representation of light penetration can be provided by
 a 61 waveband formulation.
 Unfortunately, such a model is very computationally expensive.
 Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this formulation in which
+Thus, \cite{lengaigne.menkes.ea_CD07} have constructed a simplified version of this formulation in which
 visible light is split into three wavebands: blue (400-500 nm), green (500-600 nm) and red (600-700nm).
 For each wave-band, the chlorophyll-dependent attenuation coefficient is fitted to the coefficients computed from
 the full spectral model of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}),
+the full spectral model of \cite{morel_JGR88} (as modified by \cite{morel.maritorena_JGR01}),
 assuming the same power-law relationship.
 As shown in \autoref{fig:traqsr_irradiance}, this formulation, called RGB (Red-Green-Blue),
 …
 \item[\np{nn\_chdta}~\forcode{= 2}]
   same as previous case except that a vertical profile of chlorophyl is used.
   Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value;
+  Following \cite{morel.berthon_LO89}, the profile is computed from the local surface chlorophyll value;
 \item[\np{ln\_qsr\_bio}~\forcode{= .true.}]
   simulated time varying chlorophyll by TOP biogeochemical model.
 …
 waveband Morel (1988) formulation (black) for a chlorophyll concentration of
       (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$.
       From \citet{Lengaigne_al_CD07}.
+      From \citet{lengaigne.menkes.ea_CD07}.
+    }
   \end{center}
 …
     \caption{
       \protect\label{fig:geothermal}
       Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
       It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}.
+      Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{emile-geay.madec_OS09}.
+      It is inferred from the age of the sea floor and the formulae of \citet{stein.stein_N92}.
+    }
   \end{center}
 …
 This is the default option in \NEMO, and it is implemented using the masking technique.
 However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling.
 This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{Stein_Stein_Nat92}),
+This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{stein.stein_N92}),
 but it warms systematically the ocean and acts on the densest water masses.
 Taking this flux into account in a global ocean model increases the deepest overturning cell
 (\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}.
+(\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{emile-geay.madec_OS09}.
 Options are defined through the  \ngn{namtra\_bbc} namelist variables.
 …
 the \np{nn\_geoflx\_cst}, which is also a namelist parameter.
 When \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is introduced which is provided in
 the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:geothermal}) \citep{Emile-Geay_Madec_OS09}.
+the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:geothermal}) \citep{emile-geay.madec_OS09}.
 % ================================================================
 …
 sometimes over a thickness much larger than the thickness of the observed gravity plume.
 A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of
 a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved.
 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{Beckmann_Doscher1997},
+a sill \citep{willebrand.barnier.ea_PO01}, and the thickness of the plume is not resolved.
+The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{beckmann.doscher_JPO97},
 is to allow a direct communication between two adjacent bottom cells at different levels,
 whenever the densest water is located above the less dense water.
 …
 In the current implementation of the BBL, only the tracers are modified, not the velocities.
 Furthermore, it only connects ocean bottom cells, and therefore does not include all the improvements introduced by
 \citet{Campin_Goosse_Tel99}.
+\citet{campin.goosse_T99}.
 % -------------------------------------------------------------------------------------------------------------
 …
 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, and
 $A_l^\sigma$ the lateral diffusivity in the BBL.
 Following \citet{Beckmann_Doscher1997}, the latter is prescribed with a spatial dependence,
+Following \citet{beckmann.doscher_JPO97}, the latter is prescribed with a spatial dependence,
 \ie in the conditional form
 \begin{equation}
 …
 \np{nn\_bbl\_adv}~\forcode{= 1}:
 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step
 (see black arrow in \autoref{fig:bbl}) \citep{Beckmann_Doscher1997}.
+(see black arrow in \autoref{fig:bbl}) \citep{beckmann.doscher_JPO97}.
 It is a \textit{conditional advection}, that is, advection is allowed only
 if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$) and
 …
 \np{nn\_bbl\_adv}~\forcode{= 2}:
 the downslope velocity is chosen to be proportional to $\Delta \rho$,
 the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}.
+the density difference between the higher cell and lower cell densities \citep{campin.goosse_T99}.
 The advection is allowed only  if dense water overlies less dense water on the slope
 (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$).
 …
 The parameter $\gamma$ should take a different value for each bathymetric step, but for simplicity,
 and because no direct estimation of this parameter is available, a uniform value has been assumed.
 The possible values for $\gamma$ range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}.
+The possible values for $\gamma$ range between 1 and $10~s$ \citep{campin.goosse_T99}.
 Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using the upwind scheme.
 …
 In the vicinity of these walls, $\gamma$ takes large values (equivalent to a time scale of a few days) whereas
 it is zero in the interior of the model domain.
 The second case corresponds to the use of the robust diagnostic method \citep{Sarmiento1982}.
+The second case corresponds to the use of the robust diagnostic method \citep{sarmiento.bryan_JGR82}.
 It allows us to find the velocity field consistent with the model dynamics whilst
 having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$).
 …
 only below the mixed layer (defined either on a density or $S_o$ criterion).
 It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here
 \citep{Madec_al_JPO96}.
+\citep{madec.delecluse.ea_JPO96}.
 For generating \ifile{resto}, see the documentation for the DMP tool provided with the source code under
 …
 Options are defined through the \ngn{namdom} namelist variables.
 The general framework for tracer time stepping is a modified leap-frog scheme \citep{Leclair_Madec_OM09},
+The general framework for tracer time stepping is a modified leap-frog scheme \citep{leclair.madec_OM09},
 \ie a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:STP_mLF}):
 \begin{equation}
 …
 Nonlinearities of the EOS are of major importance, in particular influencing the circulation through
 determination of the static stability below the mixed layer,
 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{Roquet_JPO2015}.
 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983}) or
 TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real ocean circulation is attempted
 \citep{Roquet_JPO2015}.
+thus controlling rates of exchange between the atmosphere and the ocean interior \citep{roquet.madec.ea_JPO15}.
+Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{fofonoff.millard_bk83}) or
+TEOS-10 \citep{ioc.iapso_bk10} standards should be used anytime a simulation of the real ocean circulation is attempted
+\citep{roquet.madec.ea_JPO15}.
 The use of TEOS-10 is highly recommended because
 \textit{(i)}   it is the new official EOS,
 …
 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature and
 practical salinity for EOS-980, both variables being more suitable for use as model variables
 \citep{TEOS10, Graham_McDougall_JPO13}.
+\citep{ioc.iapso_bk10, graham.mcdougall_JPO13}.
 EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility.
 For process studies, it is often convenient to use an approximation of the EOS.
 To that purposed, a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available.
+To that purposed, a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is also available.
 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density.
 …
 This is a sensible choice for the reference density used in a Boussinesq ocean climate model, as,
 with the exception of only a small percentage of the ocean,
 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}.
+density in the World Ocean varies by no more than 2$\%$ from that value \citep{gill_bk82}.
 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which
 …
 \begin{description}
 \item[\np{nn\_eos}~\forcode{= -1}]
   the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.
+  the polyTEOS10-bsq equation of seawater \citep{roquet.madec.ea_OM15} is used.
   The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,
   but it is optimized for a boussinesq fluid and the polynomial expressions have simpler and
 …
   use in ocean models.
   Note that a slightly higher precision polynomial form is now used replacement of
   the TEOS-10 rational function approximation for hydrographic data analysis \citep{TEOS10}.
+  the TEOS-10 rational function approximation for hydrographic data analysis \citep{ioc.iapso_bk10}.
   A key point is that conservative state variables are used:
   Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \deg{C}, notation: $\Theta$).
   The pressure in decibars is approximated by the depth in meters.
   With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant.
   It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}.
+  It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{ioc.iapso_bk10}.
   Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$.
   In particular, the initial state deined by the user have to be given as \textit{Conservative} Temperature and
 …
   The pressure in decibars is approximated by the depth in meters.
   With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, salinity and
   pressure \citep{UNESCO1983}.
+  pressure \citep{fofonoff.millard_bk83}.
   Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which
   is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value.
 \item[\np{nn\_eos}~\forcode{= 1}]
   a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen,
+  a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is chosen,
   the coefficients of which has been optimized to fit the behavior of TEOS10
   (Roquet, personal comm.) (see also \citet{Roquet_JPO2015}).
+  (Roquet, personal comm.) (see also \citet{roquet.madec.ea_JPO15}).
   It provides a simplistic linear representation of both cabbeling and thermobaricity effects which
   is enough for a proper treatment of the EOS in theoretical studies \citep{Roquet_JPO2015}.
+  is enough for a proper treatment of the EOS in theoretical studies \citep{roquet.madec.ea_JPO15}.
   With such an equation of state there is no longer a distinction between
   \textit{conservative} and \textit{potential} temperature,
 …
 \label{subsec:TRA_fzp}
 The freezing point of seawater is a function of salinity and pressure \citep{UNESCO1983}:
+The freezing point of seawater is a function of salinity and pressure \citep{fofonoff.millard_bk83}:
 \begin{equation}
   \label{eq:tra_eos_fzp}

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex

-                      r10442
+                      r11123
 a dependency between the vertical eddy coefficients and the local Richardson number
 (\ie the ratio of stratification to vertical shear).
 Following \citet{Pacanowski_Philander_JPO81}, the following formulation has been implemented:
+Following \citet{pacanowski.philander_JPO81}, the following formulation has been implemented:
 \[
   % \label{eq:zdfric}
 …
 The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}.
 Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to
 the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{Lermusiaux2001}.
+the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{lermusiaux_JMS01}.
 % -------------------------------------------------------------------------------------------------------------
 …
 a prognostic equation for $\bar{e}$, the turbulent kinetic energy,
 and a closure assumption for the turbulent length scales.
 This turbulent closure model has been developed by \citet{Bougeault1989} in the atmospheric case,
 adapted by \citet{Gaspar1990} for the oceanic case, and embedded in OPA, the ancestor of NEMO,
 by \citet{Blanke1993} for equatorial Atlantic simulations.
 Since then, significant modifications have been introduced by \citet{Madec1998} in both the implementation and
+This turbulent closure model has been developed by \citet{bougeault.lacarrere_MWR89} in the atmospheric case,
+adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of NEMO,
+by \citet{blanke.delecluse_JPO93} for equatorial Atlantic simulations.
+Since then, significant modifications have been introduced by \citet{madec.delecluse.ea_NPM98} in both the implementation and
 the formulation of the mixing length scale.
 The time evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical shear,
 its destruction through stratification, its vertical diffusion, and its dissipation of \citet{Kolmogorov1942} type:
+its destruction through stratification, its vertical diffusion, and its dissipation of \citet{kolmogorov_IANS42} type:
 \begin{equation}
   \label{eq:zdftke_e}
 …
 $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients.
 The constants $C_k =  0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with
 vertical mixing at any depth \citep{Gaspar1990}.
+vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}.
 They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}.
 $P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function of the local Richardson number, $R_i$:
+$P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$:
 \begin{align*}
   % \label{eq:prt}
 …
 At the sea surface, the value of $\bar{e}$ is prescribed from the wind stress field as
 $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn\_ebb} namelist parameter.
 The default value of $e_{bb}$ is 3.75. \citep{Gaspar1990}), however a much larger value can be used when
+The default value of $e_{bb}$ is 3.75. \citep{gaspar.gregoris.ea_JGR90}), however a much larger value can be used when
 taking into account the surface wave breaking (see below Eq. \autoref{eq:ZDF_Esbc}).
 The bottom value of TKE is assumed to be equal to the value of the level just above.
 …
 the numerical scheme does not ensure its positivity.
 To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} namelist parameter).
 Following \citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$.
 This allows the subsequent formulations to match that of \citet{Gargett1984} for the diffusion in
+Following \citet{gaspar.gregoris.ea_JGR90}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$.
+This allows the subsequent formulations to match that of \citet{gargett_JMR84} for the diffusion in
 the thermocline and deep ocean :  $K_\rho = 10^{-3} / N$.
 In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical instabilities associated with
 …
 For computational efficiency, the original formulation of the turbulent length scales proposed by
 \citet{Gaspar1990} has been simplified.
+\citet{gaspar.gregoris.ea_JGR90} has been simplified.
 Four formulations are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist parameter.
 The first two are based on the following first order approximation \citep{Blanke1993}:
+The first two are based on the following first order approximation \citep{blanke.delecluse_JPO93}:
 \begin{equation}
   \label{eq:tke_mxl0_1}
 …
 The resulting length scale is bounded by the distance to the surface or to the bottom
 (\np{nn\_mxl}\forcode{ = 0}) or by the local vertical scale factor (\np{nn\_mxl}\forcode{ = 1}).
 \citet{Blanke1993} notice that this simplification has two major drawbacks:
+\citet{blanke.delecluse_JPO93} notice that this simplification has two major drawbacks:
 it makes no sense for locally unstable stratification and the computation no longer uses all
 the information contained in the vertical density profile.
 To overcome these drawbacks, \citet{Madec1998} introduces the \np{nn\_mxl}\forcode{ = 2..3} cases,
+To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{ = 2..3} cases,
 which add an extra assumption concerning the vertical gradient of the computed length scale.
 So, the length scales are first evaluated as in \autoref{eq:tke_mxl0_1} and then bounded such that:
 …
 \autoref{eq:tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than
 the variations of depth.
 It provides a better approximation of the \citet{Gaspar1990} formulation while being much less
+It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less
 time consuming.
 In particular, it allows the length scale to be limited not only by the distance to the surface or
 …
 In the \np{nn\_mxl}\forcode{ = 2} case, the dissipation and mixing length scales take the same value:
 $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the \np{nn\_mxl}\forcode{ = 3} case,
 the dissipation and mixing turbulent length scales are give as in \citet{Gaspar1990}:
+the dissipation and mixing turbulent length scales are give as in \citet{gaspar.gregoris.ea_JGR90}:
 \[
   % \label{eq:tke_mxl_gaspar}
 …
 Usually the surface scale is given by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant and
 $z_o$ the roughness parameter of the surface.
 Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} leads to a 0.04~m, the default value of \np{rn\_mxl0}.
+Assuming $z_o=0.1$~m \citep{craig.banner_JPO94} leads to a 0.04~m, the default value of \np{rn\_mxl0}.
 In the ocean interior a minimum length scale is set to recover the molecular viscosity when
 $\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
 …
 %-----------------------------------------------------------------------%
 Following \citet{Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified to
+Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to
 include the effect of surface wave breaking energetics.
 This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow.
 The \citet{Mellor_Blumberg_JPO04} modifications acts on surface length scale and TKE values and
+The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and
 air-sea drag coefficient.
 The latter concerns the bulk formulea and is not discussed here.
 Following \citet{Craig_Banner_JPO94}, the boundary condition on surface TKE value is :
+Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is :
 \begin{equation}
   \label{eq:ZDF_Esbc}
   \bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
 \end{equation}
 where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of proportionality which depends on the ''wave age'',
 ranging from 57 for mature waves to 146 for younger waves \citep{Mellor_Blumberg_JPO04}.
+where $\alpha_{CB}$ is the \citet{craig.banner_JPO94} constant of proportionality which depends on the ''wave age'',
+ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}.
 The boundary condition on the turbulent length scale follows the Charnock's relation:
 \begin{equation}
 …
 \end{equation}
 where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant.
 \citet{Mellor_Blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by
 \citet{Stacey_JPO99} citing observation evidence, and
+\citet{mellor.blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by
+\citet{stacey_JPO99} citing observation evidence, and
 $\alpha_{CB} = 100$ the Craig and Banner's value.
 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$,
 …
 Although LC have nothing to do with convection, the circulation pattern is rather similar to
 so-called convective rolls in the atmospheric boundary layer.
 The detailed physics behind LC is described in, for example, \citet{Craik_Leibovich_JFM76}.
+The detailed physics behind LC is described in, for example, \citet{craik.leibovich_JFM76}.
 The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and
 wind drift currents.
 Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by
 \citep{Axell_JGR02} for a $k-\epsilon$ turbulent closure.
+\citep{axell_JGR02} for a $k-\epsilon$ turbulent closure.
 The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in
 an extra source terms of TKE, $P_{LC}$.
 …
 \forcode{.true.} in the namtke namelist.
 By making an analogy with the characteristic convective velocity scale (\eg, \citet{D'Alessio_al_JPO98}),
+By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}),
 $P_{LC}$ is assumed to be :
 \[
 …
 With no information about the wave field, $w_{LC}$ is assumed to be proportional to
 the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module
 \footnote{Following \citet{Li_Garrett_JMR93}, the surface Stoke drift velocity may be expressed as
+\footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as
   $u_s =  0.016 \,|U_{10m}|$.
   Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of
 …
   \end{cases}
 \]
 where $c_{LC} = 0.15$ has been chosen by \citep{Axell_JGR02} as a good compromise to fit LES data.
+where $c_{LC} = 0.15$ has been chosen by \citep{axell_JGR02} as a good compromise to fit LES data.
 The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second.
 The value of $c_{LC}$ is set through the \np{rn\_lc} namelist parameter,
 having in mind that it should stay between 0.15 and 0.54 \citep{Axell_JGR02}.
+having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}.
 The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations:
 …
 produce mixed-layer depths that are too shallow during summer months and windy conditions.
 This bias is particularly acute over the Southern Ocean.
 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{Rodgers_2014}.
+To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}.
 The parameterization is an empirical one, \ie not derived from theoretical considerations,
 but rather is meant to account for observed processes that affect the density structure of
 …
 (first line in \autoref{eq:PE_zdf}).
 To do so a special care have to be taken for both the time and space discretization of
 the TKE equation \citep{Burchard_OM02,Marsaleix_al_OM08}.
+the TKE equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}.
 Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows how
 …
 The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations:
 one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale,
 $\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}.
+$\psi$ \citep{umlauf.burchard_JMR03, umlauf.burchard_CSR05}.
 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$,
 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:GLS} allows to recover a number of
 well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982}, $k$-$\epsilon$ \citep{Rodi_1987},
 $k$-$\omega$ \citep{Wilcox_1988} among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}).
+well-known turbulent closures ($k$-$kl$ \citep{mellor.yamada_RG82}, $k$-$\epsilon$ \citep{rodi_JGR87},
+$k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}).
 The GLS scheme is given by the following set of equations:
 \begin{equation}
 …
     \begin{tabular}{ccccc}
       &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\
       % & \citep{Mellor_Yamada_1982} &  \citep{Rodi_1987}       & \citep{Wilcox_1988} &                 \\
+      % & \citep{mellor.yamada_RG82} &  \citep{rodi_JGR87}       & \citep{wilcox_AJ88} &                 \\
       \hline
       \hline
 …
 the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length) value near physical boundaries
 (logarithmic boundary layer law).
 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{Galperin_al_JAS88},
 or by \citet{Kantha_Clayson_1994} or one of the two functions suggested by \citet{Canuto_2001}
+$C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88},
+or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01}
 (\np{nn\_stab\_func}\forcode{ = 0..3}, resp.).
 The value of $C_{0\mu}$ depends of the choice of the stability function.
 …
 Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp.
 As for TKE closure, the wave effect on the mixing is considered when
 \np{ln\_crban}\forcode{ = .true.} \citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}.
+\np{ln\_crban}\forcode{ = .true.} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}.
 The \np{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and
 \np{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}.
 …
 almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy.
 With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$.
 A value of $c_{lim} = 0.53$ is often used \citep{Galperin_al_JAS88}.
 \cite{Umlauf_Burchard_CSR05} show that the value of the clipping factor is of crucial importance for
+A value of $c_{lim} = 0.53$ is often used \citep{galperin.kantha.ea_JAS88}.
+\cite{umlauf.burchard_CSR05} show that the value of the clipping factor is of crucial importance for
 the entrainment depth predicted in stably stratified situations,
 and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes.
 …
 The time and space discretization of the GLS equations follows the same energetic consideration as for
 the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{Burchard_OM02}.
 Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}.
+the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}.
+Examples of performance of the 4 turbulent closure scheme can be found in \citet{warner.sherwood.ea_OM05}.
 % -------------------------------------------------------------------------------------------------------------
 …
 the water column, but only until the density structure becomes neutrally stable
 (\ie until the mixed portion of the water column has \textit{exactly} the density of the water just below)
 \citep{Madec_al_JPO91}.
+\citep{madec.delecluse.ea_JPO91}.
 The associated algorithm is an iterative process used in the following way (\autoref{fig:npc}):
 starting from the top of the ocean, the first instability is found.
 …
 the algorithm used in \NEMO converges for any profile in a number of iterations which is less than
 the number of vertical levels.
 This property is of paramount importance as pointed out by \citet{Killworth1989}:
+This property is of paramount importance as pointed out by \citet{killworth_iprc89}:
 it avoids the existence of permanent and unrealistic static instabilities at the sea surface.
 This non-penetrative convective algorithm has been proved successful in studies of the deep water formation in
 the north-western Mediterranean Sea \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.
+the north-western Mediterranean Sea \citep{madec.delecluse.ea_JPO91, madec.chartier.ea_DAO91, madec.crepon_iprc91}.
 The current implementation has been modified in order to deal with any non linear equation of seawater
 …
 In this case, the vertical eddy mixing coefficients are assigned very large values
 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable
 (\ie when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar_PhD97, Lazar_al_JPO99}.
+(\ie when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}.
 This is done either on tracers only (\np{nn\_evdm}\forcode{ = 0}) or
 on both momentum and tracers (\np{nn\_evdm}\forcode{ = 1}).
 …
 Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$.
 This removes a potential source of divergence of odd and even time step in
 a leapfrog environment \citep{Leclair_PhD2010} (see \autoref{sec:STP_mLF}).
+a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:STP_mLF}).
 % -------------------------------------------------------------------------------------------------------------
 …
 The former condition leads to salt fingering and the latter to diffusive convection.
 Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the ocean.
 \citet{Merryfield1999} include a parameterisation of such phenomena in a global ocean model and show that
+\citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that
 it leads to relatively minor changes in circulation but exerts significant regional influences on
 temperature and salinity.
 …
     \caption{
       \protect\label{fig:zdfddm}
       From \citet{Merryfield1999} :
+      From \citet{merryfield.holloway.ea_JPO99} :
       (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in regions of salt fingering.
       Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$;
 …
 The factor 0.7 in \autoref{eq:zdfddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx  0.7$ of
 buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{McDougall_Taylor_JMR84}).
 Following  \citet{Merryfield1999}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
+buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{mcdougall.taylor_JMR84}).
+Following  \citet{merryfield.holloway.ea_JPO99}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
 To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by
 …
 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean,
 and setting $r = H / \tau$, where $H$ is the ocean depth.
 Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}.
+Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}.
 A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used in quasi-geostrophic models.
 One may consider the linear friction as an approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$
 (\citet{Gill1982}, Eq. 9.6.6).
+(\citet{gill_bk82}, Eq. 9.6.6).
 For example, with a drag coefficient $C_D = 0.002$, a typical speed of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$,
 and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.
 …
 internal waves breaking and other short time scale currents.
 A typical value of the drag coefficient is $C_D = 10^{-3} $.
 As an example, the CME experiment \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and
 $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} uses $C_D = 1.4\;10^{-3}$ and
+As an example, the CME experiment \citep{treguier_JGR92} uses $C_D = 10^{-3}$ and
+$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{killworth_JPO92} uses $C_D = 1.4\;10^{-3}$ and
 $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.
 The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} namelist parameters).
 …
 Options are defined through the  \ngn{namzdf\_tmx} namelist variables.
 The parameterization of tidal mixing follows the general formulation for the vertical eddy diffusivity proposed by
 \citet{St_Laurent_al_GRL02} and first introduced in an OGCM by \citep{Simmons_al_OM04}.
+\citet{st-laurent.simmons.ea_GRL02} and first introduced in an OGCM by \citep{simmons.jayne.ea_OM04}.
 In this formulation an additional vertical diffusivity resulting from internal tide breaking,
 $A^{vT}_{tides}$ is expressed as a function of $E(x,y)$,
 …
 with the remaining $1-q$ radiating away as low mode internal waves and
 contributing to the background internal wave field.
 A value of $q=1/3$ is typically used \citet{St_Laurent_al_GRL02}.
+A value of $q=1/3$ is typically used \citet{st-laurent.simmons.ea_GRL02}.
 The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical.
 It is implemented as a simple exponential decaying upward away from the bottom,
 with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter,
 with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04},
+with a typical value of $500\,m$) \citep{st-laurent.nash_DSR04},
 \[
   % \label{eq:Fz}
 …
 the unrepresented internal waves induced by the tidal flow over rough topography in a stratified ocean.
 In the current version of \NEMO, the map is built from the output of
 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.
+the barotropic global ocean tide model MOG2D-G \citep{carrere.lyard_GRL03}.
 This model provides the dissipation associated with internal wave energy for the M2 and K1 tides component
 (\autoref{fig:ZDF_M2_K1_tmx}).
 …
 The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$.
 Its global mean value is $1.1$ TW,
 in agreement with independent estimates \citep{Egbert_Ray_Nat00, Egbert_Ray_JGR01}.
+in agreement with independent estimates \citep{egbert.ray_N00, egbert.ray_JGR01}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
     \caption{
       \protect\label{fig:ZDF_M2_K1_tmx}
       (a) M2 and (b) K1 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$).
+      (a) M2 and (b) K1 internal wave drag energy from \citet{carrere.lyard_GRL03} ($W/m^2$).
+    }
   \end{center}
 …
 When \np{ln\_tmx\_itf}\forcode{ = .true.}, the two key parameters $q$ and $F(z)$ are adjusted following
 the parameterisation developed by \citet{Koch-Larrouy_al_GRL07}:
+the parameterisation developed by \citet{koch-larrouy.madec.ea_GRL07}:
 First, the Indonesian archipelago is a complex geographic region with a series of
 …
 Second, the vertical structure function, $F(z)$, is no more associated with a bottom intensification of the mixing,
 but with a maximum of energy available within the thermocline.
 \citet{Koch-Larrouy_al_GRL07} have suggested that the vertical distribution of
+\citet{koch-larrouy.madec.ea_GRL07} have suggested that the vertical distribution of
 the energy dissipation proportional to $N^2$ below the core of the thermocline and to $N$ above.
 The resulting $F(z)$ is:
 …
 Introduced in a regional OGCM, the parameterization improves the water mass characteristics in
 the different Indonesian seas, suggesting that the horizontal and vertical distributions of
 the mixing are adequately prescribed \citep{Koch-Larrouy_al_GRL07, Koch-Larrouy_al_OD08a, Koch-Larrouy_al_OD08b}.
+the mixing are adequately prescribed \citep{koch-larrouy.madec.ea_GRL07, koch-larrouy.madec.ea_OD08*a, koch-larrouy.madec.ea_OD08*b}.
 Note also that such a parameterisation has a significant impact on the behaviour of
 global coupled GCMs \citep{Koch-Larrouy_al_CD10}.
+global coupled GCMs \citep{koch-larrouy.lengaigne.ea_CD10}.
 % ================================================================
 …
 The parameterization of mixing induced by breaking internal waves is a generalization of
 the approach originally proposed by \citet{St_Laurent_al_GRL02}.
+the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}.
 A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,
 and the resulting diffusivity is obtained as
 …
 the energy available for mixing.
 If the \np{ln\_mevar} namelist parameter is set to false, the mixing efficiency is taken as constant and
 equal to 1/6 \citep{Osborn_JPO80}.
+equal to 1/6 \citep{osborn_JPO80}.
 In the opposite (recommended) case, $R_f$ is instead a function of
 the turbulence intensity parameter $Re_b = \frac{ \epsilon}{\nu \, N^2}$,
 with $\nu$ the molecular viscosity of seawater, following the model of \cite{Bouffard_Boegman_DAO2013} and
 the implementation of \cite{de_lavergne_JPO2016_efficiency}.
+with $\nu$ the molecular viscosity of seawater, following the model of \cite{bouffard.boegman_DAO13} and
+the implementation of \cite{de-lavergne.madec.ea_JPO16}.
 Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when
 the mixing efficiency is constant.
 …
 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary
 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice.
 This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014},
 is implemented as in \cite{de_lavergne_JPO2016_efficiency}.
+This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14},
+is implemented as in \cite{de-lavergne.madec.ea_JPO16}.
 The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$,
 …
 $h_{cri}$ is related to the large-scale topography of the ocean (etopo2) and
 $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of
 the abyssal hill topography \citep{Goff_JGR2010} and the latitude.
+the abyssal hill topography \citep{goff_JGR10} and the latitude.
 % ================================================================

NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex

-                      r10442
+                      r11123
 horizontal kinetic energy and/or potential enstrophy of horizontally non-divergent flow,
 and variance of temperature and salinity will be conserved in the absence of dissipative effects and forcing.
 \citet{Arakawa1966} has first pointed out the advantage of this approach.
+\citet{arakawa_JCP66} has first pointed out the advantage of this approach.
 He showed that if integral constraints on energy are maintained,
 the computation will be free of the troublesome "non linear" instability originally pointed out by
 \citet{Phillips1959}.
+\citet{phillips_TAMS59}.
 A consistent formulation of the energetic properties is also extremely important in carrying out
 long-term numerical simulations for an oceanographic model.
 Such a formulation avoids systematic errors that accumulate with time \citep{Bryan1997}.
+Such a formulation avoids systematic errors that accumulate with time \citep{bryan_JCP97}.
 The general philosophy of OPA which has led to the discrete formulation presented in {\S}II.2 and II.3 is to
 …
 Note that in some very specific cases as passive tracer studies, the positivity of the advective scheme is required.
 In that case, and in that case only, the advective scheme used for passive tracer is a flux correction scheme
 \citep{Marti1992, Levy1996, Levy1998}.
+\citep{Marti1992?, Levy1996?, Levy1998?}.
 % -------------------------------------------------------------------------------------------------------------

NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex

r10601	r11123
272	272	and their propagation and accumulation cause uncertainty in final simulation reproducibility on
273	273	different numbers of processors.
274		To avoid so, based on \citet{~~He_Ding_JSC~~01} review of different technics,
	274	To avoid so, based on \citet{he.ding_JS01} review of different technics,
275	275	we use a so called self-compensated summation method.
276	276	The idea is to estimate the roundoff error, store it in a buffer, and then add it back in the next addition.

NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex

-                      r10544
+                      r11123
 If further, an approximative conservation of heat and salt contents is sufficient for the problem solved,
 then it is sufficient to solve a linearized version of \autoref{eq:PE_ssh},
 which still allows to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}.
+which still allows to take into account freshwater fluxes applied at the ocean surface \citep{roullet.madec_JGR00}.
 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost.
 The filtering of EGWs in models with a free surface is usually a matter of discretisation of
 the temporal derivatives,
 using a split-explicit method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92} or
 the implicit scheme \citep{Dukowicz1994} or
 the addition of a filtering force in the momentum equation \citep{Roullet_Madec_JGR00}.
+using a split-explicit method \citep{killworth.webb.ea_JPO91, zhang.endoh_JGR92} or
+the implicit scheme \citep{dukowicz.smith_JGR94} or
+the addition of a filtering force in the momentum equation \citep{roullet.madec_JGR00}.
 With the present release, \NEMO offers the choice between
 an explicit free surface (see \autoref{subsec:DYN_spg_exp}) or
 a split-explicit scheme strongly inspired the one proposed by \citet{Shchepetkin_McWilliams_OM05}
+a split-explicit scheme strongly inspired the one proposed by \citet{shchepetkin.mcwilliams_OM05}
 (see \autoref{subsec:DYN_spg_ts}).
 …
 cannot be easily treated in a global model without filtering.
 A solution consists of introducing an appropriate coordinate transformation that
 shifts the singular point onto land \citep{Madec_Imbard_CD96, Murray_JCP96}.
+shifts the singular point onto land \citep{madec.imbard_CD96, murray_JCP96}.
 As a consequence, it is important to solve the primitive equations in various curvilinear coordinate systems.
 An efficient way of introducing an appropriate coordinate transform can be found when using a tensorial formalism.
 …
 Ocean modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth approximation),
 with preservation of the local vertical. Here we give the simplified equations for this particular case.
 The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics.
+The general case is detailed by \citet{eiseman.stone_SR80} in their survey of the conservation laws of fluid dynamics.
 Let $(i,j,k)$ be a set of orthogonal curvilinear coordinates on
 …
 In order to satisfy two or more constrains one can even be tempted to mixed these coordinate systems, as in
 HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at
 the ocean bottom) \citep{Chassignet_al_JPO03} or
+the ocean bottom) \citep{chassignet.smith.ea_JPO03} or
 OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and $\sigma$-coordinate elsewhere)
 \citep{Madec_al_JPO96} among others.
+\citep{madec.delecluse.ea_JPO96} among others.
 In fact one is totally free to choose any space and time vertical coordinate by
 …
 the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through
 \autoref{eq:PE_s}.
 This so-called \textit{generalised vertical coordinate} \citep{Kasahara_MWR74} is in fact
+This so-called \textit{generalised vertical coordinate} \citep{kasahara_MWR74} is in fact
 an Arbitrary Lagrangian--Eulerian (ALE) coordinate.
 Indeed, choosing an expression for $s$ is an arbitrary choice that determines
 which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and
 which part will be used to move them (Lagrangian part).
 The coordinate is also sometime referenced as an adaptive coordinate \citep{Hofmeister_al_OM09},
+The coordinate is also sometime referenced as an adaptive coordinate \citep{hofmeister.burchard.ea_OM10},
 since the coordinate system is adapted in the course of the simulation.
 Its most often used implementation is via an ALE algorithm,
 in which a pure lagrangian step is followed by regridding and remapping steps,
 the later step implicitly embedding the vertical advection
 \citep{Hirt_al_JCP74, Chassignet_al_JPO03, White_al_JCP09}.
 Here we follow the \citep{Kasahara_MWR74} strategy:
+\citep{hirt.amsden.ea_JCP74, chassignet.smith.ea_JPO03, white.adcroft.ea_JCP09}.
+Here we follow the \citep{kasahara_MWR74} strategy:
 a regridding step (an update of the vertical coordinate) followed by an eulerian step with
 an explicit computation of vertical advection relative to the moving s-surfaces.
 …
       (b) $z$-coordinate in non-linear free surface case ;
       (c) re-scaled height coordinate
       (become popular as the \zstar-coordinate \citep{Adcroft_Campin_OM04}).
+      (become popular as the \zstar-coordinate \citep{adcroft.campin_OM04}).
+    }
   \end{center}
 …
 In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account.
 These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site.
+These coordinates systems is presented in a report \citep{levier.treguier.ea_rpt07} available on the \NEMO web site.
 The \zstar coordinate approach is an unapproximated, non-linear free surface implementation which allows one to
 deal with large amplitude free-surface variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}.
+deal with large amplitude free-surface variations relative to the vertical resolution \citep{adcroft.campin_OM04}.
 In the \zstar formulation,
 the variation of the column thickness due to sea-surface undulations is not concentrated in the surface level,
 …
 The quasi -horizontal nature of the coordinate surfaces also facilitates the implementation of
 neutral physics parameterizations in \zstar models using the same techniques as in $z$-models
 (see Chapters 13-16 of \cite{Griffies_Bk04}) for a discussion of neutral physics in $z$-models,
+(see Chapters 13-16 of \cite{griffies_bk04}) for a discussion of neutral physics in $z$-models,
 as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO).
 …
 The response to such a velocity field often leads to numerical dispersion effects.
 One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of
 a full step one \cite{Pacanowski_Gnanadesikan_MWR98}.
+a full step one \cite{pacanowski.gnanadesikan_MWR98}.
 Another solution is to introduce a terrain-following coordinate system (hereafter $s$-coordinate).
 …
 introduces a truncation error that is not present in a $z$-model.
 In the special case of a $\sigma$-coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$),
 \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error.
+\citet{haney_JPO91} and \citet{beckmann.haidvogel_JPO93} have given estimates of the magnitude of this truncation error.
 It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state,
 and the finite difference scheme.
 …
 The large-scale slopes require high horizontal resolution, and the computational cost becomes prohibitive.
 This problem can be at least partially overcome by mixing $s$-coordinate and
 step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec_al_JPO96}.
+step-like representation of bottom topography \citep{gerdes_JGR93*a,gerdes_JGR93*b,madec.delecluse.ea_JPO96}.
 However, the definition of the model domain vertical coordinate becomes then a non-trivial thing for
 a realistic bottom topography:
 …
 In contrast, the ocean will stay at rest in a $z$-model.
 As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below
 the strongly stratified portion of the water column (\ie the main thermocline) \citep{Madec_al_JPO96}.
+the strongly stratified portion of the water column (\ie the main thermocline) \citep{madec.delecluse.ea_JPO96}.
 An alternate solution consists of rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces
 (see \autoref{subsec:PE_ldf}).
 …
 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}).
 The $s$-coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from
+The $s$-coordinates introduced here \citep{lott.madec.ea_OM90,madec.delecluse.ea_JPO96} differ mainly in two aspects from
 similar models:
 it allows a representation of bottom topography with mixed full or partial step-like/terrain following topography;
 …
 \label{subsec:PE_zco_tilde}
 The \ztilde -coordinate has been developed by \citet{Leclair_Madec_OM11}.
+The \ztilde -coordinate has been developed by \citet{leclair.madec_OM11}.
 It is available in \NEMO since the version 3.4.
 Nevertheless, it is currently not robust enough to be used in all possible configurations.
 …
 The resulting lateral diffusive and dissipative operators are of second order.
 Observations show that lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces
 (or more precisely neutral surfaces \cite{McDougall1987}) rather than across them.
+(or more precisely neutral surfaces \cite{mcdougall_JPO87}) rather than across them.
 As the slope of neutral surfaces is small in the ocean, a common approximation is to assume that
 the `lateral' direction is the horizontal, \ie the lateral mixing is performed along geopotential surfaces.
 …
 both horizontal and isoneutral operators have no effect on mean (\ie large scale) potential energy whereas
 potential energy is a main source of turbulence (through baroclinic instabilities).
 \citet{Gent1990} have proposed a parameterisation of mesoscale eddy-induced turbulence which
+\citet{gent.mcwilliams_JPO90} have proposed a parameterisation of mesoscale eddy-induced turbulence which
 associates an eddy-induced velocity to the isoneutral diffusion.
 Its mean effect is to reduce the mean potential energy of the ocean.
 …
 There are not all available in \NEMO. For active tracers (temperature and salinity) the main ones are:
 Laplacian and bilaplacian operators acting along geopotential or iso-neutral surfaces,
 \citet{Gent1990} parameterisation, and various slightly diffusive advection schemes.
+\citet{gent.mcwilliams_JPO90} parameterisation, and various slightly diffusive advection schemes.
 For momentum, the main ones are: Laplacian and bilaplacian operators acting along geopotential surfaces,
 and UBS advection schemes when flux form is chosen for the momentum advection.
 …
 the rotation between geopotential and $s$-surfaces,
 while it is only an approximation for the rotation between isoneutral and $z$- or $s$-surfaces.
 Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:PE_iso_tensor} \citep{Cox1987}.
+Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:PE_iso_tensor} \citep{cox_OM87}.
 First, the horizontal contribution of the dianeutral mixing is neglected since the ratio between iso and
 dia-neutral diffusive coefficients is known to be several orders of magnitude smaller than unity.
 …
 \subsubsection{Eddy induced velocity}
 When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used,
+When the \textit{eddy induced velocity} parametrisation (eiv) \citep{gent.mcwilliams_JPO90} is used,
 an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers:
 \[
 …
 \ie on a $f$- or $\beta$-plane, not on the sphere.
 It is also a very good approximation in vicinity of the Equator in
 a geographical coordinate system \citep{Lengaigne_al_JGR03}.
+a geographical coordinate system \citep{lengaigne.madec.ea_JGR03}.
 \subsubsection{lateral bilaplacian momentum diffusive operator}

NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex

-                      r10544
+                      r11123
 In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account.
 These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site.
+These coordinates systems is presented in a report \citep{levier.treguier.ea_rpt07} available on the \NEMO web site.
 \colorbox{yellow}{  end of to be updated}
 …
 which imposes a very small time step when an explicit time stepping is used.
 Two methods are proposed to allow a longer time step for the three-dimensional equations:
 the filtered free surface, which is a modification of the continuous equations %(see \autoref{eq:PE_flt}),
+the filtered free surface, which is a modification of the continuous equations %(see \autoref{eq:PE_flt?}),
 and the split-explicit free surface described below.
 The extra term introduced in the filtered method is calculated implicitly,
 …
 \nlst{namdom}
 %--------------------------------------------------------------------------------------------------------------
 The split-explicit free surface formulation used in OPA follows the one proposed by \citet{Griffies2004}.
+The split-explicit free surface formulation used in OPA follows the one proposed by \citet{Griffies2004?}.
 The general idea is to solve the free surface equation with a small time step,
 while the three dimensional prognostic variables are solved with a longer time step that
 …
       \protect\label{fig:DYN_dynspg_ts}
       Schematic of the split-explicit time stepping scheme for the barotropic and baroclinic modes,
       after \citet{Griffies2004}.
+      after \citet{Griffies2004?}.
       Time increases to the right.
       Baroclinic time steps are denoted by $t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$.
 …
 The split-explicit formulation has a damping effect on external gravity waves,
 which is weaker than the filtered free surface but still significant as shown by \citet{Levier2007} in
+which is weaker than the filtered free surface but still significant as shown by \citet{levier.treguier.ea_rpt07} in
 the case of an analytical barotropic Kelvin wave.
 …
 \label{subsec:DYN_spg_flt}
 The filtered formulation follows the \citet{Roullet2000} implementation.
+The filtered formulation follows the \citet{Roullet2000?} implementation.
 The extra term introduced in the equations (see {\S}I.2.2) is solved implicitly.
 The elliptic solvers available in the code are documented in \autoref{chap:MISC}.
 The amplitude of the extra term is given by the namelist variable \np{rnu}.
 The default value is 1, as recommended by \citet{Roullet2000}
+The default value is 1, as recommended by \citet{Roullet2000?}
 \colorbox{red}{\np{rnu}\forcode{ = 1} to be suppressed from namelist !}
 …
 In the non-linear free surface formulation, the variations of volume are fully taken into account.
 This option is presented in a report \citep{Levier2007} available on the NEMO web site.
+This option is presented in a report \citep{levier.treguier.ea_rpt07} available on the NEMO web site.
 The three time-stepping methods (explicit, split-explicit and filtered) are the same as in
 \autoref{DYN_spg_linear} except that the ocean depth is now time-dependent.

NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex

-                      r10501
+                      r11123
 The time stepping used for processes other than diffusion is the well-known leapfrog scheme
 \citep{Mesinger_Arakawa_Bk76}.
+\citep{mesinger.arakawa_bk76}.
 This scheme is widely used for advection processes in low-viscosity fluids.
 It is a time centred scheme, \ie the RHS in \autoref{eq:STP} is evaluated at time step $t$, the now time step.
 …
 To prevent it, the leapfrog scheme is often used in association with a Robert-Asselin time filter
 (hereafter the LF-RA scheme).
 This filter, first designed by \citet{Robert_JMSJ66} and more comprehensively studied by \citet{Asselin_MWR72},
+This filter, first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72},
 is a kind of laplacian diffusion in time that mixes odd and even time steps:
 \begin{equation}
 …
 $\gamma$ is initialized as \np{rn\_atfp} (namelist parameter).
 Its default value is \np{rn\_atfp}~\forcode{= 10.e-3} (see \autoref{sec:STP_mLF}),
 causing only a weak dissipation of high frequency motions (\citep{Farge1987}).
+causing only a weak dissipation of high frequency motions (\citep{farge-coulombier_phd87}).
 The addition of a time filter degrades the accuracy of the calculation from second to first order.
 However, the second order truncation error is proportional to $\gamma$, which is small compared to 1.
 …
 This is diffusive in time and conditionally stable.
 The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies_Bk04}:
+The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}:
 \begin{equation}
   \label{eq:STP_euler_stability}
 …
 $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms,
 therefore a special adaptation of the Gauss elimination procedure is used to find the solution
 (see for example \citet{Richtmyer1967}).
+(see for example \citet{richtmyer.morton_bk67}).
 % -------------------------------------------------------------------------------------------------------------
 …
     \caption{
       \protect\label{fig:TimeStep_flowchart}
       Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}.
+      Sketch of the leapfrog time stepping sequence in \NEMO from \citet{leclair.madec_OM09}.
       The use of a semi -implicit computation of the hydrostatic pressure gradient requires the tracer equation to
       be stepped forward prior to the momentum equation.
 …
 \label{sec:STP_mLF}
 Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to
+Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to
 ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter.
 The modifications affect both the forcing and filtering treatments in the LF-RA scheme.
 …
 The change in the forcing formulation given by \autoref{eq:STP_forcing} (see \autoref{fig:MLF_forcing})
 has a significant effect:
 the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}.
+the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}.
 % forcing seen by the model....
 This property improves the LF-RA scheme in two respects.
 …
 (last term in \autoref{eq:STP_RA} compared to \autoref{eq:STP_asselin}).
 Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme,
 the modified formulation becomes conservative \citep{Leclair_Madec_OM09}.
+the modified formulation becomes conservative \citep{leclair.madec_OM09}.
 Second, the LF-RA becomes a truly quasi -second order scheme.
 Indeed, \autoref{eq:STP_forcing} used in combination with a careful treatment of static instability

NEMO/trunk/doc/latex/NEMO/subfiles/introduction.tex

-                      r10544
+                      r11123
 The ocean component of \NEMO has been developed from the legacy of the OPA model, release 8.2,
 described in \citet{Madec1998}.
+described in \citet{madec.delecluse.ea_NPM98}.
 This model has been used for a wide range of applications, both regional or global, as a forced ocean model and
 as a model coupled with the sea-ice and/or the atmosphere.
 …
 Within the \NEMO system the ocean model is interactively coupled with a sea ice model (SI$^3$) and
 a biogeochemistry model (PISCES).
 Interactive coupling to Atmospheric models is possible via the OASIS coupler \citep{OASIS2006}.
+Interactive coupling to Atmospheric models is possible via the \href{https://portal.enes.org/oasis}{OASIS coupler}.
 Two-way nesting is also available through an interface to the AGRIF package
 (Adaptative Grid Refinement in \fortran) \citep{Debreu_al_CG2008}.
+(Adaptative Grid Refinement in \fortran) \citep{debreu.vouland.ea_CG08}.
 % Needs to be reviewed
 %The interface code for coupling to an alternative sea ice model (CICE, \citet{Hunke2008}) has now been upgraded so
 …
 The lateral Laplacian and biharmonic viscosity and diffusion can be rotated following
 a geopotential or neutral direction.
 There is an optional eddy induced velocity \citep{Gent1990} with a space and time variable coefficient
 \citet{Treguier1997}.
+There is an optional eddy induced velocity \citep{gent.mcwilliams_JPO90} with a space and time variable coefficient
+\citet{treguier.held.ea_JPO97}.
 The model has vertical harmonic viscosity and diffusion with a space and time variable coefficient,
 with options to compute the coefficients with \citet{Blanke1993}, \citet{Pacanowski_Philander_JPO81}, or
 \citet{Umlauf_Burchard_JMS03} mixing schemes.
+with options to compute the coefficients with \citet{blanke.delecluse_JPO93}, \citet{pacanowski.philander_JPO81}, or
+\citet{umlauf.burchard_JMR03} mixing schemes.
 %%gm    To be put somewhere else ....
 …
 NEMO/OPA, like all research tools, is in perpetual evolution.
 The present document describes the OPA version include in the release 3.4 of NEMO.
 This release differs significantly from version 8, documented in \citet{Madec1998}. \\
+This release differs significantly from version 8, documented in \citet{madec.delecluse.ea_NPM98}. \\
 The main modifications from OPA v8 and NEMO/OPA v3.2 are :
 …
 \item
   introduction of partial step representation of bottom topography
   \citep{Barnier_al_OD06, Le_Sommer_al_OM09, Penduff_al_OS07};
+  \citep{barnier.madec.ea_OD06, le-sommer.penduff.ea_OM09, penduff.le-sommer.ea_OS07};
 \item
   partial reactivation of a terrain-following vertical coordinate ($s$- and hybrid $s$-$z$) with
 …
   additional advection schemes for tracers;
 \item
   implementation of the AGRIF package (Adaptative Grid Refinement in \fortran) \citep{Debreu_al_CG2008};
+  implementation of the AGRIF package (Adaptative Grid Refinement in \fortran) \citep{debreu.vouland.ea_CG08};
 \item
   online diagnostics : tracers trend in the mixed layer and vorticity balance;
 …
   RGB light penetration and optional use of ocean color
 \item
   major changes in the TKE schemes: it now includes a Langmuir cell parameterization \citep{Axell_JGR02},
   the \citet{Mellor_Blumberg_JPO04} surface wave breaking parameterization, and has a time discretization which
   is energetically consistent with the ocean model equations \citep{Burchard_OM02, Marsaleix_al_OM08};
+  major changes in the TKE schemes: it now includes a Langmuir cell parameterization \citep{axell_JGR02},
+  the \citet{mellor.blumberg_JPO04} surface wave breaking parameterization, and has a time discretization which
+  is energetically consistent with the ocean model equations \citep{burchard_OM02, marsaleix.auclair.ea_OM08};
 \item
   tidal mixing parametrisation (bottom intensification) + Indonesian specific tidal mixing
   \citep{Koch-Larrouy_al_GRL07};
+  \citep{koch-larrouy.madec.ea_GRL07};
 \item
   introduction of LIM-3, the new Louvain-la-Neuve sea-ice model
   (C-grid rheology and new thermodynamics including bulk ice salinity)
   \citep{Vancoppenolle_al_OM09a, Vancoppenolle_al_OM09b}
+  \citep{vancoppenolle.fichefet.ea_OM09*a, vancoppenolle.fichefet.ea_OM09*b}
 \end{itemize}
 …
 \item
   introduction of a modified leapfrog-Asselin filter time stepping scheme
   \citep{Leclair_Madec_OM09};
 \item
   additional scheme for iso-neutral mixing \citep{Griffies_al_JPO98}, although it is still a "work in progress";
 \item
   a rewriting of the bottom boundary layer scheme, following \citet{Campin_Goosse_Tel99};
 \item
   addition of a Generic Length Scale vertical mixing scheme, following \citet{Umlauf_Burchard_JMS03};
+  \citep{leclair.madec_OM09};
+\item
+  additional scheme for iso-neutral mixing \citep{griffies.gnanadesikan.ea_JPO98}, although it is still a "work in progress";
+\item
+  a rewriting of the bottom boundary layer scheme, following \citet{campin.goosse_T99};
+\item
+  addition of a Generic Length Scale vertical mixing scheme, following \citet{umlauf.burchard_JMR03};
 \item
   addition of the atmospheric pressure as an external forcing on both ocean and sea-ice dynamics;
 \item
   addition of a diurnal cycle on solar radiation \citep{Bernie_al_CD07};
+  addition of a diurnal cycle on solar radiation \citep{bernie.guilyardi.ea_CD07};
 \item
   river runoffs added through a non-zero depth, and having its own temperature and salinity;
 …
   coupling interface adjusted for WRF atmospheric model;
 \item
   C-grid ice rheology now available fro both LIM-2 and LIM-3 \citep{Bouillon_al_OM09};
+  C-grid ice rheology now available fro both LIM-2 and LIM-3 \citep{bouillon.maqueda.ea_OM09};
 \item
   LIM-3 ice-ocean momentum coupling applied to LIM-2;
 …
 \begin{itemize}
 \item finalisation of above iso-neutral mixing \citep{Griffies_al_JPO98}";
+\item finalisation of above iso-neutral mixing \citep{griffies.gnanadesikan.ea_JPO98}";
 \item "Neptune effect" parametrisation;
 \item horizontal pressure gradient suitable for s-coordinate;

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

Legend:

NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.sty

NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.tex

NEMO/trunk/doc/latex/NEMO/subfiles

NEMO/trunk/doc/latex/NEMO/subfiles/annex_A.tex

NEMO/trunk/doc/latex/NEMO/subfiles/annex_B.tex

NEMO/trunk/doc/latex/NEMO/subfiles/annex_D.tex

NEMO/trunk/doc/latex/NEMO/subfiles/annex_E.tex

NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_CONFIG.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIU.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_STO.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_misc.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex

NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex

NEMO/trunk/doc/latex/NEMO/subfiles/introduction.tex

Download in other formats: