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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9407 – NEMO

Changeset 9407


Ignore:
Timestamp:
2018-03-15T17:40:35+01:00 (5 years ago)
Author:
nicolasmartin
Message:

Complete refactoring of cross-referencing

  • Use of \autoref instead of simple \ref for contextual text depending on target type
  • creation of few prefixes for marker to identify the type reference: apdx|chap|eq|fig|sec|subsec|tab
Location:
branches/2017/dev_merge_2017/DOC
Files:
3 added
25 edited
5 moved

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_merge_2017/DOC/HTML_htlatex.sh

    r9394 r9407  
    11#!/bin/bash 
    22 
    3 latex     -shell-escape NEMO_manual 
    4 makeindex               NEMO_manual 
    5 bibtex                  NEMO_manual 
    6 latex     -shell-escape NEMO_manual 
     3./inc/clean.sh 
     4./inc/build.sh 
    75 
    86mkdir -p html_htlatex 
    9 htlatex NEMO_manual "htlatex,2" "" "-dhtml_htlatex/" "-shell-escape" 
     7cd tex_main 
     8htlatex NEMO_manual "NEMO_manual,2" "" "-d../html_htlatex/" "-shell-escape" 
     9cd - 
    1010 
    1111exit 0 

  • branches/2017/dev_merge_2017/DOC/HTML_latex2html.sh

    r9394 r9407  
    11#!/bin/bash 
    22 
     3./inc/clean.sh 
     4./inc/build.sh 
     5 
     6cd tex_main 
    37sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s/^/%/' \ 
    4    texfiles/chapters/*.tex 
    5 sed -i    '30,${s#\\subfile{#\\include{#g}' \ 
     8   ../tex_sub/*.tex 
     9sed -i    's#\\subfile{#\\include{#g' \ 
    610   NEMO_manual.tex 
    7  
    8 latex     -shell-escape    NEMO_manual 
    9 makeindex                  NEMO_manual 
    10 bibtex                     NEMO_manual 
    11  
    12 latex2html -local_icons -no_footnode -split 4 -link 2 -mkdir -dir html_LaTeX2HTML   \ 
    13             $*                                                                      \ 
     11latex2html -local_icons -no_footnode -split 4 -link 2 -mkdir -dir ../html_LaTeX2HTML   \ 
     12            $*                                                                         \ 
    1413   NEMO_manual 
    15  
    1614sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s/^%//' \ 
    17    texfiles/chapters/*.tex 
    18 sed -i    '30,${s#\\include{#\\subfile{#g}' \ 
     15   ../tex_sub/*.tex 
     16sed -i    's#\\include{#\\subfile{#g' \ 
    1917   NEMO_manual.tex 
     18cd - 
    2019 
    2120exit 0 

  • branches/2017/dev_merge_2017/DOC/inc/build.sh

    r9394 r9407  
    1111latex       ${latex_opts}           ${latex_file} 
    1212 
    13 pdflatex    ${latex_opts}           ${latex_file} 
    14  
    15 mv ${latex_file}.pdf .. 
    16  
    1713cd - 
    1814 

  • branches/2017/dev_merge_2017/DOC/inc/clean.sh

    r9394 r9407  
    11#!/bin/bash 
    22 
    3 rm -f $( ls -1 tex_main/NEMO_* | egrep -v "\.(bib|ist|sty|tex)$" ) 
     3rm -f $( ls -1 tex_main/NEMO_* | egrep -v "\.(bib|cfg|ist|sty|tex)$" ) 
    44#rm -rf _minted-* 
    55#rm -rf html* 

  • branches/2017/dev_merge_2017/DOC/tex_main/NEMO_manual.cfg

    r9394 r9407  
    1 \Preamble{html} 
     1\Preamble{xhtml,mathml} 
     2 
     3\Configure{@HEAD}{% 
     4\HCode{<script type="text/javascript" 
     5   src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=MML_CHTML"> 
     6 </script>\Hnewline}} 
    27 
    38\begin{document} 

  • branches/2017/dev_merge_2017/DOC/tex_sub/annex_A.tex

    r9393 r9407  
    66% ================================================================ 
    77\chapter{Curvilinear $s-$Coordinate Equations} 
    8 \label{Apdx_A} 
     8\label{apdx:A} 
    99\minitoc 
    1010 
    … …  
    1616% ================================================================ 
    1717\section{Chain rule for $s-$coordinates} 
    18 \label{Apdx_A_continuity} 
     18\label{sec:A_continuity} 
    1919 
    2020In order to establish the set of Primitive Equation in curvilinear $s$-coordinates 
    2121($i.e.$ an orthogonal curvilinear coordinate in the horizontal and an Arbitrary Lagrangian  
    2222Eulerian (ALE) coordinate in the vertical), we start from the set of equations established  
    23 in \S\ref{PE_zco_Eq} for the special case $k = z$ and thus $e_3 = 1$, and we introduce  
     23in \autoref{subsec:PE_zco_Eq} for the special case $k = z$ and thus $e_3 = 1$, and we introduce  
    2424an arbitrary vertical coordinate $a = a(i,j,z,t)$. Let us define a new vertical scale factor by  
    2525$e_3 = \partial z / \partial s$ (which now depends on $(i,j,z,t)$) and the horizontal  
    2626slope of $s-$surfaces by : 
    27 \begin{equation} \label{Apdx_A_s_slope} 
     27\begin{equation} \label{apdx:A_s_slope} 
    2828\sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s  
    2929\quad \text{and} \quad  
    … …  
    3333The chain rule to establish the model equations in the curvilinear $s-$coordinate  
    3434system is: 
    35 \begin{equation} \label{Apdx_A_s_chain_rule} 
     35\begin{equation} \label{apdx:A_s_chain_rule} 
    3636\begin{aligned} 
    3737&\left. {\frac{\partial \bullet }{\partial t}} \right|_z  = 
    … …  
    5454In particular applying the time derivative chain rule to $z$ provides the expression  
    5555for $w_s$,  the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate: 
    56 \begin{equation} \label{Apdx_A_w_in_s} 
     56\begin{equation} \label{apdx:A_w_in_s} 
    5757w_s   =  \left.   \frac{\partial z }{\partial t}   \right|_s  
    5858            = \frac{\partial z}{\partial s} \; \frac{\partial s}{\partial t}  
    … …  
    6565% ================================================================ 
    6666\section{Continuity equation in $s-$coordinates} 
    67 \label{Apdx_A_continuity} 
    68  
    69 Using (\ref{Apdx_A_s_chain_rule}) and the fact that the horizontal scale factors  
     67\label{sec:A_continuity} 
     68 
     69Using (\autoref{apdx:A_s_chain_rule}) and the fact that the horizontal scale factors  
    7070$e_1$ and $e_2$ do not depend on the vertical coordinate, the divergence of  
    7171the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows 
    … …  
    131131Introducing the dia-surface velocity component, $\omega $, defined as  
    132132the volume flux across the moving $s$-surfaces per unit horizontal area: 
    133 \begin{equation} \label{Apdx_A_w_s} 
     133\begin{equation} \label{apdx:A_w_s} 
    134134\omega  = w - w_s - \sigma _1 \,u - \sigma _2 \,v    \\ 
    135135\end{equation} 
    136 with $w_s$ given by \eqref{Apdx_A_w_in_s}, we obtain the expression for  
     136with $w_s$ given by \autoref{apdx:A_w_in_s}, we obtain the expression for  
    137137the divergence of the velocity in the curvilinear $s-$coordinate system: 
    138138\begin{subequations}  
    … …  
    167167\end{subequations} 
    168168 
    169 As a result, the continuity equation \eqref{Eq_PE_continuity} in the  
     169As a result, the continuity equation \autoref{eq:PE_continuity} in the  
    170170$s-$coordinates is: 
    171 \begin{equation} \label{Apdx_A_sco_Continuity} 
     171\begin{equation} \label{apdx:A_sco_Continuity} 
    172172\frac{1}{e_3 } \frac{\partial e_3}{\partial t}  
    173173+ \frac{1}{e_1 \,e_2 \,e_3 }\left[  
    … …  
    184184% ================================================================ 
    185185\section{Momentum equation in $s-$coordinate} 
    186 \label{Apdx_A_momentum} 
     186\label{sec:A_momentum} 
    187187 
    188188Here we only consider the first component of the momentum equation,  
    … …  
    193193$\bullet$ \textbf{Total derivative in vector invariant form} 
    194194 
    195 Let us consider \eqref{Eq_PE_dyn_vect}, the first component of the momentum  
     195Let us consider \autoref{eq:PE_dyn_vect}, the first component of the momentum  
    196196equation in the vector invariant form. Its total $z-$coordinate time derivative,  
    197197$\left. \frac{D u}{D t} \right|_z$ can be transformed as follows in order to obtain  
    … …  
    258258\end{subequations} 
    259259% 
    260 Applying the time derivative chain rule (first equation of (\ref{Apdx_A_s_chain_rule})) 
    261 to $u$ and using (\ref{Apdx_A_w_in_s}) provides the expression of the last term  
     260Applying the time derivative chain rule (first equation of (\autoref{apdx:A_s_chain_rule})) 
     261to $u$ and using (\autoref{apdx:A_w_in_s}) provides the expression of the last term  
    262262of the right hand side, 
    263263\begin{equation*} {\begin{array}{*{20}l}  
    … …  
    269269leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative,  
    270270$i.e.$ the total $s-$coordinate time derivative : 
    271 \begin{align} \label{Apdx_A_sco_Dt_vect} 
     271\begin{align} \label{apdx:A_sco_Dt_vect} 
    272272\left. \frac{D u}{D t} \right|_s  
    273273  = \left. {\frac{\partial u }{\partial t}} \right|_s        
    … …  
    285285 
    286286Let us start from the total time derivative in the curvilinear $s-$coordinate system  
    287 we have just establish. Following the procedure used to establish (\ref{Eq_PE_flux_form}),  
     287we have just establish. Following the procedure used to establish (\autoref{eq:PE_flux_form}),  
    288288it can be transformed into : 
    289289%\begin{subequations}  
    … …  
    356356which leads to the $s-$coordinate flux formulation of the total $s-$coordinate time derivative,  
    357357$i.e.$ the total $s-$coordinate time derivative in flux form : 
    358 \begin{flalign}\label{Apdx_A_sco_Dt_flux} 
     358\begin{flalign}\label{apdx:A_sco_Dt_flux} 
    359359\left. \frac{D u}{D t} \right|_s   = \frac{1}{e_3}  \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s   
    360360           + \left.  \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right)    \right|_s 
    … …  
    365365It has the same form as in the $z-$coordinate but for the vertical scale factor  
    366366that has appeared inside the time derivative which comes from the modification  
    367 of (\ref{Apdx_A_sco_Continuity}), the continuity equation. 
     367of (\autoref{apdx:A_sco_Continuity}), the continuity equation. 
    368368 
    369369$\ $\newline    % force a new ligne 
    … …  
    381381\end{equation*} 
    382382Applying similar manipulation to the second component and replacing  
    383 $\sigma _1$ and $\sigma _2$ by their expression \eqref{Apdx_A_s_slope}, it comes: 
    384 \begin{equation} \label{Apdx_A_grad_p} 
     383$\sigma _1$ and $\sigma _2$ by their expression \autoref{apdx:A_s_slope}, it comes: 
     384\begin{equation} \label{apdx:A_grad_p} 
    385385\begin{split} 
    386386 -\frac{1}{\rho _o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z 
    … …  
    394394\end{equation} 
    395395 
    396 An additional term appears in (\ref{Apdx_A_grad_p}) which accounts for the  
     396An additional term appears in (\autoref{apdx:A_grad_p}) which accounts for the  
    397397tilt of $s-$surfaces with respect to geopotential $z-$surfaces. 
    398398 
    … …  
    408408\end{equation*} 
    409409Therefore, $p$ and $p_h'$ are linked through: 
    410 \begin{equation} \label{Apdx_A_pressure} 
     410\begin{equation} \label{apdx:A_pressure} 
    411411   p = \rho_o \; p_h' + g \, ( z + \eta ) 
    412412\end{equation} 
    … …  
    416416\end{equation*} 
    417417 
    418 Substituing \eqref{Apdx_A_pressure} in \eqref{Apdx_A_grad_p} and using the definition of  
     418Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p} and using the definition of  
    419419the density anomaly it comes the expression in two parts: 
    420 \begin{equation} \label{Apdx_A_grad_p} 
     420\begin{equation} \label{apdx:A_grad_p} 
    421421\begin{split} 
    422422 -\frac{1}{\rho _o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z 
    … …  
    430430\end{equation} 
    431431This formulation of the pressure gradient is characterised by the appearance of a term depending on the  
    432 the sea surface height only (last term on the right hand side of expression \eqref{Apdx_A_grad_p}). 
     432the sea surface height only (last term on the right hand side of expression \autoref{apdx:A_grad_p}). 
    433433This term will be loosely termed \textit{surface pressure gradient} 
    434434whereas the first term will be termed the  
    … …  
    445445The coriolis and forcing terms as well as the the vertical physics remain unchanged  
    446446as they involve neither time nor space derivatives. The form of the lateral physics is  
    447 discussed in appendix~\ref{Apdx_B}. 
     447discussed in \autoref{apdx:B}. 
    448448 
    449449 
    … …  
    455455solved by the model has the same mathematical expression as the one in a curvilinear  
    456456$z-$coordinate, except for the pressure gradient term : 
    457 \begin{subequations} \label{Apdx_A_dyn_vect} 
    458 \begin{multline} \label{Apdx_A_PE_dyn_vect_u} 
     457\begin{subequations} \label{apdx:A_dyn_vect} 
     458\begin{multline} \label{apdx:A_PE_dyn_vect_u} 
    459459 \frac{\partial u}{\partial t}= 
    460460   +   \left( {\zeta +f} \right)\,v                                     
    … …  
    465465   +   D_u^{\vect{U}}  +   F_u^{\vect{U}} 
    466466\end{multline} 
    467 \begin{multline} \label{Apdx_A_dyn_vect_v} 
     467\begin{multline} \label{apdx:A_dyn_vect_v} 
    468468\frac{\partial v}{\partial t}= 
    469469   -   \left( {\zeta +f} \right)\,u    
    … …  
    477477whereas the flux form momentum equation differ from it by the formulation of both 
    478478the time derivative and the pressure gradient term  : 
    479 \begin{subequations} \label{Apdx_A_dyn_flux} 
    480 \begin{multline} \label{Apdx_A_PE_dyn_flux_u} 
     479\begin{subequations} \label{apdx:A_dyn_flux} 
     480\begin{multline} \label{apdx:A_PE_dyn_flux_u} 
    481481 \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t} = 
    482482   \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right)  
    … …  
    487487   +   D_u^{\vect{U}}  +   F_u^{\vect{U}} 
    488488\end{multline} 
    489 \begin{multline} \label{Apdx_A_dyn_flux_v} 
     489\begin{multline} \label{apdx:A_dyn_flux_v} 
    490490 \frac{1}{e_3}\frac{\partial \left(  e_3\,v  \right) }{\partial t}= 
    491491   -  \nabla \cdot \left(   {{\rm {\bf U}}\,v}   \right)  
    … …  
    499499Both formulation share the same hydrostatic pressure balance expressed in terms of 
    500500hydrostatic pressure and density anomalies, $p_h'$ and $d=( \frac{\rho}{\rho_o}-1 )$: 
    501 \begin{equation} \label{Apdx_A_dyn_zph} 
     501\begin{equation} \label{apdx:A_dyn_zph} 
    502502\frac{\partial p_h'}{\partial k} = - d \, g \, e_3 
    503503\end{equation} 
    … …  
    516516% ================================================================ 
    517517\section{Tracer equation} 
    518 \label{Apdx_A_tracer} 
     518\label{sec:A_tracer} 
    519519 
    520520The tracer equation is obtained using the same calculation as for the continuity  
    521521equation and then regrouping the time derivative terms in the left hand side : 
    522522 
    523 \begin{multline} \label{Apdx_A_tracer} 
     523\begin{multline} \label{apdx:A_tracer} 
    524524 \frac{1}{e_3} \frac{\partial \left(  e_3 T  \right)}{\partial t}  
    525525   = -\frac{1}{e_1 \,e_2 \,e_3}  

  • branches/2017/dev_merge_2017/DOC/tex_sub/annex_B.tex

    r9393 r9407  
    55% ================================================================ 
    66\chapter{Appendix B : Diffusive Operators} 
    7 \label{Apdx_B} 
     7\label{apdx:B} 
    88\minitoc 
    99 
    … …  
    1616% ================================================================ 
    1717\section{Horizontal/Vertical $2^{nd}$ order tracer diffusive operators} 
    18 \label{Apdx_B_1} 
     18\label{sec:B_1} 
    1919 
    2020\subsubsection*{In z-coordinates} 
    2121In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator 
    2222is given by: 
    23 \begin{eqnarray} \label{Apdx_B1} 
     23\begin{eqnarray} \label{apdx:B1} 
    2424 &D^T = \frac{1}{e_1 \, e_2}      \left[ 
    2525  \left. \frac{\partial}{\partial i} \left(  \frac{e_2}{e_1}A^{lT} \;\left. \frac{\partial T}{\partial i} \right|_z   \right)   \right|_z      \right. 
    … …  
    3131\subsubsection*{In generalized vertical coordinates} 
    3232In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and 
    33 $\sigma_2$ by \eqref{Apdx_A_s_slope} and the vertical/horizontal ratio of diffusion 
     33$\sigma_2$ by \autoref{apdx:A_s_slope} and the vertical/horizontal ratio of diffusion 
    3434coefficient by $\epsilon = A^{vT} / A^{lT}$. The diffusion operator is given by: 
    3535 
    36 \begin{equation} \label{Apdx_B2} 
     36\begin{equation} \label{apdx:B2} 
    3737D^T = \left. \nabla \right|_s \cdot 
    3838           \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\ 
    … …  
    5656\end{subequations} 
    5757 
    58 Equation \eqref{Apdx_B2} is obtained from \eqref{Apdx_B1} without any 
     58Equation \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any 
    5959additional assumption. Indeed, for the special case $k=z$ and thus $e_3 =1$, 
    60 we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in Appendix~\ref{Apdx_A} 
    61 and use \eqref{Apdx_A_s_slope} and \eqref{Apdx_A_s_chain_rule}. 
     60we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} 
     61and use \autoref{apdx:A_s_slope} and \autoref{apdx:A_s_chain_rule}. 
    6262Since no cross horizontal derivative $\partial _i \partial _j $ appears in 
    63 \eqref{Apdx_B1}, the ($i$,$z$) and ($j$,$z$) planes are independent. 
     63\autoref{apdx:B1}, the ($i$,$z$) and ($j$,$z$) planes are independent. 
    6464The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$) 
    6565transformation without any loss of generality: 
    … …  
    139139% ================================================================ 
    140140\section{Iso/Diapycnal $2^{nd}$ order tracer diffusive operators} 
    141 \label{Apdx_B_2} 
     141\label{sec:B_2} 
    142142 
    143143\subsubsection*{In z-coordinates} 
    … …  
    147147formulated, takes the following form \citep{Redi_JPO82}: 
    148148 
    149 \begin{equation} \label{Apdx_B3} 
     149\begin{equation} \label{apdx:B3} 
    150150\textbf {A}_{\textbf I} = \frac{A^{lT}}{\left( {1+a_1 ^2+a_2 ^2} \right)} 
    151151\left[ {{\begin{array}{*{20}c} 
    … …  
    166166In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean, so 
    167167$\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}: 
    168 \begin{subequations} \label{Apdx_B4} 
    169 \begin{equation} \label{Apdx_B4a} 
     168\begin{subequations} \label{apdx:B4} 
     169\begin{equation} \label{apdx:B4a} 
    170170{\textbf{A}_{\textbf{I}}} \approx A^{lT}\;\Re\;\text{where} \;\Re = 
    171171\left[ {{\begin{array}{*{20}c} 
    … …  
    176176\end{equation} 
    177177and the iso/dianeutral diffusive operator in $z$-coordinates is then 
    178 \begin{equation}\label{Apdx_B4b} 
     178\begin{equation}\label{apdx:B4b} 
    179179 D^T = \left. \nabla \right|_z \cdot 
    180180           \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_z T  \right]. \\ 
    … …  
    183183 
    184184 
    185 Physically, the full tensor \eqref{Apdx_B3} 
     185Physically, the full tensor \autoref{apdx:B3} 
    186186represents strong isoneutral diffusion on a plane parallel to the isoneutral 
    187187surface and weak dianeutral diffusion perpendicular to this plane. 
    188 However, the approximate `weak-slope' tensor \eqref{Apdx_B4a} represents strong 
     188However, the approximate `weak-slope' tensor \autoref{apdx:B4a} represents strong 
    189189diffusion along the isoneutral surface, with weak 
    190190\emph{vertical}  diffusion -- the principal axes of the tensor are no 
    191191longer orthogonal. This simplification also decouples 
    192192the ($i$,$z$) and ($j$,$z$) planes of the tensor. The weak-slope operator therefore takes the same 
    193 form, \eqref{Apdx_B4}, as \eqref{Apdx_B2}, the diffusion operator for geopotential 
     193form, \autoref{apdx:B4}, as \autoref{apdx:B2}, the diffusion operator for geopotential 
    194194diffusion written in non-orthogonal $i,j,s$-coordinates. Written out 
    195195explicitly, 
    196196 
    197 \begin{multline} \label{Apdx_B_ldfiso} 
     197\begin{multline} \label{apdx:B_ldfiso} 
    198198 D^T=\frac{1}{e_1 e_2 }\left\{ 
    199199 {\;\frac{\partial }{\partial i}\left[ {A_h \left( {\frac{e_2}{e_1}\frac{\partial T}{\partial i}-a_1 \frac{e_2}{e_3}\frac{\partial T}{\partial k}} \right)} \right]} 
    … …  
    203203 
    204204 
    205 The isopycnal diffusion operator \eqref{Apdx_B4}, 
    206 \eqref{Apdx_B_ldfiso} conserves tracer quantity and dissipates its 
    207 square. The demonstration of the first property is trivial as \eqref{Apdx_B4} is the divergence 
     205The isopycnal diffusion operator \autoref{apdx:B4}, 
     206\autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its 
     207square. The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence 
    208208of fluxes. Let us demonstrate the second one: 
    209209\begin{equation*} 
    … …  
    233233\subsubsection*{In generalized vertical coordinates} 
    234234 
    235 Because the weak-slope operator \eqref{Apdx_B4}, \eqref{Apdx_B_ldfiso} is decoupled 
     235Because the weak-slope operator \autoref{apdx:B4}, \autoref{apdx:B_ldfiso} is decoupled 
    236236in the ($i$,$z$) and ($j$,$z$) planes, it may be transformed into 
    237 generalized $s$-coordinates in the same way as \eqref{Apdx_B_1} was transformed into 
    238 \eqref{Apdx_B_2}. The resulting operator then takes the simple form 
    239  
    240 \begin{equation} \label{Apdx_B_ldfiso_s} 
     237generalized $s$-coordinates in the same way as \autoref{sec:B_1} was transformed into 
     238\autoref{sec:B_2}. The resulting operator then takes the simple form 
     239 
     240\begin{equation} \label{apdx:B_ldfiso_s} 
    241241D^T = \left. \nabla \right|_s \cdot 
    242242           \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\ 
    … …  
    258258\end{equation*} 
    259259 
    260 To prove  \eqref{Apdx_B5}  by direct re-expression of \eqref{Apdx_B_ldfiso} is 
     260To prove  \autoref{apdx:B5}  by direct re-expression of \autoref{apdx:B_ldfiso} is 
    261261straightforward, but laborious. An easier way is first to note (by reversing the 
    262 derivation of \eqref{Apdx_B_2} from \eqref{Apdx_B_1} ) that the 
     262derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that the 
    263263weak-slope operator may be \emph{exactly} reexpressed in  
    264264non-orthogonal $i,j,\rho$-coordinates as 
    265265 
    266 \begin{equation} \label{Apdx_B5} 
     266\begin{equation} \label{apdx:B5} 
    267267D^T = \left. \nabla \right|_\rho \cdot 
    268268           \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_\rho T  \right] \\ 
    … …  
    274274\end{equation} 
    275275Then direct transformation from $i,j,\rho$-coordinates to 
    276 $i,j,s$-coordinates gives \eqref{Apdx_B_ldfiso_s} immediately. 
     276$i,j,s$-coordinates gives \autoref{apdx:B_ldfiso_s} immediately. 
    277277 
    278278Note that the weak-slope approximation is only made in 
    … …  
    282282the  $s$-surfaces, in the same way as the transformation of 
    283283horizontal/vertical Laplacian diffusion in $z$-coordinates, 
    284 \eqref{Apdx_B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces. 
     284\autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces. 
    285285 
    286286 
    … …  
    289289% ================================================================ 
    290290\section{Lateral/Vertical momentum diffusive operators} 
    291 \label{Apdx_B_3} 
     291\label{sec:B_3} 
    292292 
    293293The second order momentum diffusion operator (Laplacian) in the $z$-coordinate 
    294 is found by applying \eqref{Eq_PE_lap_vector}, the expression for the Laplacian 
     294is found by applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian 
    295295of a vector,  to the horizontal velocity vector : 
    296296\begin{align*} 
    … …  
    329329\end{array} }} \right) 
    330330\end{align*} 
    331 Using \eqref{Eq_PE_div}, the definition of the horizontal divergence, the third 
     331Using \autoref{eq:PE_div}, the definition of the horizontal divergence, the third 
    332332componant of the second vector is obviously zero and thus : 
    333333\begin{equation*} 
    … …  
    336336 
    337337Note that this operator ensures a full separation between the vorticity and horizontal 
    338 divergence fields (see Appendix~\ref{Apdx_C}). It is only equal to a Laplacian 
     338divergence fields (see \autoref{apdx:C}). It is only equal to a Laplacian 
    339339applied to each component in Cartesian coordinates, not on the sphere. 
    340340 
    341341The horizontal/vertical second order (Laplacian type) operator used to diffuse 
    342342horizontal momentum in the $z$-coordinate therefore takes the following form : 
    343 \begin{equation} \label{Apdx_B_Lap_U} 
     343\begin{equation} \label{apdx:B_Lap_U} 
    344344 {\textbf{D}}^{\textbf{U}} = 
    345345     \nabla _h \left( {A^{lm}\;\chi } \right) 
    … …  
    360360\end{align*} 
    361361 
    362 Note Bene: introducing a rotation in \eqref{Apdx_B_Lap_U} does not lead to a 
     362Note Bene: introducing a rotation in \autoref{apdx:B_Lap_U} does not lead to a 
    363363useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate. 
    364364Similarly, we did not found an expression of practical use for the geopotential 
    365365horizontal/vertical Laplacian operator in the $s$-coordinate. Generally, 
    366 \eqref{Apdx_B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is 
     366\autoref{apdx:B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is 
    367367a Laplacian diffusion is applied on momentum along the coordinate directions. 
    368368\end{document} 

  • branches/2017/dev_merge_2017/DOC/tex_sub/annex_C.tex

    r9393 r9407  
    55% ================================================================ 
    66\chapter{Discrete Invariants of the Equations} 
    7 \label{Apdx_C} 
     7\label{apdx:C} 
    88\minitoc 
    99 
    … …  
    2020% ================================================================ 
    2121\section{Introduction / Notations} 
    22 \label{Apdx_C.0} 
     22\label{sec:C.0} 
    2323 
    2424Notation used in this appendix in the demonstations : 
    … …  
    6969\end{flalign*} 
    7070that is in a more compact form : 
    71 \begin{flalign} \label{Eq_Q2_flux} 
     71\begin{flalign} \label{eq:Q2_flux} 
    7272\partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right) 
    7373=&                   \int_D { \frac{Q}{e_3}  \partial_t \left( e_3 \, Q \right) dv }   
    … …  
    8383\end{flalign*} 
    8484that is in a more compact form : 
    85 \begin{flalign} \label{Eq_Q2_vect} 
     85\begin{flalign} \label{eq:Q2_vect} 
    8686\partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right) 
    8787=& \int_D {         Q   \,\partial_t Q  \;dv }   
    … …  
    9494% ================================================================ 
    9595\section{Continuous conservation} 
    96 \label{Apdx_C.1} 
     96\label{sec:C.1} 
    9797 
    9898 
    … …  
    104104Let us first establish those constraint in the continuous world. 
    105105The total energy ($i.e.$ kinetic plus potential energies) is conserved : 
    106 \begin{flalign} \label{Eq_Tot_Energy} 
     106\begin{flalign} \label{eq:Tot_Energy} 
    107107  \partial_t \left( \int_D \left( \frac{1}{2} {\textbf{U}_h}^2 +  \rho \, g \, z \right) \;dv \right)  = & 0 
    108108\end{flalign} 
    … …  
    114114The transformation for the advection term depends on whether  
    115115the vector invariant form or the flux form is used for the momentum equation. 
    116 Using \eqref{Eq_Q2_vect} and introducing \eqref{Apdx_A_dyn_vect} in \eqref{Eq_Tot_Energy}  
     116Using \autoref{eq:Q2_vect} and introducing \autoref{apdx:A_dyn_vect} in \autoref{eq:Tot_Energy}  
    117117for the former form and 
    118 Using \eqref{Eq_Q2_flux} and introducing \eqref{Apdx_A_dyn_flux} in \eqref{Eq_Tot_Energy}  
     118Using \autoref{eq:Q2_flux} and introducing \autoref{apdx:A_dyn_flux} in \autoref{eq:Tot_Energy}  
    119119for the latter form  leads to: 
    120120 
    121 \begin{subequations} \label{E_tot} 
     121\begin{subequations} \label{eq:E_tot} 
    122122 
    123123advection term (vector invariant form): 
    124 \begin{equation} \label{E_tot_vect_vor} 
     124\begin{equation} \label{eq:E_tot_vect_vor} 
    125125\int\limits_D  \zeta \; \left( \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0   \\ 
    126126\end{equation} 
    127127% 
    128 \begin{equation} \label{E_tot_vect_adv} 
     128\begin{equation} \label{eq:E_tot_vect_adv} 
    129129   \int\limits_D  \textbf{U}_h \cdot \nabla_h \left( \frac{{\textbf{U}_h}^2}{2} \right)     dv  
    130130+ \int\limits_D  \textbf{U}_h \cdot \nabla_z \textbf{U}_h  \;dv     
    … …  
    133133 
    134134advection term (flux form): 
    135 \begin{equation} \label{E_tot_flux_metric} 
     135\begin{equation} \label{eq:E_tot_flux_metric} 
    136136\int\limits_D  \frac{1} {e_1 e_2 } \left( v \,\partial_i e_2 - u \,\partial_j e_1  \right)\;  
    137137 \left(  \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0   \\ 
    138138\end{equation} 
    139139 
    140 \begin{equation} \label{E_tot_flux_adv} 
     140\begin{equation} \label{eq:E_tot_flux_adv} 
    141141   \int\limits_D \textbf{U}_h \cdot     \left(                 {{\begin{array} {*{20}c} 
    142142\nabla \cdot \left( \textbf{U}\,u \right) \hfill \\ 
    … …  
    146146 
    147147coriolis term 
    148 \begin{equation} \label{E_tot_cor} 
     148\begin{equation} \label{eq:E_tot_cor} 
    149149\int\limits_D  f   \; \left( \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0   \\ 
    150150\end{equation} 
    151151 
    152152pressure gradient: 
    153 \begin{equation} \label{E_tot_pg} 
     153\begin{equation} \label{eq:E_tot_pg} 
    154154   - \int\limits_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv  
    155155= - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv 
    … …  
    171171 
    172172Vector invariant form: 
    173 \begin{subequations} \label{E_tot_vect} 
    174 \begin{equation} \label{E_tot_vect_vor} 
     173\begin{subequations} \label{eq:E_tot_vect} 
     174\begin{equation} \label{eq:E_tot_vect_vor} 
    175175\int\limits_D   \textbf{U}_h \cdot \text{VOR} \;dv   = 0   \\ 
    176176\end{equation} 
    177 \begin{equation} \label{E_tot_vect_adv} 
     177\begin{equation} \label{eq:E_tot_vect_adv} 
    178178   \int\limits_D  \textbf{U}_h \cdot \text{KEG}  \;dv  
    179179+ \int\limits_D  \textbf{U}_h \cdot \text{ZAD}  \;dv     
    180180-  \int\limits_D { \frac{{\textbf{U}_h}^2}{2} \frac{1}{e_3} \partial_t e_3 \;dv }   = 0   \\ 
    181181\end{equation} 
    182 \begin{equation} \label{E_tot_pg} 
     182\begin{equation} \label{eq:E_tot_pg} 
    183183   - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv  
    184184= - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv 
    … …  
    188188 
    189189Flux form: 
    190 \begin{subequations} \label{E_tot_flux} 
    191 \begin{equation} \label{E_tot_flux_metric} 
     190\begin{subequations} \label{eq:E_tot_flux} 
     191\begin{equation} \label{eq:E_tot_flux_metric} 
    192192\int\limits_D  \textbf{U}_h \cdot \text {COR} \;  dv   = 0   \\ 
    193193\end{equation} 
    194 \begin{equation} \label{E_tot_flux_adv} 
     194\begin{equation} \label{eq:E_tot_flux_adv} 
    195195   \int\limits_D \textbf{U}_h \cdot \text{ADV}   \;dv  
    196196+   \frac{1}{2} \int\limits_D {  {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t e_3  \;dv } =\;0  \\ 
    197197\end{equation} 
    198 \begin{equation} \label{E_tot_pg} 
     198\begin{equation} \label{eq:E_tot_pg} 
    199199   - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv  
    200200= - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv 
    … …  
    207207 
    208208 
    209 \eqref{E_tot_pg} is the balance between the conversion KE to PE and PE to KE.  
    210 Indeed the left hand side of \eqref{E_tot_pg} can be transformed as follows: 
     209\autoref{eq:E_tot_pg} is the balance between the conversion KE to PE and PE to KE.  
     210Indeed the left hand side of \autoref{eq:E_tot_pg} can be transformed as follows: 
    211211\begin{flalign*} 
    212212\partial_t  \left( \int\limits_D { \rho \, g \, z  \;dv} \right)  
    … …  
    221221\end{flalign*} 
    222222where the last equality is obtained by noting that the brackets is exactly the expression of $w$,  
    223 the vertical velocity referenced to the fixe $z$-coordinate system (see \eqref{Apdx_A_w_s}).  
     223the vertical velocity referenced to the fixe $z$-coordinate system (see \autoref{apdx:A_w_s}).  
    224224  
    225 The left hand side of \eqref{E_tot_pg} can be transformed as follows: 
     225The left hand side of \autoref{eq:E_tot_pg} can be transformed as follows: 
    226226\begin{flalign*} 
    227227- \int\limits_D  \left. \nabla p \right|_z & \cdot \textbf{U}_h \;dv   
    … …  
    325325% ================================================================ 
    326326\section{Discrete total energy conservation: vector invariant form} 
    327 \label{Apdx_C.1} 
     327\label{sec:C.1} 
    328328 
    329329% ------------------------------------------------------------------------------------------------------------- 
    … …  
    331331% ------------------------------------------------------------------------------------------------------------- 
    332332\subsection{Total energy conservation} 
    333 \label{Apdx_C_KE+PE} 
    334  
    335 The discrete form of the total energy conservation, \eqref{Eq_Tot_Energy}, is given by: 
     333\label{subsec:C_KE+PE} 
     334 
     335The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by: 
    336336\begin{flalign*} 
    337337\partial_t  \left(  \sum\limits_{i,j,k} \biggl\{ \frac{u^2}{2} \,b_u + \frac{v^2}{2}\, b_v +  \rho \, g \, z_t \,b_t  \biggr\} \right) &=0  \\ 
    338338\end{flalign*} 
    339339which in vector invariant forms, it leads to: 
    340 \begin{equation} \label{KE+PE_vect_discrete}   \begin{split} 
     340\begin{equation} \label{eq:KE+PE_vect_discrete}   \begin{split} 
    341341                        \sum\limits_{i,j,k} \biggl\{   u\,                        \partial_t u         \;b_u  
    342342                                                              + v\,                        \partial_t v          \;b_v  \biggr\} 
    … …  
    348348 
    349349Substituting the discrete expression of the time derivative of the velocity either in vector invariant, 
    350 leads to the discrete equivalent of the four equations \eqref{E_tot_flux}.  
     350leads to the discrete equivalent of the four equations \autoref{eq:E_tot_flux}.  
    351351 
    352352% ------------------------------------------------------------------------------------------------------------- 
    … …  
    354354% ------------------------------------------------------------------------------------------------------------- 
    355355\subsection{Vorticity term (coriolis + vorticity part of the advection)} 
    356 \label{Apdx_C_vor} 
     356\label{subsec:C_vor} 
    357357 
    358358Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$), or   
    … …  
    364364% ------------------------------------------------------------------------------------------------------------- 
    365365\subsubsection{Vorticity term with ENE scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 
    366 \label{Apdx_C_vorENE}  
     366\label{subsec:C_vorENE}  
    367367 
    368368For the ENE scheme, the two components of the vorticity term are given by : 
    … …  
    401401% ------------------------------------------------------------------------------------------------------------- 
    402402\subsubsection{Vorticity term with EEN scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 
    403 \label{Apdx_C_vorEEN}  
     403\label{subsec:C_vorEEN}  
    404404 
    405405With the EEN scheme, the vorticity terms are represented as:  
    406 \begin{equation} \label{Eq_dynvor_een} 
     406\begin{equation} \label{eq:dynvor_een} 
    407407\left\{ {    \begin{aligned} 
    408408 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}}  
    … …  
    415415$i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$, 
    416416and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by:  
    417 \begin{equation} \label{Q_triads} 
     417\begin{equation} \label{eq:Q_triads} 
    418418_i^j \mathbb{Q}^{i_p}_{j_p} 
    419419= \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right) 
    … …  
    471471% ------------------------------------------------------------------------------------------------------------- 
    472472\subsubsection{Gradient of kinetic energy / Vertical advection} 
    473 \label{Apdx_C_zad}  
     473\label{subsec:C_zad}  
    474474 
    475475The change of Kinetic Energy (KE) due to the vertical advection is exactly  
    … …  
    480480   +   \frac{1}{2} \int_D {  \frac{{\textbf{U}_h}^2}{e_3} \partial_t ( e_3) \;dv }  \\ 
    481481\end{equation*} 
    482 Indeed, using successively \eqref{DOM_di_adj} ($i.e.$ the skew symmetry  
     482Indeed, using successively \autoref{eq:DOM_di_adj} ($i.e.$ the skew symmetry  
    483483property of the $\delta$ operator) and the continuity equation, then  
    484 \eqref{DOM_di_adj} again, then the commutativity of operators  
    485 $\overline {\,\cdot \,}$ and $\delta$, and finally \eqref{DOM_mi_adj}  
     484\autoref{eq:DOM_di_adj} again, then the commutativity of operators  
     485$\overline {\,\cdot \,}$ and $\delta$, and finally \autoref{eq:DOM_mi_adj}  
    486486($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator)  
    487487applied in the horizontal and vertical directions, it becomes: 
    … …  
    536536% 
    537537\intertext{The first term provides the discrete expression for the vertical advection of momentum (ZAD),  
    538 while the second term corresponds exactly to \eqref{KE+PE_vect_discrete}, therefore:} 
     538while the second term corresponds exactly to \autoref{eq:KE+PE_vect_discrete}, therefore:} 
    539539\equiv&                   \int\limits_D \textbf{U}_h \cdot \text{ZAD} \;dv   
    540540           + \frac{1}{2} \int_D { {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t  (e_3)  \;dv }    &&&\\ 
    … …  
    568568\end{flalign*} 
    569569which is (over-)satified by defining the vertical scale factor as follows: 
    570 \begin{flalign} \label{e3u-e3v} 
     570\begin{flalign} \label{eq:e3u-e3v} 
    571571e_{3u} = \frac{1}{e_{1u}\,e_{2u}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,i+1/2}    \\ 
    572572e_{3v} = \frac{1}{e_{1v}\,e_{2v}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,j+1/2}  
    … …  
    580580% ------------------------------------------------------------------------------------------------------------- 
    581581\subsection{Pressure gradient term} 
    582 \label{Apdx_C.1.4} 
     582\label{subsec:C.1.4} 
    583583 
    584584\gmcomment{ 
    585585A pressure gradient has no contribution to the evolution of the vorticity as the  
    586586curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally  
    587 on a C-grid with 2nd order finite differences (property \eqref{Eq_DOM_curl_grad}).  
     587on a C-grid with 2nd order finite differences (property \autoref{eq:DOM_curl_grad}).  
    588588} 
    589589 
    … …  
    611611% 
    612612\allowdisplaybreaks 
    613 \intertext{Using successively  \eqref{DOM_di_adj}, $i.e.$ the skew symmetry property of  
    614 the $\delta$ operator, \eqref{Eq_wzv}, the continuity equation, \eqref{Eq_dynhpg_sco},  
     613\intertext{Using successively \autoref{eq:DOM_di_adj}, $i.e.$ the skew symmetry property of  
     614the $\delta$ operator, \autoref{eq:wzv}, the continuity equation, \autoref{dynhpg_sco},  
    615615the hydrostatic equation in the $s$-coordinate, and $\delta_{k+1/2} \left[ z_t \right] \equiv e_{3w} $, 
    616616which comes from the definition of $z_t$, it becomes: } 
    … …  
    657657% 
    658658\end{flalign*} 
    659 The first term is exactly the first term of the right-hand-side of \eqref{KE+PE_vect_discrete}. 
     659The first term is exactly the first term of the right-hand-side of \autoref{eq:KE+PE_vect_discrete}. 
    660660It remains to demonstrate that the last term, which is obviously a discrete analogue of  
    661 $\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to the last term of \eqref{KE+PE_vect_discrete}. 
     661$\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to the last term of \autoref{eq:KE+PE_vect_discrete}. 
    662662In other words, the following property must be satisfied: 
    663663\begin{flalign*} 
    … …  
    733733% ================================================================ 
    734734\section{Discrete total energy conservation: flux form} 
    735 \label{Apdx_C.1} 
     735\label{sec:C.1} 
    736736 
    737737% ------------------------------------------------------------------------------------------------------------- 
    … …  
    739739% ------------------------------------------------------------------------------------------------------------- 
    740740\subsection{Total energy conservation} 
    741 \label{Apdx_C_KE+PE} 
    742  
    743 The discrete form of the total energy conservation, \eqref{Eq_Tot_Energy}, is given by: 
     741\label{subsec:C_KE+PE} 
     742 
     743The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by: 
    744744\begin{flalign*} 
    745745\partial_t \left(  \sum\limits_{i,j,k} \biggl\{ \frac{u^2}{2} \,b_u + \frac{v^2}{2}\, b_v +  \rho \, g \, z_t \,b_t  \biggr\} \right) &=0  \\ 
    … …  
    763763% ------------------------------------------------------------------------------------------------------------- 
    764764\subsection{Coriolis and advection terms: flux form} 
    765 \label{Apdx_C.1.3} 
     765\label{subsec:C.1.3} 
    766766 
    767767% ------------------------------------------------------------------------------------------------------------- 
    … …  
    769769% ------------------------------------------------------------------------------------------------------------- 
    770770\subsubsection{Coriolis plus ``metric'' term} 
    771 \label{Apdx_C.1.3.1}  
     771\label{subsec:C.1.3.1}  
    772772 
    773773In flux from the vorticity term reduces to a Coriolis term in which the Coriolis  
    … …  
    783783Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form.  
    784784It therefore conserves the total KE. The derivation is the same as for the  
    785 vorticity term in the vector invariant form (\S\ref{Apdx_C_vor}). 
     785vorticity term in the vector invariant form (\autoref{subsec:C_vor}). 
    786786 
    787787% ------------------------------------------------------------------------------------------------------------- 
    … …  
    789789% ------------------------------------------------------------------------------------------------------------- 
    790790\subsubsection{Flux form advection} 
    791 \label{Apdx_C.1.3.2}  
     791\label{subsec:C.1.3.2}  
    792792 
    793793The flux form operator of the momentum advection is evaluated using a  
    … …  
    797797the horizontal kinetic energy, that is : 
    798798 
    799 \begin{equation} \label{Apdx_C_ADV_KE_flux} 
     799\begin{equation} \label{eq:C_ADV_KE_flux} 
    800800 -  \int_D \textbf{U}_h \cdot     \left(                 {{\begin{array} {*{20}c} 
    801801\nabla \cdot \left( \textbf{U}\,u \right) \hfill \\ 
    … …  
    856856which is the discrete form of  
    857857$ \frac{1}{2} \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv $.  
    858 \eqref{Apdx_C_ADV_KE_flux} is thus satisfied. 
     858\autoref{eq:C_ADV_KE_flux} is thus satisfied. 
    859859 
    860860 
    … …  
    877877% ================================================================ 
    878878\section{Discrete enstrophy conservation} 
    879 \label{Apdx_C.1} 
     879\label{sec:C.1} 
    880880 
    881881 
    … …  
    884884% ------------------------------------------------------------------------------------------------------------- 
    885885\subsubsection{Vorticity term with ENS scheme  (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 
    886 \label{Apdx_C_vorENS}  
     886\label{subsec:C_vorENS}  
    887887 
    888888In the ENS scheme, the vorticity term is descretized as follows: 
    889 \begin{equation} \label{Eq_dynvor_ens} 
     889\begin{equation} \label{eq:dynvor_ens} 
    890890\left\{   \begin{aligned} 
    891891+\frac{1}{e_{1u}} & \overline{q}^{\,i}  & {\overline{ \overline{\left( e_{1v}\,e_{3v}\;  v \right) } } }^{\,i, j+1/2}    \\ 
    … …  
    896896The scheme does not allow but the conservation of the total kinetic energy but the conservation  
    897897of $q^2$, the potential enstrophy for a horizontally non-divergent flow ($i.e.$ when $\chi$=$0$).  
    898 Indeed, using the symmetry or skew symmetry properties of the operators (Eqs \eqref{DOM_mi_adj}  
    899 and \eqref{DOM_di_adj}), it can be shown that: 
    900 \begin{equation} \label{Apdx_C_1.1} 
     898Indeed, using the symmetry or skew symmetry properties of the operators ( \autoref{eq:DOM_mi_adj}  
     899and \autoref{eq:DOM_di_adj}), it can be shown that: 
     900\begin{equation} \label{eq:C_1.1} 
    901901\int_D {q\,\;{\textbf{k}}\cdot \frac{1} {e_3} \nabla \times \left( {e_3 \, q \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} \equiv 0 
    902902\end{equation} 
    903903where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element. Indeed, using  
    904 \eqref{Eq_dynvor_ens}, the discrete form of the right hand side of \eqref{Apdx_C_1.1}  
     904\autoref{eq:dynvor_ens}, the discrete form of the right hand side of \autoref{eq:C_1.1}  
    905905can be transformed as follow: 
    906906\begin{flalign*}  
    … …  
    944944% ------------------------------------------------------------------------------------------------------------- 
    945945\subsubsection{Vorticity Term with EEN scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 
    946 \label{Apdx_C_vorEEN}  
     946\label{subsec:C_vorEEN}  
    947947 
    948948With the EEN scheme, the vorticity terms are represented as:  
    949 \begin{equation} \label{Eq_dynvor_een} 
     949\begin{equation} \label{eq:dynvor_een} 
    950950\left\{ {    \begin{aligned} 
    951951 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}}  
    … …  
    958958$i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$, 
    959959and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by:  
    960 \begin{equation} \label{Q_triads} 
     960\begin{equation} \label{eq:Q_triads} 
    961961_i^j \mathbb{Q}^{i_p}_{j_p} 
    962962= \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right) 
    … …  
    968968Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2} $,  
    969969similar manipulation can be done for the 3 others. The discrete form of the right hand  
    970 side of \eqref{Apdx_C_1.1} applied to this triad only can be transformed as follow: 
     970side of \autoref{eq:C_1.1} applied to this triad only can be transformed as follow: 
    971971 
    972972\begin{flalign*}  
    … …  
    10171017% ================================================================ 
    10181018\section{Conservation properties on tracers} 
    1019 \label{Apdx_C.2} 
     1019\label{sec:C.2} 
    10201020 
    10211021 
    … …  
    10331033% ------------------------------------------------------------------------------------------------------------- 
    10341034\subsection{Advection term} 
    1035 \label{Apdx_C.2.1} 
     1035\label{subsec:C.2.1} 
    10361036 
    10371037conservation of a tracer, $T$: 
    … …  
    11001100% ================================================================ 
    11011101\section{Conservation properties on lateral momentum physics} 
    1102 \label{Apdx_dynldf_properties} 
     1102\label{sec:dynldf_properties} 
    11031103 
    11041104 
    … …  
    11141114 
    11151115These properties of the horizontal diffusion operator are a direct consequence  
    1116 of properties \eqref{Eq_DOM_curl_grad} and \eqref{Eq_DOM_div_curl}.  
     1116of properties \autoref{eq:DOM_curl_grad} and \autoref{eq:DOM_div_curl}.  
    11171117When the vertical curl of the horizontal diffusion of momentum (discrete sense)  
    11181118is taken, the term associated with the horizontal gradient of the divergence is  
    … …  
    11231123% ------------------------------------------------------------------------------------------------------------- 
    11241124\subsection{Conservation of potential vorticity} 
    1125 \label{Apdx_C.3.1} 
     1125\label{subsec:C.3.1} 
    11261126 
    11271127The lateral momentum diffusion term conserves the potential vorticity : 
    … …  
    11431143   \right\}     \\  
    11441144% 
    1145 \intertext{Using \eqref{DOM_di_adj}, it follows:} 
     1145\intertext{Using \autoref{eq:DOM_di_adj}, it follows:} 
    11461146% 
    11471147\equiv& \sum\limits_{i,j,k}  
    … …  
    11571157% ------------------------------------------------------------------------------------------------------------- 
    11581158\subsection{Dissipation of horizontal kinetic energy} 
    1159 \label{Apdx_C.3.2} 
     1159\label{subsec:C.3.2} 
    11601160 
    11611161The lateral momentum diffusion term dissipates the horizontal kinetic energy: 
    … …  
    12091209% ------------------------------------------------------------------------------------------------------------- 
    12101210\subsection{Dissipation of enstrophy} 
    1211 \label{Apdx_C.3.3} 
     1211\label{subsec:C.3.3} 
    12121212 
    12131213The lateral momentum diffusion term dissipates the enstrophy when the eddy  
    … …  
    12231223             + \delta_{j+1/2} \left[  \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta  \right]   \right]      \right\}   &&&\\  
    12241224% 
    1225 \intertext{Using \eqref{DOM_di_adj}, it follows:} 
     1225\intertext{Using \autoref{eq:DOM_di_adj}, it follows:} 
    12261226% 
    12271227&\quad \equiv  - A^{\,lm} \sum\limits_{i,j,k}  
    … …  
    12341234% ------------------------------------------------------------------------------------------------------------- 
    12351235\subsection{Conservation of horizontal divergence} 
    1236 \label{Apdx_C.3.4} 
     1236\label{subsec:C.3.4} 
    12371237 
    12381238When the horizontal divergence of the horizontal diffusion of momentum  
    12391239(discrete sense) is taken, the term associated with the vertical curl of the  
    1240 vorticity is zero locally, due to \eqref{Eq_DOM_div_curl}.  
     1240vorticity is zero locally, due to \autoref{eq:DOM_div_curl}.  
    12411241The resulting term conserves the $\chi$ and dissipates $\chi^2$  
    12421242when the eddy coefficients are horizontally uniform. 
    … …  
    12511251           + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}}  \delta_{j+1/2} \left[ \chi \right]  \right]    \right\}    \\  
    12521252% 
    1253 \intertext{Using \eqref{DOM_di_adj}, it follows:} 
     1253\intertext{Using \autoref{eq:DOM_di_adj}, it follows:} 
    12541254% 
    12551255&\equiv \sum\limits_{i,j,k}  
    … …  
    12631263% ------------------------------------------------------------------------------------------------------------- 
    12641264\subsection{Dissipation of horizontal divergence variance} 
    1265 \label{Apdx_C.3.5} 
     1265\label{subsec:C.3.5} 
    12661266 
    12671267\begin{flalign*} 
    … …  
    12771277   \right\} \;   e_{1t}\,e_{2t}\,e_{3t}    \\  
    12781278% 
    1279 \intertext{Using \eqref{DOM_di_adj}, it turns out to be:} 
     1279\intertext{Using \autoref{eq:DOM_di_adj}, it turns out to be:} 
    12801280% 
    12811281&\equiv - A^{\,lm} \sum\limits_{i,j,k} 
    … …  
    12891289% ================================================================ 
    12901290\section{Conservation properties on vertical momentum physics} 
    1291 \label{Apdx_C_4} 
     1291\label{sec:C_4} 
    12921292 
    12931293As for the lateral momentum physics, the continuous form of the vertical diffusion  
    … …  
    14611461% ================================================================ 
    14621462\section{Conservation properties on tracer physics} 
    1463 \label{Apdx_C.5} 
     1463\label{sec:C.5} 
    14641464 
    14651465The numerical schemes used for tracer subgridscale physics are written such  
    … …  
    14731473% ------------------------------------------------------------------------------------------------------------- 
    14741474\subsection{Conservation of tracers} 
    1475 \label{Apdx_C.5.1} 
     1475\label{subsec:C.5.1} 
    14761476 
    14771477constraint of conservation of tracers: 
    … …  
    15071507% ------------------------------------------------------------------------------------------------------------- 
    15081508\subsection{Dissipation of tracer variance} 
    1509 \label{Apdx_C.5.2} 
     1509\label{subsec:C.5.2} 
    15101510 
    15111511constraint on the dissipation of tracer variance: 

  • branches/2017/dev_merge_2017/DOC/tex_sub/annex_D.tex

    r9393 r9407  
    55% ================================================================ 
    66\chapter{Coding Rules} 
    7 \label{Apdx_D} 
     7\label{apdx:D} 
    88\minitoc 
    99 
    … …  
    4747% ================================================================ 
    4848\section{Program structure} 
    49 \label{Apdx_D_structure} 
     49\label{sec:D_structure} 
    5050 
    5151Each program begins with a set of headline comments containing : 
    … …  
    7676% ================================================================ 
    7777\section{Coding conventions} 
    78 \label{Apdx_D_coding} 
     78\label{sec:D_coding} 
    7979 
    8080- Use of the universal language \textsc{Fortran} 90, and try to avoid obsolescent 
    … …  
    107107% ================================================================ 
    108108\section{Naming conventions} 
    109 \label{Apdx_D_naming} 
     109\label{sec:D_naming} 
    110110 
    111111The purpose of the naming conventions is to use prefix letters to classify  
    … …  
    116116 
    117117%--------------------------------------------------TABLE-------------------------------------------------- 
    118 \begin{table}[htbp]  \label{Tab_VarName} 
     118\begin{table}[htbp]  \label{tab:VarName} 
    119119\begin{center} 
    120120\begin{tabular}{|p{45pt}|p{35pt}|p{45pt}|p{40pt}|p{40pt}|p{40pt}|p{40pt}|p{40pt}|} 
    … …  
    187187\hline 
    188188\end{tabular} 
    189 \label{tab1} 
     189\label{tab:tab1} 
    190190\end{center} 
    191191\end{table} 
    … …  
    201201% ================================================================ 
    202202%\section{Program structure} 
    203 %label{Apdx_D_structure} 
     203%abel{sec:Apdx_D_structure} 
    204204 
    205205%To be done.... 

  • branches/2017/dev_merge_2017/DOC/tex_sub/annex_E.tex

    r9393 r9407  
    55% ================================================================ 
    66\chapter{Note on some algorithms} 
    7 \label{Apdx_E} 
     7\label{apdx:E} 
    88\minitoc 
    99 
    … …  
    2020% ------------------------------------------------------------------------------------------------------------- 
    2121\section{Upstream Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})} 
    22 \label{TRA_adv_ubs} 
     22\label{sec:TRA_adv_ubs} 
    2323 
    2424The UBS advection scheme is an upstream biased third order scheme based on  
    … …  
    2626QUICK scheme (Quadratic Upstream Interpolation for Convective  
    2727Kinematics). For example, in the $i$-direction : 
    28 \begin{equation} \label{Eq_tra_adv_ubs2} 
     28\begin{equation} \label{eq:tra_adv_ubs2} 
    2929\tau _u^{ubs} = \left\{  \begin{aligned} 
    3030  & \tau _u^{cen4} + \frac{1}{12} \,\tau"_i     & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ 
    … …  
    3333\end{equation} 
    3434or equivalently, the advective flux is 
    35 \begin{equation} \label{Eq_tra_adv_ubs2} 
     35\begin{equation} \label{eq:tra_adv_ubs2} 
    3636U_{i+1/2} \ \tau _u^{ubs}  
    3737=U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2} 
    … …  
    6161scheme when \np{ln\_traadv\_ubs}\forcode{ = .true.}. 
    6262 
    63 For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds  
     63For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds  
    6464to a second order centred scheme is evaluated using the \textit{now} velocity  
    6565(centred in time) while the second term which is the diffusive part of the scheme,  
    … …  
    6767by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. UBS and QUICK  
    6868schemes only differ by one coefficient. Substituting 1/6 with 1/8 in  
    69 (\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.  
     69(\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.  
    7070This option is not available through a namelist parameter, since the 1/6  
    7171coefficient is hard coded. Nevertheless it is quite easy to make the  
    … …  
    8787eight-order accurate conventional scheme. 
    8888 
    89 NB 3 : It is straight forward to rewrite \eqref{Eq_tra_adv_ubs} as follows: 
    90 \begin{equation} \label{Eq_tra_adv_ubs2} 
     89NB 3 : It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows: 
     90\begin{equation} \label{eq:tra_adv_ubs2} 
    9191\tau _u^{ubs} = \left\{  \begin{aligned} 
    9292   & \tau _u^{cen4} + \frac{1}{12} \tau"_i      & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ 
    … …  
    9595\end{equation} 
    9696or equivalently  
    97 \begin{equation} \label{Eq_tra_adv_ubs2} 
     97\begin{equation} \label{eq:tra_adv_ubs2} 
    9898\begin{split} 
    9999e_{2u} e_{3u}\,u_{i+1/2} \ \tau _u^{ubs}  
    … …  
    102102\end{split} 
    103103\end{equation} 
    104 \eqref{Eq_tra_adv_ubs2} has several advantages. First it clearly evidence that  
     104\autoref{eq:tra_adv_ubs2} has several advantages. First it clearly evidence that  
    105105the UBS scheme is based on the fourth order scheme to which is added an  
    106106upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order  
    107107part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order  
    108 part as stated above using \eqref{Eq_tra_adv_ubs}. Third, the diffusive term is  
     108part as stated above using \autoref{eq:tra_adv_ubs}. Third, the diffusive term is  
    109109in fact a biharmonic operator with a eddy coefficient with is simply proportional  
    110110to the velocity. 
    111111 
    112112laplacian diffusion: 
    113 \begin{equation} \label{Eq_tra_ldf_lap} 
     113\begin{equation} \label{eq:tra_ldf_lap} 
    114114\begin{split} 
    115115D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta _i  
    … …  
    124124 
    125125bilaplacian: 
    126 \begin{equation} \label{Eq_tra_ldf_lap} 
     126\begin{equation} \label{eq:tra_ldf_lap} 
    127127\begin{split} 
    128128D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\ 
    … …  
    136136$i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$ 
    137137it comes : 
    138 \begin{equation} \label{Eq_tra_ldf_lap} 
     138\begin{equation} \label{eq:tra_ldf_lap} 
    139139\begin{split} 
    140140D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\ 
    … …  
    146146\end{equation} 
    147147if the velocity is uniform ($i.e.$ $|u|=cst$) then the diffusive flux is 
    148 \begin{equation} \label{Eq_tra_ldf_lap} 
     148\begin{equation} \label{eq:tra_ldf_lap} 
    149149\begin{split} 
    150150F_u^{lT} = - \frac{1}{12} 
    … …  
    157157beurk....  reverte the logic: starting from the diffusive part of the advective flux it comes: 
    158158 
    159 \begin{equation} \label{Eq_tra_adv_ubs2} 
     159\begin{equation} \label{eq:tra_adv_ubs2} 
    160160\begin{split} 
    161161F_u^{lT} 
    … …  
    166166 
    167167sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$): 
    168 \begin{equation} \label{Eq_tra_adv_ubs2} 
     168\begin{equation} \label{eq:tra_adv_ubs2} 
    169169\begin{split} 
    170170F_u^{lT} 
    … …  
    175175 
    176176sol 2 coefficient at u-point: split $|u|$ into $\sqrt{|u|}$ and $e_{1T}$ into $\sqrt{e_{1u}}$ 
    177 \begin{equation} \label{Eq_tra_adv_ubs2} 
     177\begin{equation} \label{eq:tra_adv_ubs2} 
    178178\begin{split} 
    179179F_u^{lT} 
    … …  
    189189% ------------------------------------------------------------------------------------------------------------- 
    190190\section{Leapfrog energetic} 
    191 \label{LF} 
     191\label{sec:LF} 
    192192 
    193193We adopt the following semi-discrete notation for time derivative. Given the values of a variable $q$ at successive time step, the time derivation and averaging operators at the mid time step are: 
    194 \begin{subequations} \label{dt_mt} 
     194\begin{subequations} \label{eq:dt_mt} 
    195195\begin{align} 
    196196 \delta _{t+\rdt/2} [q]     &=  \  \ \,   q^{t+\rdt}  - q^{t}     \\ 
    … …  
    202202, respectively.  
    203203 
    204 The Leap-frog time stepping given by \eqref{Eq_DOM_nxt} can be defined as: 
    205 \begin{equation} \label{LF} 
     204The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as: 
     205\begin{equation} \label{eq:LF} 
    206206   \frac{\partial q}{\partial t}  
    207207         \equiv \frac{1}{\rdt} \overline{ \delta _{t+\rdt/2}[q]}^{\,t}  
    208208      =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt} 
    209209\end{equation}  
    210 Note that \eqref{LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$  
     210Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$  
    211211as it can be found sometime in literature.  
    212212The leap-Frog time stepping is a second order centered scheme. As such it respects  
    213213the quadratic invariant in integral forms, $i.e.$ the following continuous property, 
    214 \begin{equation} \label{Energy} 
     214\begin{equation} \label{eq:Energy} 
    215215\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}  
    216216   =\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt}  
    … …  
    252252scheme, but is formulated within the \NEMO framework ($i.e.$ using scale  
    253253factors rather than grid-size and having a position of $T$-points that is not  
    254 necessary in the middle of vertical velocity points, see Fig.~\ref{Fig_zgr_e3}). 
    255  
    256 In the formulation \eqref{Eq_tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO,  
     254necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}). 
     255 
     256In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO,  
    257257the off-diagonal terms of the small angle diffusion tensor contain several double  
    258258spatial averages of a gradient, for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$.  
    … …  
    263263In other word, the operator applied to a tracer does not warranties the decrease of  
    264264its global average variance. To circumvent this, we have introduced a smoothing of  
    265 the slopes of the iso-neutral surfaces (see \S\ref{LDF}). Nevertheless, this technique  
     265the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). Nevertheless, this technique  
    266266works fine for $T$ and $S$ as they are active tracers ($i.e.$ they enter the computation  
    267267of density), but it does not work for a passive tracer.   \citep{Griffies_al_JPO98} introduce  
    … …  
    270270with a derivative in the same direction by considering triads. For example in the  
    271271(\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows: 
    272 \begin{equation} \label{Gf_triads} 
     272\begin{equation} \label{eq:Gf_triads} 
    273273_i^k \mathbb{T}_{i_p}^{k_p} (T) 
    274274= \frac{1}{4} \ {b_u}_{\,i+i_p}^{\,k}  \  A_i^k    \left(   
    … …  
    282282$A_i^k$ is the lateral eddy diffusivity coefficient defined at $T$-point, 
    283283and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad : 
    284 \begin{equation} \label{Gf_slopes} 
     284\begin{equation} \label{eq:Gf_slopes} 
    285285_i^k \mathbb{R}_{i_p}^{k_p}  
    286286=\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} \ \frac  
    … …  
    288288{\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] } 
    289289\end{equation} 
    290 Note that in \eqref{Gf_slopes} we use the ratio $\alpha / \beta$ instead of  
     290Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of  
    291291multiplying the temperature derivative by $\alpha$ and the salinity derivative  
    292292by $\beta$. This is more efficient as the ratio $\alpha / \beta$ can to be  
    293293evaluated directly. 
    294294 
    295 Note that in \eqref{Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of  
     295Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of  
    296296${b_{uw}}_{\,i+i_p}^{\,k+k_p}$. This choice has been motivated by the decrease  
    297297of tracer variance and the presence of partial cell at the ocean bottom  
    298 (see Appendix~\ref{Apdx_Gf_operator}). 
     298(see \autoref{apdx:Gf_operator}). 
    299299 
    300300%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    301 \begin{figure}[!ht] \label{Fig_ISO_triad} 
    302 \begin{center} 
     301\begin{figure}[!ht] \begin{center} 
    303302\includegraphics[width=0.70\textwidth]{Fig_ISO_triad} 
    304 \caption{  \protect\label{Fig_ISO_triad}    
     303\caption{  \protect\label{fig:ISO_triad}    
    305304Triads used in the Griffies's like iso-neutral diffision scheme for  
    306305$u$-component (upper panel) and $w$-component (lower panel).} 
    … …  
    311310The four iso-neutral fluxes associated with the triads are defined at $T$-point.  
    312311They take the following expression : 
    313 \begin{flalign} \label{Gf_fluxes} 
     312\begin{flalign} \label{eq:Gf_fluxes} 
    314313\begin{split} 
    315314{_i^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)  
    … …  
    322321 
    323322The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the  
    324 sum of the fluxes that cross the $u$- and $w$-face (Fig.~\ref{Fig_ISO_triad}): 
    325 \begin{flalign} \label{Eq_iso_flux}  
     323sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}): 
     324\begin{flalign} \label{eq:iso_flux}  
    326325\textbf{F}_{iso}(T)  
    327326&\equiv  \sum_{\substack{i_p,\,k_p}}  
    … …  
    353352resulting in a iso-neutral diffusion tendency on temperature given by the divergence  
    354353of the sum of all the four triad fluxes : 
    355 \begin{equation} \label{Gf_operator} 
     354\begin{equation} \label{eq:Gf_operator} 
    356355D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{   
    357356       \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]  
    … …  
    365364\item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator  
    366365recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction : 
    367 \begin{equation} \label{Gf_property1a} 
     366\begin{equation} \label{eq:Gf_property1a} 
    368367D_l^T = \frac{1}{b_T}  \ \delta_{i}  
    369368   \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right]  
    … …  
    388387\item[$\bullet$ pure iso-neutral operator]  The iso-neutral flux of locally referenced  
    389388potential density is zero, $i.e.$ 
    390 \begin{align} \label{Gf_property2} 
     389\begin{align} \label{eq:Gf_property2} 
    391390\begin{matrix} 
    392391&{_i^k {\mathbb{F}_u}_{i_p}^{k_p} (\rho)}  
    … …  
    398397\end{matrix} 
    399398\end{align} 
    400 This result is trivially obtained using the \eqref{Gf_triads} applied to $T$ and $S$  
    401 and the definition of the triads' slopes \eqref{Gf_slopes}. 
     399This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$  
     400and the definition of the triads' slopes \autoref{eq:Gf_slopes}. 
    402401 
    403402\item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the  
    404403total tracer content, $i.e.$ 
    405 \begin{equation} \label{Gf_property1} 
     404\begin{equation} \label{eq:Gf_property1} 
    406405\sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0 
    407406\end{equation} 
    … …  
    411410\item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does  
    412411not increase the total tracer variance, $i.e.$ 
    413 \begin{equation} \label{Gf_property1} 
     412\begin{equation} \label{eq:Gf_property1} 
    414413\sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0 
    415414\end{equation} 
    416 The property is demonstrated in the Appendix~\ref{Apdx_Gf_operator}. It is a  
     415The property is demonstrated in the \autoref{apdx:Gf_operator}. It is a  
    417416key property for a diffusion term. It means that the operator is also a dissipation  
    418417term, $i.e.$ it is a sink term for the square of the quantity on which it is applied.  
    … …  
    422421\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint,  
    423422$i.e.$ 
    424 \begin{equation} \label{Gf_property1} 
     423\begin{equation} \label{eq:Gf_property1} 
    425424\sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}  
    426425\end{equation} 
    … …  
    428427operator. We just have to apply the same routine. This properties can be demonstrated  
    429428quite easily in a similar way the "non increase of tracer variance" property has been  
    430 proved (see Appendix~\ref{Apdx_Gf_operator}). 
     429proved (see \autoref{apdx:Gf_operator}). 
    431430\end{description} 
    432431 
    … …  
    442441eddy induced velocity, the formulation of which depends on the slopes of iso- 
    443442neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used  
    444 here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo}  
    445 is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo} 
    446 + \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates.  
     443here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo}  
     444is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} 
     445+ \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.  
    447446 
    448447The eddy induced velocity is given by:  
    449 \begin{equation} \label{Eq_eiv_v} 
     448\begin{equation} \label{eq:eiv_v} 
    450449\begin{split} 
    451450 u^* & = - \frac{1}{e_2\,e_{3}}          \;\partial_k \left( e_2 \, A_e \; r_i  \right)    
    … …  
    467466A traditional way to implement this additional advection is to add it to the eulerian  
    468467velocity prior to compute the tracer advection. This allows us to take advantage of  
    469 all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not just  
     468all the advection schemes offered for the tracers (see \autoref{sec:TRA_adv}) and not just  
    470469a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers  
    471470where \emph{positivity} of the advection scheme is of paramount importance.  
    472 % give here the expression using the triads. It is different from the one given in \eqref{Eq_ldfeiv} 
     471% give here the expression using the triads. It is different from the one given in \autoref{eq:ldfeiv} 
    473472% see just below a copy of this equation: 
    474 %\begin{equation} \label{Eq_ldfeiv} 
     473%\begin{equation} \label{eq:ldfeiv} 
    475474%\begin{split} 
    476475% u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\ 
    … …  
    479478%\end{split} 
    480479%\end{equation} 
    481 \begin{equation} \label{Eq_eiv_vd}   
     480\begin{equation} \label{eq:eiv_vd}   
    482481\textbf{F}_{eiv}^T   \equiv   \left( \begin{aligned}                                 
    483482 \sum_{\substack{i_p,\,k_p}} & 
    … …  
    491490\end{equation} 
    492491 
    493 \ref{Griffies_JPO98} introduces another way to implement the eddy induced advection,  
     492\citep{Griffies_JPO98} introduces another way to implement the eddy induced advection,  
    494493the so-called skew form. It is based on a transformation of the advective fluxes  
    495494using the non-divergent nature of the eddy induced velocity.  
    … …  
    522521and since the eddy induces velocity field is no-divergent, we end up with the skew  
    523522form of the eddy induced advective fluxes: 
    524 \begin{equation} \label{Eq_eiv_skew_continuous} 
     523\begin{equation} \label{eq:eiv_skew_continuous} 
    525524\textbf{F}_{eiv}^T = \begin{pmatrix}  
    526525           {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\ 
    … …  
    529528\end{equation} 
    530529The tendency associated with eddy induced velocity is then simply the divergence  
    531 of the \eqref{Eq_eiv_skew_continuous} fluxes. It naturally conserves the tracer  
     530of the \autoref{eq:eiv_skew_continuous} fluxes. It naturally conserves the tracer  
    532531content, as it is expressed in flux form and, as the advective form, it preserve the  
    533 tracer variance. Another interesting property of \eqref{Eq_eiv_skew_continuous}  
     532tracer variance. Another interesting property of \autoref{eq:eiv_skew_continuous}  
    534533form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral  
    535534diffusion and eddy induced velocity terms: 
    536 \begin{flalign} \label{Eq_eiv_skew+eiv_continuous} 
     535\begin{flalign} \label{eq:eiv_skew+eiv_continuous} 
    537536\textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=  
    538537\begin{pmatrix}  
    … …  
    554553has been used to reduce the computational time \citep{Griffies_JPO98}, but it is  
    555554not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to  
    556 choose a discret form of  \eqref{Eq_eiv_skew_continuous} which is consistent with the  
    557 iso-neutral operator \eqref{Gf_operator}. Using the slopes \eqref{Gf_slopes}  
     555choose a discret form of  \autoref{eq:eiv_skew_continuous} which is consistent with the  
     556iso-neutral operator \autoref{eq:Gf_operator}. Using the slopes \autoref{eq:Gf_slopes}  
    558557and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient), 
    559558the resulting discret form is given by: 
    560 \begin{equation} \label{Eq_eiv_skew}   
     559\begin{equation} \label{eq:eiv_skew}   
    561560\textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left( \begin{aligned}                                 
    562561 \sum_{\substack{i_p,\,k_p}} & 
    … …  
    569568\end{aligned}   \right) 
    570569\end{equation} 
    571 Note that \eqref{Eq_eiv_skew} is valid in $z$-coordinate with or without partial cells.  
     570Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells.  
    572571In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces  
    573572must be added to $\mathbb{R}$ for the discret form to be exact.  
    … …  
    575574Such a choice of discretisation is consistent with the iso-neutral operator as it uses the  
    576575same definition for the slopes. It also ensures the conservation of the tracer variance  
    577 (see Appendix \ref{Apdx_eiv_skew}), $i.e.$ it does not include a diffusive component  
     576(see Appendix \autoref{apdx:eiv_skew}), $i.e.$ it does not include a diffusive component  
    578577but is a "pure" advection term. 
    579578 
    … …  
    586585% ================================================================ 
    587586\subsection{Discrete invariants of the iso-neutral diffrusion} 
    588 \label{Apdx_Gf_operator} 
     587\label{subsec:Gf_operator} 
    589588 
    590589Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.  
    … …  
    596595\int_D  D_l^T \; T \;dv   \leq 0 
    597596\end{align*} 
    598 The discrete form of its left hand side is obtained using \eqref{Eq_iso_flux} 
     597The discrete form of its left hand side is obtained using \autoref{eq:iso_flux} 
    599598 
    600599\begin{align*} 
    … …  
    673672% 
    674673\allowdisplaybreaks 
    675 \intertext{Then outing in factor the triad in each of the four terms of the summation and substituting the triads by their expression given in \eqref{Gf_triads}. It becomes: } 
     674\intertext{Then outing in factor the triad in each of the four terms of the summation and substituting the triads by their expression given in \autoref{eq:Gf_triads}. It becomes: } 
    676675% 
    677676&\equiv -\sum_{i,k} 
    … …  
    739738% ================================================================ 
    740739\subsection{Discrete invariants of the skew flux formulation} 
    741 \label{Apdx_eiv_skew} 
     740\label{subsec:eiv_skew} 
    742741 
    743742 
    … …  
    750749\int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0 
    751750\end{align*} 
    752 The discrete form of its left hand side is obtained using \eqref{Eq_eiv_skew} 
     751The discrete form of its left hand side is obtained using \autoref{eq:eiv_skew} 
    753752\begin{align*} 
    754753 \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\; 

  • branches/2017/dev_merge_2017/DOC/tex_sub/annex_iso.tex

    r9393 r9407  
    44% Iso-neutral diffusion : 
    55% ================================================================ 
    6 \chapter{Iso-neutral diffusion and eddy advection using triads} 
    7 \label{sec:triad} 
     6\chapter[Iso-Neutral Diffusion and Eddy Advection using Triads] 
     7         {\texorpdfstring{Iso-Neutral Diffusion and\\ Eddy Advection using Triads}{Iso-Neutral Diffusion and Eddy Advection using Triads}} 
     8\label{apdx:triad} 
    89\minitoc 
    910\pagebreak 
    … …  
    1819of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes.  
    1920If the namelist logical \np{ln\_traldf\_iso} is set true,  
    20 the filtered version of Cox's original scheme (the Standard scheme) is employed (\S\ref{LDF_slp}).  
     21the filtered version of Cox's original scheme (the Standard scheme) is employed (\autoref{sec:LDF_slp}).  
    2122In the present implementation of the Griffies scheme,  
    2223the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false. 
    2324 
    2425Values of iso-neutral diffusivity and GM coefficient are set as 
    25 described in \S\ref{LDF_coef}. Note that when GM fluxes are used,  
     26described in \autoref{sec:LDF_coef}. Note that when GM fluxes are used,  
    2627the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS,  
    2728even though the eddy advection is accomplished by means of the skew fluxes. 
    … …  
    3031The options specific to the Griffies scheme include: 
    3132\begin{description}[font=\normalfont] 
    32 \item[\np{ln\_triad\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then 
     33\item[\np{ln\_triad\_iso}] See \autoref{sec:taper}. If this is set false (the default), then 
    3334  `iso-neutral' mixing is accomplished within the surface mixed-layer 
    3435  along slopes linearly decreasing with depth from the value immediately below 
    35   the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}).  
    36   This is the same treatment as used in the default implementation \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}.   
     36  the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}).  
     37  This is the same treatment as used in the default implementation \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}.   
    3738  Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced  
    3839  to ensure no vertical buoyancy flux, giving an almost pure 
    3940  horizontal diffusive tracer flux within the mixed layer. This is similar to 
    40   the tapering suggested by \citet{Gerdes1991}. See \S\ref{sec:triad:Gerdes-taper} 
    41 \item[\np{ln\_botmix\_triad}] See \S\ref{sec:triad:iso_bdry}.  
     41  the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper} 
     42\item[\np{ln\_botmix\_triad}] See \autoref{sec:iso_bdry}.  
    4243  If this is set false (the default) then the lateral diffusive fluxes 
    4344  associated with triads partly masked by topography are neglected.  
    … …  
    5354 
    5455\section{Triad formulation of iso-neutral diffusion} 
    55 \label{sec:triad:iso} 
     56\label{sec:iso} 
    5657We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98},  
    5758but formulated within the \NEMO framework, using scale factors rather than grid-sizes. 
    … …  
    6061The iso-neutral second order tracer diffusive operator for small 
    6162angles between iso-neutral surfaces and geopotentials is given by 
    62 \eqref{Eq_PE_iso_tensor}: 
    63 \begin{subequations} \label{eq:triad:PE_iso_tensor} 
     63\autoref{eq:PE_iso_tensor}: 
     64\begin{subequations} \label{eq:PE_iso_tensor} 
    6465  \begin{equation} 
    6566    D^{lT}=-\Div\vect{f}^{lT}\equiv 
    … …  
    7273  \end{equation} 
    7374  \begin{equation} 
    74     \label{eq:triad:PE_iso_tensor:c} 
     75    \label{eq:PE_iso_tensor:c} 
    7576    \mbox{with}\quad \;\;\Re = 
    7677    \begin{pmatrix} 
    … …  
    9293%  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 
    9394% \end{array} }} \right) 
    94  Here \eqref{Eq_PE_iso_slopes}  
     95 Here \autoref{eq:PE_iso_slopes}  
    9596\begin{align*} 
    9697  r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} 
    … …  
    108109space; we write 
    109110\begin{equation} 
    110   \label{eq:triad:Fijk} 
     111  \label{eq:Fijk} 
    111112  \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right). 
    112113\end{equation} 
    … …  
    117118 
    118119The off-diagonal terms of the small angle diffusion tensor 
    119 \eqref{Eq_PE_iso_tensor}, \eqref{eq:triad:PE_iso_tensor:c} produce skew-fluxes along the 
     120\autoref{eq:PE_iso_tensor}, \autoref{eq:PE_iso_tensor:c} produce skew-fluxes along the 
    120121$i$- and $j$-directions resulting from the vertical tracer gradient: 
    121122\begin{align} 
    122   \label{eq:triad:i13c} 
     123  \label{eq:i13c} 
    123124  f_{13}=&+\Alt r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+\Alt r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\ 
    124125\intertext{and in the k-direction resulting from the lateral tracer gradients} 
    125   \label{eq:triad:i31c} 
     126  \label{eq:i31c} 
    126127 f_{31}+f_{32}=& \Alt r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\Alt r_2\frac{1}{e_1}\frac{\partial T}{\partial i} 
    127128\end{align} 
    … …  
    130131component of the small angle diffusion tensor is 
    131132\begin{equation} 
    132   \label{eq:triad:i33c} 
     133  \label{eq:i33c} 
    133134  f_{33}=-\Alt(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}. 
    134135\end{equation} 
    … …  
    141142 
    142143There is no natural discretization for the $i$-component of the 
    143 skew-flux, \eqref{eq:triad:i13c}, as 
     144skew-flux, \autoref{eq:i13c}, as 
    144145although it must be evaluated at $u$-points, it involves vertical 
    145146gradients (both for the tracer and the slope $r_1$), defined at 
    146 $w$-points. Similarly, the vertical skew flux, \eqref{eq:triad:i31c}, is evaluated at 
     147$w$-points. Similarly, the vertical skew flux, \autoref{eq:i31c}, is evaluated at 
    147148$w$-points but involves horizontal gradients defined at $u$-points. 
    148149 
    149150\subsection{Standard discretization} 
    150151The straightforward approach to discretize the lateral skew flux 
    151 \eqref{eq:triad:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 
    152 into OPA, \eqref{Eq_tra_ldf_iso}, is to calculate a mean vertical 
     152\autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 
     153into OPA, \autoref{eq:tra_ldf_iso}, is to calculate a mean vertical 
    153154gradient at the $u$-point from the average of the four surrounding 
    154155vertical tracer gradients, and multiply this by a mean slope at the 
    … …  
    159160$e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with 
    160161the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer 
    161 gradient, is then \eqref{Eq_tra_ldf_iso} 
     162gradient, is then \autoref{eq:tra_ldf_iso} 
    162163\begin{equation*} 
    163164  \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k 
    … …  
    181182operator to a tracer does not guarantee the decrease of its 
    182183global-average variance. To correct this, we introduced a smoothing of 
    183 the slopes of the iso-neutral surfaces (see \S\ref{LDF}). This 
     184the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). This 
    184185technique works for $T$ and $S$ in so far as they are active tracers 
    185186($i.e.$ they enter the computation of density), but it does not work 
    … …  
    194195\begin{figure}[tb] \begin{center} 
    195196    \includegraphics[width=1.05\textwidth]{Fig_GRIFF_triad_fluxes} 
    196     \caption{ \protect\label{fig:triad:ISO_triad} 
     197    \caption{ \protect\label{fig:ISO_triad} 
    197198      (a) Arrangement of triads $S_i$ and tracer gradients to 
    198199           give lateral tracer flux from box $i,k$ to $i+1,k$ 
    … …  
    205206slope calculated from the lateral density gradient across the $u$-point 
    206207divided by the vertical density gradient at the same $w$-point as the 
    207 tracer gradient. See Fig.~\ref{fig:triad:ISO_triad}a, where the thick lines 
     208tracer gradient. See \autoref{fig:ISO_triad}a, where the thick lines 
    208209denote the tracer gradients, and the thin lines the corresponding 
    209210triads, with slopes $s_1, \dotsc s_4$. The total area-integrated 
    210211skew-flux from tracer cell $i,k$ to $i+1,k$ 
    211212\begin{multline} 
    212   \label{eq:triad:i13} 
     213  \label{eq:i13} 
    213214  \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_1 s_1 
    214215  \delta _{k+\frac{1}{2}} \left[ T^{i+1} 
    … …  
    225226stencil, and disallows the two-point computational modes. 
    226227 
    227  The vertical skew flux \eqref{eq:triad:i31c} from tracer cell $i,k$ to $i,k+1$ at the 
    228 $w$-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{fig:triad:ISO_triad}b) 
     228 The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at the 
     229$w$-point $i,k+\hhalf$ is constructed similarly (\autoref{fig:ISO_triad}b) 
    229230by multiplying lateral tracer gradients from each of the four 
    230231surrounding $u$-points by the appropriate triad slope: 
    231232\begin{multline} 
    232   \label{eq:triad:i31} 
     233  \label{eq:i31} 
    233234  \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  \Alts_i^{k+1} a_{1}' 
    234235  s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1} 
    … …  
    241242(appearing in both the vertical and lateral gradient), and the $u$- and 
    242243$w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the 
    243 triad as follows (see also Fig.~\ref{fig:triad:ISO_triad}): 
    244 \begin{equation} 
    245   \label{eq:triad:R} 
     244triad as follows (see also \autoref{fig:ISO_triad}): 
     245\begin{equation} 
     246  \label{eq:R} 
    246247  _i^k \mathbb{R}_{i_p}^{k_p} 
    247248  =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} 
    … …  
    258259\begin{figure}[tb] \begin{center} 
    259260    \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells} 
    260     \caption{   \protect\label{fig:triad:qcells} 
     261    \caption{   \protect\label{fig:qcells} 
    261262    Triad notation for quarter cells. $T$-cells are inside 
    262263      boxes, while the  $i+\half,k$ $u$-cell is shaded in green and the 
    … …  
    265266% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    266267 
    267 Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter 
     268Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter 
    268269cell that is the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ $u$-cell and the $i,k+k_p$ $w$-cell.  
    269 Expressing the slopes $s_i$ and $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation,  
     270Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation,  
    270271we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$.  
    271272Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$)  
    … …  
    276277of a unique triad, and we notate these areas, similarly to the triad slopes,  
    277278as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$,  
    278 where $e.g.$ in \eqref{eq:triad:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,  
    279 and in \eqref{eq:triad:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 
     279where $e.g.$ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,  
     280and in \autoref{eq:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 
    280281 
    281282\subsection{Full triad fluxes} 
    … …  
    287288form 
    288289\begin{equation} 
    289   \label{eq:triad:i11} 
     290  \label{eq:i11} 
    290291  \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} = 
    291292  - \left( \Alts_i^{k+1} a_{1} + \Alts_i^{k+1} a_{2} + \Alts_i^k 
    … …  
    293294  \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}}, 
    294295\end{equation} 
    295 where the areas $a_i$ are as in \eqref{eq:triad:i13}. In this case, 
    296 separating the total lateral flux, the sum of \eqref{eq:triad:i13} and 
    297 \eqref{eq:triad:i11}, into triad components, a lateral tracer 
     296where the areas $a_i$ are as in \autoref{eq:i13}. In this case, 
     297separating the total lateral flux, the sum of \autoref{eq:i13} and 
     298\autoref{eq:i11}, into triad components, a lateral tracer 
    298299flux 
    299300\begin{equation} 
    300   \label{eq:triad:latflux-triad} 
     301  \label{eq:latflux-triad} 
    301302  _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p} 
    302303  \left( 
    … …  
    312313density flux associated with each triad separately disappears. 
    313314\begin{equation} 
    314   \label{eq:triad:latflux-rho} 
     315  \label{eq:latflux-rho} 
    315316  {\mathbb{F}_u}_{i_p}^{k_p} (\rho)=-\alpha _i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (S)=0 
    316317\end{equation} 
    … …  
    319320to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. 
    320321 
    321 The squared slope $r_1^2$ in the expression \eqref{eq:triad:i33c} for the 
     322The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the 
    322323$_{33}$ component is also expressed in terms of area-weighted 
    323324squared triad slopes, so the area-integrated vertical flux from tracer 
    324325cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is 
    325326\begin{equation} 
    326   \label{eq:triad:i33} 
     327  \label{eq:i33} 
    327328  \left( F_w^{33} \right) _i^{k+\frac{1}{2}} = 
    328329    - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2 
    … …  
    332333\end{equation} 
    333334where the areas $a'$ and slopes $s'$ are the same as in 
    334 \eqref{eq:triad:i31}. 
    335 Then, separating the total vertical flux, the sum of \eqref{eq:triad:i31} and 
    336 \eqref{eq:triad:i33}, into triad components,  a vertical flux 
     335\autoref{eq:i31}. 
     336Then, separating the total vertical flux, the sum of \autoref{eq:i31} and 
     337\autoref{eq:i33}, into triad components,  a vertical flux 
    337338\begin{align} 
    338   \label{eq:triad:vertflux-triad} 
     339  \label{eq:vertflux-triad} 
    339340  _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) 
    340341  &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} 
    … …  
    345346  \right) \\ 
    346347  &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right) 
    347    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:triad:vertflux-triad2} 
     348   {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2} 
    348349\end{align} 
    349350may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ 
    … …  
    355356fluxes. 
    356357 
    357 We can explicitly identify (Fig.~\ref{fig:triad:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of 
     358We can explicitly identify (\autoref{fig:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of 
    358359the $u$-fluxes and $w$-fluxes in 
    359 \eqref{eq:triad:i31}, \eqref{eq:triad:i13}, \eqref{eq:triad:i11} \eqref{eq:triad:i33} and 
    360 Fig.~\ref{fig:triad:ISO_triad} to  write out the iso-neutral fluxes at $u$- and 
     360\autoref{eq:i31}, \autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and 
     361\autoref{fig:ISO_triad} to  write out the iso-neutral fluxes at $u$- and 
    361362$w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: 
    362 %(Fig.~\ref{Fig_ISO_triad}): 
    363 \begin{flalign} \label{Eq_iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
     363%(\autoref{fig:ISO_triad}): 
     364\begin{flalign} \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
    364365  \sum_{\substack{i_p,\,k_p}} 
    365366  \begin{pmatrix} 
    … …  
    371372 
    372373\subsection{Ensuring the scheme does not increase tracer variance} 
    373 \label{sec:triad:variance} 
     374\label{subsec:variance} 
    374375 
    375376We now require that this operator should not increase the 
    … …  
    397398  &= -T_{i+i_p-1/2}^k{\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \quad + \quad  T_{i+i_p+1/2}^k 
    398399  {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ 
    399   &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:triad:dvar_iso_i} 
     400  &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:dvar_iso_i} 
    400401 \end{aligned} 
    401402\end{multline} 
    … …  
    404405$i,k+k_p+\half$ (below) of 
    405406\begin{equation} 
    406 \label{eq:triad:dvar_iso_k} 
     407\label{eq:dvar_iso_k} 
    407408  _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    408409\end{equation} 
    409410The total variance tendency driven by the triad is the sum of these 
    410411two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and 
    411 $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:triad:latflux-triad} and 
    412 \eqref{eq:triad:vertflux-triad}, it is 
     412$_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \autoref{eq:latflux-triad} and 
     413\autoref{eq:vertflux-triad}, it is 
    413414\begin{multline*} 
    414415-\Alts_i^k\left \{ 
    … …  
    430431to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by 
    431432\begin{equation} 
    432   \label{eq:triad:V-A} 
     433  \label{eq:V-A} 
    433434  _i^k\mathbb{V}_{i_p}^{k_p} 
    434435  ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k} 
    … …  
    437438the variance tendency reduces to the perfect square 
    438439\begin{equation} 
    439   \label{eq:triad:perfect-square} 
     440  \label{eq:perfect-square} 
    440441  -\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
    441442  \left( 
    … …  
    445446  \right)^2\leq 0. 
    446447\end{equation} 
    447 Thus, the constraint \eqref{eq:triad:V-A} ensures that the fluxes (\ref{eq:triad:latflux-triad}, \ref{eq:triad:vertflux-triad}) associated 
     448Thus, the constraint \autoref{eq:V-A} ensures that the fluxes (\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated 
    448449with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase 
    449450the net variance. Since the total fluxes are sums of such fluxes from 
    … …  
    452453increase. 
    453454 
    454 The expression \eqref{eq:triad:V-A} can be interpreted as a discretization 
     455The expression \autoref{eq:V-A} can be interpreted as a discretization 
    455456of the global integral 
    456457\begin{equation} 
    457   \label{eq:triad:cts-var} 
     458  \label{eq:cts-var} 
    458459  \frac{\partial}{\partial t}\int\!\half T^2\, dV = 
    459460  \int\!\mathbf{F}\cdot\nabla T\, dV, 
    … …  
    480481cells, defined in terms of the distances between $T$, $u$,$f$ and 
    481482$w$-points. This is the natural discretization of 
    482 \eqref{eq:triad:cts-var}. The \NEMO model, however, operates with scale 
     483\autoref{eq:cts-var}. The \NEMO model, however, operates with scale 
    483484factors instead of grid sizes, and scale factors for the quarter 
    484485cells are not defined. Instead, therefore we simply choose 
    485486\begin{equation} 
    486   \label{eq:triad:V-NEMO} 
     487  \label{eq:V-NEMO} 
    487488  _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k, 
    488489\end{equation} 
    … …  
    492493$i+1,k$ reduces to the classical form 
    493494\begin{equation} 
    494   \label{eq:triad:lat-normal} 
     495  \label{eq:lat-normal} 
    495496-\overline\Alts_{\,i+1/2}^k\; 
    496497\frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    … …  
    500501In fact if the diffusive coefficient is defined at $u$-points, so that 
    501502we employ $\Alts_{i+i_p}^k$ instead of  $\Alts_i^k$ in the definitions of the 
    502 triad fluxes \eqref{eq:triad:latflux-triad} and \eqref{eq:triad:vertflux-triad}, 
     503triad fluxes \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, 
    503504we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above. 
    504505 
    … …  
    506507The iso-neutral fluxes at $u$- and 
    507508$w$-points are the sums of the triad fluxes that cross the $u$- and 
    508 $w$-faces \eqref{Eq_iso_flux}: 
    509 \begin{subequations}\label{eq:triad:alltriadflux} 
    510   \begin{flalign}\label{eq:triad:vect_isoflux} 
     509$w$-faces \autoref{eq:iso_flux}: 
     510\begin{subequations}\label{eq:alltriadflux} 
     511  \begin{flalign}\label{eq:vect_isoflux} 
    511512    \vect{F}_{\mathrm{iso}}(T) &\equiv 
    512513    \sum_{\substack{i_p,\,k_p}} 
    … …  
    517518    \end{pmatrix}, 
    518519  \end{flalign} 
    519   where \eqref{eq:triad:latflux-triad}: 
     520  where \autoref{eq:latflux-triad}: 
    520521  \begin{align} 
    521     \label{eq:triad:triadfluxu} 
     522    \label{eq:triadfluxu} 
    522523    _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{ 
    523524      \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    … …  
    534535      -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ 
    535536      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    536     \right),\label{eq:triad:triadfluxw} 
     537    \right),\label{eq:triadfluxw} 
    537538  \end{align} 
    538   with \eqref{eq:triad:V-NEMO} 
     539  with \autoref{eq:V-NEMO} 
    539540  \begin{equation} 
    540     \label{eq:triad:V-NEMO2} 
     541    \label{eq:V-NEMO2} 
    541542    _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k. 
    542543  \end{equation} 
    543544\end{subequations} 
    544545 
    545  The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 
     546 The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 
    546547each tracer point: 
    547 \begin{equation} \label{eq:triad:iso_operator} D_l^T = \frac{1}{b_T} 
     548\begin{equation} \label{eq:iso_operator} D_l^T = \frac{1}{b_T} 
    548549  \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k 
    549550        {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[ 
    … …  
    554555\begin{description} 
    555556\item[$\bullet$ horizontal diffusion] The discretization of the 
    556   diffusion operator recovers \eqref{eq:triad:lat-normal} the traditional five-point Laplacian in 
     557  diffusion operator recovers \autoref{eq:lat-normal} the traditional five-point Laplacian in 
    557558  the limit of flat iso-neutral direction : 
    558   \begin{equation} \label{eq:triad:iso_property0} D_l^T = \frac{1}{b_T} \ 
     559  \begin{equation} \label{eq:iso_property0} D_l^T = \frac{1}{b_T} \ 
    559560    \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; 
    560561      \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad 
    … …  
    564565\item[$\bullet$ implicit treatment in the vertical] Only tracer values 
    565566  associated with a single water column appear in the expression 
    566   \eqref{eq:triad:i33} for the $_{33}$ fluxes, vertical fluxes driven by 
     567  \autoref{eq:i33} for the $_{33}$ fluxes, vertical fluxes driven by 
    567568  vertical gradients. This is of paramount importance since it means 
    568569  that a time-implicit algorithm can be used to solve the vertical 
    … …  
    582583\item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of 
    583584  locally referenced potential density is zero. See 
    584   \eqref{eq:triad:latflux-rho} and \eqref{eq:triad:vertflux-triad2}. 
     585  \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}. 
    585586 
    586587\item[$\bullet$ conservation of tracer] The iso-neutral diffusion 
    587588  conserves tracer content, $i.e.$ 
    588   \begin{equation} \label{eq:triad:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ 
     589  \begin{equation} \label{eq:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ 
    589590      b_T \right\} = 0 
    590591  \end{equation} 
    … …  
    594595\item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion 
    595596  does not increase the tracer variance, $i.e.$ 
    596   \begin{equation} \label{eq:triad:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T 
     597  \begin{equation} \label{eq:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T 
    597598      \ b_T \right\} \leq 0 
    598599  \end{equation} 
    599600  The property is demonstrated in 
    600   \S\ref{sec:triad:variance} above. It is a key property for a diffusion 
     601  \autoref{subsec:variance} above. It is a key property for a diffusion 
    601602  term. It means that it is also a dissipation term, $i.e.$ it 
    602603  dissipates the square of the quantity on which it is applied.  It 
    … …  
    607608\item[$\bullet$ self-adjoint operator] The iso-neutral diffusion 
    608609  operator is self-adjoint, $i.e.$ 
    609   \begin{equation} \label{eq:triad:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T 
     610  \begin{equation} \label{eq:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T 
    610611      \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
    611612  \end{equation} 
    … …  
    614615  routine. This property can be demonstrated similarly to the proof of 
    615616  the `no increase of tracer variance' property. The contribution by a 
    616   single triad towards the left hand side of \eqref{eq:triad:iso_property3}, can 
    617   be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{eq:triad:dvar_iso_i} 
    618   and \eqref{eq:triad:dvar_iso_k}. This results in a term similar to 
    619   \eqref{eq:triad:perfect-square}, 
    620 \begin{equation} 
    621   \label{eq:triad:TScovar} 
     617  single triad towards the left hand side of \autoref{eq:iso_property3}, can 
     618  be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} 
     619  and \autoref{eq:dvar_iso_k}. This results in a term similar to 
     620  \autoref{eq:perfect-square}, 
     621\begin{equation} 
     622  \label{eq:TScovar} 
    622623  - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
    623624  \left( 
    … …  
    634635This is symmetrical in $T $ and $S$, so exactly the same term arises 
    635636from the discretization of this triad's contribution towards the 
    636 RHS of \eqref{eq:triad:iso_property3}. 
     637RHS of \autoref{eq:iso_property3}. 
    637638\end{description} 
    638639 
    639 \subsection{Treatment of the triads at the boundaries}\label{sec:triad:iso_bdry} 
     640\subsection{Treatment of the triads at the boundaries}\label{sec:iso_bdry} 
    640641The triad slope can only be defined where both the grid boxes centred at 
    641642the end of the arms exist. Triads that would poke up 
    642643through the upper ocean surface into the atmosphere, or down into the 
    643 ocean floor, must be masked out. See Fig.~\ref{fig:triad:bdry_triads}. Surface layer triads 
     644ocean floor, must be masked out. See \autoref{fig:bdry_triads}. Surface layer triads 
    644645$\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and 
    645646$\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be 
    646 specified above the ocean surface are masked (Fig.~\ref{fig:triad:bdry_triads}a): this ensures that lateral 
     647specified above the ocean surface are masked (\autoref{fig:bdry_triads}a): this ensures that lateral 
    647648tracer gradients produce no flux through the ocean surface. However, 
    648649to prevent surface noise, it is customary to retain the $_{11}$ contributions towards 
    … …  
    650651$\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer 
    651652fluxes. Similar comments apply to triads that would intersect the 
    652 ocean floor (Fig.~\ref{fig:triad:bdry_triads}b). Note that both near bottom 
     653ocean floor (\autoref{fig:bdry_triads}b). Note that both near bottom 
    653654triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 
    654655$\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ 
    … …  
    665666\begin{figure}[h] \begin{center} 
    666667    \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads} 
    667     \caption{  \protect\label{fig:triad:bdry_triads} 
     668    \caption{  \protect\label{fig:bdry_triads} 
    668669      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer 
    669670      points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad 
    … …  
    678679      or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point 
    679680      is masked. The associated lateral fluxes (grey-black dashed 
    680       line) are masked if \np{botmix\_triad}\forcode{ = .false.}, but left 
    681       unmasked, giving bottom mixing, if \np{botmix\_triad}\forcode{ = .true.}} 
     681      line) are masked if \protect\np{botmix\_triad}\forcode{ = .false.}, but left 
     682      unmasked, giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.}} 
    682683 \end{center} \end{figure} 
    683684% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    684685 
    685 \subsection{ Limiting of the slopes within the interior}\label{sec:triad:limit} 
    686 As discussed in \S\ref{LDF_slp_iso}, iso-neutral slopes relative to 
     686\subsection{ Limiting of the slopes within the interior}\label{sec:limit} 
     687As discussed in \autoref{subsec:LDF_slp_iso}, iso-neutral slopes relative to 
    687688geopotentials must be bounded everywhere, both for consistency with the small-slope 
    688689approximation and for numerical stability \citep{Cox1987, 
    … …  
    692693It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to 
    693694geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to 
    694 geopotentials) \eqref{Eq_PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate 
     695geopotentials) \autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate 
    695696surfaces, so we require 
    696697\begin{equation*} 
    … …  
    700701Each individual triad slope 
    701702 \begin{equation} 
    702    \label{eq:triad:Rtilde} 
     703   \label{eq:Rtilde} 
    703704_i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p}  + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    704705 \end{equation} 
    … …  
    709710is always downwards, and so acts to reduce gravitational potential energy. 
    710711 
    711 \subsection{Tapering within the surface mixed layer}\label{sec:triad:taper} 
     712\subsection{Tapering within the surface mixed layer}\label{sec:taper} 
    712713Additional tapering of the iso-neutral fluxes is necessary within the 
    713714surface mixed layer. When the Griffies triads are used, we offer two 
    714715options for this. 
    715716 
    716 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper} 
     717\subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:lintaper} 
    717718This is the option activated by the default choice 
    718719\np{ln\_triad\_iso}\forcode{ = .false.}. Slopes $\tilde{r}_i$ relative to 
    719720geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the 
    720 surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values 
     721surface, as described in option (c) of \autoref{fig:eiv_slp}, to values 
    721722\begin{subequations} 
    722723  \begin{equation} 
    723    \label{eq:triad:rmtilde} 
     724   \label{eq:rmtilde} 
    724725     \rMLt = 
    725726  -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h, 
    … …  
    728729adjusted to 
    729730  \begin{equation} 
    730    \label{eq:triad:rm} 
     731   \label{eq:rm} 
    731732 \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h. 
    732733  \end{equation} 
    733734\end{subequations} 
    734735Thus the diffusion operator within the mixed layer is given by: 
    735 \begin{equation} \label{eq:triad:iso_tensor_ML} 
     736\begin{equation} \label{eq:iso_tensor_ML} 
    736737D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    737738\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 
    … …  
    745746mixed-layer and in isopycnal layers immediately below, in the 
    746747thermocline. It is consistent with the way the $\tilde{r}_i$ are 
    747 tapered within the mixed layer (see \S\ref{sec:triad:taperskew} below) 
     748tapered within the mixed layer (see \autoref{sec:taperskew} below) 
    748749so as to ensure a uniform GM eddy-induced velocity throughout the 
    749750mixed layer. However, it gives a downwards density flux and so acts so 
    750751as to reduce potential energy in the same way as does the slope 
    751 limiting discussed above in \S\ref{sec:triad:limit}. 
     752limiting discussed above in \autoref{sec:limit}. 
    752753  
    753 As in \S\ref{sec:triad:limit} above, the tapering 
    754 \eqref{eq:triad:rmtilde} is applied separately to each triad 
     754As in \autoref{sec:limit} above, the tapering 
     755\autoref{eq:rmtilde} is applied separately to each triad 
    755756$_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the 
    756757$_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume 
    757758$z$-coordinates in the following; the conversion from 
    758759$\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described 
    759 above by \eqref{eq:triad:Rtilde}. 
     760above by \autoref{eq:Rtilde}. 
    760761\begin{enumerate} 
    761762\item Mixed-layer depth is defined so as to avoid including regions of weak 
    762763vertical stratification in the slope definition. 
    763764 At each $i,j$ (simplified to $i$ in 
    764 Fig.~\ref{fig:triad:MLB_triad}), we define the mixed-layer by setting 
     765\autoref{fig:MLB_triad}), we define the mixed-layer by setting 
    765766the vertical index of the tracer point immediately below the mixed 
    766767layer, $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) 
    … …  
    768769${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is 
    769770the tracer gridbox within which the depth reaches 10~m. See the left 
    770 side of Fig.~\ref{fig:triad:MLB_triad}. We use the $k_{10}$-gridbox 
     771side of \autoref{fig:MLB_triad}. We use the $k_{10}$-gridbox 
    771772instead of the surface gridbox to avoid problems e.g.\ with thin 
    772773daytime mixed-layers. Currently we use the same 
    … …  
    784785${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are 
    785786representative of the thermocline. The four basal triads defined in the bottom part 
    786 of Fig.~\ref{fig:triad:MLB_triad} are then 
     787of \autoref{fig:MLB_triad} are then 
    787788\begin{align} 
    788789  {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &= 
    789  {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}, \label{eq:triad:Rbase} 
     790 {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}, \label{eq:Rbase} 
    790791\\ 
    791792\intertext{with e.g.\ the green triad} 
    … …  
    797798so it is this depth 
    798799\begin{equation} 
    799   \label{eq:triad:zbase} 
     800  \label{eq:zbase} 
    800801  {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2} 
    801802\end{equation} 
    802803(one gridbox deeper than the 
    803804diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper 
    804 the slopes in \eqref{eq:triad:rmtilde}. 
     805the slopes in \autoref{eq:rmtilde}. 
    805806\item Finally, we calculate the adjusted triads 
    806807${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within the mixed 
    … …  
    815816\intertext{and more generally} 
    816817 {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p} &= 
    817 \frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}.\label{eq:triad:RML} 
     818\frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}.\label{eq:RML} 
    818819\end{align} 
    819820\end{enumerate} 
    … …  
    822823\begin{figure}[h] 
    823824%  \fcapside { 
    824     \caption{\protect\label{fig:triad:MLB_triad} Definition of 
     825    \caption{\protect\label{fig:MLB_triad} Definition of 
    825826      mixed-layer depth and calculation of linearly tapered 
    826827      triads. The figure shows a water column at a given $i,j$ 
    … …  
    846847 
    847848\subsubsection{Additional truncation of skew iso-neutral flux components} 
    848 \label{sec:triad:Gerdes-taper} 
     849\label{subsec:Gerdes-taper} 
    849850The alternative option is activated by setting \np{ln\_triad\_iso} = 
    850851  true. This retains the same tapered slope $\rML$  described above for the 
    … …  
    853854replaces the $\rML$ in the skew term by 
    854855\begin{equation} 
    855   \label{eq:triad:rm*} 
     856  \label{eq:rm*} 
    856857  \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i, 
    857858\end{equation} 
    858859giving a ML diffusive operator 
    859 \begin{equation} \label{eq:triad:iso_tensor_ML2} 
     860\begin{equation} \label{eq:iso_tensor_ML2} 
    860861D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    861862\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 
    … …  
    887888% Skew flux formulation for Eddy Induced Velocity : 
    888889% ================================================================ 
    889 \section{Eddy induced advection formulated as a skew flux}\label{sec:triad:skew-flux} 
    890  
    891 \subsection{Continuous skew flux formulation}\label{sec:triad:continuous-skew-flux} 
     890\section{Eddy induced advection formulated as a skew flux}\label{sec:skew-flux} 
     891 
     892\subsection{Continuous skew flux formulation}\label{sec:continuous-skew-flux} 
    892893 
    893894 When Gent and McWilliams's [1990] diffusion is used, 
    … …  
    895896eddy induced velocity, the formulation of which depends on the slopes of iso- 
    896897neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used 
    897 here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} 
    898 is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo} 
    899 + \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates. 
     898here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo} 
     899is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} 
     900+ \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. 
    900901 
    901902The eddy induced velocity is given by: 
    902 \begin{subequations} \label{eq:triad:eiv} 
    903 \begin{equation}\label{eq:triad:eiv_v} 
     903\begin{subequations} \label{eq:eiv} 
     904\begin{equation}\label{eq:eiv_v} 
    904905\begin{split} 
    905906 u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\ 
    … …  
    910911\end{equation} 
    911912where the streamfunctions $\psi_i$ are given by 
    912 \begin{equation} \label{eq:triad:eiv_psi} 
     913\begin{equation} \label{eq:eiv_psi} 
    913914\begin{split} 
    914915\psi_1 & = A_{e} \; \tilde{r}_1,   \\ 
    … …  
    924925default implementation, where \np{ln\_traldf\_triad} is set 
    925926false. This allows us to take advantage of all the advection schemes 
    926 offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$ 
     927offered for the tracers (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ 
    927928order advection scheme. This is particularly useful for passive 
    928929tracers where \emph{positivity} of the advection scheme is of 
    … …  
    962963and since the eddy induced velocity field is non-divergent, we end up with the skew 
    963964form of the eddy induced advective fluxes per unit area in $ijk$ space: 
    964 \begin{equation} \label{eq:triad:eiv_skew_ijk} 
     965\begin{equation} \label{eq:eiv_skew_ijk} 
    965966\textbf{F}_\mathrm{eiv}^T = \begin{pmatrix} 
    966967           {+ e_{2} \, \psi_1  \; \partial_k T}   \\ 
    … …  
    969970\end{equation} 
    970971The total fluxes per unit physical area are then 
    971 \begin{equation}\label{eq:triad:eiv_skew_physical} 
     972\begin{equation}\label{eq:eiv_skew_physical} 
    972973\begin{split} 
    973974 f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\ 
    … …  
    977978\end{split} 
    978979\end{equation} 
    979 Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the 
     980Note that  \autoref{eq:eiv_skew_physical} takes the same form whatever the 
    980981vertical coordinate, though of course the slopes 
    981 $\tilde{r}_i$ which define the $\psi_i$ in \eqref{eq:triad:eiv_psi} are relative to geopotentials. 
     982$\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative to geopotentials. 
    982983The tendency associated with eddy induced velocity is then simply the convergence 
    983 of the fluxes (\ref{eq:triad:eiv_skew_ijk}, \ref{eq:triad:eiv_skew_physical}), so 
    984 \begin{equation} \label{eq:triad:skew_eiv_conv} 
     984of the fluxes (\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so 
     985\begin{equation} \label{eq:skew_eiv_conv} 
    985986\frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[ 
    986987  \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right) 
    … …  
    995996 
    996997\subsection{Discrete skew flux formulation} 
    997 The skew fluxes in (\ref{eq:triad:eiv_skew_physical}, \ref{eq:triad:eiv_skew_ijk}), like the off-diagonal terms 
    998 (\ref{eq:triad:i13c}, \ref{eq:triad:i31c}) of the small angle diffusion tensor, are best 
    999 expressed in terms of the triad slopes, as in Fig.~\ref{fig:triad:ISO_triad} 
    1000 and Eqs~(\ref{eq:triad:i13}, \ref{eq:triad:i31}); but now in terms of the triad slopes 
     998The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}), like the off-diagonal terms 
     999(\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor, are best 
     1000expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad} 
     1001and (\autoref{eq:i13}, \autoref{eq:i31}); but now in terms of the triad slopes 
    10011002$\tilde{\mathbb{R}}$ relative to geopotentials instead of the 
    10021003$\mathbb{R}$ relative to coordinate surfaces. The discrete form of 
    1003 \eqref{eq:triad:eiv_skew_ijk} using the slopes \eqref{eq:triad:R} and 
     1004\autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and 
    10041005defining $A_e$ at $T$-points is then given by: 
    10051006 
    10061007 
    1007 \begin{subequations}\label{eq:triad:allskewflux} 
    1008   \begin{flalign}\label{eq:triad:vect_skew_flux} 
     1008\begin{subequations}\label{eq:allskewflux} 
     1009  \begin{flalign}\label{eq:vect_skew_flux} 
    10091010    \vect{F}_{\mathrm{eiv}}(T) &\equiv 
    10101011    \sum_{\substack{i_p,\,k_p}} 
    … …  
    10161017  \end{flalign} 
    10171018  where the skew flux in the $i$-direction associated with a given 
    1018   triad is (\ref{eq:triad:latflux-triad}, \ref{eq:triad:triadfluxu}): 
     1019  triad is (\autoref{eq:latflux-triad}, \autoref{eq:triadfluxu}): 
    10191020  \begin{align} 
    1020     \label{eq:triad:skewfluxu} 
     1021    \label{eq:skewfluxu} 
    10211022    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{ 
    10221023      \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    … …  
    10241025      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, 
    10251026   \\ 
    1026     \intertext{and \eqref{eq:triad:triadfluxw} in the $k$-direction, changing the sign 
    1027       to be consistent with \eqref{eq:triad:eiv_skew_ijk}:} 
     1027    \intertext{and \autoref{eq:triadfluxw} in the $k$-direction, changing the sign 
     1028      to be consistent with \autoref{eq:eiv_skew_ijk}:} 
    10281029    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T) 
    10291030    &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}} 
    1030      {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:triad:skewfluxw} 
     1031     {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:skewfluxw} 
    10311032  \end{align} 
    10321033\end{subequations} 
    … …  
    10401041include a diffusive component but is a `pure' advection term. This can 
    10411042be seen 
    1042 %either from Appendix \ref{Apdx_eiv_skew} or 
     1043%either from Appendix \autoref{apdx:eiv_skew} or 
    10431044by considering the 
    10441045fluxes associated with a given triad slope 
    10451046$_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following 
    1046 \S\ref{sec:triad:variance} and \eqref{eq:triad:dvar_iso_i}, the 
     1047\autoref{subsec:variance} and \autoref{eq:dvar_iso_i}, the 
    10471048associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ 
    10481049drives a net rate of change of variance, summed over the two 
    10491050$T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 
    10501051\begin{equation} 
    1051 \label{eq:triad:dvar_eiv_i} 
     1052\label{eq:dvar_eiv_i} 
    10521053  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], 
    10531054\end{equation} 
    … …  
    10551056$T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of 
    10561057\begin{equation} 
    1057 \label{eq:triad:dvar_eiv_k} 
     1058\label{eq:dvar_eiv_k} 
    10581059  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    10591060\end{equation} 
    1060 Inspection of the definitions (\ref{eq:triad:skewfluxu}, \ref{eq:triad:skewfluxw}) 
    1061 shows that these two variance changes (\ref{eq:triad:dvar_eiv_i}, \ref{eq:triad:dvar_eiv_k}) 
     1061Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw}) 
     1062shows that these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k}) 
    10621063sum to zero. Hence the two fluxes associated with each triad make no 
    10631064net contribution to the variance budget. 
    … …  
    10721073For the change in gravitational PE driven by the $k$-flux is 
    10731074\begin{align} 
    1074   \label{eq:triad:vert_densityPE} 
     1075  \label{eq:vert_densityPE} 
    10751076  g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) 
    10761077  &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k 
    … …  
    10781079    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\ 
    10791080\intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from 
    1080   \eqref{eq:triad:skewfluxw}, gives} 
     1081  \autoref{eq:skewfluxw}, gives} 
    10811082% and separating out 
    10821083% $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$, 
    … …  
    10901091\end{align} 
    10911092using the definition of the triad slope $\rtriad{R}$, 
    1092 \eqref{eq:triad:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ 
     1093\autoref{eq:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ 
    10931094\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of  $-\alpha_i^k \delta_{k+ 
    10941095  k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$. 
    … …  
    10961097Where the coordinates slope, the $i$-flux gives a PE change 
    10971098\begin{multline} 
    1098   \label{eq:triad:lat_densityPE} 
     1099  \label{eq:lat_densityPE} 
    10991100 g \delta_{i+i_p}[z_T^k] 
    11001101\left[ 
    … …  
    11061107\frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}, 
    11071108\end{multline} 
    1108 (using \eqref{eq:triad:skewfluxu}) and so the total PE change 
    1109 \eqref{eq:triad:vert_densityPE} + \eqref{eq:triad:lat_densityPE} associated with the triad fluxes is 
     1109(using \autoref{eq:skewfluxu}) and so the total PE change 
     1110\autoref{eq:vert_densityPE} + \autoref{eq:lat_densityPE} associated with the triad fluxes is 
    11101111\begin{multline} 
    1111   \label{eq:triad:tot_densityPE} 
     1112  \label{eq:tot_densityPE} 
    11121113  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) + 
    11131114g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\ 
    … …  
    11191120\beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative. 
    11201121 
    1121 \subsection{Treatment of the triads at the boundaries}\label{sec:triad:skew_bdry} 
     1122\subsection{Treatment of the triads at the boundaries}\label{sec:skew_bdry} 
    11221123Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes 
    11231124are masked at the boundaries in exactly the same way as are the triad 
    11241125slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as 
    1125 described in \S\ref{sec:triad:iso_bdry} and 
    1126 Fig.~\ref{fig:triad:bdry_triads}. Thus surface layer triads 
     1126described in \autoref{sec:iso_bdry} and 
     1127\autoref{fig:bdry_triads}. Thus surface layer triads 
    11271128$\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are 
    11281129masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ 
    … …  
    11321133no effect on the eddy-induced skew-fluxes. 
    11331134 
    1134 \subsection{Limiting of the slopes within the interior}\label{sec:triad:limitskew} 
     1135\subsection{Limiting of the slopes within the interior}\label{sec:limitskew} 
    11351136Presently, the iso-neutral slopes $\tilde{r}_i$ relative 
    11361137to geopotentials are limited to be less than $1/100$, exactly as in 
    1137 calculating the iso-neutral diffusion, \S \ref{sec:triad:limit}. Each 
     1138calculating the iso-neutral diffusion, \S \autoref{sec:limit}. Each 
    11381139individual triad \rtriadt{R} is so limited. 
    11391140 
    1140 \subsection{Tapering within the surface mixed layer}\label{sec:triad:taperskew} 
     1141\subsection{Tapering within the surface mixed layer}\label{sec:taperskew} 
    11411142The slopes $\tilde{r}_i$ relative to 
    11421143geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the 
    1143 surface \eqref{eq:triad:rmtilde}, as described in \S\ref{sec:triad:lintaper}. This is 
    1144 option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the 
     1144surface \autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}. This is 
     1145option (c) of \autoref{fig:eiv_slp}. This linear tapering for the 
    11451146slopes used to calculate the eddy-induced fluxes is 
    11461147unaffected by the value of \np{ln\_triad\_iso}. 
    … …  
    11481149The justification for this linear slope tapering is that, for $A_e$ 
    11491150that is constant or varies only in the horizontal (the most commonly 
    1150 used options in \NEMO: see \S\ref{LDF_coef}), it is 
     1151used options in \NEMO: see \autoref{sec:LDF_coef}), it is 
    11511152equivalent to a horizontal eiv (eddy-induced velocity) that is uniform 
    1152 within the mixed layer \eqref{eq:triad:eiv_v}. This ensures that the 
     1153within the mixed layer \autoref{eq:eiv_v}. This ensures that the 
    11531154eiv velocities do not restratify the mixed layer \citep{Treguier1997, 
    11541155  Danabasoglu_al_2008}. Equivantly, in terms 
    … …  
    11581159horizontal flux convergence is relatively insignificant within the mixed-layer). 
    11591160 
    1160 \subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag} 
     1161\subsection{Streamfunction diagnostics}\label{sec:sfdiag} 
    11611162Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, diagnosed 
    11621163mean eddy-induced velocities are output. Each time step, 
    … …  
    11641165$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ 
    11651166(integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table 
    1166 \ref{Tab_cell}) respectively. We follow \citep{Griffies_Bk04} and 
     1167\autoref{tab:cell}) respectively. We follow \citep{Griffies_Bk04} and 
    11671168calculate the streamfunction at a given $uw$-point from the 
    11681169surrounding four triads according to: 
    11691170\begin{equation} 
    1170   \label{eq:triad:sfdiagi} 
     1171  \label{eq:sfdiagi} 
    11711172  {\psi_1}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}} 
    11721173  {A_e}_{i+1/2-i_p}^{k+1/2-k_p}\:\triadd{i+1/2-i_p}{k+1/2-k_p}{R}{i_p}{k_p}. 
    … …  
    11741175The streamfunction $\psi_1$ is calculated similarly at $vw$ points. 
    11751176The eddy-induced velocities are then calculated from the 
    1176 straightforward discretisation of \eqref{eq:triad:eiv_v}: 
    1177 \begin{equation}\label{eq:triad:eiv_v_discrete} 
     1177straightforward discretisation of \autoref{eq:eiv_v}: 
     1178\begin{equation}\label{eq:eiv_v_discrete} 
    11781179\begin{split} 
    11791180 {u^*}_{i+1/2}^{k} & = - \frac{1}{{e_{3u}}_{i}^{k}}\left({\psi_1}_{i+1/2}^{k+1/2}-{\psi_1}_{i+1/2}^{k+1/2}\right),   \\ 

  • branches/2017/dev_merge_2017/DOC/tex_sub/chap_ASM.tex

    r9393 r9407  
    55% ================================================================ 
    66\chapter{Apply Assimilation Increments (ASM)} 
    7 \label{ASM} 
     7\label{chap:ASM} 
    88 
    99Authors: D. Lea,  M. Martin, K. Mogensen, A. Weaver, ...   % do we keep 
    … …  
    2626 
    2727\section{Direct initialization} 
    28 \label{ASM_DI} 
     28\label{sec:ASM_DI} 
    2929 
    3030Direct initialization (DI) refers to the instantaneous correction 
    … …  
    3333 
    3434\section{Incremental analysis updates} 
    35 \label{ASM_IAU} 
     35\label{sec:ASM_IAU} 
    3636 
    3737Rather than updating the model state directly with the analysis increment, 
    … …  
    8383\end{eqnarray} 
    8484where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even.  
    85 The weights described by \eqref{eq:F2_i} provide a  
     85The weights described by \autoref{eq:F2_i} provide a  
    8686smoother transition of the analysis trajectory from one assimilation cycle  
    87 to the next than that described by \eqref{eq:F1_i}. 
     87to the next than that described by \autoref{eq:F1_i}. 
    8888 
    8989%========================================================================== 
    9090% Divergence damping description %%% 
    9191\section{Divergence damping initialisation} 
    92 \label{ASM_details} 
     92\label{sec:ASM_details} 
    9393 
    9494The velocity increments may be initialized by the iterative application of  
    … …  
    110110                       +\delta _j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right). 
    111111\end{equation} 
    112 By the application of \eqref{eq:asm_dmp} and \eqref{eq:asm_dmp} the divergence is filtered 
     112By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered 
    113113in each iteration, and the vorticity is left unchanged. In the presence of coastal boundaries 
    114114with zero velocity increments perpendicular to the coast the divergence is strongly damped. 
    … …  
    125125 
    126126\section{Implementation details} 
    127 \label{ASM_details} 
     127\label{sec:ASM_details} 
    128128 
    129129Here we show an example \ngn{namasm} namelist and the header of an example assimilation  

  • branches/2017/dev_merge_2017/DOC/tex_sub/chap_CONFIG.tex

    r9394 r9407  
    55% ================================================================ 
    66\chapter{Configurations} 
    7 \label{CFG} 
     7\label{chap:CFG} 
    88\minitoc 
    99 
    … …  
    1515% ================================================================ 
    1616\section{Introduction} 
    17 \label{CFG_intro} 
     17\label{sec:CFG_intro} 
    1818 
    1919 
    … …  
    3333% ================================================================ 
    3434\section{C1D: 1D Water column model (\protect\key{c1d}) } 
    35 \label{CFG_c1d} 
     35\label{sec:CFG_c1d} 
    3636 
    3737$\ $\newline 
    … …  
    8181% ================================================================ 
    8282\section{ORCA family: global ocean with tripolar grid} 
    83 \label{CFG_orca} 
     83\label{sec:CFG_orca} 
    8484 
    8585The ORCA family is a series of global ocean configurations that are run together with  
    … …  
    9595\begin{figure}[!t]   \begin{center} 
    9696\includegraphics[width=0.98\textwidth]{Fig_ORCA_NH_mesh} 
    97 \caption{  \protect\label{Fig_MISC_ORCA_msh}      
     97\caption{  \protect\label{fig:MISC_ORCA_msh}      
    9898ORCA mesh conception. The departure from an isotropic Mercator grid start poleward of 20\degN. 
    9999The two "north pole" are the foci of a series of embedded ellipses (blue curves)  
    … …  
    108108% ------------------------------------------------------------------------------------------------------------- 
    109109\subsection{ORCA tripolar grid} 
    110 \label{CFG_orca_grid} 
     110\label{subsec:CFG_orca_grid} 
    111111 
    112112The ORCA grid is a tripolar is based on the semi-analytical method of \citet{Madec_Imbard_CD96}.  
    … …  
    116116computing the associated set of mesh meridians, and projecting the resulting mesh onto the sphere.  
    117117The set of mesh parallels used is a series of embedded ellipses which foci are the two mesh north  
    118 poles (Fig.~\ref{Fig_MISC_ORCA_msh}). The resulting mesh presents no loss of continuity in  
     118poles (\autoref{fig:MISC_ORCA_msh}). The resulting mesh presents no loss of continuity in  
    119119either the mesh lines or the scale factors, or even the scale factor derivatives over the whole  
    120120ocean domain, as the mesh is not a composite mesh.  
    … …  
    123123\includegraphics[width=1.0\textwidth]{Fig_ORCA_NH_msh05_e1_e2} 
    124124\includegraphics[width=0.80\textwidth]{Fig_ORCA_aniso} 
    125 \caption {  \protect\label{Fig_MISC_ORCA_e1e2} 
     125\caption {  \protect\label{fig:MISC_ORCA_e1e2} 
    126126\textit{Top}: Horizontal scale factors ($e_1$, $e_2$) and  
    127127\textit{Bottom}: ratio of anisotropy ($e_1 / e_2$) 
    … …  
    141141the Gulf Stream) and keeping the smallest scale factor in the northern hemisphere larger  
    142142than the smallest one in the southern hemisphere. 
    143 The resulting mesh is shown in Fig.~\ref{Fig_MISC_ORCA_msh} and \ref{Fig_MISC_ORCA_e1e2}  
     143The resulting mesh is shown in \autoref{fig:MISC_ORCA_msh} and \autoref{fig:MISC_ORCA_e1e2}  
    144144for a half a degree grid (ORCA\_R05). 
    145145The smallest ocean scale factor is found in along  Antarctica, while the ratio of anisotropy remains close to one except near the Victoria Island  
    … …  
    150150% ------------------------------------------------------------------------------------------------------------- 
    151151\subsection{ORCA pre-defined resolution} 
    152 \label{CFG_orca_resolution} 
     152\label{subsec:CFG_orca_resolution} 
    153153 
    154154 
    … …  
    156156horizontal resolution. The value of the resolution is given by the resolution at the Equator  
    157157expressed in degrees. Each of configuration is set through the \textit{domain\_cfg} domain configuration file,  
    158 which sets the grid size and configuration name parameters. The NEMO System Team provides only ORCA2 domain input file "\ifile{ORCA\_R2\_zps\_domcfg}" file  (Tab. \ref{Tab_ORCA}). 
     158which sets the grid size and configuration name parameters. The NEMO System Team provides only ORCA2 domain input file "\ifile{ORCA\_R2\_zps\_domcfg}" file  (Tab. \autoref{tab:ORCA}). 
    159159 
    160160 
    … …  
    175175\hline   \hline 
    176176\end{tabular} 
    177 \caption{ \protect\label{Tab_ORCA}    
     177\caption{ \protect\label{tab:ORCA}    
    178178Domain size of ORCA family configurations. 
    179179The flag for configurations of ORCA family need to be set in \textit{domain\_cfg} file. } 
    … …  
    196196For ORCA\_R1 and R025, setting the configuration key to 75 allows to use 75 vertical levels,  
    197197otherwise 46 are used. In the other ORCA configurations, 31 levels are used  
    198 (see Tab.~\ref{Tab_orca_zgr} \sfcomment{HERE I need to put new table for ORCA2 values} and Fig.~\ref{Fig_zgr}). 
     198(see \autoref{tab:orca_zgr} \sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}). 
    199199 
    200200Only the ORCA\_R2 is provided with all its input files in the \NEMO distribution.  
    … …  
    204204 
    205205This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m)  
    206 in the upper 150m (see Tab.~\ref{Tab_orca_zgr} and Fig.~\ref{Fig_zgr}).  
     206in the upper 150m (see \autoref{tab:orca_zgr} and \autoref{fig:zgr}).  
    207207The bottom topography and the coastlines are derived from the global atlas of Smith and Sandwell (1997).  
    208 The default forcing uses the boundary forcing from \citet{Large_Yeager_Rep04} (see \S\ref{SBC_blk_core}),  
     208The default forcing uses the boundary forcing from \citet{Large_Yeager_Rep04} (see \autoref{subsec:SBC_blk_core}),  
    209209which was developed for the purpose of running global coupled ocean-ice simulations  
    210210without an interactive atmosphere. This \citet{Large_Yeager_Rep04} dataset is available  
    … …  
    222222% ------------------------------------------------------------------------------------------------------------- 
    223223\section{GYRE family: double gyre basin } 
    224 \label{CFG_gyre} 
     224\label{sec:CFG_gyre} 
    225225 
    226226The GYRE configuration \citep{Levy_al_OM10} has been built to simulate 
    … …  
    234234The domain geometry is a closed rectangular basin on the $\beta$-plane centred  
    235235at $\sim$ 30\degN and rotated by 45\deg, 3180~km long, 2120~km wide  
    236 and 4~km deep (Fig.~\ref{Fig_MISC_strait_hand}).  
     236and 4~km deep (\autoref{fig:MISC_strait_hand}).  
    237237The domain is bounded by vertical walls and by a flat bottom. The configuration is  
    238238meant to represent an idealized North Atlantic or North Pacific basin.  
    … …  
    257257Obviously, the namelist parameters have to be adjusted to the chosen resolution, see the Configurations  
    258258pages on the NEMO web site (Using NEMO\/Configurations) . 
    259 In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}\forcode{ = 31}) (Fig.~\ref{Fig_zgr}). 
     259In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}\forcode{ = 31}) (\autoref{fig:zgr}). 
    260260 
    261261The GYRE configuration is also used in benchmark test as it is very simple to increase  
    … …  
    267267\begin{figure}[!t]   \begin{center} 
    268268\includegraphics[width=1.0\textwidth]{Fig_GYRE} 
    269 \caption{  \protect\label{Fig_GYRE}    
     269\caption{  \protect\label{fig:GYRE}    
    270270Snapshot of relative vorticity at the surface of the model domain  
    271271in GYRE R9, R27 and R54. From \citet{Levy_al_OM10}.} 
    … …  
    277277% ------------------------------------------------------------------------------------------------------------- 
    278278\section{AMM: atlantic margin configuration} 
    279 \label{MISC_config_AMM} 
     279\label{sec:MISC_config_AMM} 
    280280 
    281281The AMM, Atlantic Margins Model, is a regional model covering the 

  • branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIA.tex

    r9394 r9407  
    55% ================================================================ 
    66\chapter{Output and Diagnostics (IOM, DIA, TRD, FLO)} 
    7 \label{DIA} 
     7\label{chap:DIA} 
    88\minitoc 
    99 
    … …  
    1515% ================================================================ 
    1616\section{Old model output (default)} 
    17 \label{DIA_io_old} 
     17\label{sec:DIA_io_old} 
    1818 
    1919The model outputs are of three types: the restart file, the output listing,  
    … …  
    5656% ================================================================ 
    5757\section{Standard model output (IOM)} 
    58 \label{DIA_iom} 
     58\label{sec:DIA_iom} 
    5959 
    6060 
    … …  
    595595 
    596596\subsection{XML reference tables} 
    597 \label{IOM_xmlref} 
     597\label{subsec:IOM_xmlref} 
    598598 
    599599(1) Simple computation: directly define the computation when refering to the variable in the file definition. 
    … …  
    998998\subsection{CF metadata standard compliance} 
    999999 
    1000 Output from the XIOS-1.0 IO server is compliant with \href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of the CF metadata standard. Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of section \ref{IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. 
     1000Output from the XIOS-1.0 IO server is compliant with \href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of the CF metadata standard. Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of section \autoref{subsec:IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. 
    10011001 
    10021002Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by the namelist parameter \np{ln\_cfmeta} in the \ngn{namrun} namelist. This must be set to true if these metadata are to be included in the output files. 
    … …  
    10071007% ================================================================ 
    10081008\section{NetCDF4 support (\protect\key{netcdf4})} 
    1009 \label{DIA_iom} 
     1009\label{sec:DIA_iom} 
    10101010 
    10111011Since version 3.3, support for NetCDF4 chunking and (loss-less) compression has 
    … …  
    10701070respectively in the mono-processor case (i.e. global domain of {\small\tt 182x149x31}). 
    10711071An illustration of the potential space savings that NetCDF4 chunking and compression 
    1072 provides is given in table \ref{Tab_NC4} which compares the results of two short 
     1072provides is given in table \autoref{tab:NC4} which compares the results of two short 
    10731073runs of the ORCA2\_LIM reference configuration with a 4x2 mpi partitioning. Note 
    10741074the variation in the compression ratio achieved which reflects chiefly the dry to wet  
    … …  
    11061106ORCA2\_2d\_grid\_W\_0007.nc & 4416 & 1368 & 70\%\\ 
    11071107\end{tabular} 
    1108 \caption{   \protect\label{Tab_NC4}  
     1108\caption{   \protect\label{tab:NC4}  
    11091109Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression} 
    11101110\end{table} 
    … …  
    11261126% ------------------------------------------------------------------------------------------------------------- 
    11271127\section{Tracer/Dynamics trends  (\protect\ngn{namtrd})} 
    1128 \label{DIA_trd} 
     1128\label{sec:DIA_trd} 
    11291129 
    11301130%------------------------------------------namtrd---------------------------------------------------- 
    … …  
    11661166% ------------------------------------------------------------------------------------------------------------- 
    11671167\section{FLO: On-Line Floats trajectories (\protect\key{floats})} 
    1168 \label{FLO} 
     1168\label{sec:FLO} 
    11691169%--------------------------------------------namflo------------------------------------------------------- 
    11701170\forfile{../namelists/namflo}  
    … …  
    12741274% ------------------------------------------------------------------------------------------------------------- 
    12751275\section{Harmonic analysis of tidal constituents (\protect\key{diaharm}) } 
    1276 \label{DIA_diag_harm} 
     1276\label{sec:DIA_diag_harm} 
    12771277 
    12781278%------------------------------------------namdia_harm---------------------------------------------------- 
    … …  
    13161316% ------------------------------------------------------------------------------------------------------------- 
    13171317\section{Transports across sections (\protect\key{diadct}) } 
    1318 \label{DIA_diag_dct} 
     1318\label{sec:DIA_diag_dct} 
    13191319 
    13201320%------------------------------------------namdct---------------------------------------------------- 
    … …  
    14671467% ================================================================ 
    14681468\section{Diagnosing the steric effect in sea surface height} 
    1469 \label{DIA_steric} 
     1469\label{sec:DIA_steric} 
    14701470 
    14711471 
    … …  
    15001500 
    15011501A non-Boussinesq fluid conserves mass. It satisfies the following relations: 
    1502 \begin{equation} \label{Eq_MV_nBq}  
     1502\begin{equation} \label{eq:MV_nBq}  
    15031503\begin{split}  
    15041504\mathcal{M} &=  \mathcal{V}  \;\bar{\rho}       \\ 
    … …  
    15071507\end{equation} 
    15081508Temporal changes in total mass is obtained from the density conservation equation : 
    1509 \begin{equation}  \label{Eq_Co_nBq} 
     1509\begin{equation}  \label{eq:Co_nBq} 
    15101510\frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} ) = \left. \frac{\textit{emp}}{e_3}\right|_\textit{surface} 
    15111511\end{equation} 
    … …  
    15131513exchanges with the other media of the Earth system (atmosphere, sea-ice, land).  
    15141514Its global averaged leads to the total mass change  
    1515 \begin{equation}  \label{Eq_Mass_nBq} 
     1515\begin{equation}  \label{eq:Mass_nBq} 
    15161516\partial_t \mathcal{M} = \mathcal{A} \;\overline{\textit{emp}} 
    15171517\end{equation} 
    15181518where $\overline{\textit{emp}}=\int_S \textit{emp}\,ds$ is the net mass flux  
    15191519through the ocean surface. 
    1520 Bringing \eqref{Eq_Mass_nBq} and the time derivative of \eqref{Eq_MV_nBq}  
     1520Bringing \autoref{eq:Mass_nBq} and the time derivative of \autoref{eq:MV_nBq}  
    15211521together leads to the evolution equation of the mean sea level 
    1522 \begin{equation} \label{Eq_ssh_nBq} 
     1522\begin{equation} \label{eq:ssh_nBq} 
    15231523  \partial_t \bar{\eta} =  \frac{\overline{\textit{emp}}}{ \bar{\rho}}  
    15241524               - \frac{\mathcal{V}}{\mathcal{A}}  \;\frac{\partial_t \bar{\rho} }{\bar{\rho}} 
    15251525\end{equation} 
    1526 The first term in equation \eqref{Eq_ssh_nBq} alters sea level by adding or  
     1526The first term in equation \autoref{eq:ssh_nBq} alters sea level by adding or  
    15271527subtracting mass from the ocean.  
    15281528The second term arises from temporal changes in the global mean  
    … …  
    15311531In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$  
    15321532appears multiplied by the gravity ($i.e.$ in the hydrostatic balance of the primitive Equations).  
    1533 In particular, the mass conservation equation, \eqref{Eq_Co_nBq}, degenerates into  
     1533In particular, the mass conservation equation, \autoref{eq:Co_nBq}, degenerates into  
    15341534the incompressibility equation: 
    1535 \begin{equation}  \label{Eq_Co_Bq} 
     1535\begin{equation}  \label{eq:Co_Bq} 
    15361536\frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} ) =  \left. \frac{\textit{emp}}{\rho_o \,e_3}\right|_ \textit{surface} 
    15371537\end{equation} 
    15381538and the global average of this equation now gives the temporal change of the total volume, 
    1539 \begin{equation}  \label{Eq_V_Bq} 
     1539\begin{equation}  \label{eq:V_Bq} 
    15401540  \partial_t \mathcal{V} =   \mathcal{A} \;\frac{\overline{\textit{emp}}}{\rho_o}  
    15411541\end{equation} 
    … …  
    15531553by the Boussinesq model, via the steric contribution to the sea level, $\eta_s$, a spatially  
    15541554uniform variable, as follows: 
    1555 \begin{equation}  \label{Eq_M_Bq} 
     1555\begin{equation}  \label{eq:M_Bq} 
    15561556   \mathcal{M}_o  =  \mathcal{M} + \rho_o \,\eta_s \,\mathcal{A}  
    15571557\end{equation} 
    … …  
    15591559the ocean surface is converted into a mean change in sea level. Introducing the total density  
    15601560anomaly, $\mathcal{D}= \int_D d_a \,dv$, where $d_a= (\rho -\rho_o ) / \rho_o$   
    1561 is the density anomaly used in \NEMO (cf. \S\ref{TRA_eos}) in \eqref{Eq_M_Bq} 
     1561is the density anomaly used in \NEMO (cf. \autoref{subsec:TRA_eos}) in \autoref{eq:M_Bq} 
    15621562leads to a very simple form for the steric height: 
    1563 \begin{equation}  \label{Eq_steric_Bq} 
     1563\begin{equation}  \label{eq:steric_Bq} 
    15641564   \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D}  
    15651565\end{equation} 
    … …  
    15811581(wetting and drying of grid point is not allowed).  
    15821582   
    1583 Third, the discretisation of \eqref{Eq_steric_Bq} depends on the type of free surface 
     1583Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface 
    15841584which is considered. In the non linear free surface case, $i.e.$ \key{vvl} defined, it is 
    15851585given by 
    1586 \begin{equation}  \label{Eq_discrete_steric_Bq} 
     1586\begin{equation}  \label{eq:discrete_steric_Bq} 
    15871587   \eta_s =  - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t} e_{2t} e_{3t} } 
    15881588                  { \sum_{i,\,j,\,k}         e_{1t} e_{2t} e_{3t} }  
    … …  
    15901590whereas in the linear free surface, the volume above the \textit{z=0} surface must be explicitly taken  
    15911591into account to better approximate the total ocean mass and thus the steric sea level: 
    1592 \begin{equation}  \label{Eq_discrete_steric_Bq} 
     1592\begin{equation}  \label{eq:discrete_steric_Bq} 
    15931593   \eta_s =  - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} d_a\; e_{1t}e_{2t} \eta } 
    15941594                     {\sum_{i,\,j,\,k} e_{1t}e_{2t}e_{3t} + \sum_{i,\,j}           e_{1t}e_{2t} \eta }  
    … …  
    16081608In AR5 outputs, the thermosteric sea level is demanded. It is steric sea level due to  
    16091609changes in ocean density arising just from changes in temperature. It is given by: 
    1610 \begin{equation}  \label{Eq_thermosteric_Bq} 
     1610\begin{equation}  \label{eq:thermosteric_Bq} 
    16111611   \eta_s = - \frac{1}{\mathcal{A}} \int_D d_a(T,S_o,p_o) \,dv 
    16121612\end{equation} 
    … …  
    16221622% ------------------------------------------------------------------------------------------------------------- 
    16231623\section{Other diagnostics (\protect\key{diahth}, \protect\key{diaar5})} 
    1624 \label{DIA_diag_others} 
     1624\label{sec:DIA_diag_others} 
    16251625 
    16261626 
    … …  
    16581658as well as for the World Ocean. The sub-basin decomposition requires an input file  
    16591659(\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask  
    1660 been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}).  
     1660been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}).  
    16611661 
    16621662%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    16631663\begin{figure}[!t]     \begin{center} 
    16641664\includegraphics[width=1.0\textwidth]{Fig_mask_subasins} 
    1665 \caption{   \protect\label{Fig_mask_subasins} 
     1665\caption{   \protect\label{fig:mask_subasins} 
    16661666Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute 
    16671667the heat and salt transports as well as the meridional stream-function: Atlantic basin (red),  
    … …  
    16811681A series of diagnostics has been added in the \mdl{diaar5}.  
    16821682They corresponds to outputs that are required for AR5 simulations (CMIP5) 
    1683 (see also Section \ref{DIA_steric} for one of them).  
     1683(see also  \autoref{sec:DIA_steric} for one of them).  
    16841684Activating those outputs requires to define the \key{diaar5} CPP key. 
    16851685 

  • branches/2017/dev_merge_2017/DOC/tex_sub/chap_DIU.tex

    r9394 r9407  
    66% ================================================================ 
    77\chapter{Diurnal SST Models (DIU)} 
    8 \label{DIU} 
     8\label{chap:DIU} 
    99 
    1010\minitoc 
    … …  
    5454%=============================================================== 
    5555\section{Warm layer model} 
    56 \label{warm_layer_sec} 
     56\label{sec:warm_layer_sec} 
    5757%=============================================================== 
    5858 
    … …  
    6262\frac{\partial{\Delta T_{\rm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p 
    6363\nu}-\frac{(\nu+1)ku^*_{w}f(L_a)\Delta T}{D_T\Phi\!\left(\frac{D_T}{L}\right)} \mbox{,} 
    64 \label{ecmwf1} \\ 
    65 L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{ecmwf2} 
     64\label{eq:ecmwf1} \\ 
     65L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:ecmwf2} 
    6666\end{eqnarray} 
    6767where $\Delta T_{\rm{wl}}$ is the temperature difference between the top of the warm 
    6868layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal. In 
    69 equation (\ref{ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion 
     69equation (\autoref{eq:ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion 
    7070coefficient of water, $\kappa=0.4$ is von K\'{a}rm\'{a}n's constant, $c_p$ is the heat 
    7171capacity at constant pressure of sea water, $\rho_w$ is the 
    … …  
    8181$u^*_{w} = u_{10}\sqrt{\frac{C_d\rho_a}{\rho_w}}$, where $C_d$ is 
    8282the drag coefficient, and $\rho_a$ is the density of air.  The symbol $Q$ in equation 
    83 (\ref{ecmwf1}) is the instantaneous total thermal energy 
     83(\autoref{eq:ecmwf1}) is the instantaneous total thermal energy 
    8484flux into 
    8585the diurnal layer, $i.e.$ 
    8686\begin{equation} 
    87 Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,} \label{e_flux_eqn} 
     87Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,} \label{eq:e_flux_eqn} 
    8888\end{equation} 
    8989where $Q_{\rm{h}}$ is the sensible and latent heat flux, $Q_{\rm{lw}}$ is the long 
    9090wave flux, and $Q_{\rm{sol}}$ is the solar flux absorbed 
    9191within the diurnal warm layer. For $Q_{\rm{sol}}$ the 9 term 
    92 representation of \citet{Gentemann_al_JGR09} is used.  In equation \ref{ecmwf1} 
     92representation of \citet{Gentemann_al_JGR09} is used.  In equation \autoref{eq:ecmwf1} 
    9393the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$, where $L_a=0.3$\footnote{This 
    9494is a global average value, more accurately $L_a$ could be computed as 
    … …  
    1031034\zeta^2}{1+3\zeta+0.25\zeta^2} &(\zeta \ge 0) \\ 
    104104                                    (1 - 16\zeta)^{-\frac{1}{2}} & (\zeta < 0) \mbox{,} 
    105                                     \end{array} \right. \label{stab_func_eqn} 
     105                                    \end{array} \right. \label{eq:stab_func_eqn} 
    106106\end{equation} 
    107107where $\zeta=\frac{D_T}{L}$.  It is clear that the first derivative of 
    108 (\ref{stab_func_eqn}), and thus of (\ref{ecmwf1}), 
    109 is discontinuous at $\zeta=0$ ($i.e.$ $Q\rightarrow0$ in equation (\ref{ecmwf2})). 
     108(\autoref{eq:stab_func_eqn}), and thus of (\autoref{eq:ecmwf1}), 
     109is discontinuous at $\zeta=0$ ($i.e.$ $Q\rightarrow0$ in equation (\autoref{eq:ecmwf2})). 
    110110 
    111 The two terms on the right hand side of (\ref{ecmwf1}) represent different processes. 
     111The two terms on the right hand side of (\autoref{eq:ecmwf1}) represent different processes. 
    112112The first term is simply the diabatic heating or cooling of the 
    113113diurnal warm 
    … …  
    121121 
    122122\section{Cool skin model} 
    123 \label{cool_skin_sec} 
     123\label{sec:cool_skin_sec} 
    124124 
    125125%=============================================================== 
    … …  
    131131Saunders equation for the cool skin temperature difference $\Delta T_{\rm{cs}}$ becomes 
    132132\begin{equation} 
    133 \label{sunders_eqn} 
     133\label{eq:sunders_eqn} 
    134134\Delta T_{\rm{cs}}=\frac{Q_{\rm{ns}}\delta}{k_t} \mbox{,} 
    135135\end{equation} 
    … …  
    138138skin layer and is given by 
    139139\begin{equation} 
    140 \label{sunders_thick_eqn} 
     140\label{eq:sunders_thick_eqn} 
    141141\delta=\frac{\lambda \mu}{u^*_{w}} \mbox{,} 
    142142\end{equation} 
    … …  
    144144proportionality which \citet{Saunders_JAS82} suggested varied between 5 and 10. 
    145145 
    146 The value of $\lambda$ used in equation (\ref{sunders_thick_eqn}) is that of 
     146The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of 
    147147\citet{Artale_al_JGR02}, 
    148148which is shown in \citet{Tu_Tsuang_GRL05} to outperform a number of other 
    149149parametrisations at both low and high wind speeds. Specifically, 
    150150\begin{equation} 
    151 \label{artale_lambda_eqn} 
     151\label{eq:artale_lambda_eqn} 
    152152\lambda = \frac{ 8.64\times10^4 u^*_{w} k_t }{ \rho c_p h \mu \gamma }\mbox{,} 
    153153\end{equation} 
    … …  
    155155$\gamma$ is a dimensionless function of wind speed $u$: 
    156156\begin{equation} 
    157 \label{artale_gamma_eqn} 
     157\label{eq:artale_gamma_eqn} 
    158158\gamma = \left\{ \begin{matrix} 
    159159                     0.2u+0.5\mbox{,} & u \le 7.5\,\mbox{ms}^{-1} \\ 

  • branches/2017/dev_merge_2017/DOC/tex_sub/chap_DOM.tex

    r9393 r9407  
    55% ================================================================ 
    66\chapter{Space Domain (DOM)} 
    7 \label{DOM} 
     7\label{chap:DOM} 
    88\minitoc 
    99 
    … …  
    2020$\ $\newline    % force a new line 
    2121 
    22 Having defined the continuous equations in Chap.~\ref{PE} and chosen a time  
    23 discretization Chap.~\ref{STP}, we need to choose a discretization on a grid,  
     22Having defined the continuous equations in \autoref{chap:PE} and chosen a time  
     23discretization \autoref{chap:STP}, we need to choose a discretization on a grid,  
    2424and numerical algorithms. In the present chapter, we provide a general description  
    2525of the staggered grid used in \NEMO, and other information relevant to the main  
    … …  
    3232% ================================================================ 
    3333\section{Fundamentals of the discretisation} 
    34 \label{DOM_basics} 
     34\label{sec:DOM_basics} 
    3535 
    3636% ------------------------------------------------------------------------------------------------------------- 
    … …  
    3838% ------------------------------------------------------------------------------------------------------------- 
    3939\subsection{Arrangement of variables} 
    40 \label{DOM_cell} 
     40\label{subsec:DOM_cell} 
    4141 
    4242%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    4343\begin{figure}[!tb]    \begin{center} 
    4444\includegraphics[width=0.90\textwidth]{Fig_cell} 
    45 \caption{ \protect\label{Fig_cell}     
     45\caption{ \protect\label{fig:cell}     
    4646Arrangement of variables. $t$ indicates scalar points where temperature,  
    4747salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$)  
    … …  
    5656space directions. The arrangement of variables is the same in all directions.  
    5757It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector  
    58 points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}).  
     58points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:cell}).  
    5959This is the generalisation to three dimensions of the well-known ``C'' grid in  
    6060Arakawa's classification \citep{Mesinger_Arakawa_Bk76}. The relative and  
    … …  
    6666by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$.  
    6767The grid-points are located at integer or integer and a half value of $(i,j,k)$ as  
    68 indicated on Table \ref{Tab_cell}. In all the following, subscripts $u$, $v$, $w$,  
     68indicated on \autoref{tab:cell}. In all the following, subscripts $u$, $v$, $w$,  
    6969$f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale  
    7070factors are defined. Each scale factor is defined as the local analytical value  
    71 provided by \eqref{Eq_scale_factors}. As a result, the mesh on which partial  
     71provided by \autoref{eq:scale_factors}. As a result, the mesh on which partial  
    7272derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, and  
    7373$\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity.  
    … …  
    7878from their analytical expression. This preserves the symmetry of the discrete set  
    7979of equations and therefore satisfies many of the continuous properties (see  
    80 Appendix~\ref{Apdx_C}). A similar, related remark can be made about the domain  
     80\autoref{apdx:C}). A similar, related remark can be made about the domain  
    8181size: when needed, an area, volume, or the total ocean depth must be evaluated  
    82 as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section).  
     82as the sum of the relevant scale factors (see \autoref{eq:DOM_bar}) in the next section).  
    8383 
    8484%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    9595fw & $i+1/2$   & $j+1/2$   & $k+1/2$   \\ \hline 
    9696\end{tabular} 
    97 \caption{ \protect\label{Tab_cell} 
     97\caption{ \protect\label{tab:cell} 
    9898Location of grid-points as a function of integer or integer and a half value of the column,  
    9999line or level. This indexing is only used for the writing of the semi-discrete equation.  
    100100In the code, the indexing uses integer values only and has a reverse direction  
    101 in the vertical (see \S\ref{DOM_Num_Index})} 
     101in the vertical (see \autoref{subsec:DOM_Num_Index})} 
    102102\end{center} 
    103103\end{table} 
    … …  
    108108% ------------------------------------------------------------------------------------------------------------- 
    109109\subsection{Discrete operators} 
    110 \label{DOM_operators} 
     110\label{subsec:DOM_operators} 
    111111 
    112112Given the values of a variable $q$ at adjacent points, the differencing and  
    113113averaging operators at the midpoint between them are: 
    114 \begin{subequations} \label{Eq_di_mi} 
     114\begin{subequations} \label{eq:di_mi} 
    115115\begin{align} 
    116116 \delta _i [q]       &=  \  \    q(i+1/2)  - q(i-1/2)    \\ 
    … …  
    120120 
    121121Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and  
    122 $k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a  
     122$k+1/2$. Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a  
    123123variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$-  
    124124and $w$-points while its Laplacien is defined at $t$-point. These operators have  
    125125the following discrete forms in the curvilinear $s$-coordinate system: 
    126 \begin{equation} \label{Eq_DOM_grad} 
     126\begin{equation} \label{eq:DOM_grad} 
    127127\nabla q\equiv    \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,\mathbf{i} 
    128128      +  \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,\mathbf{j} 
    129129      +  \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,\mathbf{k} 
    130130\end{equation} 
    131 \begin{multline} \label{Eq_DOM_lap} 
     131\begin{multline} \label{eq:DOM_lap} 
    132132\Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
    133133       \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right] 
    … …  
    136136\end{multline} 
    137137 
    138 Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$  
     138Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$  
    139139defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$,  
    140140and $f$-points, and its divergence defined at $t$-points: 
    141 \begin{eqnarray}  \label{Eq_DOM_curl} 
     141\begin{eqnarray}  \label{eq:DOM_curl} 
    142142 \nabla \times {\rm{\bf A}}\equiv & 
    143143      \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \mathbf{i} \\  
    … …  
    145145 +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \mathbf{k} 
    146146 \end{eqnarray} 
    147 \begin{eqnarray} \label{Eq_DOM_div} 
     147\begin{eqnarray} \label{eq:DOM_div} 
    148148\nabla \cdot \rm{\bf A} \equiv  
    149149    \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] 
    … …  
    153153The vertical average over the whole water column denoted by an overbar becomes  
    154154for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area): 
    155 \begin{equation} \label{DOM_bar} 
     155\begin{equation} \label{eq:DOM_bar} 
    156156\bar q   =         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk}  
    157157      \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } 
    … …  
    163163 
    164164In continuous form, the following properties are satisfied: 
    165 \begin{equation} \label{Eq_DOM_curl_grad} 
     165\begin{equation} \label{eq:DOM_curl_grad} 
    166166\nabla \times \nabla q ={\rm {\bf {0}}} 
    167167\end{equation} 
    168 \begin{equation} \label{Eq_DOM_div_curl} 
     168\begin{equation} \label{eq:DOM_div_curl} 
    169169\nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0 
    170170\end{equation} 
    … …  
    181181operators, $i.e.$ 
    182182\begin{align}  
    183 \label{DOM_di_adj} 
     183\label{eq:DOM_di_adj} 
    184184\sum\limits_i { a_i \;\delta _i \left[ b \right]}  
    185185   &\equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} }      \\ 
    186 \label{DOM_mi_adj} 
     186\label{eq:DOM_mi_adj} 
    187187\sum\limits_i { a_i \;\overline b^{\,i}}  
    188188   & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} }  
    … …  
    192192$\delta_i^*=\delta_{i+1/2}$ and  
    193193${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively.  
    194 These two properties will be used extensively in the Appendix~\ref{Apdx_C} to  
     194These two properties will be used extensively in the \autoref{apdx:C} to  
    195195demonstrate integral conservative properties of the discrete formulation chosen. 
    196196 
    … …  
    199199% ------------------------------------------------------------------------------------------------------------- 
    200200\subsection{Numerical indexing} 
    201 \label{DOM_Num_Index} 
     201\label{subsec:DOM_Num_Index} 
    202202 
    203203%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    204204\begin{figure}[!tb]  \begin{center} 
    205205\includegraphics[width=0.90\textwidth]{Fig_index_hor} 
    206 \caption{   \protect\label{Fig_index_hor}     
     206\caption{   \protect\label{fig:index_hor}     
    207207Horizontal integer indexing used in the \textsc{Fortran} code. The dashed area indicates  
    208208the cell in which variables contained in arrays have the same $i$- and $j$-indices} 
    … …  
    211211 
    212212The array representation used in the \textsc{Fortran} code requires an integer  
    213 indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is  
     213indexing while the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is  
    214214associated with the use of integer values for $t$-points and both integer and  
    215215integer and a half values for all the other points. Therefore a specific integer  
    … …  
    222222% ----------------------------------- 
    223223\subsubsection{Horizontal indexing} 
    224 \label{DOM_Num_Index_hor} 
    225  
    226 The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}.  
     224\label{subsec:DOM_Num_Index_hor} 
     225 
     226The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}.  
    227227For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point  
    228 (northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}).  
     228(northward $v$-point) have the same index (see the dashed area in \autoref{fig:index_hor}).  
    229229A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. 
    230230 
    … …  
    233233% ----------------------------------- 
    234234\subsubsection{Vertical indexing} 
    235 \label{DOM_Num_Index_vertical} 
     235\label{subsec:DOM_Num_Index_vertical} 
    236236 
    237237In the vertical, the chosen indexing requires special attention since the  
    238238$k$-axis is re-orientated downward in the \textsc{Fortran} code compared  
    239 to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}.  
     239to the indexing used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}.  
    240240The sea surface corresponds to the $w$-level $k=1$ which is the same index  
    241 as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$)  
     241as $t$-level just below (\autoref{fig:index_vert}). The last $w$-level ($k=jpk$)  
    242242either corresponds to the ocean floor or is inside the bathymetry while the last  
    243 $t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that  
     243$t$-level is always inside the bathymetry (\autoref{fig:index_vert}). Note that  
    244244for an increasing $k$ index, a $w$-point and the $t$-point just below have the  
    245245same $k$ index, in opposition to what is done in the horizontal plane where  
    246246it is the $t$-point and the nearest velocity points in the direction of the horizontal  
    247247axis that have the same $i$ or $j$ index (compare the dashed area in  
    248 Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are  
     248\autoref{fig:index_hor} and \autoref{fig:index_vert}). Since the scale factors are  
    249249chosen to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran}  
    250250code \emph{before all the vertical derivatives} of the discrete equations given in  
    … …  
    254254\begin{figure}[!pt]    \begin{center} 
    255255\includegraphics[width=.90\textwidth]{Fig_index_vert} 
    256 \caption{ \protect\label{Fig_index_vert}      
     256\caption{ \protect\label{fig:index_vert}      
    257257Vertical integer indexing used in the \textsc{Fortran } code. Note that  
    258258the $k$-axis is orientated downward. The dashed area indicates the cell in  
    … …  
    265265% ----------------------------------- 
    266266\subsubsection{Domain size} 
    267 \label{DOM_size} 
     267\label{subsec:DOM_size} 
    268268 
    269269The total size of the computational domain is set by the parameters \np{jpiglo},  
    … …  
    273273%%% 
    274274Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when the code is  
    275 run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \S\ref{LBC_mpp}). 
     275run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). 
    276276 
    277277 
    … …  
    282282% ================================================================ 
    283283\section{Needed fields} 
    284 \label{DOM_fields} 
     284\label{sec:DOM_fields} 
    285285The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined  
    286286by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.  
    287287The grid-points are located at integer or integer and a half values of as indicated  
    288 in Table~\ref{Tab_cell}. The associated scale factors are defined using the   
    289 analytical first derivative of the transformation \eqref{Eq_scale_factors}.  
     288in \autoref{tab:cell}. The associated scale factors are defined using the   
     289analytical first derivative of the transformation \autoref{eq:scale_factors}.  
    290290Necessary fields for configuration definition are: \\ 
    291291Geographic position : 
    … …  
    316316% ------------------------------------------------------------------------------------------------------------- 
    317317%\subsection{List of needed fields to build DOMAIN} 
    318 %\label{DOM_fields_list} 
     318%\label{subsec:DOM_fields_list} 
    319319 
    320320 
    … …  
    323323% ================================================================ 
    324324\section{Horizontal grid mesh (\protect\mdl{domhgr})} 
    325 \label{DOM_hgr} 
     325\label{sec:DOM_hgr} 
    326326 
    327327% ------------------------------------------------------------------------------------------------------------- 
    … …  
    329329% ------------------------------------------------------------------------------------------------------------- 
    330330\subsection{Coordinates and scale factors} 
    331 \label{DOM_hgr_coord_e} 
     331\label{subsec:DOM_hgr_coord_e} 
    332332 
    333333The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined  
    334334by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.  
    335335The grid-points are located at integer or integer and a half values of as indicated  
    336 in Table~\ref{Tab_cell}. The associated scale factors are defined using the  
    337 analytical first derivative of the transformation \eqref{Eq_scale_factors}. These  
     336in \autoref{tab:cell}. The associated scale factors are defined using the  
     337analytical first derivative of the transformation \autoref{eq:scale_factors}. These  
    338338definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which  
    339339provide the horizontal and vertical meshes, respectively. This section deals with  
    … …  
    343343analytical expressions of the longitude $\lambda$ and  latitude $\varphi$ as a  
    344344function of  $(i,j)$. The horizontal scale factors are calculated using  
    345 \eqref{Eq_scale_factors}. For example, when the longitude and latitude are  
     345\autoref{eq:scale_factors}. For example, when the longitude and latitude are  
    346346function of a single value ($i$ and $j$, respectively) (geographical configuration  
    347347of the mesh), the horizontal mesh definition reduces to define the wanted  
    … …  
    382382allowing the user to set arbitrary jumps in thickness between adjacent layers)  
    383383\citep{Treguier1996}. An example of the effect of such a choice is shown in  
    384 Fig.~\ref{Fig_zgr_e3}. 
     384\autoref{fig:zgr_e3}. 
    385385%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    386386\begin{figure}[!t]     \begin{center} 
    387387\includegraphics[width=0.90\textwidth]{Fig_zgr_e3} 
    388 \caption{ \protect\label{Fig_zgr_e3}     
     388\caption{ \protect\label{fig:zgr_e3}     
    389389Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,  
    390390and (b) analytically derived grid-point position and scale factors.  
    … …  
    401401% ------------------------------------------------------------------------------------------------------------- 
    402402\subsection{Choice of horizontal grid} 
    403 \label{DOM_hgr_msh_choice} 
     403\label{subsec:DOM_hgr_msh_choice} 
    404404 
    405405 
    … …  
    408408% ------------------------------------------------------------------------------------------------------------- 
    409409\subsection{Output grid files} 
    410 \label{DOM_hgr_files} 
     410\label{subsec:DOM_hgr_files} 
    411411 
    412412All the arrays relating to a particular ocean model configuration (grid-point  
    … …  
    426426% ================================================================ 
    427427\section{Vertical grid (\protect\mdl{domzgr})} 
    428 \label{DOM_zgr} 
     428\label{sec:DOM_zgr} 
    429429%-----------------------------------------nam_zgr & namdom------------------------------------------- 
    430430%\forfile{../namelists/namzgr}  
    … …  
    444444\begin{figure}[!tb]    \begin{center} 
    445445\includegraphics[width=1.0\textwidth]{Fig_z_zps_s_sps} 
    446 \caption{  \protect\label{Fig_z_zps_s_sps}    
     446\caption{  \protect\label{fig:z_zps_s_sps}    
    447447The ocean bottom as seen by the model:  
    448448(a) $z$-coordinate with full step,  
    … …  
    451451(d) hybrid $s-z$ coordinate,  
    452452(e) hybrid $s-z$ coordinate with partial step, and  
    453 (f) same as (e) but in the non-linear free surface (\np{ln\_linssh}\forcode{ = .false.}).  
     453(f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}).  
    454454Note that the non-linear free surface can be used with any of the  
    4554555 coordinates (a) to (e).} 
    … …  
    460460must be done once of all at the beginning of an experiment. It is not intended as an  
    461461option which can be enabled or disabled in the middle of an experiment. Three main  
    462 choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step  
     462choices are offered (\autoref{fig:z_zps_s_sps}a to c): $z$-coordinate with full step  
    463463bathymetry (\np{ln\_zco}\forcode{ = .true.}), $z$-coordinate with partial step bathymetry  
    464464(\np{ln\_zps}\forcode{ = .true.}), or generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}).  
    465465Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate  
    466 (Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). By default a non-linear free surface is used: 
     466(\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e). By default a non-linear free surface is used: 
    467467the coordinate follow the time-variation of the free surface so that the transformation is time dependent:  
    468 $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln\_linssh}\forcode{ = .true.}),  
     468$z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}f). When a linear free surface is assumed (\np{ln\_linssh}\forcode{ = .true.}),  
    469469the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface  
    470470(in other words, the top of the ocean in not a rigid-lid).  
    … …  
    513513% ------------------------------------------------------------------------------------------------------------- 
    514514\subsection{Meter bathymetry} 
    515 \label{DOM_bathy} 
     515\label{subsec:DOM_bathy} 
    516516 
    517517Three options are possible for defining the bathymetry, according to the  
    … …  
    541541This is unnecessary when the ocean is forced by fixed atmospheric conditions,  
    542542so these seas can be removed from the ocean domain. The user has the option  
    543 to set the bathymetry in closed seas to zero (see \S\ref{MISC_closea}), but the  
     543to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), but the  
    544544code has to be adapted to the user's configuration.  
    545545 
    … …  
    549549\subsection[$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and ref. coordinate] 
    550550            {$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and reference coordinate} 
    551 \label{DOM_zco} 
     551\label{subsec:DOM_zco} 
    552552 
    553553%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    554554\begin{figure}[!tb]    \begin{center} 
    555555\includegraphics[width=0.90\textwidth]{Fig_zgr} 
    556 \caption{ \protect\label{Fig_zgr}     
     556\caption{ \protect\label{fig:zgr}     
    557557Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level functions for  
    558558(a) T-point depth and (b) the associated scale factor as computed  
    559 from \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.} 
     559from \autoref{eq:DOM_zgr_ana} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate.} 
    560560\end{center}   \end{figure} 
    561561%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    563563The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$  
    564564and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on  
    565 Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the  
     565\autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the  
    566566ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the  
    567567additional $t$-point at $jk=jpk$ is below the sea floor and is not used.  
    … …  
    579579near the ocean surface. The following function is proposed as a standard for a  
    580580$z$-coordinate (with either full or partial steps):  
    581 \begin{equation} \label{DOM_zgr_ana} 
     581\begin{equation} \label{eq:DOM_zgr_ana} 
    582582\begin{split} 
    583583 z_0 (k)    &= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\  
    … …  
    588588expression allows us to define a nearly uniform vertical location of levels at the  
    589589ocean top and bottom with a smooth hyperbolic tangent transition in between  
    590 (Fig.~\ref{Fig_zgr}). 
     590(\autoref{fig:zgr}). 
    591591 
    592592If the ice shelf cavities are opened (\np{ln\_isfcav}\forcode{ = .true.}), the definition of $z_0$ is the same.  
    593593However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 
    594 \begin{equation} \label{DOM_zgr_ana} 
     594\begin{equation} \label{eq:DOM_zgr_ana} 
    595595\begin{split} 
    596596 e_3^T(k) &= z_W (k+1) - z_W (k)   \\ 
    … …  
    605605surface (bottom) layers and a depth which varies from 0 at the sea surface to a  
    606606minimum of $-5000~m$. This leads to the following conditions: 
    607 \begin{equation} \label{DOM_zgr_coef} 
     607\begin{equation} \label{eq:DOM_zgr_coef} 
    608608\begin{split} 
    609609 e_3 (1+1/2)      &=10. \\  
    … …  
    616616With the choice of the stretching $h_{cr} =3$ and the number of levels  
    617617\jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in  
    618 \eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is  
     618\autoref{eq:DOM_zgr_ana} have been determined such that \autoref{eq:DOM_zgr_coef} is  
    619619satisfied, through an optimisation procedure using a bisection method. For the first  
    620620standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$,  
    621621$h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and  
    622 scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and  
    623 given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters  
     622scale factors as a function of the model levels are shown in \autoref{fig:zgr} and  
     623given in \autoref{tab:orca_zgr}. Those values correspond to the parameters  
    624624\np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist.  
    625625 
    … …  
    67567531 &  \textbf{5250.23}& 5000.00 &   \textbf{500.56} & 500.33 \\ \hline 
    676676\end{tabular} \end{center}  
    677 \caption{ \protect\label{Tab_orca_zgr}    
     677\caption{ \protect\label{tab:orca_zgr}    
    678678Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed  
    679 from \eqref{DOM_zgr_ana} using the coefficients given in \eqref{DOM_zgr_coef}} 
     679from \autoref{eq:DOM_zgr_ana} using the coefficients given in \autoref{eq:DOM_zgr_coef}} 
    680680\end{table} 
    681681%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    685685% ------------------------------------------------------------------------------------------------------------- 
    686686\subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}\forcode{ = .true.})} 
    687 \label{DOM_zps} 
     687\label{subsec:DOM_zps} 
    688688%--------------------------------------------namdom------------------------------------------------------- 
    689689\forfile{../namelists/namdom}  
    … …  
    717717% ------------------------------------------------------------------------------------------------------------- 
    718718\subsection{$S$-coordinate (\protect\np{ln\_sco}\forcode{ = .true.})} 
    719 \label{DOM_sco} 
     719\label{subsec:DOM_sco} 
    720720%------------------------------------------nam_zgr_sco--------------------------------------------------- 
    721721%\forfile{../namelists/namzgr_sco}  
    … …  
    726726function or its derivative, respectively: 
    727727 
    728 \begin{equation} \label{DOM_sco_ana} 
     728\begin{equation} \label{eq:DOM_sco_ana} 
    729729\begin{split} 
    730730 z(k)       &= h(i,j) \; z_0(k)  \\ 
    … …  
    737737surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean  
    738738depth, since a mixed step-like and bottom-following representation of the  
    739 topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f). 
     739topography can be used (\autoref{fig:z_zps_s_sps}d-e) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}f). 
    740740The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects  
    741741the sea bed and becomes a pseudo z-coordinate.  
    … …  
    764764\end{equation} 
    765765 
    766 \begin{equation} \label{DOM_sco_function} 
     766\begin{equation} \label{eq:DOM_sco_function} 
    767767\begin{split} 
    768768C(s)  &=  \frac{ \left[   \tanh{ \left( \theta \, (s+b) \right)}  
    … …  
    784784\begin{figure}[!ht]    \begin{center} 
    785785\includegraphics[width=1.0\textwidth]{Fig_sco_function} 
    786 \caption{  \protect\label{Fig_sco_function}    
     786\caption{  \protect\label{fig:sco_function}    
    787787Examples of the stretching function applied to a seamount; from left to right:  
    788788surface, surface and bottom, and bottom intensified resolutions} 
    … …  
    794794are the surface and bottom control parameters such that $0\leqslant \theta \leqslant 20$, and  
    795795$0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom  
    796 increase of the vertical resolution (Fig.~\ref{Fig_sco_function}). 
     796increase of the vertical resolution (\autoref{fig:sco_function}). 
    797797 
    798798Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows  
    … …  
    807807The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 
    808808 
    809 \begin{equation} \label{DOM_gamma_deriv} 
     809\begin{equation} \label{eq:DOM_gamma_deriv} 
    810810\gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right) 
    811811\end{equation} 
    812812 
    813813Where: 
    814 \begin{equation} \label{DOM_gamma} 
     814\begin{equation} \label{eq:DOM_gamma} 
    815815f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} \quad \text{ and } \quad \sigma = \frac{k}{n-1}  
    816816\end{equation} 
    … …  
    821821and bottom depths. The bottom cell depth in this example is given as a function of water depth: 
    822822 
    823 \begin{equation} \label{DOM_zb} 
     823\begin{equation} \label{eq:DOM_zb} 
    824824Z_b= h a + b 
    825825\end{equation} 
    … …  
    831831   \includegraphics[width=1.0\textwidth]{FIG_DOM_compare_coordinates_surface} 
    832832        \caption{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 100m for a idealised bathymetry that goes from 50m to 5500m depth. For clarity every third coordinate surface is shown.} 
    833     \label{fig_compare_coordinates_surface} 
     833    \label{fig:fig_compare_coordinates_surface} 
    834834\end{figure} 
    835835%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    845845% ------------------------------------------------------------------------------------------------------------- 
    846846\subsection{$Z^*$- or $S^*$-coordinate (\protect\np{ln\_linssh}\forcode{ = .false.}) } 
    847 \label{DOM_zgr_star} 
     847\label{subsec:DOM_zgr_star} 
    848848 
    849849This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site.  
    … …  
    855855% ------------------------------------------------------------------------------------------------------------- 
    856856\subsection{Level bathymetry and mask} 
    857 \label{DOM_msk} 
     857\label{subsec:DOM_msk} 
    858858 
    859859Whatever the vertical coordinate used, the model offers the possibility of  
    … …  
    892892 
    893893Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with  
    894 the numerical indexing used (\S~\ref{DOM_Num_Index}). Nevertheless, $wmask$ are required  
     894the numerical indexing used (\autoref{subsec:DOM_Num_Index}). Nevertheless, $wmask$ are required  
    895895with oceean cavities to deal with the top boundary (ice shelf/ocean interface)  
    896896exactly in the same way as for the bottom boundary.  
    … …  
    900900case of an east-west cyclical boundary condition, \textit{mbathy} has its last  
    901901column equal to the second one and its first column equal to the last but one  
    902 (and so too the mask arrays) (see \S~\ref{LBC_jperio}). 
     902(and so too the mask arrays) (see \autoref{fig:LBC_jperio}). 
    903903 
    904904 
    … …  
    907907% ================================================================ 
    908908\section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
    909 \label{DTA_tsd} 
     909\label{sec:DTA_tsd} 
    910910%-----------------------------------------namtsd------------------------------------------- 
    911911\forfile{../namelists/namtsd}  
    … …  
    918918\item[\np{ln\_tsd\_init}\forcode{ = .true.}] use a T and S input files that can be given on the model grid itself or  
    919919on their native input data grid. In the latter case, the data will be interpolated on-the-fly both in the  
    920 horizontal and the vertical to the model grid (see \S~\ref{SBC_iof}). The information relative to the  
     920horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}). The information relative to the  
    921921input files are given in the \np{sn\_tem} and \np{sn\_sal} structures.  
    922922The computation is done in the \mdl{dtatsd} module. 

  • branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex

    r9394 r9407  
    55% ================================================================ 
    66\chapter{Ocean Dynamics (DYN)} 
    7 \label{DYN} 
     7\label{chap:DYN} 
    88\minitoc 
    99 
    … …  
    1111$\ $\newline      %force an empty line 
    1212 
    13 Using the representation described in Chapter \ref{DOM}, several semi-discrete  
     13Using the representation described in \autoref{chap:DOM}, several semi-discrete  
    1414space forms of the dynamical equations are available depending on the vertical  
    1515coordinate used and on the conservation properties of the vorticity term. In all  
    … …  
    3636inputs (surface wind stress calculation using bulk formulae, estimation of mixing  
    3737coefficients) that are carried out in modules SBC, LDF and ZDF and are described  
    38 in Chapters \ref{SBC}, \ref{LDF} and \ref{ZDF}, respectively.  
     38in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.  
    3939 
    4040In the present chapter we also describe the diagnostic equations used to compute  
    … …  
    5151The user has the option of extracting and outputting each tendency term from the 
    52523D momentum equations (\key{trddyn} defined), as described in  
    53 Chap.\ref{MISC}.  Furthermore, the tendency terms associated with the 2D  
     53\autoref{chap:MISC}.  Furthermore, the tendency terms associated with the 2D  
    5454barotropic vorticity balance (when \key{trdvor} is defined) can be derived from the  
    55553D terms. 
    … …  
    6464% ================================================================ 
    6565\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} 
    66 \label{DYN_divcur_wzv} 
     66\label{sec:DYN_divcur_wzv} 
    6767 
    6868%-------------------------------------------------------------------------------------------------------------- 
    … …  
    7070%-------------------------------------------------------------------------------------------------------------- 
    7171\subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
    72 \label{DYN_divcur} 
     72\label{subsec:DYN_divcur} 
    7373 
    7474The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows: 
    75 \begin{equation} \label{Eq_divcur_cur} 
     75\begin{equation} \label{eq:divcur_cur} 
    7676\zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta _{i+1/2} \left[ {e_{2v}\;v} \right] 
    7777                          -\delta _{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) 
    … …  
    7979 
    8080The horizontal divergence is defined at a $T$-point. It is given by: 
    81 \begin{equation} \label{Eq_divcur_div} 
     81\begin{equation} \label{eq:divcur_div} 
    8282\chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
    8383      \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u} \right] 
    … …  
    102102%-------------------------------------------------------------------------------------------------------------- 
    103103\subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
    104 \label{DYN_sshwzv} 
     104\label{subsec:DYN_sshwzv} 
    105105 
    106106The sea surface height is given by : 
    107 \begin{equation} \label{Eq_dynspg_ssh} 
     107\begin{equation} \label{eq:dynspg_ssh} 
    108108\begin{aligned} 
    109109\frac{\partial \eta }{\partial t} 
    … …  
    117117expressed in Kg/m$^2$/s (which is equal to mm/s), and $\rho _w$=1,035~Kg/m$^3$  
    118118is the reference density of sea water (Boussinesq approximation). If river runoff is  
    119 expressed as a surface freshwater flux (see \S\ref{SBC}) then \textit{emp} can be  
     119expressed as a surface freshwater flux (see \autoref{chap:SBC}) then \textit{emp} can be  
    120120written as the evaporation minus precipitation, minus the river runoff.  
    121121The sea-surface height is evaluated using exactly the same time stepping scheme  
    122 as the tracer equation \eqref{Eq_tra_nxt}:  
     122as the tracer equation \autoref{eq:tra_nxt}:  
    123123a leapfrog scheme in combination with an Asselin time filter, $i.e.$ the velocity appearing  
    124 in \eqref{Eq_dynspg_ssh} is centred in time (\textit{now} velocity).  
     124in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity).  
    125125This is of paramount importance. Replacing $T$ by the number $1$ in the tracer equation and summing 
    126126over the water column must lead to the sea surface height equation otherwise tracer content 
    … …  
    129129The vertical velocity is computed by an upward integration of the horizontal  
    130130divergence starting at the bottom, taking into account the change of the thickness of the levels : 
    131 \begin{equation} \label{Eq_wzv} 
     131\begin{equation} \label{eq:wzv} 
    132132\left\{   \begin{aligned} 
    133133&\left. w \right|_{k_b-1/2} \quad= 0    \qquad \text{where } k_b \text{ is the level just above the sea floor }   \\ 
    … …  
    141141of the level thicknesses, re-orientated downward. 
    142142\gmcomment{not sure of this...  to be modified with the change in emp setting} 
    143 In the case of a linear free surface, the time derivative in \eqref{Eq_wzv} disappears. 
     143In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears. 
    144144The upper boundary condition applies at a fixed level $z=0$. The top vertical velocity  
    145145is thus equal to the divergence of the barotropic transport ($i.e.$ the first term in the 
    146 right-hand-side of \eqref{Eq_dynspg_ssh}). 
     146right-hand-side of \autoref{eq:dynspg_ssh}). 
    147147 
    148148Note also that whereas the vertical velocity has the same discrete  
    … …  
    150150in the second case, $w$ is the velocity normal to the $s$-surfaces.  
    151151Note also that the $k$-axis is re-orientated downwards in the \textsc{fortran} code compared  
    152 to the indexing used in the semi-discrete equations such as \eqref{Eq_wzv}  
    153 (see  \S\ref{DOM_Num_Index_vertical}).  
     152to the indexing used in the semi-discrete equations such as \autoref{eq:wzv}  
     153(see \autoref{subsec:DOM_Num_Index_vertical}).  
    154154 
    155155 
    … …  
    158158% ================================================================ 
    159159\section{Coriolis and advection: vector invariant form} 
    160 \label{DYN_adv_cor_vect} 
     160\label{sec:DYN_adv_cor_vect} 
    161161%-----------------------------------------nam_dynadv---------------------------------------------------- 
    162162\forfile{../namelists/namdyn_adv}  
    … …  
    171171time (\textit{now} velocity).  
    172172At the lateral boundaries either free slip, no slip or partial slip boundary  
    173 conditions are applied following Chap.\ref{LBC}. 
     173conditions are applied following \autoref{chap:LBC}. 
    174174 
    175175% ------------------------------------------------------------------------------------------------------------- 
    … …  
    177177% ------------------------------------------------------------------------------------------------------------- 
    178178\subsection{Vorticity term (\protect\mdl{dynvor})} 
    179 \label{DYN_vor} 
     179\label{subsec:DYN_vor} 
    180180%------------------------------------------nam_dynvor---------------------------------------------------- 
    181181\forfile{../namelists/namdyn_vor}  
    … …  
    188188the relative vorticity term and horizontal kinetic energy for the planetary vorticity  
    189189term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent  
    190 flow and horizontal kinetic energy (EEN scheme) (see  Appendix~\ref{Apdx_C_vorEEN}). In the  
     190flow and horizontal kinetic energy (EEN scheme) (see \autoref{subsec:C_vorEEN}). In the  
    191191case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the  
    192192consistency of vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{ = .true.}). 
    … …  
    198198%------------------------------------------------------------- 
    199199\subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 
    200 \label{DYN_vor_ens} 
     200\label{subsec:DYN_vor_ens} 
    201201 
    202202In the enstrophy conserving case (ENS scheme), the discrete formulation of the  
    … …  
    204204($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent  
    205205flow ($i.e.$ $\chi$=$0$), but does not conserve the total kinetic energy. It is given by: 
    206 \begin{equation} \label{Eq_dynvor_ens} 
     206\begin{equation} \label{eq:dynvor_ens} 
    207207\left\{  
    208208\begin{aligned} 
    … …  
    219219%------------------------------------------------------------- 
    220220\subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 
    221 \label{DYN_vor_ene} 
     221\label{subsec:DYN_vor_ene} 
    222222 
    223223The kinetic energy conserving scheme (ENE scheme) conserves the global  
    224224kinetic energy but not the global enstrophy. It is given by: 
    225 \begin{equation} \label{Eq_dynvor_ene} 
     225\begin{equation} \label{eq:dynvor_ene} 
    226226\left\{   \begin{aligned} 
    227227{+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) 
    … …  
    236236%------------------------------------------------------------- 
    237237\subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) } 
    238 \label{DYN_vor_mix} 
     238\label{subsec:DYN_vor_mix} 
    239239 
    240240For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the  
    241 two previous schemes is used. It consists of the ENS scheme (\ref{Eq_dynvor_ens})  
    242 for the relative vorticity term, and of the ENE scheme (\ref{Eq_dynvor_ene}) applied  
     241two previous schemes is used. It consists of the ENS scheme (\autoref{eq:dynvor_ens})  
     242for the relative vorticity term, and of the ENE scheme (\autoref{eq:dynvor_ene}) applied  
    243243to the planetary vorticity term. 
    244 \begin{equation} \label{Eq_dynvor_mix} 
     244\begin{equation} \label{eq:dynvor_mix} 
    245245\left\{ {     \begin{aligned} 
    246246 {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i}  
    … …  
    259259%------------------------------------------------------------- 
    260260\subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) } 
    261 \label{DYN_vor_een} 
     261\label{subsec:DYN_vor_een} 
    262262 
    263263In both the ENS and ENE schemes, it is apparent that the combination of $i$ and $j$  
    … …  
    277277The idea is to get rid of the double averaging by considering triad combinations of vorticity.  
    278278It is noteworthy that this solution is conceptually quite similar to the one proposed by 
    279 \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see App.\ref{Apdx_C}). 
     279\citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). 
    280280 
    281281The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified  
    282282for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.  
    283283First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:  
    284 \begin{equation} \label{Eq_pot_vor} 
     284\begin{equation} \label{eq:pot_vor} 
    285285q  = \frac{\zeta +f} {e_{3f} } 
    286286\end{equation} 
    287 where the relative vorticity is defined by (\ref{Eq_divcur_cur}), the Coriolis parameter  
     287where the relative vorticity is defined by (\autoref{eq:divcur_cur}), the Coriolis parameter  
    288288is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:  
    289 \begin{equation} \label{Eq_een_e3f} 
     289\begin{equation} \label{eq:een_e3f} 
    290290e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} 
    291291\end{equation} 
    … …  
    294294\begin{figure}[!ht]    \begin{center} 
    295295\includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad} 
    296 \caption{ \protect\label{Fig_DYN_een_triad}   
     296\caption{ \protect\label{fig:DYN_een_triad}   
    297297Triads used in the energy and enstrophy conserving scheme (een) for  
    298298$u$-component (upper panel) and $v$-component (lower panel).} 
    … …  
    300300%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    301301 
    302 A key point in \eqref{Eq_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.  
     302A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.  
    303303It uses the sum of masked t-point vertical scale factor divided either  
    304304by the sum of the four t-point masks (\np{nn\_een\_e3f}\forcode{ = 1}),  
    … …  
    312312Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as  
    313313the following triad combinations of the neighbouring potential vorticities defined at f-points  
    314 (Fig.~\ref{Fig_DYN_een_triad}):  
    315 \begin{equation} \label{Q_triads} 
     314(\autoref{fig:DYN_een_triad}):  
     315\begin{equation} \label{eq:Q_triads} 
    316316_i^j \mathbb{Q}^{i_p}_{j_p} 
    317317= \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right) 
    … …  
    320320 
    321321Finally, the vorticity terms are represented as:  
    322 \begin{equation} \label{Eq_dynvor_een} 
     322\begin{equation} \label{eq:dynvor_een} 
    323323\left\{ { 
    324324\begin{aligned} 
    … …  
    333333This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes.  
    334334It conserves both total energy and potential enstrophy in the limit of horizontally  
    335 nondivergent flow ($i.e.$ $\chi$=$0$) (see  Appendix~\ref{Apdx_C_vorEEN}).  
     335nondivergent flow ($i.e.$ $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).  
    336336Applied to a realistic ocean configuration, it has been shown that it leads to a significant  
    337337reduction of the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}.  
    … …  
    344344%-------------------------------------------------------------------------------------------------------------- 
    345345\subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})} 
    346 \label{DYN_keg} 
    347  
    348 As demonstrated in Appendix~\ref{Apdx_C}, there is a single discrete formulation  
     346\label{subsec:DYN_keg} 
     347 
     348As demonstrated in \autoref{apdx:C}, there is a single discrete formulation  
    349349of the kinetic energy gradient term that, together with the formulation chosen for  
    350350the vertical advection (see below), conserves the total kinetic energy: 
    351 \begin{equation} \label{Eq_dynkeg} 
     351\begin{equation} \label{eq:dynkeg} 
    352352\left\{ \begin{aligned} 
    353353 -\frac{1}{2 \; e_{1u} }  & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   \\ 
    … …  
    360360%-------------------------------------------------------------------------------------------------------------- 
    361361\subsection{Vertical advection term (\protect\mdl{dynzad}) } 
    362 \label{DYN_zad} 
     362\label{subsec:DYN_zad} 
    363363 
    364364The discrete formulation of the vertical advection, together with the formulation  
    365365chosen for the gradient of kinetic energy (KE) term, conserves the total kinetic  
    366366energy. Indeed, the change of KE due to the vertical advection is exactly  
    367 balanced by the change of KE due to the gradient of KE (see Appendix~\ref{Apdx_C}). 
    368 \begin{equation} \label{Eq_dynzad} 
     367balanced by the change of KE due to the gradient of KE (see \autoref{apdx:C}). 
     368\begin{equation} \label{eq:dynzad} 
    369369\left\{     \begin{aligned} 
    370370-\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2}  \;\delta _{k+1/2} \left[ u \right]\  }^{\,k}  \\ 
    … …  
    377377Note that in this case, a similar split-explicit time stepping should be used on  
    378378vertical advection of tracer to ensure a better stability,  
    379 an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \S\ref{TRA_adv_tvd}). 
     379an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). 
    380380 
    381381 
    … …  
    384384% ================================================================ 
    385385\section{Coriolis and advection: flux form} 
    386 \label{DYN_adv_cor_flux} 
     386\label{sec:DYN_adv_cor_flux} 
    387387%------------------------------------------nam_dynadv---------------------------------------------------- 
    388388\forfile{../namelists/namdyn_adv}  
    … …  
    394394appearing in their expressions is centred in time (\textit{now} velocity). At the  
    395395lateral boundaries either free slip, no slip or partial slip boundary conditions  
    396 are applied following Chap.\ref{LBC}. 
     396are applied following \autoref{chap:LBC}. 
    397397 
    398398 
    … …  
    401401%-------------------------------------------------------------------------------------------------------------- 
    402402\subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) } 
    403 \label{DYN_cor_flux} 
     403\label{subsec:DYN_cor_flux} 
    404404 
    405405In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis  
    406406parameter has been modified to account for the "metric" term. This altered  
    407407Coriolis parameter is thus discretised at $f$-points. It is given by:  
    408 \begin{multline} \label{Eq_dyncor_metric} 
     408\begin{multline} \label{eq:dyncor_metric} 
    409409f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\ 
    410410   \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right]   
    … …  
    412412\end{multline}  
    413413 
    414 Any of the (\ref{Eq_dynvor_ens}), (\ref{Eq_dynvor_ene}) and (\ref{Eq_dynvor_een})  
     414Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een})  
    415415schemes can be used to compute the product of the Coriolis parameter and the  
    416 vorticity. However, the energy-conserving scheme  (\ref{Eq_dynvor_een}) has  
     416vorticity. However, the energy-conserving scheme  (\autoref{eq:dynvor_een}) has  
    417417exclusively been used to date. This term is evaluated using a leapfrog scheme,  
    418418$i.e.$ the velocity is centred in time (\textit{now} velocity). 
    … …  
    422422%-------------------------------------------------------------------------------------------------------------- 
    423423\subsection{Flux form advection term (\protect\mdl{dynadv}) } 
    424 \label{DYN_adv_flux} 
     424\label{subsec:DYN_adv_flux} 
    425425 
    426426The discrete expression of the advection term is given by : 
    427 \begin{equation} \label{Eq_dynadv} 
     427\begin{equation} \label{eq:dynadv} 
    428428\left\{  
    429429\begin{aligned} 
    … …  
    454454%------------------------------------------------------------- 
    455455\subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 
    456 \label{DYN_adv_cen2} 
     456\label{subsec:DYN_adv_cen2} 
    457457 
    458458In the centered $2^{nd}$ order formulation, the velocity is evaluated as the  
    459459mean of the two neighbouring points : 
    460 \begin{equation} \label{Eq_dynadv_cen2} 
     460\begin{equation} \label{eq:dynadv_cen2} 
    461461\left\{     \begin{aligned} 
    462462 u_T^{cen2} &=\overline u^{i }       \quad &  u_F^{cen2} &=\overline u^{j+1/2}  \quad &  u_{uw}^{cen2} &=\overline u^{k+1/2}   \\ 
    … …  
    475475%------------------------------------------------------------- 
    476476\subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 
    477 \label{DYN_adv_ubs} 
     477\label{subsec:DYN_adv_ubs} 
    478478 
    479479The UBS advection scheme is an upstream biased third order scheme based on  
    480480an upstream-biased parabolic interpolation. For example, the evaluation of  
    481481$u_T^{ubs} $ is done as follows: 
    482 \begin{equation} \label{Eq_dynadv_ubs} 
     482\begin{equation} \label{eq:dynadv_ubs} 
    483483u_T^{ubs} =\overline u ^i-\;\frac{1}{6}   \begin{cases} 
    484484      u"_{i-1/2}&    \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  \geqslant 0$ }    \\ 
    … …  
    498498The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$  
    499499order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and  
    500 $u_{vw}^{ubs}$ in \eqref{Eq_dynadv_cen2} are used. UBS is diffusive and is  
     500$u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used. UBS is diffusive and is  
    501501associated with vertical mixing of momentum. \gmcomment{ gm  pursue the  
    502502sentence:Since vertical mixing of momentum is a source term of the TKE equation...  } 
    503503 
    504 For stability reasons,  the first term in (\ref{Eq_dynadv_ubs}), which corresponds  
     504For stability reasons, the first term in (\autoref{eq:dynadv_ubs}), which corresponds  
    505505to a second order centred scheme, is evaluated using the \textit{now} velocity  
    506506(centred in time), while the second term, which is the diffusion part of the scheme,  
    … …  
    510510Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics)  
    511511schemes only differ by one coefficient. Replacing $1/6$ by $1/8$ in  
    512 (\ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.  
     512(\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.  
    513513This option is not available through a namelist parameter, since the $1/6$ coefficient  
    514514is hard coded. Nevertheless it is quite easy to make the substitution in the 
    … …  
    526526% ================================================================ 
    527527\section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
    528 \label{DYN_hpg} 
     528\label{sec:DYN_hpg} 
    529529%------------------------------------------nam_dynhpg--------------------------------------------------- 
    530530\forfile{../namelists/namdyn_hpg}  
    … …  
    547547%-------------------------------------------------------------------------------------------------------------- 
    548548\subsection{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 
    549 \label{DYN_hpg_zco} 
     549\label{subsec:DYN_hpg_zco} 
    550550 
    551551The hydrostatic pressure can be obtained by integrating the hydrostatic equation  
    … …  
    556556 
    557557for $k=km$ (surface layer, $jk=1$ in the code) 
    558 \begin{equation} \label{Eq_dynhpg_zco_surf} 
     558\begin{equation} \label{eq:dynhpg_zco_surf} 
    559559\left\{ \begin{aligned} 
    560560               \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k=km}  
    … …  
    566566 
    567567for $1<k<km$ (interior layer) 
    568 \begin{equation} \label{Eq_dynhpg_zco} 
     568\begin{equation} \label{eq:dynhpg_zco} 
    569569\left\{ \begin{aligned} 
    570570               \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k}  
    … …  
    577577\end{equation}  
    578578 
    579 Note that the $1/2$ factor in (\ref{Eq_dynhpg_zco_surf}) is adequate because of  
     579Note that the $1/2$ factor in (\autoref{eq:dynhpg_zco_surf}) is adequate because of  
    580580the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface  
    581581level ($z=0$). Note also that in case of variable volume level (\key{vvl} defined), the  
    582 surface pressure gradient is included in \eqref{Eq_dynhpg_zco_surf} and \eqref{Eq_dynhpg_zco}  
     582surface pressure gradient is included in \autoref{eq:dynhpg_zco_surf} and \autoref{eq:dynhpg_zco}  
    583583through the space and time variations of the vertical scale factor $e_{3w}$. 
    584584 
    … …  
    587587%-------------------------------------------------------------------------------------------------------------- 
    588588\subsection{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 
    589 \label{DYN_hpg_zps} 
     589\label{subsec:DYN_hpg_zps} 
    590590 
    591591With partial bottom cells, tracers in horizontally adjacent cells generally live at  
    … …  
    596596Apart from this modification, the horizontal hydrostatic pressure gradient evaluated  
    597597in the $z$-coordinate with partial step is exactly as in the pure $z$-coordinate case.  
    598 As explained in detail in section \S\ref{TRA_zpshde}, the nonlinearity of pressure  
     598As explained in detail in section \autoref{sec:TRA_zpshde}, the nonlinearity of pressure  
    599599effects in the equation of state is such that it is better to interpolate temperature and  
    600600salinity vertically before computing the density. Horizontal gradients of temperature  
    601601and salinity are needed for the TRA modules, which is the reason why the horizontal  
    602602gradients of density at the deepest model level are computed in module \mdl{zpsdhe}  
    603 located in the TRA directory and described in \S\ref{TRA_zpshde}. 
     603located in the TRA directory and described in \autoref{sec:TRA_zpshde}. 
    604604 
    605605%-------------------------------------------------------------------------------------------------------------- 
    … …  
    607607%-------------------------------------------------------------------------------------------------------------- 
    608608\subsection{$S$- and $Z$-$S$-coordinates} 
    609 \label{DYN_hpg_sco} 
     609\label{subsec:DYN_hpg_sco} 
    610610 
    611611Pressure gradient formulations in an $s$-coordinate have been the subject of a vast  
    … …  
    615615 
    616616$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.}) 
    617 \begin{equation} \label{Eq_dynhpg_sco} 
     617\begin{equation} \label{eq:dynhpg_sco} 
    618618\left\{ \begin{aligned} 
    619619 - \frac{1}                   {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  p^h  \right]  
    … …  
    625625 
    626626Where the first term is the pressure gradient along coordinates, computed as in  
    627 \eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of  
     627\autoref{eq:dynhpg_zco_surf} - \autoref{eq:dynhpg_zco}, and $z_T$ is the depth of  
    628628the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point  
    629629($e_{3w}$). 
    … …  
    637637(\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development) 
    638638 
    639 Note that expression \eqref{Eq_dynhpg_sco} is commonly used when the variable volume formulation is 
     639Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is 
    640640activated (\key{vvl}) because in that case, even with a flat bottom, the coordinate surfaces are not 
    641641horizontal but follow the free surface \citep{Levier2007}. The pressure jacobian scheme 
    … …  
    648648 
    649649\subsection{Ice shelf cavity} 
    650 \label{DYN_hpg_isf} 
     650\label{subsec:DYN_hpg_isf} 
    651651Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and 
    652652 the pressure gradient due to the ocean load. If cavity opened (\np{ln\_isfcav}\forcode{ = .true.}) these 2 terms can be 
    … …  
    658658This top pressure is constant over time. A detailed description of this method is described in \citet{Losch2008}.\\ 
    659659 
    660 $\bullet$ The ocean load is computed using the expression \eqref{Eq_dynhpg_sco} described in \ref{DYN_hpg_sco}.  
     660$\bullet$ The ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in \autoref{subsec:DYN_hpg_sco}.  
    661661 
    662662%-------------------------------------------------------------------------------------------------------------- 
    … …  
    664664%-------------------------------------------------------------------------------------------------------------- 
    665665\subsection{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})} 
    666 \label{DYN_hpg_imp} 
     666\label{subsec:DYN_hpg_imp} 
    667667 
    668668The default time differencing scheme used for the horizontal pressure gradient is  
    … …  
    680680$\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}): 
    681681 
    682 \begin{equation} \label{Eq_dynhpg_lf} 
     682\begin{equation} \label{eq:dynhpg_lf} 
    683683\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 
    684684   -\frac{1}{\rho _o \,e_{1u} }\delta _{i+1/2} \left[ {p_h^t } \right] 
    … …  
    686686 
    687687$\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}): 
    688 \begin{equation} \label{Eq_dynhpg_imp} 
     688\begin{equation} \label{eq:dynhpg_imp} 
    689689\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 
    690690   -\frac{1}{4\,\rho _o \,e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2\,p_h^t +p_h^{t-\rdt}  \right] 
    691691\end{equation} 
    692692 
    693 The semi-implicit time scheme \eqref{Eq_dynhpg_imp} is made possible without  
     693The semi-implicit time scheme \autoref{eq:dynhpg_imp} is made possible without  
    694694significant additional computation since the density can be updated to time level  
    695695$t+\rdt$ before computing the horizontal hydrostatic pressure gradient. It can  
    696696be easily shown that the stability limit associated with the hydrostatic pressure  
    697 gradient doubles using \eqref{Eq_dynhpg_imp} compared to that using the  
    698 standard leapfrog scheme \eqref{Eq_dynhpg_lf}. Note that \eqref{Eq_dynhpg_imp}  
     697gradient doubles using \autoref{eq:dynhpg_imp} compared to that using the  
     698standard leapfrog scheme \autoref{eq:dynhpg_lf}. Note that \autoref{eq:dynhpg_imp}  
    699699is equivalent to applying a time filter to the pressure gradient to eliminate high  
    700 frequency IGWs. Obviously, when using \eqref{Eq_dynhpg_imp}, the doubling of  
     700frequency IGWs. Obviously, when using \autoref{eq:dynhpg_imp}, the doubling of  
    701701the time-step is achievable only if no other factors control the time-step, such as  
    702702the stability limits associated with advection or diffusion. 
    … …  
    708708compute the hydrostatic pressure gradient (whatever the formulation) is evaluated  
    709709as follows: 
    710 \begin{equation} \label{Eq_rho_flt} 
     710\begin{equation} \label{eq:rho_flt} 
    711711   \rho^t = \rho( \widetilde{T},\widetilde {S},z_t) 
    712712 \quad     \text{with}  \quad  
    … …  
    722722% ================================================================ 
    723723\section{Surface pressure gradient (\protect\mdl{dynspg})} 
    724 \label{DYN_spg} 
     724\label{sec:DYN_spg} 
    725725%-----------------------------------------nam_dynspg---------------------------------------------------- 
    726726\forfile{../namelists/namdyn_spg}  
    … …  
    730730 
    731731Options are defined through the \ngn{namdyn\_spg} namelist variables. 
    732 The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}).  
     732The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:PE_hor_pg}).  
    733733The main distinction is between the fixed volume case (linear free surface) and the variable volume case  
    734 (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface})  
     734(nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\autoref{subsec:PE_free_surface})  
    735735the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case  
    736 (\S\ref{PE_free_surface}).  
     736(\autoref{subsec:PE_free_surface}).  
    737737With both linear and nonlinear free surface, external gravity waves are allowed in the equations,  
    738738which imposes a very small time step when an explicit time stepping is used.  
    739739Two methods are proposed to allow a longer time step for the three-dimensional equations:  
    740 the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}),  
     740the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt}),  
    741741and the split-explicit free surface described below.  
    742742The extra term introduced in the filtered method is calculated implicitly,  
    … …  
    745745 
    746746The form of the surface pressure gradient term depends on how the user wants to handle  
    747 the fast external gravity waves that are a solution of the analytical equation (\S\ref{PE_hor_pg}).  
     747the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:PE_hor_pg}).  
    748748Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx): 
    749749an explicit formulation which requires a small time step ; 
    … …  
    761761%-------------------------------------------------------------------------------------------------------------- 
    762762\subsection{Explicit free surface (\protect\key{dynspg\_exp})} 
    763 \label{DYN_spg_exp} 
     763\label{subsec:DYN_spg_exp} 
    764764 
    765765In the explicit free surface formulation (\key{dynspg\_exp} defined), the model time step  
    … …  
    767767The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time), 
    768768is thus simply given by : 
    769 \begin{equation} \label{Eq_dynspg_exp} 
     769\begin{equation} \label{eq:dynspg_exp} 
    770770\left\{ \begin{aligned} 
    771771 - \frac{1}{e_{1u}\,\rho_o} \;   \delta _{i+1/2} \left[  \,\rho \,\eta\,  \right]   \\ 
    … …  
    782782%-------------------------------------------------------------------------------------------------------------- 
    783783\subsection{Split-explicit free surface (\protect\key{dynspg\_ts})} 
    784 \label{DYN_spg_ts} 
     784\label{subsec:DYN_spg_ts} 
    785785%------------------------------------------namsplit----------------------------------------------------------- 
    786786%\forfile{../namelists/namsplit} 
    … …  
    792792equation and the associated barotropic velocity equations with a smaller time  
    793793step than $\rdt$, the time step used for the three dimensional prognostic  
    794 variables (Fig.~\ref{Fig_DYN_dynspg_ts}).  
     794variables (\autoref{fig:DYN_dynspg_ts}).  
    795795The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) 
    796796 is provided through the \np{nn\_baro} namelist parameter as:  
    … …  
    802802%%% 
    803803The barotropic mode solves the following equations: 
    804 \begin{subequations} \label{Eq_BT} 
    805   \begin{equation}     \label{Eq_BT_dyn} 
     804\begin{subequations} \label{eq:BT} 
     805  \begin{equation}     \label{eq:BT_dyn} 
    806806\frac{\partial {\rm \overline{{\bf U}}_h} }{\partial t}= 
    807807 -f\;{\rm {\bf k}}\times {\rm \overline{{\bf U}}_h}  
    … …  
    809809  \end{equation} 
    810810 
    811   \begin{equation} \label{Eq_BT_ssh} 
     811  \begin{equation} \label{eq:BT_ssh} 
    812812\frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]+P-E 
    813813  \end{equation} 
    814814\end{subequations} 
    815 where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. The third term on the right hand side of \eqref{Eq_BT_dyn} represents the bottom stress (see section \S\ref{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}. AB3-AM4 coefficients used in \NEMO follow the second-order accurate, "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08} (see their figure 12, lower left).  
     815where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress (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}. AB3-AM4 coefficients used in \NEMO follow the second-order accurate, "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08} (see their figure 12, lower left).  
    816816 
    817817%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   > 
    818818\begin{figure}[!t]    \begin{center} 
    819819\includegraphics[width=0.7\textwidth]{Fig_DYN_dynspg_ts} 
    820 \caption{  \protect\label{Fig_DYN_dynspg_ts} 
     820\caption{  \protect\label{fig:DYN_dynspg_ts} 
    821821Schematic of the split-explicit time stepping scheme for the external  
    822822and internal modes. Time increases to the right. In this particular exemple,  
    … …  
    827827The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged  
    828828transports to advect tracers. 
    829 a) Forward time integration: \np{ln\_bt\_fw}\forcode{ = .true.},  \np{ln\_bt\_av}\forcode{ = .true.}. 
    830 b) Centred time integration: \np{ln\_bt\_fw}\forcode{ = .false.}, \np{ln\_bt\_av}\forcode{ = .true.}. 
    831 c) Forward time integration with no time filtering (POM-like scheme): \np{ln\_bt\_fw}\forcode{ = .true.}, \np{ln\_bt\_av}\forcode{ = .false.}. } 
     829a) Forward time integration: \protect\np{ln\_bt\_fw}\forcode{ = .true.},  \protect\np{ln\_bt\_av}\forcode{ = .true.}. 
     830b) Centred time integration: \protect\np{ln\_bt\_fw}\forcode{ = .false.}, \protect\np{ln\_bt\_av}\forcode{ = .true.}. 
     831c) Forward time integration with no time filtering (POM-like scheme): \protect\np{ln\_bt\_fw}\forcode{ = .true.}, \protect\np{ln\_bt\_av}\forcode{ = .false.}. } 
    832832\end{center}    \end{figure} 
    833833%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   > 
    834834 
    835835In the default case (\np{ln\_bt\_fw}\forcode{ = .true.}), the external mode is integrated  
    836 between \textit{now} and  \textit{after} baroclinic time-steps (Fig.~\ref{Fig_DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic  
     836between \textit{now} and  \textit{after} baroclinic time-steps (\autoref{fig:DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic  
    837837quantities (\np{ln\_bt\_av}\forcode{ = .true.}). In that case, the integration is extended slightly beyond  \textit{after} time step to provide time filtered quantities.  
    838838These are used for the subsequent initialization of the barotropic mode in the following baroclinic step.  
    … …  
    850850at \textit{now} time step. This implies to change the traditional order of computations in \NEMO: most of momentum   
    851851trends (including the barotropic mode calculation) updated first, tracers' after. This \textit{de facto} makes semi-implicit hydrostatic  
    852 pressure gradient (see section \S\ref{DYN_hpg_imp}) and time splitting not compatible.  
     852pressure gradient (see section \autoref{subsec:DYN_hpg_imp}) and time splitting not compatible.  
    853853Advective barotropic velocities are obtained by using a secondary set of filtering weights, uniquely defined from the filter  
    854854coefficients 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. 
    … …  
    872872scheme using the small barotropic time step $\rdt$. We have  
    873873 
    874 \begin{equation} \label{DYN_spg_ts_eta} 
     874\begin{equation} \label{eq:DYN_spg_ts_eta} 
    875875\eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1}) 
    876876   = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right]  
    877877\end{equation} 
    878 \begin{multline} \label{DYN_spg_ts_u} 
     878\begin{multline} \label{eq:DYN_spg_ts_u} 
    879879\textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1})  \\ 
    880880   = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n})  
    … …  
    886886and $U^{(b)}$ denotes the baroclinic time at which the vertically integrated forcing $\textbf{M}(\tau)$ (note that this forcing includes the surface freshwater forcing), the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$, and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over a single cycle. This is also the time  
    887887that sets the barotropic time steps via  
    888 \begin{equation} \label{DYN_spg_ts_t} 
     888\begin{equation} \label{eq:DYN_spg_ts_t} 
    889889t_n=\tau+n\rdt    
    890890\end{equation} 
    891891with $n$ an integer. The density scaled surface pressure is evaluated via  
    892 \begin{equation} \label{DYN_spg_ts_ps} 
     892\begin{equation} \label{eq:DYN_spg_ts_ps} 
    893893p_s^{(b)}(\tau,t_{n}) = \begin{cases} 
    894894   g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o  &      \text{non-linear case} \\ 
    … …  
    897897\end{equation} 
    898898To get started, we assume the following initial conditions  
    899 \begin{equation} \label{DYN_spg_ts_eta} 
     899\begin{equation} \label{eq:DYN_spg_ts_eta} 
    900900\begin{split} 
    901901\eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)} 
    … …  
    905905\end{equation} 
    906906with  
    907 \begin{equation} \label{DYN_spg_ts_etaF} 
     907\begin{equation} \label{eq:DYN_spg_ts_etaF} 
    908908 \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n}) 
    909909\end{equation} 
    910910the time averaged surface height taken from the previous barotropic cycle. Likewise,  
    911 \begin{equation} \label{DYN_spg_ts_u} 
     911\begin{equation} \label{eq:DYN_spg_ts_u} 
    912912\textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)}   \\ 
    913913\\ 
    … …  
    915915\end{equation} 
    916916with  
    917 \begin{equation} \label{DYN_spg_ts_u} 
     917\begin{equation} \label{eq:DYN_spg_ts_u} 
    918918 \overline{\textbf{U}^{(b)}(\tau)}  
    919919   = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n}) 
    … …  
    922922 
    923923Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ , the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at baroclinic time $\tau + \rdt \tau$  
    924 \begin{equation} \label{DYN_spg_ts_u} 
     924\begin{equation} \label{eq:DYN_spg_ts_u} 
    925925\textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)}  
    926926   = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n}) 
    … …  
    928928The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using the following form  
    929929 
    930 \begin{equation} \label{DYN_spg_ts_ssh} 
     930\begin{equation} \label{eq:DYN_spg_ts_ssh} 
    931931\eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right]   
    932932\end{equation} 
    … …  
    935935  
    936936In general, some form of time filter is needed to maintain integrity of the surface  
    937 height field due to the leap-frog splitting mode in equation \ref{DYN_spg_ts_ssh}. We  
     937height field due to the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}. We  
    938938have tried various forms of such filtering, with the following method discussed in  
    939939\cite{Griffies_al_MWR01} chosen due to its stability and reasonably good maintenance of  
    940 tracer conservation properties (see Section ??)  
    941  
    942 \begin{equation} \label{DYN_spg_ts_sshf} 
     940tracer conservation properties (see ??)  
     941 
     942\begin{equation} \label{eq:DYN_spg_ts_sshf} 
    943943\eta^{F}(\tau-\Delta) =  \overline{\eta^{(b)}(\tau)}  
    944944\end{equation} 
    945945Another approach tried was  
    946946 
    947 \begin{equation} \label{DYN_spg_ts_sshf2} 
     947\begin{equation} \label{eq:DYN_spg_ts_sshf2} 
    948948\eta^{F}(\tau-\Delta) = \eta(\tau)  
    949949   + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt) 
    … …  
    953953which is useful since it isolates all the time filtering aspects into the term multiplied  
    954954by $\alpha$. This isolation allows for an easy check that tracer conservation is exact when  
    955 eliminating tracer and surface height time filtering (see Section ?? for more complete discussion). However, in the general case with a non-zero $\alpha$, the filter \ref{DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.  
     955eliminating tracer and surface height time filtering (see ?? for more complete discussion). However, in the general case with a non-zero $\alpha$, the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.  
    956956 
    957957}            %%end gm comment (copy of griffies book) 
    … …  
    964964%-------------------------------------------------------------------------------------------------------------- 
    965965\subsection{Filtered free surface (\protect\key{dynspg\_flt})} 
    966 \label{DYN_spg_fltp} 
     966\label{subsec:DYN_spg_fltp} 
    967967 
    968968The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation.  
    969 The extra term introduced in the equations (see \S\ref{PE_free_surface}) is solved implicitly.  
    970 The elliptic solvers available in the code are documented in \S\ref{MISC}. 
     969The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly.  
     970The elliptic solvers available in the code are documented in \autoref{chap:MISC}. 
    971971 
    972972%% gm %%======>>>>   given here the discrete eqs provided to the solver 
    973973\gmcomment{               %%% copy from chap-model basics  
    974 \begin{equation} \label{Eq_spg_flt} 
     974\begin{equation} \label{eq:spg_flt} 
    975975\frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} 
    976976- g \nabla \left( \tilde{\rho} \ \eta \right)  
    … …  
    980980$\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$  
    981981represents the collected contributions of the Coriolis, hydrostatic pressure gradient,  
    982 non-linear and viscous terms in \eqref{Eq_PE_dyn}. 
     982non-linear and viscous terms in \autoref{eq:PE_dyn}. 
    983983}   %end gmcomment 
    984984 
    … …  
    990990% ================================================================ 
    991991\section{Lateral diffusion term and operators (\protect\mdl{dynldf})} 
    992 \label{DYN_ldf} 
     992\label{sec:DYN_ldf} 
    993993%------------------------------------------nam_dynldf---------------------------------------------------- 
    994994\forfile{../namelists/namdyn_ldf}  
    … …  
    999999(rotated or not) or biharmonic operators. The coefficients may be constant  
    10001000or spatially variable; the description of the coefficients is found in the chapter  
    1001 on lateral physics (Chap.\ref{LDF}). The lateral diffusion of momentum is  
     1001on lateral physics (\autoref{chap:LDF}). The lateral diffusion of momentum is  
    10021002evaluated using a forward scheme, $i.e.$ the velocity appearing in its expression  
    10031003is the \textit{before} velocity in time, except for the pure vertical component  
    10041004that appears when a tensor of rotation is used. This latter term is solved  
    1005 implicitly together with the vertical diffusion term (see \S\ref{STP})  
     1005implicitly together with the vertical diffusion term (see \autoref{chap:STP})  
    10061006 
    10071007At the lateral boundaries either free slip, no slip or partial slip boundary  
    1008 conditions are applied according to the user's choice (see Chap.\ref{LBC}). 
     1008conditions are applied according to the user's choice (see \autoref{chap:LBC}). 
    10091009 
    10101010\gmcomment{ 
    … …  
    10251025\subsection[Iso-level laplacian (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})] 
    10261026            {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 
    1027 \label{DYN_ldf_lap} 
     1027\label{subsec:DYN_ldf_lap} 
    10281028 
    10291029For lateral iso-level diffusion, the discrete operator is:  
    1030 \begin{equation} \label{Eq_dynldf_lap} 
     1030\begin{equation} \label{eq:dynldf_lap} 
    10311031\left\{ \begin{aligned} 
    10321032 D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta _{i+1/2} \left[ {A_T^{lm}  
    … …  
    10401040\end{equation}  
    10411041 
    1042 As explained in \S\ref{PE_ldf}, this formulation (as the gradient of a divergence  
     1042As explained in \autoref{subsec:PE_ldf}, this formulation (as the gradient of a divergence  
    10431043and curl of the vorticity) preserves symmetry and ensures a complete  
    10441044separation between the vorticity and divergence parts of the momentum diffusion.  
    … …  
    10491049\subsection[Rotated laplacian (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})] 
    10501050            {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 
    1051 \label{DYN_ldf_iso} 
     1051\label{subsec:DYN_ldf_iso} 
    10521052 
    10531053A rotation of the lateral momentum diffusion operator is needed in several cases:  
    … …  
    10611061constraints on the stress tensor such as symmetry. The resulting discrete  
    10621062representation is: 
    1063 \begin{equation} \label{Eq_dyn_ldf_iso} 
     1063\begin{equation} \label{eq:dyn_ldf_iso} 
    10641064\begin{split} 
    10651065 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ 
    … …  
    11111111diffusion operator acts and the surface of computation ($z$- or $s$-surfaces).  
    11121112The way these slopes are evaluated is given in the lateral physics chapter  
    1113 (Chap.\ref{LDF}). 
     1113(\autoref{chap:LDF}). 
    11141114 
    11151115%-------------------------------------------------------------------------------------------------------------- 
    … …  
    11181118\subsection[Iso-level bilaplacian (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})] 
    11191119            {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 
    1120 \label{DYN_ldf_bilap} 
     1120\label{subsec:DYN_ldf_bilap} 
    11211121 
    11221122The lateral fourth order operator formulation on momentum is obtained by  
    1123 applying \eqref{Eq_dynldf_lap} twice. It requires an additional assumption on  
     1123applying \autoref{eq:dynldf_lap} twice. It requires an additional assumption on  
    11241124boundary conditions: the first derivative term normal to the coast depends on  
    11251125the free or no-slip lateral boundary conditions chosen, while the third  
    1126 derivative terms normal to the coast are set to zero (see Chap.\ref{LBC}). 
     1126derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}). 
    11271127%%% 
    11281128\gmcomment{add a remark on the the change in the position of the coefficient} 
    … …  
    11331133% ================================================================ 
    11341134\section{Vertical diffusion term (\protect\mdl{dynzdf})} 
    1135 \label{DYN_zdf} 
     1135\label{sec:DYN_zdf} 
    11361136%----------------------------------------------namzdf------------------------------------------------------ 
    11371137\forfile{../namelists/namzdf}  
    … …  
    11451145scheme (\np{ln\_zdfexp}\forcode{ = .true.}) using a time splitting technique  
    11461146(\np{nn\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme  
    1147 (\np{ln\_zdfexp}\forcode{ = .false.}) (see \S\ref{STP}). Note that namelist variables  
     1147(\np{ln\_zdfexp}\forcode{ = .false.}) (see \autoref{chap:STP}). Note that namelist variables  
    11481148\np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.  
    11491149 
    11501150The formulation of the vertical subgrid scale physics is the same whatever  
    11511151the vertical coordinate is. The vertical diffusion operators given by  
    1152 \eqref{Eq_PE_zdf} take the following semi-discrete space form: 
    1153 \begin{equation} \label{Eq_dynzdf} 
     1152\autoref{eq:PE_zdf} take the following semi-discrete space form: 
     1153\begin{equation} \label{eq:dynzdf} 
    11541154\left\{   \begin{aligned} 
    11551155D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta _k \left[ \frac{A_{uw}^{vm} }{e_{3uw} } 
    … …  
    11621162where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and  
    11631163diffusivity coefficients. The way these coefficients are evaluated  
    1164 depends on the vertical physics used (see \S\ref{ZDF}). 
     1164depends on the vertical physics used (see \autoref{chap:ZDF}). 
    11651165 
    11661166The surface boundary condition on momentum is the stress exerted by  
    11671167the wind. At the surface, the momentum fluxes are prescribed as the boundary  
    11681168condition on the vertical turbulent momentum fluxes, 
    1169 \begin{equation} \label{Eq_dynzdf_sbc} 
     1169\begin{equation} \label{eq:dynzdf_sbc} 
    11701170\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 
    11711171    = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v } 
    … …  
    11771177is small (when no mixed layer scheme is used) the surface stress enters only  
    11781178the top model level, as a body force. The surface wind stress is calculated  
    1179 in the surface module routines (SBC, see Chap.\ref{SBC}) 
     1179in the surface module routines (SBC, see \autoref{chap:SBC}) 
    11801180 
    11811181The turbulent flux of momentum at the bottom of the ocean is specified through  
    1182 a bottom friction parameterisation (see \S\ref{ZDF_bfr}) 
     1182a bottom friction parameterisation (see \autoref{sec:ZDF_bfr}) 
    11831183 
    11841184% ================================================================ 
    … …  
    11861186% ================================================================ 
    11871187\section{External forcings} 
    1188 \label{DYN_forcing} 
     1188\label{sec:DYN_forcing} 
    11891189 
    11901190Besides the surface and bottom stresses (see the above section) which are  
    … …  
    11921192may enter the dynamical equations by affecting the surface pressure gradient.  
    11931193 
    1194 (1) When \np{ln\_apr\_dyn}\forcode{ = .true.} (see \S\ref{SBC_apr}), the atmospheric pressure is taken  
     1194(1) When \np{ln\_apr\_dyn}\forcode{ = .true.} (see \autoref{sec:SBC_apr}), the atmospheric pressure is taken  
    11951195into account when computing the surface pressure gradient. 
    11961196 
    1197 (2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ = .true.} (see \S\ref{SBC_tide}),  
     1197(2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ = .true.} (see \autoref{sec:SBC_tide}),  
    11981198the tidal potential is taken into account when computing the surface pressure gradient. 
    11991199 
    … …  
    12091209% ================================================================ 
    12101210\section{Time evolution term (\protect\mdl{dynnxt})} 
    1211 \label{DYN_nxt} 
     1211\label{sec:DYN_nxt} 
    12121212 
    12131213%----------------------------------------------namdom---------------------------------------------------- 
    … …  
    12181218The general framework for dynamics time stepping is a leap-frog scheme,  
    12191219$i.e.$ a three level centred time scheme associated with an Asselin time filter  
    1220 (cf. Chap.\ref{STP}). The scheme is applied to the velocity, except when using  
    1221 the flux form of momentum advection (cf. \S\ref{DYN_adv_cor_flux}) in the variable  
     1220(cf. \autoref{chap:STP}). The scheme is applied to the velocity, except when using  
     1221the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux}) in the variable  
    12221222volume case (\key{vvl} defined), where it has to be applied to the thickness  
    1223 weighted velocity (see \S\ref{Apdx_A_momentum})   
     1223weighted velocity (see \autoref{sec:A_momentum})   
    12241224 
    12251225$\bullet$ vector invariant form or linear free surface (\np{ln\_dynhpg\_vec}\forcode{ = .true.} ; \key{vvl} not defined): 
    1226 \begin{equation} \label{Eq_dynnxt_vec} 
     1226\begin{equation} \label{eq:dynnxt_vec} 
    12271227\left\{   \begin{aligned} 
    12281228&u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t     \\ 
    … …  
    12321232 
    12331233$\bullet$ flux form and nonlinear free surface (\np{ln\_dynhpg\_vec}\forcode{ = .false.} ; \key{vvl} defined): 
    1234 \begin{equation} \label{Eq_dynnxt_flux} 
     1234\begin{equation} \label{eq:dynnxt_flux} 
    12351235\left\{   \begin{aligned} 
    12361236&\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t     \\ 

  • branches/2017/dev_merge_2017/DOC/tex_sub/chap_LBC.tex

    r9393 r9407  
    55% ================================================================ 
    66\chapter{Lateral Boundary Condition (LBC)} 
    7 \label{LBC} 
     7\label{chap:LBC} 
    88\minitoc 
    99 
    … …  
    1818% ================================================================ 
    1919\section{Boundary condition at the coast (\protect\np{rn\_shlat})} 
    20 \label{LBC_coast} 
     20\label{sec:LBC_coast} 
    2121%--------------------------------------------nam_lbc------------------------------------------------------- 
    2222\forfile{../namelists/namlbc}  
    2323%-------------------------------------------------------------------------------------------------------------- 
    2424 
    25 %The lateral ocean boundary conditions contiguous to coastlines are Neumann conditions for heat and salt (no flux across boundaries) and Dirichlet conditions for momentum (ranging from free-slip to "strong" no-slip). They are handled automatically by the mask system (see \S\ref{DOM_msk}).  
    26  
    27 %OPA allows land and topography grid points in the computational domain due to the presence of continents or islands, and includes the use of a full or partial step representation of bottom topography. The computation is performed over the whole domain, i.e. we do not try to restrict the computation to ocean-only points. This choice has two motivations. Firstly, working on ocean only grid points overloads the code and harms the code readability. Secondly, and more importantly, it drastically reduces the vector portion of the computation, leading to a dramatic increase of CPU time requirement on vector computers.  The current section describes how the masking affects the computation of the various terms of the equations with respect to the boundary condition at solid walls. The process of defining which areas are to be masked is described in \S\ref{DOM_msk}. 
     25%The lateral ocean boundary conditions contiguous to coastlines are Neumann conditions for heat and salt (no flux across boundaries) and Dirichlet conditions for momentum (ranging from free-slip to "strong" no-slip). They are handled automatically by the mask system (see \autoref{subsec:DOM_msk}).  
     26 
     27%OPA allows land and topography grid points in the computational domain due to the presence of continents or islands, and includes the use of a full or partial step representation of bottom topography. The computati