# source:NEMO/trunk/doc/latex/NEMO/subfiles/chap_LDF.tex@10414

Last change on this file since 10414 was 10414, checked in by nicolasmartin, 2 years ago
• Comment \label commands on maths environments for unreferenced equations and adapt the unnumbered math container accordingly (mainly switch to shortanded LateX syntax with $...$)
• Add a code trick to build subfile with its own bibliography
• Fix right path for main LaTeX document in first line of subfiles (\documentclass[...]{subfiles})
• Rename abstract_foreword.tex to foreword.tex
• Fix some non-ASCII codes inserted here or there in LaTeX (0-9]*) • Made a first iteration on the indentation and alignement within math, figure and table environments to improve source code readability File size: 27.0 KB Line 1\documentclass[../main/NEMO_manual]{subfiles} 2 3\begin{document} 4 5% ================================================================ 6% Chapter Lateral Ocean Physics (LDF) 7% ================================================================ 8\chapter{Lateral Ocean Physics (LDF)} 9\label{chap:LDF} 10 11\minitoc 12 13\newpage 14 15The lateral physics terms in the momentum and tracer equations have been described in \autoref{eq:PE_zdf} and 16their discrete formulation in \autoref{sec:TRA_ldf} and \autoref{sec:DYN_ldf}). 17In this section we further discuss each lateral physics option. 18Choosing one lateral physics scheme means for the user defining, 19(1) the type of operator used (laplacian or bilaplacian operators, or no lateral mixing term); 20(2) the direction along which the lateral diffusive fluxes are evaluated 21(model level, geopotential or isopycnal surfaces); and 22(3) the space and time variations of the eddy coefficients. 23These three aspects of the lateral diffusion are set through namelist parameters 24(see the \textit{\ngn{nam\_traldf}} and \textit{\ngn{nam\_dynldf}} below). 25Note that this chapter describes the standard implementation of iso-neutral tracer mixing, 26and Griffies's implementation, which is used if \np{traldf\_grif}\forcode{ = .true.}, 27is described in Appdx\autoref{apdx:triad} 28 29%-----------------------------------nam_traldf - nam_dynldf-------------------------------------------- 30 31\nlst{namtra_ldf} 32 33\nlst{namdyn_ldf} 34%-------------------------------------------------------------------------------------------------------------- 35 36 37% ================================================================ 38% Direction of lateral Mixing 39% ================================================================ 40\section{Direction of lateral mixing (\protect\mdl{ldfslp})} 41\label{sec:LDF_slp} 42 43%%% 44\gmcomment{ 45 we should emphasize here that the implementation is a rather old one. 46 Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme. 47} 48 49A direction for lateral mixing has to be defined when the desired operator does not act along the model levels. 50This occurs when (a) horizontal mixing is required on tracer or momentum 51(\np{ln\_traldf\_hor} or \np{ln\_dynldf\_hor}) in s- or mixed s-z- coordinates, 52and (b) isoneutral mixing is required whatever the vertical coordinate is. 53This direction of mixing is defined by its slopes in the \textbf{i}- and \textbf{j}-directions at the face of 54the cell of the quantity to be diffused. 55For a tracer, this leads to the following four slopes: 56r_{1u}, r_{1w}, r_{2v}, r_{2w} (see \autoref{eq:tra_ldf_iso}), 57while for momentum the slopes are r_{1t}, r_{1uw}, r_{2f}, r_{2uw} for u and 58r_{1f}, r_{1vw}, r_{2t}, r_{2vw} for v. 59 60%gm% add here afigure of the slope in i-direction 61 62\subsection{Slopes for tracer geopotential mixing in the s-coordinate} 63 64In s-coordinates, geopotential mixing (i.e. horizontal mixing) r_1 and r_2 are the slopes between 65the geopotential and computational surfaces. 66Their discrete formulation is found by locally solving \autoref{eq:tra_ldf_iso} when 67the diffusive fluxes in the three directions are set to zero and T is assumed to be horizontally uniform, 68i.e. a linear function of z_T, the depth of a T-point. 69%gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient} 70 71\begin{equation} 72 \label{eq:ldfslp_geo} 73 \begin{aligned} 74 r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)} 75 \;\delta_{i+1/2}[z_t] 76 &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] \ \ \ \\ 77 r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)} 78 \;\delta_{j+1/2} [z_t] 79 &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] \ \ \ \\ 80 r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2} 81 &\approx \frac{1}{e_{1w}}\; \delta_{i+1/2}[z_{uw}\\ 82 r_{2w} &= \frac{1}{e_{2w}}\;\overline{\overline{\delta_{j+1/2}[z_t]}}^{\,j,\,k+1/2} 83 &\approx \frac{1}{e_{2w}}\; \delta_{j+1/2}[z_{vw}] 84 \end{aligned} 85\end{equation} 86 87%gm% caution I'm not sure the simplification was a good idea! 88 89These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}\forcode{ = .true.}rue, 90and either \np{ln\_traldf\_hor}\forcode{ = .true.} or \np{ln\_dynldf\_hor}\forcode{ = .true.}. 91 92\subsection{Slopes for tracer iso-neutral mixing} 93\label{subsec:LDF_slp_iso} 94 95In iso-neutral mixing r_1 and r_2 are the slopes between the iso-neutral and computational surfaces. 96Their formulation does not depend on the vertical coordinate used. 97Their discrete formulation is found using the fact that the diffusive fluxes of 98locally referenced potential density (i.e. in situ density) vanish. 99So, substituting T by \rho in \autoref{eq:tra_ldf_iso} and setting the diffusive fluxes in 100the three directions to zero leads to the following definition for the neutral slopes: 101 102\begin{equation} 103 \label{eq:ldfslp_iso} 104 \begin{split} 105 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac{\delta_{i+1/2}[\rho]} 106 {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,i+1/2,\,k}} \\ 107 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac{\delta_{j+1/2}\left[\rho \right]} 108 {\overline{\overline{\delta_{k+1/2}[\rho]}}^{\,j+1/2,\,k}} \\ 109 r_{1w} &= \frac{e_{3w}}{e_{1w}}\; 110 \frac{\overline{\overline{\delta_{i+1/2}[\rho]}}^{\,i,\,k+1/2}} 111 {\delta_{k+1/2}[\rho]} \\ 112 r_{2w} &= \frac{e_{3w}}{e_{2w}}\; 113 \frac{\overline{\overline{\delta_{j+1/2}[\rho]}}^{\,j,\,k+1/2}} 114 {\delta_{k+1/2}[\rho]} 115 \end{split} 116\end{equation} 117 118%gm% rewrite this as the explanation is not very clear !!! 119%In practice, \autoref{eq:ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:ldfslp_iso} using the in situ density, \rho, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the in situ density with depth. 120 121%By definition, neutral surfaces are tangent to the local in situ density \citep{McDougall1987}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters). 122 123%In the z-coordinate, the derivative of the \autoref{eq:ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} (T-level, which is the same as the u- and v-levels), so the in situ density can be used for its evaluation. 124 125As the mixing is performed along neutral surfaces, the gradient of \rho in \autoref{eq:ldfslp_iso} has to 126be evaluated at the same local pressure 127(which, in decibars, is approximated by the depth in meters in the model). 128Therefore \autoref{eq:ldfslp_iso} cannot be used as such, 129but further transformation is needed depending on the vertical coordinate used: 130 131\begin{description} 132 133\item[z-coordinate with full step: ] 134 in \autoref{eq:ldfslp_iso} the densities appearing in the i and j derivatives are taken at the same depth, 135 thus the in situ density can be used. 136 This is not the case for the vertical derivatives: \delta_{k+1/2}[\rho] is replaced by -\rho N^2/g, 137 where N^2 is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{McDougall1987} 138 (see \autoref{subsec:TRA_bn2}). 139 140\item[z-coordinate with partial step: ] 141 this case is identical to the full step case except that at partial step level, 142 the \emph{horizontal} density gradient is evaluated as described in \autoref{sec:TRA_zpshde}. 143 144\item[s- or hybrid s-z- coordinate: ] 145 in the current release of \NEMO, iso-neutral mixing is only employed for s-coordinates if 146 the Griffies scheme is used (\np{traldf\_grif}\forcode{ = .true.}; 147 see Appdx \autoref{apdx:triad}). 148 In other words, iso-neutral mixing will only be accurately represented with a linear equation of state 149 (\np{nn\_eos}\forcode{ = 1..2}). 150 In the case of a "true" equation of state, the evaluation of i and j derivatives in \autoref{eq:ldfslp_iso} 151 will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes. 152 153%gm% 154 Note: The solution for s-coordinate passes trough the use of different (and better) expression for 155 the constraint on iso-neutral fluxes. 156 Following \citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral diffusive flux of 157 locally referenced potential density, we stay in the T-S plane and consider the balance between 158 the neutral direction diffusive fluxes of potential temperature and salinity: 159 \[ 160 \alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S) 161
162  % gm{  where vector F is ....}
163
164This constraint leads to the following definition for the slopes:
165
166$167 % \label{eq:ldfslp_iso2} 168 \begin{split} 169 r_{1u} &= \frac{e_{3u}}{e_{1u}}\; \frac 170 {\alpha_u \;\delta_{i+1/2}[T] - \beta_u \;\delta_{i+1/2}[S]} 171 {\alpha_u \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,i+1/2,\,k} 172 -\beta_u \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,i+1/2,\,k} } \\ 173 r_{2v} &= \frac{e_{3v}}{e_{2v}}\; \frac 174 {\alpha_v \;\delta_{j+1/2}[T] - \beta_v \;\delta_{j+1/2}[S]} 175 {\alpha_v \;\overline{\overline{\delta_{k+1/2}[T]}}^{\,j+1/2,\,k} 176 -\beta_v \;\overline{\overline{\delta_{k+1/2}[S]}}^{\,j+1/2,\,k} } \\ 177 r_{1w} &= \frac{e_{3w}}{e_{1w}}\; \frac 178 {\alpha_w \;\overline{\overline{\delta_{i+1/2}[T]}}^{\,i,\,k+1/2} 179 -\beta_w \;\overline{\overline{\delta_{i+1/2}[S]}}^{\,i,\,k+1/2} } 180 {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]} \\ 181 r_{2w} &= \frac{e_{3w}}{e_{2w}}\; \frac 182 {\alpha_w \;\overline{\overline{\delta_{j+1/2}[T]}}^{\,j,\,k+1/2} 183 -\beta_w \;\overline{\overline{\delta_{j+1/2}[S]}}^{\,j,\,k+1/2} } 184 {\alpha_w \;\delta_{k+1/2}[T] - \beta_w \;\delta_{k+1/2}[S]} \\ 185 \end{split} 186$
187where $\alpha$ and $\beta$, the thermal expansion and saline contraction coefficients introduced in
188\autoref{subsec:TRA_bn2}, have to be evaluated at the three velocity points.
189In order to save computation time, they should be approximated by the mean of their values at $T$-points
190(for example in the case of $\alpha$:
191$\alpha_u=\overline{\alpha_T}^{i+1/2}$$\alpha_v=\overline{\alpha_T}^{j+1/2}$ and
192$\alpha_w=\overline{\alpha_T}^{k+1/2}$).
193
194Note that such a formulation could be also used in the $z$-coordinate and $z$-coordinate with partial steps cases.
195
196\end{description}
197
198This implementation is a rather old one.
199It is similar to the one proposed by Cox [1987], except for the background horizontal diffusion.
200Indeed, the Cox implementation of isopycnal diffusion in GFDL-type models requires
201a minimum background horizontal diffusion for numerical stability reasons.
202To overcome this problem, several techniques have been proposed in which the numerical schemes of
203the ocean model are modified \citep{Weaver_Eby_JPO97, Griffies_al_JPO98}.
204Griffies's scheme is now available in \NEMO if \np{traldf\_grif\_iso} is set true; see Appdx \autoref{apdx:triad}.
205Here, another strategy is presented \citep{Lazar_PhD97}:
206a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of
207grid point noise generated by the iso-neutral diffusion operator (\autoref{fig:LDF_ZDF1}).
208This allows an iso-neutral diffusion scheme without additional background horizontal mixing.
209This technique can be viewed as a diffusion operator that acts along large-scale
210(2~$\Delta$x) \gmcomment{2deltax doesnt seem very large scale} iso-neutral surfaces.
211The diapycnal diffusion required for numerical stability is thus minimized and its net effect on the flow is quite small when compared to the effect of an horizontal background mixing.
212
213Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,
214contrary to the \citet{Griffies_al_JPO98} operator which has that property.
215
216%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
217\begin{figure}[!ht]
218  \begin{center}
219    \includegraphics[width=0.70\textwidth]{Fig_LDF_ZDF1}
220    \caption {
221      \protect\label{fig:LDF_ZDF1}
222      averaging procedure for isopycnal slope computation.
223    }
224  \end{center}
225\end{figure}
226%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
227
229%First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes.
230%Second, numerical stability issues also require a bound on slopes.
231%Third, the question of boundary condition specified on slopes...
232
233%from griffies: chapter 13.1....
234
235
236
237% In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04},
238% the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly
239% to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the
240% surface motivates this flattening of isopycnals near the surface).
241
242For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also be bounded by
243$1/100$ everywhere.
244This constraint is applied in a piecewise linear fashion, increasing from zero at the surface to
245$1/100$ at $70$ metres and thereafter decreasing to zero at the bottom of the ocean
246(the fact that the eddies "feel" the surface motivates this flattening of isopycnals near the surface).
247
248%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
249\begin{figure}[!ht]
250  \begin{center}
251    \includegraphics[width=0.70\textwidth]{Fig_eiv_slp}
252    \caption{
253      \protect\label{fig:eiv_slp}
254      Vertical profile of the slope used for lateral mixing in the mixed layer:
255      \textit{(a)} in the real ocean the slope is the iso-neutral slope in the ocean interior,
256      which has to be adjusted at the surface boundary
257      (i.e. it must tend to zero at the surface since there is no mixing across the air-sea interface:
258      wall boundary condition).
259      Nevertheless, the profile between the surface zero value and the interior iso-neutral one is unknown,
260      and especially the value at the base of the mixed layer;
261      \textit{(b)} profile of slope using a linear tapering of the slope near the surface and
262      imposing a maximum slope of 1/100;
263      \textit{(c)} profile of slope actually used in \NEMO: a linear decrease of the slope from
264      zero at the surface to its ocean interior value computed just below the mixed layer.
265      Note the huge change in the slope at the base of the mixed layer between \textit{(b)} and \textit{(c)}.
266    }
267  \end{center}
268\end{figure}
269%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
270
271\colorbox{yellow}{add here a discussion about the flattening of the slopes, vs tapering the coefficient.}
272
273\subsection{Slopes for momentum iso-neutral mixing}
274
275The iso-neutral diffusion operator on momentum is the same as the one used on tracers but
276applied to each component of the velocity separately
277(see \autoref{eq:dyn_ldf_iso} in section~\autoref{subsec:DYN_ldf_iso}).
278The slopes between the surface along which the diffusion operator acts and the surface of computation
279($z$- or $s$-surfaces) are defined at $T$-, $f$-, and \textit{uw}- points for the $u$-component, and $T$-, $f$- and
280\textit{vw}- points for the $v$-component.
281They are computed from the slopes used for tracer diffusion,
282$i.e.$ \autoref{eq:ldfslp_geo} and \autoref{eq:ldfslp_iso} :
283
284285 % \label{eq:ldfslp_dyn} 286 \begin{aligned} 287 &r_{1t}\ \ = \overline{r_{1u}}^{\,i} &&& r_{1f}\ \ &= \overline{r_{1u}}^{\,i+1/2} \\ 288 &r_{2f} \ \ = \overline{r_{2v}}^{\,j+1/2} &&& r_{2t}\ &= \overline{r_{2v}}^{\,j} \\ 289 &r_{1uw} = \overline{r_{1w}}^{\,i+1/2} &&\ \ \text{and} \ \ & r_{1vw}&= \overline{r_{1w}}^{\,j+1/2} \\ 290 &r_{2uw}= \overline{r_{2w}}^{\,j+1/2} &&& r_{2vw}&= \overline{r_{2w}}^{\,j+1/2}\\ 291 \end{aligned} 292
293
294The major issue remaining is in the specification of the boundary conditions.
295The same boundary conditions are chosen as those used for lateral diffusion along model level surfaces,
296$i.e.$ using the shear computed along the model levels and with no additional friction at the ocean bottom
297(see \autoref{sec:LBC_coast}).
298
299
300% ================================================================
301% Lateral Mixing Operator
302% ================================================================
303\section{Lateral mixing operators (\protect\mdl{traldf}, \protect\mdl{traldf}) }
304\label{sec:LDF_op}
305
306
307
308% ================================================================
309% Lateral Mixing Coefficients
310% ================================================================
311\section{Lateral mixing coefficient (\protect\mdl{ldftra}, \protect\mdl{ldfdyn}) }
312\label{sec:LDF_coef}
313
314Introducing a space variation in the lateral eddy mixing coefficients changes the model core memory requirement,
315adding up to four extra three-dimensional arrays for the geopotential or isopycnal second order operator applied to
316momentum.
317Six CPP keys control the space variation of eddy coefficients: three for momentum and three for tracer.
318The three choices allow:
319a space variation in the three space directions (\key{traldf\_c3d}\key{dynldf\_c3d}),
320in the horizontal plane (\key{traldf\_c2d}\key{dynldf\_c2d}),
321or in the vertical only (\key{traldf\_c1d}\key{dynldf\_c1d}).
322The default option is a constant value over the whole ocean on both momentum and tracers.
323
324The number of additional arrays that have to be defined and the gridpoint position at which
325they are defined depend on both the space variation chosen and the type of operator used.
326The resulting eddy viscosity and diffusivity coefficients can be a function of more than one variable.
327Changes in the computer code when switching from one option to another have been minimized by
328introducing the eddy coefficients as statement functions
329(include file \hf{ldftra\_substitute} and \hf{ldfdyn\_substitute}).
330The functions are replaced by their actual meaning during the preprocessing step (CPP).
331The specification of the space variation of the coefficient is made in \mdl{ldftra} and \mdl{ldfdyn},
332or more precisely in include files \hf{traldf\_cNd} and \hf{dynldf\_cNd}, with N=1, 2 or 3.
333The user can modify these include files as he/she wishes.
334The way the mixing coefficient are set in the reference version can be briefly described as follows:
335
336\subsubsection{Constant mixing coefficients (default option)}
337When none of the \key{dynldf\_...} and \key{traldf\_...} keys are defined,
338a constant value is used over the whole ocean for momentum and tracers,
339which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist parameters.
340
341\subsubsection{Vertically varying mixing coefficients (\protect\key{traldf\_c1d} and \key{dynldf\_c1d})}
342The 1D option is only available when using the $z$-coordinate with full step.
343Indeed in all the other types of vertical coordinate,
344the depth is a 3D function of (\textbf{i},\textbf{j},\textbf{k}) and therefore,
345introducing depth-dependent mixing coefficients will require 3D arrays.
346In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced in which
347the surface value is \np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value,
348and the transition takes place around z=300~m with a width of 300~m
349($i.e.$ both the depth and the width of the inflection point are set to 300~m).
350This profile is hard coded in file \hf{traldf\_c1d}, but can be easily modified by users.
351
352\subsubsection{Horizontally varying mixing coefficients (\protect\key{traldf\_c2d} and \protect\key{dynldf\_c2d})}
353By default the horizontal variation of the eddy coefficient depends on the local mesh size and
354the type of operator used:
355\begin{equation}
356  \label{eq:title}
357  A_l = \left\{
358    \begin{aligned}
359      & \frac{\max(e_1,e_2)}{e_{max}} A_o^l           & \text{for laplacian operator } \\
360      & \frac{\max(e_1,e_2)^{3}}{e_{max}^{3}} A_o^l          & \text{for bilaplacian operator }
361    \end{aligned}
362  \right.
363\end{equation}
364where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain,
365and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer) namelist parameter.
366This variation is intended to reflect the lesser need for subgrid scale eddy mixing where
367the grid size is smaller in the domain.
368It was introduced in the context of the DYNAMO modelling project \citep{Willebrand_al_PO01}.
369Note that such a grid scale dependance of mixing coefficients significantly increase the range of stability of
370model configurations presenting large changes in grid pacing such as global ocean models.
371Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to
372large coefficient compare to the smallest grid size (see \autoref{sec:STP_forward_imp}),
373especially when using a bilaplacian operator.
374
375Other formulations can be introduced by the user for a given configuration.
376For example, in the ORCA2 global ocean model (see Configurations),
377the laplacian viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and
378decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}.
379This modification can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}.
380Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of
381ORCA2 and ORCA05 (see \&namcfg namelist).
382
383\subsubsection{Space varying mixing coefficients (\protect\key{traldf\_c3d} and \key{dynldf\_c3d})}
384
385The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases,
386$i.e.$ a hyperbolic tangent variation with depth associated with a grid size dependence of
387the magnitude of the coefficient.
388
389\subsubsection{Space and time varying mixing coefficients}
390
391There is no default specification of space and time varying mixing coefficient.
392The only case available is specific to the ORCA2 and ORCA05 global ocean configurations.
393It provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and
394eddy induced velocity (ORCA05) that depends on the local growth rate of baroclinic instability.
395This specification is actually used when an ORCA key and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined.
396
397The following points are relevant when the eddy coefficient varies spatially:
398
399(1) the momentum diffusion operator acting along model level surfaces is written in terms of curl and
400divergent components of the horizontal current (see \autoref{subsec:PE_ldf}).
401Although the eddy coefficient could be set to different values in these two terms,
402this option is not currently available.
403
404(2) with an horizontally varying viscosity, the quadratic integral constraints on enstrophy and on the square of
405the horizontal divergence for operators acting along model-surfaces are no longer satisfied
406(\autoref{sec:dynldf_properties}).
407
408(3) for isopycnal diffusion on momentum or tracers, an additional purely horizontal background diffusion with
409uniform coefficient can be added by setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0},
410a background horizontal eddy viscosity or diffusivity coefficient
411(namelist parameters whose default values are $0$).
412However, the technique used to compute the isopycnal slopes is intended to get rid of such a background diffusion,
413since it introduces spurious diapycnal diffusion (see \autoref{sec:LDF_slp}).
414
415(4) when an eddy induced advection term is used (\key{traldf\_eiv}),
416$A^{eiv}$, the eddy induced coefficient has to be defined.
417Its space variations are controlled by the same CPP variable as for the eddy diffusivity coefficient
418($i.e.$ \key{traldf\_cNd}).
419
420(5) the eddy coefficient associated with a biharmonic operator must be set to a \emph{negative} value.
421
422(6) it is possible to use both the laplacian and biharmonic operators concurrently.
423
424(7) it is possible to run without explicit lateral diffusion on momentum
425(\np{ln\_dynldf\_lap}\forcode{ = .?.}\np{ln\_dynldf\_bilap}\forcode{ = .false.}).
426This is recommended when using the UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{ = .true.},
427see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.
428
429% ================================================================
430% Eddy Induced Mixing
431% ================================================================
433\label{sec:LDF_eiv}
434
435%%gm  from Triad appendix  : to be incorporated....
436\gmcomment{
437  Values of iso-neutral diffusivity and GM coefficient are set as described in \autoref{sec:LDF_coef}.
438  If none of the keys \key{traldf\_cNd}, N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and
439  GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and \np{rn\_aeiv\_0}.
440  If 2D-varying coefficients are set with \key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal
441  scale factor according to \autoref{eq:title}
442  \footnote{
443    Except in global ORCA  $0.5^{\circ}$ runs with \key{traldf\_eiv},
444    where $A_l$ is set like $A_e$ but with a minimum vale of $100\;\mathrm{m}^2\;\mathrm{s}^{-1}$
445  }.
446  In idealised setups with \key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv} is set in
447  the global configurations with \key{traldf\_c2d}, a horizontally varying $A_e$ is instead set from
448  the Held-Larichev parameterisation
449  \footnote{
450    In this case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further reduced by a factor $|f/f_{20}|$,
451    where $f_{20}$ is the value of $f$ at $20^{\circ}$~N
452  } (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored unless it is zero.
453}
454
455When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),
457the formulation of which depends on the slopes of iso-neutral surfaces.
458Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces,
459$i.e.$ \autoref{eq:ldfslp_geo} is used in $z$-coordinates,
460and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $s$-coordinates.
461The eddy induced velocity is given by:
462\begin{equation}
463  \label{eq:ldfeiv}
464  \begin{split}
465    u^* & = \frac{1}{e_{2u}e_{3u}}\; \delta_k \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right]\\
466    v^* & = \frac{1}{e_{1u}e_{3v}}\; \delta_k \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right]\\
467    w^* & = \frac{1}{e_{1w}e_{2w}}\; \left\{ \delta_i \left[e_{2u} \, A_{uw}^{eiv} \; \overline{r_{1w}}^{\,i+1/2} \right] + \delta_j \left[e_{1v} \, A_{vw}^{eiv} \; \overline{r_{2w}}^{\,j+1/2} \right] \right\} \\
468  \end{split}
469\end{equation}
470where $A^{eiv}$ is the eddy induced velocity coefficient whose value is set through \np{rn\_aeiv},
471a \textit{nam\_traldf} namelist parameter.
472The three components of the eddy induced velocity are computed and
474This has been preferred to a separate computation of the advective trends associated with the eiv velocity,
475since it allows us to take advantage of all the advection schemes offered for the tracers
476(see \autoref{sec:TRA_adv}) and not just the $2^{nd}$ order advection scheme as in