Lateral Mixing Coefficient (ldftra.F90, ldfdyn.F90)

Introducing a space variation in the lateral eddy mixing coefficients changes the model core memory requirement, adding up to four extra three-dimensional arrays for the geopotential or isopycnal second order operator applied to momentum. Six CPP keys control the space variation of eddy coefficients: three for momentum and three for tracer. The three choices allow: a space variation in the three space directions (key_ traldf_c3d, key_ dynldf_c3d), in the horizontal plane (key_ traldf_c2d, key_ dynldf_c2d), or in the vertical only (key_ traldf_c1d, key_ dynldf_c1d). The default option is a constant value over the whole ocean on both momentum and tracers.

The number of additional arrays that have to be defined and the gridpoint position at which they are defined depend on both the space variation chosen and the type of operator used. The resulting eddy viscosity and diffusivity coefficients can be a function of more than one variable. Changes in the computer code when switching from one option to another have been minimized by introducing the eddy coefficients as statement functions (include file ldftra_substitute.h90 and ldfdyn_substitute.h90). The functions are replaced by their actual meaning during the preprocessing step (CPP). The specification of the space variation of the coefficient is made in ldftra.F90 and ldfdyn.F90, or more precisely in include files traldf_cNd.h90 and dynldf_cNd.h90, with N=1, 2 or 3. The user can modify these include files as he/she wishes. The way the mixing coefficient are set in the reference version can be briefly described as follows:

Constant Mixing Coefficients (default option)

When none of the key_dynldf_... and key_traldf_... keys are defined, a constant value is used over the whole ocean for momentum and tracers, which is specified through the rn_ahm_0_lap and rn_aht_0 namelist parameters.

Vertically varying Mixing Coefficients (key_ traldf_c1d and key_ dynldf_c1d)

The 1D option is only available when using the $z$-coordinate with full step. Indeed in all the other types of vertical coordinate, the depth is a 3D function of (i,j,k) and therefore, introducing depth-dependent mixing coefficients will require 3D arrays. In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced in which the surface value is rn_aht_0 (rn_ahm_0_lap), the bottom value is 1/4 of the surface value, and the transition takes place around z=300 m with a width of 300 m ($i.e.$ both the depth and the width of the inflection point are set to 300 m). This profile is hard coded in file traldf_c1d.h90, but can be easily modified by users.

Horizontally Varying Mixing Coefficients (key_ traldf_c2d and key_ dynldf_c2d)

By default the horizontal variation of the eddy coefficient depends on the local mesh size and the type of operator used:

$$A_l = \left\{ \begin{aligned}& \frac{\max(e_1,e_2)}{e_{max}} A_... ...}} A_o^l & \text{for bilaplacian operator } \end{aligned} \right.$$

where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain, and $A_o^l$ is the rn_ahm_0_lap (momentum) or rn_aht_0 (tracer) namelist parameter. This variation is intended to reflect the lesser need for subgrid scale eddy mixing where the grid size is smaller in the domain. It was introduced in the context of the DYNAMO modelling project [Willebrand et al., 2001]. Note that such a grid scale dependance of mixing coefficients significantly increase the range of stability of model configurations presenting large changes in grid pacing such as global ocean models. Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to large coefficient compare to the smallest grid size (see §3.3), especially when using a bilaplacian operator.

Other formulations can be introduced by the user for a given configuration. For example, in the ORCA2 global ocean model (see Configurations), the laplacian viscosity operator uses rn_ahm_0_lap = 4.10$^4$ m$^2$/s poleward of 20 north and south and decreases linearly to rn_aht_0 = 2.10$^3$ m$^2$/s at the equator [Delecluse and Madec, 2000, Madec et al., 1996]. This modification can be found in routine ldf_dyn_c2d_orca defined in ldfdyn_c2d.F90. Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of ORCA2 and ORCA05 (see &namcfg namelist).

Space Varying Mixing Coefficients (key_ traldf_c3d and key_ dynldf_c3d)

The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases, $i.e.$ a hyperbolic tangent variation with depth associated with a grid size dependence of the magnitude of the coefficient.

Space and Time Varying Mixing Coefficients

There are no default specifications of space and time varying mixing coefficient. One available case is specific to the ORCA2 and ORCA05 global ocean configurations. It provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and eddy induced velocity (ORCA05) that depends on the local growth rate of baroclinic instability. This specification is actually used when an ORCA key and both key_ traldf_eiv and key_ traldf_c2d are defined.

Smagorinsky viscosity (key_ dynldf_c3d and key_ dynldf_smag)

The key_ dynldf_smag key activates a 3D, time-varying viscosity that depends on the resolved motions. Following [Smagorinsky, 1993] the viscosity coefficient is set proportional to a local deformation rate based on the horizontal shear and tension, namely:

$$A_{m_{Smag}} = \left(\frac{{\sf CM_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert$$ (9.2)

where the deformation rate $\vert{D}\vert$ is given by

$$\vert{D}\vert=\sqrt{\left({\frac{\partial{u}} {\partial{x}}} -{\f... ...\frac{\partial{u}} {\partial{y}}} +{\frac{\partial{v}} {\partial{x}}}\right)^2}$$ (9.3)

and $L$ is the local gridscale given by:

$$L^2 = \frac{2{e_1}^2 {e_2}^2}{\left ( {e_1}^2 + {e_2}^2 \right )}$$ (9.4)

[Griffies and Hallberg, 2000] suggest values in the range 2.2 to 4.0 of the coefficient $\sf CM_{Smag}$ for oceanic flows. This value is set via the rn_cmsmag_1 namelist parameter. An additional parameter: rn_cmsh is included in NEMO for experimenting with the contribution of the shear term. A value of 1.0 (the default) calculates the deformation rate as above; a value of 0.0 will discard the shear term entirely.

For numerical stability, the calculated viscosity is bounded according to the following:

$${\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_ahm\_m\_lap\right ) \geq A_{m_{Smag}} \geq rn\_ahm\_0\_lap$$ (9.5)

with both parameters for the upper and lower bounds being provided via the indicated namelist parameters.

When $ln\_dynldf\_bilap = .true.$, a biharmonic version of the Smagorinsky viscosity is also available which sets a coefficient for the biharmonic viscosity as:

$$B_{m_{Smag}} = - \left(\frac{{\sf CM_{bSmag}}}{\pi}\right)^2 {L^4\over 8}\vert{D}\vert$$ (9.6)

which is bounded according to:

$${\rm MAX}\left (-{ L^4\over {64\Delta{t}}}, rn\_ahm\_m\_blp\right ) \leq B_{m_{Smag}} \leq rn\_ahm\_0\_blp$$ (9.7)

Note the reversal of the inequalities here because NEMO requires the biharmonic coefficients as negative numbers. $\sf CM_{bSmag}$ is set via the rn_cmsmag_2 namelist parameter and the bounding values have corresponding entries in the namelist too.

The current implementation in NEMO also allows for 3D, time-varying diffusivities to be set using the Smagorinsky approach. Users should note that this option is not recommended for many applications since diffusivities will tend to be largest near boundaries (where shears are greatest) leading to spurious upwellings ([Griffies, 2004], chapter 18.3.4). Nevertheless the option is there for those wishing to experiment. This choice requires both key_ traldf_c3d and key_ traldf_smag and uses the rn_chsmag ( ${\sf CH_{Smag}}$), rn_smsh and rn_aht_m namelist parameters in an analogous way to rn_cmsmag_1, rn_cmsh and rn_ahm_m_lap (see above) to set the diffusion coefficient:

$$A_{h_{Smag}} = \left(\frac{{\sf CH_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert$$ (9.8)

For numerical stability, the calculated diffusivity is bounded according to the following:

$${\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_aht\_m\right ) \geq A_{h_{Smag}} \geq rn\_aht\_0$$ (9.9)

$$

The following points are relevant when the eddy coefficient varies spatially:

(1) the momentum diffusion operator acting along model level surfaces is written in terms of curl and divergent components of the horizontal current (see §2.5.2). Although the eddy coefficient could be set to different values in these two terms, this option is not currently available.

(2) with an horizontally varying viscosity, the quadratic integral constraints on enstrophy and on the square of the horizontal divergence for operators acting along model-surfaces are no longer satisfied (Appendix C.7).

(3) for isopycnal diffusion on momentum or tracers, an additional purely horizontal background diffusion with uniform coefficient can be added by setting a non zero value of rn_ahmb_0 or rn_ahtb_0, a background horizontal eddy viscosity or diffusivity coefficient (namelist parameters whose default values are 0). However, the technique used to compute the isopycnal slopes is intended to get rid of such a background diffusion, since it introduces spurious diapycnal diffusion (see §9.2).

(4) when an eddy induced advection term is used (key_ traldf_eiv), $A^{eiv}$, the eddy induced coefficient has to be defined. Its space variations are controlled by the same CPP variable as for the eddy diffusivity coefficient ($i.e.$ key_traldf_cNd).

(5) the eddy coefficient associated with a biharmonic operator must be set to a negative value.

(6) it is possible to use both the laplacian and biharmonic operators concurrently.

(7) it is possible to run without explicit lateral diffusion on momentum (ln_dynldf_lap = ln_dynldf_bilap = false). This is recommended when using the UBS advection scheme on momentum (ln_dynadv_ubs = true, see 6.3.2) and can be useful for testing purposes.