source: trunk/DOC/BETA/Chapters/Chap_MISC.tex @ 817

% ================================================================
% 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. 
% also same problem as soon as important straits are not properly resolve by the mesh. For example in ORCA 1/4\r{}, this occors in several straits of the Indonesian archipelagos (Ombai, Lombok...).
% We describe here the three methods that can be used in NEMO to handle such unproperly resolved straits. 
%The communication consists of mixing tracers and mass/volume between non-adjacent water columns. 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 have been made to set them in a generic way. Nevertheless, 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 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]{./Figures/Fig_Gibraltar.pdf}
\includegraphics[width=0.80\textwidth]{./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 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}

%--------------------------------------------namcro--------------------------------------------------------
\namdisplay{namcro} 
%--------------------------------------------------------------------------------------------------------------

\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 zoom option (improperly because it is not associated with a change in model resolution) is a quite simple function that allows to 
perform a simulation over a sub-domain of an already defined configuration 
(i.e. without defining a new set of mesh, initial state and forcings). This 
option can be useful for testing the user setting of surface boundary 
conditions, or the initial ocean state of a huge ocean model configuration 
while having a small central memory requirement. It can also be used to 
easily test specific physics in a sub-domain (for example, test of the 
coupling between sea-ice and ocean model over the Arctic or Antarctic ocean 
as they are set in the global ocean version of OPA \citep{Madec1996}. 
In standard, this option does not include any specific treatment for 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 those 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 configuration, the user set \jp{jpiglo}=\jp{jpidta} \jp{jpjglo}=\jp{jpjdta }and \jp{jpk}=\jp{jpkdta}. When running the model 
over a sub-domain, the user have 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 dimension 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)

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!ht] \label{Fig_LBC_zoom}  \begin{center}
\includegraphics[width=0.90\textwidth]{./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 is a stand alone water column based on 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 parameterizations of unresolved turbulence (wind steering, langmuir circulation, skin layers) ; \textit{(c)} to compare the behaviour of different mixing vertical scheme  ; \textit{(d)} to perform sensitivity study to the vertical diffusion on a particular point of the ocean global domain ; \textit{(d)} to access to specific diagnostics, aside from the standard model variables,  because of having small in core memory requirement.

The methodology is based on the use of the zoom functionality (see \S\ref{MISC_zoom}), and the addition of some specific routines. There is no need of defining a new set of mesh, bathymetry, initial state and forcing, as 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} 
%--------------------------------------------------------------------------------------------------------------

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 consideration (less than a 
few hours), this is usually requires a large elapse time. In order overcome 
this problem, \citet{Bryan1984} introduces a technique that allows to accelerate 
the spin up to the equilibrium. It consists in using 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 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. Both are settled from \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 analysed 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 there is a wider range of instability. This technique is efficient for 
searching 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 keep this technique when using 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 
$\Delta \tilde {t}$ implies that the heat and salt contents are no more 
conserved due to the vertical coupling of the ocean level through both 
advection and diffusion.


% ================================================================
% Model optimisation, Control Print and Benchmark
% ================================================================
\section{Model optimisation, Control Print and Benchmark}
\label{MISC_opt}
%--------------------------------------------namdom-------------------------------------------------------
\namdisplay{nam_ctl} 
%--------------------------------------------------------------------------------------------------------------

Three points to be described here:

$\bullet$ Vector and memory optimisation:

 \key{vectopt\_loop} enable the internal loop collapse, a very efficient way to increase the length of vector 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 by incresing the CPU time requirement, as it suppress intermediate computation saved in in-core memory. This possibility have not been intensively used. In fact up to now, it only concern the TKE physics, in which,  when \key{vectopt\_memory} is defined, the coefficient used for horizontal smoothing of $A_v^T$ and $A_v^m$ are no more computed once 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 inner domain averaged trends in all TRA, DYN LDF and ZDF modules. Very useful to detect the origin of an undesired change in model results. 

2- also \np{ln\_ctl} but using the nictl. njctl. namelist parameters to check the origin of differences between mono and multi processor

3- \key{esopa} (to be rename key\_nemo) : also a option of model management. When defined, this key force the activation of all options and CPP keys. For example, all the tracer and momentum advection scheme are called ! There is therefore no physical meaning associated with model results. In fact, this option allows both the compilator and the model run to go through all the Fortran lines of the model. This allows to check is there is obvious compilation or running bugs for CPP options, and running bugs for namelist options.

5- last digit comparison (\np{nbit\_cmp}). In MPP simulation, the computation of sum of the whole domain is performed as the sum over all processors of the sum of the inner domain of each processor. This double sum nevr give exactly the same result as a single sum over the whole domain, due to truncature differences. The "bit comparison" option has been introduced in order to be able to check that mono-processor and multi-processor give exactly the same results. 


$\bullet$  Benchmark (\np{nbench}). This option, defines a benchmark run based on GYRE configuration in which the resolution remains the same whatever the domain size is. This allows to run very large model domain by just changing the domain size (\jp{jpiglo}, \jp{jpjglo}) without adjusting neither the time-step nor the physical parametrisations. 


% ================================================================
% 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 to compute $\partial_t \psi$, the time evolution of the 
barotropic streamfunction. $\partial_t \psi$ is solution of an elliptic 
equation (I.2.4) for which two solvers are available, a 
Successive-Over-Relaxation (SOR) or a preconditioned conjugate gradient 
(PCG) \citep{Madec1988, Madec1990}. The PCG is a very efficient 
method for solving elliptic equations on vector computers. It is a fast and 
rather easy to use method, which is an attractive feature 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 the search of an optimal parameter as in the SOR method. 
Nevertheless, the SOR has been kept because it is a linear solver, a very 
useful property when using the adjoint model of OPA.

The surface pressure gradient is computed in \mdl{dynspg}. The 
default option is the use of PCG or SOR 
depending on \np{nsolv} (namelist parameter). 
At each time step the time derivative of the barotropic streamfunction is 
the solution of (II.2.3). 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 equivalent equation (II.2.3) 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} 
    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 further discussion]:

initialisation (evaluate a first guess from former 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 
accelerates significantly the convergence, but it has to be adjusted 
empirically for each model domain, except for an uniform grid where an 
analytical expression for \textit{$\omega $ }can be found \citep{Richtmyer1967}. This 
parameter is defined as \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 highly recommended 
to use a $\epsilon $ smaller or equal to $10^{-2}$ as larger values may leads 
to numerically induced basin scale barotropic oscillations, and to use a two 
or three order of magnitude smaller value in eddy resolving configuration. 
The precision of the solver is not only a numerical parameter, but 
influences the physics. Indeed, the zero change of kinetic energy due to the 
work of surface pressure forces in rigid-lid approximation is only achieved 
at the precision demanded on the solver ({\S}~C.1-e). The precision is 
specified by setting \textit{eps} (\textbf{namelist parameter}). In addition, two other tests are used to stop the iterative 
algorithm. They concern the number of iterations and the module 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 while 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 (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 for $\epsilon = 10^{-2}$, increasing to a few 
thousand for $\epsilon = 10^{-12}$.

The vectorization of the SOR algorithm is not straightforward. (VII.5.4) 
contains two linear recurrences on $i$ and $j$. This inhibits the vectorisation 
({\S}~IV.2-a). Therefore (VII.5.4a) 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}

If the three equations in (VII.5.6) are solved inside the same do-loop, 
(VII.5.4a) and (VII.5.6) are strictly equivalent. In the model they are 
solved successively over the whole domain. The convergence is slower but 
(VII.5.6a) and (VII.5.6b) are vector loops on $i$-index (the inner loop) and 
(VII.5.6c) is adapted to Cray vector computers by using a routine from the 
Cray library (namely the FOLR routine) to solve the first-order linear 
recurrence. 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 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 (VII.5.1) 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 the solution in the following 
iterative way: assuming that $\textbf{x}^n$ is obtained, $\textbf{x}^{n+1}$ is searched under the form $\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 \,{\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 choice of the descent vector $\textbf{d}^n$ depends 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 base 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, 
both the PCG algorithm and the whole model computation are stopped when the 
number of iteration, $nmax$, or the module of the right hand side exceeds a 
specified value (see {\S}~VII.5-a for further discussion). The precision and 
the maximum number of iteration are specified by setting $eps$ and $nmax$ (\textbf{namelist parameters}).

It can be demonstrated that the 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 all the more fast as the matrix is closed to 
identity (i.e. the eigen values are closed 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} 
\textbf{Q}^{-1} \textbf{A x} = \textbf{Q}^{-1} \textbf{b}
\end{equation}

The same algorithm can be used to solve (VII.5.7) 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 OPA, \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 (VII.5.8) is in 
fact given by (VII.5.2) and thus the matrix and right hand side are computed 
independently from the solver used.

% -------------------------------------------------------------------------------------------------------------
%       FETI
% -------------------------------------------------------------------------------------------------------------
\subsection{FETI}
\label{MISC_solfeti}

FETI is a powerfull solver that was developed by Marc Guyon (ref !!!).  It as been just converted from Fortan 77 to 90, but never successfully tested after that. Since then, nobody has found enough time to further investigate the implmenntationof FETI and debug it.

Its main advantaged is to save a lot of CPU time compared to SOR or PCG algorithm. Nevertheless, its main drawbacks is that the solution is depends 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 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 done to 
determine the barotropic streamfunction with the correct boundary 
conditions. This is detailed below.

The model does not recognise itself the islands which must be defined using 
bathymetry information, i.e. $mbathy$ array, equals $-1$ over the first island, $-2$ over 
the second, ... , $-N$ over the $N^{th}$ island (see {\S}~VII.2-b). The model 
determines the position of island grid-points and defines a closed contour 
around each island which will be used to compute the circulation around each 
island. The closed contour is formed of the ocean grid-points which are the 
closest to the island.

First, the island barotropic streamfunctions $\psi_n$ are computed 
using the PCG (they are solutions of (VII.5.1) with the right-hand side 
equals 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 define a 
specific right hand side for (VII.5.1)~). The precision of this computation 
can be very high since it is done once. The absolute precision, $epsisl$, and the 
maximum number of iteration, $nmisl$, are the \textbf{namelist parameters} used for that 
computation. Their typical values are $epsisl = 10^{-10}$ and $nmisl = 4000$. Then the island matrix A is computed from (I.2.8) and reversed. At 
each time step, $\psi_0$, the solution of (I.2.4) with $\psi_0 = 0$ along all coastlines, is computed using either SOR or PCG. It has to 
be noted that the first guess of this computation is defined as in (VII.5.3) 
except that $\partial_t \psi_0$ is used, not $\partial_t \psi$. Indeed we are 
computing $\partial_t \psi_0$ which is usually far 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$ using (I.2.9) 
in order 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 not 
to introduce changes in the model results (even at the computer precision) 
between a run performed with several restarts and the same run performed in 
one time. It has to be noticed that this requires that the restart file 
contains two consecutive time steps for all the prognostic variables, and 
that it is save in the same binary format as the one used by the computer to 
calculate (in particular, 32 bits binary IEEE format for this file must not 
be used). The output listing and file(s) are defined but should be checked 
and eventually adapted to the user's needs. The output listing is stored in 
the $ocean.output$ file. The information are printed all the way through 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 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 standard and available on-line diagnostics are described in {\S}III-12-c.

The default output 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 the graphic outputs from the 
model's results.
}

% -------------------------------------------------------------------------------------------------------------
%       Tracer/Dynamics Trends
% -------------------------------------------------------------------------------------------------------------
\subsection{Tracer/Dynamics Trends (\key{trdlmd}, \key{diatrdtra})}
\label{MISC_tratrd}

%to be udated this corresponds to OPA8
When \key{diatrddyn} and/or \key{diatrddyn} cpp 
variables are defined, each trends of the dynamics and/or temperature and 
salinity time evolution equations are stored in three-dimensional arrays 
just after their computation (i.e. at the end of each $dyn\cdots .F$ 
and/or $tra\cdots .F$ routines). These trends are then used in diagnostic 
routines $diadyn.F$ and $diatra.F$ respectively. In standard, these routines check the basin averaged properties of the momentum and tracer equations every \textit{ntrd } time-steps (\textbf{namelist parameter}). These routines 
are given as an example; they must be adapted by the user to his/her 
desiderata.

These two options imply the definition of several arrays in the in core 
memory, increasing quite sensitively the code memory requirements (search 
for \key{diatrddyn} or 
\key{diatrdtra}in \textit{common.h} file)

% -------------------------------------------------------------------------------------------------------------
%       On-line Floats trajectories
% -------------------------------------------------------------------------------------------------------------
\subsection{On-line Floats trajectories}
\label{FLO}
%--------------------------------------------namflo-------------------------------------------------------
\namdisplay{namflo} 
%--------------------------------------------------------------------------------------------------------------

\colorbox{yellow}{a description is to be added here}

% -------------------------------------------------------------------------------------------------------------
%       Other Diagnostics
% -------------------------------------------------------------------------------------------------------------
\subsection{Other Diagnostics}
\label{MISC_diag_others}

%To be updated  this corresponds to OPA 8

Aside from the standard model variables, some other diagnostics are computed 
on-line or can be added in the model. The available ready-to-add diagnostics 
can be found in DIA. Among the available diagnostics one can 
quote:

- 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})