% ================================================================ % Chapter Ñ Miscellaneous Topics % ================================================================ \chapter{Miscellaneous Topics (xxx)} \label{MISC} \minitoc % ================================================================ % Representation of Unresolved Straits % ================================================================ \section{Representation of Unresolved Straits} \label{MISC_strait} % In climate modeling, it is often necessary to allow water masses that are separated by land to exchange properties. This situation arises in models when the grid mesh is too coarse to resolve narrow passageways that in reality provide crucial connections between water masses. For example, coarse grid spacing typically closes off the Mediterranean from the Atlantic at the Straits of Gibraltar. In this case, it is important for climate models to include the effects of salty water entering the Atlantic from the Mediterranean. Likewise, it is important for the Mediterranean to replenish its supply of water from the Atlantic to balance the net evaporation occurring over the Mediterranean region. % Same problem occurs as soon as important straits are not properly resolved by the mesh. For example it arises in ORCA 1/4\r{} in several straits of the Indonesian archipelago (Ombai, Lombok...). % We describe here the three methods that can be used in \NEMO to handle such improperly resolved straits. %The problem is resolved here by allowing the mixing of tracers and mass/volume between non-adjacent water columns at nominated regions within the model. Momentum is not mixed. The scheme conserves total tracer content, and total volume (the latter in $z*$- or $s*$-coordinate), and maintains compatibility between the tracer and mass/volume budgets. % Note: such changes are so specific to a given configuration that no attempt has been made to set them in a generic way. Nevertheless, an example of how they can be set up is given in the ORCA2 configuration (\key{ORCA_R2}). % ------------------------------------------------------------------------------------------------------------- % Hand made geometry changes % ------------------------------------------------------------------------------------------------------------- \subsection{Hand made geometry changes} \label{MISC_strait_hand} $\bullet$ reduced scale factor, also called partially open face (Hallberg, personnal communication 2006) %also called partially open cell or partial closed cells $\bullet$ increase of the viscous boundary layer thickness by local increase of the fmask value at the coast \colorbox{yellow}{Add a short description of scale factor changes staff and fmask increase} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!tbp] \label{Fig_MISC_strait_hand} \begin{center} \includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_Gibraltar.pdf} \includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_Gibraltar2.pdf} \caption {Example of the Gibraltar strait defined in a 1\r{} x 1\r{} mesh. \textit{Top}: using partially open cells. The meridional scale factor at $v$-point is reduced on both sides of the strait to account for the real width of the strait (about 20 km). Note that the scale factors of the strait $T$-point remains unchanged. \textit{Bottom}: using viscous boundary layers. The four fmask parameters along the strait coastlines are set to a value larger than 4, $i.e.$ "strong" no-slip case (see Fig.\ref{Fig_LBC_shlat}) creating a large viscous boundary layer that allows a reduced transport through the strait.} \end{center} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> % ------------------------------------------------------------------------------------------------------------- % Cross Land Advection % ------------------------------------------------------------------------------------------------------------- \subsection{Cross Land Advection (\textit{tracla} module)} \label{MISC_strait_cla} %--------------------------------------------namcla-------------------------------------------------------- \namdisplay{namcla} %-------------------------------------------------------------------------------------------------------------- \colorbox{yellow}{Add a short description of CLA staff here or in lateral boundary condition chapter?} % ================================================================ % Closed seas % ================================================================ \section{Closed seas} \label{MISC_closea} % ================================================================ % Sub-Domain Functionality (\textit{nizoom, njzoom}, namelist parameters) % ================================================================ \section{Sub-Domain Functionality (\jp{jpizoom}, \jp{jpjzoom})} \label{MISC_zoom} The sub-domain functionality, also improperly called the zoom option (improperly because it is not associated with a change in model resolution) is a quite simple function that allows a simulation over a sub-domain of an already defined configuration ($i.e.$ without defining a new mesh, initial state and forcings). This option can be useful for testing the user settings of surface boundary conditions, or the initial ocean state of a huge ocean model configuration while having a small computer memory requirement. It can also be used to easily test specific physics in a sub-domain (for example, see \citep{Madec1996} for a test of the coupling used in the global ocean version of OPA between sea-ice and ocean model over the Arctic or Antarctic ocean, using a sub-domain). In the standard model, this option does not include any specific treatment for the ocean boundaries of the sub-domain: they are considered as artificial vertical walls. Nevertheless, it is quite easy to add a restoring term toward a climatology in the vicinity of such boundaries (see \S\ref{TRA_dmp}). In order to easily define a sub-domain over which the computation can be performed, the dimension of all input arrays (ocean mesh, bathymetry, forcing, initial state, ...) are defined as \jp{jpidta}, \jp{jpjdta} and \jp{jpkdta} (\mdl{par\_oce} module), while the computational domain is defined through \jp{jpiglo}, \jp{jpjglo} and \jp{jpk} (\mdl{par\_oce} module). When running the model over the whole domain, the user sets \jp{jpiglo}=\jp{jpidta} \jp{jpjglo}=\jp{jpjdta} and \jp{jpk}=\jp{jpkdta}. When running the model over a sub-domain, the user has to provide the size of the sub-domain, (\jp{jpiglo}, \jp{jpjglo}, \jp{jpkglo}), and the indices of the south western corner as \jp{jpizoom} and \jp{jpjzoom} in the \mdl{par\_oce} module (Fig.~\ref{Fig_LBC_zoom}). Note that a third set of dimensions exist, \jp{jpi}, \jp{jpj} and \jp{jpk} which is actually used to perform the computation. It is set by default to \jp{jpi}=\jp{jpjglo} and \jp{jpj}=\jp{jpjglo}, except for massively parallel computing where the computational domain is laid out on local processor memories following a 2D horizontal splitting. % (see {\S}IV.2-c) ref to the section to be updated %>>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[!ht] \label{Fig_LBC_zoom} \begin{center} \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_zoom.pdf} \caption {Position of a model domain compared to the data input domain when the zoom functionality is used.} \end{center} \end{figure} %>>>>>>>>>>>>>>>>>>>>>>>>>>>> % ================================================================ % 1D model functionality % ================================================================ \section{Water column model: 1D model (\key{cfg\_1d})} \label{MISC_1D} The 1D model option simulates a stand alone water column within the 3D \NEMO system. It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers or a biogeochemical model. It is set up by defining the \key{cfg\_1d} CPP key. This 1D model is a very useful tool \textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes ; \textit{(b)} to investigate suitable parameterisations of unresolved turbulence (wind steering, langmuir circulation, skin layers) ; \textit{(c)} to compare the behaviour of different vertical mixing schemes ; \textit{(d)} to perform sensitivity studies on the vertical diffusion at a particular point of the ocean global domain ; \textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model. The methodology is based on the use of the zoom functionality (see \S\ref{MISC_zoom}), with some extra routines. There is no need to define a new mesh, bathymetry, initial state or forcing, since the 1D model will use those of the configuration it is a zoom of. % ================================================================ % Accelerating the Convergence % ================================================================ \section{Accelerating the Convergence (\np{nn\_acc} = 1)} \label{MISC_acc} %--------------------------------------------namdom------------------------------------------------------- \namdisplay{namdom} %-------------------------------------------------------------------------------------------------------------- % Steven update just bellow... not sure it is what I vant to say %Searching an equilibrium state with an ocean model requires repeated very long time %integrations (each of a few thousand years for a global model). Due to the size of Searching an equilibrium state with an ocean model requires very long time integration (a few thousand years for a global model). Due to the size of the time step required for numerical stability (less than a few hours), this usually requires a large elapsed time. In order to overcome this problem, \citet{Bryan1984} introduces a technique that is intended to accelerate the spin up to equilibrium. It uses a larger time step in the thermodynamic evolution equations than in the dynamic evolution equations. It does not affect the equilibrium solution but modifies the trajectory to reach it. The acceleration of convergence option is used when \np{nn\_acc}=1. In that case, $\Delta t=rdt$ is the time step of dynamics while $\widetilde{\Delta t}=rdttra$ is the tracer time-step. Their values are set from the \np{rdt} and \np{rdttra} namelist parameters. The set of prognostic equations to solve becomes: \begin{equation} \label{Eq_acc} \begin{split} \frac{\partial \textbf{U}_h }{\partial t} &\equiv \frac{\textbf{U}_h ^{t+1}-\textbf{U}_h^{t-1} }{2\Delta t} = \ldots \\ \frac{\partial T}{\partial t} &\equiv \frac{T^{t+1}-T^{t-1}}{2 \widetilde{\Delta t}} = \ldots \\ \frac{\partial S}{\partial t} &\equiv \frac{S^{t+1} -S^{t-1}}{2 \widetilde{\Delta t}} = \ldots \\ \end{split} \end{equation} \citet{Bryan1984} has examined the consequences of this distorted physics. Free waves have a slower phase speed, their meridional structure is slightly modified, and the growth rate of baroclinically unstable waves is reduced but with a wider range of instability. This technique is efficient for searching for an equilibrium state in coarse resolution models. However its application is not suitable for many oceanic problems: it cannot be used for transient or time evolving problems (in particular, it is very questionable to use this technique when there is a seasonal cycle in the forcing fields), and it cannot be used in high-resolution models where baroclinically unstable processes are important. Moreover, the vertical variation of $\widetilde{ \Delta t}$ implies that the heat and salt contents are no longer conserved due to the vertical coupling of the ocean level through both advection and diffusion. % first mention of depth varying tilda-delta-t! and namelist parameter associated to that are to be described % ================================================================ % Model optimisation, Control Print and Benchmark % ================================================================ \section{Model Optimisation, Control Print and Benchmark} \label{MISC_opt} %--------------------------------------------namctl------------------------------------------------------- \namdisplay{namctl} %-------------------------------------------------------------------------------------------------------------- % \gmcomment{why not make these bullets into subsections?} Three issues to be described here: $\bullet$ Vector and memory optimisation: \key{vectopt\_loop} enables internal loop collapse, a very efficient way to increase the length of vector calculations and thus speed up the model on vector computers. % Add here also one word on NPROMA technique that has been found useless, since compiler have made significant progress during the last decade. % Add also one word on NEC specific optimisation (Novercheck option for example) \key{vectopt\_memory} has been introduced in order to reduce the memory requirement of the model. This is obviously done at the cost of increasing the CPU time requirement, since it suppress intermediate computations which would have been saved in in-core memory. This feature has not been intensively used. In fact, currently it is only implemented for the TKE physics, in which, when \key{vectopt\_memory} is defined, the coefficients used for horizontal smoothing of $A_v^T$ and $A_v^m$ are no longer computed once and for all. This reduces the memory requirement by three 2D arrays. $\bullet$ Control print %: describe here 4 things: 1- \np{ln\_ctl} : compute and print the trends averaged over the interior domain in all TRA, DYN, LDF and ZDF modules. This option is very helpful when diagnosing the origin of an undesired change in model results. 2- also \np{ln\_ctl} but using the nictl and njctl namelist parameters to check the source of differences between mono and multi processor runs. 3- \key{esopa} (to be rename key\_nemo) : which is another option for model management. When defined, this key forces the activation of all options and CPP keys. For example, all tracer and momentum advection schemes are called! There is therefore no physical meaning associated with the model results. However, this option forces both the compiler and the model to run through all the \textsc{Fortran} lines of the model. This allows the user to check for obvious compilation or execution errors with all CPP options, and errors in namelist options. 3- \key{esopa} (to be rename key\_nemo) : which is another option for model management. When defined, this key forces the activation of all options and CPP keys. For example, all tracer and momentum advection schemes are called! There is therefore no physical meaning associated with the model results. However, this option forces both the compiler and the model to run through all the Fortran lines of the model. This allows the user to check for obvious compilation or execution errors with all CPP options, and errors in namelist options. 4- last digit comparison (\np{nbit\_cmp}). In an MPP simulation, the computation of a sum over the whole domain is performed as the summation over all processors of each of their sums over their interior domains. This double sum never gives exactly the same result as a single sum over the whole domain, due to truncation differences. The "bit comparison" option has been introduced in order to be able to check that mono-processor and multi-processor runs give exactly the same results. $\bullet$ Benchmark (\np{nbench}). This option defines a benchmark run based on a GYRE configuration in which the resolution remains the same whatever the domain size. This allows a very large model domain to be used, just by changing the domain size (\jp{jpiglo}, \jp{jpjglo}) and without adjusting either the time-step or the physical parameterisations. % ================================================================ % Elliptic solvers (SOL) % ================================================================ \section{Elliptic solvers (SOL directory)} \label{MISC_sol} %--------------------------------------------namdom------------------------------------------------------- \namdisplay{namsol} %-------------------------------------------------------------------------------------------------------------- The computation of the surface pressure gradient with a rigid lid assumption requires the computation of $\partial_t \psi$, the time evolution of the barotropic streamfunction. $\partial_t \psi$ is the solution of an elliptic equation \eqref{Eq_PE_psi} for which three solvers are available: a Successive-Over-Relaxation scheme (SOR), a preconditioned conjugate gradient scheme(PCG) \citep{Madec1988, Madec1990} and a Finite Elements Tearing and Interconnecting scheme (FETI) \citep{Guyon_al_EP99, Guyon_al_CalPar99}. The PCG is a very efficient method for solving elliptic equations on vector computers. It is a fast and rather easy method to use; which are attractive features for a large number of ocean situations (variable bottom topography, complex coastal geometry, variable grid spacing, islands, open or cyclic boundaries, etc ...). It does not require a search for an optimal parameter as in the SOR method. However, the SOR has been retained because it is a linear solver, which is a very useful property when using the adjoint model of \NEMO. The FETI solver is a very efficient method on massively parallel computers. However, it has never been used since OPA 8.0. The current version in \NEMO should not even successfully go through the compilation phase. The surface pressure gradient is computed in \mdl{dynspg}. The solver used (PCG or SOR) depending on the value of \np{nsolv} (namelist parameter). At each time step the time derivative of the barotropic streamfunction is the solution of \eqref{Eq_psi_total}. Introducing the following coefficients: \begin{equation} \begin{aligned} &C_{i,j}^{NS} &&= \frac{e_{2v}(i,j)}{\left( H_v (i,j) e_{1v} (i,j)\right)}\\ &C_{i,j}^{EW} &&= \frac{e_{1u}(i,j)}{\left( H_u (i,j) e_{2u} (i,j)\right)}\\ &B_{i,j} &&= \delta_i \left( e_{2v}M_v \right) - \delta_j \left( e_{1u}M_u \right)\\ \end{aligned} \end{equation} the five-point finite difference equation \eqref{Eq_psi_total} can be rewritten as: \begin{multline} \label{Eq_solmat} C_{i+1,j}^{NS} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j} + \; C_{i,j+1}^{EW}{ \left( \frac{\partial \psi}{\partial t } \right) }_{i,j+1} + \; C_{i,j} ^{NS} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1,j} + \; C_{i,j} ^{EW}{ \left( \frac{\partial \psi}{\partial t } \right) }_{i,j-1}\\ -\left(C_{i+1,j}^{NS} + C_{i,j+1}^{EW} + C_{i,j}^{NS} + C_{i,j}^{EW} \right) { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} = B_{i,j} \end{multline} \eqref{Eq_solmat} is a linear symmetric system of equations. All the elements of the corresponding matrix \textbf{A} vanish except those of five diagonals. With the natural ordering of the grid points (i.e. from west to east and from south to north), the structure of \textbf{A} is block-tridiagonal with tridiagonal or diagonal blocks. \textbf{A} is a positive-definite symmetric matrix of size $(jpi \cdot jpj)^2$, and \textbf{B}, the right hand side of \eqref{Eq_solmat}, is a vector. % ------------------------------------------------------------------------------------------------------------- % Successive Over Relaxation % ------------------------------------------------------------------------------------------------------------- \subsection{Successive Over Relaxation \np{nsolv}=2} \label{MISC_solsor} Let us introduce the four cardinal coefficients: $A_{i,j}^S = C_{i,j}^{NS}/D_{i,j}\,$, $A_{i,j}^W = C_{i,j}^{EW}/D_{i,j}\,$, $A_{i,j}^E = C_{i,j+1}^{EW}/D_{i,j}\,$ and $A_{i,j}^N = C_{i+1,j}^{NS}/D_{i,j}\,$, and define $\tilde B_{i,j} = B_{i,j}/D_{i,j}\,$, where $D_{i,j} = C_{i,j}^{NS}+ C_{i+1,j}^{NS} + C_{i,j}^{EW} + C_{i,j+1}^{EW} $ (i.e. the diagonal of \textbf{A}). (VII.5.1) can be rewritten as: \begin{multline} \label{Eq_solmat_p} A_{i,j}^{N} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j } +\,A_{i,j}^{E} { \left( \frac{\partial \psi}{\partial t } \right) }_{i ,j+1} \\ +\,A_{i,j}^{S} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1 ,j } +\,A_{i,j}^{W} { \left( \frac{\partial \psi}{\partial t } \right) }_{i ,j-1} - { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} = \tilde B_{i,j} \end{multline} The SOR method used is an iterative method. Its algorithm can be summarised as follows [see \citet{Haltiner1980} for a further discussion]: initialisation (evaluate a first guess from previous time step computations) \begin{equation} \left( \frac{\partial \psi}{\partial t } \right)_{i,j}^0 = 2{\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^t} - {\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^{t-1}} \end{equation} iteration $n$, from $n=0$ until convergence, do : \begin{multline} R_{i,j}^n = A_{i,j}^{N} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j}^n +\,A_{i,j}^{E} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j+1} ^n \\ +\,A_{i,j}^{S} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1,j} ^{n+1} +\,A_{i,j}^{W} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j-1} ^{n+1} - { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j}^n - \tilde B_{i,j} \end{multline} \begin{equation} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} ^{n+1} = { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j} ^{n} + \omega \;R_{i,j}^n \end{equation} where \textit{$\omega $ }satisfies $1\leq \omega \leq 2$. An optimal value exists for \textit{$\omega$} which significantly accelerates the convergence, but it has to be adjusted empirically for each model domain (except for a uniform grid where an analytical expression for \textit{$\omega$} can be found \citep{Richtmyer1967}). The value of $\omega$ is set using \textit{sor}, a \textbf{namelist} parameter. The convergence test is of the form: \begin{equation} \delta = \frac{\sum\limits_{i,j}{R_{i,j}^n}{R_{i,j}^n}}{\sum\limits_{i,j}{ \tilde B_{i,j}^n}{\tilde B_{i,j}^n}} \leq \epsilon \end{equation} where $\epsilon$ is the absolute precision that is required. It is recommended that a value smaller or equal to $10^{-3}$ is used for $\epsilon$ since larger values may lead to numerically induced basin scale barotropic oscillations. In fact, for an eddy resolving configuration or in a filtered free surface case, a value three orders of magnitude smaller than this should be used. The precision is specified by setting \textit{eps} (\textbf{namelist parameter}). In addition, two other tests are used to halt the iterative algorithm. They involve the number of iterations and the modulus of the right hand side. If the former exceeds a specified value, \textit{nmax} (\textbf{namelist parameter}), or the latter is greater than $10^{15}$, the whole model computation is stopped and the last computed time step fields are saved in the standard output file. In both cases, this usually indicates that there is something wrong in the model configuration (an error in the mesh, the initial state, the input forcing, or the magnitude of the time step or of the mixing coefficients). A typical value of $nmax$ is a few hundred when $\epsilon = 10^{-6}$, increasing to a few thousand when $\epsilon = 10^{-12}$. The vectorization of the SOR algorithm is not straightforward. The scheme contains two linear recurrences on $i$ and $j$. This inhibits the vectorisation. Therefore it has been rewritten as: \begin{multline} R_{i,j}^n = A_{i,j}^{N} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j}^n +\,A_{i,j}^{E} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j+1} ^n \\ +\,A_{i,j}^{S} { \left( \frac{\partial \psi}{\partial t } \right) }_{i-1,j} ^{n} +\,A_{i,j}^{W} { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j-1} ^{n} - { \left( \frac{\partial \psi}{\partial t } \right) }_{i,j}^n - \tilde B_{i,j} \end{multline} \begin{equation} R_{i,j}^n = R_{i,j}^n - \omega \;A_{i,j}^{S}\; R_{i,j-1}^n \end{equation} \begin{equation} R_{i,j}^n = R_{i,j}^n - \omega \;A_{i,j}^{W}\; R_{i-1,j}^n \end{equation} The SOR method is very flexible and can be used under a wide range of conditions, including irregular boundaries, interior boundary points, etc. Proofs of convergence, etc. may be found in the standard numerical methods texts for partial differential equations. % ------------------------------------------------------------------------------------------------------------- % Preconditioned Conjugate Gradient % ------------------------------------------------------------------------------------------------------------- \subsection{Preconditioned Conjugate Gradient} \label{MISC_solpcg} \begin{flushright} (\textit{nbsfs=1}, \textbf{namelist parameter}) \end{flushright} \textbf{A} is a definite positive symmetric matrix, thus solving the linear system \eqref{Eq_solmat} is equivalent to the minimisation of a quadratic functional: \begin{equation*} \textbf{Ax} = \textbf{b} \leftrightarrow \textbf{x} =\text{inf}_{y} \,\phi (\textbf{y}) \quad , \qquad \phi (\textbf{y}) = 1/2 \langle \textbf{Ay},\textbf{y}\rangle - \langle \textbf{b},\textbf{y} \rangle \end{equation*} where $\langle , \rangle$ is the canonical dot product. The idea of the conjugate gradient method is to search for the solution in the following iterative way: assuming that $\textbf{x}^n$ has been obtained, $\textbf{x}^{n+1}$ is found from $\textbf {x}^{n+1}={\textbf {x}}^n+\alpha^n{\textbf {d}}^n$ which satisfies: \begin{equation*} {\textbf{ x}}^{n+1}=\text{inf} _{{\textbf{ y}}={\textbf{ x}}^n+\alpha^n \,{\textbf{ d}}^n} \,\phi ({\textbf{ y}})\;\;\Leftrightarrow \;\;\frac{d\phi }{d\alpha}=0 \end{equation*} and expressing $\phi (\textbf{y})$ as a function of \textit{$\alpha $}, we obtain the value that minimises the functional: \begin{equation*} \alpha ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{ A d}^n, \textbf{d}^n} \rangle \end{equation*} where $\textbf{r}^n = \textbf{b}-\textbf{A x}^n = \textbf{A} (\textbf{x}-\textbf{x}^n)$ is the error at rank $n$. The descent vector $\textbf{d}^n$ s chosen to be dependent on the error: $\textbf{d}^n = \textbf{r}^n + \beta^n \,\textbf{d}^{n-1}$. $\beta ^n$ is searched such that the descent vectors form an orthogonal basis for the dot product linked to \textbf{A}. Expressing the condition $\langle \textbf{A d}^n, \textbf{d}^{n-1} \rangle = 0$ the value of $\beta ^n$ is found: $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{r}^{n-1}, \textbf{r}^{n-1}} \rangle$. As a result, the errors $ \textbf{r}^n$ form an orthogonal base for the canonic dot product while the descent vectors $\textbf{d}^n$ form an orthogonal base for the dot product linked to \textbf{A}. The resulting algorithm is thus the following one: initialisation : \begin{equation*} \textbf{x}^0 = \left( \frac{\partial \psi}{\partial t } \right)_{i,j}^0 = 2{\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^t} - {\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^{t-1}} , \text{the initial guess } \end{equation*} \begin{equation*} \textbf{r}^0 = \textbf{d}^0 = \textbf{b} - \textbf{A x}^0 \end{equation*} \begin{equation*} \gamma_0 = \langle{ \textbf{r}^0 , \textbf{r}^0} \rangle \end{equation*} iteration $n,$ from $n=0$ until convergence, do : \begin{equation} \begin{split} \text{z}^n& = \textbf{A d}^n \\ \alpha_n &= \gamma_n \langle{ \textbf{z}^n , \textbf{d}^n} \rangle \\ \textbf{x}^{n+1} &= \textbf{x}^n + \alpha_n \,\textbf{d}^n \\ \textbf{r}^{n+1} &= \textbf{r}^n - \alpha_n \,\textbf{z}^n \\ \gamma_{n+1} &= \langle{ \textbf{r}^{n+1} , \textbf{r}^{n+1}} \rangle \\ \beta_{n+1} &= \gamma_{n+1}/\gamma_{n} \\ \textbf{d}^{n+1} &= \textbf{r}^{n+1} + \beta_{n+1}\; \textbf{d}^{n}\\ \end{split} \end{equation} The convergence test is: \begin{equation} \delta = \gamma_{n}\; / \langle{ \textbf{b} , \textbf{b}} \rangle \leq \epsilon \end{equation} where $\epsilon $ is the absolute precision that is required. As for the SOR algorithm, the whole model computation is stopped when the number of iterations, $nmax$, or the modulus of the right hand side of the convergence equation exceeds a specified value (see \S\ref{MISC_solsor} for a further discussion). The required precision and the maximum number of iterations allowed are specified by setting $eps$ and $nmax$ (\textbf{namelist parameters}). It can be demonstrated that the above algorithm is optimal, provides the exact solution in a number of iterations equal to the size of the matrix, and that the convergence rate is faster as the matrix is closer to the identity matrix, $i.e.$ its eigenvalues are closer to 1. Therefore, it is more efficient to solve a better conditioned system which has the same solution. For that purpose, we introduce a preconditioning matrix \textbf{Q} which is an approximation of \textbf{A} but much easier to invert than \textbf{A}, and solve the system: \begin{equation} \label{Eq_pmat} \textbf{Q}^{-1} \textbf{A x} = \textbf{Q}^{-1} \textbf{b} \end{equation} The same algorithm can be used to solve \eqref{Eq_pmat} if instead of the canonical dot product the following one is used: ${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, and if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ are substituted to \textbf{b} and \textbf{A} \citep{Madec1988}. In \NEMO, \textbf{Q} is chosen as the diagonal of \textbf{ A}, i.e. the simplest form for \textbf{Q} so that it can be easily inverted. In this case, the discrete formulation of \eqref{Eq_pmat} is in fact given by \eqref{Eq_solmat_p} and thus the matrix and right hand side are computed independently from the solver used. % ------------------------------------------------------------------------------------------------------------- % FETI % ------------------------------------------------------------------------------------------------------------- \subsection{FETI} \label{MISC_solfeti} FETI is a powerful solver that was developed by Marc Guyon \citep{Guyon_al_EP99, Guyon_al_CalPar99}. It has been converted from Fortan 77 to 90, but never successfully tested after that. Its main advantage is to save a lot of CPU time when compared to the SOR or PCG algorithms. However, its main drawback is that the solution is dependent on the domain decomposition chosen. Using a different number of processors, the solution is the same at the precision required, but not the same at the computer precision. This make it hard to debug. % ------------------------------------------------------------------------------------------------------------- % Boundary Conditions --- Islands % ------------------------------------------------------------------------------------------------------------- \subsection{Boundary Conditions --- Islands (\key{islands})} \label{MISC_solisl} The boundary condition used for both recommended solvers is that the time derivative of the barotropic streamfunction is zero along all the coastlines. When islands are present in the model domain, additional computations must be performed to determine the barotropic streamfunction with the correct boundary conditions. This is detailed below. The model does not have specialised code for islands. They must instead be identified to the solvers by the user via bathymetry information, i.e. the $mbathy$ array should equal $-1$ over the first island, $-2$ over the second, ... , $-N$ over the $N^{th}$ island. The model determines the position of island grid-points and defines a closed contour around each island which is used to compute the circulation around each one. The closed contour is formed from the ocean grid-points which are the closest to the island. First, the island barotropic streamfunctions $\psi_n$ are computed using the SOR or PCG scheme (they are solutions of \eqref{Eq_solmat} with the right-hand side equal to zero and with $\psi_n = 1$ along the island $n$ and $\psi_n = 0$ along the other coastlines) (Note that specifying $1$ as boundary condition on an island for $\psi$ is equivalent to defining a specific right hand side for \eqref{Eq_solmat}). The requested precision of this computation can be very high since it is only performed once. The absolute precision, $epsisl$, and the maximum number of iterations allowed, $nmisl$, are the \textbf{namelist parameters} used for this computation. Their typical values are $epsisl = 10^{-10}$ and $nmisl = 4000$. Then the island matrix A is computed from \eqref{Eq_PE_psi_isl_matrix} and inverted. At each time step, $\psi_0$, the solution of \eqref{Eq_solmat} with $\psi = 0$ along all coastlines, is computed using either SOR or PCG. (It should be noted that the first guess of this computation remains as usual except that $\partial_t \psi_0$ is used, instead of $\partial_t \psi$. Indeed we are computing $\partial_t \psi_0$ which is usually very different from $\partial_t \psi$.) Then, it is easy to find the time evolution of the barotropic streamfunction on each island and to deduce $\partial_t \psi$, and to use it to compute the surface pressure gradient. Note that the value of the barotropic streamfunction itself is also computed as the time integration of $\partial_t \psi$ for further diagnostics. % ================================================================ % Diagnostics % ================================================================ \section{Diagnostics } \label{MISC_diag} % ------------------------------------------------------------------------------------------------------------- % Standard Model Output % ------------------------------------------------------------------------------------------------------------- \subsection{Standard Model Output (default option or \key{dimg})} \label{MISC_iom} %to be updated with Seb documentation on the IO \amtcomment{ The model outputs are of three types: the restart file, the output listing, and the output file(s). The restart file is used internally by the code when the user wants to start the model with initial conditions defined by a previous simulation. It contains all the information that is necessary in order for there to be no changes in the model results (even at the computer precision) between a run performed with several restarts and the same run performed in one step. It should be noted that this requires that the restart file contain two consecutive time steps for all the prognostic variables, and that it is saved in the same binary format as the one used by the computer that is to read it (in particular, 32 bits binary IEEE format must not be used for this file). The output listing and file(s) are predefined but should be checked and eventually adapted to the user's needs. The output listing is stored in the $ocean.output$ file. The information is printed from within the code on the logical unit $numout$. To locate these prints, use the UNIX command "$grep -i numout^*$" in the source code directory. In the standard configuration, the user will find the model results in two output files containing values for every time-step where output is demanded: a \textbf{VO} file containing all the three dimensional fields in logical unit $numwvo$, and a \textbf{SO} file containing all the two dimensional (horizontal) fields in logical unit $numwso$. These outputs are defined in the $diawri.F$ routine. The default and user-selectable diagnostics are described in {\S}III-12-c. The default output (for all output files apart from the restart file) is 32 bits binary IEEE format, compatible with the Vairmer software (see the Climate Modelling and global Change team WEB server at CERFACS: http://www.cerfacs.fr). The model's reference directory also contains a visualisation tool based on \textbf{NCAR Graphics} (http://ngwww.ucar.edu). If a user has access to the NCAR software, he or she can copy the \textbf{LODMODEL/UTILS/OPADRA} directory from the reference and, following the \textbf{README}, create graphical outputs from the model's results. } % ------------------------------------------------------------------------------------------------------------- % Tracer/Dynamics Trends % ------------------------------------------------------------------------------------------------------------- \subsection[Tracer/Dynamics Trends (\key{trdlmd}, \textbf{key\_diatrd...})] {Tracer/Dynamics Trends (\key{trdlmd}, \key{diatrdtra}, \key{diatrddyn})} \label{MISC_tratrd} %to be udated this corresponds to OPA8 When \key{diatrddyn} and/or \key{diatrddyn} cpp variables are defined, each trend of the dynamics and/or temperature and salinity time evolution equations is stored in three-dimensional arrays just after their computation ($i.e.$ at the end of each $dyn\cdots .F90$ and/or $tra\cdots .F90$ routine). These trends are then used in diagnostic routines $diadyn.F90$ and $diatra.F90$ respectively. In the standard model, these routines check the basin averaged properties of the momentum and tracer equations every \textit{ntrd } time-steps (\textbf{namelist parameter}). These routines are supplied as an example; they must be adapted by the user to his/her requirements. These two options imply the creation of several extra arrays in the in-core memory, increasing quite seriously the code memory requirements. % ------------------------------------------------------------------------------------------------------------- % On-line Floats trajectories % ------------------------------------------------------------------------------------------------------------- \subsection{On-line Floats trajectories} \label{FLO} %--------------------------------------------namflo------------------------------------------------------- \namdisplay{namflo} %-------------------------------------------------------------------------------------------------------------- The on-line computation of floats adevected either by the three dimensional velocity field or constraint to remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project. The algorithm used is based on the work of \cite{Blanke_Raynaud_JPO97}. (see also the web site describing the off-line use of this marvellous diagnostic tool (http://stockage.univ-brest.fr/~grima/Ariane/). % ------------------------------------------------------------------------------------------------------------- % Other Diagnostics % ------------------------------------------------------------------------------------------------------------- \subsection{Other Diagnostics} \label{MISC_diag_others} %To be updated this mainly corresponds to OPA 8 Aside from the standard model variables, other diagnostics can be computed on-line or can be added to the model. The available ready-to-add diagnostics routines can be found in directory DIA. Among the available diagnostics are: - the mixed layer depth (based on a density criterion) (\mdl{diamxl}) - the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diamxl}) - the depth of the 20\r{}C isotherm (\mdl{diahth}) - the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth}) - the meridional heat and salt transports and their decomposition (\mdl{diamfl}) - the surface pressure (\mdl{diaspr})