Index: trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 5118)
+++ trunk/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 6289)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter � Miscellaneous Topics
+% Chapter ——— Miscellaneous Topics
 % ================================================================
 \chapter{Miscellaneous Topics}
@@ -34,12 +34,6 @@
 has been made to set them in a generic way. However, examples of how 
 they can be set up is given in the ORCA 2\deg and 0.5\deg configurations. For example, 
-for details of implementation in ORCA2, search:
-\vspace{-10pt}  
-\begin{alltt}  
-\tiny    
-\begin{verbatim}
-IF( cp_cfg == "orca" .AND. jp_cfg == 2 )
-\end{verbatim}  
-\end{alltt}
+for details of implementation in ORCA2, search: 
+\texttt{ IF( cp\_cfg == "orca" .AND. jp\_cfg == 2 ) }
 
 % -------------------------------------------------------------------------------------------------------------
@@ -80,17 +74,4 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-% -------------------------------------------------------------------------------------------------------------
-% Cross Land Advection 
-% -------------------------------------------------------------------------------------------------------------
-\subsection{Cross Land Advection (\mdl{tracla})}
-\label{MISC_strait_cla}
-%--------------------------------------------namcla--------------------------------------------------------
-\namdisplay{namcla} 
-%--------------------------------------------------------------------------------------------------------------
-
-\colorbox{yellow}{Add a short description of CLA staff here or in lateral boundary condition chapter?}
-Options are defined through the  \ngn{namcla} namelist variables.
-
-%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.  
 
 % ================================================================
@@ -208,56 +189,4 @@
 
 % ================================================================
-% Accelerating the Convergence 
-% ================================================================
-\section{Accelerating the Convergence (\np{nn\_acc} = 1)}
-\label{MISC_acc}
-%--------------------------------------------namdom-------------------------------------------------------
-\namdisplay{namdom} 
-%--------------------------------------------------------------------------------------------------------------
-
-Searching an equilibrium state with an global ocean model requires a very long time 
-integration period (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 tracer evolution equations than in the momentum evolution 
-equations. It does not affect the equilibrium solution but modifies the 
-trajectory to reach it.
-
-Options are defined through the  \ngn{namdom} namelist variables.
-The acceleration of convergence option is used when \np{nn\_acc}=1. In that case, 
-$\rdt=rn\_rdt$ is the time step of dynamics while $\widetilde{\rdt}=rdttra$ is the 
-tracer time-step. the former is set from the \np{rn\_rdt} namelist parameter while the latter
-is computed using a hyperbolic tangent profile and the following namelist parameters : 
-\np{rn\_rdtmin}, \np{rn\_rdtmax} and \np{rn\_rdth}. Those three parameters correspond 
-to the surface value the deep ocean value and the depth at which the transition occurs, respectively. 
-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\rdt} = \ldots \\ 
-\frac{\partial T}{\partial t} &\equiv \frac{T^{t+1}-T^{t-1}}{2 \widetilde{\rdt}} = \ldots \\ 
-\frac{\partial S}{\partial t} &\equiv \frac{S^{t+1} -S^{t-1}}{2 \widetilde{\rdt}} = \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{ \rdt}$ 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. Therefore \np{rn\_rdtmin} = \np{rn\_rdtmax} should be
-a more clever choice.
-
-
-% ================================================================
 % Accuracy and Reproducibility
 % ================================================================
@@ -370,13 +299,6 @@
 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! 
-Therefore the model results have no physical meaning. 
-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.
-
-4- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of 
+%%gm   to be removed both here and in the code
+3- last digit comparison (\np{nn\_bit\_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 
@@ -386,4 +308,5 @@
 %THIS is to be updated with the mpp_sum_glo  introduced in v3.3
 % nn_bit_cmp  today only check that the nn_cla = 0 (no cross land advection)
+%%gm end
 
 $\bullet$  Benchmark (\np{nn\_bench}). This option defines a benchmark run based on 
@@ -393,247 +316,8 @@
 or the physical parameterisations. 
 
-
-% ================================================================
-% Elliptic solvers (SOL)
-% ================================================================
-\section{Elliptic solvers (SOL)}
-\label{MISC_sol}
-%--------------------------------------------namdom-------------------------------------------------------
-\namdisplay{namsol} 
-%--------------------------------------------------------------------------------------------------------------
-
-When the filtered sea surface height option is used, the surface pressure gradient is 
-computed in \mdl{dynspg\_flt}. The force added in the momentum equation is solved implicitely.
-It is thus solution of an elliptic equation \eqref{Eq_PE_flt} for which two solvers are available: 
-a Successive-Over-Relaxation scheme (SOR) and a preconditioned conjugate gradient 
-scheme(PCG) \citep{Madec_al_OM88, Madec_PhD90}. The solver is selected trough the
-the value of \np{nn\_solv}   \ngn{namsol} namelist variable. 
-
-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, 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.
-
-At each time step, the time derivative of the sea surface height at time step $t+1$ 
-(or equivalently the divergence of the \textit{after} barotropic transport) that appears 
-in the filtering forced is the solution of the elliptic equation obtained from the horizontal 
-divergence of the vertical summation of \eqref{Eq_PE_flt}. 
-Introducing the following coefficients:
-\begin{equation}  \label{Eq_sol_matrix}
-\begin{aligned}
-&c_{i,j}^{NS}  &&= {2 \rdt }^2 \; \frac{H_v (i,j) \; e_{1v} (i,j)}{e_{2v}(i,j)}              \\
-&c_{i,j}^{EW} &&= {2 \rdt }^2 \; \frac{H_u (i,j) \; e_{2u} (i,j)}{e_{1u}(i,j)}            \\
-&b_{i,j} &&= \delta_i \left[ e_{2u}M_u \right] - \delta_j \left[ e_{1v}M_v \right]\ ,   \\
-\end{aligned}
-\end{equation}
-the resulting five-point finite difference equation is given by:
-\begin{equation}  \label{Eq_solmat}
-\begin{split}
-       c_{i+1,j}^{NS} D_{i+1,j}  + \;  c_{i,j+1}^{EW} D_{i,j+1}   
-  +   c_{i,j}    ^{NS} D_{i-1,j}   + \;  c_{i,j}    ^{EW} D_{i,j-1}                                          &    \\
-  -    \left(c_{i+1,j}^{NS} + c_{i,j+1}^{EW} + c_{i,j}^{NS} + c_{i,j}^{EW} \right)   D_{i,j}  &=  b_{i,j}
-\end{split}
-\end{equation}
-\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.
-
-Note that in the linear free surface case, the depth that appears in \eqref{Eq_sol_matrix}
-does not vary with time, and thus the matrix can be computed once for all. In non-linear free surface 
-(\key{vvl} defined) the matrix have to be updated at each time step.
-
-% -------------------------------------------------------------------------------------------------------------
-%       Successive Over Relaxation
-% -------------------------------------------------------------------------------------------------------------
-\subsection{Successive Over Relaxation (\np{nn\_solv}=2, \mdl{solsor})}
-\label{MISC_solsor}
-
-Let us introduce the four cardinal coefficients: 
-\begin{align*}
-a_{i,j}^S &= c_{i,j    }^{NS}/d_{i,j}     &\qquad  a_{i,j}^W &= c_{i,j}^{EW}/d_{i,j}       \\
-a_{i,j}^E &= c_{i,j+1}^{EW}/d_{i,j}    &\qquad   a_{i,j}^N &= c_{i+1,j}^{NS}/d_{i,j}   
-\end{align*}
-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 the matrix). \eqref{Eq_solmat} can be rewritten as:
-\begin{equation}  \label{Eq_solmat_p}
-\begin{split}
-a_{i,j}^{N}  D_{i+1,j} +\,a_{i,j}^{E}  D_{i,j+1} +\, a_{i,j}^{S}  D_{i-1,j} +\,a_{i,j}^{W} D_{i,j-1}  -  D_{i,j} = \tilde{b}_{i,j}
-\end{split}
-\end{equation}
-with $\tilde b_{i,j} = b_{i,j}/d_{i,j}$. \eqref{Eq_solmat_p} is the equation actually solved 
-with the SOR method. This method used is an iterative one. 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}
-D_{i,j}^0 = 2 \, D_{i,j}^t - D_{i,j}^{t-1}
-\end{equation}
-iteration $n$, from $n=0$ until convergence, do :
-\begin{equation} \label{Eq_sor_algo}
-\begin{split}
-R_{i,j}^n  = &a_{i,j}^{N} D_{i+1,j}^n       +\,a_{i,j}^{E}  D_{i,j+1} ^n	    	
-			+\, a_{i,j}^{S}  D_{i-1,j} ^{n+1}+\,a_{i,j}^{W} D_{i,j-1} ^{n+1}
-                 -  D_{i,j}^n - \tilde{b}_{i,j}                                           \\
-D_{i,j} ^{n+1}  = &D_{i,j} ^{n}   + \omega \;R_{i,j}^n     
-\end{split}
-\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 \np{rn\_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^{-6}$ is used for $\epsilon$ since larger 
-values may lead to numerically induced basin scale barotropic oscillations. 
-The precision is specified by setting \np{rn\_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, \np{nn\_max} (\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 a abort.nc NetCDF 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 
-$nn\_max$ 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. 
-\eqref{Eq_sor_algo} can be been rewritten as:
-\begin{equation} 
-\begin{split}
-R_{i,j}^n
-= &a_{i,j}^{N}  D_{i+1,j}^n +\,a_{i,j}^{E}  D_{i,j+1} ^n
- +\,a_{i,j}^{S}  D_{i-1,j} ^{n}+\,_{i,j}^{W} D_{i,j-1} ^{n} -  D_{i,j}^n - \tilde{b}_{i,j}      \\
-R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{S}\; R_{i,j-1}^n                                             \\
-R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{W}\; R_{i-1,j}^n
-\end{split}
-\end{equation}
-This technique slightly increases the number of iteration required to reach the convergence,
-but this is largely compensated by the gain obtained by the suppression of the recurrences.
-
-Another technique have been chosen, the so-called red-black SOR. It consist in solving successively 
-\eqref{Eq_sor_algo} for odd and even grid points. It also slightly reduced the convergence rate
-but allows the vectorisation. In addition, and this is the reason why it has been chosen, it is able to handle the north fold boundary condition used in ORCA configuration ($i.e.$ tri-polar global ocean mesh).
-
-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  (\np{nn\_solv}=1, \mdl{solpcg}) }
-\label{MISC_solpcg}
-
-\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*} 
-\begin{split}
-\textbf{x}^0 &= D_{i,j}^0   = 2 D_{i,j}^t - D_{i,j}^{t-1}       \quad, \text{the initial guess }     \\
-\textbf{r}^0 &= \textbf{d}^0 = \textbf{b} - \textbf{A x}^0       \\
-\gamma_0 &= \langle{ \textbf{r}^0 , \textbf{r}^0} \rangle
-\end{split}
-\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, \np{nn\_max}, 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 
-\np{rn\_eps} and \np{nn\_max} (\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{Madec_al_OM88}. 
-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.
-
-% ================================================================
-
-
-
-
-
-
-
-
-
-
-
-
+% ================================================================
+
+
+
+
+
