# Changeset 10368

Ignore:
Timestamp:
2018-12-03T12:45:01+01:00 (23 months ago)
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10365, see #2133

Location:
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE
Files:
74 edited

Unmodified
Removed
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/field_def_nemo-oce.xml

 r9990
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/field_def_nemo-pisces.xml

 r10345
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_pisces_ref

 r10345 !                         !  ln_ligand enabled wfep       =  0.2      ! FeP sinking speed ldocp      =  1.E-5    ! Phyto ligand production per unit doc ldocz      =  1.E-5    ! Zoo ligand production per unit doc ldocp      =  1.E-4    ! Phyto ligand production per unit doc ldocz      =  1.E-4    ! Zoo ligand production per unit doc lthet      =  0.5      ! Proportional loss of ligands due to Fe uptake !                         !  ln_p5z enabled resrat2    =  0.005    ! exsudation rate of mesozooplankton mzrat2     =  0.03     ! mesozooplankton mortality rate xprefc     =  1.       ! mesozoo preference for diatoms xprefp     =  0.3      ! mesozoo preference for nanophyto. xprefz     =  1.       ! mesozoo preference for microzoo. xprefpoc   =  0.3      ! mesozoo preference for poc xpref2d    =  1.       ! mesozoo preference for diatoms xpref2n    =  0.3      ! mesozoo preference for nanophyto. xpref2z    =  1.       ! mesozoo preference for microzoo. xpref2c    =  0.3      ! mesozoo preference for poc xthresh2zoo = 1E-8     ! zoo feeding threshold for mesozooplankton xthresh2dia = 1E-8     ! diatoms feeding threshold for mesozooplankton xkgraz2    =  20.E-6   ! half saturation constant for meso grazing epsher2    =  0.35     ! Efficicency of Mesozoo growth epsher2min =  0.35     ! Minimum efficiency of mesozoo growth sigma2     =  0.6      ! Fraction of mesozoo excretion as DOM unass2     =  0.3      ! non assimilated fraction of P by mesozoo xkgraz2     =  20.E-6   ! half sturation constant for meso grazing epsher2     =  0.5      ! Efficicency of Mesozoo growth epsher2min  =  0.2     ! Minimum efficiency of mesozoo growth ssigma2     =  0.5     ! Fraction excreted as semi-labile DOM srespir2    =  0.2     ! Active respiration resrat     =  0.03     ! exsudation rate of zooplankton mzrat      =  0.004    ! zooplankton mortality rate xpref2c    =  0.1      ! Microzoo preference for POM xpref2p    =  1.       ! Microzoo preference for Nanophyto xpref2d    =  0.5      ! Microzoo preference for Diatoms xprefc     =  0.1      ! Microzoo preference for POM xprefn     =  1.       ! Microzoo preference for Nanophyto xprefd     =  0.6      ! Microzoo preference for Diatoms xthreshdia =  1.E-8    ! Diatoms feeding threshold for microzooplankton xthreshphy =  1.E-8    ! Nanophyto feeding threshold for microzooplankton xkgraz     =  20.E-6   ! half sturation constant for grazing epsher     =  0.3      ! Efficiency of microzoo growth epshermin  =  0.3      ! Minimum efficiency of microzoo growth sigma1     =  0.6      ! Fraction of microzoo excretion as DOM unass      =  0.3      ! non assimilated fraction of phyto by zoo xkgraz     =  20.E-6   ! half sturation constant for grazing epsher     =  0.5      ! Efficiency of microzoo growth epshermin  =  0.2      ! Minimum efficiency of microzoo growth ssigma     =  0.5      ! Fraction excreted as semi-labile DOM srespir    =  0.2      ! Active respiration ln_fechem =  .false.   ! complex iron chemistry ( T/F ) ln_ligvar =  .false.   ! variable ligand concentration ln_fecolloid = .false. ! variable colloidal fraction xlam1        =  0.005  ! scavenging rate of Iron xlamdust     =  150.0  ! Scavenging rate of dust ligand       =  0.6E-9 ! Ligands concentration kfep         =  0.     ! Nanoparticle formation rate constant kfep         =  0.01   ! Nanoparticle formation rate constant / !----------------------------------------------------------------------- fep_rats    =  1.       ! Fep/Fer ratio from sed sources fep_rath    =  1.       ! Fep/Fer ratio from sed hydro sources rdustfep    =  0.0      ! Fraction of dust that is FeP lgw_rath    =  0.5      ! Weak ligand ratio from sed hydro sources / !----------------------------------------------------------------------- rfep        =  0.001    ! Dissolution rate of FeP rlgw        =  1.       ! Lifetime (years) of weak ligands rlgw        =  100.     ! Lifetime (years) of weak ligands rlig        =  1.E-4    ! Remin ligand production per unit C prlgw       =  1.E-4    ! Photolysis of weak ligand rlgs        =  1000.    ! Lifetime (years) of strong ligands rlgs        =  1.       ! Lifetime (years) of strong ligands / !-----------------------------------------------------------------------
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_ref

 r10345 &namzdf        !   vertical physics manager                             (default: NO selection) !----------------------------------------------------------------------- !                       ! adaptive-implicit vertical advection ln_zad_Aimp = .false.      !  Courant number dependent scheme (Shchepetkin 2015) ! !                       ! type of vertical closure (required) ln_zdfcst   = .false.      !  constant mixing !                       !  buffer blocking send or immediate non-blocking sends, resp. nn_buffer   =   0       !  size in bytes of exported buffer ('B' case), 0 no exportation ln_nnogather =  .false.  !  activate code to avoid mpi_allgather use at the northfold ln_nnogather =  .true.  !  activate code to avoid mpi_allgather use at the northfold jpni        =   0       !  jpni   number of processors following i (set automatically if < 1) jpnj        =   0       !  jpnj   number of processors following j (set automatically if < 1)
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/abstract_foreword.tex

 r9393 {\small The ocean engine of NEMO (Nucleus for European Modelling of the Ocean) is a primitive equation model adapted to regional and global ocean circulation problems. It is intended to be a flexible tool for studying the ocean and its interactions with the others components of the earth climate system over a wide range of space and time scales. Prognostic variables are the three-dimensional velocity field, a non-linear sea surface height, the \textit{Conservative} Temperature and the \textit{Absolute} Salinity. In the horizontal direction, the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step $z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a three-dimensional Arakawa C-type grid. Various physical choices are available to describe ocean physics, including TKE, and GLS vertical physics. Within NEMO, the ocean is interfaced with a sea-ice model (LIM or CICE), passive tracer and biogeochemical models (TOP) and, via the OASIS coupler, with several atmospheric general circulation models. It also support two-way grid embedding via the AGRIF software. The ocean engine of NEMO (Nucleus for European Modelling of the Ocean) is a primitive equation model adapted to regional and global ocean circulation problems. It is intended to be a flexible tool for studying the ocean and its interactions with the others components of the earth climate system over a wide range of space and time scales. Prognostic variables are the three-dimensional velocity field, a non-linear sea surface height, the \textit{Conservative} Temperature and the \textit{Absolute} Salinity. In the horizontal direction, the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step $z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a three-dimensional Arakawa C-type grid. Various physical choices are available to describe ocean physics, including TKE, and GLS vertical physics. Within NEMO, the ocean is interfaced with a sea-ice model (LIM or CICE), passive tracer and biogeochemical models (TOP) and, via the OASIS coupler, with several atmospheric general circulation models. It also support two-way grid embedding via the AGRIF software. } \chapter*{Disclaimer} Like all components of NEMO, the ocean component is developed under the \href{http://www.cecill.info/}{CECILL license}, which is a French adaptation of the GNU GPL (General Public License). Anyone may use it freely for research purposes, and is encouraged to communicate back to the NEMO team its own developments and improvements. The model and the present document have been made available as a service to the community. We cannot certify that the code and its manual are free of errors. Bugs are inevitable and some have undoubtedly survived the testing phase. Users are encouraged to bring them to our attention. The author assumes no responsibility for problems, errors, or incorrect usage of NEMO. Like all components of NEMO, the ocean component is developed under the \href{http://www.cecill.info/}{CECILL license}, which is a French adaptation of the GNU GPL (General Public License). Anyone may use it freely for research purposes, and is encouraged to communicate back to the NEMO team its own developments and improvements. The model and the present document have been made available as a service to the community. We cannot certify that the code and its manual are free of errors. Bugs are inevitable and some have undoubtedly survived the testing phase. Users are encouraged to bring them to our attention. The author assumes no responsibility for problems, errors, or incorrect usage of NEMO. \vspace{1cm}
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_A.tex

 r9414 In order to establish the set of Primitive Equation in curvilinear $s$-coordinates ($i.e.$ an orthogonal curvilinear coordinate in the horizontal and an Arbitrary Lagrangian Eulerian (ALE) coordinate in the vertical), we start from the set of equations established in \autoref{subsec:PE_zco_Eq} for the special case $k = z$ and thus $e_3 = 1$, and we introduce an arbitrary vertical coordinate $a = a(i,j,z,t)$. Let us define a new vertical scale factor by $e_3 = \partial z / \partial s$ (which now depends on $(i,j,z,t)$) and the horizontal slope of $s-$surfaces by : ($i.e.$ an orthogonal curvilinear coordinate in the horizontal and an Arbitrary Lagrangian Eulerian (ALE) coordinate in the vertical), we start from the set of equations established in \autoref{subsec:PE_zco_Eq} for the special case $k = z$ and thus $e_3 = 1$, and we introduce an arbitrary vertical coordinate $a = a(i,j,z,t)$. Let us define a new vertical scale factor by $e_3 = \partial z / \partial s$ (which now depends on $(i,j,z,t)$) and the horizontal slope of $s-$surfaces by: \label{apdx:A_s_slope} \sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s The chain rule to establish the model equations in the curvilinear $s-$coordinate system is: The chain rule to establish the model equations in the curvilinear $s-$coordinate system is: \label{apdx:A_s_chain_rule} \begin{aligned} In particular applying the time derivative chain rule to $z$ provides the expression for $w_s$,  the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate: In particular applying the time derivative chain rule to $z$ provides the expression for $w_s$, the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate: \label{apdx:A_w_in_s} w_s   =  \left.   \frac{\partial z }{\partial t}   \right|_s \label{sec:A_continuity} Using (\autoref{apdx:A_s_chain_rule}) and the fact that the horizontal scale factors $e_1$ and $e_2$ do not depend on the vertical coordinate, the divergence of the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows in order to obtain its expression in the curvilinear $s-$coordinate system: Using (\autoref{apdx:A_s_chain_rule}) and the fact that the horizontal scale factors $e_1$ and $e_2$ do not depend on the vertical coordinate, the divergence of the velocity relative to the ($i$,$j$,$z$) coordinate system is transformed as follows in order to obtain its expression in the curvilinear $s-$coordinate system: \begin{subequations} \end{subequations} Here, $w$ is the vertical velocity relative to the $z-$coordinate system. Introducing the dia-surface velocity component, $\omega$, defined as the volume flux across the moving $s$-surfaces per unit horizontal area: Here, $w$ is the vertical velocity relative to the $z-$coordinate system. Introducing the dia-surface velocity component, $\omega$, defined as the volume flux across the moving $s$-surfaces per unit horizontal area: \label{apdx:A_w_s} \omega  = w - w_s - \sigma _1 \,u - \sigma _2 \,v    \\ with $w_s$ given by \autoref{apdx:A_w_in_s}, we obtain the expression for the divergence of the velocity in the curvilinear $s-$coordinate system: with $w_s$ given by \autoref{apdx:A_w_in_s}, we obtain the expression for the divergence of the velocity in the curvilinear $s-$coordinate system: \begin{subequations} \begin{align*} {\begin{array}{*{20}l} \end{subequations} As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is: As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is: \label{apdx:A_sco_Continuity} \frac{1}{e_3 } \frac{\partial e_3}{\partial t} +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0 A additional term has appeared that take into account the contribution of the time variation of the vertical coordinate to the volume budget. A additional term has appeared that take into account the contribution of the time variation of the vertical coordinate to the volume budget. \label{sec:A_momentum} Here we only consider the first component of the momentum equation, Here we only consider the first component of the momentum equation, the generalization to the second one being straightforward. $\bullet$ \textbf{Total derivative in vector invariant form} Let us consider \autoref{eq:PE_dyn_vect}, the first component of the momentum equation in the vector invariant form. Its total $z-$coordinate time derivative, $\left. \frac{D u}{D t} \right|_z$ can be transformed as follows in order to obtain Let us consider \autoref{eq:PE_dyn_vect}, the first component of the momentum equation in the vector invariant form. Its total $z-$coordinate time derivative, $\left. \frac{D u}{D t} \right|_z$ can be transformed as follows in order to obtain its expression in the curvilinear $s-$coordinate system: \end{subequations} % Applying the time derivative chain rule (first equation of (\autoref{apdx:A_s_chain_rule})) to $u$ and using (\autoref{apdx:A_w_in_s}) provides the expression of the last term of the right hand side, Applying the time derivative chain rule (first equation of (\autoref{apdx:A_s_chain_rule})) to $u$ and using (\autoref{apdx:A_w_in_s}) provides the expression of the last term of the right hand side, \begin{equation*} {\begin{array}{*{20}l} w_s  \;\frac{\partial u}{\partial s} \end{array} } \end{equation*} leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative, leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative, $i.e.$ the total $s-$coordinate time derivative : \begin{align} \label{apdx:A_sco_Dt_vect} + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} \end{align} Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in $z-$ and $s-$coordinates. This is not the case for the flux form as shown in next paragraph. Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in $z-$ and $s-$coordinates. This is not the case for the flux form as shown in next paragraph. $\$\newline    % force a new ligne $\bullet$ \textbf{Total derivative in flux form} Let us start from the total time derivative in the curvilinear $s-$coordinate system we have just establish. Following the procedure used to establish (\autoref{eq:PE_flux_form}), it can be transformed into : Let us start from the total time derivative in the curvilinear $s-$coordinate system we have just establish. Following the procedure used to establish (\autoref{eq:PE_flux_form}), it can be transformed into : %\begin{subequations} \begin{align*} {\begin{array}{*{20}l} \end{subequations} which leads to the $s-$coordinate flux formulation of the total $s-$coordinate time derivative, $i.e.$ the total $s-$coordinate time derivative in flux form : $i.e.$ the total $s-$coordinate time derivative in flux form: \begin{flalign}\label{apdx:A_sco_Dt_flux} \left. \frac{D u}{D t} \right|_s   = \frac{1}{e_3}  \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s \end{flalign} which is the total time derivative expressed in the curvilinear $s-$coordinate system. It has the same form as in the $z-$coordinate but for the vertical scale factor that has appeared inside the time derivative which comes from the modification of (\autoref{apdx:A_sco_Continuity}), the continuity equation. It has the same form as in the $z-$coordinate but for the vertical scale factor that has appeared inside the time derivative which comes from the modification of (\autoref{apdx:A_sco_Continuity}), the continuity equation. $\$\newline    % force a new ligne \end{split} \end{equation*} Applying similar manipulation to the second component and replacing $\sigma _1$ and $\sigma _2$ by their expression \autoref{apdx:A_s_slope}, it comes: Applying similar manipulation to the second component and replacing $\sigma _1$ and $\sigma _2$ by their expression \autoref{apdx:A_s_slope}, it comes: \label{apdx:A_grad_p_1} \begin{split} An additional term appears in (\autoref{apdx:A_grad_p_1}) which accounts for the tilt of $s-$surfaces with respect to geopotential $z-$surfaces. As in $z$-coordinate, the horizontal pressure gradient can be split in two parts following \citet{Marsaleix_al_OM08}. Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$, and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$. An additional term appears in (\autoref{apdx:A_grad_p_1}) which accounts for the tilt of $s-$surfaces with respect to geopotential $z-$surfaces. As in $z$-coordinate, the horizontal pressure gradient can be split in two parts following \citet{Marsaleix_al_OM08}. Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$, and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$. The pressure is then given by: \begin{equation*} \end{equation*} Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and using the definition of the density anomaly it comes the expression in two parts: Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and using the definition of the density anomaly it comes the expression in two parts: \label{apdx:A_grad_p_2} \begin{split} \end{split} This formulation of the pressure gradient is characterised by the appearance of a term depending on the the sea surface height only (last term on the right hand side of expression \autoref{apdx:A_grad_p_2}). This term will be loosely termed \textit{surface pressure gradient} whereas the first term will be termed the \textit{hydrostatic pressure gradient} by analogy to the $z$-coordinate formulation. In fact, the the true surface pressure gradient is $1/\rho_o \nabla (\rho \eta)$, and $\eta$ is implicitly included in the computation of $p_h'$ through the upper bound of the vertical integration. This formulation of the pressure gradient is characterised by the appearance of a term depending on the sea surface height only (last term on the right hand side of expression \autoref{apdx:A_grad_p_2}). This term will be loosely termed \textit{surface pressure gradient} whereas the first term will be termed the \textit{hydrostatic pressure gradient} by analogy to the $z$-coordinate formulation. In fact, the true surface pressure gradient is $1/\rho_o \nabla (\rho \eta)$, and $\eta$ is implicitly included in the computation of $p_h'$ through the upper bound of the vertical integration. $\$\newline    % force a new ligne $\bullet$ \textbf{The other terms of the momentum equation} The coriolis and forcing terms as well as the the vertical physics remain unchanged as they involve neither time nor space derivatives. The form of the lateral physics is discussed in \autoref{apdx:B}. The coriolis and forcing terms as well as the the vertical physics remain unchanged as they involve neither time nor space derivatives. The form of the lateral physics is discussed in \autoref{apdx:B}. $\bullet$ \textbf{Full momentum equation} To sum up, in a curvilinear $s$-coordinate system, the vector invariant momentum equation solved by the model has the same mathematical expression as the one in a curvilinear $z-$coordinate, except for the pressure gradient term : To sum up, in a curvilinear $s$-coordinate system, the vector invariant momentum equation solved by the model has the same mathematical expression as the one in a curvilinear $z-$coordinate, except for the pressure gradient term: \begin{subequations} \label{apdx:A_dyn_vect} \begin{multline} \label{apdx:A_PE_dyn_vect_u} \end{multline} \end{subequations} whereas the flux form momentum equation differ from it by the formulation of both the time derivative and the pressure gradient term  : whereas the flux form momentum equation differs from it by the formulation of both the time derivative and the pressure gradient term: \begin{subequations} \label{apdx:A_dyn_flux} \begin{multline} \label{apdx:A_PE_dyn_flux_u} It is important to realize that the change in coordinate system has only concerned the position on the vertical. It has not affected (\textbf{i},\textbf{j},\textbf{k}), the orthogonal curvilinear set of unit vectors. ($u$,$v$) are always horizontal velocities so that their evolution is driven by \emph{horizontal} forces, in particular the pressure gradient. By contrast, $\omega$ is not $w$, the third component of the velocity, but the dia-surface velocity component, $i.e.$ the volume flux across the moving $s$-surfaces per unit horizontal area. It is important to realize that the change in coordinate system has only concerned the position on the vertical. It has not affected (\textbf{i},\textbf{j},\textbf{k}), the orthogonal curvilinear set of unit vectors. ($u$,$v$) are always horizontal velocities so that their evolution is driven by \emph{horizontal} forces, in particular the pressure gradient. By contrast, $\omega$ is not $w$, the third component of the velocity, but the dia-surface velocity component, $i.e.$ the volume flux across the moving $s$-surfaces per unit horizontal area. \label{sec:A_tracer} The tracer equation is obtained using the same calculation as for the continuity equation and then regrouping the time derivative terms in the left hand side : The tracer equation is obtained using the same calculation as for the continuity equation and then regrouping the time derivative terms in the left hand side : \begin{multline} \label{apdx:A_tracer} The expression for the advection term is a straight consequence of (A.4), the expression of the 3D divergence in the $s-$coordinates established above. The expression for the advection term is a straight consequence of (A.4), the expression of the 3D divergence in the $s-$coordinates established above. \end{document}
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_B.tex

 r9407 \subsubsection*{In z-coordinates} In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator is given by: In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator is given by: \begin{eqnarray} \label{apdx:B1} &D^T = \frac{1}{e_1 \, e_2}      \left[ \subsubsection*{In generalized vertical coordinates} In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by \autoref{apdx:A_s_slope} and the vertical/horizontal ratio of diffusion coefficient by $\epsilon = A^{vT} / A^{lT}$. The diffusion operator is given by: In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and $\sigma_2$ by \autoref{apdx:A_s_slope} and the vertical/horizontal ratio of diffusion coefficient by $\epsilon = A^{vT} / A^{lT}$. The diffusion operator is given by: \label{apdx:B2} \end{subequations} Equation \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption. Indeed, for the special case $k=z$ and thus $e_3 =1$, we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and use \autoref{apdx:A_s_slope} and \autoref{apdx:A_s_chain_rule}. Since no cross horizontal derivative $\partial _i \partial _j$ appears in \autoref{apdx:B1}, the ($i$,$z$) and ($j$,$z$) planes are independent. The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$) transformation without any loss of generality: Equation \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption. Indeed, for the special case $k=z$ and thus $e_3 =1$, we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and use \autoref{apdx:A_s_slope} and \autoref{apdx:A_s_chain_rule}. Since no cross horizontal derivative $\partial _i \partial _j$ appears in \autoref{apdx:B1}, the ($i$,$z$) and ($j$,$z$) planes are independent. The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$) transformation without any loss of generality: \begin{subequations} \subsubsection*{In z-coordinates} The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in the ($i$,$j$,$k$) curvilinear coordinate system in which the equations of the ocean circulation model are formulated, takes the following form \citep{Redi_JPO82}: The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in the ($i$,$j$,$k$) curvilinear coordinate system in which the equations of the ocean circulation model are formulated, takes the following form \citep{Redi_JPO82}: \label{apdx:B3} \end{array} }} \right] where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials: where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials: \begin{equation*} a_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \end{equation*} In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean, so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}: In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean, so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}: \begin{subequations} \label{apdx:B4} \label{apdx:B4a} Physically, the full tensor \autoref{apdx:B3} represents strong isoneutral diffusion on a plane parallel to the isoneutral surface and weak dianeutral diffusion perpendicular to this plane. However, the approximate weak-slope' tensor \autoref{apdx:B4a} represents strong diffusion along the isoneutral surface, with weak \emph{vertical}  diffusion -- the principal axes of the tensor are no longer orthogonal. This simplification also decouples the ($i$,$z$) and ($j$,$z$) planes of the tensor. The weak-slope operator therefore takes the same form, \autoref{apdx:B4}, as \autoref{apdx:B2}, the diffusion operator for geopotential diffusion written in non-orthogonal $i,j,s$-coordinates. Written out explicitly, Physically, the full tensor \autoref{apdx:B3} represents strong isoneutral diffusion on a plane parallel to the isoneutral surface and weak dianeutral diffusion perpendicular to this plane. However, the approximate weak-slope' tensor \autoref{apdx:B4a} represents strong diffusion along the isoneutral surface, with weak \emph{vertical} diffusion -- the principal axes of the tensor are no longer orthogonal. This simplification also decouples the ($i$,$z$) and ($j$,$z$) planes of the tensor. The weak-slope operator therefore takes the same form, \autoref{apdx:B4}, as \autoref{apdx:B2}, the diffusion operator for geopotential diffusion written in non-orthogonal $i,j,s$-coordinates. Written out explicitly, \begin{multline} \label{apdx:B_ldfiso} The isopycnal diffusion operator \autoref{apdx:B4}, \autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square. The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence of fluxes. Let us demonstrate the second one: \autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square. The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence of fluxes. Let us demonstrate the second one: \begin{equation*} \iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv \end{subequations} \addtocounter{equation}{-1} the property becomes obvious. the property becomes obvious. \subsubsection*{In generalized vertical coordinates} Because the weak-slope operator \autoref{apdx:B4}, \autoref{apdx:B_ldfiso} is decoupled in the ($i$,$z$) and ($j$,$z$) planes, it may be transformed into generalized $s$-coordinates in the same way as \autoref{sec:B_1} was transformed into \autoref{sec:B_2}. The resulting operator then takes the simple form Because the weak-slope operator \autoref{apdx:B4}, \autoref{apdx:B_ldfiso} is decoupled in the ($i$,$z$) and ($j$,$z$) planes, it may be transformed into generalized $s$-coordinates in the same way as \autoref{sec:B_1} was transformed into \autoref{sec:B_2}. The resulting operator then takes the simple form \label{apdx:B_ldfiso_s} where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to $s$-coordinate surfaces: where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to $s$-coordinate surfaces: \begin{equation*} r_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial s}} \right)^{-1} \end{equation*} To prove  \autoref{apdx:B5}  by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious. An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as To prove \autoref{apdx:B5} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious. An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as \label{apdx:B5} \end{array} }} \right). Then direct transformation from $i,j,\rho$-coordinates to $i,j,s$-coordinates gives \autoref{apdx:B_ldfiso_s} immediately. Note that the weak-slope approximation is only made in transforming from the (rotated,orthogonal) isoneutral axes to the non-orthogonal $i,j,\rho$-coordinates. The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the  $s$-surfaces, in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates, Then direct transformation from $i,j,\rho$-coordinates to $i,j,s$-coordinates gives \autoref{apdx:B_ldfiso_s} immediately. Note that the weak-slope approximation is only made in transforming from the (rotated,orthogonal) isoneutral axes to the non-orthogonal $i,j,\rho$-coordinates. The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the $s$-surfaces, in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates, \autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces. \label{sec:B_3} The second order momentum diffusion operator (Laplacian) in the $z$-coordinate is found by applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector,  to the horizontal velocity vector : The second order momentum diffusion operator (Laplacian) in the $z$-coordinate is found by applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector, to the horizontal velocity vector: \begin{align*} \Delta {\textbf{U}}_h \end{array} }} \right) \end{align*} Using \autoref{eq:PE_div}, the definition of the horizontal divergence, the third componant of the second vector is obviously zero and thus : Using \autoref{eq:PE_div}, the definition of the horizontal divergence, the third componant of the second vector is obviously zero and thus : \begin{equation*} \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) \end{equation*} Note that this operator ensures a full separation between the vorticity and horizontal divergence fields (see \autoref{apdx:C}). It is only equal to a Laplacian applied to each component in Cartesian coordinates, not on the sphere. The horizontal/vertical second order (Laplacian type) operator used to diffuse horizontal momentum in the $z$-coordinate therefore takes the following form : Note that this operator ensures a full separation between the vorticity and horizontal divergence fields (see \autoref{apdx:C}). It is only equal to a Laplacian applied to each component in Cartesian coordinates, not on the sphere. The horizontal/vertical second order (Laplacian type) operator used to diffuse horizontal momentum in the $z$-coordinate therefore takes the following form: \label{apdx:B_Lap_U} {\textbf{D}}^{\textbf{U}} = \end{align*} Note Bene: introducing a rotation in \autoref{apdx:B_Lap_U} does not lead to a useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate. Similarly, we did not found an expression of practical use for the geopotential horizontal/vertical Laplacian operator in the $s$-coordinate. Generally, \autoref{apdx:B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is a Laplacian diffusion is applied on momentum along the coordinate directions. Note Bene: introducing a rotation in \autoref{apdx:B_Lap_U} does not lead to a useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate. Similarly, we did not found an expression of practical use for the geopotential horizontal/vertical Laplacian operator in the $s$-coordinate. Generally, \autoref{apdx:B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is a Laplacian diffusion is applied on momentum along the coordinate directions. \end{document}
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_C.tex

 r9414 \label{sec:C.0} Notation used in this appendix in the demonstations : Notation used in this appendix in the demonstations: fluxes at the faces of a $T$-box: $dv=e_1\,e_2\,e_3 \,di\,dj\,dk$  is the volume element, with only $e_3$ that depends on time. $D$ and $S$ are the ocean domain volume and surface, respectively. No wetting/drying is allow ($i.e.$ $\frac{\partial S}{\partial t} = 0$) Let $k_s$ and $k_b$ be the ocean surface and bottom, resp. $D$ and $S$ are the ocean domain volume and surface, respectively. No wetting/drying is allow ($i.e.$ $\frac{\partial S}{\partial t} = 0$). Let $k_s$ and $k_b$ be the ocean surface and bottom, resp. ($i.e.$ $s(k_s) = \eta$ and $s(k_b)=-H$, where $H$ is the bottom depth). \begin{flalign*} =  \int_D { \frac{1}{e_3} \partial_t \left( e_3 \, Q \right) dv } =0 \end{equation*} equation of evolution of $Q$ written as the time evolution of the vertical content of $Q$ like for tracers, or momentum in flux form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when : equation of evolution of $Q$ written as the time evolution of the vertical content of $Q$ like for tracers, or momentum in flux form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when: \begin{flalign*} \partial_t \left(   \int_D{ \frac{1}{2} \,Q^2\;dv }   \right) -   \frac{1}{2} \int_D {  \frac{Q^2}{e_3} \partial_t (e_3) \;dv } \end{flalign} equation of evolution of $Q$ written as the time evolution of $Q$ like for momentum in vector invariant form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when : equation of evolution of $Q$ written as the time evolution of $Q$ like for momentum in vector invariant form, the quadratic quantity $\frac{1}{2}Q^2$ is conserved when: \begin{flalign*} \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right) +  \int_D { \frac{1}{2} Q^2 \, \partial_t e_3  \;e_1e_2\;di\,dj\,dk } \\ \end{flalign*} that is in a more compact form : that is in a more compact form: \begin{flalign} \label{eq:Q2_vect} \partial_t \left( \int_D {\frac{1}{2} Q^2\;dv} \right) The discretization of pimitive equation in $s$-coordinate ($i.e.$ time and space varying vertical coordinate) must be chosen so that the discrete equation of the model satisfy integral constrains on energy and enstrophy. The discretization of pimitive equation in $s$-coordinate ($i.e.$ time and space varying vertical coordinate) must be chosen so that the discrete equation of the model satisfy integral constrains on energy and enstrophy. Let us first establish those constraint in the continuous world. The total energy ($i.e.$ kinetic plus potential energies) is conserved : The total energy ($i.e.$ kinetic plus potential energies) is conserved: \begin{flalign} \label{eq:Tot_Energy} \partial_t \left( \int_D \left( \frac{1}{2} {\textbf{U}_h}^2 +  \rho \, g \, z \right) \;dv \right)  = & 0 \end{flalign} under the following assumptions: no dissipation, no forcing (wind, buoyancy flux, atmospheric pressure variations), mass conservation, and closed domain. This equation can be transformed to obtain several sub-equalities. The transformation for the advection term depends on whether the vector invariant form or the flux form is used for the momentum equation. Using \autoref{eq:Q2_vect} and introducing \autoref{apdx:A_dyn_vect} in \autoref{eq:Tot_Energy} for the former form and Using \autoref{eq:Q2_flux} and introducing \autoref{apdx:A_dyn_flux} in \autoref{eq:Tot_Energy} for the latter form  leads to: under the following assumptions: no dissipation, no forcing (wind, buoyancy flux, atmospheric pressure variations), mass conservation, and closed domain. This equation can be transformed to obtain several sub-equalities. The transformation for the advection term depends on whether the vector invariant form or the flux form is used for the momentum equation. Using \autoref{eq:Q2_vect} and introducing \autoref{apdx:A_dyn_vect} in \autoref{eq:Tot_Energy} for the former form and using \autoref{eq:Q2_flux} and introducing \autoref{apdx:A_dyn_flux} in \autoref{eq:Tot_Energy} for the latter form leads to: \begin{subequations} \label{eq:E_tot} Substituting the discrete expression of the time derivative of the velocity either in vector invariant, leads to the discrete equivalent of the four equations \autoref{eq:E_tot_flux}. leads to the discrete equivalent of the four equations \autoref{eq:E_tot_flux}. % ------------------------------------------------------------------------------------------------------------- \label{subsec:C_vor} Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$), or the planetary ($q=f/e_{3f}$), or the total potential vorticity ($q=(\zeta +f) /e_{3f}$). Two discretisation of the vorticity term (ENE and EEN) allows the conservation of the kinetic energy. Let $q$, located at $f$-points, be either the relative ($q=\zeta / e_{3f}$), or the planetary ($q=f/e_{3f}$), or the total potential vorticity ($q=(\zeta +f) /e_{3f}$). Two discretisation of the vorticity term (ENE and EEN) allows the conservation of the kinetic energy. % ------------------------------------------------------------------------------------------------------------- %       Vorticity Term with ENE scheme \label{subsec:C_vorENE} For the ENE scheme, the two components of the vorticity term are given by : For the ENE scheme, the two components of the vorticity term are given by: \begin{equation*} - e_3 \, q \;{\textbf{k}}\times {\textbf {U}}_h    \equiv \end{equation*} This formulation does not conserve the enstrophy but it does conserve the total kinetic energy. Indeed, the kinetic energy tendency associated to the vorticity term and averaged over the ocean domain can be transformed as follows: This formulation does not conserve the enstrophy but it does conserve the total kinetic energy. Indeed, the kinetic energy tendency associated to the vorticity term and averaged over the ocean domain can be transformed as follows: \begin{flalign*} &\int\limits_D -  \left(  e_3 \, q \;\textbf{k} \times \textbf{U}_h  \right) \cdot \textbf{U}_h  \;  dv &&  \\ \end{aligned}   } \right. where the indices $i_p$ and $j_p$ take the following value: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$, where the indices $i_p$ and $j_p$ take the following value: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$, and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: \tag{\ref{eq:Q_triads}} This formulation does conserve the total kinetic energy. Indeed, This formulation does conserve the total kinetic energy. Indeed, \begin{flalign*} &\int\limits_D - \textbf{U}_h \cdot   \left(  \zeta \;\textbf{k} \times \textbf{U}_h  \right)  \;  dv &&  \\ \label{subsec:C_zad} The 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~: The 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~: \begin{equation*} \int_D \textbf{U}_h \cdot \frac{1}{e_3 } \omega \partial_k \textbf{U}_h \;dv +   \frac{1}{2} \int_D {  \frac{{\textbf{U}_h}^2}{e_3} \partial_t ( e_3) \;dv }  \\ \end{equation*} Indeed, using successively \autoref{eq:DOM_di_adj} ($i.e.$ the skew symmetry property of the $\delta$ operator) and the continuity equation, then \autoref{eq:DOM_di_adj} again, then the commutativity of operators $\overline {\,\cdot \,}$ and $\delta$, and finally \autoref{eq:DOM_mi_adj} ($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator) Indeed, using successively \autoref{eq:DOM_di_adj} ($i.e.$ the skew symmetry property of the $\delta$ operator) and the continuity equation, then \autoref{eq:DOM_di_adj} again, then the commutativity of operators $\overline {\,\cdot \,}$ and $\delta$, and finally \autoref{eq:DOM_mi_adj} ($i.e.$ the symmetry property of the $\overline {\,\cdot \,}$ operator) applied in the horizontal and vertical directions, it becomes: \begin{flalign*} \end{flalign*} There is two main points here. First, the satisfaction of this property links the choice of the discrete formulation of the vertical advection and of the horizontal gradient of KE. Choosing one imposes the other. For example KE can also be discretized as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$. This leads to the following expression for the vertical advection: There is two main points here. First, the satisfaction of this property links the choice of the discrete formulation of the vertical advection and of the horizontal gradient of KE. Choosing one imposes the other. For example KE can also be discretized as $1/2\,({\overline u^{\,i}}^2 + {\overline v^{\,j}}^2)$. This leads to the following expression for the vertical advection: \begin{equation*} \frac{1} {e_3 }\; \omega\; \partial_k \textbf{U}_h \end{array}} } \right) \end{equation*} a formulation that requires an additional horizontal mean in contrast with the one used in NEMO. Nine velocity points have to be used instead of 3. a formulation that requires an additional horizontal mean in contrast with the one used in NEMO. Nine velocity points have to be used instead of 3. This is the reason why it has not been chosen. Second, as soon as the chosen $s$-coordinate depends on time, an extra constraint arises on the time derivative of the volume at $u$- and $v$-points: Second, as soon as the chosen $s$-coordinate depends on time, an extra constraint arises on the time derivative of the volume at $u$- and $v$-points: \begin{flalign*} e_{1u}\,e_{2u}\,\partial_t (e_{3u}) =\overline{ e_{1t}\,e_{2t}\;\partial_t (e_{3t}) }^{\,i+1/2}    \\ \gmcomment{ A pressure gradient has no contribution to the evolution of the vorticity as the curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally on a C-grid with 2nd order finite differences (property \autoref{eq:DOM_curl_grad}). A pressure gradient has no contribution to the evolution of the vorticity as the curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally on a C-grid with 2nd order finite differences (property \autoref{eq:DOM_curl_grad}). } When the equation of state is linear ($i.e.$ when an advection-diffusion equation for density can be derived from those of temperature and salinity) the change of KE due to the work of pressure forces is balanced by the change of potential energy due to buoyancy forces: When the equation of state is linear ($i.e.$ when an advection-diffusion equation for density can be derived from those of temperature and salinity) the change of KE due to the work of pressure forces is balanced by the change of potential energy due to buoyancy forces: \begin{equation*} - \int_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv \end{equation*} This property can be satisfied in a discrete sense for both $z$- and $s$-coordinates. Indeed, defining the depth of a $T$-point, $z_t$, as the sum of the vertical scale factors at $w$-points starting from the surface, the work of pressure forces can be written as: This property can be satisfied in a discrete sense for both $z$- and $s$-coordinates. Indeed, defining the depth of a $T$-point, $z_t$, as the sum of the vertical scale factors at $w$-points starting from the surface, the work of pressure forces can be written as: \begin{flalign*} &- \int_D  \left. \nabla p \right|_z \cdot \textbf{U}_h \;dv \end{flalign*} The first term is exactly the first term of the right-hand-side of \autoref{eq:KE+PE_vect_discrete}. It remains to demonstrate that the last term, which is obviously a discrete analogue of $\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to the last term of \autoref{eq:KE+PE_vect_discrete}. It remains to demonstrate that the last term, which is obviously a discrete analogue of $\int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to the last term of \autoref{eq:KE+PE_vect_discrete}. In other words, the following property must be satisfied: \begin{flalign*} \end{flalign*} Let introduce $p_w$ the pressure at $w$-point such that $\delta_k [p_w] = - \rho \,g\,e_{3t}$. Let introduce $p_w$ the pressure at $w$-point such that $\delta_k [p_w] = - \rho \,g\,e_{3t}$. The right-hand-side of the above equation can be transformed as follows: Note that this property strongly constrains the discrete expression of both the depth of $T-$points and of the term added to the pressure gradient in the $s$-coordinate. Nevertheless, it is almost never satisfied since a linear equation of state is rarely used. Note that this property strongly constrains the discrete expression of both the depth of $T-$points and of the term added to the pressure gradient in the $s$-coordinate. Nevertheless, it is almost never satisfied since a linear equation of state is rarely used. \end{flalign*} Substituting the discrete expression of the time derivative of the velocity either in vector invariant or in flux form, leads to the discrete equivalent of the Substituting the discrete expression of the time derivative of the velocity either in vector invariant or in flux form, leads to the discrete equivalent of the ???? \label{subsec:C.3.3} In flux from the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the metric'' term. This altered Coriolis parameter is discretised at an f-point. It is given by: In flux from the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the metric'' term. This altered Coriolis parameter is discretised at an f-point. It is given by: \begin{equation*} f+\frac{1} {e_1 e_2 } \left( v \frac{\partial e_2 } {\partial i} - u \frac{\partial e_1 } {\partial j}\right)\; \end{equation*} Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form. It therefore conserves the total KE. The derivation is the same as for the vorticity term in the vector invariant form (\autoref{subsec:C_vor}). Either the ENE or EEN scheme is then applied to obtain the vorticity term in flux form. It therefore conserves the total KE. The derivation is the same as for the vorticity term in the vector invariant form (\autoref{subsec:C_vor}). % ------------------------------------------------------------------------------------------------------------- \label{subsec:C.3.4} The flux form operator of the momentum advection is evaluated using a centered second order finite difference scheme. Because of the flux form, the discrete operator does not contribute to the global budget of linear momentum. Because of the centered second order scheme, it conserves the horizontal kinetic energy, that is : The flux form operator of the momentum advection is evaluated using a centered second order finite difference scheme. Because of the flux form, the discrete operator does not contribute to the global budget of linear momentum. Because of the centered second order scheme, it conserves the horizontal kinetic energy, that is: \label{eq:C_ADV_KE_flux} Let us first consider the first term of the scalar product ($i.e.$ just the the terms associated with the i-component of the advection) : Let us first consider the first term of the scalar product ($i.e.$ just the the terms associated with the i-component of the advection): \begin{flalign*} &  - \int_D u \cdot \nabla \cdot \left(   \textbf{U}\,u   \right) \; dv   \\ \biggl\{     \left(   \frac{1}{e_{3t}} \frac{\partial e_{3t}}{\partial t}   \right) \; b_t     \biggr\}    &&& \\ \end{flalign*} Applying similar manipulation applied to the second term of the scalar product leads to : Applying similar manipulation applied to the second term of the scalar product leads to: \begin{equation*} -  \int_D \textbf{U}_h \cdot     \left(                 {{\begin{array} {*{20}c} \biggl\{     \left(   \frac{1}{e_{3t}} \frac{\partial e_{3t}}{\partial t}   \right) \; b_t     \biggr\} \end{equation*} which is the discrete form of $\frac{1}{2} \int_D u \cdot \nabla \cdot \left( \textbf{U}\,u \right) \; dv$. which is the discrete form of $\frac{1}{2} \int_D u \cdot \nabla \cdot \left( \textbf{U}\,u \right) \; dv$. \autoref{eq:C_ADV_KE_flux} is thus satisfied. When the UBS scheme is used to evaluate the flux form momentum advection, the discrete operator does not contribute to the global budget of linear momentum (flux form). The horizontal kinetic energy is not conserved, but forced to decay ($i.e.$ the scheme is diffusive). When the UBS scheme is used to evaluate the flux form momentum advection, the discrete operator does not contribute to the global budget of linear momentum (flux form). The horizontal kinetic energy is not conserved, but forced to decay ($i.e.$ the scheme is diffusive). The scheme does not allow but the conservation of the total kinetic energy but the conservation of $q^2$, the potential enstrophy for a horizontally non-divergent flow ($i.e.$ when $\chi$=$0$). Indeed, using the symmetry or skew symmetry properties of the operators ( \autoref{eq:DOM_mi_adj} and \autoref{eq:DOM_di_adj}), it can be shown that: The scheme does not allow but the conservation of the total kinetic energy but the conservation of $q^2$, the potential enstrophy for a horizontally non-divergent flow ($i.e.$ when $\chi$=$0$). Indeed, using the symmetry or skew symmetry properties of the operators ( \autoref{eq:DOM_mi_adj} and \autoref{eq:DOM_di_adj}), it can be shown that: \label{eq:C_1.1} \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 where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element. Indeed, using \autoref{eq:dynvor_ens}, the discrete form of the right hand side of \autoref{eq:C_1.1} can be transformed as follow: where $dv=e_1\,e_2\,e_3 \; di\,dj\,dk$ is the volume element. Indeed, using \autoref{eq:dynvor_ens}, the discrete form of the right hand side of \autoref{eq:C_1.1} can be transformed as follow: \begin{flalign*} &\int_D q \,\; \textbf{k} \cdot \frac{1} {e_3 } \nabla \times \end{aligned}   } \right. where the indices $i_p$ and $k_p$ take the following value: where the indices $i_p$ and $k_p$ take the following values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$, and the vorticity triads, ${^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $T$-point, are given by: This formulation does conserve the potential enstrophy for a horizontally non-divergent flow ($i.e.$ $\chi=0$). Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2}$, similar manipulation can be done for the 3 others. The discrete form of the right hand side of \autoref{eq:C_1.1} applied to this triad only can be transformed as follow: Let consider one of the vorticity triad, for example ${^{i}_j}\mathbb{Q}^{+1/2}_{+1/2}$, similar manipulation can be done for the 3 others. The discrete form of the right hand side of \autoref{eq:C_1.1} applied to this triad only can be transformed as follow: \begin{flalign*} All the numerical schemes used in NEMO are written such that the tracer content is conserved by the internal dynamics and physics (equations in flux form). For advection, only the CEN2 scheme ($i.e.$ $2^{nd}$ order finite different scheme) conserves the global variance of tracer. Nevertheless the other schemes ensure that the global variance decreases ($i.e.$ they are at least slightly diffusive). For diffusion, all the schemes ensure the decrease of the total tracer variance, except the iso-neutral operator. There is generally no strict conservation of mass, as the equation of state is non linear with respect to $T$ and $S$. In practice, the mass is conserved to a very high accuracy. All the numerical schemes used in NEMO are written such that the tracer content is conserved by the internal dynamics and physics (equations in flux form). For advection, only the CEN2 scheme ($i.e.$ $2^{nd}$ order finite different scheme) conserves the global variance of tracer. Nevertheless the other schemes ensure that the global variance decreases ($i.e.$ they are at least slightly diffusive). For diffusion, all the schemes ensure the decrease of the total tracer variance, except the iso-neutral operator. There is generally no strict conservation of mass, as the equation of state is non linear with respect to $T$ and $S$. In practice, the mass is conserved to a very high accuracy. % ------------------------------------------------------------------------------------------------------------- %       Advection Term Whatever the advection scheme considered it conserves of the tracer content as all the scheme are written in flux form. Indeed,  let $T$ be the tracer and $\tau_u$, $\tau_v$, and $\tau_w$ its interpolated values at velocity point (whatever the interpolation is), Whatever the advection scheme considered it conserves of the tracer content as all the scheme are written in flux form. Indeed, let $T$ be the tracer and its $\tau_u$, $\tau_v$, and $\tau_w$ interpolated values at velocity point (whatever the interpolation is), the conservation of the tracer content due to the advection tendency is obtained as follows: \begin{flalign*} \end{flalign*} The conservation of the variance of tracer due to the advection tendency can be achieved only with the CEN2 scheme, $i.e.$ when $\tau_u= \overline T^{\,i+1/2}$, $\tau_v= \overline T^{\,j+1/2}$, and $\tau_w= \overline T^{\,k+1/2}$. The conservation of the variance of tracer due to the advection tendency can be achieved only with the CEN2 scheme, $i.e.$ when $\tau_u= \overline T^{\,i+1/2}$, $\tau_v= \overline T^{\,j+1/2}$, and $\tau_w= \overline T^{\,k+1/2}$. It can be demonstarted as follows: \begin{flalign*} The discrete formulation of the horizontal diffusion of momentum ensures the conservation of potential vorticity and the horizontal divergence, and the dissipation of the square of these quantities ($i.e.$ enstrophy and the variance of the horizontal divergence) as well as the dissipation of the horizontal kinetic energy. In particular, when the eddy coefficients are horizontally uniform, it ensures a complete separation of vorticity and horizontal divergence fields, so that diffusion (dissipation) of vorticity (enstrophy) does not generate horizontal divergence (variance of the horizontal divergence) and \textit{vice versa}. These properties of the horizontal diffusion operator are a direct consequence of properties \autoref{eq:DOM_curl_grad} and \autoref{eq:DOM_div_curl}. When the vertical curl of the horizontal diffusion of momentum (discrete sense) is taken, the term associated with the horizontal gradient of the divergence is locally zero. The discrete formulation of the horizontal diffusion of momentum ensures the conservation of potential vorticity and the horizontal divergence, and the dissipation of the square of these quantities ($i.e.$ enstrophy and the variance of the horizontal divergence) as well as the dissipation of the horizontal kinetic energy. In particular, when the eddy coefficients are horizontally uniform, it ensures a complete separation of vorticity and horizontal divergence fields, so that diffusion (dissipation) of vorticity (enstrophy) does not generate horizontal divergence (variance of the horizontal divergence) and \textit{vice versa}. These properties of the horizontal diffusion operator are a direct consequence of properties \autoref{eq:DOM_curl_grad} and \autoref{eq:DOM_div_curl}. When the vertical curl of the horizontal diffusion of momentum (discrete sense) is taken, the term associated with the horizontal gradient of the divergence is locally zero. % ------------------------------------------------------------------------------------------------------------- \label{subsec:C.6.1} The lateral momentum diffusion term conserves the potential vorticity : The lateral momentum diffusion term conserves the potential vorticity: \begin{flalign*} &\int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times \label{subsec:C.6.3} The lateral momentum diffusion term dissipates the enstrophy when the eddy coefficients are horizontally uniform: The lateral momentum diffusion term dissipates the enstrophy when the eddy coefficients are horizontally uniform: \begin{flalign*} &\int\limits_D  \zeta \; \textbf{k} \cdot \nabla \times \label{subsec:C.6.4} When the horizontal divergence of the horizontal diffusion of momentum (discrete sense) is taken, the term associated with the vertical curl of the vorticity is zero locally, due to \autoref{eq:DOM_div_curl}. The resulting term conserves the $\chi$ and dissipates $\chi^2$ when the eddy coefficients are horizontally uniform. When the horizontal divergence of the horizontal diffusion of momentum (discrete sense) is taken, the term associated with the vertical curl of the vorticity is zero locally, due to \autoref{eq:DOM_div_curl}. The resulting term conserves the $\chi$ and dissipates $\chi^2$ when the eddy coefficients are horizontally uniform. \begin{flalign*} & \int\limits_D  \nabla_h \cdot \label{sec:C.7} As for the lateral momentum physics, the continuous form of the vertical diffusion of momentum satisfies several integral constraints. The first two are associated with the conservation of momentum and the dissipation of horizontal kinetic energy: As for the lateral momentum physics, the continuous form of the vertical diffusion of momentum satisfies several integral constraints. The first two are associated with the conservation of momentum and the dissipation of horizontal kinetic energy: \begin{align*} \int\limits_D   \frac{1} {e_3 }\; \frac{\partial } {\partial k} \end{align*} The first property is obvious. The second results from: The first property is obvious. The second results from: \begin{flalign*} \int\limits_D \end{flalign*} The vorticity is also conserved. Indeed: The vorticity is also conserved. Indeed: \begin{flalign*} \int \limits_D \end{flalign*} If the vertical diffusion coefficient is uniform over the whole domain, the enstrophy is dissipated, $i.e.$ If the vertical diffusion coefficient is uniform over the whole domain, the enstrophy is dissipated, $i.e.$ \begin{flalign*} \int\limits_D \zeta \, \textbf{k} \cdot \nabla \times &\left[  \frac{A_u^{\,vm}} {e_{3uw}} \delta_{k+1/2} \left[ \delta_{j+1/2} \left[ e_{1u}\,u \right] \right]  \right]  \biggr\}  &&\\ \end{flalign*} Using the fact that the vertical diffusion coefficients are uniform, and that in $z$-coordinate, the vertical scale factors do not depend on $i$ and $j$ so that: $e_{3f} =e_{3u} =e_{3v} =e_{3t}$ and $e_{3w} =e_{3uw} =e_{3vw}$, it follows: Using the fact that the vertical diffusion coefficients are uniform, and that in $z$-coordinate, the vertical scale factors do not depend on $i$ and $j$ so that: $e_{3f} =e_{3u} =e_{3v} =e_{3t}$ and $e_{3w} =e_{3uw} =e_{3vw}$, it follows: \begin{flalign*} \equiv A^{\,vm} \sum\limits_{i,j,k} \zeta \;\delta_k \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right) \right)\; dv = 0    &&&\\ \end{flalign*} and the square of the horizontal divergence decreases ($i.e.$ the horizontal divergence is dissipated) if the vertical diffusion coefficient is uniform over the whole domain: and the square of the horizontal divergence decreases ($i.e.$ the horizontal divergence is dissipated) if the vertical diffusion coefficient is uniform over the whole domain: \begin{flalign*} \label{sec:C.8} The numerical schemes used for tracer subgridscale physics are written such that the heat and salt contents are conserved (equations in flux form). Since a flux form is used to compute the temperature and salinity, the quadratic form of these quantities ($i.e.$ their variance) globally tends to diminish. The numerical schemes used for tracer subgridscale physics are written such that the heat and salt contents are conserved (equations in flux form). Since a flux form is used to compute the temperature and salinity, the quadratic form of these quantities ($i.e.$ their variance) globally tends to diminish. As for the advection term, there is conservation of mass only if the Equation Of Seawater is linear.
• ## NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_D.tex

 r9407 \newpage $\$\newline    % force a new ligne This appendix some on going consideration on algorithms used or planned to be used in \NEMO. This appendix some on going consideration on algorithms used or planned to be used in \NEMO. $\$\newline    % force a new ligne \label{sec:TRA_adv_ubs} The UBS advection scheme is an upstream biased third order scheme based on an upstream-biased parabolic interpolation. It is also known as Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics). For example, in the $i$-direction : The UBS advection scheme is an upstream biased third order scheme based on an upstream-biased parabolic interpolation. It is also known as Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics). For example, in the $i$-direction: \label{eq:tra_adv_ubs2} \tau _u^{ubs} = \left\{  \begin{aligned} - \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. By choosing this expression for $\tau "$ we consider a fourth order approximation of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$). where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. By choosing this expression for $\tau "$ we consider a fourth order approximation of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$). Alternative choice: introduce the scale factors: This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}. It is a relatively good compromise between accuracy and smoothness. It is not a \emph{positive} scheme meaning false extrema are permitted but the amplitude of such are significantly reduced over the centred second order method. Nevertheless it is not recommended to apply it to a passive tracer that requires positivity. The intrinsic diffusion of UBS makes its use risky in the vertical direction where the control of artificial diapycnal fluxes is of paramount importance. It has therefore been preferred to evaluate the vertical flux using the TVD scheme when \np{ln\_traadv\_ubs}\forcode{ = .true.}. For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds to a second order centred scheme is evaluated using the \textit{now} velocity (centred in time) while the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity (forward in time. This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. UBS and QUICK schemes only differ by one coefficient. Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can be controlled by vertical advection (not vertical diffusion which is usually solved using an implicit scheme). Computer time can be saved by using a time-splitting technique on vertical advection. This possibility have been implemented and validated in ORCA05-L301. It is not currently offered in the current reference version. NB 2 : In a forthcoming release four options will be proposed for the vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme , or  \textit{(b)} a TVD scheme, or  \textit{(c)} an interpolation based on conservative parabolic splines following \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme. NB 3 : It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows: This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}. It is a relatively good compromise between accuracy and smoothness. It is not a \emph{positive} scheme meaning false extrema are permitted but the amplitude of such are significantly reduced over the centred second order method. Nevertheless it is not recommended to apply it to a passive tracer that requires positivity. The intrinsic diffusion of UBS makes its use risky in the vertical direction where the control of artificial diapycnal fluxes is of paramount importance. It has therefore been preferred to evaluate the vertical flux using the TVD scheme when \np{ln\_traadv\_ubs}\forcode{ = .true.}. For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds to a second order centred scheme is evaluated using the \textit{now} velocity (centred in time) while the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity (forward in time). This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. UBS and QUICK schemes only differ by one coefficient. Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme. NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can be controlled by vertical advection (not vertical diffusion which is usually solved using an implicit scheme). Computer time can be saved by using a time-splitting technique on vertical advection. This possibility have been implemented and validated in ORCA05-L301. It is not currently offered in the current reference version. NB 2: In a forthcoming release four options will be proposed for the vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme, or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme. NB 3: It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows: \label{eq:tra_adv_ubs2} \tau _u^{ubs} = \left\{  \begin{aligned} \end{split} \autoref{eq:tra_adv_ubs2} has several advantages. First it clearly evidence that the UBS scheme is based on the fourth order scheme to which is added an upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order part as stated above using \autoref{eq:tra_adv_ubs}. Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient with is simply proportional to the velocity. \autoref{eq:tra_adv_ubs2} has several advantages. First it clearly evidences that the UBS scheme is based on the fourth order scheme to which is added an upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order part as stated above using \autoref{eq:tra_adv_ubs}. Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient which is simply proportional to the velocity. laplacian diffusion: with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$, $i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$ it comes : it comes: \label{eq:tra_ldf_lap} \begin{split} \end{split} if the velocity is uniform ($i.e.$ $|u|=cst$) and choosing $\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]$ if the velocity is uniform ($i.e.$ $|u|=cst$) and choosing $\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]$ sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$): \label{sec:LF} We adopt the following semi-discrete notation for time derivative. Given the values of a variable $q$ at successive time step, the time derivation and averaging operators at the mid time step are: We adopt the following semi-discrete notation for time derivative. Given the values of a variable $q$ at successive time step, the time derivation and averaging operators at the mid time step are: \begin{subequations} \label{eq:dt_mt} \begin{align} \end{align} \end{subequations} As for space operator, the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\rdt/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$ , respectively. As for space operator, the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\rdt/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$, respectively. The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as: =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt} Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$ as it can be found sometime in literature. The leap-Frog time stepping is a second order centered scheme. As such it respects the quadratic invariant in integral forms, $i.e.$ the following continuous property, Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$ as it can be found sometimes in literature. The leap-Frog time stepping is a second order centered scheme. As such it respects the quadratic invariant in integral forms, $i.e.$ the following continuous property, \label{eq:Energy} \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} =  \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) , is satisfied in discrete form. Indeed, is satisfied in discrete form. Indeed, \begin{split} \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) \end{split} NB here pb of boundary condition when applying the adjoin! In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition (equivalently of the boundary value of the integration by part). In time this boundary condition is not physical and \textbf{add something here!!!} NB here pb of boundary condition when applying the adjoint! In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition (equivalently of the boundary value of the integration by part). In time this boundary condition is not physical and \textbf{add something here!!!} \subsection{Griffies iso-neutral diffusion operator} Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} scheme, but is formulated within the \NEMO framework ($i.e.$ using scale factors rather than grid-size and having a position of $T$-points that is not necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}). In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO, the off-diagonal terms of the small angle diffusion tensor contain several double spatial averages of a gradient, for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$. It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer allows for the presence of grid point oscillation structures that will be invisible to the operator. These structures are \textit{computational modes}. They will not be damped by the iso-neutral operator, and even possibly amplified by it. In other word, the operator applied to a tracer does not warranties the decrease of its global average variance. To circumvent this, we have introduced a smoothing of the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers ($i.e.$ they enter the computation of density), but it does not work for a passive tracer.   \citep{Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that nicely solve the problem. The idea is to get rid of combinations of an averaged in one direction combined with a derivative in the same direction by considering triads. For example in the (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows: Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} scheme, but is formulated within the \NEMO framework ($i.e.$ using scale factors rather than grid-size and having a position of $T$-points that is not necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}). In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO, the off-diagonal terms of the small angle diffusion tensor contain several double spatial averages of a gradient, for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$. It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer allows for the presence of grid point oscillation structures that will be invisible to the operator. These structures are \textit{computational modes}. They will not be damped by the iso-neutral operator, and even possibly amplified by it. In other word, the operator applied to a tracer does not warranties the decrease of its global average variance. To circumvent this, we have introduced a smoothing of the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers ($i.e.$ they enter the computation of density), but it does not work for a passive tracer. \citep{Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that nicely solve the problem. The idea is to get rid of combinations of an averaged in one direction combined with a derivative in the same direction by considering triads. For example in the (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows: \label{eq:Gf_triads} _i^k \mathbb{T}_{i_p}^{k_p} (T) \right) where the indices $i_p$ and $k_p$ define the four triads and take the following value: $i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$, $b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells, where the indices $i_p$ and $k_p$ define the four triads and take the following value: $i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$, $b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells, $A_i^k$ is the lateral eddy diffusivity coefficient defined at $T$-point, and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad : and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad: \label{eq:Gf_slopes} _i^k \mathbb{R}_{i_p}^{k_p} {\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] } Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of multiplying the temperature derivative by $\alpha$ and the salinity derivative by $\beta$. This is more efficient as the ratio $\alpha / \beta$ can to be evaluated directly. Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$. This choice has been motivated by the decrease of tracer variance and the presence of partial cell at the ocean bottom (see \autoref{apdx:Gf_operator}). Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of multiplying the temperature derivative by $\alpha$ and the salinity derivative by $\beta$. This is more efficient as the ratio $\alpha / \beta$ can to be evaluated directly. Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$. This choice has been motivated by the decrease of tracer variance and the presence of partial cell at the ocean bottom (see \autoref{apdx:Gf_operator}). %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!ht] \begin{center} \includegraphics[width=0.70\textwidth]{Fig_ISO_triad} \caption{  \protect\label{fig:ISO_triad} Triads used in the Griffies's like iso-neutral diffision scheme for $u$-component (upper panel) and $w$-component (lower panel).} \caption{  \protect\label{fig:ISO_triad} Triads used in the Griffies's like iso-neutral diffision scheme for $u$-component (upper panel) and $w$-component (lower panel).} \end{center} \end{figure} The four iso-neutral fluxes associated with the triads are defined at $T$-point. They take the following expression : They take the following expression: \begin{flalign} \label{eq:Gf_fluxes} \begin{split} \end{flalign} The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}): The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}): \begin{flalign} \label{eq:iso_flux} \textbf{F}_{iso}(T) %    \end{pmatrix} \end{flalign} resulting in a iso-neutral diffusion tendency on temperature given by the divergence of the sum of all the four triad fluxes : resulting in a iso-neutral diffusion tendency on temperature given by the divergence of the sum of all the four triad fluxes: \label{eq:Gf_operator} D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{ where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. This expression of the iso-neutral diffusion has been chosen in order to satisfy the following six properties: This expression of the iso-neutral diffusion has been chosen in order to satisfy the following six properties: \begin{description} \item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction : \item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction: \label{eq:Gf_property1a} D_l^T = \frac{1}{b_T}  \ \delta_{i} \item[$\bullet$ implicit treatment in the vertical]  In the diagonal term associated with the vertical divergence of the iso-neutral fluxes (i.e. the term associated with a second order vertical derivative) appears only tracer values associated with a single water column. This is of paramount importance since it means that the implicit in time algorithm for solving the vertical diffusion equation can be used to evaluate this term. It is a necessity since the vertical eddy diffusivity associated with this term, \item[$\bullet$ implicit treatment in the vertical] In the diagonal term associated with the vertical divergence of the iso-neutral fluxes (i.e. the term associated with a second order vertical derivative) appears only tracer values associated with a single water column. This is of paramount importance since it means that the implicit in time algorithm for solving the vertical diffusion equation can be used to evaluate this term. It is a necessity since the vertical eddy diffusivity associated with this term, \sum_{\substack{i_p, \,k_p}} \left\{ can be quite large. \item[$\bullet$ pure iso-neutral operator]  The iso-neutral flux of locally referenced potential density is zero, $i.e.$ \item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of locally referenced potential density is zero, $i.e.$ \begin{align} \label{eq:Gf_property2} \begin{matrix} \end{matrix} \end{align} This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ and the definition of the triads' slopes \autoref{eq:Gf_slopes}. \item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the total tracer content, $i.e.$ This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ and the definition of the triads' slopes \autoref{eq:Gf_slopes}. \item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the total tracer content, $i.e.$ \label{eq:Gf_property1} \sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0 This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form. \item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does not increase the total tracer variance, $i.e.$ This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form. \item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does not increase the total tracer variance, $i.e.$ \label{eq:Gf_property1} \sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0 The property is demonstrated in the \autoref{apdx:Gf_operator}. It is a key property for a diffusion term. It means that the operator is also a dissipation term, $i.e.$ it is a sink term for the square of the quantity on which it is applied. It therfore ensure that, when the diffusivity coefficient is large enough, the field on which it is applied become free of grid-point noise. \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint, $i.e.$ The property is demonstrated in the \autoref{apdx:Gf_operator}. It is a key property for a diffusion term. It means that the operator is also a dissipation term, $i.e.$ it is a sink term for the square of the quantity on which it is applied. It therfore ensures that, when the diffusivity coefficient is large enough, the field on which it is applied become free of grid-point noise. \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint, $i.e.$ \label{eq:Gf_property1} \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} In other word, there is no needs to develop a specific routine from the adjoint of this operator. We just have to apply the same routine. This properties can be demonstrated quite easily in a similar way the "non increase of tracer variance" property has been proved (see \autoref{apdx:Gf_operator}). In other word, there is no needs to develop a specific routine from the adjoint of this operator. We just have to apply the same routine. This properties can be demonstrated quite easily in a similar way the "non increase of tracer variance" property has been proved (see \autoref{apdx:Gf_operator}). \end{description} \subsection{Eddy induced velocity and skew flux formulation} When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), an additional advection term is added. The associated velocity is the so called eddy induced velocity, the formulation of which depends on the slopes of iso- neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo} is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), an additional advection term is added. The associated velocity is the so called eddy induced velocity, the formulation of which depends on the slopes of iso-neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo} is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. The eddy induced velocity is given by: \end{split} where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces. where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces. %%gm wrong: to be modified with 2 2D streamfunctions In other words, the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, which is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^* = \textbf{k} \times \nabla \phi$ In other words, the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, which is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^* = \textbf{k} \times \nabla \phi$. %%end gm A traditional way to implement this additional advection is to add it to the eulerian velocity prior to compute the tracer advection. This allows us to take advantage of all the advection schemes offered for the tracers (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers where \emph{positivity} of the advection scheme is of paramount importance. A traditional way to implement this additional advection is to add it to the eulerian velocity prior to compute the tracer advection. This allows us to take advantage of all the advection schemes offered for the tracers (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers where \emph{positivity} of the advection scheme is of paramount importance. % give here the expression using the triads. It is different from the one given in \autoref{eq:ldfeiv} % see just below a copy of this equation: \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form. It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows: \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form. It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows: \begin{flalign*} \begin{split} \end{split} \end{flalign*} and since the eddy induces velocity field is no-divergent, we end up with the skew form of the eddy induced advective fluxes: and since the eddy induces velocity field is no-divergent, we end up with the skew form of the eddy induced advective fluxes: \label{eq:eiv_skew_continuous} \textbf{F}_{eiv}^T = \begin{pmatrix} \end{pmatrix} The tendency associated with eddy induced velocity is then simply the divergence of the \autoref{eq:eiv_skew_continuous} fluxes. It naturally conserves the tracer content, as it is expressed in flux form and, as the advective form, it preserve the tracer variance. Another interesting property of \autoref{eq:eiv_skew_continuous} form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms: The tendency associated with eddy induced velocity is then simply the divergence of the \autoref{eq:eiv_skew_continuous} fluxes. It naturally conserves the tracer content, as it is expressed in flux form and, as the advective form, it preserves the tracer variance. Another interesting property of \autoref{eq:eiv_skew_continuous} form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms: \begin{flalign} \label{eq:eiv_skew+eiv_continuous} \textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &= \end{pmatrix} \end{flalign} The horizontal component reduces to the one use for an horizontal laplacian operator and the vertical one keep the same complexity, but not more. This property has been used to reduce the computational time \citep{Griffies_JPO98}, but it is not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to choose a discret form of  \autoref{eq:eiv_skew_continuous} which is consistent with the iso-neutral operator \autoref{eq:Gf_operator}. Using the slopes \autoref{eq:Gf_slopes} and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient), the resulting discret form is given by: The horizontal component reduces to the one use for an horizontal laplacian operator and the vertical one keeps the same complexity, but not more. This property has been used to reduce the computational time \citep{Griffies_JPO98}, but it is not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which is consistent with the iso-neutral operator \autoref{eq:Gf_operator}. Using the slopes \autoref{eq:Gf_slopes} and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient), the resulting discret form is given by: \label{eq:eiv_skew} \textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left( \begin{aligned} Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells. In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to $\mathbb{R}$ for the discret form to be exact. Such a choice of discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes. It also ensures the conservation of the tracer variance (see Appendix \autoref{apdx:eiv_skew}), $i.e.$ it does not include a diffusive component but is a "pure" advection term. In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to $\mathbb{R}$ for the discret form to be exact. Such a choice of discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes. It also ensures the conservation of the tracer variance (see Appendix \autoref{apdx:eiv_skew}), $i.e.$ it does not include a diffusive component but is a "pure" advection term. This part will be moved in an Appendix. The continuous property to be demonstrated is : The continuous property to be demonstrated is: \begin{align*} \int_D  D_l^T \; T \;dv   \leq 0 % \allowdisplaybreaks \intertext{The summation is done over all $i$ and $k$ indices, it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to regroup all the terms of the summation by triad at a ($i$,$k$) point. In other words, we regroup all the terms in the neighbourhood  that contain a triad at the same ($i$,$k$) indices. It becomes: } \intertext{The summation is done over all $i$ and $k$ indices, it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to regroup all the terms of the summation by triad at a ($i$,$k$) point. In other words, we regroup all the terms in the neighbourhood that contain a triad at the same ($i$,$k$) indices. It becomes: } % &\equiv -\sum_{i,k} % \allowdisplaybreaks \intertext{Then outing in factor the triad in each of the four terms of the summation and substituting the triads by their expression given in \autoref{eq:Gf_triads}. It becomes: } \intertext{Then outing in factor the triad in each of the four terms of the summation and substituting the triads by their expression given in \autoref{eq:Gf_triads}. It becomes: } % &\equiv -\sum_{i,k} The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities. Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$, then the previous demonstration would have let to: Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$, then the previous demonstration would have let to: \begin{align*} \int_D  S \; D_l^T  \;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\ &\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\} \end{align*} This means that the iso-neutral operator is self-adjoint. There is no need to develop a specific to obtain it. This means that the iso-neutral operator is self-adjoint. There is no need to develop a specific to obtain it. This have to be moved in an Appendix. The continuous property to be demonstrated is : The continuous property to be demonstrated is: \begin{align*} \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0 \end{matrix} \end{align*} The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs, they cancel out. Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ they cancel out with the neighbouring grid points. Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the $i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the tracer is preserved by the discretisation of the skew fluxes. The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs, they cancel out. Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ they cancel out with the neighbouring grid points. Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the $i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the tracer is preserved by the discretisation of the skew fluxes. \end{document}