Changeset 11543


Ignore:
Timestamp:
2019-09-13T15:57:52+02:00 (13 months ago)
Author:
nicolasmartin
Message:

Implementation of convention for labelling references + files renaming
Now each reference is supposed to have the information of the chapter in its name
to identify quickly which file contains the reference (\label{$prefix:$chap_...)

Rename the appendices from 'annex_' to 'apdx_' to conform with the prefix used in labels (apdx:...)
Suppress the letter numbering

Location:
NEMO/trunk/doc/latex
Files:
19 edited
7 moved

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/doc/latex/NEMO/main/appendices.tex

    r11330 r11543  
    11 
    2 \subfile{../subfiles/annex_A}             %% Generalised vertical coordinate 
    3 \subfile{../subfiles/annex_B}             %% Diffusive operator 
    4 \subfile{../subfiles/annex_C}             %% Discrete invariants of the eqs. 
    5 \subfile{../subfiles/annex_iso}            %% Isoneutral diffusion using triads 
    6 \subfile{../subfiles/annex_DOMAINcfg}     %% Brief notes on DOMAINcfg 
     2\subfile{../subfiles/apdx_s_coord}      %% Generalised vertical coordinate 
     3\subfile{../subfiles/apdx_diff_opers}   %% Diffusive operators 
     4\subfile{../subfiles/apdx_invariants}   %% Discrete invariants of the eqs. 
     5\subfile{../subfiles/apdx_triads}       %% Isoneutral diffusion using triads 
     6\subfile{../subfiles/apdx_DOMAINcfg}    %% Brief notes on DOMAINcfg 
    77 
    88%% Not included 
     
    1010%\subfile{../subfiles/chap_DIU} 
    1111%\subfile{../subfiles/chap_conservation} 
    12 %\subfile{../subfiles/annex_E}            %% Notes on some on going staff 
    13  
     12%\subfile{../subfiles/apdx_algos}   %% Notes on some on going staff 
  • NEMO/trunk/doc/latex/NEMO/main/chapters.tex

    r11522 r11543  
    1313\subfile{../subfiles/chap_STO}            %% Stochastic param. 
    1414\subfile{../subfiles/chap_misc}           %% Miscellaneous topics 
    15 \subfile{../subfiles/chap_CONFIG}         %% Predefined configurations 
     15\subfile{../subfiles/chap_cfgs}           %% Predefined configurations 
    1616 
    1717%% Not included 
  • NEMO/trunk/doc/latex/NEMO/main/introduction.tex

    r11522 r11543  
    5656 
    5757\begin{description} 
    58 \item [\nameref{chap:PE}] presents the equations and their assumptions, the vertical coordinates used, 
     58\item [\nameref{chap:MB}] presents the equations and their assumptions, the vertical coordinates used, 
    5959and the subgrid scale physics. 
    6060The equations are written in a curvilinear coordinate system, with a choice of vertical coordinates 
     
    6363Dimensional units in the meter, kilogram, second (MKS) international system are used throughout. 
    6464The following chapters deal with the discrete equations. 
    65 \item [\nameref{chap:STP}] presents the model time stepping environment. 
     65\item [\nameref{chap:TD}] presents the model time stepping environment. 
    6666it is a three level scheme in which the tendency terms of the equations are evaluated either 
    6767centered in time, or forward, or backward depending of the nature of the term. 
     
    123123\item [\nameref{chap:ASM}] describes how increments produced by 
    124124data \textbf{A}s\textbf{S}i\textbf{M}ilation may be applied to the model equations. 
     125\item [\nameref{chap:STO}] 
    125126\item [\nameref{chap:MISC}] (including solvers) 
    126 \item [\nameref{chap:CFG}] provides finally a brief introduction to 
     127\item [\nameref{chap:CFGS}] provides finally a brief introduction to 
    127128the pre-defined model configurations 
    128129(water column model \texttt{C1D}, ORCA and GYRE families of configurations). 
     
    133134 
    134135\begin{description} 
    135 \item [\nameref{apdx:s_coord}] 
    136 \item [\nameref{apdx:diff_oper}] 
    137 \item [\nameref{apdx:invariants}] 
    138 \item [\nameref{apdx:triads}] 
    139 \item [\nameref{apdx:DOMAINcfg}] 
    140 \item [\nameref{apdx:coding}] 
     136\item [\nameref{apdx:SCOORD}] 
     137\item [\nameref{apdx:DIFFOPERS}] 
     138\item [\nameref{apdx:INVARIANTS}] 
     139\item [\nameref{apdx:TRIADS}] 
     140\item [\nameref{apdx:DOMCFG}] 
     141\item [\nameref{apdx:CODING}] 
    141142\end{description} 
  • NEMO/trunk/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex

    r11529 r11543  
    66% ================================================================ 
    77\chapter{A brief guide to the DOMAINcfg tool} 
    8 \label{apdx:DOMAINcfg} 
     8\label{apdx:DOMCFG} 
    99 
    1010\chaptertoc 
     
    121121The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and 
    122122$gdepw_0$ for $t$- and $w$-points, respectively. See \autoref{sec:DOMCFG_sco} for the 
    123 S-coordinate options.  As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of 
     123S-coordinate options.  As indicated on \autoref{fig:DOM_index_vert} \jp{jpk} is the number of 
    124124$w$-levels.  $gdepw_0(1)$ is the ocean surface.  There are at most \jp{jpk}-1 $t$-points 
    125125inside the ocean, the additional $t$-point at $jk = jpk$ is below the sea floor and is not 
     
    421421The depth field $h$ is not necessary the ocean depth, 
    422422since a mixed step-like and bottom-following representation of the topography can be used 
    423 (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). 
     423(\autoref{fig:DOM_z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:DOM_z_zps_s_sps}). 
    424424The namelist parameter \np{rn\_rmax} determines the slope at which 
    425425the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 
     
    436436\[ 
    437437  z = s_{min} + C (s) (H - s_{min}) 
    438   % \label{eq:SH94_1} 
     438  % \label{eq:DOMCFG_SH94_1} 
    439439\] 
    440440 
     
    458458         + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)} 
    459459                        {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)} 
    460  \label{eq:SH94_2} 
     460 \label{eq:DOMCFG_SH94_2} 
    461461\] 
    462462 
     
    466466    \includegraphics[width=\textwidth]{Fig_sco_function} 
    467467    \caption{ 
    468       \protect\label{fig:sco_function} 
     468      \protect\label{fig:DOMCFG_sco_function} 
    469469      Examples of the stretching function applied to a seamount; 
    470470      from left to right: surface, surface and bottom, and bottom intensified resolutions 
     
    478478bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. 
    479479$b$ has been designed to allow surface and/or bottom increase of the vertical resolution 
    480 (\autoref{fig:sco_function}). 
     480(\autoref{fig:DOMCFG_sco_function}). 
    481481 
    482482Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in 
     
    486486\begin{equation} 
    487487  z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 
    488   % \label{eq:z} 
     488  % \label{eq:DOMCFG_z} 
    489489\end{equation} 
    490490 
     
    524524    For clarity every third coordinate surface is shown. 
    525525  } 
    526   \label{fig:fig_compare_coordinates_surface} 
     526  \label{fig:DOMCFG_fig_compare_coordinates_surface} 
    527527\end{figure} 
    528528 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • NEMO/trunk/doc/latex/NEMO/subfiles/apdx_algos.tex

    r11529 r11543  
    66% ================================================================ 
    77\chapter{Note on some algorithms} 
    8 \label{apdx:E} 
     8\label{apdx:ALGOS} 
    99 
    1010\chaptertoc 
     
    1212\newpage 
    1313 
    14 This appendix some on going consideration on algorithms used or planned to be used in \NEMO.  
     14This appendix some on going consideration on algorithms used or planned to be used in \NEMO. 
    1515 
    1616% ------------------------------------------------------------------------------------------------------------- 
    17 %        UBS scheme   
     17%        UBS scheme 
    1818% ------------------------------------------------------------------------------------------------------------- 
    1919\section{Upstream Biased Scheme (UBS) (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})} 
     
    4545a constant i-grid spacing ($\Delta i=1$). 
    4646 
    47 Alternative choice: introduce the scale factors:   
     47Alternative choice: introduce the scale factors: 
    4848$\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta_i \left[ \frac{e_{2u} e_{3u} }{e_{1u} }\delta_{i+1/2}[\tau] \right]$. 
    4949 
     
    5454It is not a \emph{positive} scheme meaning false extrema are permitted but 
    5555the amplitude of such are significantly reduced over the centred second order method. 
    56 Nevertheless it is not recommended to apply it to a passive tracer that requires positivity.  
     56Nevertheless it is not recommended to apply it to a passive tracer that requires positivity. 
    5757 
    5858The intrinsic diffusion of UBS makes its use risky in the vertical direction where 
     
    6161\np{ln\_traadv\_ubs}\forcode{ = .true.}. 
    6262 
    63 For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds to 
     63For stability reasons, in \autoref{eq:TRA_adv_ubs}, the first term which corresponds to 
    6464a second order centred scheme is evaluated using the \textit{now} velocity (centred in time) while 
    6565the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity 
     
    6767This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. 
    6868UBS and QUICK schemes only differ by one coefficient. 
    69 Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
     69Substituting 1/6 with 1/8 in (\autoref{eq:TRA_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
    7070This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. 
    7171Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme. 
     
    7575Computer time can be saved by using a time-splitting technique on vertical advection. 
    7676This possibility have been implemented and validated in ORCA05-L301. 
    77 It is not currently offered in the current reference version.  
     77It is not currently offered in the current reference version. 
    7878 
    7979NB 2: In a forthcoming release four options will be proposed for the vertical component used in the UBS scheme. 
     
    8383The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme. 
    8484 
    85 NB 3: It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows: 
     85NB 3: It is straight forward to rewrite \autoref{eq:TRA_adv_ubs} as follows: 
    8686\begin{equation} 
    8787  \label{eq:tra_adv_ubs2} 
     
    9393  \right. 
    9494\end{equation} 
    95 or equivalently  
     95or equivalently 
    9696\begin{equation} 
    9797  \label{eq:tra_adv_ubs2} 
     
    102102  \end{split} 
    103103\end{equation} 
    104 \autoref{eq:tra_adv_ubs2} has several advantages. 
     104\autoref{eq:TRA_adv_ubs2} has several advantages. 
    105105First it clearly evidences that the UBS scheme is based on the fourth order scheme to which 
    106106is added an upstream biased diffusive term. 
    107107Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, 
    108 not only the $2^{th}$ order part as stated above using \autoref{eq:tra_adv_ubs}. 
     108not only the $2^{th}$ order part as stated above using \autoref{eq:TRA_adv_ubs}. 
    109109Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient which 
    110110is simply proportional to the velocity. 
     
    134134  \end{split} 
    135135\end{equation} 
    136 with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$,  
     136with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$, 
    137137\ie\ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$ 
    138138it comes: 
     
    189189 
    190190% ------------------------------------------------------------------------------------------------------------- 
    191 %        Leap-Frog energetic   
     191%        Leap-Frog energetic 
    192192% ------------------------------------------------------------------------------------------------------------- 
    193193\section{Leapfrog energetic} 
     
    214214  \equiv \frac{1}{\rdt} \overline{ \delta_{t+\rdt/2}[q]}^{\,t} 
    215215  =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt} 
    216 \]  
     216\] 
    217217Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, 
    218218not $2\rdt$ as it can be found sometimes in literature. 
     
    226226\] 
    227227is satisfied in discrete form. 
    228 Indeed,  
     228Indeed, 
    229229\[ 
    230230  \begin{split} 
     
    240240\] 
    241241NB here pb of boundary condition when applying the adjoint! 
    242 In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition  
     242In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition 
    243243(equivalently of the boundary value of the integration by part). 
    244244In time this boundary condition is not physical and \textbf{add something here!!!} 
    245245 
    246246% ================================================================ 
    247 % Iso-neutral diffusion :  
     247% Iso-neutral diffusion : 
    248248% ================================================================ 
    249249 
     
    251251 
    252252% ================================================================ 
    253 % Griffies' iso-neutral diffusion operator :  
     253% Griffies' iso-neutral diffusion operator : 
    254254% ================================================================ 
    255255\subsection{Griffies iso-neutral diffusion operator} 
     
    258258but is formulated within the \NEMO\ framework 
    259259(\ie\ using scale factors rather than grid-size and having a position of $T$-points that 
    260 is not necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}). 
    261  
    262 In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO, 
     260is not necessary in the middle of vertical velocity points, see \autoref{fig:DOM_zgr_e3}). 
     261 
     262In the formulation \autoref{eq:TRA_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO, 
    263263the off-diagonal terms of the small angle diffusion tensor contain several double spatial averages of a gradient, 
    264264for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$. 
     
    318318%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    319319 
    320 The four iso-neutral fluxes associated with the triads are defined at $T$-point.  
     320The four iso-neutral fluxes associated with the triads are defined at $T$-point. 
    321321They take the following expression: 
    322322\begin{flalign*} 
     
    332332 
    333333The resulting iso-neutral fluxes at $u$- and $w$-points are then given by 
    334 the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}): 
     334the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:TRIADS_ISO_triad}): 
    335335\begin{flalign} 
    336336  \label{eq:iso_flux} 
     
    369369    + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]   \right\} 
    370370\end{equation} 
    371 where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.  
     371where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. 
    372372 
    373373This expression of the iso-neutral diffusion has been chosen in order to satisfy the following six properties: 
     
    448448 
    449449% ================================================================ 
    450 % Skew flux formulation for Eddy Induced Velocity :  
     450% Skew flux formulation for Eddy Induced Velocity : 
    451451% ================================================================ 
    452452\subsection{Eddy induced velocity and skew flux formulation} 
     
    457457the formulation of which depends on the slopes of iso-neutral surfaces. 
    458458Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, 
    459 \ie\ \autoref{eq:ldfslp_geo} is used in $z$-coordinate, 
    460 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.  
    461  
    462 The eddy induced velocity is given by:  
     459\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinate, 
     460and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $z^*$ or $s$-coordinates. 
     461 
     462The eddy induced velocity is given by: 
    463463\begin{equation} 
    464464  \label{eq:eiv_v} 
     
    484484(see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme. 
    485485This is particularly useful for passive tracers where 
    486 \emph{positivity} of the advection scheme is of paramount importance.  
    487 % give here the expression using the triads. It is different from the one given in \autoref{eq:ldfeiv} 
     486\emph{positivity} of the advection scheme is of paramount importance. 
     487% give here the expression using the triads. It is different from the one given in \autoref{eq:LDF_eiv} 
    488488% see just below a copy of this equation: 
    489489%\begin{equation} \label{eq:ldfeiv} 
     
    593593  \right) 
    594594\end{equation} 
    595 Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells.  
     595Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells. 
    596596In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to 
    597 $\mathbb{R}$ for the discret form to be exact.  
     597$\mathbb{R}$ for the discret form to be exact. 
    598598 
    599599Such a choice of discretisation is consistent with the iso-neutral operator as 
     
    604604$\ $\newpage      %force an empty line 
    605605% ================================================================ 
    606 % Discrete Invariants of the iso-neutral diffrusion  
     606% Discrete Invariants of the iso-neutral diffrusion 
    607607% ================================================================ 
    608608\subsection{Discrete invariants of the iso-neutral diffrusion} 
    609609\label{subsec:Gf_operator} 
    610610 
    611 Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane.  
     611Demonstration of the decrease of the tracer variance in the (\textbf{i},\textbf{j}) plane. 
    612612 
    613613This part will be moved in an Appendix. 
     
    617617  \int_D  D_l^T \; T \;dv   \leq 0 
    618618\] 
    619 The discrete form of its left hand side is obtained using \autoref{eq:iso_flux} 
     619The discrete form of its left hand side is obtained using \autoref{eq:TRIADS_iso_flux} 
    620620 
    621621\begin{align*} 
     
    740740        \right\} 
    741741        \quad   \leq 0 
    742 \end{align*}  
     742\end{align*} 
    743743The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities. 
    744744 
     
    764764                             % 
    765765                           &\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\} 
    766 \end{align*}  
     766\end{align*} 
    767767This means that the iso-neutral operator is self-adjoint. 
    768768There is no need to develop a specific to obtain it. 
     
    776776\label{subsec:eiv_skew} 
    777777 
    778 Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.  
     778Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. 
    779779 
    780780This have to be moved in an Appendix. 
     
    830830    &{\ \ \;_i^k  \mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}] 
    831831    &\Bigr\}  \\ 
    832   \end{matrix}    
     832  \end{matrix} 
    833833\end{align*} 
    834834The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs, 
    835 they cancel out.  
     835they cancel out. 
    836836Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. 
    837837The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and 
  • NEMO/trunk/doc/latex/NEMO/subfiles/apdx_diff_opers.tex

    r11529 r11543  
    55% Chapter Appendix B : Diffusive Operators 
    66% ================================================================ 
    7 \chapter{Appendix B : Diffusive Operators} 
    8 \label{apdx:B} 
     7\chapter{Diffusive Operators} 
     8\label{apdx:DIFFOPERS} 
    99 
    1010\chaptertoc 
     
    1616% ================================================================ 
    1717\section{Horizontal/Vertical $2^{nd}$ order tracer diffusive operators} 
    18 \label{sec:B_1} 
     18\label{sec:DIFFOPERS_1} 
    1919 
    2020\subsubsection*{In z-coordinates} 
     
    2222In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator is given by: 
    2323\begin{align} 
    24   \label{apdx:B1} 
     24  \label{eq:DIFFOPERS_1} 
    2525  &D^T = \frac{1}{e_1 \, e_2}      \left[ 
    2626    \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. 
     
    3232\subsubsection*{In generalized vertical coordinates} 
    3333 
    34 In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by \autoref{apdx:A_s_slope} and 
     34In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by \autoref{eq:SCOORD_s_slope} and 
    3535the vertical/horizontal ratio of diffusion coefficient by $\epsilon = A^{vT} / A^{lT}$. 
    3636The diffusion operator is given by: 
    3737 
    3838\begin{equation} 
    39   \label{apdx:B2} 
     39  \label{eq:DIFFOPERS_2} 
    4040  D^T = \left. \nabla \right|_s \cdot 
    4141  \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\ 
     
    5454  \begin{array}{*{20}l} 
    5555    D^T= \frac{1}{e_1\,e_2\,e_3 } & \left\{ \quad \quad \frac{\partial }{\partial i}  \left. \left[  e_2\,e_3 \, A^{lT} 
    56                                \left( \  \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s  
     56                               \left( \  \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s 
    5757                                       -\frac{\sigma_1 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\ 
    5858        &  \quad \  +   \            \left.   \frac{\partial }{\partial j}  \left. \left[  e_1\,e_3 \, A^{lT} 
    59                                \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s  
     59                               \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s 
    6060                                       -\frac{\sigma_2 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\ 
    61         &  \quad \  +   \           \left.  e_1\,e_2\, \frac{\partial }{\partial s}  \left[ A^{lT} \; \left(  
    62                      -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s  
    63                      -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s  
     61        &  \quad \  +   \           \left.  e_1\,e_2\, \frac{\partial }{\partial s}  \left[ A^{lT} \; \left( 
     62                     -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s 
     63                     -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s 
    6464                          +\left( \varepsilon +\sigma_1^2+\sigma_2 ^2 \right) \; \frac{1}{e_3 } \; \frac{\partial T}{\partial s} \right) \; \right] \;  \right\} . 
    6565  \end{array} 
     
    6767\end{align*} 
    6868 
    69 \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption. 
     69\autoref{eq:DIFFOPERS_2} is obtained from \autoref{eq:DIFFOPERS_1} without any additional assumption. 
    7070Indeed, for the special case $k=z$ and thus $e_3 =1$, 
    71 we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and 
    72 use \autoref{apdx:A_s_slope} and \autoref{apdx:A_s_chain_rule}. 
    73 Since no cross horizontal derivative $\partial _i \partial _j $ appears in \autoref{apdx:B1}, 
     71we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:SCOORD} and 
     72use \autoref{eq:SCOORD_s_slope} and \autoref{eq:SCOORD_s_chain_rule}. 
     73Since no cross horizontal derivative $\partial _i \partial _j $ appears in \autoref{eq:DIFFOPERS_1}, 
    7474the ($i$,$z$) and ($j$,$z$) planes are independent. 
    7575The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$) transformation without 
     
    160160% ================================================================ 
    161161\section{Iso/Diapycnal $2^{nd}$ order tracer diffusive operators} 
    162 \label{sec:B_2} 
     162\label{sec:DIFFOPERS_2} 
    163163 
    164164\subsubsection*{In z-coordinates} 
     
    170170 
    171171\begin{equation} 
    172   \label{apdx:B3} 
     172  \label{eq:DIFFOPERS_3} 
    173173  \textbf {A}_{\textbf I} = \frac{A^{lT}}{\left( {1+a_1 ^2+a_2 ^2} \right)} 
    174174  \left[ {{ 
     
    193193 
    194194In practice, $\epsilon$ is small and isopycnal slopes are generally less than $10^{-2}$ in the ocean, 
    195 so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}. Keeping leading order terms\footnote{Apart from the (1,0)  
     195so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}. Keeping leading order terms\footnote{Apart from the (1,0) 
    196196and (0,1) elements which are set to zero. See \citet{griffies_bk04}, section 14.1.4.1 for a discussion of this point.}: 
    197197\begin{subequations} 
    198   \label{apdx:B4} 
     198  \label{eq:DIFFOPERS_4} 
    199199  \begin{equation} 
    200     \label{apdx:B4a} 
     200    \label{eq:DIFFOPERS_4a} 
    201201    {\textbf{A}_{\textbf{I}}} \approx A^{lT}\;\Re\;\text{where} \;\Re = 
    202202    \left[ {{ 
     
    210210  and the iso/dianeutral diffusive operator in $z$-coordinates is then 
    211211  \begin{equation} 
    212     \label{apdx:B4b} 
     212    \label{eq:DIFFOPERS_4b} 
    213213    D^T = \left. \nabla \right|_z \cdot 
    214214    \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_z T  \right]. \\ 
     
    216216\end{subequations} 
    217217 
    218 Physically, the full tensor \autoref{apdx:B3} represents strong isoneutral diffusion on a plane parallel to 
     218Physically, the full tensor \autoref{eq:DIFFOPERS_3} represents strong isoneutral diffusion on a plane parallel to 
    219219the isoneutral surface and weak dianeutral diffusion perpendicular to this plane. 
    220220However, 
    221 the approximate `weak-slope' tensor \autoref{apdx:B4a} represents strong diffusion along the isoneutral surface, 
     221the approximate `weak-slope' tensor \autoref{eq:DIFFOPERS_4a} represents strong diffusion along the isoneutral surface, 
    222222with weak \emph{vertical} diffusion -- the principal axes of the tensor are no longer orthogonal. 
    223223This simplification also decouples the ($i$,$z$) and ($j$,$z$) planes of the tensor. 
    224 The weak-slope operator therefore takes the same form, \autoref{apdx:B4}, as \autoref{apdx:B2}, 
     224The weak-slope operator therefore takes the same form, \autoref{eq:DIFFOPERS_4}, as \autoref{eq:DIFFOPERS_2}, 
    225225the diffusion operator for geopotential diffusion written in non-orthogonal $i,j,s$-coordinates. 
    226226Written out explicitly, 
    227227 
    228228\begin{multline} 
    229   \label{apdx:B_ldfiso} 
     229  \label{eq:DIFFOPERS_ldfiso} 
    230230  D^T=\frac{1}{e_1 e_2 }\left\{ 
    231231    {\;\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]} 
     
    234234\end{multline} 
    235235 
    236 The isopycnal diffusion operator \autoref{apdx:B4}, 
    237 \autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square. 
    238 As \autoref{apdx:B4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero  
     236The isopycnal diffusion operator \autoref{eq:DIFFOPERS_4}, 
     237\autoref{eq:DIFFOPERS_ldfiso} conserves tracer quantity and dissipates its square. 
     238As \autoref{eq:DIFFOPERS_4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero 
    239239(as it is when $A_h$ is zero at the boundary). Let us demonstrate the second one: 
    240240\[ 
     
    256256             j}-a_2 \frac{\partial T}{\partial k}} \right)^2} 
    257257             +\varepsilon \left(\frac{\partial T}{\partial k}\right) ^2\right]      \\ 
    258            & \geq 0 .  
     258           & \geq 0 . 
    259259  \end{array} 
    260260             } 
     
    265265\subsubsection*{In generalized vertical coordinates} 
    266266 
    267 Because the weak-slope operator \autoref{apdx:B4}, 
    268 \autoref{apdx:B_ldfiso} is decoupled in the ($i$,$z$) and ($j$,$z$) planes, 
     267Because the weak-slope operator \autoref{eq:DIFFOPERS_4}, 
     268\autoref{eq:DIFFOPERS_ldfiso} is decoupled in the ($i$,$z$) and ($j$,$z$) planes, 
    269269it may be transformed into generalized $s$-coordinates in the same way as 
    270 \autoref{sec:B_1} was transformed into \autoref{sec:B_2}. 
     270\autoref{sec:DIFFOPERS_1} was transformed into \autoref{sec:DIFFOPERS_2}. 
    271271The resulting operator then takes the simple form 
    272272 
    273273\begin{equation} 
    274   \label{apdx:B_ldfiso_s} 
     274  \label{eq:DIFFOPERS_ldfiso_s} 
    275275  D^T = \left. \nabla \right|_s \cdot 
    276276  \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\ 
     
    295295\] 
    296296 
    297 To prove \autoref{apdx:B_ldfiso_s} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious. 
    298 An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that 
     297To prove \autoref{eq:DIFFOPERS_ldfiso_s} by direct re-expression of \autoref{eq:DIFFOPERS_ldfiso} is straightforward, but laborious. 
     298An easier way is first to note (by reversing the derivation of \autoref{sec:DIFFOPERS_2} from \autoref{sec:DIFFOPERS_1} ) that 
    299299the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as 
    300300 
    301301\begin{equation} 
    302   \label{apdx:B5} 
     302  \label{eq:DIFFOPERS_5} 
    303303  D^T = \left. \nabla \right|_\rho \cdot 
    304304  \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_\rho T  \right] \\ 
     
    312312\end{equation} 
    313313Then direct transformation from $i,j,\rho$-coordinates to $i,j,s$-coordinates gives 
    314 \autoref{apdx:B_ldfiso_s} immediately. 
     314\autoref{eq:DIFFOPERS_ldfiso_s} immediately. 
    315315 
    316316Note that the weak-slope approximation is only made in transforming from 
     
    318318The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the $s$-surfaces, 
    319319in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates in 
    320 \autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces. 
     320\autoref{sec:DIFFOPERS_1} onto $s$-coordinates is exact, however steep the $s$-surfaces. 
    321321 
    322322 
     
    325325% ================================================================ 
    326326\section{Lateral/Vertical momentum diffusive operators} 
    327 \label{sec:B_3} 
     327\label{sec:DIFFOPERS_3} 
    328328 
    329329The second order momentum diffusion operator (Laplacian) in $z$-coordinates is found by 
    330 applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector, 
     330applying \autoref{eq:MB_lap_vector}, the expression for the Laplacian of a vector, 
    331331to the horizontal velocity vector: 
    332332\begin{align*} 
     
    371371  }} \right) 
    372372\end{align*} 
    373 Using \autoref{eq:PE_div}, the definition of the horizontal divergence, 
     373Using \autoref{eq:MB_div}, the definition of the horizontal divergence, 
    374374the third component of the second vector is obviously zero and thus : 
    375375\[ 
    376   \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) .  
     376  \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) . 
    377377\] 
    378378 
    379379Note that this operator ensures a full separation between 
    380 the vorticity and horizontal divergence fields (see \autoref{apdx:C}). 
     380the vorticity and horizontal divergence fields (see \autoref{apdx:INVARIANTS}). 
    381381It is only equal to a Laplacian applied to each component in Cartesian coordinates, not on the sphere. 
    382382 
     
    384384the $z$-coordinate therefore takes the following form: 
    385385\begin{equation} 
    386   \label{apdx:B_Lap_U} 
     386  \label{eq:DIFFOPERS_Lap_U} 
    387387  { 
    388388    \textbf{D}}^{\textbf{U}} = 
     
    404404\end{align*} 
    405405 
    406 Note Bene: introducing a rotation in \autoref{apdx:B_Lap_U} does not lead to 
     406Note Bene: introducing a rotation in \autoref{eq:DIFFOPERS_Lap_U} does not lead to 
    407407a useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate. 
    408408Similarly, we did not found an expression of practical use for 
    409409the geopotential horizontal/vertical Laplacian operator in the $s$-coordinate. 
    410 Generally, \autoref{apdx:B_Lap_U} is used in both $z$- and $s$-coordinate systems, 
     410Generally, \autoref{eq:DIFFOPERS_Lap_U} is used in both $z$- and $s$-coordinate systems, 
    411411that is a Laplacian diffusion is applied on momentum along the coordinate directions. 
    412412 
  • NEMO/trunk/doc/latex/NEMO/subfiles/apdx_invariants.tex

    r11529 r11543  
    66% ================================================================ 
    77\chapter{Discrete Invariants of the Equations} 
    8 \label{apdx:C} 
     8\label{apdx:INVARIANTS} 
    99 
    1010\chaptertoc 
     
    2121% ================================================================ 
    2222\section{Introduction / Notations} 
    23 \label{sec:C.0} 
     23\label{sec:INVARIANTS_0} 
    2424 
    2525Notation used in this appendix in the demonstations: 
     
    7272that is in a more compact form : 
    7373\begin{flalign} 
    74   \label{eq:Q2_flux} 
     74  \label{eq:INVARIANTS_Q2_flux} 
    7575  \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right) 
    7676  =&                   \int_D { \frac{Q}{e_3}  \partial_t \left( e_3 \, Q \right) dv } 
     
    8787that is in a more compact form: 
    8888\begin{flalign} 
    89   \label{eq:Q2_vect} 
     89  \label{eq:INVARIANTS_Q2_vect} 
    9090  \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right) 
    9191  =& \int_D {         Q   \,\partial_t Q  \;dv } 
     
    9797% ================================================================ 
    9898\section{Continuous conservation} 
    99 \label{sec:C.1} 
     99\label{sec:INVARIANTS_1} 
    100100 
    101101The discretization of pimitive equation in $s$-coordinate (\ie\ time and space varying vertical coordinate) 
     
    105105The total energy (\ie\ kinetic plus potential energies) is conserved: 
    106106\begin{flalign} 
    107   \label{eq:Tot_Energy} 
     107  \label{eq:INVARIANTS_Tot_Energy} 
    108108  \partial_t \left( \int_D \left( \frac{1}{2} {\textbf{U}_h}^2 +  \rho \, g \, z \right) \;dv \right)  = & 0 
    109109\end{flalign} 
     
    114114The transformation for the advection term depends on whether the vector invariant form or 
    115115the flux form is used for the momentum equation. 
    116 Using \autoref{eq:Q2_vect} and introducing \autoref{apdx:A_dyn_vect} in 
    117 \autoref{eq:Tot_Energy} for the former form and 
    118 using \autoref{eq:Q2_flux} and introducing \autoref{apdx:A_dyn_flux} in 
    119 \autoref{eq:Tot_Energy} for the latter form leads to: 
    120  
    121 % \label{eq:E_tot} 
     116Using \autoref{eq:INVARIANTS_Q2_vect} and introducing \autoref{eq:SCOORD_dyn_vect} in 
     117\autoref{eq:INVARIANTS_Tot_Energy} for the former form and 
     118using \autoref{eq:INVARIANTS_Q2_flux} and introducing \autoref{eq:SCOORD_dyn_flux} in 
     119\autoref{eq:INVARIANTS_Tot_Energy} for the latter form leads to: 
     120 
     121% \label{eq:INVARIANTS_E_tot} 
    122122advection term (vector invariant form): 
    123123\[ 
    124   % \label{eq:E_tot_vect_vor_1} 
     124  % \label{eq:INVARIANTS_E_tot_vect_vor_1} 
    125125  \int\limits_D  \zeta \; \left( \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0   \\ 
    126126\] 
    127127% 
    128128\[ 
    129   % \label{eq:E_tot_vect_adv_1} 
     129  % \label{eq:INVARIANTS_E_tot_vect_adv_1} 
    130130  \int\limits_D  \textbf{U}_h \cdot \nabla_h \left( \frac{{\textbf{U}_h}^2}{2} \right)     dv 
    131131  + \int\limits_D  \textbf{U}_h \cdot \nabla_z \textbf{U}_h  \;dv 
     
    134134advection term (flux form): 
    135135\[ 
    136   % \label{eq:E_tot_flux_metric} 
     136  % \label{eq:INVARIANTS_E_tot_flux_metric} 
    137137  \int\limits_D  \frac{1} {e_1 e_2 } \left( v \,\partial_i e_2 - u \,\partial_j e_1  \right)\; 
    138138  \left(  \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0 
    139139\] 
    140140\[ 
    141   % \label{eq:E_tot_flux_adv} 
     141  % \label{eq:INVARIANTS_E_tot_flux_adv} 
    142142  \int\limits_D \textbf{U}_h \cdot     \left(                 {{ 
    143143        \begin{array} {*{20}c} 
     
    150150coriolis term 
    151151\[ 
    152   % \label{eq:E_tot_cor} 
     152  % \label{eq:INVARIANTS_E_tot_cor} 
    153153  \int\limits_D  f   \; \left( \textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv   = 0 
    154154\] 
    155155pressure gradient: 
    156156\[ 
    157   % \label{eq:E_tot_pg_1} 
     157  % \label{eq:INVARIANTS_E_tot_pg_1} 
    158158  - \int\limits_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv 
    159159  = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv 
     
    173173 
    174174Vector invariant form: 
    175 % \label{eq:E_tot_vect} 
    176 \[ 
    177   % \label{eq:E_tot_vect_vor_2} 
     175% \label{eq:INVARIANTS_E_tot_vect} 
     176\[ 
     177  % \label{eq:INVARIANTS_E_tot_vect_vor_2} 
    178178  \int\limits_D   \textbf{U}_h \cdot \text{VOR} \;dv   = 0 
    179179\] 
    180180\[ 
    181   % \label{eq:E_tot_vect_adv_2} 
     181  % \label{eq:INVARIANTS_E_tot_vect_adv_2} 
    182182  \int\limits_D  \textbf{U}_h \cdot \text{KEG}  \;dv 
    183183  + \int\limits_D  \textbf{U}_h \cdot \text{ZAD}  \;dv 
     
    185185\] 
    186186\[ 
    187   % \label{eq:E_tot_pg_2} 
     187  % \label{eq:INVARIANTS_E_tot_pg_2} 
    188188  - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv 
    189189  = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv 
     
    193193Flux form: 
    194194\begin{subequations} 
    195   \label{eq:E_tot_flux} 
     195  \label{eq:INVARIANTS_E_tot_flux} 
    196196  \[ 
    197     % \label{eq:E_tot_flux_metric_2} 
     197    % \label{eq:INVARIANTS_E_tot_flux_metric_2} 
    198198    \int\limits_D  \textbf{U}_h \cdot \text {COR} \;  dv   = 0 
    199199  \] 
    200200  \[ 
    201     % \label{eq:E_tot_flux_adv_2} 
     201    % \label{eq:INVARIANTS_E_tot_flux_adv_2} 
    202202    \int\limits_D \textbf{U}_h \cdot \text{ADV}   \;dv 
    203203    +   \frac{1}{2} \int\limits_D {  {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t e_3  \;dv } =\;0 
    204204  \] 
    205205  \begin{equation} 
    206     \label{eq:E_tot_pg_3} 
     206    \label{eq:INVARIANTS_E_tot_pg_3} 
    207207    - \int\limits_D  \textbf{U}_h \cdot (\text{HPG}+ \text{SPG}) \;dv 
    208208    = - \int\limits_D \nabla \cdot \left( \rho \,\textbf {U} \right)\;g\;z\;\;dv 
     
    211211\end{subequations} 
    212212 
    213 \autoref{eq:E_tot_pg_3} is the balance between the conversion KE to PE and PE to KE. 
    214 Indeed the left hand side of \autoref{eq:E_tot_pg_3} can be transformed as follows: 
     213\autoref{eq:INVARIANTS_E_tot_pg_3} is the balance between the conversion KE to PE and PE to KE. 
     214Indeed the left hand side of \autoref{eq:INVARIANTS_E_tot_pg_3} can be transformed as follows: 
    215215\begin{flalign*} 
    216216  \partial_t  \left( \int\limits_D { \rho \, g \, z  \;dv} \right) 
     
    225225\end{flalign*} 
    226226where the last equality is obtained by noting that the brackets is exactly the expression of $w$, 
    227 the vertical velocity referenced to the fixe $z$-coordinate system (see \autoref{apdx:A_w_s}). 
    228  
    229 The left hand side of \autoref{eq:E_tot_pg_3} can be transformed as follows: 
     227the vertical velocity referenced to the fixe $z$-coordinate system (see \autoref{eq:SCOORD_w_s}). 
     228 
     229The left hand side of \autoref{eq:INVARIANTS_E_tot_pg_3} can be transformed as follows: 
    230230\begin{flalign*} 
    231231  - \int\limits_D  \left. \nabla p \right|_z & \cdot \textbf{U}_h \;dv 
     
    326326% ================================================================ 
    327327\section{Discrete total energy conservation: vector invariant form} 
    328 \label{sec:C.2} 
     328\label{sec:INVARIANTS_2} 
    329329 
    330330% ------------------------------------------------------------------------------------------------------------- 
     
    332332% ------------------------------------------------------------------------------------------------------------- 
    333333\subsection{Total energy conservation} 
    334 \label{subsec:C_KE+PE_vect} 
    335  
    336 The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by: 
     334\label{subsec:INVARIANTS_KE+PE_vect} 
     335 
     336The discrete form of the total energy conservation, \autoref{eq:INVARIANTS_Tot_Energy}, is given by: 
    337337\begin{flalign*} 
    338338  \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 
     
    340340which in vector invariant forms, it leads to: 
    341341\begin{equation} 
    342   \label{eq:KE+PE_vect_discrete} 
     342  \label{eq:INVARIANTS_KE+PE_vect_discrete} 
    343343  \begin{split} 
    344344    \sum\limits_{i,j,k} \biggl\{   u\,                        \partial_t u         \;b_u 
     
    352352 
    353353Substituting the discrete expression of the time derivative of the velocity either in vector invariant, 
    354 leads to the discrete equivalent of the four equations \autoref{eq:E_tot_flux}. 
     354leads to the discrete equivalent of the four equations \autoref{eq:INVARIANTS_E_tot_flux}. 
    355355 
    356356% ------------------------------------------------------------------------------------------------------------- 
     
    358358% ------------------------------------------------------------------------------------------------------------- 
    359359\subsection{Vorticity term (coriolis + vorticity part of the advection)} 
    360 \label{subsec:C_vor} 
     360\label{subsec:INVARIANTS_vor} 
    361361 
    362362Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$), 
     
    367367% ------------------------------------------------------------------------------------------------------------- 
    368368\subsubsection{Vorticity term with ENE scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 
    369 \label{subsec:C_vorENE} 
     369\label{subsec:INVARIANTS_vorENE} 
    370370 
    371371For the ENE scheme, the two components of the vorticity term are given by: 
     
    407407% ------------------------------------------------------------------------------------------------------------- 
    408408\subsubsection{Vorticity term with EEN scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 
    409 \label{subsec:C_vorEEN_vect} 
     409\label{subsec:INVARIANTS_vorEEN_vect} 
    410410 
    411411With the EEN scheme, the vorticity terms are represented as: 
    412412\begin{equation} 
    413   \tag{\ref{eq:dynvor_een}} 
     413  \label{eq:INVARIANTS_dynvor_een} 
    414414  \left\{ { 
    415415      \begin{aligned} 
     
    424424and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: 
    425425\begin{equation} 
    426   \tag{\ref{eq:Q_triads}} 
     426  \label{eq:INVARIANTS_Q_triads} 
    427427  _i^j \mathbb{Q}^{i_p}_{j_p} 
    428428  = \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) 
     
    479479% ------------------------------------------------------------------------------------------------------------- 
    480480\subsubsection{Gradient of kinetic energy / Vertical advection} 
    481 \label{subsec:C_zad} 
     481\label{subsec:INVARIANTS_zad} 
    482482 
    483483The change of Kinetic Energy (KE) due to the vertical advection is exactly balanced by the change of KE due to the horizontal gradient of KE~: 
     
    542542  % 
    543543  \intertext{The first term provides the discrete expression for the vertical advection of momentum (ZAD), 
    544     while the second term corresponds exactly to \autoref{eq:KE+PE_vect_discrete}, therefore:} 
     544    while the second term corresponds exactly to \autoref{eq:INVARIANTS_KE+PE_vect_discrete}, therefore:} 
    545545  \equiv&                   \int\limits_D \textbf{U}_h \cdot \text{ZAD} \;dv 
    546546  + \frac{1}{2} \int_D { {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t  (e_3)  \;dv }    &&&\\ 
     
    578578which is (over-)satified by defining the vertical scale factor as follows: 
    579579\begin{flalign*} 
    580   % \label{eq:e3u-e3v} 
     580  % \label{eq:INVARIANTS_e3u-e3v} 
    581581  e_{3u} = \frac{1}{e_{1u}\,e_{2u}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,i+1/2}    \\ 
    582582  e_{3v} = \frac{1}{e_{1v}\,e_{2v}}\;\overline{ e_{1t}^{ }\,e_{2t}^{ }\,e_{3t}^{ } }^{\,j+1/2} 
     
    590590% ------------------------------------------------------------------------------------------------------------- 
    591591\subsection{Pressure gradient term} 
    592 \label{subsec:C.2.6} 
     592\label{subsec:INVARIANTS_2.6} 
    593593 
    594594\gmcomment{ 
     
    622622  \allowdisplaybreaks 
    623623  \intertext{Using successively \autoref{eq:DOM_di_adj}, \ie\ the skew symmetry property of 
    624     the $\delta$ operator, \autoref{eq:wzv}, the continuity equation, \autoref{eq:dynhpg_sco}, 
     624    the $\delta$ operator, \autoref{eq:DYN_wzv}, the continuity equation, \autoref{eq:DYN_hpg_sco}, 
    625625    the hydrostatic equation in the $s$-coordinate, and $\delta_{k+1/2} \left[ z_t \right] \equiv e_{3w} $, 
    626626    which comes from the definition of $z_t$, it becomes: } 
     
    667667  % 
    668668\end{flalign*} 
    669 The first term is exactly the first term of the right-hand-side of \autoref{eq:KE+PE_vect_discrete}. 
     669The first term is exactly the first term of the right-hand-side of \autoref{eq:INVARIANTS_KE+PE_vect_discrete}. 
    670670It remains to demonstrate that the last term, 
    671671which is obviously a discrete analogue of $\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to 
    672 the last term of \autoref{eq:KE+PE_vect_discrete}. 
     672the last term of \autoref{eq:INVARIANTS_KE+PE_vect_discrete}. 
    673673In other words, the following property must be satisfied: 
    674674\begin{flalign*} 
     
    735735% ================================================================ 
    736736\section{Discrete total energy conservation: flux form} 
    737 \label{sec:C.3} 
     737\label{sec:INVARIANTS_3} 
    738738 
    739739% ------------------------------------------------------------------------------------------------------------- 
     
    741741% ------------------------------------------------------------------------------------------------------------- 
    742742\subsection{Total energy conservation} 
    743 \label{subsec:C_KE+PE_flux} 
    744  
    745 The discrete form of the total energy conservation, \autoref{eq:Tot_Energy}, is given by: 
     743\label{subsec:INVARIANTS_KE+PE_flux} 
     744 
     745The discrete form of the total energy conservation, \autoref{eq:INVARIANTS_Tot_Energy}, is given by: 
    746746\begin{flalign*} 
    747747  \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  \\ 
     
    765765% ------------------------------------------------------------------------------------------------------------- 
    766766\subsection{Coriolis and advection terms: flux form} 
    767 \label{subsec:C.3.2} 
     767\label{subsec:INVARIANTS_3.2} 
    768768 
    769769% ------------------------------------------------------------------------------------------------------------- 
     
    771771% ------------------------------------------------------------------------------------------------------------- 
    772772\subsubsection{Coriolis plus ``metric'' term} 
    773 \label{subsec:C.3.3} 
     773\label{subsec:INVARIANTS_3.3} 
    774774 
    775775In flux from the vorticity term reduces to a Coriolis term in which 
     
    786786Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form. 
    787787It therefore conserves the total KE. 
    788 The derivation is the same as for the vorticity term in the vector invariant form (\autoref{subsec:C_vor}). 
     788The derivation is the same as for the vorticity term in the vector invariant form (\autoref{subsec:INVARIANTS_vor}). 
    789789 
    790790% ------------------------------------------------------------------------------------------------------------- 
     
    792792% ------------------------------------------------------------------------------------------------------------- 
    793793\subsubsection{Flux form advection} 
    794 \label{subsec:C.3.4} 
     794\label{subsec:INVARIANTS_3.4} 
    795795 
    796796The flux form operator of the momentum advection is evaluated using 
     
    800800 
    801801\begin{equation} 
    802   \label{eq:C_ADV_KE_flux} 
     802  \label{eq:INVARIANTS_ADV_KE_flux} 
    803803  -  \int_D \textbf{U}_h \cdot     \left(                 {{ 
    804804        \begin{array} {*{20}c} 
     
    863863\] 
    864864which is the discrete form of $ \frac{1}{2} \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv $. 
    865 \autoref{eq:C_ADV_KE_flux} is thus satisfied. 
     865\autoref{eq:INVARIANTS_ADV_KE_flux} is thus satisfied. 
    866866 
    867867When the UBS scheme is used to evaluate the flux form momentum advection, 
     
    873873% ================================================================ 
    874874\section{Discrete enstrophy conservation} 
    875 \label{sec:C.4} 
     875\label{sec:INVARIANTS_4} 
    876876 
    877877% ------------------------------------------------------------------------------------------------------------- 
     
    879879% ------------------------------------------------------------------------------------------------------------- 
    880880\subsubsection{Vorticity term with ENS scheme  (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 
    881 \label{subsec:C_vorENS} 
     881\label{subsec:INVARIANTS_vorENS} 
    882882 
    883883In the ENS scheme, the vorticity term is descretized as follows: 
    884884\begin{equation} 
    885   \tag{\ref{eq:dynvor_ens}} 
     885  \label{eq:INVARIANTS_dynvor_ens} 
    886886  \left\{ 
    887887    \begin{aligned} 
     
    898898it can be shown that: 
    899899\begin{equation} 
    900   \label{eq:C_1.1} 
     900  \label{eq:INVARIANTS_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. 
    904 Indeed, using \autoref{eq:dynvor_ens}, 
    905 the discrete form of the right hand side of \autoref{eq:C_1.1} can be transformed as follow: 
     904Indeed, using \autoref{eq:DYN_vor_ens}, 
     905the discrete form of the right hand side of \autoref{eq:INVARIANTS_1.1} can be transformed as follow: 
    906906\begin{flalign*} 
    907907  &\int_D q \,\; \textbf{k} \cdot \frac{1} {e_3 } \nabla \times 
     
    948948% ------------------------------------------------------------------------------------------------------------- 
    949949\subsubsection{Vorticity Term with EEN scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 
    950 \label{subsec:C_vorEEN} 
     950\label{subsec:INVARIANTS_vorEEN} 
    951951 
    952952With the EEN scheme, the vorticity terms are represented as: 
    953953\begin{equation} 
    954   \tag{\ref{eq:dynvor_een}} 
     954  \label{eq:INVARIANTS_dynvor_een} 
    955955  \left\{ { 
    956956      \begin{aligned} 
     
    966966and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: 
    967967\begin{equation} 
    968   \tag{\ref{eq:Q_triads}} 
     968  \tag{\ref{eq:INVARIANTS_Q_triads}} 
    969969  _i^j \mathbb{Q}^{i_p}_{j_p} 
    970970  = \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) 
     
    975975Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2} $, 
    976976similar manipulation can be done for the 3 others. 
    977 The discrete form of the right hand side of \autoref{eq:C_1.1} applied to 
     977The discrete form of the right hand side of \autoref{eq:INVARIANTS_1.1} applied to 
    978978this triad only can be transformed as follow: 
    979979 
     
    10211021% ================================================================ 
    10221022\section{Conservation properties on tracers} 
    1023 \label{sec:C.5} 
     1023\label{sec:INVARIANTS_5} 
    10241024 
    10251025All the numerical schemes used in \NEMO\ are written such that the tracer content is conserved by 
     
    10371037% ------------------------------------------------------------------------------------------------------------- 
    10381038\subsection{Advection term} 
    1039 \label{subsec:C.5.1} 
     1039\label{subsec:INVARIANTS_5.1} 
    10401040 
    10411041conservation of a tracer, $T$: 
     
    11031103% ================================================================ 
    11041104\section{Conservation properties on lateral momentum physics} 
    1105 \label{sec:dynldf_properties} 
     1105\label{sec:INVARIANTS_dynldf_properties} 
    11061106 
    11071107The discrete formulation of the horizontal diffusion of momentum ensures 
     
    11241124% ------------------------------------------------------------------------------------------------------------- 
    11251125\subsection{Conservation of potential vorticity} 
    1126 \label{subsec:C.6.1} 
     1126\label{subsec:INVARIANTS_6.1} 
    11271127 
    11281128The lateral momentum diffusion term conserves the potential vorticity: 
     
    11581158% ------------------------------------------------------------------------------------------------------------- 
    11591159\subsection{Dissipation of horizontal kinetic energy} 
    1160 \label{subsec:C.6.2} 
     1160\label{subsec:INVARIANTS_6.2} 
    11611161 
    11621162The lateral momentum diffusion term dissipates the horizontal kinetic energy: 
     
    12101210% ------------------------------------------------------------------------------------------------------------- 
    12111211\subsection{Dissipation of enstrophy} 
    1212 \label{subsec:C.6.3} 
     1212\label{subsec:INVARIANTS_6.3} 
    12131213 
    12141214The lateral momentum diffusion term dissipates the enstrophy when the eddy coefficients are horizontally uniform: 
     
    12341234% ------------------------------------------------------------------------------------------------------------- 
    12351235\subsection{Conservation of horizontal divergence} 
    1236 \label{subsec:C.6.4} 
     1236\label{subsec:INVARIANTS_6.4} 
    12371237 
    12381238When the horizontal divergence of the horizontal diffusion of momentum (discrete sense) is taken, 
     
    12611261% ------------------------------------------------------------------------------------------------------------- 
    12621262\subsection{Dissipation of horizontal divergence variance} 
    1263 \label{subsec:C.6.5} 
     1263\label{subsec:INVARIANTS_6.5} 
    12641264 
    12651265\begin{flalign*} 
     
    12871287% ================================================================ 
    12881288\section{Conservation properties on vertical momentum physics} 
    1289 \label{sec:C.7} 
     1289\label{sec:INVARIANTS_7} 
    12901290 
    12911291As for the lateral momentum physics, 
     
    14581458% ================================================================ 
    14591459\section{Conservation properties on tracer physics} 
    1460 \label{sec:C.8} 
     1460\label{sec:INVARIANTS_8} 
    14611461 
    14621462The numerical schemes used for tracer subgridscale physics are written such that 
     
    14701470% ------------------------------------------------------------------------------------------------------------- 
    14711471\subsection{Conservation of tracers} 
    1472 \label{subsec:C.8.1} 
     1472\label{subsec:INVARIANTS_8.1} 
    14731473 
    14741474constraint of conservation of tracers: 
     
    15031503% ------------------------------------------------------------------------------------------------------------- 
    15041504\subsection{Dissipation of tracer variance} 
    1505 \label{subsec:C.8.2} 
     1505\label{subsec:INVARIANTS_8.2} 
    15061506 
    15071507constraint on the dissipation of tracer variance: 
  • NEMO/trunk/doc/latex/NEMO/subfiles/apdx_s_coord.tex

    r11529 r11543  
    77% ================================================================ 
    88\chapter{Curvilinear $s-$Coordinate Equations} 
    9 \label{apdx:A} 
     9\label{apdx:SCOORD} 
    1010 
    1111\chaptertoc 
     
    2828% ================================================================ 
    2929\section{Chain rule for $s-$coordinates} 
    30 \label{sec:A_chain} 
     30\label{sec:SCOORD_chain} 
    3131 
    3232In order to establish the set of Primitive Equation in curvilinear $s$-coordinates 
    3333(\ie\ an orthogonal curvilinear coordinate in the horizontal and 
    3434an Arbitrary Lagrangian Eulerian (ALE) coordinate in the vertical), 
    35 we start from the set of equations established in \autoref{subsec:PE_zco_Eq} for 
     35we start from the set of equations established in \autoref{subsec:MB_zco_Eq} for 
    3636the special case $k = z$ and thus $e_3 = 1$, 
    3737and we introduce an arbitrary vertical coordinate $a = a(i,j,z,t)$. 
     
    3939the horizontal slope of $s-$surfaces by: 
    4040\begin{equation} 
    41   \label{apdx:A_s_slope} 
     41  \label{eq:SCOORD_s_slope} 
    4242  \sigma_1 =\frac{1}{e_1 } \; \left. {\frac{\partial z}{\partial i}} \right|_s 
    4343  \quad \text{and} \quad 
     
    4646 
    4747The model fields (e.g. pressure $p$) can be viewed as functions of $(i,j,z,t)$ (e.g. $p(i,j,z,t)$) or as 
    48 functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of  
    49 these fields.  Any ``infinitesimal'' change in $\bullet$ can be written in two forms:  
    50 \begin{equation} 
    51   \label{apdx:A_s_infin_changes} 
     48functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of 
     49these fields.  Any ``infinitesimal'' change in $\bullet$ can be written in two forms: 
     50\begin{equation} 
     51  \label{eq:SCOORD_s_infin_changes} 
    5252  \begin{aligned} 
    53     & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t}  
    54                 + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t}  
    55                 + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t}  
     53    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t} 
     54                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t} 
     55                + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t} 
    5656                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,s} , \\ 
    57     & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t}  
    58                 + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t}  
    59                 + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t}  
     57    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t} 
     58                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t} 
     59                + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t} 
    6060                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,z} . 
    6161  \end{aligned} 
     
    6363Using the first form and considering a change $\delta i$ with $j, z$ and $t$ held constant, shows that 
    6464\begin{equation} 
    65   \label{apdx:A_s_chain_rule} 
     65  \label{eq:SCOORD_s_chain_rule} 
    6666      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,z,t}  = 
    6767      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,s,t} 
    68     + \left. {\frac{\partial s       }{\partial i}} \right|_{j,z,t} \;  
    69       \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} .     
    70 \end{equation} 
    71 The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces,  
    72 (\autoref{apdx:A_s_slope}), by applying the second of (\autoref{apdx:A_s_infin_changes}) with $\bullet$ set to  
     68    + \left. {\frac{\partial s       }{\partial i}} \right|_{j,z,t} \; 
     69      \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} . 
     70\end{equation} 
     71The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces, 
     72(\autoref{eq:SCOORD_s_slope}), by applying the second of (\autoref{eq:SCOORD_s_infin_changes}) with $\bullet$ set to 
    7373$s$ and $j, t$ held constant 
    7474\begin{equation} 
    75 \label{apdx:a_delta_s} 
    76 \delta s|_{j,t} =  
    77          \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t}  
     75\label{eq:SCOORD_delta_s} 
     76\delta s|_{j,t} = 
     77         \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t} 
    7878       + \delta z \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} . 
    7979\end{equation} 
    8080Choosing to look at a direction in the $(i,z)$ plane in which $\delta s = 0$ and using 
    81 (\autoref{apdx:A_s_slope}) we obtain  
    82 \begin{equation} 
    83 \left. \frac{ \partial s }{\partial i} \right|_{j,z,t} =   
     81(\autoref{eq:SCOORD_s_slope}) we obtain 
     82\begin{equation} 
     83\left. \frac{ \partial s }{\partial i} \right|_{j,z,t} = 
    8484         -  \left. \frac{ \partial z }{\partial i} \right|_{j,s,t} \; 
    8585            \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} 
    8686    = - \frac{e_1 }{e_3 }\sigma_1  . 
    87 \label{apdx:a_ds_di_z} 
    88 \end{equation} 
    89 Another identity, similar in form to (\autoref{apdx:a_ds_di_z}), can be derived 
    90 by choosing $\bullet$ to be $s$ and using the second form of (\autoref{apdx:A_s_infin_changes}) to consider 
     87\label{eq:SCOORD_ds_di_z} 
     88\end{equation} 
     89Another identity, similar in form to (\autoref{eq:SCOORD_ds_di_z}), can be derived 
     90by choosing $\bullet$ to be $s$ and using the second form of (\autoref{eq:SCOORD_s_infin_changes}) to consider 
    9191changes in which $i , j$ and $s$ are constant. This shows that 
    9292\begin{equation} 
    93 \label{apdx:A_w_in_s} 
    94 w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} =   
     93\label{eq:SCOORD_w_in_s} 
     94w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} = 
    9595- \left. \frac{ \partial z }{\partial s} \right|_{i,j,t} 
    96   \left. \frac{ \partial s }{\partial t} \right|_{i,j,z}  
    97   = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} .  
    98 \end{equation} 
    99  
    100 In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is  
    101 usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish  
     96  \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} 
     97  = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} . 
     98\end{equation} 
     99 
     100In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is 
     101usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish 
    102102the model equations in the curvilinear $s-$coordinate system are: 
    103103\begin{equation} 
    104   \label{apdx:A_s_chain_rule} 
     104  \label{eq:SCOORD_s_chain_rule} 
    105105  \begin{aligned} 
    106106    &\left. {\frac{\partial \bullet }{\partial t}} \right|_z  = 
    107     \left. {\frac{\partial \bullet }{\partial t}} \right|_s  
     107    \left. {\frac{\partial \bullet }{\partial t}} \right|_s 
    108108    + \frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial t} , \\ 
    109109    &\left. {\frac{\partial \bullet }{\partial i}} \right|_z  = 
    110110    \left. {\frac{\partial \bullet }{\partial i}} \right|_s 
    111111    +\frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial i}= 
    112     \left. {\frac{\partial \bullet }{\partial i}} \right|_s  
     112    \left. {\frac{\partial \bullet }{\partial i}} \right|_s 
    113113    -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} , \\ 
    114114    &\left. {\frac{\partial \bullet }{\partial j}} \right|_z  = 
    115     \left. {\frac{\partial \bullet }{\partial j}} \right|_s  
     115    \left. {\frac{\partial \bullet }{\partial j}} \right|_s 
    116116    + \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}= 
    117     \left. {\frac{\partial \bullet }{\partial j}} \right|_s  
     117    \left. {\frac{\partial \bullet }{\partial j}} \right|_s 
    118118    - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} , \\ 
    119119    &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} . 
     
    126126% ================================================================ 
    127127\section{Continuity equation in $s-$coordinates} 
    128 \label{sec:A_continuity} 
    129  
    130 Using (\autoref{apdx:A_s_chain_rule}) and 
     128\label{sec:SCOORD_continuity} 
     129 
     130Using (\autoref{eq:SCOORD_s_chain_rule}) and 
    131131the fact that the horizontal scale factors $e_1$ and $e_2$ do not depend on the vertical coordinate, 
    132132the divergence of the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows in order to 
     
    189189\end{subequations} 
    190190 
    191 Here, $w$ is the vertical velocity relative to the $z-$coordinate system.  
    192 Using the first form of (\autoref{apdx:A_s_infin_changes})  
    193 and the definitions (\autoref{apdx:A_s_slope}) and (\autoref{apdx:A_w_in_s}) for $\sigma_1$, $\sigma_2$ and  $w_s$, 
     191Here, $w$ is the vertical velocity relative to the $z-$coordinate system. 
     192Using the first form of (\autoref{eq:SCOORD_s_infin_changes}) 
     193and the definitions (\autoref{eq:SCOORD_s_slope}) and (\autoref{eq:SCOORD_w_in_s}) for $\sigma_1$, $\sigma_2$ and  $w_s$, 
    194194one can show that the vertical velocity, $w_p$ of a point 
    195 moving with the horizontal velocity of the fluid along an $s$ surface is given by  
    196 \begin{equation} 
    197 \label{apdx:A_w_p} 
     195moving with the horizontal velocity of the fluid along an $s$ surface is given by 
     196\begin{equation} 
     197\label{eq:SCOORD_w_p} 
    198198\begin{split} 
    199199w_p  = & \left. \frac{ \partial z }{\partial t} \right|_s 
    200      + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s  
     200     + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s 
    201201     + \frac{v}{e_2} \left. \frac{ \partial z }{\partial j} \right|_s \\ 
    202202     = & w_s + u \sigma_1 + v \sigma_2 . 
    203 \end{split}      
     203\end{split} 
    204204\end{equation} 
    205205 The vertical velocity across this surface is denoted by 
    206206\begin{equation} 
    207   \label{apdx:A_w_s} 
    208   \omega  = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v )  .  
    209 \end{equation} 
    210 Hence  
    211 \begin{equation} 
    212 \frac{1}{e_3 } \frac{\partial}{\partial s}   \left[  w -  u\;\sigma_1  - v\;\sigma_2  \right] =  
    213 \frac{1}{e_3 } \frac{\partial}{\partial s} \left[  \omega + w_s \right] =  
    214    \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s}  
    215  + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] =  
    216    \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s  
    217 \end{equation} 
    218  
    219 Using (\autoref{apdx:A_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain  
     207  \label{eq:SCOORD_w_s} 
     208  \omega  = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v )  . 
     209\end{equation} 
     210Hence 
     211\begin{equation} 
     212\frac{1}{e_3 } \frac{\partial}{\partial s}   \left[  w -  u\;\sigma_1  - v\;\sigma_2  \right] = 
     213\frac{1}{e_3 } \frac{\partial}{\partial s} \left[  \omega + w_s \right] = 
     214   \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s} 
     215 + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] = 
     216   \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s 
     217\end{equation} 
     218 
     219Using (\autoref{eq:SCOORD_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain 
    220220our final expression for the divergence of the velocity in the curvilinear $s-$coordinate system: 
    221221\begin{equation} 
     
    228228\end{equation} 
    229229 
    230 As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is: 
    231 \begin{equation} 
    232   \label{apdx:A_sco_Continuity} 
     230As a result, the continuity equation \autoref{eq:MB_PE_continuity} in the $s-$coordinates is: 
     231\begin{equation} 
     232  \label{eq:SCOORD_sco_Continuity} 
    233233  \frac{1}{e_3 } \frac{\partial e_3}{\partial t} 
    234234  + \frac{1}{e_1 \,e_2 \,e_3 }\left[ 
     
    245245% ================================================================ 
    246246\section{Momentum equation in $s-$coordinate} 
    247 \label{sec:A_momentum} 
     247\label{sec:SCOORD_momentum} 
    248248 
    249249Here we only consider the first component of the momentum equation, 
     
    252252$\bullet$ \textbf{Total derivative in vector invariant form} 
    253253 
    254 Let us consider \autoref{eq:PE_dyn_vect}, the first component of the momentum equation in the vector invariant form. 
     254Let us consider \autoref{eq:MB_dyn_vect}, the first component of the momentum equation in the vector invariant form. 
    255255Its total $z-$coordinate time derivative, 
    256256$\left. \frac{D u}{D t} \right|_z$ can be transformed as follows in order to obtain 
     
    272272        +  w \;\frac{\partial u}{\partial z}      \\ 
    273273        % 
    274       \intertext{introducing the chain rule (\autoref{apdx:A_s_chain_rule}) } 
     274      \intertext{introducing the chain rule (\autoref{eq:SCOORD_s_chain_rule}) } 
    275275      % 
    276276      &= \left. {\frac{\partial u }{\partial t}} \right|_z 
     
    306306        \; \frac{\partial u}{\partial s} .  \\ 
    307307        % 
    308       \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{apdx:A_w_s}) } 
     308      \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{eq:SCOORD_w_s}) } 
    309309      % 
    310310      &= \left. {\frac{\partial u }{\partial t}} \right|_z 
     
    317317\end{subequations} 
    318318% 
    319 Applying the time derivative chain rule (first equation of (\autoref{apdx:A_s_chain_rule})) to $u$ and 
    320 using (\autoref{apdx:A_w_in_s}) provides the expression of the last term of the right hand side, 
     319Applying the time derivative chain rule (first equation of (\autoref{eq:SCOORD_s_chain_rule})) to $u$ and 
     320using (\autoref{eq:SCOORD_w_in_s}) provides the expression of the last term of the right hand side, 
    321321\[ 
    322322  { 
     
    331331\ie\ the total $s-$coordinate time derivative : 
    332332\begin{align} 
    333   \label{apdx:A_sco_Dt_vect} 
     333  \label{eq:SCOORD_sco_Dt_vect} 
    334334  \left. \frac{D u}{D t} \right|_s 
    335335  = \left. {\frac{\partial u }{\partial t}} \right|_s 
    336336  - \left. \zeta \right|_s \;v 
    337337  + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s 
    338   + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} .  
     338  + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} . 
    339339\end{align} 
    340340Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in 
     
    345345 
    346346Let us start from the total time derivative in the curvilinear $s-$coordinate system we have just establish. 
    347 Following the procedure used to establish (\autoref{eq:PE_flux_form}), it can be transformed into : 
     347Following the procedure used to establish (\autoref{eq:MB_flux_form}), it can be transformed into : 
    348348% \begin{subequations} 
    349349\begin{align*} 
     
    367367\end{align*} 
    368368% 
    369 Introducing the vertical scale factor inside the horizontal derivative of the first two terms  
     369Introducing the vertical scale factor inside the horizontal derivative of the first two terms 
    370370(\ie\ the horizontal divergence), it becomes : 
    371371\begin{align*} 
     
    373373  \begin{array}{*{20}l} 
    374374    % \begin{align*} {\begin{array}{*{20}l} 
    375     %     {\begin{array}{*{20}l} \left. \frac{D u}{D t} \right|_s   
     375    %     {\begin{array}{*{20}l} \left. \frac{D u}{D t} \right|_s 
    376376    &= \left. {\frac{\partial u }{\partial t}} \right|_s 
    377377    &+ \frac{1}{e_1\,e_2\,e_3}  \left(  \frac{\partial( e_2 e_3 \,u^2 )}{\partial i} 
     
    398398     % 
    399399    \intertext {Introducing a more compact form for the divergence of the momentum fluxes, 
    400     and using (\autoref{apdx:A_sco_Continuity}), the $s-$coordinate continuity equation, 
     400    and using (\autoref{eq:SCOORD_sco_Continuity}), the $s-$coordinate continuity equation, 
    401401    it becomes : } 
    402402  % 
     
    410410  } 
    411411\end{align*} 
    412 which leads to the $s-$coordinate flux formulation of the total $s-$coordinate time derivative,  
     412which leads to the $s-$coordinate flux formulation of the total $s-$coordinate time derivative, 
    413413\ie\ the total $s-$coordinate time derivative in flux form: 
    414414\begin{flalign} 
    415   \label{apdx:A_sco_Dt_flux} 
     415  \label{eq:SCOORD_sco_Dt_flux} 
    416416  \left. \frac{D u}{D t} \right|_s   = \frac{1}{e_3}  \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s 
    417417  + \left.  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)    \right|_s 
     
    422422It has the same form as in the $z-$coordinate but for 
    423423the vertical scale factor that has appeared inside the time derivative which 
    424 comes from the modification of (\autoref{apdx:A_sco_Continuity}), 
     424comes from the modification of (\autoref{eq:SCOORD_sco_Continuity}), 
    425425the continuity equation. 
    426426 
     
    437437\] 
    438438Applying similar manipulation to the second component and 
    439 replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it becomes: 
    440 \begin{equation} 
    441   \label{apdx:A_grad_p_1} 
     439replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{eq:SCOORD_s_slope}, it becomes: 
     440\begin{equation} 
     441  \label{eq:SCOORD_grad_p_1} 
    442442  \begin{split} 
    443443    -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z 
     
    451451\end{equation} 
    452452 
    453 An additional term appears in (\autoref{apdx:A_grad_p_1}) which accounts for 
     453An additional term appears in (\autoref{eq:SCOORD_grad_p_1}) which accounts for 
    454454the tilt of $s-$surfaces with respect to geopotential $z-$surfaces. 
    455455 
     
    467467Therefore, $p$ and $p_h'$ are linked through: 
    468468\begin{equation} 
    469   \label{apdx:A_pressure} 
     469  \label{eq:SCOORD_pressure} 
    470470  p = \rho_o \; p_h' + \rho_o \, g \, ( \eta - z ) 
    471471\end{equation} 
     
    475475\] 
    476476 
    477 Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and 
     477Substituing \autoref{eq:SCOORD_pressure} in \autoref{eq:SCOORD_grad_p_1} and 
    478478using the definition of the density anomaly it becomes an expression in two parts: 
    479479\begin{equation} 
    480   \label{apdx:A_grad_p_2} 
     480  \label{eq:SCOORD_grad_p_2} 
    481481  \begin{split} 
    482482    -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z 
     
    491491This formulation of the pressure gradient is characterised by the appearance of 
    492492a term depending on the sea surface height only 
    493 (last term on the right hand side of expression \autoref{apdx:A_grad_p_2}). 
     493(last term on the right hand side of expression \autoref{eq:SCOORD_grad_p_2}). 
    494494This term will be loosely termed \textit{surface pressure gradient} whereas 
    495495the first term will be termed the \textit{hydrostatic pressure gradient} by analogy to 
     
    502502The coriolis and forcing terms as well as the the vertical physics remain unchanged as 
    503503they involve neither time nor space derivatives. 
    504 The form of the lateral physics is discussed in \autoref{apdx:B}. 
     504The form of the lateral physics is discussed in \autoref{apdx:DIFFOPERS}. 
    505505 
    506506$\bullet$ \textbf{Full momentum equation} 
     
    510510the one in a curvilinear $z-$coordinate, except for the pressure gradient term: 
    511511\begin{subequations} 
    512   \label{apdx:A_dyn_vect} 
     512  \label{eq:SCOORD_dyn_vect} 
    513513  \begin{multline} 
    514     \label{apdx:A_PE_dyn_vect_u} 
     514    \label{eq:SCOORD_PE_dyn_vect_u} 
    515515    \frac{\partial u}{\partial t}= 
    516516    +   \left( {\zeta +f} \right)\,v 
     
    522522  \end{multline} 
    523523  \begin{multline} 
    524     \label{apdx:A_dyn_vect_v} 
     524    \label{eq:SCOORD_dyn_vect_v} 
    525525    \frac{\partial v}{\partial t}= 
    526526    -   \left( {\zeta +f} \right)\,u 
     
    535535the formulation of both the time derivative and the pressure gradient term: 
    536536\begin{subequations} 
    537   \label{apdx:A_dyn_flux} 
     537  \label{eq:SCOORD_dyn_flux} 
    538538  \begin{multline} 
    539     \label{apdx:A_PE_dyn_flux_u} 
     539    \label{eq:SCOORD_PE_dyn_flux_u} 
    540540    \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t} = 
    541541    - \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right) 
     
    547547  \end{multline} 
    548548  \begin{multline} 
    549     \label{apdx:A_dyn_flux_v} 
     549    \label{eq:SCOORD_dyn_flux_v} 
    550550    \frac{1}{e_3}\frac{\partial \left(  e_3\,v  \right) }{\partial t}= 
    551551    -  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,v}   \right) 
     
    554554    -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right) 
    555555    -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j} 
    556     +  D_v^{\vect{U}}  +   F_v^{\vect{U}} .  
     556    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} . 
    557557  \end{multline} 
    558558\end{subequations} 
     
    560560hydrostatic pressure and density anomalies, $p_h'$ and $d=( \frac{\rho}{\rho_o}-1 )$: 
    561561\begin{equation} 
    562   \label{apdx:A_dyn_zph} 
     562  \label{eq:SCOORD_dyn_zph} 
    563563  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 . 
    564564\end{equation} 
     
    569569in particular the pressure gradient. 
    570570By contrast, $\omega$ is not $w$, the third component of the velocity, but the dia-surface velocity component, 
    571 \ie\ the volume flux across the moving $s$-surfaces per unit horizontal area.  
     571\ie\ the volume flux across the moving $s$-surfaces per unit horizontal area. 
    572572 
    573573 
     
    576576% ================================================================ 
    577577\section{Tracer equation} 
    578 \label{sec:A_tracer} 
     578\label{sec:SCOORD_tracer} 
    579579 
    580580The tracer equation is obtained using the same calculation as for the continuity equation and then 
     
    582582 
    583583\begin{multline} 
    584   \label{apdx:A_tracer} 
     584  \label{eq:SCOORD_tracer} 
    585585  \frac{1}{e_3} \frac{\partial \left(  e_3 T  \right)}{\partial t} 
    586586  = -\frac{1}{e_1 \,e_2 \,e_3} 
     
    591591\end{multline} 
    592592 
    593 The expression for the advection term is a straight consequence of (\autoref{apdx:A_sco_Continuity}), 
    594 the expression of the 3D divergence in the $s-$coordinates established above.  
     593The expression for the advection term is a straight consequence of (\autoref{eq:SCOORD_sco_Continuity}), 
     594the expression of the 3D divergence in the $s-$coordinates established above. 
    595595 
    596596\biblio 
  • NEMO/trunk/doc/latex/NEMO/subfiles/apdx_triads.tex

    r11529 r11543  
    1616% ================================================================ 
    1717\chapter{Iso-Neutral Diffusion and Eddy Advection using Triads} 
    18 \label{apdx:triad} 
     18\label{apdx:TRIADS} 
    1919 
    2020\chaptertoc 
     
    4545\begin{description} 
    4646\item[\np{ln\_triad\_iso}] 
    47   See \autoref{sec:taper}. 
     47  See \autoref{sec:TRIADS_taper}. 
    4848  If this is set false (the default), 
    4949  then `iso-neutral' mixing is accomplished within the surface mixed-layer along slopes linearly decreasing with 
    50   depth from the value immediately below the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}). 
     50  depth from the value immediately below the mixed-layer to zero (flat) at the surface (\autoref{sec:TRIADS_lintaper}). 
    5151  This is the same treatment as used in the default implementation 
    52   \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}. 
     52  \autoref{subsec:LDF_slp_iso}; \autoref{fig:LDF_eiv_slp}. 
    5353  Where \np{ln\_triad\_iso} is set true, 
    5454  the vertical skew flux is further reduced to ensure no vertical buoyancy flux, 
    5555  giving an almost pure horizontal diffusive tracer flux within the mixed layer. 
    56   This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:Gerdes-taper} 
     56  This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:TRIADS_Gerdes-taper} 
    5757\item[\np{ln\_botmix\_triad}] 
    58   See \autoref{sec:iso_bdry}.  
     58  See \autoref{sec:TRIADS_iso_bdry}. 
    5959  If this is set false (the default) then the lateral diffusive fluxes 
    60   associated with triads partly masked by topography are neglected.  
    61   If it is set true, however, then these lateral diffusive fluxes are applied,  
     60  associated with triads partly masked by topography are neglected. 
     61  If it is set true, however, then these lateral diffusive fluxes are applied, 
    6262  giving smoother bottom tracer fields at the cost of introducing diapycnal mixing. 
    6363\item[\np{rn\_sw\_triad}] 
     
    7171 
    7272\section{Triad formulation of iso-neutral diffusion} 
    73 \label{sec:iso} 
     73\label{sec:TRIADS_iso} 
    7474 
    7575We have implemented into \NEMO\ a scheme inspired by \citet{griffies.gnanadesikan.ea_JPO98}, 
     
    7979 
    8080The iso-neutral second order tracer diffusive operator for small angles between 
    81 iso-neutral surfaces and geopotentials is given by \autoref{eq:iso_tensor_1}: 
     81iso-neutral surfaces and geopotentials is given by \autoref{eq:TRIADS_iso_tensor_1}: 
    8282\begin{subequations} 
    83   \label{eq:iso_tensor_1} 
     83  \label{eq:TRIADS_iso_tensor_1} 
    8484  \begin{equation} 
    8585    D^{lT}=-\nabla \cdot\vect{f}^{lT}\equiv 
     
    9292  \end{equation} 
    9393  \begin{equation} 
    94     \label{eq:iso_tensor_2} 
     94    \label{eq:TRIADS_iso_tensor_2} 
    9595    \mbox{with}\quad \;\;\Re = 
    9696    \begin{pmatrix} 
     
    113113%  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 
    114114% \end{array} }} \right) 
    115 Here \autoref{eq:PE_iso_slopes}  
     115Here \autoref{eq:MB_iso_slopes} 
    116116\begin{align*} 
    117117  r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} 
     
    128128We will find it useful to consider the fluxes per unit area in $i,j,k$ space; we write 
    129129\[ 
    130   % \label{eq:Fijk} 
     130  % \label{eq:TRIADS_Fijk} 
    131131  \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right). 
    132132\] 
     
    136136 
    137137The off-diagonal terms of the small angle diffusion tensor 
    138 \autoref{eq:iso_tensor_1}, \autoref{eq:iso_tensor_2} produce skew-fluxes along 
     138\autoref{eq:TRIADS_iso_tensor_1}, \autoref{eq:TRIADS_iso_tensor_2} produce skew-fluxes along 
    139139the $i$- and $j$-directions resulting from the vertical tracer gradient: 
    140140\begin{align} 
    141   \label{eq:i13c} 
     141  \label{eq:TRIADS_i13c} 
    142142  f_{13}=&+{A^{lT}} r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+{A^{lT}} r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\ 
    143143  \intertext{and in the k-direction resulting from the lateral tracer gradients} 
    144   \label{eq:i31c} 
     144  \label{eq:TRIADS_i31c} 
    145145  f_{31}+f_{32}=& {A^{lT}} r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+{A^{lT}} r_2\frac{1}{e_1}\frac{\partial T}{\partial i} 
    146146\end{align} 
     
    148148The vertical diffusive flux associated with the $_{33}$ component of the small angle diffusion tensor is 
    149149\begin{equation} 
    150   \label{eq:i33c} 
     150  \label{eq:TRIADS_i33c} 
    151151  f_{33}=-{A^{lT}}(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}. 
    152152\end{equation} 
     
    157157The following description will describe the fluxes on the $i$-$k$ plane. 
    158158 
    159 There is no natural discretization for the $i$-component of the skew-flux, \autoref{eq:i13c}, 
     159There is no natural discretization for the $i$-component of the skew-flux, \autoref{eq:TRIADS_i13c}, 
    160160as although it must be evaluated at $u$-points, 
    161161it involves vertical gradients (both for the tracer and the slope $r_1$), defined at $w$-points. 
    162 Similarly, the vertical skew flux, \autoref{eq:i31c}, 
     162Similarly, the vertical skew flux, \autoref{eq:TRIADS_i31c}, 
    163163is evaluated at $w$-points but involves horizontal gradients defined at $u$-points. 
    164164 
     
    166166 
    167167The straightforward approach to discretize the lateral skew flux 
    168 \autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA, 
    169 \autoref{eq:tra_ldf_iso}, is to calculate a mean vertical gradient at the $u$-point from 
     168\autoref{eq:TRIADS_i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA, 
     169\autoref{eq:TRA_ldf_iso}, is to calculate a mean vertical gradient at the $u$-point from 
    170170the average of the four surrounding vertical tracer gradients, and multiply this by a mean slope at the $u$-point, 
    171171calculated from the averaged surrounding vertical density gradients. 
    172172The total area-integrated skew-flux (flux per unit area in $ijk$ space) from tracer cell $i,k$ to $i+1,k$, 
    173173noting that the $e_{{3}_{i+1/2}^k}$ in the area $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with 
    174 the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer gradient, is then \autoref{eq:tra_ldf_iso} 
     174the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer gradient, is then \autoref{eq:TRA_ldf_iso} 
    175175\[ 
    176176  \left(F_u^{13} \right)_{i+\frac{1}{2}}^k = {A}_{i+\frac{1}{2}}^k 
     
    205205    \includegraphics[width=\textwidth]{Fig_GRIFF_triad_fluxes} 
    206206    \caption{ 
    207       \protect\label{fig:ISO_triad} 
     207      \protect\label{fig:TRIADS_ISO_triad} 
    208208      (a) Arrangement of triads $S_i$ and tracer gradients to 
    209209      give lateral tracer flux from box $i,k$ to $i+1,k$ 
     
    217217the corresponding `triad' slope calculated from the lateral density gradient across the $u$-point divided by 
    218218the vertical density gradient at the same $w$-point as the tracer gradient. 
    219 See \autoref{fig:ISO_triad}a, where the thick lines denote the tracer gradients, 
     219See \autoref{fig:TRIADS_ISO_triad}a, where the thick lines denote the tracer gradients, 
    220220and the thin lines the corresponding triads, with slopes $s_1, \dotsc s_4$. 
    221221The total area-integrated skew-flux from tracer cell $i,k$ to $i+1,k$ 
    222222\begin{multline} 
    223   \label{eq:i13} 
     223  \label{eq:TRIADS_i13} 
    224224  \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = {A}_{i+1}^k a_1 s_1 
    225225  \delta_{k+\frac{1}{2}} \left[ T^{i+1} 
     
    235235This discretization gives a much closer stencil, and disallows the two-point computational modes. 
    236236 
    237 The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at 
    238 the $w$-point $i,k+\frac{1}{2}$ is constructed similarly (\autoref{fig:ISO_triad}b) by 
     237The vertical skew flux \autoref{eq:TRIADS_i31c} from tracer cell $i,k$ to $i,k+1$ at 
     238the $w$-point $i,k+\frac{1}{2}$ is constructed similarly (\autoref{fig:TRIADS_ISO_triad}b) by 
    239239multiplying lateral tracer gradients from each of the four surrounding $u$-points by the appropriate triad slope: 
    240240\begin{multline} 
    241   \label{eq:i31} 
     241  \label{eq:TRIADS_i31} 
    242242  \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  {A}_i^{k+1} a_{1}' 
    243243  s_{1}' \delta_{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1} 
     
    250250(appearing in both the vertical and lateral gradient), 
    251251and the $u$- and $w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the triad as follows 
    252 (see also \autoref{fig:ISO_triad}): 
    253 \begin{equation} 
    254   \label{eq:R} 
     252(see also \autoref{fig:TRIADS_ISO_triad}): 
     253\begin{equation} 
     254  \label{eq:TRIADS_R} 
    255255  _i^k \mathbb{R}_{i_p}^{k_p} 
    256256  =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} 
     
    269269    \includegraphics[width=\textwidth]{Fig_GRIFF_qcells} 
    270270    \caption{ 
    271       \protect\label{fig:qcells} 
     271      \protect\label{fig:TRIADS_qcells} 
    272272      Triad notation for quarter cells. $T$-cells are inside boxes, 
    273273      while the  $i+\fractext{1}{2},k$ $u$-cell is shaded in green and 
     
    278278% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    279279 
    280 Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter cell that is 
     280Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:TRIADS_qcells}) with the quarter cell that is 
    281281the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ $u$-cell and the $i,k+k_p$ $w$-cell. 
    282 Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation, 
     282Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:TRIADS_i13} and \autoref{eq:TRIADS_i31} in this notation, 
    283283we have \eg\ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. 
    284284Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to 
     
    289289and we notate these areas, similarly to the triad slopes, 
    290290as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, 
    291 where \eg\ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, 
    292 and in \autoref{eq:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 
     291where \eg\ in \autoref{eq:TRIADS_i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, 
     292and in \autoref{eq:TRIADS_i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 
    293293 
    294294\subsection{Full triad fluxes} 
     
    299299tracer cell $i,k$ to $i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the form 
    300300\begin{equation} 
    301   \label{eq:i11} 
     301  \label{eq:TRIADS_i11} 
    302302  \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} = 
    303303  - \left( {A}_i^{k+1} a_{1} + {A}_i^{k+1} a_{2} + {A}_i^k 
     
    305305  \frac{\delta_{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}}, 
    306306\end{equation} 
    307 where the areas $a_i$ are as in \autoref{eq:i13}. 
    308 In this case, separating the total lateral flux, the sum of \autoref{eq:i13} and \autoref{eq:i11}, 
     307where the areas $a_i$ are as in \autoref{eq:TRIADS_i13}. 
     308In this case, separating the total lateral flux, the sum of \autoref{eq:TRIADS_i13} and \autoref{eq:TRIADS_i11}, 
    309309into triad components, a lateral tracer flux 
    310310\begin{equation} 
    311   \label{eq:latflux-triad} 
     311  \label{eq:TRIADS_latflux-triad} 
    312312  _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - {A}_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p} 
    313313  \left( 
     
    322322the lateral density flux associated with each triad separately disappears. 
    323323\begin{equation} 
    324   \label{eq:latflux-rho} 
     324  \label{eq:TRIADS_latflux-rho} 
    325325  {\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 
    326326\end{equation} 
     
    328328tracer cell $i,k$ to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. 
    329329 
    330 The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the $_{33}$ component is also expressed in 
     330The squared slope $r_1^2$ in the expression \autoref{eq:TRIADS_i33c} for the $_{33}$ component is also expressed in 
    331331terms of area-weighted squared triad slopes, 
    332332so the area-integrated vertical flux from tracer cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is 
    333333\begin{equation} 
    334   \label{eq:i33} 
     334  \label{eq:TRIADS_i33} 
    335335  \left( F_w^{33} \right) _i^{k+\frac{1}{2}} = 
    336336  - \left( {A}_i^{k+1} a_{1}' s_{1}'^2 
     
    339339    + {A}_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right], 
    340340\end{equation} 
    341 where the areas $a'$ and slopes $s'$ are the same as in \autoref{eq:i31}. 
    342 Then, separating the total vertical flux, the sum of \autoref{eq:i31} and \autoref{eq:i33}, 
     341where the areas $a'$ and slopes $s'$ are the same as in \autoref{eq:TRIADS_i31}. 
     342Then, separating the total vertical flux, the sum of \autoref{eq:TRIADS_i31} and \autoref{eq:TRIADS_i33}, 
    343343into triad components, a vertical flux 
    344344\begin{align} 
    345   \label{eq:vertflux-triad} 
     345  \label{eq:TRIADS_vertflux-triad} 
    346346  _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) 
    347347  &= {A}_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} 
     
    352352    \right) \\ 
    353353  &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right) 
    354     {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2} 
     354    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:TRIADS_vertflux-triad2} 
    355355\end{align} 
    356356may be associated with each triad. 
     
    361361tracer cell $i,k$ to $i,k+1$ must also vanish since it is a sum of four such triad fluxes. 
    362362 
    363 We can explicitly identify (\autoref{fig:qcells}) the triads associated with the $s_i$, $a_i$, 
    364 and $s'_i$, $a'_i$ used in the definition of the $u$-fluxes and $w$-fluxes in \autoref{eq:i31}, 
    365 \autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and \autoref{fig:ISO_triad} to write out 
     363We can explicitly identify (\autoref{fig:TRIADS_qcells}) the triads associated with the $s_i$, $a_i$, 
     364and $s'_i$, $a'_i$ used in the definition of the $u$-fluxes and $w$-fluxes in \autoref{eq:TRIADS_i31}, 
     365\autoref{eq:TRIADS_i13}, \autoref{eq:TRIADS_i11} \autoref{eq:TRIADS_i33} and \autoref{fig:TRIADS_ISO_triad} to write out 
    366366the iso-neutral fluxes at $u$- and $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: 
    367 %(\autoref{fig:ISO_triad}): 
     367%(\autoref{fig:TRIADS_ISO_triad}): 
    368368\begin{flalign} 
    369   \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
     369  \label{eq:TRIADS_iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
    370370  \sum_{\substack{i_p,\,k_p}} 
    371371  \begin{pmatrix} 
     
    376376 
    377377\subsection{Ensuring the scheme does not increase tracer variance} 
    378 \label{subsec:variance} 
     378\label{subsec:TRIADS_variance} 
    379379 
    380380We now require that this operator should not increase the globally-integrated tracer variance. 
     
    400400    &= -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 
    401401    {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ 
    402     &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:dvar_iso_i} 
     402    &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:TRIADS_dvar_iso_i} 
    403403  \end{aligned} 
    404404\end{multline} 
     
    406406the $T$-points $i,k+k_p-\fractext{1}{2}$ (above) and $i,k+k_p+\fractext{1}{2}$ (below) of 
    407407\begin{equation} 
    408   \label{eq:dvar_iso_k} 
     408  \label{eq:TRIADS_dvar_iso_k} 
    409409  _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    410410\end{equation} 
    411411The total variance tendency driven by the triad is the sum of these two. 
    412412Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with 
    413 \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, it is 
     413\autoref{eq:TRIADS_latflux-triad} and \autoref{eq:TRIADS_vertflux-triad}, it is 
    414414\begin{multline*} 
    415415  -{A}_i^k\left \{ 
     
    430430be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by 
    431431\begin{equation} 
    432   \label{eq:V-A} 
     432  \label{eq:TRIADS_V-A} 
    433433  _i^k\mathbb{V}_{i_p}^{k_p} 
    434434  ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k} 
     
    437437the variance tendency reduces to the perfect square 
    438438\begin{equation} 
    439   \label{eq:perfect-square} 
     439  \label{eq:TRIADS_perfect-square} 
    440440  -{A}_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
    441441  \left( 
     
    445445  \right)^2\leq 0. 
    446446\end{equation} 
    447 Thus, the constraint \autoref{eq:V-A} ensures that the fluxes 
    448 (\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated with 
     447Thus, the constraint \autoref{eq:TRIADS_V-A} ensures that the fluxes 
     448(\autoref{eq:TRIADS_latflux-triad}, \autoref{eq:TRIADS_vertflux-triad}) associated with 
    449449a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase the net variance. 
    450450Since the total fluxes are sums of such fluxes from the various triads, this constraint, applied to all triads, 
    451451is sufficient to ensure that the globally integrated variance does not increase. 
    452452 
    453 The expression \autoref{eq:V-A} can be interpreted as a discretization of the global integral 
    454 \begin{equation} 
    455   \label{eq:cts-var} 
     453The expression \autoref{eq:TRIADS_V-A} can be interpreted as a discretization of the global integral 
     454\begin{equation} 
     455  \label{eq:TRIADS_cts-var} 
    456456  \frac{\partial}{\partial t}\int\!\fractext{1}{2} T^2\, dV = 
    457457  \int\!\mathbf{F}\cdot\nabla T\, dV, 
     
    477477\citet{griffies.gnanadesikan.ea_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells, 
    478478defined in terms of the distances between $T$, $u$,$f$ and $w$-points. 
    479 This is the natural discretization of \autoref{eq:cts-var}. 
     479This is the natural discretization of \autoref{eq:TRIADS_cts-var}. 
    480480The \NEMO\ model, however, operates with scale factors instead of grid sizes, 
    481481and scale factors for the quarter cells are not defined. 
    482482Instead, therefore we simply choose 
    483483\begin{equation} 
    484   \label{eq:V-NEMO} 
     484  \label{eq:TRIADS_V-NEMO} 
    485485  _i^k\mathbb{V}_{i_p}^{k_p}=\fractext{1}{4} {b_u}_{i+i_p}^k, 
    486486\end{equation} 
     
    489489the lateral flux from tracer cell $i,k$ to $i+1,k$ reduces to the classical form 
    490490\begin{equation} 
    491   \label{eq:lat-normal} 
     491  \label{eq:TRIADS_lat-normal} 
    492492  -\overline{A}_{\,i+1/2}^k\; 
    493493  \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
     
    497497In fact if the diffusive coefficient is defined at $u$-points, 
    498498so that we employ ${A}_{i+i_p}^k$ instead of  ${A}_i^k$ in the definitions of the triad fluxes 
    499 \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, 
     499\autoref{eq:TRIADS_latflux-triad} and \autoref{eq:TRIADS_vertflux-triad}, 
    500500we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above. 
    501501 
     
    503503 
    504504The iso-neutral fluxes at $u$- and $w$-points are the sums of the triad fluxes that 
    505 cross the $u$- and $w$-faces \autoref{eq:iso_flux}: 
     505cross the $u$- and $w$-faces \autoref{eq:TRIADS_iso_flux}: 
    506506\begin{subequations} 
    507   % \label{eq:alltriadflux} 
     507  % \label{eq:TRIADS_alltriadflux} 
    508508  \begin{flalign*} 
    509     % \label{eq:vect_isoflux} 
     509    % \label{eq:TRIADS_vect_isoflux} 
    510510    \vect{F}_{\mathrm{iso}}(T) &\equiv 
    511511    \sum_{\substack{i_p,\,k_p}} 
     
    515515    \end{pmatrix}, 
    516516  \end{flalign*} 
    517   where \autoref{eq:latflux-triad}: 
     517  where \autoref{eq:TRIADS_latflux-triad}: 
    518518  \begin{align} 
    519     \label{eq:triadfluxu} 
     519    \label{eq:TRIADS_triadfluxu} 
    520520    _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - {A}_i^k{ 
    521521                                          \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}} 
     
    532532                                          -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ 
    533533                                          \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    534                                           \right),\label{eq:triadfluxw} 
     534                                          \right),\label{eq:TRIADS_triadfluxw} 
    535535  \end{align} 
    536   with \autoref{eq:V-NEMO} 
     536  with \autoref{eq:TRIADS_V-NEMO} 
    537537  \[ 
    538     % \label{eq:V-NEMO2} 
     538    % \label{eq:TRIADS_V-NEMO2} 
    539539    _i^k{\mathbb{V}}_{i_p}^{k_p}=\fractext{1}{4} {b_u}_{i+i_p}^k. 
    540540  \] 
    541541\end{subequations} 
    542542 
    543 The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 
     543The divergence of the expression \autoref{eq:TRIADS_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 
    544544each tracer point: 
    545545\[ 
    546   % \label{eq:iso_operator} 
     546  % \label{eq:TRIADS_iso_operator} 
    547547  D_l^T = \frac{1}{b_T} 
    548548  \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k 
     
    555555\item[$\bullet$ horizontal diffusion] 
    556556  The discretization of the diffusion operator recovers the traditional five-point Laplacian 
    557   \autoref{eq:lat-normal} in the limit of flat iso-neutral direction: 
     557  \autoref{eq:TRIADS_lat-normal} in the limit of flat iso-neutral direction: 
    558558  \[ 
    559     % \label{eq:iso_property0} 
     559    % \label{eq:TRIADS_iso_property0} 
    560560    D_l^T = \frac{1}{b_T} \ 
    561561    \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; 
     
    565565 
    566566\item[$\bullet$ implicit treatment in the vertical] 
    567   Only tracer values associated with a single water column appear in the expression \autoref{eq:i33} for 
     567  Only tracer values associated with a single water column appear in the expression \autoref{eq:TRIADS_i33} for 
    568568  the $_{33}$ fluxes, vertical fluxes driven by vertical gradients. 
    569569  This is of paramount importance since it means that a time-implicit algorithm can be used to 
     
    582582\item[$\bullet$ pure iso-neutral operator] 
    583583  The iso-neutral flux of locally referenced potential density is zero. 
    584   See \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}. 
     584  See \autoref{eq:TRIADS_latflux-rho} and \autoref{eq:TRIADS_vertflux-triad2}. 
    585585 
    586586\item[$\bullet$ conservation of tracer] 
    587587  The iso-neutral diffusion conserves tracer content, \ie 
    588588  \[ 
    589     % \label{eq:iso_property1} 
     589    % \label{eq:TRIADS_iso_property1} 
    590590    \sum_{i,j,k} \left\{ D_l^T \      b_T \right\} = 0 
    591591  \] 
     
    595595  The iso-neutral diffusion does not increase the tracer variance, \ie 
    596596  \[ 
    597     % \label{eq:iso_property2} 
     597    % \label{eq:TRIADS_iso_property2} 
    598598    \sum_{i,j,k} \left\{ T \ D_l^T      \ b_T \right\} \leq 0 
    599599  \] 
    600   The property is demonstrated in \autoref{subsec:variance} above. 
     600  The property is demonstrated in \autoref{subsec:TRIADS_variance} above. 
    601601  It is a key property for a diffusion term. 
    602602  It means that it is also a dissipation term, 
     
    608608  The iso-neutral diffusion operator is self-adjoint, \ie 
    609609  \begin{equation} 
    610     \label{eq:iso_property3} 
     610    \label{eq:TRIADS_iso_property3} 
    611611    \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
    612612  \end{equation} 
     
    614614  We just have to apply the same routine. 
    615615  This property can be demonstrated similarly to the proof of the `no increase of tracer variance' property. 
    616   The contribution by a single triad towards the left hand side of \autoref{eq:iso_property3}, 
    617   can be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} and \autoref{eq:dvar_iso_k}. 
    618   This results in a term similar to \autoref{eq:perfect-square}, 
     616  The contribution by a single triad towards the left hand side of \autoref{eq:TRIADS_iso_property3}, 
     617  can be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:TRIADS_dvar_iso_i} and \autoref{eq:TRIADS_dvar_iso_k}. 
     618  This results in a term similar to \autoref{eq:TRIADS_perfect-square}, 
    619619  \[ 
    620     % \label{eq:TScovar} 
     620    % \label{eq:TRIADS_TScovar} 
    621621    - {A}_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
    622622    \left( 
     
    632632  \] 
    633633This is symmetrical in $T $ and $S$, so exactly the same term arises from 
    634 the discretization of this triad's contribution towards the RHS of \autoref{eq:iso_property3}. 
     634the discretization of this triad's contribution towards the RHS of \autoref{eq:TRIADS_iso_property3}. 
    635635\end{description} 
    636636 
    637637\subsection{Treatment of the triads at the boundaries} 
    638 \label{sec:iso_bdry} 
     638\label{sec:TRIADS_iso_bdry} 
    639639 
    640640The triad slope can only be defined where both the grid boxes centred at the end of the arms exist. 
    641641Triads that would poke up through the upper ocean surface into the atmosphere, 
    642642or down into the ocean floor, must be masked out. 
    643 See \autoref{fig:bdry_triads}. 
     643See \autoref{fig:TRIADS_bdry_triads}. 
    644644Surface layer triads \triad{i}{1}{R}{1/2}{-1/2} (magenta) and \triad{i+1}{1}{R}{-1/2}{-1/2} (blue) that 
    645 require density to be specified above the ocean surface are masked (\autoref{fig:bdry_triads}a): 
     645require density to be specified above the ocean surface are masked (\autoref{fig:TRIADS_bdry_triads}a): 
    646646this ensures that lateral tracer gradients produce no flux through the ocean surface. 
    647647However, to prevent surface noise, it is customary to retain the $_{11}$ contributions towards 
    648648the lateral triad fluxes \triad[u]{i}{1}{F}{1/2}{-1/2} and \triad[u]{i+1}{1}{F}{-1/2}{-1/2}; 
    649649this drives diapycnal tracer fluxes. 
    650 Similar comments apply to triads that would intersect the ocean floor (\autoref{fig:bdry_triads}b). 
     650Similar comments apply to triads that would intersect the ocean floor (\autoref{fig:TRIADS_bdry_triads}b). 
    651651Note that both near bottom triad slopes \triad{i}{k}{R}{1/2}{1/2} and \triad{i+1}{k}{R}{-1/2}{1/2} are masked when 
    652652either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, \ie\ the $i,k+1$ $u$-point is masked. 
     
    662662    \includegraphics[width=\textwidth]{Fig_GRIFF_bdry_triads} 
    663663    \caption{ 
    664       \protect\label{fig:bdry_triads} 
     664      \protect\label{fig:TRIADS_bdry_triads} 
    665665      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots), 
    666666      and $i+1/2,1$ $u$-point (blue square). 
     
    683683 
    684684\subsection{ Limiting of the slopes within the interior} 
    685 \label{sec:limit} 
     685\label{sec:TRIADS_limit} 
    686686 
    687687As discussed in \autoref{subsec:LDF_slp_iso}, 
     
    693693It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to geopotentials 
    694694(here the $\sigma_i$ are the slopes of the coordinate surfaces relative to geopotentials) 
    695 \autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate surfaces, so we require 
     695\autoref{eq:MB_slopes_eiv} rather than the slope $r_i$ relative to coordinate surfaces, so we require 
    696696\[ 
    697697  |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01. 
     
    700700Each individual triad slope 
    701701\begin{equation} 
    702   \label{eq:Rtilde} 
     702  \label{eq:TRIADS_Rtilde} 
    703703  _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}} 
    704704\end{equation} 
     
    711711 
    712712\subsection{Tapering within the surface mixed layer} 
    713 \label{sec:taper} 
     713\label{sec:TRIADS_taper} 
    714714 
    715715Additional tapering of the iso-neutral fluxes is necessary within the surface mixed layer. 
     
    717717 
    718718\subsubsection{Linear slope tapering within the surface mixed layer} 
    719 \label{sec:lintaper} 
     719\label{sec:TRIADS_lintaper} 
    720720 
    721721This is the option activated by the default choice \np{ln\_triad\_iso}\forcode{ = .false.}. 
    722722Slopes $\tilde{r}_i$ relative to geopotentials are tapered linearly from their value immediately below 
    723 the mixed layer to zero at the surface, as described in option (c) of \autoref{fig:eiv_slp}, to values 
    724 \begin{equation} 
    725   \label{eq:rmtilde} 
     723the mixed layer to zero at the surface, as described in option (c) of \autoref{fig:LDF_eiv_slp}, to values 
     724\begin{equation} 
     725  \label{eq:TRIADS_rmtilde} 
    726726  \rMLt = -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h, 
    727727\end{equation} 
    728728and then the $r_i$ relative to vertical coordinate surfaces are appropriately adjusted to 
    729729\[ 
    730   % \label{eq:rm} 
     730  % \label{eq:TRIADS_rm} 
    731731  \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h. 
    732732\] 
    733733Thus the diffusion operator within the mixed layer is given by: 
    734734\[ 
    735   % \label{eq:iso_tensor_ML} 
     735  % \label{eq:TRIADS_iso_tensor_ML} 
    736736  D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    737737  \mbox{with}\quad \;\;\Re =\left( {{ 
     
    747747in isopycnal layers immediately below, in the thermocline. 
    748748It is consistent with the way the $\tilde{r}_i$ are tapered within the mixed layer 
    749 (see \autoref{sec:taperskew} below) so as to ensure a uniform GM eddy-induced velocity throughout the mixed layer. 
     749(see \autoref{sec:TRIADS_taperskew} below) so as to ensure a uniform GM eddy-induced velocity throughout the mixed layer. 
    750750However, it gives a downwards density flux and so acts so as to reduce potential energy in the same way as 
    751 does the slope limiting discussed above in \autoref{sec:limit}. 
    752   
    753 As in \autoref{sec:limit} above, the tapering \autoref{eq:rmtilde} is applied separately to 
     751does the slope limiting discussed above in \autoref{sec:TRIADS_limit}. 
     752 
     753As in \autoref{sec:TRIADS_limit} above, the tapering \autoref{eq:TRIADS_rmtilde} is applied separately to 
    754754each triad $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. 
    755755For clarity, we assume $z$-coordinates in the following; 
    756756the conversion from $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as 
    757 described above by \autoref{eq:Rtilde}. 
     757described above by \autoref{eq:TRIADS_Rtilde}. 
    758758\begin{enumerate} 
    759759\item 
    760760  Mixed-layer depth is defined so as to avoid including regions of weak vertical stratification in 
    761761  the slope definition. 
    762   At each $i,j$ (simplified to $i$ in \autoref{fig:MLB_triad}), 
     762  At each $i,j$ (simplified to $i$ in \autoref{fig:TRIADS_MLB_triad}), 
    763763  we define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer, 
    764764  $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) such that 
    765765  the potential density ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, 
    766766  where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10~m. 
    767   See the left side of \autoref{fig:MLB_triad}. 
     767  See the left side of \autoref{fig:TRIADS_MLB_triad}. 
    768768  We use the $k_{10}$-gridbox instead of the surface gridbox to avoid problems \eg\ with thin daytime mixed-layers. 
    769769  Currently we use the same $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is used to 
     
    776776  This is to ensure that the vertical density gradients associated with 
    777777  these basal triad slopes ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are representative of the thermocline. 
    778   The four basal triads defined in the bottom part of \autoref{fig:MLB_triad} are then 
     778  The four basal triads defined in the bottom part of \autoref{fig:TRIADS_MLB_triad} are then 
    779779  \begin{align*} 
    780780    {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &= 
    781781                                                       {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}, 
    782                                                        % \label{eq:Rbase} 
     782                                                       % \label{eq:TRIADS_Rbase} 
    783783    \\ 
    784784    \intertext{with \eg\ the green triad} 
     
    789789the $w$-point $i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point, so it is this depth 
    790790\[ 
    791   % \label{eq:zbase} 
     791  % \label{eq:TRIADS_zbase} 
    792792  {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2} 
    793793\] 
    794794one gridbox deeper than the diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper the slopes in 
    795 \autoref{eq:rmtilde}. 
     795\autoref{eq:TRIADS_rmtilde}. 
    796796\item 
    797797  Finally, we calculate the adjusted triads ${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within 
     
    805805    {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p} &= 
    806806                                                       \frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}. 
    807                                                        % \label{eq:RML} 
     807                                                       % \label{eq:TRIADS_RML} 
    808808  \end{align*} 
    809809\end{enumerate} 
     
    813813%  \fcapside { 
    814814  \caption{ 
    815     \protect\label{fig:MLB_triad} 
     815    \protect\label{fig:TRIADS_MLB_triad} 
    816816    Definition of mixed-layer depth and calculation of linearly tapered triads. 
    817817    The figure shows a water column at a given $i,j$ (simplified to $i$), with the ocean surface at the top. 
     
    836836 
    837837\subsubsection{Additional truncation of skew iso-neutral flux components} 
    838 \label{subsec:Gerdes-taper} 
     838\label{subsec:TRIADS_Gerdes-taper} 
    839839 
    840840The alternative option is activated by setting \np{ln\_triad\_iso} = true. 
     
    843843but replaces the $\rML$ in the skew term by 
    844844\begin{equation} 
    845   \label{eq:rm*} 
     845  \label{eq:TRIADS_rm*} 
    846846  \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i, 
    847847\end{equation} 
    848848giving a ML diffusive operator 
    849849\[ 
    850   % \label{eq:iso_tensor_ML2} 
     850  % \label{eq:TRIADS_iso_tensor_ML2} 
    851851  D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    852852  \mbox{with}\quad \;\;\Re =\left( {{ 
     
    881881% ================================================================ 
    882882\section{Eddy induced advection formulated as a skew flux} 
    883 \label{sec:skew-flux} 
     883\label{sec:TRIADS_skew-flux} 
    884884 
    885885\subsection{Continuous skew flux formulation} 
    886 \label{sec:continuous-skew-flux} 
     886\label{sec:TRIADS_continuous-skew-flux} 
    887887 
    888888When Gent and McWilliams's [1990] diffusion is used, an additional advection term is added. 
     
    890890the formulation of which depends on the slopes of iso-neutral surfaces. 
    891891Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, 
    892 \ie\ \autoref{eq:ldfslp_geo} is used in $z$-coordinate, 
    893 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. 
     892\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinate, 
     893and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $z^*$ or $s$-coordinates. 
    894894 
    895895The eddy induced velocity is given by: 
    896896\begin{subequations} 
    897   % \label{eq:eiv} 
     897  % \label{eq:TRIADS_eiv} 
    898898  \begin{equation} 
    899     \label{eq:eiv_v} 
     899    \label{eq:TRIADS_eiv_v} 
    900900    \begin{split} 
    901901      u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\ 
     
    907907  where the streamfunctions $\psi_i$ are given by 
    908908  \begin{equation} 
    909     \label{eq:eiv_psi} 
     909    \label{eq:TRIADS_eiv_psi} 
    910910    \begin{split} 
    911911      \psi_1 & = A_{e} \; \tilde{r}_1,   \\ 
     
    958958we end up with the skew form of the eddy induced advective fluxes per unit area in $ijk$ space: 
    959959\begin{equation} 
    960   \label{eq:eiv_skew_ijk} 
     960  \label{eq:TRIADS_eiv_skew_ijk} 
    961961  \textbf{F}_\mathrm{eiv}^T = 
    962962  \begin{pmatrix} 
     
    967967The total fluxes per unit physical area are then 
    968968\begin{equation} 
    969   \label{eq:eiv_skew_physical} 
     969  \label{eq:TRIADS_eiv_skew_physical} 
    970970  \begin{split} 
    971971    f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\ 
     
    974974\end{split} 
    975975\end{equation} 
    976 Note that \autoref{eq:eiv_skew_physical} takes the same form whatever the vertical coordinate, 
    977 though of course the slopes $\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative to 
     976Note that \autoref{eq:TRIADS_eiv_skew_physical} takes the same form whatever the vertical coordinate, 
     977though of course the slopes $\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:TRIADS_eiv_psi} are relative to 
    978978geopotentials. 
    979979The tendency associated with eddy induced velocity is then simply the convergence of the fluxes 
    980 (\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so 
    981 \[ 
    982   % \label{eq:skew_eiv_conv} 
     980(\autoref{eq:TRIADS_eiv_skew_ijk}, \autoref{eq:TRIADS_eiv_skew_physical}), so 
     981\[ 
     982  % \label{eq:TRIADS_skew_eiv_conv} 
    983983  \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[ 
    984984    \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right) 
     
    993993\subsection{Discrete skew flux formulation} 
    994994 
    995 The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}), 
    996 like the off-diagonal terms (\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor, 
    997 are best expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad} and 
    998 (\autoref{eq:i13}, \autoref{eq:i31}); 
     995The skew fluxes in (\autoref{eq:TRIADS_eiv_skew_physical}, \autoref{eq:TRIADS_eiv_skew_ijk}), 
     996like the off-diagonal terms (\autoref{eq:TRIADS_i13c}, \autoref{eq:TRIADS_i31c}) of the small angle diffusion tensor, 
     997are best expressed in terms of the triad slopes, as in \autoref{fig:TRIADS_ISO_triad} and 
     998(\autoref{eq:TRIADS_i13}, \autoref{eq:TRIADS_i31}); 
    999999but now in terms of the triad slopes $\tilde{\mathbb{R}}$ relative to geopotentials instead of 
    10001000the $\mathbb{R}$ relative to coordinate surfaces. 
    1001 The discrete form of \autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and 
     1001The discrete form of \autoref{eq:TRIADS_eiv_skew_ijk} using the slopes \autoref{eq:TRIADS_R} and 
    10021002defining $A_e$ at $T$-points is then given by: 
    10031003 
    10041004\begin{subequations} 
    1005   % \label{eq:allskewflux} 
     1005  % \label{eq:TRIADS_allskewflux} 
    10061006  \begin{flalign*} 
    1007     % \label{eq:vect_skew_flux} 
     1007    % \label{eq:TRIADS_vect_skew_flux} 
    10081008    \vect{F}_{\mathrm{eiv}}(T) &\equiv    \sum_{\substack{i_p,\,k_p}} 
    10091009    \begin{pmatrix} 
     
    10121012    \end{pmatrix}, 
    10131013  \end{flalign*} 
    1014   where the skew flux in the $i$-direction associated with a given triad is (\autoref{eq:latflux-triad}, 
    1015   \autoref{eq:triadfluxu}): 
     1014  where the skew flux in the $i$-direction associated with a given triad is (\autoref{eq:TRIADS_latflux-triad}, 
     1015  \autoref{eq:TRIADS_triadfluxu}): 
    10161016  \begin{align} 
    1017     \label{eq:skewfluxu} 
     1017    \label{eq:TRIADS_skewfluxu} 
    10181018    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \fractext{1}{4} {A_e}_i^k{ 
    10191019                                          \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
     
    10211021                                          \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, \\ 
    10221022    \intertext{ 
    1023     and \autoref{eq:triadfluxw} in the $k$-direction, changing the sign 
    1024     to be consistent with \autoref{eq:eiv_skew_ijk}: 
     1023    and \autoref{eq:TRIADS_triadfluxw} in the $k$-direction, changing the sign 
     1024    to be consistent with \autoref{eq:TRIADS_eiv_skew_ijk}: 
    10251025    } 
    10261026    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T) 
    10271027                                        &= -\fractext{1}{4} {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}} 
    1028                                           {_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} 
     1028                                          {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:TRIADS_skewfluxw} 
    10291029  \end{align} 
    10301030\end{subequations} 
     
    10381038This can be seen %either from Appendix \autoref{apdx:eiv_skew} or 
    10391039by considering the fluxes associated with a given triad slope $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. 
    1040 For, following \autoref{subsec:variance} and \autoref{eq:dvar_iso_i}, 
     1040For, following \autoref{subsec:TRIADS_variance} and \autoref{eq:TRIADS_dvar_iso_i}, 
    10411041the associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ drives a net rate of change of variance, 
    10421042summed over the two $T$-points $i+i_p-\fractext{1}{2},k$ and $i+i_p+\fractext{1}{2},k$, of 
    10431043\begin{equation} 
    1044   \label{eq:dvar_eiv_i} 
     1044  \label{eq:TRIADS_dvar_eiv_i} 
    10451045  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], 
    10461046\end{equation} 
     
    10481048the $T$-points $i,k+k_p-\fractext{1}{2}$ (above) and $i,k+k_p+\fractext{1}{2}$ (below) of 
    10491049\begin{equation} 
    1050   \label{eq:dvar_eiv_k} 
     1050  \label{eq:TRIADS_dvar_eiv_k} 
    10511051  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    10521052\end{equation} 
    1053 Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw}) shows that 
    1054 these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k}) sum to zero. 
     1053Inspection of the definitions (\autoref{eq:TRIADS_skewfluxu}, \autoref{eq:TRIADS_skewfluxw}) shows that 
     1054these two variance changes (\autoref{eq:TRIADS_dvar_eiv_i}, \autoref{eq:TRIADS_dvar_eiv_k}) sum to zero. 
    10551055Hence the two fluxes associated with each triad make no net contribution to the variance budget. 
    10561056 
     
    10641064For the change in gravitational PE driven by the $k$-flux is 
    10651065\begin{align} 
    1066   \label{eq:vert_densityPE} 
     1066  \label{eq:TRIADS_vert_densityPE} 
    10671067  g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) 
    10681068  &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k 
    10691069    {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k 
    10701070    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\ 
    1071   \intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from \autoref{eq:skewfluxw}, gives} 
     1071  \intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from \autoref{eq:TRIADS_skewfluxw}, gives} 
    10721072  % and separating out 
    10731073  % $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$, 
     
    10801080    \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}}, 
    10811081\end{align} 
    1082 using the definition of the triad slope $\rtriad{R}$, \autoref{eq:R} to 
     1082using the definition of the triad slope $\rtriad{R}$, \autoref{eq:TRIADS_R} to 
    10831083express $-\alpha _i^k\delta_{i+ i_p}[T^k]+\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of 
    10841084$-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$. 
     
    10861086Where the coordinates slope, the $i$-flux gives a PE change 
    10871087\begin{multline} 
    1088   \label{eq:lat_densityPE} 
     1088  \label{eq:TRIADS_lat_densityPE} 
    10891089  g \delta_{i+i_p}[z_T^k] 
    10901090  \left[ 
     
    10961096  \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}}, 
    10971097\end{multline} 
    1098 (using \autoref{eq:skewfluxu}) and so the total PE change \autoref{eq:vert_densityPE} + 
    1099 \autoref{eq:lat_densityPE} associated with the triad fluxes is 
     1098(using \autoref{eq:TRIADS_skewfluxu}) and so the total PE change \autoref{eq:TRIADS_vert_densityPE} + 
     1099\autoref{eq:TRIADS_lat_densityPE} associated with the triad fluxes is 
    11001100\begin{multline*} 
    1101   % \label{eq:tot_densityPE} 
     1101  % \label{eq:TRIADS_tot_densityPE} 
    11021102  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) + 
    11031103  g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\ 
     
    11101110 
    11111111\subsection{Treatment of the triads at the boundaries} 
    1112 \label{sec:skew_bdry} 
    1113  
    1114 Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries  
    1115 in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes,  
    1116 as described in \autoref{sec:iso_bdry} and \autoref{fig:bdry_triads}.  
    1117 Thus surface layer triads $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are masked,  
    1118 and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when  
    1119 either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, \ie\ the $i,k+1$ $u$-point is masked.  
     1112\label{sec:TRIADS_skew_bdry} 
     1113 
     1114Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries 
     1115in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes, 
     1116as described in \autoref{sec:TRIADS_iso_bdry} and \autoref{fig:TRIADS_bdry_triads}. 
     1117Thus surface layer triads $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are masked, 
     1118and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when 
     1119either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, \ie\ the $i,k+1$ $u$-point is masked. 
    11201120The namelist parameter \np{ln\_botmix\_triad} has no effect on the eddy-induced skew-fluxes. 
    11211121 
    11221122\subsection{Limiting of the slopes within the interior} 
    1123 \label{sec:limitskew} 
    1124  
    1125 Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$,  
    1126 exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:limit}.  
     1123\label{sec:TRIADS_limitskew} 
     1124 
     1125Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$, 
     1126exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:TRIADS_limit}. 
    11271127Each individual triad \rtriadt{R} is so limited. 
    11281128 
    11291129\subsection{Tapering within the surface mixed layer} 
    1130 \label{sec:taperskew} 
    1131  
    1132 The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R})  
    1133 are always tapered linearly from their value immediately below the mixed layer to zero at the surface  
    1134 \autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}.  
    1135 This is option (c) of \autoref{fig:eiv_slp}.  
    1136 This linear tapering for the slopes used to calculate the eddy-induced fluxes is unaffected by  
     1130\label{sec:TRIADS_taperskew} 
     1131 
     1132The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R}) 
     1133are always tapered linearly from their value immediately below the mixed layer to zero at the surface 
     1134\autoref{eq:TRIADS_rmtilde}, as described in \autoref{sec:TRIADS_lintaper}. 
     1135This is option (c) of \autoref{fig:LDF_eiv_slp}. 
     1136This linear tapering for the slopes used to calculate the eddy-induced fluxes is unaffected by 
    11371137the value of \np{ln\_triad\_iso}. 
    11381138 
     
    11401140the horizontal (the most commonly used options in \NEMO: see \autoref{sec:LDF_coef}), 
    11411141it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer 
    1142 \autoref{eq:eiv_v}. 
     1142\autoref{eq:TRIADS_eiv_v}. 
    11431143This ensures that the eiv velocities do not restratify the mixed layer \citep{treguier.held.ea_JPO97,danabasoglu.ferrari.ea_JC08}. 
    11441144Equivantly, in terms of the skew-flux formulation we use here, 
     
    11481148 
    11491149\subsection{Streamfunction diagnostics} 
    1150 \label{sec:sfdiag} 
     1150\label{sec:TRIADS_sfdiag} 
    11511151 
    11521152Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, 
     
    11541154Each time step, streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at 
    11551155$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$) 
    1156 points (see Table \autoref{tab:cell}) respectively. 
     1156points (see Table \autoref{tab:DOM_cell}) respectively. 
    11571157We follow \citep{griffies_bk04} and calculate the streamfunction at a given $uw$-point from 
    11581158the surrounding four triads according to: 
    11591159\[ 
    1160   % \label{eq:sfdiagi} 
     1160  % \label{eq:TRIADS_sfdiagi} 
    11611161  {\psi_1}_{i+1/2}^{k+1/2}={\fractext{1}{4}}\sum_{\substack{i_p,\,k_p}} 
    11621162  {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}. 
    11631163\] 
    11641164The streamfunction $\psi_1$ is calculated similarly at $vw$ points. 
    1165 The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:eiv_v}: 
    1166 \[ 
    1167   % \label{eq:eiv_v_discrete} 
     1165The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:TRIADS_eiv_v}: 
     1166\[ 
     1167  % \label{eq:TRIADS_eiv_v_discrete} 
    11681168  \begin{split} 
    11691169    {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),   \\ 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_ASM.tex

    r11537 r11543  
    5353additional tendency terms to the prognostic equations: 
    5454\begin{align*} 
    55   % \label{eq:wa_traj_iau} 
     55  % \label{eq:ASM_wa_traj_iau} 
    5656  {\mathbf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\mathbf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\mathbf x}^{a} 
    5757\end{align*} 
     
    6666The first function (namelist option \np{niaufn}=0) employs constant weights, 
    6767\begin{align} 
    68   \label{eq:F1_i} 
     68  \label{eq:ASM_F1_i} 
    6969  F^{(1)}_{i} 
    7070  =\left\{ 
     
    8080with the weighting reduced linearly to a small value at the window end-points: 
    8181\begin{align} 
    82   \label{eq:F2_i} 
     82  \label{eq:ASM_F2_i} 
    8383  F^{(2)}_{i} 
    8484  =\left\{ 
     
    9292\end{align} 
    9393where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. 
    94 The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from 
    95 one assimilation cycle to the next than that described by \autoref{eq:F1_i}. 
     94The weights described by \autoref{eq:ASM_F2_i} provide a smoother transition of the analysis trajectory from 
     95one assimilation cycle to the next than that described by \autoref{eq:ASM_F1_i}. 
    9696 
    9797%========================================================================== 
     
    106106 
    107107\begin{equation} 
    108   \label{eq:asm_dmp} 
     108  \label{eq:ASM_dmp} 
    109109  \left\{ 
    110110    \begin{aligned} 
     
    120120 
    121121\[ 
    122   % \label{eq:asm_div} 
     122  % \label{eq:ASM_div} 
    123123  \chi^{n-1}_I = \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
    124124  \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u^{n-1}_I} \right] 
     
    126126\] 
    127127 
    128 By the application of \autoref{eq:asm_dmp} the divergence is filtered in each iteration, 
     128By the application of \autoref{eq:ASM_dmp} the divergence is filtered in each iteration, 
    129129and the vorticity is left unchanged. 
    130130In the presence of coastal boundaries with zero velocity increments perpendicular to the coast 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_DIA.tex

    r11537 r11543  
    696696\begin{table} 
    697697  \scriptsize 
    698   \begin{tabularx}{\textwidth}{|X|c|c|c|} 
     698  \begin{tabular}{|l|c|c|} 
    699699    \hline 
    700700    tag ids affected by automatic definition of some of their attributes & 
    701701    name attribute                                                       & 
    702     attribute value                      \\ 
     702    attribute value                                                      \\ 
    703703    \hline 
    704704    \hline 
    705705    field\_definition                                                    & 
    706706    freq\_op                                                             & 
    707     \np{rn\_rdt}                         \\ 
     707    \np{rn\_rdt}                                                         \\ 
    708708    \hline 
    709709    SBC                                                                  & 
    710710    freq\_op                                                             & 
    711     \np{rn\_rdt} $\times$ \np{nn\_fsbc}  \\ 
     711    \np{rn\_rdt} $\times$ \np{nn\_fsbc}                                  \\ 
    712712    \hline 
    713713    ptrc\_T                                                              & 
    714714    freq\_op                                                             & 
    715     \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\ 
     715    \np{rn\_rdt} $\times$ \np{nn\_dttrc}                                 \\ 
    716716    \hline 
    717717    diad\_T                                                              & 
    718718    freq\_op                                                             & 
    719     \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\ 
     719    \np{rn\_rdt} $\times$ \np{nn\_dttrc}                                 \\ 
    720720    \hline 
    721721    EqT, EqU, EqW                                                        & 
    722722    jbegin, ni,                                                          & 
    723     according to the grid                \\ 
    724     & 
     723    according to the grid                                                \\ 
     724                                                                         & 
    725725    name\_suffix                                                         & 
    726     \\ 
     726                                                                         \\ 
    727727    \hline 
    728728    TAO, RAMA and PIRATA moorings                                        & 
    729729    zoom\_ibegin, zoom\_jbegin,                                          & 
    730     according to the grid                \\ 
    731     & 
     730    according to the grid                                                \\ 
     731                                                                         & 
    732732    name\_suffix                                                         & 
    733     \\ 
    734     \hline 
    735   \end{tabularx} 
     733                                                                         \\ 
     734    \hline 
     735  \end{tabular} 
    736736\end{table} 
    737737 
     
    739739 
    740740\subsection{XML reference tables} 
    741 \label{subsec:IOM_xmlref} 
     741\label{subsec:DIA_IOM_xmlref} 
    742742 
    743743\begin{enumerate} 
     
    13361336the CF metadata standard. 
    13371337Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of 
    1338 section \autoref{subsec:IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. 
     1338section \autoref{subsec:DIA_IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. 
    13391339 
    13401340Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by 
     
    14071407the mono-processor case (\ie\ global domain of {\small\ttfamily 182x149x31}). 
    14081408An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in 
    1409 table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with 
     1409table \autoref{tab:DIA_NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with 
    14101410a 4x2 mpi partitioning. 
    14111411Note the variation in the compression ratio achieved which reflects chiefly the dry to wet volume ratio of 
     
    14471447  \end{tabular} 
    14481448  \caption{ 
    1449     \protect\label{tab:NC4} 
     1449    \protect\label{tab:DIA_NC4} 
    14501450    Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression 
    14511451  } 
     
    15151515\section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})] 
    15161516{FLO: On-Line Floats trajectories (\protect\key{floats})} 
    1517 \label{sec:FLO} 
     1517\label{sec:DIA_FLO} 
    15181518%--------------------------------------------namflo------------------------------------------------------- 
    15191519 
     
    18471847    \mathcal{V} &=  \mathcal{A}  \;\bar{\eta} 
    18481848  \end{split} 
    1849   \label{eq:MV_nBq} 
     1849  \label{eq:DIA_MV_nBq} 
    18501850\end{equation} 
    18511851 
     
    18551855  \frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} ) 
    18561856  = \left. \frac{\textit{emp}}{e_3}\right|_\textit{surface} 
    1857   \label{eq:Co_nBq} 
     1857  \label{eq:DIA_Co_nBq} 
    18581858\end{equation} 
    18591859 
     
    18641864\begin{equation} 
    18651865  \partial_t \mathcal{M} = \mathcal{A} \;\overline{\textit{emp}} 
    1866   \label{eq:Mass_nBq} 
     1866  \label{eq:DIA_Mass_nBq} 
    18671867\end{equation} 
    18681868 
    18691869where $\overline{\textit{emp}} = \int_S \textit{emp}\,ds$ is the net mass flux through the ocean surface. 
    1870 Bringing \autoref{eq:Mass_nBq} and the time derivative of \autoref{eq:MV_nBq} together leads to 
     1870Bringing \autoref{eq:DIA_Mass_nBq} and the time derivative of \autoref{eq:DIA_MV_nBq} together leads to 
    18711871the evolution equation of the mean sea level 
    18721872 
     
    18741874  \partial_t \bar{\eta} = \frac{\overline{\textit{emp}}}{ \bar{\rho}} 
    18751875  - \frac{\mathcal{V}}{\mathcal{A}} \;\frac{\partial_t \bar{\rho} }{\bar{\rho}} 
    1876   \label{eq:ssh_nBq} 
     1876  \label{eq:DIA_ssh_nBq} 
    18771877\end{equation} 
    18781878 
    1879 The first term in equation \autoref{eq:ssh_nBq} alters sea level by adding or subtracting mass from the ocean. 
     1879The first term in equation \autoref{eq:DIA_ssh_nBq} alters sea level by adding or subtracting mass from the ocean. 
    18801880The second term arises from temporal changes in the global mean density; \ie\ from steric effects. 
    18811881 
    18821882In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$ appears multiplied by 
    18831883the gravity (\ie\ in the hydrostatic balance of the primitive Equations). 
    1884 In particular, the mass conservation equation, \autoref{eq:Co_nBq}, degenerates into the incompressibility equation: 
     1884In particular, the mass conservation equation, \autoref{eq:DIA_Co_nBq}, degenerates into the incompressibility equation: 
    18851885 
    18861886\[ 
    18871887  \frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} ) = \left. \frac{\textit{emp}}{\rho_o \,e_3}\right|_ \textit{surface} 
    1888   % \label{eq:Co_Bq} 
     1888  % \label{eq:DIA_Co_Bq} 
    18891889\] 
    18901890 
     
    18931893\[ 
    18941894  \partial_t \mathcal{V} = \mathcal{A} \;\frac{\overline{\textit{emp}}}{\rho_o} 
    1895   % \label{eq:V_Bq} 
     1895  % \label{eq:DIA_V_Bq} 
    18961896\] 
    18971897 
     
    19121912\begin{equation} 
    19131913  \mathcal{M}_o = \mathcal{M} + \rho_o \,\eta_s \,\mathcal{A} 
    1914   \label{eq:M_Bq} 
     1914  \label{eq:DIA_M_Bq} 
    19151915\end{equation} 
    19161916 
     
    19191919Introducing the total density anomaly, $\mathcal{D}= \int_D d_a \,dv$, 
    19201920where $d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO\ (cf. \autoref{subsec:TRA_eos}) 
    1921 in \autoref{eq:M_Bq} leads to a very simple form for the steric height: 
     1921in \autoref{eq:DIA_M_Bq} leads to a very simple form for the steric height: 
    19221922 
    19231923\begin{equation} 
    19241924  \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D} 
    1925   \label{eq:steric_Bq} 
     1925  \label{eq:DIA_steric_Bq} 
    19261926\end{equation} 
    19271927 
     
    19431943(wetting and drying of grid point is not allowed). 
    19441944 
    1945 Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered. 
     1945Third, the discretisation of \autoref{eq:DIA_steric_Bq} depends on the type of free surface which is considered. 
    19461946In the non linear free surface case, \ie\ \np{ln\_linssh}\forcode{=.true.}, it is given by 
    19471947 
    19481948\[ 
    19491949  \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t} e_{2t} e_{3t} }{ \sum_{i,\,j,\,k}       e_{1t} e_{2t} e_{3t} } 
    1950   % \label{eq:discrete_steric_Bq_nfs} 
     1950  % \label{eq:DIA_discrete_steric_Bq_nfs} 
    19511951\] 
    19521952 
     
    19581958  \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 } 
    19591959                  { \sum_{i,\,j,\,k}       e_{1t}e_{2t}e_{3t} + \sum_{i,\,j}       e_{1t}e_{2t} \eta } 
    1960   % \label{eq:discrete_steric_Bq_fs} 
     1960  % \label{eq:DIA_discrete_steric_Bq_fs} 
    19611961\] 
    19621962 
     
    19781978\[ 
    19791979  \eta_s = - \frac{1}{\mathcal{A}} \int_D d_a(T,S_o,p_o) \,dv 
    1980   % \label{eq:thermosteric_Bq} 
     1980  % \label{eq:DIA_thermosteric_Bq} 
    19811981\] 
    19821982 
     
    20142014    \includegraphics[width=\textwidth]{Fig_mask_subasins} 
    20152015    \caption{ 
    2016       \protect\label{fig:mask_subasins} 
     2016      \protect\label{fig:DIA_mask_subasins} 
    20172017      Decomposition of the World Ocean (here ORCA2) into sub-basin used in to 
    20182018      compute the heat and salt transports as well as the meridional stream-function: 
     
    20452045Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean. 
    20462046The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays, 
    2047 the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}). 
     2047the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:DIA_mask_subasins}). 
    20482048 
    20492049%------------------------------------------namptr----------------------------------------- 
     
    20932093\[ 
    20942094  C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}} 
    2095   % \label{eq:CFL} 
     2095  % \label{eq:DIA_CFL} 
    20962096\] 
    20972097 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex

    r11537 r11543  
    2929    {\em 
    3030      Compatibility changes Major simplification has moved many of the options to external domain configuration tools. 
    31       (see \autoref{apdx:DOMAINcfg}) 
     31      (see \autoref{apdx:DOMCFG}) 
    3232    }                                                                                            \\ 
    3333    {\em 3.x} & {\em Rachid Benshila, Gurvan Madec \& S\'{e}bastien Masson} & 
     
    3838\newpage 
    3939 
    40 Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretisation \autoref{chap:STP}, 
     40Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD}, 
    4141we need to choose a grid for spatial discretisation and related numerical algorithms. 
    4242In the present chapter, we provide a general description of the staggered grid used in \NEMO, 
     
    6060    \includegraphics[width=\textwidth]{Fig_cell} 
    6161    \caption{ 
    62       \protect\label{fig:cell} 
     62      \protect\label{fig:DOM_cell} 
    6363      Arrangement of variables. 
    6464      $t$ indicates scalar points where temperature, salinity, density, pressure and 
     
    7676The arrangement of variables is the same in all directions. 
    7777It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in 
    78 the centre of each face of the cells (\autoref{fig:cell}). 
     78the centre of each face of the cells (\autoref{fig:DOM_cell}). 
    7979This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 
    8080\citep{mesinger.arakawa_bk76}. 
     
    8484The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that 
    8585gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    86 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}. 
     86The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}. 
    8787In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 
    8888the grid-point where the scale factors are defined. 
    89 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. 
     89Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. 
    9090As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 
    9191$\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 
     
    9595centred finite difference approximation, not from their analytical expression. 
    9696This preserves the symmetry of the discrete set of equations and therefore satisfies many of 
    97 the continuous properties (see \autoref{apdx:C}). 
     97the continuous properties (see \autoref{apdx:INVARIANTS}). 
    9898A similar, related remark can be made about the domain size: 
    9999when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors 
     
    123123    \end{tabular} 
    124124    \caption{ 
    125       \protect\label{tab:cell} 
     125      \protect\label{tab:DOM_cell} 
    126126      Location of grid-points as a function of integer or integer and a half value of the column, line or level. 
    127127      This indexing is only used for the writing of the semi -discrete equations. 
     
    145145secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 
    146146(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. 
    147 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. 
     147An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. 
    148148%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    149149\begin{figure}[!t] 
     
    151151    \includegraphics[width=\textwidth]{Fig_zgr_e3} 
    152152    \caption{ 
    153       \protect\label{fig:zgr_e3} 
     153      \protect\label{fig:DOM_zgr_e3} 
    154154      Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
    155155      and (b) analytically derived grid-point position and scale factors. 
     
    174174the midpoint between them are: 
    175175\begin{alignat*}{2} 
    176   % \label{eq:di_mi} 
     176  % \label{eq:DOM_di_mi} 
    177177  \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\ 
    178178  \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2 
     
    180180 
    181181Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 
    182 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at a $t$-point has 
     182Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has 
    183183its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $t$-point. 
    184184These operators have the following discrete forms in the curvilinear $s$-coordinates system: 
     
    198198\end{multline*} 
    199199 
    200 Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 
     200Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 
    201201vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and 
    202202its divergence defined at $t$-points: 
     
    255255In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 
    256256$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 
    257 These two properties will be used extensively in the \autoref{apdx:C} to 
     257These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to 
    258258demonstrate integral conservative properties of the discrete formulation chosen. 
    259259 
     
    269269    \includegraphics[width=\textwidth]{Fig_index_hor} 
    270270    \caption{ 
    271       \protect\label{fig:index_hor} 
     271      \protect\label{fig:DOM_index_hor} 
    272272      Horizontal integer indexing used in the \fortran code. 
    273273      The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices 
     
    290290\label{subsec:DOM_Num_Index_hor} 
    291291 
    292 The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}. 
     292The indexing in the horizontal plane has been chosen as shown in \autoref{fig:DOM_index_hor}. 
    293293For an increasing $i$ index ($j$ index), 
    294294the $t$-point and the eastward $u$-point (northward $v$-point) have the same index 
    295 (see the dashed area in \autoref{fig:index_hor}). 
     295(see the dashed area in \autoref{fig:DOM_index_hor}). 
    296296A $t$-point and its nearest north-east $f$-point have the same $i$-and $j$-indices. 
    297297 
     
    306306given in \autoref{subsec:DOM_cell}. 
    307307The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below 
    308 (\autoref{fig:index_vert}). 
     308(\autoref{fig:DOM_index_vert}). 
    309309The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 
    310 the last $t$-level is always outside the ocean domain (\autoref{fig:index_vert}). 
     310the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). 
    311311Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index 
    312312(\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), 
    313313in contrast to the indexing on the horizontal plane where the $t$-point has the same index as 
    314314the nearest velocity points in the positive direction of the respective horizontal axis index 
    315 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). 
     315(compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). 
    316316Since the scale factors are chosen to be strictly positive, 
    317317a \textit{minus sign} is included in the \fortran implementations of 
     
    324324    \includegraphics[width=\textwidth]{Fig_index_vert} 
    325325    \caption{ 
    326       \protect\label{fig:index_vert} 
     326      \protect\label{fig:DOM_index_vert} 
    327327      Vertical integer indexing used in the \fortran code. 
    328328      Note that the $k$-axis is oriented downward. 
     
    363363the model domain itself can be altered by runtime selections. 
    364364The code previously used to perform vertical discretisation has been incorporated into an external tool 
    365 (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMAINcfg}. 
     365(\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 
    366366 
    367367The next subsections summarise the parameter and fields related to the configuration of the whole model domain. 
     
    418418The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to 
    419419the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$, 
    420 evaluated at the values as specified in \autoref{tab:cell} for the respective grid-point position. 
     420evaluated at the values as specified in \autoref{tab:DOM_cell} for the respective grid-point position. 
    421421The calculation of the values of the horizontal scale factor arrays in general additionally involves 
    422422partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$, 
     
    485485    \includegraphics[width=\textwidth]{Fig_z_zps_s_sps} 
    486486    \caption{ 
    487       \protect\label{fig:z_zps_s_sps} 
     487      \protect\label{fig:DOM_z_zps_s_sps} 
    488488      The ocean bottom as seen by the model: 
    489489      (a) $z$-coordinate with full step, 
     
    510510By default a non-linear free surface is used (\np{ln\_linssh} set to \forcode{=.false.} in \nam{dom}): 
    511511the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 
    512 $z(i,j,k,t)$ (\eg\ \autoref{fig:z_zps_s_sps}f). 
     512$z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 
    513513When a linear free surface is assumed (\np{ln\_linssh} set to \forcode{=.true.} in \nam{dom}), 
    514514the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface 
     
    527527\medskip 
    528528The decision on these choices must be made when the \np{cn\_domcfg} file is constructed. 
    529 Three main choices are offered (\autoref{fig:z_zps_s_sps}a-c): 
     529Three main choices are offered (\autoref{fig:DOM_z_zps_s_sps}a-c): 
    530530 
    531531\begin{itemize} 
     
    536536 
    537537Additionally, hybrid combinations of the three main coordinates are available: 
    538 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e). 
     538$s-z$ or $s-zps$ coordinate (\autoref{fig:DOM_z_zps_s_sps}d and \autoref{fig:DOM_z_zps_s_sps}e). 
    539539 
    540540A further choice related to vertical coordinate concerns 
     
    678678\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})] 
    679679{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
    680 \label{sec:DTA_tsd} 
     680\label{sec:DOM_DTA_tsd} 
    681681%-----------------------------------------namtsd------------------------------------------- 
    682682\nlst{namtsd} 
     
    697697  Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 
    698698  The default version sets horizontally uniform T and profiles as used in the GYRE configuration 
    699   (see \autoref{sec:CFG_gyre}). 
     699  (see \autoref{sec:CFGS_gyre}). 
    700700\end{description} 
    701701 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex

    r11537 r11543  
    6565%           Horizontal divergence and relative vorticity 
    6666%-------------------------------------------------------------------------------------------------------------- 
    67 \subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})] 
    68 {Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
     67\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
    6968\label{subsec:DYN_divcur} 
    7069 
    7170The vorticity is defined at an $f$-point (\ie\ corner point) as follows: 
    7271\begin{equation} 
    73   \label{eq:divcur_cur} 
     72  \label{eq:DYN_divcur_cur} 
    7473  \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right] 
    7574      -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) 
     
    7978It is given by: 
    8079\[ 
    81   % \label{eq:divcur_div} 
     80  % \label{eq:DYN_divcur_div} 
    8281  \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
    8382  \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right] 
     
    102101%           Sea Surface Height evolution 
    103102%-------------------------------------------------------------------------------------------------------------- 
    104 \subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})] 
    105 {Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
     103\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
    106104\label{subsec:DYN_sshwzv} 
    107105 
    108106The sea surface height is given by: 
    109107\begin{equation} 
    110   \label{eq:dynspg_ssh} 
     108  \label{eq:DYN_spg_ssh} 
    111109  \begin{aligned} 
    112110    \frac{\partial \eta }{\partial t} 
     
    123121\textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. 
    124122The sea-surface height is evaluated using exactly the same time stepping scheme as 
    125 the tracer equation \autoref{eq:tra_nxt}: 
     123the tracer equation \autoref{eq:TRA_nxt}: 
    126124a leapfrog scheme in combination with an Asselin time filter, 
    127 \ie\ the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity). 
     125\ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity). 
    128126This is of paramount importance. 
    129127Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to 
     
    134132taking into account the change of the thickness of the levels: 
    135133\begin{equation} 
    136   \label{eq:wzv} 
     134  \label{eq:DYN_wzv} 
    137135  \left\{ 
    138136    \begin{aligned} 
     
    148146re-orientated downward. 
    149147\gmcomment{not sure of this...  to be modified with the change in emp setting} 
    150 In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears. 
     148In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears. 
    151149The upper boundary condition applies at a fixed level $z=0$. 
    152150The top vertical velocity is thus equal to the divergence of the barotropic transport 
    153 (\ie\ the first term in the right-hand-side of \autoref{eq:dynspg_ssh}). 
     151(\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}). 
    154152 
    155153Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates, 
     
    157155in the second case, $w$ is the velocity normal to the $s$-surfaces. 
    158156Note also that the $k$-axis is re-orientated downwards in the \fortran code compared to 
    159 the indexing used in the semi-discrete equations such as \autoref{eq:wzv} 
     157the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv} 
    160158(see \autoref{subsec:DOM_Num_Index_vertical}). 
    161159 
     
    183181%        Vorticity term 
    184182% ------------------------------------------------------------------------------------------------------------- 
    185 \subsection[Vorticity term (\textit{dynvor.F90})] 
    186 {Vorticity term (\protect\mdl{dynvor})} 
     183\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} 
    187184\label{subsec:DYN_vor} 
    188185%------------------------------------------nam_dynvor---------------------------------------------------- 
     
    198195horizontal kinetic energy for the planetary vorticity term (MIX scheme); 
    199196or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy 
    200 (EEN scheme) (see \autoref{subsec:C_vorEEN}). 
     197(EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}). 
    201198In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of 
    202199vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{=.true.}). 
     
    206203%                 enstrophy conserving scheme 
    207204%------------------------------------------------------------- 
    208 \subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens=.true.})] 
    209 {Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{=.true.})} 
     205\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens = .true.})]{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 
    210206\label{subsec:DYN_vor_ens} 
    211207 
     
    216212It is given by: 
    217213\begin{equation} 
    218   \label{eq:dynvor_ens} 
     214  \label{eq:DYN_vor_ens} 
    219215  \left\{ 
    220216    \begin{aligned} 
     
    230226%                 energy conserving scheme 
    231227%------------------------------------------------------------- 
    232 \subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene=.true.})] 
    233 {Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{=.true.})} 
     228\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene = .true.})]{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 
    234229\label{subsec:DYN_vor_ene} 
    235230 
     
    237232It is given by: 
    238233\begin{equation} 
    239   \label{eq:dynvor_ene} 
     234  \label{eq:DYN_vor_ene} 
    240235  \left\{ 
    241236    \begin{aligned} 
     
    251246%                 mix energy/enstrophy conserving scheme 
    252247%------------------------------------------------------------- 
    253 \subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix=.true.})] 
    254 {Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{=.true.})} 
     248\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix = .true.})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})} 
    255249\label{subsec:DYN_vor_mix} 
    256250 
    257251For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used. 
    258 It consists of the ENS scheme (\autoref{eq:dynvor_ens}) for the relative vorticity term, 
    259 and of the ENE scheme (\autoref{eq:dynvor_ene}) applied to the planetary vorticity term. 
    260 \[ 
    261   % \label{eq:dynvor_mix} 
     252It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term, 
     253and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term. 
     254\[ 
     255  % \label{eq:DYN_vor_mix} 
    262256  \left\{ { 
    263257      \begin{aligned} 
     
    277271%                 energy and enstrophy conserving scheme 
    278272%------------------------------------------------------------- 
    279 \subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een=.true.})] 
    280 {Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{=.true.})} 
     273\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een = .true.})]{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 
    281274\label{subsec:DYN_vor_een} 
    282275 
     
    297290The idea is to get rid of the double averaging by considering triad combinations of vorticity. 
    298291It is noteworthy that this solution is conceptually quite similar to the one proposed by 
    299 \citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). 
     292\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}). 
    300293 
    301294The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified 
     
    303296First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: 
    304297\[ 
    305   % \label{eq:pot_vor} 
     298  % \label{eq:DYN_pot_vor} 
    306299  q  = \frac{\zeta +f} {e_{3f} } 
    307300\] 
    308 where the relative vorticity is defined by (\autoref{eq:divcur_cur}), 
     301where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}), 
    309302the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: 
    310303\begin{equation} 
    311   \label{eq:een_e3f} 
     304  \label{eq:DYN_een_e3f} 
    312305  e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} 
    313306\end{equation} 
     
    326319% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    327320 
    328 A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. 
     321A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. 
    329322It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks 
    330323(\np{nn\_een\_e3f}\forcode{=1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{=.true.}). 
     
    340333(\autoref{fig:DYN_een_triad}): 
    341334\begin{equation} 
    342   \label{eq:Q_triads} 
     335  \label{eq:DYN_Q_triads} 
    343336  _i^j \mathbb{Q}^{i_p}_{j_p} 
    344337  = \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) 
     
    348341Finally, the vorticity terms are represented as: 
    349342\begin{equation} 
    350   \label{eq:dynvor_een} 
     343  \label{eq:DYN_vor_een} 
    351344  \left\{ { 
    352345      \begin{aligned} 
     
    361354This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. 
    362355It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow 
    363 (\ie\ $\chi$=$0$) (see \autoref{subsec:C_vorEEN}). 
     356(\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}). 
    364357Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of 
    365358the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}. 
     
    371364%           Kinetic Energy Gradient term 
    372365%-------------------------------------------------------------------------------------------------------------- 
    373 \subsection[Kinetic energy gradient term (\textit{dynkeg.F90})] 
    374 {Kinetic energy gradient term (\protect\mdl{dynkeg})} 
     366\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})} 
    375367\label{subsec:DYN_keg} 
    376368 
    377 As demonstrated in \autoref{apdx:C}, 
     369As demonstrated in \autoref{apdx:INVARIANTS}, 
    378370there is a single discrete formulation of the kinetic energy gradient term that, 
    379371together with the formulation chosen for the vertical advection (see below), 
    380372conserves the total kinetic energy: 
    381373\[ 
    382   % \label{eq:dynkeg} 
     374  % \label{eq:DYN_keg} 
    383375  \left\{ 
    384376    \begin{aligned} 
     
    392384%           Vertical advection term 
    393385%-------------------------------------------------------------------------------------------------------------- 
    394 \subsection[Vertical advection term (\textit{dynzad.F90})] 
    395 {Vertical advection term (\protect\mdl{dynzad})} 
     386\subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})} 
    396387\label{subsec:DYN_zad} 
    397388 
     
    400391conserves the total kinetic energy. 
    401392Indeed, the change of KE due to the vertical advection is exactly balanced by 
    402 the change of KE due to the gradient of KE (see \autoref{apdx:C}). 
    403 \[ 
    404   % \label{eq:dynzad} 
     393the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}). 
     394\[ 
     395  % \label{eq:DYN_zad} 
    405396  \left\{ 
    406397    \begin{aligned} 
     
    439430%           Coriolis plus curvature metric terms 
    440431%-------------------------------------------------------------------------------------------------------------- 
    441 \subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})] 
    442 {Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 
     432\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 
    443433\label{subsec:DYN_cor_flux} 
    444434 
     
    447437It is given by: 
    448438\begin{multline*} 
    449   % \label{eq:dyncor_metric} 
     439  % \label{eq:DYN_cor_metric} 
    450440  f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\ 
    451441  \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right] 
     
    453443\end{multline*} 
    454444 
    455 Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to 
     445Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to 
    456446compute the product of the Coriolis parameter and the vorticity. 
    457 However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has exclusively been used to date. 
     447However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date. 
    458448This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity). 
    459449 
     
    461451%           Flux form Advection term 
    462452%-------------------------------------------------------------------------------------------------------------- 
    463 \subsection[Flux form advection term (\textit{dynadv.F90})] 
    464 {Flux form advection term (\protect\mdl{dynadv})} 
     453\subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})} 
    465454\label{subsec:DYN_adv_flux} 
    466455 
    467456The discrete expression of the advection term is given by: 
    468457\[ 
    469   % \label{eq:dynadv} 
     458  % \label{eq:DYN_adv} 
    470459  \left\{ 
    471460    \begin{aligned} 
     
    495484%                 2nd order centred scheme 
    496485%------------------------------------------------------------- 
    497 \subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2=.true.})] 
    498 {CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{=.true.})} 
     486\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2 = .true.})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 
    499487\label{subsec:DYN_adv_cen2} 
    500488 
    501489In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points: 
    502490\begin{equation} 
    503   \label{eq:dynadv_cen2} 
     491  \label{eq:DYN_adv_cen2} 
    504492  \left\{ 
    505493    \begin{aligned} 
     
    519507%                 UBS scheme 
    520508%------------------------------------------------------------- 
    521 \subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs=.true.})] 
    522 {UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{=.true.})} 
     509\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs = .true.})]{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 
    523510\label{subsec:DYN_adv_ubs} 
    524511 
     
    527514For example, the evaluation of $u_T^{ubs} $ is done as follows: 
    528515\begin{equation} 
    529   \label{eq:dynadv_ubs} 
     516  \label{eq:DYN_adv_ubs} 
    530517  u_T^{ubs} =\overline u ^i-\;\frac{1}{6} 
    531518  \begin{cases} 
     
    547534The UBS scheme is not used in all directions. 
    548535In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and 
    549 $u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used. 
     536$u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used. 
    550537UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the 
    551538sentence:Since vertical mixing of momentum is a source term of the TKE equation...  } 
    552539 
    553 For stability reasons, the first term in (\autoref{eq:dynadv_ubs}), 
     540For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}), 
    554541which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time), 
    555542while the second term, which is the diffusion part of the scheme, 
     
    559546Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by 
    560547one coefficient. 
    561 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
     548Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
    562549This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. 
    563550Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. 
     
    573560%           Hydrostatic pressure gradient term 
    574561% ================================================================ 
    575 \section[Hydrostatic pressure gradient (\textit{dynhpg.F90})] 
    576 {Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
     562\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
    577563\label{sec:DYN_hpg} 
    578564%------------------------------------------nam_dynhpg--------------------------------------------------- 
     
    596582%           z-coordinate with full step 
    597583%-------------------------------------------------------------------------------------------------------------- 
    598 \subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco=.true.})] 
    599 {Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{=.true.})} 
     584\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco = .true.})]{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 
    600585\label{subsec:DYN_hpg_zco} 
    601586 
     
    607592for $k=km$ (surface layer, $jk=1$ in the code) 
    608593\begin{equation} 
    609   \label{eq:dynhpg_zco_surf} 
     594  \label{eq:DYN_hpg_zco_surf} 
    610595  \left\{ 
    611596    \begin{aligned} 
     
    620605for $1<k<km$ (interior layer) 
    621606\begin{equation} 
    622   \label{eq:dynhpg_zco} 
     607  \label{eq:DYN_hpg_zco} 
    623608  \left\{ 
    624609    \begin{aligned} 
     
    633618\end{equation} 
    634619 
    635 Note that the $1/2$ factor in (\autoref{eq:dynhpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as 
     620Note that the $1/2$ factor in (\autoref{eq:DYN_hpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as 
    636621the vertical derivative of the scale factor at the surface level ($z=0$). 
    637622Note also that in case of variable volume level (\texttt{vvl?} defined), 
    638 the surface pressure gradient is included in \autoref{eq:dynhpg_zco_surf} and 
    639 \autoref{eq:dynhpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$. 
     623the surface pressure gradient is included in \autoref{eq:DYN_hpg_zco_surf} and 
     624\autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$. 
    640625 
    641626%-------------------------------------------------------------------------------------------------------------- 
    642627%           z-coordinate with partial step 
    643628%-------------------------------------------------------------------------------------------------------------- 
    644 \subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps=.true.})] 
    645 {Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{=.true.})} 
     629\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps = .true.})]{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 
    646630\label{subsec:DYN_hpg_zps} 
    647631 
     
    674658$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{=.true.}) 
    675659\begin{equation} 
    676   \label{eq:dynhpg_sco} 
     660  \label{eq:DYN_hpg_sco} 
    677661  \left\{ 
    678662    \begin{aligned} 
     
    686670 
    687671Where the first term is the pressure gradient along coordinates, 
    688 computed as in \autoref{eq:dynhpg_zco_surf} - \autoref{eq:dynhpg_zco}, 
     672computed as in \autoref{eq:DYN_hpg_zco_surf} - \autoref{eq:DYN_hpg_zco}, 
    689673and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point 
    690674($e_{3w}$). 
     
    698682(\np{ln\_dynhpg\_djc}\forcode{=.true.}) (currently disabled; under development) 
    699683 
    700 Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated 
     684Note that expression \autoref{eq:DYN_hpg_sco} is commonly used when the variable volume formulation is activated 
    701685(\texttt{vvl?}) because in that case, even with a flat bottom, 
    702686the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}. 
     
    712696\subsection{Ice shelf cavity} 
    713697\label{subsec:DYN_hpg_isf} 
     698 
    714699Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and 
    715700the pressure gradient due to the ocean load (\np{ln\_dynhpg\_isf}\forcode{=.true.}).\\ 
     
    722707A detailed description of this method is described in \citet{losch_JGR08}.\\ 
    723708 
    724 The pressure gradient due to ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in 
     709The pressure gradient due to ocean load is computed using the expression \autoref{eq:DYN_hpg_sco} described in 
    725710\autoref{subsec:DYN_hpg_sco}. 
    726711 
     
    728713%           Time-scheme 
    729714%-------------------------------------------------------------------------------------------------------------- 
    730 \subsection[Time-scheme (\forcode{ln_dynhpg_imp={.true.,.false.}})] 
    731 {Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{=.true.,.false.})} 
     715\subsection[Time-scheme (\forcode{ln_dynhpg_imp = .{true,false}.})]{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .\{true,false\}}.)} 
    732716\label{subsec:DYN_hpg_imp} 
    733717 
     
    748732 
    749733\begin{equation} 
    750   \label{eq:dynhpg_lf} 
     734  \label{eq:DYN_hpg_lf} 
    751735  \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 
    752736  -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right] 
     
    755739$\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{=.true.}): 
    756740\begin{equation} 
    757   \label{eq:dynhpg_imp} 
     741  \label{eq:DYN_hpg_imp} 
    758742  \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 
    759743  -\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] 
    760744\end{equation} 
    761745 
    762 The semi-implicit time scheme \autoref{eq:dynhpg_imp} is made possible without 
     746The semi-implicit time scheme \autoref{eq:DYN_hpg_imp} is made possible without 
    763747significant additional computation since the density can be updated to time level $t+\rdt$ before 
    764748computing the horizontal hydrostatic pressure gradient. 
    765749It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using 
    766 \autoref{eq:dynhpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:dynhpg_lf}. 
    767 Note that \autoref{eq:dynhpg_imp} is equivalent to applying a time filter to the pressure gradient to 
     750\autoref{eq:DYN_hpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:DYN_hpg_lf}. 
     751Note that \autoref{eq:DYN_hpg_imp} is equivalent to applying a time filter to the pressure gradient to 
    768752eliminate high frequency IGWs. 
    769 Obviously, when using \autoref{eq:dynhpg_imp}, 
     753Obviously, when using \autoref{eq:DYN_hpg_imp}, 
    770754the doubling of the time-step is achievable only if no other factors control the time-step, 
    771755such as the stability limits associated with advection or diffusion. 
     
    777761The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows: 
    778762\[ 
    779   % \label{eq:rho_flt} 
     763  % \label{eq:DYN_rho_flt} 
    780764  \rho^t = \rho( \widetilde{T},\widetilde {S},z_t) 
    781765  \quad    \text{with}  \quad 
     
    790774% Surface Pressure Gradient 
    791775% ================================================================ 
    792 \section[Surface pressure gradient (\textit{dynspg.F90})] 
    793 {Surface pressure gradient (\protect\mdl{dynspg})} 
     776\section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})} 
    794777\label{sec:DYN_spg} 
    795778%-----------------------------------------nam_dynspg---------------------------------------------------- 
     
    799782 
    800783Options are defined through the \nam{dyn\_spg} namelist variables. 
    801 The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:PE_hor_pg}). 
     784The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:MB_hor_pg}). 
    802785The main distinction is between the fixed volume case (linear free surface) and 
    803786the variable volume case (nonlinear free surface, \texttt{vvl?} is defined). 
    804 In the linear free surface case (\autoref{subsec:PE_free_surface}) 
     787In the linear free surface case (\autoref{subsec:MB_free_surface}) 
    805788the vertical scale factors $e_{3}$ are fixed in time, 
    806 while they are time-dependent in the nonlinear case (\autoref{subsec:PE_free_surface}). 
     789while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}). 
    807790With both linear and nonlinear free surface, external gravity waves are allowed in the equations, 
    808791which imposes a very small time step when an explicit time stepping is used. 
    809792Two methods are proposed to allow a longer time step for the three-dimensional equations: 
    810 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt?}), 
     793the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:MB_flt?}), 
    811794and the split-explicit free surface described below. 
    812795The extra term introduced in the filtered method is calculated implicitly, 
     
    815798 
    816799The form of the surface pressure gradient term depends on how the user wants to 
    817 handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:PE_hor_pg}). 
     800handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:MB_hor_pg}). 
    818801Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx): 
    819802an explicit formulation which requires a small time step; 
     
    829812% Explicit free surface formulation 
    830813%-------------------------------------------------------------------------------------------------------------- 
    831 \subsection[Explicit free surface (\texttt{ln\_dynspg\_exp}\forcode{=.true.})] 
    832 {Explicit free surface (\protect\np{ln\_dynspg\_exp}\forcode{=.true.})} 
     814\subsection[Explicit free surface (\texttt{ln\_dynspg\_exp}\forcode{ = .true.})]{Explicit free surface (\protect\np{ln\_dynspg\_exp}\forcode{ = .true.})} 
    833815\label{subsec:DYN_spg_exp} 
    834816 
     
    839821is thus simply given by : 
    840822\begin{equation} 
    841   \label{eq:dynspg_exp} 
     823  \label{eq:DYN_spg_exp} 
    842824  \left\{ 
    843825    \begin{aligned} 
     
    856838% Split-explict free surface formulation 
    857839%-------------------------------------------------------------------------------------------------------------- 
    858 \subsection[Split-explicit free surface (\texttt{ln\_dynspg\_ts}\forcode{=.true.})] 
    859 {Split-explicit free surface (\protect\np{ln\_dynspg\_ts}\forcode{=.true.})} 
     840\subsection[Split-explicit free surface (\texttt{ln\_dynspg\_ts}\forcode{ = .true.})]{Split-explicit free surface (\protect\np{ln\_dynspg\_ts}\forcode{ = .true.})} 
    860841\label{subsec:DYN_spg_ts} 
    861842%------------------------------------------namsplit----------------------------------------------------------- 
     
    868849The general idea is to solve the free surface equation and the associated barotropic velocity equations with 
    869850a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables 
    870 (\autoref{fig:DYN_dynspg_ts}). 
     851(\autoref{fig:DYN_spg_ts}). 
    871852The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through 
    872853the \np{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$. 
     
    879860The barotropic mode solves the following equations: 
    880861% \begin{subequations} 
    881 %  \label{eq:BT} 
    882 \begin{equation} 
    883   \label{eq:BT_dyn} 
     862%  \label{eq:DYN_BT} 
     863\begin{equation} 
     864  \label{eq:DYN_BT_dyn} 
    884865  \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}= 
    885866  -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h} 
     
    887868\end{equation} 
    888869\[ 
    889   % \label{eq:BT_ssh} 
     870  % \label{eq:DYN_BT_ssh} 
    890871  \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E 
    891872\] 
     
    893874where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes, 
    894875surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. 
    895 The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress 
    896 (see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration. 
     876The third term on the right hand side of \autoref{eq:DYN_BT_dyn} represents the bottom stress 
     877(see section \autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration. 
    897878Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm 
    898879detailed in \citet{shchepetkin.mcwilliams_OM05}. 
     
    906887    \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts} 
    907888    \caption{ 
    908       \protect\label{fig:DYN_dynspg_ts} 
     889      \protect\label{fig:DYN_spg_ts} 
    909890      Schematic of the split-explicit time stepping scheme for the external and internal modes. 
    910891      Time increases to the right. In this particular exemple, 
     
    929910In the default case (\np{ln\_bt\_fw}\forcode{=.true.}), 
    930911the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps 
    931 (\autoref{fig:DYN_dynspg_ts}a). 
     912(\autoref{fig:DYN_spg_ts}a). 
    932913To avoid aliasing of fast barotropic motions into three dimensional equations, 
    933914time filtering is eventually applied on barotropic quantities (\np{ln\_bt\_av}\forcode{=.true.}). 
     
    11011082% Filtered free surface formulation 
    11021083%-------------------------------------------------------------------------------------------------------------- 
    1103 \subsection[Filtered free surface (\texttt{dynspg\_flt?})] 
    1104 {Filtered free surface (\protect\texttt{dynspg\_flt?})} 
     1084\subsection[Filtered free surface (\texttt{dynspg\_flt?})]{Filtered free surface (\protect\texttt{dynspg\_flt?})} 
    11051085\label{subsec:DYN_spg_fltp} 
    11061086 
    11071087The filtered formulation follows the \citet{roullet.madec_JGR00} implementation. 
    1108 The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly. 
     1088The extra term introduced in the equations (see \autoref{subsec:MB_free_surface}) is solved implicitly. 
    11091089The elliptic solvers available in the code are documented in \autoref{chap:MISC}. 
    11101090 
     
    11121092\gmcomment{               %%% copy from chap-model basics 
    11131093  \[ 
    1114     % \label{eq:spg_flt} 
     1094    % \label{eq:DYN_spg_flt} 
    11151095    \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}} 
    11161096    - g \nabla \left( \tilde{\rho} \ \eta \right) 
     
    11201100  $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, 
    11211101  and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 
    1122   non-linear and viscous terms in \autoref{eq:PE_dyn}. 
     1102  non-linear and viscous terms in \autoref{eq:MB_dyn}. 
    11231103}   %end gmcomment 
    11241104 
     
    11301110% Lateral diffusion term 
    11311111% ================================================================ 
    1132 \section[Lateral diffusion term and operators (\textit{dynldf.F90})] 
    1133 {Lateral diffusion term and operators (\protect\mdl{dynldf})} 
     1112\section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})} 
    11341113\label{sec:DYN_ldf} 
    11351114%------------------------------------------nam_dynldf---------------------------------------------------- 
     
    11451124\ie\ the velocity appearing in its expression is the \textit{before} velocity in time, 
    11461125except for the pure vertical component that appears when a tensor of rotation is used. 
    1147 This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}). 
     1126This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}). 
    11481127 
    11491128At the lateral boundaries either free slip, 
     
    11651144 
    11661145% ================================================================ 
    1167 \subsection[Iso-level laplacian (\forcode{ln_dynldf_lap=.true.})] 
    1168 {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{=.true.})} 
     1146\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap = .true.})]{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 
    11691147\label{subsec:DYN_ldf_lap} 
    11701148 
    11711149For lateral iso-level diffusion, the discrete operator is: 
    11721150\begin{equation} 
    1173   \label{eq:dynldf_lap} 
     1151  \label{eq:DYN_ldf_lap} 
    11741152  \left\{ 
    11751153    \begin{aligned} 
     
    11841162\end{equation} 
    11851163 
    1186 As explained in \autoref{subsec:PE_ldf}, 
     1164As explained in \autoref{subsec:MB_ldf}, 
    11871165this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and 
    11881166ensures a complete separation between the vorticity and divergence parts of the momentum diffusion. 
     
    11911169%           Rotated laplacian operator 
    11921170%-------------------------------------------------------------------------------------------------------------- 
    1193 \subsection[Rotated laplacian (\forcode{ln_dynldf_iso=.true.})] 
    1194 {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{=.true.})} 
     1171\subsection[Rotated laplacian (\forcode{ln_dynldf_iso = .true.})]{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 
    11951172\label{subsec:DYN_ldf_iso} 
    11961173 
     
    12061183The resulting discrete representation is: 
    12071184\begin{equation} 
    1208   \label{eq:dyn_ldf_iso} 
     1185  \label{eq:DYN_ldf_iso} 
    12091186  \begin{split} 
    12101187    D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ 
     
    12501227%           Iso-level bilaplacian operator 
    12511228%-------------------------------------------------------------------------------------------------------------- 
    1252 \subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap=.true.})] 
    1253 {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{=.true.})} 
     1229\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap = .true.})]{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 
    12541230\label{subsec:DYN_ldf_bilap} 
    12551231 
    1256 The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:dynldf_lap} twice. 
     1232The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice. 
    12571233It requires an additional assumption on boundary conditions: 
    12581234the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen, 
     
    12651241%           Vertical diffusion term 
    12661242% ================================================================ 
    1267 \section[Vertical diffusion term (\textit{dynzdf.F90})] 
    1268 {Vertical diffusion term (\protect\mdl{dynzdf})} 
     1243\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})} 
    12691244\label{sec:DYN_zdf} 
    12701245%----------------------------------------------namzdf------------------------------------------------------ 
     
    12791254(\np{ln\_zdfexp}\forcode{=.true.}) using a time splitting technique (\np{nn\_zdfexp} $>$ 1) or 
    12801255$(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}\forcode{=.false.}) 
    1281 (see \autoref{chap:STP}). 
     1256(see \autoref{chap:TD}). 
    12821257Note that namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 
    12831258 
    12841259The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is. 
    1285 The vertical diffusion operators given by \autoref{eq:PE_zdf} take the following semi-discrete space form: 
    1286 \[ 
    1287   % \label{eq:dynzdf} 
     1260The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form: 
     1261\[ 
     1262  % \label{eq:DYN_zdf} 
    12881263  \left\{ 
    12891264    \begin{aligned} 
     
    13031278the vertical turbulent momentum fluxes, 
    13041279\begin{equation} 
    1305   \label{eq:dynzdf_sbc} 
     1280  \label{eq:DYN_zdf_sbc} 
    13061281  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 
    13071282  = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v } 
     
    13161291 
    13171292The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation 
    1318 (see \autoref{sec:ZDF_bfr}) 
     1293(see \autoref{sec:ZDF_drg}) 
    13191294 
    13201295% ================================================================ 
     
    13471322\section{Wetting and drying } 
    13481323\label{sec:DYN_wetdry} 
     1324 
    13491325There are two main options for wetting and drying code (wd): 
    13501326(a) an iterative limiter (il) and (b) a directional limiter (dl). 
     
    13951371%   Iterative limiters 
    13961372%----------------------------------------------------------------------------------------- 
    1397 \subsection[Directional limiter (\textit{wet\_dry.F90})] 
    1398 {Directional limiter (\mdl{wet\_dry})} 
     1373\subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})} 
    13991374\label{subsec:DYN_wd_directional_limiter} 
     1375 
    14001376The principal idea of the directional limiter is that 
    14011377water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn\_wdmin1}). 
     
    14351411%----------------------------------------------------------------------------------------- 
    14361412 
    1437 \subsection[Iterative limiter (\textit{wet\_dry.F90})] 
    1438 {Iterative limiter (\mdl{wet\_dry})} 
     1413\subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})} 
    14391414\label{subsec:DYN_wd_iterative_limiter} 
    14401415 
    1441 \subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})] 
    1442 {Iterative flux limiter (\mdl{wet\_dry})} 
    1443 \label{subsubsec:DYN_wd_il_spg_limiter} 
     1416\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})} 
     1417\label{subsec:DYN_wd_il_spg_limiter} 
    14441418 
    14451419The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry'' 
     
    14491423 
    14501424The continuity equation for the total water depth in a column 
    1451 \begin{equation} \label{dyn_wd_continuity} 
    1452  \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 . 
     1425\begin{equation} 
     1426  \label{eq:DYN_wd_continuity} 
     1427  \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 . 
    14531428\end{equation} 
    14541429can be written in discrete form  as 
    14551430 
    1456 \begin{align} \label{dyn_wd_continuity_2} 
    1457 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 
    1458 &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\ 
    1459 &= \mathrm{zzflx}_{i,j} . 
     1431\begin{align} 
     1432  \label{eq:DYN_wd_continuity_2} 
     1433  \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 
     1434  &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\ 
     1435  &= \mathrm{zzflx}_{i,j} . 
    14601436\end{align} 
    14611437 
     
    14701446(zzflxp) and fluxes that are into the cell (zzflxn).  Clearly 
    14711447 
    1472 \begin{equation} \label{dyn_wd_zzflx_p_n_1} 
    1473 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} . 
     1448\begin{equation} 
     1449  \label{eq:DYN_wd_zzflx_p_n_1} 
     1450  \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} . 
    14741451\end{equation} 
    14751452 
     
    14821459$\mathrm{zcoef}_{i,j}^{(m)}$ such that: 
    14831460 
    1484 \begin{equation} \label{dyn_wd_continuity_coef} 
    1485 \begin{split} 
    1486 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\ 
    1487 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j} 
    1488 \end{split} 
     1461\begin{equation} 
     1462  \label{eq:DYN_wd_continuity_coef} 
     1463  \begin{split} 
     1464    \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\ 
     1465    \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j} 
     1466  \end{split} 
    14891467\end{equation} 
    14901468 
     
    14941472The iteration is initialised by setting 
    14951473 
    1496 \begin{equation} \label{dyn_wd_zzflx_initial} 
    1497 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 
     1474\begin{equation} 
     1475  \label{eq:DYN_wd_zzflx_initial} 
     1476  \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 
    14981477\end{equation} 
    14991478 
    15001479The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the 
    15011480cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell 
    1502 times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this 
     1481times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this 
    15031482condition is 
    15041483 
    1505 \begin{equation} \label{dyn_wd_continuity_if} 
    1506 h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) . 
    1507 \end{equation} 
    1508  
    1509 Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum 
     1484\begin{equation} 
     1485  \label{eq:DYN_wd_continuity_if} 
     1486  h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) . 
     1487\end{equation} 
     1488 
     1489Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum 
    15101490outward flux that can be allowed and still maintain the minimum wet depth: 
    15111491 
    1512 \begin{equation} \label{dyn_wd_max_flux} 
    1513 \begin{split} 
    1514 \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 
    1515 \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] 
    1516 \end{split} 
     1492\begin{equation} 
     1493  \label{eq:DYN_wd_max_flux} 
     1494  \begin{split} 
     1495    \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 
     1496    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] 
     1497  \end{split} 
    15171498\end{equation} 
    15181499 
    15191500Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is 
    1520 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an 
     1501this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an 
    15211502expression for the coefficient needed to multiply the outward flux at this cell in order 
    15221503to avoid drying. 
    15231504 
    1524 \begin{equation} \label{dyn_wd_continuity_nxtcoef} 
    1525 \begin{split} 
    1526 \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 
    1527 \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 
    1528 \end{split} 
     1505\begin{equation} 
     1506  \label{eq:DYN_wd_continuity_nxtcoef} 
     1507  \begin{split} 
     1508    \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 
     1509    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 
     1510  \end{split} 
    15291511\end{equation} 
    15301512 
     
    15451527%      Surface pressure gradients 
    15461528%---------------------------------------------------------------------------------------- 
    1547 \subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})] 
    1548 {Modification of surface pressure gradients (\mdl{dynhpg})} 
    1549 \label{subsubsec:DYN_wd_il_spg} 
     1529\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})} 
     1530\label{subsec:DYN_wd_il_spg} 
    15501531 
    15511532At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the 
     
    15601541neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals 
    15611542variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid 
    1562 column.  The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}. 
     1543column.  The three possible combinations are illustrated in figure \autoref{fig:DYN_WAD_dynhpg}. 
    15631544 
    15641545%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    15651546\begin{figure}[!ht] \begin{center} 
    15661547\includegraphics[width=\textwidth]{Fig_WAD_dynhpg} 
    1567 \caption{ \label{Fig_WAD_dynhpg} 
    1568 Illustrations of the three possible combinations of the logical variables controlling the 
    1569 limiting of the horizontal pressure gradient in wetting and drying regimes} 
    1570 \end{center}\end{figure} 
     1548\caption{ 
     1549  \label{fig:DYN_WAD_dynhpg} 
     1550  Illustrations of the three possible combinations of the logical variables controlling the 
     1551  limiting of the horizontal pressure gradient in wetting and drying regimes} 
     1552\end{center} 
     1553\end{figure} 
    15711554%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    15721555 
     
    15761559of the topography at the two points: 
    15771560 
    1578 \begin{equation} \label{dyn_ll_tmp1} 
    1579 \begin{split} 
    1580 \mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\ 
     1561\begin{equation} 
     1562  \label{eq:DYN_ll_tmp1} 
     1563  \begin{split} 
     1564    \mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\ 
    15811565                     & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\ 
    1582 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\ 
    1583 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\ 
    1584 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 } 
    1585 \end{split} 
     1566                     & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\ 
     1567                     & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\ 
     1568                     & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 } 
     1569  \end{split} 
    15861570\end{equation} 
    15871571 
     
    15901574at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ 
    15911575 
    1592 \begin{equation} \label{dyn_ll_tmp2} 
    1593 \begin{split} 
    1594 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\ 
    1595 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\ 
    1596 & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) . 
    1597 \end{split} 
     1576\begin{equation} 
     1577  \label{eq:DYN_ll_tmp2} 
     1578  \begin{split} 
     1579    \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\ 
     1580    & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\ 
     1581    & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) . 
     1582  \end{split} 
    15981583\end{equation} 
    15991584 
     
    16111596conditions. 
    16121597 
    1613 \subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})] 
    1614 {Additional considerations (\mdl{usrdef\_zgr})} 
    1615 \label{subsubsec:WAD_additional} 
     1598\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})} 
     1599\label{subsec:DYN_WAD_additional} 
    16161600 
    16171601In the very shallow water where wetting and drying occurs the parametrisation of 
     
    16261610%      The WAD test cases 
    16271611%---------------------------------------------------------------------------------------- 
    1628 \subsection[The WAD test cases (\textit{usrdef\_zgr.F90})] 
    1629 {The WAD test cases (\mdl{usrdef\_zgr})} 
    1630 \label{WAD_test_cases} 
     1612\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})} 
     1613\label{subsec:DYN_WAD_test_cases} 
    16311614 
    16321615See the WAD tests MY\_DOC documention for details of the WAD test cases. 
     
    16371620% Time evolution term 
    16381621% ================================================================ 
    1639 \section[Time evolution term (\textit{dynnxt.F90})] 
    1640 {Time evolution term (\protect\mdl{dynnxt})} 
     1622\section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})} 
    16411623\label{sec:DYN_nxt} 
    16421624 
     
    16481630Options are defined through the \nam{dom} namelist variables. 
    16491631The general framework for dynamics time stepping is a leap-frog scheme, 
    1650 \ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:STP}). 
     1632\ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}). 
    16511633The scheme is applied to the velocity, except when 
    16521634using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux}) 
    16531635in the variable volume case (\texttt{vvl?} defined), 
    1654 where it has to be applied to the thickness weighted velocity (see \autoref{sec:A_momentum}) 
     1636where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum}) 
    16551637 
    16561638$\bullet$ vector invariant form or linear free surface 
    16571639(\np{ln\_dynhpg\_vec}\forcode{=.true.} ; \texttt{vvl?} not defined): 
    16581640\[ 
    1659   % \label{eq:dynnxt_vec} 
     1641  % \label{eq:DYN_nxt_vec} 
    16601642  \left\{ 
    16611643    \begin{aligned} 
     
    16691651(\np{ln\_dynhpg\_vec}\forcode{=.false.} ; \texttt{vvl?} defined): 
    16701652\[ 
    1671   % \label{eq:dynnxt_flux} 
     1653  % \label{eq:DYN_nxt_flux} 
    16721654  \left\{ 
    16731655    \begin{aligned} 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_LBC.tex

    r11537 r11543  
    5757 
    5858\[ 
    59   % \label{eq:lbc_aaaa} 
     59  % \label{eq:LBC_aaaa} 
    6060  \frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT} 
    6161  }{e_{1u} } \; \delta_{i+1 / 2} \left[ T \right]\;\;mask_u 
     
    134134  the no-slip boundary condition, simply by multiplying it by the mask$_{f}$ : 
    135135  \[ 
    136     % \label{eq:lbc_bbbb} 
     136    % \label{eq:LBC_bbbb} 
    137137    \zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta_{i+1/2} 
    138138        \left[ {e_{2v} \,v} \right]-\delta_{j+1/2} \left[ {e_{1u} \,u} \right]} 
     
    226226The north fold boundary condition has been introduced in order to handle the north boundary of 
    227227a three-polar ORCA grid. 
    228 Such a grid has two poles in the northern hemisphere (\autoref{fig:MISC_ORCA_msh}, 
    229 and thus requires a specific treatment illustrated in \autoref{fig:North_Fold_T}. 
     228Such a grid has two poles in the northern hemisphere (\autoref{fig:CFGS_ORCA_msh}, 
     229and thus requires a specific treatment illustrated in \autoref{fig:LBC_North_Fold_T}. 
    230230Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition. 
    231231 
     
    235235    \includegraphics[width=\textwidth]{Fig_North_Fold_T} 
    236236    \caption{ 
    237       \protect\label{fig:North_Fold_T} 
     237      \protect\label{fig:LBC_North_Fold_T} 
    238238      North fold boundary with a $T$-point pivot and cyclic east-west boundary condition ($jperio=4$), 
    239239      as used in ORCA 2, 1/4, and 1/12. 
     
    256256%----------------------------------------------------------------------------------------------- 
    257257 
    258 For massively parallel processing (mpp), a domain decomposition method is used. The basic idea of the method is to split the large computation domain of a numerical experiment into several smaller domains and solve the set of equations by addressing independent local problems. Each processor has its own local memory and computes the model equation over a subdomain of the whole model domain. The subdomain boundary conditions are specified through communications between processors which are organized by explicit statements (message passing method). The present implementation is largely inspired by Guyon's work [Guyon 1995]. 
     258For massively parallel processing (mpp), a domain decomposition method is used. 
     259The basic idea of the method is to split the large computation domain of a numerical experiment into several smaller domains and 
     260solve the set of equations by addressing independent local problems. 
     261Each processor has its own local memory and computes the model equation over a subdomain of the whole model domain. 
     262The subdomain boundary conditions are specified through communications between processors which are organized by 
     263explicit statements (message passing method). 
     264The present implementation is largely inspired by Guyon's work [Guyon 1995]. 
    259265 
    260266The parallelization strategy is defined by the physical characteristics of the ocean model. 
     
    272278each processor sends to its neighbouring processors the update values of the points corresponding to 
    273279the interior overlapping area to its neighbouring sub-domain (\ie\ the innermost of the two overlapping rows). 
    274 Communications are first done according to the east-west direction and next according to the north-south direction. There is no specific communications for the corners. The communication is done through the Message Passing Interface (MPI) and requires \key{mpp\_mpi}. Use also \key{mpi2} if MPI3 is not available on your computer. 
     280Communications are first done according to the east-west direction and next according to the north-south direction. 
     281There is no specific communications for the corners. 
     282The communication is done through the Message Passing Interface (MPI) and requires \key{mpp\_mpi}. 
     283Use also \key{mpi2} if MPI3 is not available on your computer. 
    275284The data exchanges between processors are required at the very place where 
    276285lateral domain boundary conditions are set in the mono-domain computation: 
     
    285294    \includegraphics[width=\textwidth]{Fig_mpp} 
    286295    \caption{ 
    287       \protect\label{fig:mpp} 
     296      \protect\label{fig:LBC_mpp} 
    288297      Positioning of a sub-domain when massively parallel processing is used. 
    289298    } 
     
    292301%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    293302 
    294 In \NEMO, the splitting is regular and arithmetic. The total number of subdomains corresponds to the number of MPI processes allocated to \NEMO\ when the model is launched (\ie\ mpirun -np x ./nemo will automatically give x subdomains). The i-axis is divided by \np{jpni} and the j-axis by \np{jpnj}. These parameters are defined in \nam{mpp} namelist. If \np{jpni} and \np{jpnj} are < 1, they will be automatically redefined in the code to give the best domain decomposition (see bellow). 
    295  
    296 Each processor is independent and without message passing or synchronous process, programs run alone and access just its own local memory. For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil) that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. 
    297 These dimensions include the internal domain and the overlapping rows. The number of rows to exchange (known as the halo) is usually set to one (nn\_hls=1, in \mdl{par\_oce}, and must be kept to one until further notice). The whole domain dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. The relationship between the whole domain and a sub-domain is: 
    298 \[ 
    299   jpi = ( jpiglo-2\times nn\_hls + (jpni-1) ) / jpni + 2\times nn\_hls 
    300 \] 
    301 \[ 
     303In \NEMO, the splitting is regular and arithmetic. 
     304The total number of subdomains corresponds to the number of MPI processes allocated to \NEMO\ when the model is launched 
     305(\ie\ mpirun -np x ./nemo will automatically give x subdomains). 
     306The i-axis is divided by \np{jpni} and the j-axis by \np{jpnj}. 
     307These parameters are defined in \nam{mpp} namelist. 
     308If \np{jpni} and \np{jpnj} are < 1, they will be automatically redefined in the code to give the best domain decomposition 
     309(see bellow). 
     310 
     311Each processor is independent and without message passing or synchronous process, programs run alone and access just its own local memory. 
     312For this reason, 
     313the main model dimensions are now the local dimensions of the subdomain (pencil) that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. 
     314These dimensions include the internal domain and the overlapping rows. 
     315The number of rows to exchange (known as the halo) is usually set to one (nn\_hls=1, in \mdl{par\_oce}, 
     316and must be kept to one until further notice). 
     317The whole domain dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. 
     318The relationship between the whole domain and a sub-domain is: 
     319\begin{gather*} 
     320  jpi = ( jpiglo-2\times nn\_hls + (jpni-1) ) / jpni + 2\times nn\_hls \\ 
    302321  jpj = ( jpjglo-2\times nn\_hls + (jpnj-1) ) / jpnj + 2\times nn\_hls 
    303 \] 
    304  
    305 One also defines variables nldi and nlei which correspond to the internal domain bounds, and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain (\autoref{fig:mpp}). Note that since the version 4, there is no more extra-halo area as defined in \autoref{fig:mpp} so \jp{jpi} is now always equal to nlci and \jp{jpj} equal to nlcj. 
     322\end{gather*} 
     323 
     324One also defines variables nldi and nlei which correspond to the internal domain bounds, and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain (\autoref{fig:LBC_mpp}). Note that since the version 4, there is no more extra-halo area as defined in \autoref{fig:LBC_mpp} so \jp{jpi} is now always equal to nlci and \jp{jpj} equal to nlcj. 
    306325 
    307326An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$, 
    308327a global array (whole domain) by the relationship: 
    309328\[ 
    310   % \label{eq:lbc_nimpp} 
     329  % \label{eq:LBC_nimpp} 
    311330  T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k), 
    312331\] 
     
    322341 
    323342If the domain decomposition is automatically defined (when \np{jpni} and \np{jpnj} are < 1), the decomposition chosen by the model will minimise the sub-domain size (defined as $max_{all domains}(jpi \times jpj)$) and maximize the number of eliminated land subdomains. This means that no other domain decomposition (a set of \np{jpni} and \np{jpnj} values) will use less processes than $(jpni  \times  jpnj - N_{land})$ and get a smaller subdomain size. 
    324 In order to specify $N_{mpi}$ properly (minimize $N_{useless}$), you must run the model once with \np{ln\_list} activated. In this case, the model will start the initialisation phase, print the list of optimum decompositions ($N_{mpi}$, \np{jpni} and \np{jpnj}) in \texttt{ocean.output} and directly abort. The maximum value of $N_{mpi}$ tested in this list is given by $max(N_{MPI\_tasks}, \np{jpni} \times \np{jpnj})$. For example, run the model on 40 nodes with ln\_list activated and $\np{jpni} = 10000$ and $\np{jpnj} = 1$, will print the list of optimum domains decomposition from 1 to about 10000.  
    325  
    326 Processors are numbered from 0 to $N_{mpi} - 1$. Subdomains containning some ocean points are numbered first from 0 to $jpni * jpnj - N_{land} -1$. The remaining $N_{useless}$ land subdomains are numbered next, which means that, for a given (\np{jpni}, \np{jpnj}), the numbers attributed to he ocean subdomains do not vary with $N_{useless}$.  
     343In order to specify $N_{mpi}$ properly (minimize $N_{useless}$), you must run the model once with \np{ln\_list} activated. In this case, the model will start the initialisation phase, print the list of optimum decompositions ($N_{mpi}$, \np{jpni} and \np{jpnj}) in \texttt{ocean.output} and directly abort. The maximum value of $N_{mpi}$ tested in this list is given by $max(N_{MPI\_tasks}, \np{jpni} \times \np{jpnj})$. For example, run the model on 40 nodes with ln\_list activated and $\np{jpni} = 10000$ and $\np{jpnj} = 1$, will print the list of optimum domains decomposition from 1 to about 10000. 
     344 
     345Processors are numbered from 0 to $N_{mpi} - 1$. Subdomains containning some ocean points are numbered first from 0 to $jpni * jpnj - N_{land} -1$. The remaining $N_{useless}$ land subdomains are numbered next, which means that, for a given (\np{jpni}, \np{jpnj}), the numbers attributed to he ocean subdomains do not vary with $N_{useless}$. 
    327346 
    328347When land processors are eliminated, the value corresponding to these locations in the model output files is undefined. \np{ln\_mskland} must be activated in order avoid Not a Number values in output files. Note that it is better to not eliminate land processors when creating a meshmask file (\ie\ when setting a non-zero value to \np{nn\_msh}). 
     
    332351  \begin{center} 
    333352    \includegraphics[width=\textwidth]{Fig_mppini2} 
    334     \caption { 
    335       \protect\label{fig:mppini2} 
     353    \caption[Atlantic domain]{ 
     354      \protect\label{fig:LBC_mppini2} 
    336355      Example of Atlantic domain defined for the CLIPPER projet. 
    337356      Initial grid is composed of 773 x 1236 horizontal points. 
     
    374393%---------------------------------------------- 
    375394\subsection{Namelists} 
    376 \label{subsec:BDY_namelist} 
     395\label{subsec:LBC_bdy_namelist} 
    377396 
    378397The BDY module is activated by setting \np{ln\_bdy}\forcode{=.true.} . 
     
    384403In the example above, there are two boundary sets, the first of which is defined via a file and 
    385404the second is defined in the namelist. 
    386 For more details of the definition of the boundary geometry see section \autoref{subsec:BDY_geometry}. 
     405For more details of the definition of the boundary geometry see section \autoref{subsec:LBC_bdy_geometry}. 
    387406 
    388407For each boundary set a boundary condition has to be chosen for the barotropic solution 
     
    441460%---------------------------------------------- 
    442461\subsection{Flow relaxation scheme} 
    443 \label{subsec:BDY_FRS_scheme} 
     462\label{subsec:LBC_bdy_FRS_scheme} 
    444463 
    445464The Flow Relaxation Scheme (FRS) \citep{davies_QJRMS76,engedahl_T95}, 
     
    448467Given a model prognostic variable $\Phi$ 
    449468\[ 
    450   % \label{eq:bdy_frs1} 
     469  % \label{eq:LBC_bdy_frs1} 
    451470  \Phi(d) = \alpha(d)\Phi_{e}(d) + (1-\alpha(d))\Phi_{m}(d)\;\;\;\;\; d=1,N 
    452471\] 
     
    457476the prognostic equation for $\Phi$ of the form: 
    458477\[ 
    459   % \label{eq:bdy_frs2} 
     478  % \label{eq:LBC_bdy_frs2} 
    460479  -\frac{1}{\tau}\left(\Phi - \Phi_{e}\right) 
    461480\] 
    462481where the relaxation time scale $\tau$ is given by a function of $\alpha$ and the model time step $\Delta t$: 
    463482\[ 
    464   % \label{eq:bdy_frs3} 
     483  % \label{eq:LBC_bdy_frs3} 
    465484  \tau = \frac{1-\alpha}{\alpha}  \,\rdt 
    466485\] 
     
    472491The function $\alpha$ is specified as a $tanh$ function: 
    473492\[ 
    474   % \label{eq:bdy_frs4} 
     493  % \label{eq:LBC_bdy_frs4} 
    475494  \alpha(d) = 1 - \tanh\left(\frac{d-1}{2}\right),       \quad d=1,N 
    476495\] 
     
    480499%---------------------------------------------- 
    481500\subsection{Flather radiation scheme} 
    482 \label{subsec:BDY_flather_scheme} 
     501\label{subsec:LBC_bdy_flather_scheme} 
    483502 
    484503The \citet{flather_JPO94} scheme is a radiation condition on the normal, 
    485504depth-mean transport across the open boundary. 
    486505It takes the form 
    487 \begin{equation}  \label{eq:bdy_fla1} 
    488 U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right), 
     506\begin{equation} 
     507  \label{eq:LBC_bdy_fla1} 
     508  U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right), 
    489509\end{equation} 
    490510where $U$ is the depth-mean velocity normal to the boundary and $\eta$ is the sea surface height, 
     
    495515the external depth-mean normal velocity, 
    496516plus a correction term that allows gravity waves generated internally to exit the model boundary. 
    497 Note that the sea-surface height gradient in \autoref{eq:bdy_fla1} is a spatial gradient across the model boundary, 
     517Note that the sea-surface height gradient in \autoref{eq:LBC_bdy_fla1} is a spatial gradient across the model boundary, 
    498518so that $\eta_{e}$ is defined on the $T$ points with $nbr=1$ and $\eta$ is defined on the $T$ points with $nbr=2$. 
    499519$U$ and $U_{e}$ are defined on the $U$ or $V$ points with $nbr=1$, \ie\ between the two $T$ grid points. 
     
    501521%---------------------------------------------- 
    502522\subsection{Orlanski radiation scheme} 
    503 \label{subsec:BDY_orlanski_scheme} 
     523\label{subsec:LBC_bdy_orlanski_scheme} 
    504524 
    505525The Orlanski scheme is based on the algorithm described by \citep{marchesiello.mcwilliams.ea_OM01}, hereafter MMS. 
     
    507527The adaptive Orlanski condition solves a wave plus relaxation equation at the boundary: 
    508528\begin{equation} 
    509 \frac{\partial\phi}{\partial t} + c_x \frac{\partial\phi}{\partial x} + c_y \frac{\partial\phi}{\partial y} = 
    510                                                 -\frac{1}{\tau}(\phi - \phi^{ext}) 
    511 \label{eq:wave_continuous} 
     529  \label{eq:LBC_wave_continuous} 
     530  \frac{\partial\phi}{\partial t} + c_x \frac{\partial\phi}{\partial x} + c_y \frac{\partial\phi}{\partial y} = 
     531  -\frac{1}{\tau}(\phi - \phi^{ext}) 
    512532\end{equation} 
    513533 
     
    515535velocities are diagnosed from the model fields as: 
    516536 
    517 \begin{equation} \label{eq:cx} 
    518 c_x = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial x}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2} 
     537\begin{equation} 
     538  \label{eq:LBC_cx} 
     539  c_x = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial x}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2} 
    519540\end{equation} 
    520541\begin{equation} 
    521 \label{eq:cy} 
    522 c_y = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial y}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2} 
     542  \label{eq:LBC_cy} 
     543  c_y = -\frac{\partial\phi}{\partial t}\frac{\partial\phi / \partial y}{(\partial\phi /\partial x)^2 + (\partial\phi /\partial y)^2} 
    523544\end{equation} 
    524545 
    525546(As noted by MMS, this is a circular diagnosis of the phase speeds which only makes sense on a discrete grid). 
    526 Equation (\autoref{eq:wave_continuous}) is defined adaptively depending on the sign of the phase velocity normal to the boundary $c_x$. 
     547Equation (\autoref{eq:LBC_wave_continuous}) is defined adaptively depending on the sign of the phase velocity normal to the boundary $c_x$. 
    527548For $c_x$ outward, we have 
    528549 
     
    534555 
    535556\begin{equation} 
    536 \tau = \tau_{in}\,\,\,;\,\,\, c_x = c_y = 0 
    537 \label{eq:tau_in} 
     557  \label{eq:LBC_tau_in} 
     558  \tau = \tau_{in}\,\,\,;\,\,\, c_x = c_y = 0 
    538559\end{equation} 
    539560 
    540561Generally the relaxation time scale at inward propagation points (\np{rn\_time\_dmp}) is set much shorter than the time scale at outward propagation 
    541562points (\np{rn\_time\_dmp\_out}) so that the solution is constrained more strongly by the external data at inward propagation points. 
    542 See \autoref{subsec:BDY_relaxation} for detailed on the spatial shape of the scaling.\\ 
     563See \autoref{subsec:LBC_bdy_relaxation} for detailed on the spatial shape of the scaling.\\ 
    543564The ``normal propagation of oblique radiation'' or NPO approximation (called \forcode{'orlanski_npo'}) involves assuming 
    544 that $c_y$ is zero in equation (\autoref{eq:wave_continuous}), but including 
    545 this term in the denominator of equation (\autoref{eq:cx}). Both versions of the scheme are options in BDY. Equations 
    546 (\autoref{eq:wave_continuous}) - (\autoref{eq:tau_in}) correspond to equations (13) - (15) and (2) - (3) in MMS.\\ 
     565that $c_y$ is zero in equation (\autoref{eq:LBC_wave_continuous}), but including 
     566this term in the denominator of equation (\autoref{eq:LBC_cx}). Both versions of the scheme are options in BDY. Equations 
     567(\autoref{eq:LBC_wave_continuous}) - (\autoref{eq:LBC_tau_in}) correspond to equations (13) - (15) and (2) - (3) in MMS.\\ 
    547568 
    548569%---------------------------------------------- 
    549570\subsection{Relaxation at the boundary} 
    550 \label{subsec:BDY_relaxation} 
     571\label{subsec:LBC_bdy_relaxation} 
    551572 
    552573In addition to a specific boundary condition specified as \np{cn\_tra} and \np{cn\_dyn3d}, relaxation on baroclinic velocities and tracers variables are available. 
     
    564585%---------------------------------------------- 
    565586\subsection{Boundary geometry} 
    566 \label{subsec:BDY_geometry} 
     587\label{subsec:LBC_bdy_geometry} 
    567588 
    568589Each open boundary set is defined as a list of points. 
     
    615636%---------------------------------------------- 
    616637\subsection{Input boundary data files} 
    617 \label{subsec:BDY_data} 
     638\label{subsec:LBC_bdy_data} 
    618639 
    619640The data files contain the data arrays in the order in which the points are defined in the $nbi$ and $nbj$ arrays. 
     
    655676%---------------------------------------------- 
    656677\subsection{Volume correction} 
    657 \label{subsec:BDY_vol_corr} 
     678\label{subsec:LBC_bdy_vol_corr} 
    658679 
    659680There is an option to force the total volume in the regional model to be constant. 
     
    672693%---------------------------------------------- 
    673694\subsection{Tidal harmonic forcing} 
    674 \label{subsec:BDY_tides} 
     695\label{subsec:LBC_bdy_tides} 
    675696 
    676697%-----------------------------------------nambdy_tide-------------------------------------------- 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex

    r11537 r11543  
    1313\newpage 
    1414 
    15 The lateral physics terms in the momentum and tracer equations have been described in \autoref{eq:PE_zdf} and 
     15The lateral physics terms in the momentum and tracer equations have been described in \autoref{eq:MB_zdf} and 
    1616their discrete formulation in \autoref{sec:TRA_ldf} and \autoref{sec:DYN_ldf}). 
    1717In this section we further discuss each lateral physics option. 
     
    2525Note that this chapter describes the standard implementation of iso-neutral tracer mixing.  
    2626Griffies's implementation, which is used if \np{ln\_traldf\_triad}\forcode{=.true.}, 
    27 is described in \autoref{apdx:triad} 
     27is described in \autoref{apdx:TRIADS} 
    2828 
    2929%-----------------------------------namtra_ldf - namdyn_ldf-------------------------------------------- 
     
    8282the cell of the quantity to be diffused. 
    8383For a tracer, this leads to the following four slopes: 
    84 $r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \autoref{eq:tra_ldf_iso}), 
     84$r_{1u}$, $r_{1w}$, $r_{2v}$, $r_{2w}$ (see \autoref{eq:TRA_ldf_iso}), 
    8585while for momentum the slopes are  $r_{1t}$, $r_{1uw}$, $r_{2f}$, $r_{2uw}$ for $u$ and 
    8686$r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$.  
     
    9292In $s$-coordinates, geopotential mixing (\ie\ horizontal mixing) $r_1$ and $r_2$ are the slopes between 
    9393the geopotential and computational surfaces. 
    94 Their discrete formulation is found by locally solving \autoref{eq:tra_ldf_iso} when 
     94Their discrete formulation is found by locally solving \autoref{eq:TRA_ldf_iso} when 
    9595the diffusive fluxes in the three directions are set to zero and $T$ is assumed to be horizontally uniform, 
    9696\ie\ a linear function of $z_T$, the depth of a $T$-point.  
     
    9898 
    9999\begin{equation} 
    100   \label{eq:ldfslp_geo} 
     100  \label{eq:LDF_slp_geo} 
    101101  \begin{aligned} 
    102102    r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)} 
     
    125125Their discrete formulation is found using the fact that the diffusive fluxes of 
    126126locally referenced potential density (\ie\ $in situ$ density) vanish. 
    127 So, substituting $T$ by $\rho$ in \autoref{eq:tra_ldf_iso} and setting the diffusive fluxes in 
     127So, substituting $T$ by $\rho$ in \autoref{eq:TRA_ldf_iso} and setting the diffusive fluxes in 
    128128the three directions to zero leads to the following definition for the neutral slopes: 
    129129 
    130130\begin{equation} 
    131   \label{eq:ldfslp_iso} 
     131  \label{eq:LDF_slp_iso} 
    132132  \begin{split} 
    133133    r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]} 
     
    145145 
    146146%gm% rewrite this as the explanation is not very clear !!! 
    147 %In practice, \autoref{eq:ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.  
    148  
    149 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters). 
    150  
    151 %In the $z$-coordinate, the derivative of the  \autoref{eq:ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.  
    152  
    153 As the mixing is performed along neutral surfaces, the gradient of $\rho$ in \autoref{eq:ldfslp_iso} has to 
     147%In practice, \autoref{eq:LDF_slp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:LDF_slp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.  
     148 
     149%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:LDF_slp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters). 
     150 
     151%In the $z$-coordinate, the derivative of the  \autoref{eq:LDF_slp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.  
     152 
     153As the mixing is performed along neutral surfaces, the gradient of $\rho$ in \autoref{eq:LDF_slp_iso} has to 
    154154be evaluated at the same local pressure 
    155155(which, in decibars, is approximated by the depth in meters in the model). 
    156 Therefore \autoref{eq:ldfslp_iso} cannot be used as such, 
     156Therefore \autoref{eq:LDF_slp_iso} cannot be used as such, 
    157157but further transformation is needed depending on the vertical coordinate used: 
    158158 
     
    160160 
    161161\item[$z$-coordinate with full step: ] 
    162   in \autoref{eq:ldfslp_iso} the densities appearing in the $i$ and $j$ derivatives  are taken at the same depth, 
     162  in \autoref{eq:LDF_slp_iso} the densities appearing in the $i$ and $j$ derivatives  are taken at the same depth, 
    163163  thus the $in situ$ density can be used. 
    164164  This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, 
     
    173173  in the current release of \NEMO, iso-neutral mixing is only employed for $s$-coordinates if 
    174174  the Griffies scheme is used (\np{ln\_traldf\_triad}\forcode{=.true.}; 
    175   see \autoref{apdx:triad}). 
     175  see \autoref{apdx:TRIADS}). 
    176176  In other words, iso-neutral mixing will only be accurately represented with a linear equation of state 
    177177  (\np{ln\_seos}\forcode{=.true.}). 
    178   In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:ldfslp_iso} 
     178  In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:LDF_slp_iso} 
    179179  will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes. 
    180180 
     
    193193 
    194194\[ 
    195   % \label{eq:ldfslp_iso2} 
     195  % \label{eq:LDF_slp_iso2} 
    196196  \begin{split} 
    197197    r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac 
     
    230230To overcome this problem, several techniques have been proposed in which the numerical schemes of 
    231231the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}. 
    232 Griffies's scheme is now available in \NEMO\ if \np{ln\_traldf\_triad}\forcode{=.true.}; see \autoref{apdx:triad}. 
     232Griffies's scheme is now available in \NEMO\ if \np{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:TRIADS}. 
    233233Here, another strategy is presented \citep{lazar_phd97}: 
    234234a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of 
     
    280280    \includegraphics[width=\textwidth]{Fig_eiv_slp} 
    281281    \caption{ 
    282       \protect\label{fig:eiv_slp} 
     282      \protect\label{fig:LDF_eiv_slp} 
    283283      Vertical profile of the slope used for lateral mixing in the mixed layer: 
    284284      \textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior, 
     
    304304The iso-neutral diffusion operator on momentum is the same as the one used on tracers but 
    305305applied to each component of the velocity separately 
    306 (see \autoref{eq:dyn_ldf_iso} in section~\autoref{subsec:DYN_ldf_iso}). 
     306(see \autoref{eq:DYN_ldf_iso} in section~\autoref{subsec:DYN_ldf_iso}). 
    307307The slopes between the surface along which the diffusion operator acts and the surface of computation 
    308308($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the $u$-component, and $T$-, $f$- and 
    309309\textit{vw}- points for the $v$-component. 
    310310They are computed from the slopes used for tracer diffusion, 
    311 \ie\ \autoref{eq:ldfslp_geo} and \autoref{eq:ldfslp_iso}: 
     311\ie\ \autoref{eq:LDF_slp_geo} and \autoref{eq:LDF_slp_iso}: 
    312312 
    313313\[ 
    314   % \label{eq:ldfslp_dyn} 
     314  % \label{eq:LDF_slp_dyn} 
    315315  \begin{aligned} 
    316316    &r_{1t}\ \ = \overline{r_{1u}}^{\,i}       &&&    r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\ 
     
    371371 
    372372\begin{equation} 
    373   \label{eq:constantah} 
     373  \label{eq:LDF_constantah} 
    374374  A_o^l = \left\{ 
    375375    \begin{aligned} 
     
    386386 
    387387In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which 
    388 the surface value is given by \autoref{eq:constantah}, the bottom value is 1/4 of the surface value, 
     388the surface value is given by \autoref{eq:LDF_constantah}, the bottom value is 1/4 of the surface value, 
    389389and the transition takes place around z=500~m with a width of 200~m. 
    390390This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users. 
     
    396396the type of operator used: 
    397397\begin{equation} 
    398   \label{eq:title} 
     398  \label{eq:LDF_title} 
    399399  A_l = \left\{ 
    400400    \begin{aligned} 
     
    411411model configurations presenting large changes in grid spacing such as global ocean models. 
    412412Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to 
    413 large coefficient compare to the smallest grid size (see \autoref{sec:STP_forward_imp}), 
     413large coefficient compare to the smallest grid size (see \autoref{sec:TD_forward_imp}), 
    414414especially when using a bilaplacian operator. 
    415415 
     
    429429 
    430430\begin{equation} 
    431   \label{eq:flowah} 
     431  \label{eq:LDF_flowah} 
    432432  A_l = \left\{ 
    433433    \begin{aligned} 
     
    445445 
    446446\begin{equation} 
    447   \label{eq:smag1} 
     447  \label{eq:LDF_smag1} 
    448448  \begin{split} 
    449449    T_{smag}^{-1} & = \sqrt{\left( \partial_x u - \partial_y v\right)^2 + \left( \partial_y u + \partial_x v\right)^2  } \\ 
     
    455455 
    456456\begin{equation} 
    457   \label{eq:smag2} 
     457  \label{eq:LDF_smag2} 
    458458  A_{smag} = \left\{ 
    459459    \begin{aligned} 
     
    464464\end{equation} 
    465465 
    466 For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:STP_forward_imp}) so that: 
    467 \begin{equation} 
    468   \label{eq:smag3} 
     466For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:TD_forward_imp}) so that: 
     467\begin{equation} 
     468  \label{eq:LDF_smag3} 
    469469    \begin{aligned} 
    470470      & C_{min} \frac{1}{2}   \lvert U \rvert  e    < A_{smag} < C_{max} \frac{e^2}{   8\rdt}                 & \text{for laplacian operator } \\ 
     
    480480 
    481481(1) the momentum diffusion operator acting along model level surfaces is written in terms of curl and 
    482 divergent components of the horizontal current (see \autoref{subsec:PE_ldf}). 
     482divergent components of the horizontal current (see \autoref{subsec:MB_ldf}). 
    483483Although the eddy coefficient could be set to different values in these two terms, 
    484484this option is not currently available.  
     
    486486(2) with an horizontally varying viscosity, the quadratic integral constraints on enstrophy and on the square of 
    487487the horizontal divergence for operators acting along model-surfaces are no longer satisfied 
    488 (\autoref{sec:dynldf_properties}). 
     488(\autoref{sec:INVARIANTS_dynldf_properties}). 
    489489 
    490490% ================================================================ 
     
    527527the formulation of which depends on the slopes of iso-neutral surfaces. 
    528528Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, 
    529 \ie\ \autoref{eq:ldfslp_geo} is used in $z$-coordinates, 
    530 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $s$-coordinates. 
     529\ie\ \autoref{eq:LDF_slp_geo} is used in $z$-coordinates, 
     530and the sum \autoref{eq:LDF_slp_geo} + \autoref{eq:LDF_slp_iso} in $s$-coordinates. 
    531531 
    532532If isopycnal mixing is used in the standard way, \ie\ \np{ln\_traldf\_triad}\forcode{=.false.}, the eddy induced velocity is given by:  
    533533\begin{equation} 
    534   \label{eq:ldfeiv} 
     534  \label{eq:LDF_eiv} 
    535535  \begin{split} 
    536536    u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\ 
     
    554554\colorbox{yellow}{CASE \np{nn\_aei\_ijk\_t} = 21 to be added} 
    555555 
    556 In case of setting \np{ln\_traldf\_triad}\forcode{=.true.}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:triad}. 
     556In case of setting \np{ln\_traldf\_triad}\forcode{ = .true.}, a skew form of the eddy induced advective fluxes is used, which is described in \autoref{apdx:TRIADS}. 
    557557 
    558558% ================================================================ 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex

    r11435 r11543  
    691691 
    692692Examples of the weights calculated for an observation with rectangular and radial footprints are shown in 
    693 \autoref{fig:obsavgrec} and~\autoref{fig:obsavgrad}. 
     693\autoref{fig:OBS_avgrec} and~\autoref{fig:OBS_avgrad}. 
    694694 
    695695%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    698698    \includegraphics[width=\textwidth]{Fig_OBS_avg_rec} 
    699699    \caption{ 
    700       \protect\label{fig:obsavgrec} 
     700      \protect\label{fig:OBS_avgrec} 
    701701      Weights associated with each model grid box (blue lines and numbers) 
    702702      for an observation at -170.5\deg{E}, 56.0\deg{N} with a rectangular footprint of 1\deg x 1\deg. 
     
    711711    \includegraphics[width=\textwidth]{Fig_OBS_avg_rad} 
    712712    \caption{ 
    713       \protect\label{fig:obsavgrad} 
     713      \protect\label{fig:OBS_avgrad} 
    714714      Weights associated with each model grid box (blue lines and numbers) 
    715715      for an observation at -170.5\deg{E}, 56.0\deg{N} with a radial footprint with diameter 1\deg. 
     
    756756          ({\phi_{}}_{\mathrm D}  \;  - \; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\ 
    757757  \end{array} 
    758   % \label{eq:cross} 
     758  % \label{eq:OBS_cross} 
    759759\end{align*} 
    760760point in the opposite direction to the unit normal $\widehat{\mathbf k}$ 
     
    791791    \includegraphics[width=\textwidth]{Fig_ASM_obsdist_local} 
    792792    \caption{ 
    793       \protect\label{fig:obslocal} 
     793      \protect\label{fig:OBS_local} 
    794794      Example of the distribution of observations with the geographical distribution of observational data. 
    795795    } 
     
    800800This is the simplest option in which the observations are distributed according to 
    801801the domain of the grid-point parallelization. 
    802 \autoref{fig:obslocal} shows an example of the distribution of the {\em in situ} data on processors with 
     802\autoref{fig:OBS_local} shows an example of the distribution of the {\em in situ} data on processors with 
    803803a different colour for each observation on a given processor for a 4 $\times$ 2 decomposition with ORCA2. 
    804804The grid-point domain decomposition is clearly visible on the plot. 
     
    820820    \includegraphics[width=\textwidth]{Fig_ASM_obsdist_global} 
    821821    \caption{ 
    822       \protect\label{fig:obsglobal} 
     822      \protect\label{fig:OBS_global} 
    823823      Example of the distribution of observations with the round-robin distribution of observational data. 
    824824    } 
     
    830830use message passing in order to retrieve the stencil for interpolation. 
    831831The simplest distribution of the observations is to distribute them using a round-robin scheme. 
    832 \autoref{fig:obsglobal} shows the distribution of the {\em in situ} data on processors for 
     832\autoref{fig:OBS_global} shows the distribution of the {\em in situ} data on processors for 
    833833the round-robin distribution of observations with a different colour for each observation on a given processor for 
    834 a 4 $\times$ 2 decomposition with ORCA2 for the same input data as in \autoref{fig:obslocal}. 
     834a 4 $\times$ 2 decomposition with ORCA2 for the same input data as in \autoref{fig:OBS_local}. 
    835835The observations are now clearly randomly distributed on the globe. 
    836836In order to be able to perform horizontal interpolation in this case, 
     
    11181118\end{minted} 
    11191119 
    1120 \autoref{fig:obsdataplotmain} shows the main window which is launched when dataplot starts. 
     1120\autoref{fig:OBS_dataplotmain} shows the main window which is launched when dataplot starts. 
    11211121This is split into three parts. 
    11221122At the top there is a menu bar which contains a variety of drop down menus. 
     
    11541154    \includegraphics[width=\textwidth]{Fig_OBS_dataplot_main} 
    11551155    \caption{ 
    1156       \protect\label{fig:obsdataplotmain} 
     1156      \protect\label{fig:OBS_dataplotmain} 
    11571157      Main window of dataplot. 
    11581158    } 
     
    11621162 
    11631163If a profile point is clicked with the mouse button a plot of the observation and background values as 
    1164 a function of depth (\autoref{fig:obsdataplotprofile}). 
     1164a function of depth (\autoref{fig:OBS_dataplotprofile}). 
    11651165 
    11661166%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    11701170    \includegraphics[width=\textwidth]{Fig_OBS_dataplot_prof} 
    11711171    \caption{ 
    1172       \protect\label{fig:obsdataplotprofile} 
     1172      \protect\label{fig:OBS_dataplotprofile} 
    11731173      Profile plot from dataplot produced by right clicking on a point in the main window. 
    11741174    } 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex

    r11537 r11543  
    109109Next, the scheme