Index: trunk/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
 trunk/DOC/TexFiles/Biblio/Biblio.bib (revision 6140)
+++ trunk/DOC/TexFiles/Biblio/Biblio.bib (revision 6289)
@@ 154,7 +154,7 @@
author={V. Artale and D. Iudicone and R. Santoleri and V. Rupolo and S. Marullo
and F. {D'O}rtenzio},
 title={{Role of surface fluxes in ocean general circulation models using satellite
+ title={Role of surface fluxes in ocean general circulation models using satellite
sea surface temperature: Validation of and sensitivity to the forcing
 frequency of the Mediterranean thermohaline circulation}},
+ frequency of the Mediterranean thermohaline circulation},
journal=JGR,
year={2002},
@@ 504,4 +504,27 @@
}
+@article{Brankart_OM2013,
+author = "J.M. Brankart",
+title = "Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modelling ",
+journal = OM,
+year = "2013",
+volume = "66", pages = "6476",
+issn = "14635003",
+doi = "http://dx.doi.org/10.1016/j.ocemod.2013.02.004",
+url = "http://www.sciencedirect.com/science/article/pii/S1463500313000309",
+}
+
+@Article{Brankart_al_GMD2015,
+AUTHOR = {Brankart, J.M. and Candille, G. and Garnier, F. and Calone, C. and Melet, A. and Bouttier, P.A. and Brasseur, P. and Verron, J.},
+TITLE = {A generic approach to explicit simulation of uncertainty in the NEMO ocean model},
+JOURNAL = {Geoscientific Model Development},
+VOLUME = {8},
+YEAR = {2015},
+NUMBER = {5},
+PAGES = {12851297},
+URL = {http://www.geoscimodeldev.net/8/1285/2015/},
+DOI = {10.5194/gmd812852015}
+}
+
@ARTICLE{Brodeau_al_OM09,
author = {L. Brodeau and B. Barnier and A.M. Tr\'{e}guier and T. Penduff and S. Gulev},
@@ 777,6 +800,6 @@
@PHDTHESIS{Demange_PhD2014,
 author = {J. Farge},
 title = {Sch\'{e}́mas num\'{e}́riques d’advection et de propagation d’ondes de gravit\'{e}́
+ author = {J. Demange},
+ title = {Sch\'{e}́mas num\'{e}́riques d'advection et de propagation d’ondes de gravit\'{e}́
dans les mod\`{e}les de circulation oc\'{e}́anique.},
school = {Doctorat es Applied Mathematiques, Grenoble University, France},
@@ 1239,4 +1262,13 @@
pages = {10811098},
}
+@ARTICLE{Griffies_Hallberg_MWR00,
+ author = {S.M. Griffies and R.H. Hallberg},
+ title = {Biharmonic friction with a smagorinskylike viscosity for use in largescale eddypermitting ocean models},
+ journal = MWR,
+ year = {2000},
+ volume = {128},
+ pages = {2935–2946},
+ url = {http://dx.doi.org/10.1175/15200493(2000)128}
+}
@ARTICLE{Guilyardi_al_JC04,
@@ 1721,4 +1753,15 @@
volume = {9},
pages = {91112}
+}
+
+@article{Lemarie_OM2012,
+author = "F. Lemari\'{e} and L. Debreu and A.F. Shchepetkin and J.C. McWilliams",
+title = "On the stability and accuracy of the harmonic and biharmonic isoneutral mixing operators in ocean models ",
+journal = OM,
+year = "2012",
+volume = "52–53", pages = "935",
+issn = "14635003",
+doi = "http://dx.doi.org/10.1016/j.ocemod.2012.04.007",
+url = "http://www.sciencedirect.com/science/article/pii/S1463500312000674",
}
@@ 1876,15 +1919,17 @@
year = {2009},
volume = {30}, number = {23},
 pages = {8894},
+ pages = {8894},
doi = {10.1016/j.ocemod.2009.06.006},
 url = {http://dx.doi.org/}
}

@ARTICLE{Leclair_Madec_OM10s,
+ url = {http://dx.doi.org/10.1016/j.ocemod.2009.06.006}
+}
+
+@ARTICLE{Leclair_Madec_OM11,
author = {M. Leclair and G. Madec},
title = {$\tilde{z}$coordinate, an Arbitrary LagrangianEulerian coordinate separating high and low frequency},
journal = OM,
 year = {2010},
 pages = {submitted},
+ year = {2011},
+ volume = {37}, pages = {139152},
+ doi = {10.1016/j.ocemod.2011.02.001},
+ url = {http://dx.doi.org/10.1016/j.ocemod.2011.02.001}
}
@@ 1894,6 +1939,5 @@
journal = JCP,
year = {1992},
 volume = {103}
 pages = {1642}
+ volume = {103}, pages = {1642}
}
@@ 1905,6 +1949,5 @@
journal = JC,
year = {2003},
 volume = {16}, number = {20},
 pages = {33303343}
+ volume = {16}, number = {20}, pages = {33303343}
}
@@ 2184,17 +2227,4 @@
volume = {3},
pages = {120}
}


@article{Martin_Adcroft_OM10,
author = {T. Martin and A. Adcroft},
title = {Parameterizing the freshwater flux from land ice to ocean with interactive icebergs in a coupled climate model},
journal = OM,
year = {2010},
volume = {34}, number = {34},
pages = {111124},
issn = {14635003},
doi = {10.1016/j.ocemod.2010.05.001},
url = {http://dx.doi.org/10.1016/j.ocemod.2010.05.001}
}
@@ 2210,4 +2240,28 @@
doi = {10.1016/j.ocemod.2007.07.005},
url = {http://dx.doi.org/10.1016/j.ocemod.2007.07.005}
+}
+
+@Article{Marsh_GMD2015,
+AUTHOR = {R. Marsh and V. O. Ivchenko and N. Skliris and S. Alderson and G. R. Bigg and G. Madec and A. T. Blaker and Y. Aksenov and B. Sinha and A. C. Coward and J. Le Sommer and N. Merino and V. B. Zalesny},
+TITLE = {NEMOICB (v1.0): interactive icebergs in the NEMO ocean model globally configured at eddypermitting resolution},
+JOURNAL = {Geoscientific Model Development},
+VOLUME = {8},
+YEAR = {2015},
+NUMBER = {5},
+PAGES = {15471562},
+url = {HTTP://www.geoscimodeldev.net/8/1547/2015/},
+DOI = {10.5194/gmd815472015}
+}
+
+@article{Martin_Adcroft_OM10,
+author = {T. Martin and A. Adcroft},
+title = {Parameterizing the freshwater flux from land ice to ocean with interactive icebergs in a coupled climate model},
+journal = OM,
+year = {2010},
+volume = {34}, number = {34},
+pages = {111124},
+issn = {14635003},
+doi = {10.1016/j.ocemod.2010.05.001},
+url = {http://dx.doi.org/10.1016/j.ocemod.2010.05.001}
}
@@ 2580,4 +2634,16 @@
}
+@Article{Rodgers_2014,
+AUTHOR = {Rodgers, K. B. and Aumont, O. and Mikaloff Fletcher, S. E. and Plancherel, Y. and Bopp, L. and de Boyer Mont\'egut, C. and Iudicone, D. and Keeling, R. F. and Madec, G. and Wanninkhof, R.},
+TITLE = {Strong sensitivity of Southern Ocean carbon uptake and nutrient cycling to wind stirring},
+JOURNAL = {Biogeosciences},
+VOLUME = {11},
+YEAR = {2014},
+NUMBER = {15},
+PAGES = {40774098},
+URL = {HTTP://www.biogeosciences.net/11/4077/2014/},
+DOI = {10.5194/bg1140772014}
+}
+
@ARTICLE{Rodi_1987,
author = {W. Rodi},
@@ 2834,12 +2900,12 @@
@ARTICLE{Takaya_al_JGR10,
 author = {Y. Takaya and JR. Bidlot and A. C. M. Beljaars and
+ author = {Y. Takaya and JR. Bidlot and A. C. M. Beljaars and
P. A. E. M. Janssen},
 title = {Refinements to a prognostic scheme of sea surface skin temperature},
 journal = JGR,
 year = {2010},
+ title = {Refinements to a prognostic scheme of sea surface skin temperature},
+ journal = JGR,
+ year = {2010},
Volume = {115},
 Pages = {C06009},
 doi = {10.1029/2009JC005985},
+ Pages = {C06009},
+ doi = {10.1029/2009JC005985},
}
Index: trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex (revision 6289)
@@ 14,5 +14,5 @@
the earth climate system over a wide range of space and time scales.
Prognostic variables are the threedimensional velocity field, a nonlinear sea surface height,
the \textit{Conservative} Temperature and the \textit{Absolut} Salinity.
+the \textit{Conservative} Temperature and the \textit{Absolute} Salinity.
In the horizontal direction, the model uses a curvilinear orthogonal grid and in the vertical direction,
a full or partial step $z$coordinate, or $s$coordinate, or a mixture of the two.
@@ 31,5 +31,5 @@
interactions avec les autres composantes du syst\`{e}me climatique terrestre.
Les variables pronostiques sont le champ tridimensionnel de vitesse, une hauteur de la mer
lin\'{e}aire, la Ttemperature Conservative et la Salinit\'{e} Absolue.
+lin\'{e}aire, la Temp\'{e}erature Conservative et la Salinit\'{e} Absolue.
La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une
coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore
Index: trunk/DOC/TexFiles/Chapters/Annex_C.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Annex_C.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Annex_C.tex (revision 6289)
@@ 1103,5 +1103,5 @@
The discrete formulation of the horizontal diffusion of momentum ensures the
conservation of potential vorticity and the horizontal divergence, and the
dissipation of the square of these quantities (i.e. enstrophy and the
+dissipation of the square of these quantities ($i.e.$ enstrophy and the
variance of the horizontal divergence) as well as the dissipation of the
horizontal kinetic energy. In particular, when the eddy coefficients are
@@ 1127,26 +1127,17 @@
&\int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times
\Bigl[ \nabla_h \left( A^{\,lm}\;\chi \right)
  \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv = 0
\end{flalign*}
+  \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv \\
+%\end{flalign*}
%%%%%%%%%% recheck here.... (gm)
\begin{flalign*}
= \int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times
 \Bigl[ \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv &&& \\
\end{flalign*}
\begin{flalign*}
+%\begin{flalign*}
+=& \int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times
+ \Bigl[ \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv \\
+%\end{flalign*}
+%\begin{flalign*}
\equiv& \sum\limits_{i,j}
\left\{
 \delta_{i+1/2}
 \left[
 \frac {e_{2v}} {e_{1v}\,e_{3v}} \delta_i
 \left[ A_f^{\,lm} e_{3f} \zeta \right]
 \right]
 + \delta_{j+1/2}
 \left[
 \frac {e_{1u}} {e_{2u}\,e_{3u}} \delta_j
 \left[ A_f^{\,lm} e_{3f} \zeta \right]
 \right]
 \right\}
 && \\
+ \delta_{i+1/2} \left[ \frac {e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ A_f^{\,lm} e_{3f} \zeta \right] \right]
+ + \delta_{j+1/2} \left[ \frac {e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ A_f^{\,lm} e_{3f} \zeta \right] \right]
+ \right\} \\
%
\intertext{Using \eqref{DOM_di_adj}, it follows:}
@@ 1154,10 +1145,8 @@
\equiv& \sum\limits_{i,j,k}
\,\left\{
 \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i
 \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_i \left[ 1\right]
 + \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j
 \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_j \left[ 1\right]
+ \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_i \left[ 1\right]
+ + \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_j \left[ 1\right]
\right\} \quad \equiv 0
 && \\
+ \\
\end{flalign*}
@@ 1167,5 +1156,4 @@
\subsection{Dissipation of Horizontal Kinetic Energy}
\label{Apdx_C.3.2}

The lateral momentum diffusion term dissipates the horizontal kinetic energy:
@@ 1221,5 +1209,4 @@
\label{Apdx_C.3.3}

The lateral momentum diffusion term dissipates the enstrophy when the eddy
coefficients are horizontally uniform:
@@ 1228,7 +1215,7 @@
\left[ \nabla_h \left( A^{\,lm}\;\chi \right)
 \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \right]\;dv &&&\\
&= A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times
+&\quad = A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times
\left[ \nabla_h \times \left( \zeta \; \textbf{k} \right) \right]\;dv &&&\\
&\equiv A^{\,lm} \sum\limits_{i,j,k} \zeta \;e_{3f}
+&\quad \equiv A^{\,lm} \sum\limits_{i,j,k} \zeta \;e_{3f}
\left\{ \delta_{i+1/2} \left[ \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ e_{3f} \zeta \right] \right]
+ \delta_{j+1/2} \left[ \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta \right] \right] \right\} &&&\\
@@ 1236,8 +1223,7 @@
\intertext{Using \eqref{DOM_di_adj}, it follows:}
%
&\equiv  A^{\,lm} \sum\limits_{i,j,k}
+&\quad \equiv  A^{\,lm} \sum\limits_{i,j,k}
\left\{ \left( \frac{1} {e_{1v}\,e_{3v}} \delta_i \left[ e_{3f} \zeta \right] \right)^2 b_v
 + \left( \frac{1} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta \right] \right)^2 b_u \right\} &&&\\
& \leq \;0 &&&\\
+ + \left( \frac{1} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta \right] \right)^2 b_u \right\} \quad \leq \;0 &&&\\
\end{flalign*}
@@ 1250,16 +1236,16 @@
When the horizontal divergence of the horizontal diffusion of momentum
(discrete sense) is taken, the term associated with the vertical curl of the
vorticity is zero locally, due to (!!! II.1.8 !!!!!). The resulting term conserves the
$\chi$ and dissipates $\chi^2$ when the eddy coefficients are
horizontally uniform.
+vorticity is zero locally, due to \eqref{Eq_DOM_div_curl}.
+The resulting term conserves the $\chi$ and dissipates $\chi^2$
+when the eddy coefficients are horizontally uniform.
\begin{flalign*}
& \int\limits_D \nabla_h \cdot
\Bigl[ \nabla_h \left( A^{\,lm}\;\chi \right)
 \nabla_h \times \left( A^{\,lm}\;\zeta \;\textbf{k} \right) \Bigr] dv
= \int\limits_D \nabla_h \cdot \nabla_h \left( A^{\,lm}\;\chi \right) dv &&&\\
+= \int\limits_D \nabla_h \cdot \nabla_h \left( A^{\,lm}\;\chi \right) dv \\
%
&\equiv \sum\limits_{i,j,k}
\left\{ \delta_i \left[ A_u^{\,lm} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} \left[ \chi \right] \right]
 + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right] \right\} &&&\\
+ + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right] \right\} \\
%
\intertext{Using \eqref{DOM_di_adj}, it follows:}
@@ 1267,6 +1253,6 @@
&\equiv \sum\limits_{i,j,k}
 \left\{ \frac{e_{2u}\,e_{3u}} {e_{1u}} A_u^{\,lm} \delta_{i+1/2} \left[ \chi \right] \delta_{i+1/2} \left[ 1 \right]
 + \frac{e_{1v}\,e_{3v}} {e_{2v}} A_v^{\,lm} \delta_{j+1/2} \left[ \chi \right] \delta_{j+1/2} \left[ 1 \right] \right\}
 \qquad \equiv 0 &&& \\
+ + \frac{e_{1v}\,e_{3v}} {e_{2v}} A_v^{\,lm} \delta_{j+1/2} \left[ \chi \right] \delta_{j+1/2} \left[ 1 \right] \right\}
+ \quad \equiv 0 \\
\end{flalign*}
@@ 1281,5 +1267,5 @@
\left[ \nabla_h \left( A^{\,lm}\;\chi \right)
 \nabla_h \times \left( A^{\,lm}\;\zeta \;\textbf{k} \right) \right]\; dv
 = A^{\,lm} \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\; dv &&&\\
+ = A^{\,lm} \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\; dv \\
%
&\equiv A^{\,lm} \sum\limits_{i,j,k} \frac{1} {e_{1t}\,e_{2t}\,e_{3t}} \chi
@@ 1287,5 +1273,5 @@
\delta_i \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} \left[ \chi \right] \right]
+ \delta_j \left[ \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right]
 \right\} \; e_{1t}\,e_{2t}\,e_{3t} &&&\\
+ \right\} \; e_{1t}\,e_{2t}\,e_{3t} \\
%
\intertext{Using \eqref{DOM_di_adj}, it turns out to be:}
@@ 1293,7 +1279,6 @@
&\equiv  A^{\,lm} \sum\limits_{i,j,k}
\left\{ \left( \frac{1} {e_{1u}} \delta_{i+1/2} \left[ \chi \right] \right)^2 b_u
 + \left( \frac{1} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right)^2 b_v \right\} \; &&&\\
%
&\leq 0 &&&\\
+ + \left( \frac{1} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right)^2 b_v \right\}
+\quad \leq 0 \\
\end{flalign*}
@@ 1303,5 +1288,4 @@
\section{Conservation Properties on Vertical Momentum Physics}
\label{Apdx_C_4}

As for the lateral momentum physics, the continuous form of the vertical diffusion
@@ 1319,6 +1303,6 @@
\left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right)\; dv \quad &\leq 0 \\
\end{align*}
+
The first property is obvious. The second results from:

\begin{flalign*}
\int\limits_D
@@ 1359,4 +1343,5 @@
e_{1f}\,e_{2f}\,e_{3f} \; \equiv 0 && \\
\end{flalign*}
+
If the vertical diffusion coefficient is uniform over the whole domain, the
enstrophy is dissipated, $i.e.$
@@ 1366,6 +1351,6 @@
\left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right) \right)\; dv = 0 &&&\\
\end{flalign*}
+
This property is only satisfied in $z$coordinates:

\begin{flalign*}
\int\limits_D \zeta \, \textbf{k} \cdot \nabla \times
@@ 1477,10 +1462,8 @@
The numerical schemes used for tracer subgridscale physics are written such
that the heat and salt contents are conserved (equations in flux form, second
order centered finite differences). Since a flux form is used to compute the
temperature and salinity, the quadratic form of these quantities (i.e. their variance)
globally tends to diminish. As for the advection term, there is generally no strict
conservation of mass, even if in practice the mass is conserved to a very high
accuracy.
+that the heat and salt contents are conserved (equations in flux form).
+Since a flux form is used to compute the temperature and salinity,
+the quadratic form of these quantities ($i.e.$ their variance) globally tends to diminish.
+As for the advection term, there is conservation of mass only if the Equation Of Seawater is linear.
% 
Index: trunk/DOC/TexFiles/Chapters/Annex_D.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Annex_D.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Annex_D.tex (revision 6289)
@@ 120,5 +120,5 @@
\hline
public \par or \par module variable&
\textbf{m n} \par \textit{but not} \par \textbf{nn\_}&
+\textbf{m n} \par \textit{but not} \par \textbf{nn\_ np\_}&
\textbf{a b e f g h o q r} \par \textbf{t} \textit{to} \textbf{x} \par but not \par \textbf{fs rn\_}&
\textbf{l} \par \textit{but not} \par \textbf{lp ld} \par \textbf{ ll ln\_}&
@@ 156,5 +156,5 @@
\hline
parameter&
\textbf{jp}&
+\textbf{jp np\_}&
\textbf{pp}&
\textbf{lp}&
@@ 190,4 +190,8 @@
%
+N.B. Parameter here, in not only parameter in the \textsc{Fortran} acceptation, it is also used for code variables
+that are read in namelist and should never been modified during a simulation.
+It is the case, for example, for the size of a domain (jpi,jpj,jpk).
+
\newpage
% ================================================================
Index: trunk/DOC/TexFiles/Chapters/Annex_ISO.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Annex_ISO.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Annex_ISO.tex (revision 6289)
@@ 11,58 +11,47 @@
\namdisplay{namtra_ldf}
%
If the namelist variable \np{ln\_traldf\_grif} is set true (and
\key{ldfslp} is set), \NEMO updates both active and passive tracers
using the Griffies triad representation of isoneutral diffusion and
the eddyinduced advective skew (GM) fluxes. Otherwise (by default) the
filtered version of Cox's original scheme is employed
(\S\ref{LDF_slp}). In the present implementation of the Griffies
scheme, the advective skew fluxes are implemented even if
\key{traldf\_eiv} is not set.
+
+Two scheme are available to perform the isoneutral diffusion.
+If the namelist logical \np{ln\_traldf\_triad} is set true,
+\NEMO updates both active and passive tracers using the Griffies triad representation
+of isoneutral diffusion and the eddyinduced advective skew (GM) fluxes.
+If the namelist logical \np{ln\_traldf\_iso} is set true,
+the filtered version of Cox's original scheme (the Standard scheme) is employed (\S\ref{LDF_slp}).
+In the present implementation of the Griffies scheme,
+the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false.
Values of isoneutral diffusivity and GM coefficient are set as
described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd},
N=1,2,3 is set (the default), spatially constant isoneutral $A_l$ and
GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and
\np{rn\_aeiv\_0}. If 2Dvarying coefficients are set with
\key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal
scale factor according to \eqref{Eq_title} \footnote{Except in global ORCA
 $0.5^{\circ}$ runs with \key{traldf\_eiv}, where
 $A_l$ is set like $A_e$ but with a minimum vale of
 $100\;\mathrm{m}^2\;\mathrm{s}^{1}$}. In idealised setups with
\key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv}
is set in the global configurations with \key{traldf\_c2d}, a horizontally varying $A_e$ is
instead set from the HeldLarichev parameterisation\footnote{In this
 case, $A_e$ at low latitudes $\theta<20^{\circ}$ is further
 reduced by a factor $f/f_{20}$, where $f_{20}$ is the value of $f$
 at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored
unless it is zero.
+described in \S\ref{LDF_coef}. Note that when GM fluxes are used,
+the eddyadvective (GM) velocities are output for diagnostic purposes using xIOS,
+even though the eddy advection is accomplished by means of the skew fluxes.
+
The options specific to the Griffies scheme include:
\begin{description}[font=\normalfont]
\item[\np{ln\_traldf\_gdia}] Default value is false. See \S\ref{sec:triad:sfdiag}. If this is set true, timemean
 eddyadvective (GM) velocities are output for diagnostic purposes, even
 though the eddy advection is accomplished by means of the skew
 fluxes.
\item[\np{ln\_traldf\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then
+\item[\np{ln\_triad\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then
`isoneutral' mixing is accomplished within the surface mixedlayer
along slopes linearly decreasing with depth from the value immediately below
 the mixedlayer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}). This is the same
 treatment as used in the default implementation
 \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}. Where
 \np{ln\_traldf\_iso} is set true, the vertical skew flux is further
 reduced to ensure no vertical buoyancy flux, giving an almost pure
+ the mixedlayer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}).
+ This is the same treatment as used in the default implementation \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}.
+ Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced
+ to ensure no vertical buoyancy flux, giving an almost pure
horizontal diffusive tracer flux within the mixed layer. This is similar to
the tapering suggested by \citet{Gerdes1991}. See \S\ref{sec:triad:Gerdestaper}
\item[\np{ln\_traldf\_botmix}] See \S\ref{sec:triad:iso_bdry}. If this
 is set false (the default) then the lateral diffusive fluxes
 associated with triads partly masked by topography are neglected. If
 it is set true, however, then these lateral diffusive fluxes are
 applied, giving smoother bottom tracer fields at the cost of
 introducing diapycnal mixing.
+\item[\np{ln\_botmix\_triad}] See \S\ref{sec:triad:iso_bdry}.
+ If this is set false (the default) then the lateral diffusive fluxes
+ associated with triads partly masked by topography are neglected.
+ If it is set true, however, then these lateral diffusive fluxes are applied,
+ giving smoother bottom tracer fields at the cost of introducing diapycnal mixing.
+\item[\np{rn\_sw\_triad}] blah blah to be added....
+\end{description}
+The options shared with the Standard scheme include:
+\begin{description}[font=\normalfont]
+\item[\np{ln\_traldf\_msc}] blah blah to be added
+\item[\np{rn\_slpmax}] blah blah to be added
\end{description}
\section{Triad formulation of isoneutral diffusion}
\label{sec:triad:iso}
We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO
framework, using scale factors rather than gridsizes.
+We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98},
+but formulated within the \NEMO framework, using scale factors rather than gridsizes.
\subsection{The isoneutral diffusion operator}
@@ 84,13 +73,13 @@
\mbox{with}\quad \;\;\Re =
\begin{pmatrix}
 1&0&r_1\mystrut \\
 0&1&r_2\mystrut \\
 r_1&r_2&r_1 ^2+r_2 ^2\mystrut
+ 1 & 0 & r_1 \mystrut \\
+ 0 & 1 & r_2 \mystrut \\
+ r_1 & r_2 & r_1 ^2+r_2 ^2 \mystrut
\end{pmatrix}
\quad \text{and} \quad\grad T=
\begin{pmatrix}
 \frac{1}{e_1}\pd[T]{i}\mystrut \\
 \frac{1}{e_2}\pd[T]{j}\mystrut \\
 \frac{1}{e_3}\pd[T]{k}\mystrut
+ \frac{1}{e_1} \pd[T]{i} \mystrut \\
+ \frac{1}{e_2} \pd[T]{j} \mystrut \\
+ \frac{1}{e_3} \pd[T]{k} \mystrut
\end{pmatrix}.
\end{equation}
@@ 101,5 +90,5 @@
% {r_1 } \hfill & {r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
% \end{array} }} \right)
 Here \eqref{Eq_PE_iso_slopes}
+ Here \eqref{Eq_PE_iso_slopes}
\begin{align*}
r_1 &=\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i}
@@ 200,5 +189,5 @@
% the mean vertical gradient at the $u$point,
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[h] \begin{center}
+\begin{figure}[tb] \begin{center}
\includegraphics[width=1.05\textwidth]{./TexFiles/Figures/Fig_GRIFF_triad_fluxes}
\caption{ \label{fig:triad:ISO_triad}
@@ 256,17 +245,13 @@
\
\frac
 {\left(\alpha / \beta \right)_i^k \ \delta_{i + i_p}[T^k]  \delta_{i + i_p}[S^k] }
 {\left(\alpha / \beta \right)_i^k \ \delta_{k+k_p}[T^i ]  \delta_{k+k_p}[S^i ] }.
\end{equation}
In calculating the slopes of the local neutral
surfaces, the expansion coefficients $\alpha$ and $\beta$ are
evaluated at the anchor points of the triad \footnote{Note that in \eqref{eq:triad:R} we use the ratio $\alpha / \beta$
instead of multiplying the temperature derivative by $\alpha$ and the
salinity derivative by $\beta$. This is more efficient as the ratio
$\alpha / \beta$ can to be evaluated directly}, while the metrics are
calculated at the $u$ and $w$points on the arms.
+ { \alpha_i^k \ \delta_{i+i_p}[T^k]  \beta_i^k \ \delta_{i+i_p}[S^k] }
+ { \alpha_i^k \ \delta_{k+k_p}[T^i]  \beta_i^k \ \delta_{k+k_p}[S^i] }.
+\end{equation}
+In calculating the slopes of the local neutral surfaces,
+the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad,
+while the metrics are calculated at the $u$ and $w$points on the arms.
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[h] \begin{center}
+\begin{figure}[tb] \begin{center}
\includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_GRIFF_qcells}
\caption{ \label{fig:triad:qcells}
@@ 277,20 +262,17 @@
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter
cell that is the intersection of the $i,k$ $T$cell, the $i+i_p,k$
$u$cell and the $i,k+k_p$ $w$cell. Expressing the slopes $s_i$ and
$s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation, we have
e.g.\ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. Each triad slope $_i^k
\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to calculate the
lateral flux along its $u$arm, at $(i+i_p,k)$, and then again as an
$s'$ to calculate the vertical flux along its $w$arm at
$(i,k+k_p)$. Each vertical area $a_i$ used to calculate the lateral
flux and horizontal area $a'_i$ used to calculate the vertical flux
can also be identified as the area across the $u$ and $w$arms of a
unique triad, and we notate these areas, similarly to the triad
slopes, as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$,
$_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq:triad:i13}
$a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq:triad:i31}
$a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
+Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter
+cell that is the intersection of the $i,k$ $T$cell, the $i+i_p,k$ $u$cell and the $i,k+k_p$ $w$cell.
+Expressing the slopes $s_i$ and $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation,
+we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$.
+Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$)
+to calculate the lateral flux along its $u$arm, at $(i+i_p,k)$,
+and then again as an $s'$ to calculate the vertical flux along its $w$arm at $(i,k+k_p)$.
+Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used
+to calculate the vertical flux can also be identified as the area across the $u$ and $w$arms
+of a unique triad, and we notate these areas, similarly to the triad slopes,
+as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$,
+where $e.g.$ in \eqref{eq:triad:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$,
+and in \eqref{eq:triad:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
\subsection{The full triad fluxes}
@@ 667,11 +649,11 @@
or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$point is
masked. The associated lateral fluxes (greyblack dashed line) are
masked if \np{ln\_botmix\_grif}=false, but left unmasked,
giving bottom mixing, if \np{ln\_botmix\_grif}=true.

The default option \np{ln\_botmix\_grif}=false is suitable when the
+masked if \np{ln\_botmix\_triad}=false, but left unmasked,
+giving bottom mixing, if \np{ln\_botmix\_triad}=true.
+
+The default option \np{ln\_botmix\_triad}=false is suitable when the
bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}=1),
or for simple idealized problems. For setups with topography without
bbl mixing, \np{ln\_botmix\_grif}=true may be necessary.
+bbl mixing, \np{ln\_botmix\_triad}=true may be necessary.
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[h] \begin{center}
@@ 690,6 +672,6 @@
or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$point
is masked. The associated lateral fluxes (greyblack dashed
 line) are masked if \np{botmix\_grif}=.false., but left
 unmasked, giving bottom mixing, if \np{botmix\_grif}=.true.}
+ line) are masked if \np{botmix\_triad}=.false., but left
+ unmasked, giving bottom mixing, if \np{botmix\_triad}=.true.}
\end{center} \end{figure}
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ 931,5 +913,5 @@
it to the Eulerian velocity prior to computing the tracer
advection. This is implemented if \key{traldf\_eiv} is set in the
default implementation, where \np{ln\_traldf\_grif} is set
+default implementation, where \np{ln\_traldf\_triad} is set
false. This allows us to take advantage of all the advection schemes
offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$
@@ 938,5 +920,5 @@
paramount importance.
However, when \np{ln\_traldf\_grif} is set true, \NEMO instead
+However, when \np{ln\_traldf\_triad} is set true, \NEMO instead
implements eddy induced advection according to the socalled skew form
\citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes
@@ 1137,5 +1119,5 @@
and $\triadt{i+1}{k}{R}{1/2}{1/2}$ are masked when either of the
$i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$
$u$point is masked. The namelist parameter \np{ln\_botmix\_grif} has
+$u$point is masked. The namelist parameter \np{ln\_botmix\_triad} has
no effect on the eddyinduced skewfluxes.
Index: trunk/DOC/TexFiles/Chapters/Chap_DIA.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_DIA.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_DIA.tex (revision 6289)
@@ 12,5 +12,5 @@
% Old Model Output
% ================================================================
\section{Old Model Output (default or \key{dimgout})}
+\section{Old Model Output (default)}
\label{DIA_io_old}
@@ 33,9 +33,18 @@
"\textit{grep i numout}" in the source code directory.
By default, diagnostic output files are written in NetCDF format but an IEEE binary output format, called DIMG, can be chosen by defining \key{dimgout}.

Since version 3.2, when defining \key{iomput}, an I/O server has been added which provides more flexibility in the choice of the fields to be written as well as how the writing work is distributed over the processors in massively parallel computing. The complete description of the use of this I/O server is presented in the next section.

By default, if neither \key{iomput} nor \key{dimgout} are defined, NEMO produces NetCDF with the old IOIPSL library which has been kept for compatibility and its easy installation. However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2, many diagnostic options have been added presuming the use of \key{iomput}. The usefulness of the default IOIPSLbased option is expected to reduce with each new release. If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of nn\_write timesteps (namelist parameter).
+By default, diagnostic output files are written in NetCDF format.
+Since version 3.2, when defining \key{iomput}, an I/O server has been added
+which provides more flexibility in the choice of the fields to be written
+as well as how the writing work is distributed over the processors in massively parallel computing.
+A complete description of the use of this I/O server is presented in the next section.
+
+By default, \key{iomput} is not defined, NEMO produces NetCDF with the old IOIPSL library
+which has been kept for compatibility and its easy installation.
+However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2,
+many diagnostic options have been added presuming the use of \key{iomput}.
+The usefulness of the default IOIPSLbased option is expected to reduce with each new release.
+If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module
+and contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period
+of nn\_write timesteps (namelist parameter).
%\gmcomment{ % start of gmcomment
@@ 48,13 +57,19 @@
Since version 3.2, iomput is the NEMO output interface of choice. It has been designed to be simple to use, flexible and efficient. The two main purposes of iomput are:
+Since version 3.2, iomput is the NEMO output interface of choice.
+It has been designed to be simple to use, flexible and efficient.
+The two main purposes of iomput are:
\begin{enumerate}
\item The complete and flexible control of the output files through external XML files adapted by the user from standard templates.
\item To achieve high performance and scalable output through the optional distribution of all diagnostic output related tasks to dedicated processes.
+\item The complete and flexible control of the output files through external XML files
+adapted by the user from standard templates.
+\item To achieve high performance and scalable output through the optional distribution
+of all diagnostic output related tasks to dedicated processes.
\end{enumerate}
The first functionality allows the user to specify, without code changes or recompilation, aspects of the diagnostic output stream, such as:
+The first functionality allows the user to specify, without code changes or recompilation,
+aspects of the diagnostic output stream, such as:
\begin{itemize}
\item The choice of output frequencies that can be different for each file (including real months and years).
\item The choice of file contents; includes complete flexibility over which data are written in which files (the same data can be written in different files).
+\item The choice of file contents; includes complete flexibility over which data are written in which files
+(the same data can be written in different files).
\item The possibility to split output files at a chosen frequency.
\item The possibility to extract a vertical or an horizontal subdomain.
@@ 62,5 +77,7 @@
\item Control over metadata via a large XML "database" of possible output fields.
\end{itemize}
In addition, iomput allows the user to add the output of any new variable (scalar, 2D or 3D) in the code in a very easy way. All details of iomput functionalities are listed in the following subsections. Examples of the XML files that control the outputs can be found in:
+In addition, iomput allows the user to add in the code the output of any new variable (scalar, 2D or 3D)
+in a very easy way. All details of iomput functionalities are listed in the following subsections.
+Examples of the XML files that control the outputs can be found in:
\begin{alltt}
\begin{verbatim}
@@ 72,9 +89,24 @@
\end{alltt}
The second functionality targets output performance when running in parallel (\key{mpp\_mpi}). Iomput provides the possibility to specify N dedicated I/O processes (in addition to the NEMO processes) to collect and write the outputs. With an appropriate choice of N by the user, the bottleneck associated with the writing of the output files can be greatly reduced.

In version 3.6, the iom\_put interface depends on an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios1.0}{XIOS1.0} (use of revision 618 or higher is required). This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to create a single output file and therefore to bypass the rebuilding phase. Note that writing in parallel into the same NetCDF files requires that your NetCDF4 library is linked to an HDF5 library that has been correctly compiled (i.e. with the configure option $$enableparallel). Note that the files created by iomput through XIOS are incompatible with NetCDF3. All postprocesssing and visualization tools must therefore be compatible with NetCDF4 and not only NetCDF3.

Even if not using the parallel I/O functionality of NetCDF4, using N dedicated I/O servers, where N is typically much less than the number of NEMO processors, will reduce the number of output files created. This can greatly reduce the postprocessing burden usually associated with using large numbers of NEMO processors. Note that for smaller configurations, the rebuilding phase can be avoided, even without a parallelenabled NetCDF4 library, simply by employing only one dedicated I/O server.
+The second functionality targets output performance when running in parallel (\key{mpp\_mpi}).
+Iomput provides the possibility to specify N dedicated I/O processes (in addition to the NEMO processes)
+to collect and write the outputs. With an appropriate choice of N by the user,
+the bottleneck associated with the writing of the output files can be greatly reduced.
+
+In version 3.6, the iom\_put interface depends on an external code called
+\href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios1.0}{XIOS1.0}
+(use of revision 618 or higher is required). This new IO server can take advantage of
+the parallel I/O functionality of NetCDF4 to create a single output file and therefore
+to bypass the rebuilding phase. Note that writing in parallel into the same NetCDF files
+requires that your NetCDF4 library is linked to an HDF5 library that has been correctly compiled
+($i.e.$ with the configure option $$enableparallel).
+Note that the files created by iomput through XIOS are incompatible with NetCDF3.
+All postprocesssing and visualization tools must therefore be compatible with NetCDF4 and not only NetCDF3.
+
+Even if not using the parallel I/O functionality of NetCDF4, using N dedicated I/O servers,
+where N is typically much less than the number of NEMO processors, will reduce the number of output files created.
+This can greatly reduce the postprocessing burden usually associated with using large numbers of NEMO processors.
+Note that for smaller configurations, the rebuilding phase can be avoided,
+even without a parallelenabled NetCDF4 library, simply by employing only one dedicated I/O server.
\subsection{XIOS: the IO\_SERVER}
@@ 82,13 +114,39 @@
\subsubsection{Attached or detached mode?}
Iomput is based on \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS}, the io\_server developed by Yann Meurdesoif from IPSL. The behaviour of the io subsystem is controlled by settings in the external XML files listed above. Key settings in the iodef.xml file are {\tt using\_server} and the {\tt type} tag associated with each defined file. The {\tt using\_server} setting determines whether or not the server will be used in ''attached mode'' (as a library) [{\tt false}] or in ''detached mode'' (as an external executable on N additional, dedicated cpus) [{\tt true}]. The ''attached mode'' is simpler to use but much less efficient for massively parallel applications. The type of each file can be either ''multiple\_file'' or ''one\_file''.

In attached mode and if the type of file is ''multiple\_file'', then each NEMO process will also act as an IO server and produce its own set of output files. Superficially, this emulates the standard behaviour in previous versions, However, the subdomain written out by each process does not correspond to the {\tt jpi x jpj x jpk} domain actually computed by the process (although it may if {\tt jpni=1}). Instead each process will have collected and written out a number of complete longitudinal strips. If the ''one\_file'' option is chosen then all processes will collect their longitudinal strips and write (in parallel) to a single output file.

In detached mode and if the type of file is ''multiple\_file'', then each standalone XIOS process will collect data for a range of complete longitudinal strips and write to its own set of output files. If the ''one\_file'' option is chosen then all XIOS processes will collect their longitudinal strips and write (in parallel) to a single output file. Note running in detached mode requires launching a Multiple Process Multiple Data (MPMD) parallel job. The following subsection provides a typical example but the syntax will vary in different MPP environments.
+Iomput is based on \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS},
+the io\_server developed by Yann Meurdesoif from IPSL.
+The behaviour of the I/O subsystem is controlled by settings in the external XML files listed above.
+Key settings in the iodef.xml file are {\tt using\_server} and the {\tt type} tag associated with each defined file.
+The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode} (as a library)
+[{\tt false}] or in \textit{detached mode} (as an external executable on N additional, dedicated cpus) [{\tt true}].
+The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications.
+The type of each file can be either ''multiple\_file'' or ''one\_file''.
+
+In \textit{attached mode} and if the type of file is ''multiple\_file'',
+then each NEMO process will also act as an IO server and produce its own set of output files.
+Superficially, this emulates the standard behaviour in previous versions.
+However, the subdomain written out by each process does not correspond to the {\tt jpi x jpj x jpk}
+domain actually computed by the process (although it may if {\tt jpni=1}).
+Instead each process will have collected and written out a number of complete longitudinal strips.
+If the ''one\_file'' option is chosen then all processes will collect their longitudinal strips
+and write (in parallel) to a single output file.
+
+In \textit{detached mode} and if the type of file is ''multiple\_file'',
+then each standalone XIOS process will collect data for a range of complete longitudinal strips
+and write to its own set of output files. If the ''one\_file'' option is chosen then
+all XIOS processes will collect their longitudinal strips and write (in parallel) to a single output file.
+Note running in detached mode requires launching a Multiple Process Multiple Data (MPMD) parallel job.
+The following subsection provides a typical example but the syntax will vary in different MPP environments.
\subsubsection{Number of cpu used by XIOS in detached mode}
The number of cores used by the XIOS is specified when launching the model. The number of cores dedicated to XIOS should be from ~1/10 to ~1/50 of the number or cores dedicated to NEMO. Some manufacturers suggest using O($\sqrt{N}$) dedicated IO processors for N processors but this is a general recommendation and not specific to NEMO. It is difficult to provide precise recommendations because the optimal choice will depend on the particular hardware properties of the target system (parallel filesystem performance, available memory, memory bandwidth etc.) and the volume and frequency of data to be created. Here is an example of 2 cpus for the io\_server and 62 cpu for nemo using mpirun:
+The number of cores used by the XIOS is specified when launching the model.
+The number of cores dedicated to XIOS should be from ~1/10 to ~1/50 of the number or cores dedicated to NEMO.
+Some manufacturers suggest using O($\sqrt{N}$) dedicated IO processors for N processors
+but this is a general recommendation and not specific to NEMO.
+It is difficult to provide precise recommendations because the optimal choice will depend on
+the particular hardware properties of the target system (parallel filesystem performance, available memory,
+memory bandwidth etc.) and the volume and frequency of data to be created.
+Here is an example of 2 cpus for the io\_server and 62 cpu for nemo using mpirun:
\texttt{ mpirun np 62 ./nemo.exe : np 2 ./xios\_server.exe }
@@ 96,5 +154,6 @@
\subsubsection{Control of XIOS: the XIOS context in iodef.xml}
As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml. See the XML basics section below for more details on XML syntax and rules.
+As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml.
+See the XML basics section below for more details on XML syntax and rules.
\begin{tabular}{p{4cm}p{6.0cm}p{2.0cm}}
@@ 106,5 +165,6 @@
\hline
buffer\_size &
 buffer size used by XIOS to send data from NEMO to XIOS. Larger is more efficient. Note that needed/used buffer sizes are summarized at the end of the job &
+ buffer size used by XIOS to send data from NEMO to XIOS. Larger is more efficient.
+ Note that needed/used buffer sizes are summarized at the end of the job &
25000000 \\
\hline
@@ 136,10 +196,16 @@
\subsubsection{Installation}
As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with NEMO. See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance. NEMO will need to link to the compiled XIOS library. The
\href{http://www.nemoocean.eu/UsingNEMO/UserGuides/Basics/XIOSIOserverinstallationanduse}{XIOS with NEMO} guide provides an example illustration of how this can be achieved.
+As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with NEMO.
+See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance.
+NEMO will need to link to the compiled XIOS library.
+The \href{http://www.nemoocean.eu/UsingNEMO/UserGuides/Basics/XIOSIOserverinstallationanduse}{XIOS with NEMO}
+guide provides an example illustration of how this can be achieved.
\subsubsection{Add your own outputs}
It is very easy to add your own outputs with iomput. Many standard fields and diagnostics are already prepared (i.e., steps 1 to 3 below have been done) and simply need to be activated by including the required output in a file definition in iodef.xml (step 4). To add new output variables, all 4 of the following steps must be taken.
+It is very easy to add your own outputs with iomput.
+Many standard fields and diagnostics are already prepared ($i.e.$, steps 1 to 3 below have been done)
+and simply need to be activated by including the required output in a file definition in iodef.xml (step 4).
+To add new output variables, all 4 of the following steps must be taken.
\begin{description}
\item[1.] in NEMO code, add a \\
@@ 151,5 +217,7 @@
to the list of used modules in the upper part of your module.
\item[3.] in the field\_def.xml file, add the definition of your variable using the same identifier you used in the f90 code (see subsequent sections for a details of the XML syntax and rules). For example:
+\item[3.] in the field\_def.xml file, add the definition of your variable using the same identifier
+you used in the f90 code (see subsequent sections for a details of the XML syntax and rules).
+For example:
\vspace{20pt}
\begin{alltt} {{\scriptsize
@@ 165,5 +233,9 @@
\end{verbatim}
}}\end{alltt}
Note your definition must be added to the field\_group whose reference grid is consistent with the size of the array passed to iomput. The grid\_ref attribute refers to definitions set in iodef.xml which, in turn, reference grids and axes either defined in the code (iom\_set\_domain\_attr and iom\_set\_axis\_attr in iom.F90) or defined in the domain\_def.xml file. E.g.:
+Note your definition must be added to the field\_group whose reference grid is consistent
+with the size of the array passed to iomput.
+The grid\_ref attribute refers to definitions set in iodef.xml which, in turn, reference grids
+and axes either defined in the code (iom\_set\_domain\_attr and iom\_set\_axis\_attr in iom.F90)
+or defined in the domain\_def.xml file. $e.g.$:
\vspace{20pt}
\begin{alltt} {{\scriptsize
@@ 173,5 +245,6 @@
}}\end{alltt}
Note, if your array is computed within the surface module each nn\_fsbc time\_step,
add the field definition within the field\_group defined with the id ''SBC'': $<$field\_group id=''SBC''...$>$ which has been defined with the correct frequency of operations (iom\_set\_field\_attr in iom.F90)
+add the field definition within the field\_group defined with the id ''SBC'': $<$field\_group id=''SBC''...$>$
+which has been defined with the correct frequency of operations (iom\_set\_field\_attr in iom.F90)
\item[4.] add your field in one of the output files defined in iodef.xml (again see subsequent sections for syntax and rules) \\
@@ 201,5 +274,6 @@
\subsubsection{Structure of the xml file used in NEMO}
The XML file used in XIOS is structured by 7 families of tags: context, axis, domain, grid, field, file and variable. Each tag family has hierarchy of three flavors (except for context):
+The XML file used in XIOS is structured by 7 families of tags: context, axis, domain, grid, field, file and variable.
+Each tag family has hierarchy of three flavors (except for context):
\\
\begin{tabular}{p{3.0cm}p{4.5cm}p{4.5cm}}
@@ 225,7 +299,14 @@
\\
Each element may have several attributes. Some attributes are mandatory, other are optional but have a default value and other are are completely optional. Id is a special attribute used to identify an element or a group of elements. It must be unique for a kind of element. It is optional, but no reference to the corresponding element can be done if it is not defined.

The XML file is split into context tags that are used to isolate IO definition from different codes or different parts of a code. No interference is possible between 2 different contexts. Each context has its own calendar and an associated timestep. In NEMO, we used the following contexts (that can be defined in any order):\\
+Each element may have several attributes.
+Some attributes are mandatory, other are optional but have a default value and other are are completely optional.
+Id is a special attribute used to identify an element or a group of elements.
+It must be unique for a kind of element.
+It is optional, but no reference to the corresponding element can be done if it is not defined.
+
+The XML file is split into context tags that are used to isolate IO definition from different codes
+or different parts of a code. No interference is possible between 2 different contexts.
+Each context has its own calendar and an associated timestep.
+In \NEMO, we used the following contexts (that can be defined in any order):\\
\\
\begin{tabular}{p{3.0cm}p{4.5cm}p{4.5cm}}
@@ 271,5 +352,6 @@
\\
\noindent Each context tag related to NEMO (mother or child grids) is divided into 5 parts (that can be defined in any order):\\
+\noindent Each context tag related to NEMO (mother or child grids) is divided into 5 parts
+(that can be defined in any order):\\
\\
\begin{tabular}{p{3.0cm}p{4.5cm}p{4.5cm}}
@@ 305,5 +387,6 @@
\subsubsection{Nesting XML files}
The XML file can be split in different parts to improve its readability and facilitate its use. The inclusion of XML files into the main XML file can be done through the attribute src: \\
+The XML file can be split in different parts to improve its readability and facilitate its use.
+The inclusion of XML files into the main XML file can be done through the attribute src: \\
{\scriptsize \verb? ?}\\
@@ 323,5 +406,8 @@
\subsubsection{Use of inheritance}
XML extensively uses the concept of inheritance. XML has a tree based structure with a parentchild oriented relation: all children inherit attributes from parent, but an attribute defined in a child replace the inherited attribute value. Note that the special attribute ''id'' is never inherited. \\
+XML extensively uses the concept of inheritance.
+XML has a tree based structure with a parentchild oriented relation:
+all children inherit attributes from parent, but an attribute defined in a child replace the inherited attribute value.
+Note that the special attribute ''id'' is never inherited. \\
\\
example 1: Direct inheritance.
@@ 362,5 +448,8 @@
\subsubsection{Use of Groups}
Groups can be used for 2 purposes. Firstly, the group can be used to define common attributes to be shared by the elements of the group through inheritance. In the following example, we define a group of field that will share a common grid ''grid\_T\_2D''. Note that for the field ''toce'', we overwrite the grid definition inherited from the group by ''grid\_T\_3D''.
+Groups can be used for 2 purposes.
+Firstly, the group can be used to define common attributes to be shared by the elements of the group through inheritance.
+In the following example, we define a group of field that will share a common grid ''grid\_T\_2D''.
+Note that for the field ''toce'', we overwrite the grid definition inherited from the group by ''grid\_T\_3D''.
\vspace{20pt}
\begin{alltt} {{\scriptsize
@@ 375,5 +464,7 @@
}}\end{alltt}
Secondly, the group can be used to replace a list of elements. Several examples of groups of fields are proposed at the end of the file {\tt CONFIG/SHARED/field\_def.xml}. For example, a short list of the usual variables related to the U grid:
+Secondly, the group can be used to replace a list of elements.
+Several examples of groups of fields are proposed at the end of the file {\tt CONFIG/SHARED/field\_def.xml}.
+For example, a short list of the usual variables related to the U grid:
\vspace{20pt}
\begin{alltt} {{\scriptsize
@@ 399,8 +490,12 @@
\subsection{Detailed functionalities }
The file {\tt NEMOGCM/CONFIG/ORCA2\_LIM/iodef\_demo.xml} provides several examples of the use of the new functionalities offered by the XML interface of XIOS.
+The file {\tt NEMOGCM/CONFIG/ORCA2\_LIM/iodef\_demo.xml} provides several examples of the use
+of the new functionalities offered by the XML interface of XIOS.
\subsubsection{Define horizontal subdomains}
Horizontal subdomains are defined through the attributs zoom\_ibegin, zoom\_jbegin, zoom\_ni, zoom\_nj of the tag family domain. It must therefore be done in the domain part of the XML file. For example, in {\tt CONFIG/SHARED/domain\_def.xml}, we provide the following example of a definition of a 5 by 5 box with the bottom left corner at point (10,10).
+Horizontal subdomains are defined through the attributs zoom\_ibegin, zoom\_jbegin, zoom\_ni, zoom\_nj of the tag family domain.
+It must therefore be done in the domain part of the XML file.
+For example, in {\tt CONFIG/SHARED/domain\_def.xml}, we provide the following example of a definition
+of a 5 by 5 box with the bottom left corner at point (10,10).
\vspace{20pt}
\begin{alltt} {{\scriptsize
@@ 419,5 +514,10 @@
\end{verbatim}
}}\end{alltt}
Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain. The Equatorial section, the TAO, RAMA and PIRATA moorings are alredy registered in the code and can therefore be outputted without taking care of their (i,j) position in the grid. These predefined domains can be activated by the use of specific domain\_ref: ''EqT'', ''EqU'' or ''EqW'' for the equatorial sections and the mooring position for TAO, RAMA and PIRATA followed by ''T'' (for example: ''8s137eT'', ''1.5s80.5eT'' ...)
+Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain.
+The Equatorial section, the TAO, RAMA and PIRATA moorings are alredy registered in the code
+and can therefore be outputted without taking care of their (i,j) position in the grid.
+These predefined domains can be activated by the use of specific domain\_ref: ''EqT'', ''EqU'' or ''EqW''
+for the equatorial sections and the mooring position for TAO, RAMA and PIRATA followed
+by ''T'' (for example: ''8s137eT'', ''1.5s80.5eT'' ...)
\vspace{20pt}
\begin{alltt} {{\scriptsize
@@ 1116,6 +1216,5 @@
% 
\section[Tracer/Dynamics Trends (TRD)]
 {Tracer/Dynamics Trends (\key{trdtra}, \key{trddyn}, \\
 \key{trddvor}, \key{trdmld})}
+ {Tracer/Dynamics Trends (\ngn{namtrd})}
\label{DIA_trd}
@@ 1124,40 +1223,33 @@
%
When \key{trddyn} and/or \key{trddyn} CPP variables are defined, each
trend of the dynamics and/or temperature and salinity time evolution equations
is stored in threedimensional arrays just after their computation ($i.e.$ at the end
of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). Options are defined by
\ngn{namtrd} namelist variables. These trends are then
used in \mdl{trdmod} (see TRD directory) every \textit{nn\_trd } timesteps.

What is done depends on the CPP keys defined:
+Each trend of the dynamics and/or temperature and salinity time evolution equations
+can be send to \mdl{trddyn} and/or \mdl{trdtra} modules (see TRD directory) just after their computation
+($i.e.$ at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines).
+This capability is controlled by options offered in \ngn{namtrd} namelist.
+Note that the output are done with xIOS, and therefore the \key{IOM} is required.
+
+What is done depends on the \ngn{namtrd} logical set to \textit{true}:
\begin{description}
\item[\key{trddyn}, \key{trdtra}] : a check of the basin averaged properties of the momentum
and/or tracer equations is performed ;
\item[\key{trdvor}] : a vertical summation of the moment tendencies is performed,
then the curl is computed to obtain the barotropic vorticity tendencies which are output ;
\item[\key{trdmld}] : output of the tracer tendencies averaged vertically
either over the mixed layer (\np{nn\_ctls}=0),
or over a fixed number of model levels (\np{nn\_ctls}$>$1 provides the number of level),
or over a spatially varying but temporally fixed number of levels (typically the base
of the winter mixed layer) read in \ifile{ctlsurf\_idx} (\np{nn\_ctls}=1) ;
+\item[\np{ln\_glo\_trd}] : at each \np{nn\_trd} timestep a check of the basin averaged properties
+of the momentum and tracer equations is performed. This also includes a check of $T^2$, $S^2$,
+$\tfrac{1}{2} (u^2+v2)$, and potential energy time evolution equations properties ;
+\item[\np{ln\_dyn\_trd}] : each 3D trend of the evolution of the two momentum components is output ;
+\item[\np{ln\_dyn\_mxl}] : each 3D trend of the evolution of the two momentum components averaged
+ over the mixed layer is output ;
+\item[\np{ln\_vor\_trd}] : a vertical summation of the moment tendencies is performed,
+ then the curl is computed to obtain the barotropic vorticity tendencies which are output ;
+\item[\np{ln\_KE\_trd}] : each 3D trend of the Kinetic Energy equation is output ;
+\item[\np{ln\_tra\_trd}] : each 3D trend of the evolution of temperature and salinity is output ;
+\item[\np{ln\_tra\_mxl}] : each 2D trend of the evolution of temperature and salinity averaged
+ over the mixed layer is output ;
\end{description}

The units in the output file can be changed using the \np{nn\_ucf} namelist parameter.
For example, in case of salinity tendency the units are given by PSU/s/\np{nn\_ucf}.
Setting \np{nn\_ucf}=86400 ($i.e.$ the number of second in a day) provides the tendencies in PSU/d.

When \key{trdmld} is defined, two time averaging procedure are proposed.
Setting \np{ln\_trdmld\_instant} to \textit{true}, a simple time averaging is performed,
so that the resulting tendency is the contribution to the change of a quantity between
the two instantaneous values taken at the extremities of the time averaging period.
Setting \np{ln\_trdmld\_instant} to \textit{false}, a double time averaging is performed,
so that the resulting tendency is the contribution to the change of a quantity between
two \textit{time mean} values. The later option requires the use of an extra file, \ifile{restart\_mld}
(\np{ln\_trdmld\_restart}=true), to restart a run.

Note that the mixed layer tendency diagnostic can also be used on biogeochemical models
via the \key{trdtrc} and \key{trdmld\_trc} CPP keys.
+
+\textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested.
+In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl}
+are not working, and none of the option have been tested with variable volume ($i.e.$ \key{vvl} defined).
+
% 
@@ 1280,28 +1372,26 @@
\label{DIA_diag_harm}
A module is available to compute the amplitude and phase for tidal waves.
This diagnostic is actived with \key{diaharm}.

%namdia_harm
\namdisplay{namdia_harm}
%
Concerning the online Harmonic analysis, some parameters are available in namelist
\ngn{namdia\_harm} :

 \texttt{nit000\_han} is the first time step used for harmonic analysis

 \texttt{nitend\_han} is the last time step used for harmonic analysis

 \texttt{nstep\_han} is the time step frequency for harmonic analysis

 \texttt{nb\_ana} is the number of harmonics to analyse

 \texttt{tname} is an array with names of tidal constituents to analyse

\texttt{nit000\_han} and \texttt{nitend\_han} must be between \texttt{nit000} and \texttt{nitend} of the simulation.
+A module is available to compute the amplitude and phase of tidal waves.
+This online Harmonic analysis is actived with \key{diaharm}.
+Some parameters are available in namelist \ngn{namdia\_harm} :
+
+ \np{nit000\_han} is the first time step used for harmonic analysis
+
+ \np{nitend\_han} is the last time step used for harmonic analysis
+
+ \np{nstep\_han} is the time step frequency for harmonic analysis
+
+ \np{nb\_ana} is the number of harmonics to analyse
+
+ \np{tname} is an array with names of tidal constituents to analyse
+
+\np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation.
The restart capability is not implemented.
The Harmonic analysis solve this equation:
+The Harmonic analysis solve the following equation:
\begin{equation}
h_{i}  A_{0} + \sum^{nb\_ana}_{j=1}[A_{j}cos(\nu_{j}t_{j}\phi_{j})] = e_{i}
@@ 1350,6 +1440,6 @@
\label{DIA_diag_dct}
A module is available to compute the transport of volume, heat and salt through sections. This diagnostic
is actived with \key{diadct}.
+A module is available to compute the transport of volume, heat and salt through sections.
+This diagnostic is actived with \key{diadct}.
Each section is defined by the coordinates of its 2 extremities. The pathways between them are contructed
@@ 1373,9 +1463,9 @@
%
\texttt{nn\_dct}: frequency of instantaneous transports computing

\texttt{nn\_dctwri}: frequency of writing ( mean of instantaneous transports )

\texttt{nn\_debug}: debugging of the section
+\np{nn\_dct}: frequency of instantaneous transports computing
+
+\np{nn\_dctwri}: frequency of writing ( mean of instantaneous transports )
+
+\np{nn\_debug}: debugging of the section
\subsubsection{ To create a binary file containing the pathway of each section }
@@ 1508,5 +1598,5 @@
the \key{diahth} CPP key:
 the mixed layer depth (based on a density criterion, \citet{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
+ the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
 the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})
Index: trunk/DOC/TexFiles/Chapters/Chap_DIU.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_DIU.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_DIU.tex (revision 6289)
@@ 13,6 +13,6 @@
Code to produce an estimate of the diurnal warming and cooling of the sea surface skin
temperature (skin SST) is found in the DIU directory. The skin
temperature can be split into three parts:
+temperature (skin SST) is found in the DIU directory.
+The skin temperature can be split into three parts:
\begin{itemize}
\item A foundation SST which is free from diurnal warming.
@@ 25,19 +25,18 @@
Models are provided for both the warm layer, diurnal\_bulk.F90, and the cool skin,
cool\_skin.F90. Foundation SST is not considered as it can be obtained
either from the main NEMO model (i.e. from the temperature of the top few model levels)
or from
some other source. It must be noted that both the cool skin and
warm layer models produce estimates of the change in temperature ($\Delta T_{\rm{cs}}$
and $\Delta T_{\rm{wl}}$) and both must
be added to a foundation SST to obtain the true skin temperature.
+either from the main NEMO model ($i.e.$ from the temperature of the top few model levels)
+or from some other source.
+It must be noted that both the cool skin and warm layer models produce estimates of
+the change in temperature ($\Delta T_{\rm{cs}}$ and $\Delta T_{\rm{wl}}$)
+and both must be added to a foundation SST to obtain the true skin temperature.
Both the cool skin and warm layer models are controlled through the namelist `namdiu':
+Both the cool skin and warm layer models are controlled through the namelist \ngn{namdiu}:
\namdisplay{namdiu}
This namelist contains only two variables:
\begin{description}
\item[ln\_diurnal] A logical switch for turning on/off both the cool skin and warm layer.
\item[ln\_diurnal\_only] A logical switch which if .TRUE. will run the diurnal model
without the other dynamical parts of NEMO. ln\_diurnal\_only must be
.FALSE. if ln\_diurnal is .FALSE.
+\item[\np{ln\_diurnal}] A logical switch for turning on/off both the cool skin and warm layer.
+\item[\np{ln\_diurnal\_only}] A logical switch which if .TRUE. will run the diurnal model
+without the other dynamical parts of NEMO.
+\np{ln\_diurnal\_only} must be .FALSE. if \np{ln\_diurnal} is .FALSE.
\end{description}
@@ 46,11 +45,12 @@
output if they are specified in the iodef.xml file.
Initialisation is through the restart file. Specifically the code will expect the presence of the 2D variable ``Dsst'' to initialise the warm layer. The cool skin model, which is determined purely by the instantaneous fluxes, has no initialisation variable.
+Initialisation is through the restart file. Specifically the code will expect
+the presence of the 2D variable ``Dsst'' to initialise the warm layer.
+The cool skin model, which is determined purely by the instantaneous fluxes,
+has no initialisation variable.
%===============================================================

\section{Warm Layer model}
\label{warm_layer_sec}

%===============================================================
@@ 81,5 +81,5 @@
(\ref{ecmwf1}) is the instantaneous total thermal energy
flux into
the diurnal layer, i.e.
+the diurnal layer, $i.e.$
\begin{equation}
Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,} \label{e_flux_eqn}
@@ 105,5 +105,5 @@
where $\zeta=\frac{D_T}{L}$. It is clear that the first derivative of
(\ref{stab_func_eqn}), and thus of (\ref{ecmwf1}),
is discontinuous at $\zeta=0$ (i.e. $Q\rightarrow0$ in equation (\ref{ecmwf2})).
+is discontinuous at $\zeta=0$ ($i.e.$ $Q\rightarrow0$ in equation (\ref{ecmwf2})).
The two terms on the right hand side of (\ref{ecmwf1}) represent different processes.
Index: trunk/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_DOM.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_DOM.tex (revision 6289)
@@ 138,23 +138,14 @@
and $f$points, and its divergence defined at $t$points:
\begin{eqnarray} \label{Eq_DOM_curl}
 \nabla \times {\rm {\bf A}}\equiv &
+ \nabla \times {\rm{\bf A}}\equiv &
\frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] \delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right) &\ \mathbf{i} \\
+& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1 \right] \delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right) &\ \mathbf{j} \\
+& \frac{1}{e_{1f} \,e_{2f} } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2 \right] \delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right) &\ \mathbf{k}
\end{eqnarray}
\begin{equation} \label{Eq_DOM_div}
\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}\,e_{3t}}\left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right]
 +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right]
\end{equation}

In the special case of a pure $z$coordinate system, \eqref{Eq_DOM_lap} and
\eqref{Eq_DOM_div} can be simplified. In this case, the vertical scale factor
becomes a function of the single variable $k$ and thus does not depend on the
horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to:
\begin{equation*}
\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}} \left( \delta_i \left[e_{2u}\,a_1 \right]
 +\delta_j \left[e_{1v}\, a_2 \right] \right)
 +\frac{1}{e_{3t}} \delta_k \left[ a_3 \right]
\end{equation*}
+\begin{eqnarray} \label{Eq_DOM_div}
+\nabla \cdot \rm{\bf A} \equiv
+ \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right]
+ +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right]
+\end{eqnarray}
The vertical average over the whole water column denoted by an overbar becomes
@@ 183,6 +174,6 @@
Let $a$ and $b$ be two fields defined on the mesh, with value zero inside
continental area. Using integration by parts it can be shown that the differencing
operators ($\delta_i$, $\delta_j$ and $\delta_k$) are antisymmetric linear
operators, and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$,
+operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skewsymmetric linear operators,
+and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$,
$\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear
operators, $i.e.$
@@ 425,5 +416,5 @@
The choice of the grid must be consistent with the boundary conditions specified
by the parameter \np{jperio} (see {\S\ref{LBC}).
+by \np{jperio}, a parameter found in \ngn{namcfg} namelist (see {\S\ref{LBC}).
% 
@@ 505,5 +496,5 @@
If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft
(in meters) (\ifile{isf\_draft\_meter}) is needed and all the location where the isf cavity thinnest
 than \np{rn\_isfhmin} meters are grounded (ie masked).
+ than \np{rn\_isfhmin} meters are grounded ($i.e.$ masked).
After reading the bathymetry, the algorithm for vertical grid definition differs
@@ 540,5 +531,5 @@
Three options are possible for defining the bathymetry, according to the
namelist variable \np{nn\_bathy}:
+namelist variable \np{nn\_bathy} (found in \ngn{namdom} namelist):
\begin{description}
\item[\np{nn\_bathy} = 0] a flatbottom domain is defined. The total depth $z_w (jpk)$
@@ 721,5 +712,5 @@
usually 10\%, of the default thickness $e_{3t}(jk)$).
 \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }
+\gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } }
% 
@@ 749,9 +740,15 @@
depth, since a mixed steplike and bottomfollowing representation of the
topography can be used (Fig.~\ref{Fig_z_zps_s_sps}de) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f).
The namelist parameter \np{rn\_rmax} determines the slope at which the terrainfollowing coordinate intersects the sea bed and becomes a pseudo zcoordinate. The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as the minimum and maximum depths at which the terrainfollowing vertical coordinate is calculated.

Options for stretching the coordinate are provided as examples, but care must be taken to ensure that the vertical stretch used is appropriate for the application.

The original default NEMO scoordinate stretching is available if neither of the other options are specified as true (\np{ln\_sco\_SH94}~=~false and \np{ln\_sco\_SF12}~=~false.) This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
+The namelist parameter \np{rn\_rmax} determines the slope at which the terrainfollowing coordinate intersects
+the sea bed and becomes a pseudo zcoordinate.
+The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max}
+as the minimum and maximum depths at which the terrainfollowing vertical coordinate is calculated.
+
+Options for stretching the coordinate are provided as examples, but care must be taken to ensure
+that the vertical stretch used is appropriate for the application.
+
+The original default NEMO scoordinate stretching is available if neither of the other options
+are specified as true (\np{ln\_s\_SH94}~=~false and \np{ln\_s\_SF12}~=~false).
+This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
\begin{equation}
@@ 760,5 +757,7 @@
\end{equation}
where $s_{min}$ is the depth at which the scoordinate stretching starts and allows a zcoordinate to placed on top of the stretched coordinate, and z is the depth (negative down from the asea surface).
+where $s_{min}$ is the depth at which the scoordinate stretching starts and
+allows a zcoordinate to placed on top of the stretched coordinate,
+and z is the depth (negative down from the asea surface).
\begin{equation}
@@ 775,5 +774,6 @@
\end{equation}
A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_sco\_SH94}~=~true), is also available and is more commonly used for shelf seas modelling:
+A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94}
+stretching (\np{ln\_s\_SH94}~=~true), is also available and is more commonly used for shelf seas modelling:
\begin{equation}
@@ 792,10 +792,13 @@
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and
bottom control parameters such that $0\leqslant \theta \leqslant 20$, and
+where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from
+pure $\sigma$ to the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb})
+are the surface and bottom control parameters such that $0\leqslant \theta \leqslant 20$, and
$0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom
increase of the vertical resolution (Fig.~\ref{Fig_sco_function}).
Another example has been provided at version 3.5 (\np{ln\_sco\_SF12}) that allows a fixed surface resolution in an analytical terrainfollowing stretching \citet{Siddorn_Furner_OM12}. In this case the a stretching function $\gamma$ is defined such that:
+Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows
+a fixed surface resolution in an analytical terrainfollowing stretching \citet{Siddorn_Furner_OM12}.
+In this case the a stretching function $\gamma$ is defined such that:
\begin{equation}
@@ 815,5 +818,8 @@
\end{equation}
This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs}) and bottom depths. The bottom cell depth in this example is given as a function of water depth:
+This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of
+the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards
+the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs})
+and bottom depths. The bottom cell depth in this example is given as a function of water depth:
\begin{equation} \label{DOM_zb}
@@ 843,6 +849,5 @@
\label{DOM_zgr_star}
This option is described in the Report by Levier \textit{et al.} (2007), available on
the \NEMO web site.
+This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site.
%gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances
Index: trunk/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_DYN.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_DYN.tex (revision 6289)
@@ 301,10 +301,13 @@
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Note that a key point in \eqref{Eq_een_e3f} is that the averaging in the \textbf{i} and
\textbf{j} directions uses the masked vertical scale factor but is always divided by
$4$, not by the sum of the masks at the four $T$points. This preserves the continuity of
$e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and
extends by continuity the value of $e_{3f}$ into the land areas. This feature is essential for
the $z$coordinate with partial steps.
+A key point in \eqref{Eq_een_e3f} is how the averaging in the \textbf{i} and \textbf{j} directions is made.
+It uses the sum of masked tpoint vertical scale factor divided either
+by the sum of the four tpoint masks (\np{nn\_een\_e3f}~=~1),
+or just by $4$ (\np{nn\_een\_e3f}~=~true).
+The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$
+tends to zero and extends by continuity the value of $e_{3f}$ into the land areas.
+This case introduces a subgridscale topography at fpoints (with a systematic reduction of $e_{3f}$
+when a model level intercept the bathymetry) that tends to reinforce the topostrophy of the flow
+($i.e.$ the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.
Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$point as
@@ 372,4 +375,11 @@
\end{aligned} \right.
\end{equation}
+When \np{ln\_dynzad\_zts}~=~\textit{true}, a splitexplicit time stepping with 5 subtimesteps is used
+on the vertical advection term.
+This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.
+Note that in this case, a similar splitexplicit time stepping should be used on
+vertical advection of tracer to ensure a better stability,
+an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \S\ref{TRA_adv_tvd}).
+
% ================================================================
@@ 716,9 +726,17 @@
$\ $\newline %force an empty line
%%%
Options are defined through the \ngn{namdyn\_spg} namelist variables.
The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}). The main distinction is between the fixed volume case (linear free surface) and the variable volume case (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface}) the vertical scale factors $e_{3}$ are fixed in time, while they are timedependent in the nonlinear case (\S\ref{PE_free_surface}). With both linear and nonlinear free surface, external gravity waves are allowed in the equations, which imposes a very small time step when an explicit time stepping is used. Two methods are proposed to allow a longer time step for the threedimensional equations: the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}), and the splitexplicit free surface described below. The extra term introduced in the filtered method is calculated implicitly, so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.

%%%
+The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}).
+The main distinction is between the fixed volume case (linear free surface) and the variable volume case
+(nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface})
+the vertical scale factors $e_{3}$ are fixed in time, while they are timedependent in the nonlinear case
+(\S\ref{PE_free_surface}).
+With both linear and nonlinear free surface, external gravity waves are allowed in the equations,
+which imposes a very small time step when an explicit time stepping is used.
+Two methods are proposed to allow a longer time step for the threedimensional equations:
+the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}),
+and the splitexplicit free surface described below.
+The extra term introduced in the filtered method is calculated implicitly,
+so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
@@ 734,5 +752,4 @@
implicitly, so that a solver is used to compute it. As a consequence the update of the $next$
velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.

@@ 777,5 +794,6 @@
$\rdt_e = \rdt / nn\_baro$. This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}=true)
considering that the stability of the barotropic system is essentially controled by external waves propagation.
Maximum allowed Courant number is in that case time independent, and easily computed online from the input bathymetry.
+Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry.
+Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn\_bt\_cmax}.
%%%
@@ 800,5 +818,5 @@
Schematic of the splitexplicit time stepping scheme for the external
and internal modes. Time increases to the right. In this particular exemple,
a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_filt=1$) and $nn\_baro=5$.
+a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_flt=1$) and $nn\_baro=5$.
Internal mode time steps (which are also the model time steps) are denoted
by $t\rdt$, $t$ and $t+\rdt$. Variables with $k$ superscript refer to instantaneous barotropic variables,
@@ 806,7 +824,7 @@
The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged
transports to advect tracers.
a) Forward time integration: \np{ln\_bt\_fw}=true, \np{ln\_bt\_ave}=true.
b) Centred time integration: \np{ln\_bt\_fw}=false, \np{ln\_bt\_ave}=true.
c) Forward time integration with no time filtering (POMlike scheme): \np{ln\_bt\_fw}=true, \np{ln\_bt\_ave}=false. }
+a) Forward time integration: \np{ln\_bt\_fw}=true, \np{ln\_bt\_av}=true.
+b) Centred time integration: \np{ln\_bt\_fw}=false, \np{ln\_bt\_av}=true.
+c) Forward time integration with no time filtering (POMlike scheme): \np{ln\_bt\_fw}=true, \np{ln\_bt\_av}=false. }
\end{center} \end{figure}
%> > > > > > > > > > > > > > > > > > > > > > > > > > > >
@@ 814,5 +832,5 @@
In the default case (\np{ln\_bt\_fw}=true), the external mode is integrated
between \textit{now} and \textit{after} baroclinic timesteps (Fig.~\ref{Fig_DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic
quantities (\np{ln\_bt\_ave}=true). In that case, the integration is extended slightly beyond \textit{after} time step to provide time filtered quantities.
+quantities (\np{ln\_bt\_av}=true). In that case, the integration is extended slightly beyond \textit{after} time step to provide time filtered quantities.
These are used for the subsequent initialization of the barotropic mode in the following baroclinic step.
Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme,
@@ 835,5 +853,5 @@
%%%
One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_ave}=false).
+One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_av}=false).
In that case, external mode equations are continuous in time, ie they are not reinitialized when starting a new
substepping sequence. This is the method used so far in the POM model, the stability being maintained by refreshing at (almost)
@@ 987,4 +1005,18 @@
At the lateral boundaries either free slip, no slip or partial slip boundary
conditions are applied according to the user's choice (see Chap.\ref{LBC}).
+
+\gmcomment{
+Hyperviscous operators are frequently used in the simulation of turbulent flows to control
+the dissipation of unresolved small scale features.
+Their primary role is to provide strong dissipation at the smallest scale supported by the grid
+while minimizing the impact on the larger scale features.
+Hyperviscous operators are thus designed to be more scale selective than the traditional,
+physically motivated Laplace operator.
+In finite difference methods, the biharmonic operator is frequently the method of choice to achieve
+this scale selective dissipation since its damping time ($i.e.$ its spin down time)
+scale like $\lambda^{4}$ for disturbances of wavelength $\lambda$
+(so that short waves damped more rapidelly than long ones),
+whereas the Laplace operator damping time scales only like $\lambda^{2}$.
+}
% ================================================================
@@ 1156,12 +1188,19 @@
Besides the surface and bottom stresses (see the above section) which are
introduced as boundary conditions on the vertical mixing, two other forcings
enter the dynamical equations.

One is the effect of atmospheric pressure on the ocean dynamics.
Another forcing term is the tidal potential.
Both of which will be introduced into the reference version soon.

\gmcomment{atmospheric pressure is there!!!! include its description }
+introduced as boundary conditions on the vertical mixing, three other forcings
+may enter the dynamical equations by affecting the surface pressure gradient.
+
+(1) When \np{ln\_apr\_dyn}~=~true (see \S\ref{SBC_apr}), the atmospheric pressure is taken
+into account when computing the surface pressure gradient.
+
+(2) When \np{ln\_tide\_pot}~=~true and \key{tide} is defined (see \S\ref{SBC_tide}),
+the tidal potential is taken into account when computing the surface pressure gradient.
+
+(3) When \np{nn\_ice\_embd}~=~2 and LIM or CICE is used ($i.e.$ when the seaice is embedded in the ocean),
+the snowice mass is taken into account when computing the surface pressure gradient.
+
+
+\gmcomment{ missing : the lateral boundary condition !!! another external forcing
+ }
% ================================================================
@@ 1211,85 +1250,2 @@
% ================================================================
% Neptune effect
% ================================================================
\section [Neptune effect (\textit{dynnept})]
 {Neptune effect (\mdl{dynnept})}
\label{DYN_nept}

The "Neptune effect" (thus named in \citep{HollowayOM86}) is a
parameterisation of the potentially large effect of topographic form stress
(caused by eddies) in driving the ocean circulation. Originally developed for
lowresolution models, in which it was applied via a Laplacian (secondorder)
diffusionlike term in the momentum equation, it can also be applied in eddy
permitting or resolving models, in which a more scaleselective bilaplacian
(fourthorder) implementation is preferred. This mechanism has a
significant effect on boundary currents (including undercurrents), and the
upwelling of deep water near continental shelves.

The theoretical basis for the method can be found in
\citep{HollowayJPO92}, including the explanation of why form stress is not
necessarily a drag force, but may actually drive the flow.
\citep{HollowayJPO94} demonstrate the effects of the parameterisation in
the GFDLMOM model, at a horizontal resolution of about 1.8 degrees.
\citep{HollowayOM08} demonstrate the biharmonic version of the
parameterisation in a global run of the POP model, with an average horizontal
grid spacing of about 32km.

The NEMO implementation is a simplified form of that supplied by
Greg Holloway, the testing of which was described in \citep{HollowayJGR09}.
The major simplification is that a time invariant Neptune velocity
field is assumed. This is computed only once, during startup, and
made available to the rest of the code via a module. Vertical
diffusive terms are also ignored, and the model topography itself
is used, rather than a separate topographic dataset as in
\citep{HollowayOM08}. This implementation is only in the isolevel
formulation, as is the case anyway for the bilaplacian operator.

The velocity field is derived from a transport stream function given by:

\begin{equation} \label{Eq_dynnept_sf}
\psi = fL^2H
\end{equation}

where $L$ is a latitudedependant length scale given by:

\begin{equation} \label{Eq_dynnept_ls}
L = l_1 + (l_2 l_1)\left ( {1 + \cos 2\phi \over 2 } \right )
\end{equation}

where $\phi$ is latitude and $l_1$ and $l_2$ are polar and equatorial length scales respectively.
Neptune velocity components, $u^*$, $v^*$ are derived from the stremfunction as:

\begin{equation} \label{Eq_dynnept_vel}
u^* = {1\over H} {\partial \psi \over \partial y}\ \ \ ,\ \ \ v^* = {1\over H} {\partial \psi \over \partial x}
\end{equation}

\smallskip
%namdom
\namdisplay{namdyn_nept}
%
\smallskip

The Neptune effect is enabled when \np{ln\_neptsimp}=true (default=false).
\np{ln\_smooth\_neptvel} controls whether a scaleselective smoothing is applied
to the Neptune effect flow field (default=false) (this smoothing method is as
used by Holloway). \np{rn\_tslse} and \np{rn\_tslsp} are the equatorial and
polar values respectively of the lengthscale parameter $L$ used in determining
the Neptune stream function \eqref{Eq_dynnept_sf} and \eqref{Eq_dynnept_ls}.
Values at intermediate latitudes are given by a cosine fit, mimicking the
variation of the deformation radius with latitude. The default values of 12km
and 3km are those given in \citep{HollowayJPO94}, appropriate for a coarse
resolution model. The finer resolution study of \citep{HollowayOM08} increased
the values of L by a factor of $\sqrt 2$ to 17km and 4.2km, thus doubling the
stream function for a given topography.

The simple formulation for ($u^*$, $v^*$) can give unacceptably large velocities
in shallow water, and \citep{HollowayOM08} add an offset to the depth in the
denominator to control this problem. In this implementation we offer instead (at
the suggestion of G. Madec) the option of ramping down the Neptune flow field to
zero over a finite depth range. The switch \np{ln\_neptramp} activates this
option (default=false), in which case velocities at depths greater than
\np{rn\_htrmax} are unaltered, but ramp down linearly with depth to zero at a
depth of \np{rn\_htrmin} (and shallower).

% ================================================================
Index: trunk/DOC/TexFiles/Chapters/Chap_LBC.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_LBC.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_LBC.tex (revision 6289)
@@ 1,4 +1,4 @@
% ================================================================
% Chapter � Lateral Boundary Condition (LBC)
+% Chapter — Lateral Boundary Condition (LBC)
% ================================================================
\chapter{Lateral Boundary Condition (LBC) }
@@ 127,22 +127,4 @@
can be quite shallow.
The alternative numerical implementation of the noslip boundary conditions for an
arbitrary coast line of \citet{Shchepetkin1996} is also available through the
\key{noslip\_accurate} CPP key. It is based on a fourth order evaluation of the shear at the
coast which, in turn, allows a true second order scheme in the interior of the domain
($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical
scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a
technique considerably improves the quality of the numerical solution. In \NEMO, such
spectacular improvements have not been found in the halfdegree global ocean
(ORCA05), but significant reductions of numerically induced coastal upwellings were
found in an eddy resolving simulation of the Alboran Sea \citep{Olivier_PhD01}.
Nevertheless, since a noslip boundary condition is not recommended in an eddy
permitting or resolving simulation \citep{Penduff_al_OS07}, the use of this option is also
not recommended.

In practice, the noslip accurate option changes the way the curl is evaluated at the
coast (see \mdl{divcur} module), and requires the nature of each coastline grid point
(convex or concave corners, straight northsouth or eastwest coast) to be specified.
This is performed in routine \rou{dom\_msk\_nsa} in the \mdl{domask} module.
% ================================================================
@@ 204,10 +186,11 @@
% North fold (\textit{jperio = 3 }to $6)$
% 
\subsection{Northfold (\textit{jperio = 3 }to $6)$ }
+\subsection{Northfold (\textit{jperio = 3 }to $6$) }
\label{LBC_north_fold}
The north fold boundary condition has been introduced in order to handle the north
boundary of a threepolar ORCA grid. Such a grid has two poles in the northern hemisphere.
\colorbox{yellow}{to be completed...}
+boundary of a threepolar ORCA grid. Such a grid has two poles in the northern hemisphere
+(Fig.\ref{Fig_MISC_ORCA_msh}, and thus requires a specific treatment illustrated in Fig.\ref{Fig_North_Fold_T}.
+Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition.
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ 250,10 +233,9 @@
ocean model. Second order finite difference schemes lead to local discrete
operators that depend at the very most on one neighbouring point. The only
nonlocal computations concern the vertical physics (implicit diffusion, 1.5
+nonlocal computations concern the vertical physics (implicit diffusion,
turbulent closure scheme, ...) (delocalization over the whole water column),
and the solving of the elliptic equation associated with the surface pressure
gradient computation (delocalization over the whole horizontal domain).
Therefore, a pencil strategy is used for the data substructuration
\gmcomment{no idea what this means!}
: the 3D initial domain is laid out on local processor
memories following a 2D horizontal topological splitting. Each subdomain
@@ 264,16 +246,14 @@
phase starts: each processor sends to its neighbouring processors the update
values of the points corresponding to the interior overlapping area to its
neighbouring subdomain (i.e. the innermost of the two overlapping rows).
The communication is done through message passing. Usually the parallel virtual
language, PVM, is used as it is a standard language available on nearly all
MPP computers. More specific languages (i.e. computer dependant languages)
can be easily used to speed up the communication, such as SHEM on a T3E
computer. The data exchanges between processors are required at the very
+neighbouring subdomain ($i.e.$ the innermost of the two overlapping rows).
+The communication is done through the Message Passing Interface (MPI).
+The data exchanges between processors are required at the very
place where lateral domain boundary conditions are set in the monodomain
computation (\S III.10c): the lbc\_lnk routine which manages such conditions
is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP
computer (\key{mpp\_mpi} defined). It has to be pointed out that when using
the MPP version of the model, the eastwest cyclic boundary condition is done
implicitly, whilst the southsymmetric boundary condition option is not available.
+computation : the \rou{lbc\_lnk} routine (found in \mdl{lbclnk} module)
+which manages such conditions is interfaced with routines found in \mdl{lib\_mpp} module
+when running on an MPP computer ($i.e.$ when \key{mpp\_mpi} defined).
+It has to be pointed out that when using the MPP version of the model,
+the eastwest cyclic boundary condition is done implicitly,
+whilst the southsymmetric boundary condition option is not available.
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ 285,12 +265,10 @@
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
In the standard version of the OPA model, the splitting is regular and arithmetic.
 the iaxis is divided by \jp{jpni} and the jaxis by \jp{jpnj} for a number of processors
 \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in
 \mdl{par\_oce}). Each processor is independent and without message passing
 or synchronous process
 \gmcomment{how does a synchronous process relate to this?},
 programs run alone and access just its own local memory. For this reason, the
 main model dimensions are now the local dimensions of the subdomain (pencil)
+In the standard version of \NEMO, the splitting is regular and arithmetic.
+The iaxis is divided by \jp{jpni} and the jaxis by \jp{jpnj} for a number of processors
+\jp{jpnij} most often equal to $jpni \times jpnj$ (parameters set in
+ \ngn{nammpp} namelist). Each processor is independent and without message passing
+ or synchronous process, programs run alone and access just its own local memory.
+ For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil)
that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal
domain and the overlapping rows. The number of rows to exchange (known as
@@ 304,12 +282,8 @@
where \jp{jpni}, \jp{jpnj} are the number of processors following the i and jaxis.
\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and
no eastwest cyclic boundary conditions.}

One also defines variables nldi and nlei which correspond to the internal
domain bounds, and the variables nimpp and njmpp which are the position
of the (1,1) gridpoint in the global domain. An element of $T_{l}$, a local array
(subdomain) corresponds to an element of $T_{g}$, a global array
(whole domain) by the relationship:
+One also defines variables nldi and nlei which correspond to the internal domain bounds,
+and the variables nimpp and njmpp which are the position of the (1,1) gridpoint in the global domain.
+An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$,
+a global array (whole domain) by the relationship:
\begin{equation} \label{Eq_lbc_nimpp}
T_{g} (i+nimpp1,j+njmpp1,k) = T_{l} (i,j,k),
@@ 320,6 +294,5 @@
nproc. In the standard version, a processor has no more than four neighbouring
processors named nono (for north), noea (east), noso (south) and nowe (west)
and two variables, nbondi and nbondj, indicate the relative position of the processor
\colorbox{yellow}{(see Fig.IV.3)}:
+and two variables, nbondi and nbondj, indicate the relative position of the processor :
\begin{itemize}
\item nbondi = 1 an east neighbour, no west processor,
@@ 332,6 +305,4 @@
processor on its overlapping row, and sends the data issued from internal
domain corresponding to the overlapping row of the other processor.

\colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos }
@@ 343,27 +314,20 @@
global ocean where more than 50 \% of points are land points. For this reason, a
preprocessing tool can be used to choose the mpp domain decomposition with a
maximum number of only land points processors, which can then be eliminated.
(For example, the mpp\_optimiz tools, available from the DRAKKAR web site.)
+maximum number of only land points processors, which can then be eliminated (Fig. \ref{Fig_mppini2})
+(For example, the mpp\_optimiz tools, available from the DRAKKAR web site).
This optimisation is dependent on the specific bathymetry employed. The user
then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with
$jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj  jpnij$
land processors. When those parameters are specified in module \mdl{par\_oce},
+land processors. When those parameters are specified in \ngn{nammpp} namelist,
the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound,
nono, noea,...) so that the landonly processors are not taken into account.
\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp
+\gmcomment{Note that the inimpp2 routine is general so that the original inimpp
routine should be suppressed from the code.}
When land processors are eliminated, the value corresponding to these locations in
the model output files is zero. Note that this is a problem for a mesh output file written
by such a model configuration, because model users often divide by the scale factors
($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be
best not to eliminate land processors when running the model especially to write the
mesh files as outputs (when \np{nn\_msh} namelist parameter differs from 0).
%%
\gmcomment{Steven : dont understand this, no land processor means no output file
covering this part of globe; its only when files are stitched together into one that you
can leave a hole}
%%
+the model output files is undefined. Note that this is a problem for the meshmask file
+which requires to be defined over the whole domain. Therefore, user should not eliminate
+land processors when creating a meshmask file ($i.e.$ when setting a nonzero value to \np{nn\_msh}).
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ 380,360 +344,4 @@
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>

% ================================================================
% Open Boundary Conditions
% ================================================================
\section{Open Boundary Conditions (\key{obc}) (OBC)}
\label{LBC_obc}
%nam_obc 
% nobc_dta = 0 ! = 0 the obc data are equal to the initial state
% ! = 1 the obc data are read in 'obc .dta' files
% rn_dpein = 1. ! ???
% rn_dpwin = 1. ! ???
% rn_dpnin = 30. ! ???
% rn_dpsin = 1. ! ???
% rn_dpeob = 1500. ! time relaxation (days) for the east open boundary
% rn_dpwob = 15. ! " " for the west open boundary
% rn_dpnob = 150. ! " " for the north open boundary
% rn_dpsob = 15. ! " " for the south open boundary
% ln_obc_clim = .true. ! climatological obc data files (default T)
% ln_vol_cst = .true. ! total volume conserved
\namdisplay{namobc}

It is often necessary to implement a model configuration limited to an oceanic
region or a basin, which communicates with the rest of the global ocean through
''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a
computational border where the aim of the calculations is to allow the perturbations
generated inside the computational domain to leave it without deterioration of the
inner model solution. However, an open boundary also has to let information from
the outer ocean enter the model and should support inflow and outflow conditions.

The open boundary package OBC is the first open boundary option developed in
NEMO (originally in OPA8.2). It allows the user to
\begin{itemize}
\item tell the model that a boundary is ''open'' and not closed by a wall, for example
by modifying the calculation of the divergence of velocity there;
\item impose values of tracers and velocities at that boundary (values which may
be taken from a climatology): this is the``fixed OBC'' option.
\item calculate boundary values by a sophisticated algorithm combining radiation
and relaxation (``radiative OBC'' option)
\end{itemize}

Options are defined through the \ngn{namobc} namelist variables.
The package resides in the OBC directory. It is described here in four parts: the
boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at
the boundaries (module \mdl{obcdta}), the radiation algorithm involving the
namelist and module \mdl{obcrad}, and a brief presentation of boundary update
and restart files.

%
\subsection{Boundary geometry}
\label{OBC_geom}
%
First one has to realize that open boundaries may not necessarily be located
at the extremities of the computational domain. They may exist in the middle
of the domain, for example at Gibraltar Straits if one wants to avoid including
the Mediterranean in an Atlantic domain. This flexibility has been found necessary
for the CLIPPER project \citep{Treguier_al_JGR01}. Because of the complexity of the
geometry of ocean basins, it may even be necessary to have more than one
''west'' open boundary, more than one ''north'', etc. This is not possible with
the OBC option: only one open boundary of each kind, west, east, south and
north is allowed; these names refer to the grid geometry (not to the direction
of the geographical ''west'', ''east'', etc).

The open boundary geometry is set by a series of parameters in the module
\mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east}
(true if an east open boundary exists), \jp{jpieob} the $i$index along which
the eastern open boundary (eob) is located, \jp{jpjed} the $j$index at which
it starts, and \jp{jpjef} the $j$index where it ends (note $d$ is for ''d\'{e}but''
and $f$ for ''fin'' in French). Similar parameters exist for the west, south and
north cases (Table~\ref{Tab_obc_param}).


%TABLE
\begin{table}[htbp] \begin{center} \begin{tabular}{lccc}
\hline
Boundary and & Constant index & Starting index (d\'{e}but) & Ending index (fin) \\
Logical flag & & & \\
\hline
West & \jp{jpiwob} $>= 2$ & \jp{jpjwd}$>= 2$ & \jp{jpjwf}<= \np{jpjglo}1 \\
lp\_obc\_west & $i$index of a $u$ point & $j$ of a $T$ point &$j$ of a $T$ point \\
\hline
East & \jp{jpieob}$<=$\np{jpiglo}2&\jp{jpjed} $>= 2$ & \jp{jpjef}$<=$ \np{jpjglo}1 \\
 lp\_obc\_east & $i$index of a $u$ point & $j$ of a $T$ point & $j$ of a $T$ point \\
\hline
South & \jp{jpjsob} $>= 2$ & \jp{jpisd} $>= 2$ & \jp{jpisf}$<=$\np{jpiglo}1 \\
lp\_obc\_south & $j$index of a $v$ point & $i$ of a $T$ point & $i$ of a $T$ point \\
\hline
North & \jp{jpjnob} $<=$ \np{jpjglo}2& \jp{jpind} $>= 2$ & \jp{jpinf}$<=$\np{jpiglo}1 \\
lp\_obc\_north & $j$index of a $v$ point & $i$ of a $T$ point & $i$ of a $T$ point \\
\hline
\end{tabular} \end{center}
\caption{ \label{Tab_obc_param}
Names of different indices relating to the open boundaries. In the case
of a completely open ocean domain with four ocean boundaries, the parameters
take exactly the values indicated.}
\end{table}
%

The open boundaries must be along coordinate lines. On the Cgrid, the boundary
itself is along a line of normal velocity points: $v$ points for a zonal open boundary
(the south or north one), and $u$ points for a meridional open boundary (the west
or east one). Another constraint is that there still must be a row of masked points
all around the domain, as if the domain were a closed basin (unless periodic conditions
are used together with open boundary conditions). Therefore, an open boundary
cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also,
the open boundary algorithm involves calculating the normal velocity points situated
just on the boundary, as well as the tangential velocity and temperature and salinity
just outside the boundary. This means that for a west/south boundary, normal
velocities and temperature are calculated at the same index \jp{jpiwob} and
\jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is
calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is
at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob}
cannot be bigger than \jp{jpiglo}2, \jp{jpjglo}2.


The starting and ending indices are to be thought of as $T$ point indices: in many
cases they indicate the first land $T$point, at the extremity of an open boundary
(the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example
of a northern open boundary). All indices are relative to the global domain. In the
free surface case it is possible to have ``ocean corners'', that is, an open boundary
starting and ending in the ocean.

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!t] \begin{center}
\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf}
\caption{ \label{Fig_obc_north}
Localization of the North open boundary points.}
\end{center} \end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Although not compulsory, it is highly recommended that the bathymetry in the
vicinity of an open boundary follows the following rule: in the direction perpendicular
to the open line, the water depth should be constant for 4 grid points. This is in
order to ensure that the radiation condition, which involves model variables next
to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we
indicate by an $=$ symbol, the points which should have the same depth. It means
that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure
why cylindrical}. The line behind the open $T$line must be 0 in the bathymetry file
(as shown on Fig.\ref{Fig_obc_north} for example).

%
\subsection{Boundary data}
\label{OBC_data}

It is necessary to provide information at the boundaries. The simplest case is
when this information does not change in time and is equal to the initial conditions
(namelist variable \np{nn\_obcdta}=0). This is the case for the standard configuration
EEL5 with open boundaries. When (\np{nn\_obcdta}=1), open boundary information
is read from netcdf files. For convenience the input files are supposed to be similar
to the ''history'' NEMO output files, for dimension names and variable names.
Open boundary arrays must be dimensioned according to the parameters of table~
\ref{Tab_obc_param}: for example, at the western boundary, arrays have a
dimension of \jp{jpwf}\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical.

When ocean observations are used to generate the boundary data (a hydrographic
section for example, as in \citet{Treguier_al_JGR01}) it happens often that only the velocity
normal to the boundary is known, which is the reason why the initial OBC code
assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be
specified. As more and more global model solutions and ocean analysis products
become available, it will be possible to provide information about all the variables
(including the tangential velocity) so that the specification of four variables at each
boundaries will become standard. For the sea surface height, one must distinguish
between the filtered free surface case and the timesplitting or explicit treatment of
the free surface.
 In the first case, it is assumed that the user does not wish to represent high
 frequency motions such as tides. The boundary condition is thus one of zero
 normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}.
No information other than the total velocity needs to be provided at the open
boundaries in that case. In the other two cases (time splitting or explicit free surface),
the user must provide barotropic information (sea surface height and barotropic
velocities) and the use of the Flather algorithm for barotropic variables is
recommanded. However, this algorithm has not yet been fully tested and bugs
remain in NEMO v2.3. Users should read the code carefully before using it. Finally,
in the case of the rigid lid approximation the barotropic streamfunction must be
provided, as documented in \citet{Treguier_al_JGR01}). This option is no longer
recommended but remains in NEMO V2.3.

One frequently encountered case is when an open boundary domain is constructed
from a global or larger scale NEMO configuration. Assuming the domain corresponds
to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the
small domain can be created by using the following netcdf utility on the global files:
ncks F $d\;x,ib,ie$ $d\;y,jb,je$ (part of the nco series of utilities,
see their \href{http://nco.sourceforge.net}{website}).
The open boundary files can be constructed using ncks
commands, following table~\ref{Tab_obc_ind}.

%TABLE
\begin{table}[htbp] \begin{center} \begin{tabular}{lccccc}
\hline
OBC & Variable & file name & Index & Start & end \\
West & T,S & obcwest\_TS.nc & $ib$+1 & $jb$+1 & $je1$ \\
 & U & obcwest\_U.nc & $ib$+1 & $jb$+1 & $je1$ \\
 & V & obcwest\_V.nc & $ib$+1 & $jb$+1 & $je1$ \\
\hline
East & T,S & obceast\_TS.nc & $ie$1 & $jb$+1 & $je1$ \\
 & U & obceast\_U.nc & $ie$2 & $jb$+1 & $je1$ \\
 & V & obceast\_V.nc & $ie$1 & $jb$+1 & $je1$ \\
\hline
South & T,S & obcsouth\_TS.nc & $jb$+1 & $ib$+1 & $ie1$ \\
 & U & obcsouth\_U.nc & $jb$+1 & $ib$+1 & $ie1$ \\
 & V & obcsouth\_V.nc & $jb$+1 & $ib$+1 & $ie1$ \\
\hline
North & T,S & obcnorth\_TS.nc & $je$1 & $ib$+1 & $ie1$ \\
 & U & obcnorth\_U.nc & $je$1 & $ib$+1 & $ie1$ \\
 & V & obcnorth\_V.nc & $je$2 & $ib$+1 & $ie1$ \\
\hline
\end{tabular} \end{center}
\caption{ \label{Tab_obc_ind}
Requirements for creating open boundary files from a global configuration,
appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the
$i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global
configuration, starting and ending with the $j$ or $i$ indices indicated.
For example, to generate file obcnorth\_V.nc, use the command ncks
$F$ $d\;y,je2$ $d\;x,ib+1,ie1$ }
\end{table}
%

It is assumed that the open boundary files contain the variables for the period of
the model integration. If the boundary files contain one time frame, the boundary
data is held fixed in time. If the files contain 12 values, it is assumed that the input
is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim}
=true). The case of an arbitrary number of time frames is not yet implemented
correctly; the user is required to write his own code in the module \mdl{obc\_dta}
to deal with this situation.

\subsection{Radiation algorithm}
\label{OBC_rad}

The art of open boundary management consists in applying a constraint strong
enough that the inner domain "feels" the rest of the ocean, but weak enough
that perturbations are allowed to leave the domain with minimum false reflections
of energy. The constraints are specified separately at each boundary as time
scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days)
by namelist parameters such as \np{rn\_dpein}, \np{rn\_dpeob} for the eastern open
boundary for example. When both time scales are zero for a given boundary
($e.g.$ for the western boundary, \jp{lp\_obc\_west}=true, \np{rn\_dpwob}=0 and
\np{rn\_dpwin}=0) this means that the boundary in question is a ''fixed '' boundary
where the solution is set exactly by the boundary data. This is not recommended,
except in combination with increased viscosity in a ''sponge'' layer next to the
boundary in order to avoid spurious reflections.


The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output}
algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is
nonzero. It has been developed and tested in the SPEM model and its
successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an
$s$coordinate model on an Arakawa Cgrid. Although the algorithm has
been numerically successful in the CLIPPER Atlantic models, the physics
do not work as expected \citep{Treguier_al_JGR01}. Users are invited to consider
open boundary conditions (OBC hereafter) with some scepticism
\citep{Durran2001, Blayo2005}.

The first part of the algorithm calculates a phase velocity to determine
whether perturbations tend to propagate toward, or away from, the
boundary. Let us consider a model variable $\phi$.
The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$,
in the directions normal and tangential to the boundary are
\begin{equation} \label{Eq_obc_cphi}
C_{\phi x} = \frac{ \phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x}
\;\;\;\;\; \;\;\;
C_{\phi y} = \frac{ \phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}.
\end{equation}
Following \citet{Treguier_al_JGR01} and \citet{Marchesiello2001} we retain only
the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$
(but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in
the expression for $C_{\phi x}$).

The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998},
takes into account the two rows of grid points situated inside the domain
next to the boundary, and the three previous time steps ($n$, $n1$,
and $n2$). The same equation can then be discretized at the boundary at
time steps $n1$, $n$ and $n+1$ \gmcomment{since the original was three timelevel}
in order to extrapolate for the new boundary value $\phi^{n+1}$.

In the open boundary algorithm as implemented in NEMO v2.3, the new boundary
values are updated differently depending on the sign of $C_{\phi x}$. Let us take
an eastern boundary as an example. The solution for variable $\phi$ at the
boundary is given by a generalized wave equation with phase velocity $C_{\phi}$,
with the addition of a relaxation term, as:
\begin{eqnarray}
\phi_{t} & = & C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}\phi)
 \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\
\phi_{t} & = & \frac{1}{\tau_{i}} (\phi_{c}\phi)
\;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi}
\end{eqnarray}
where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary
data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio
$\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward
propagation), the radiation condition (\ref{Eq_obc_rado}) is used.
When $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is
used with a strong relaxation to climatology (usually $\tau_{i}=\np{rn\_dpein}=$1~day).
Equation (\ref{Eq_obc_radi}) is solved with a Euler timestepping scheme. As a
consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent
to a fixed boundary condition. A time scale of one day is usually a good compromise
which guarantees that the inflow conditions remain close to climatology while ensuring
numerical stability.

In the case of a western boundary located in the Eastern Atlantic, \citet{Penduff_al_JGR00}
have been able to implement the radiation algorithm without any boundary data,
using persistence from the previous time step instead. This solution has not worked
in other cases \citep{Treguier_al_JGR01}, so that the use of boundary data is recommended.
Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to
maintain a weak relaxation to climatology. The time step is usually chosen so as to
be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}).

The radiation condition is applied to the model variables: temperature, salinity,
tangential and normal velocities. For normal and tangential velocities, $u$ and $v$,
radiation is applied with phase velocities calculated from $u$ and $v$ respectively.
For the radiation of tracers, we use the phase velocity calculated from the tangential
velocity in order to avoid calculating too many independent radiation velocities and
because tangential velocities and tracers have the same position along the boundary
on a Cgrid.

\subsection{Domain decomposition (\key{mpp\_mpi})}
\label{OBC_mpp}
When \key{mpp\_mpi} is active in the code, the computational domain is divided
into rectangles that are attributed each to a different processor. The open boundary
code is ``mppcompatible'' up to a certain point. The radiation algorithm will not
work if there is an mpp subdomain boundary parallel to the open boundary at the
index of the boundary, or the grid point after (outside), or three grid points before
(inside). On the other hand, there is no problem if an mpp subdomain boundary
cuts the open boundary perpendicularly. These geometrical limitations must be
checked for by the user (there is no safeguard in the code).
The general principle for the open boundary mpp code is that loops over the open
boundaries {not sure what this means} are performed on local indices (nie0,
nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module
\mdl{obc\_ini}. Those indices have relevant values on the processors that contain
a segment of an open boundary. For processors that do not include an open
boundary segment, the indices are such that the calculations within the loops are
not performed.
\gmcomment{I dont understand most of the last few sentences}

Arrays of climatological data that are read from files are seen by all processors
and have the same dimensions for all (for instance, for the eastern boundary,
uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation
are local to each processor (uebnd(jpj,jpk,3,3) for instance). This allowed the
CLIPPER model for example, to save on memory where the eastern boundary
crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}\jp{jpjed}+1).

\subsection{Volume conservation}
\label{OBC_vol}

It is necessary to control the volume inside a domain when using open boundaries.
With fixed boundaries, it is enough to ensure that the total inflow/outflow has
reasonable values (either zero or a value compatible with an observed volume
balance). When using radiative boundary conditions it is necessary to have a
volume constraint because each open boundary works independently from the
others. The methodology used to control this volume is identical to the one
coded in the ROMS model \citep{Marchesiello2001}.


% EXTRAS
\colorbox{yellow}{Explain obc\_vol{\ldots}}

\colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}}

\colorbox{yellow}{OBC rigid lid? {\ldots}}
% ====================================================================
Index: trunk/DOC/TexFiles/Chapters/Chap_LDF.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_LDF.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_LDF.tex (revision 6289)
@@ 69,9 +69,9 @@
r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)}
\;\delta_{i+1/2}[z_t]
 &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t]
+ &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] \ \ \
\\
r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)}
\;\delta_{j+1/2} [z_t]
 &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t]
+ &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] \ \ \
\\
r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2}
@@ 448,4 +448,25 @@
\label{LDF_eiv}
+%%gm from Triad appendix : to be incorporated....
+\gmcomment{
+Values of isoneutral diffusivity and GM coefficient are set as
+described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd},
+N=1,2,3 is set (the default), spatially constant isoneutral $A_l$ and
+GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and
+\np{rn\_aeiv\_0}. If 2Dvarying coefficients are set with
+\key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal
+scale factor according to \eqref{Eq_title} \footnote{Except in global ORCA
+ $0.5^{\circ}$ runs with \key{traldf\_eiv}, where
+ $A_l$ is set like $A_e$ but with a minimum vale of
+ $100\;\mathrm{m}^2\;\mathrm{s}^{1}$}. In idealised setups with
+\key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv}
+is set in the global configurations with \key{traldf\_c2d}, a horizontally varying $A_e$ is
+instead set from the HeldLarichev parameterisation\footnote{In this
+ case, $A_e$ at low latitudes $\theta<20^{\circ}$ is further
+ reduced by a factor $f/f_{20}$, where $f_{20}$ is the value of $f$
+ at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored
+unless it is zero.
+}
+
When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),
an eddy induced tracer advection term is added, the formulation of which
Index: trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_MISC.tex (revision 6140)
+++ 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 nonadjacent 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 timestep. 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^{t1} }{2\rdt} = \ldots \\
\frac{\partial T}{\partial t} &\equiv \frac{T^{t+1}T^{t1}}{2 \widetilde{\rdt}} = \ldots \\
\frac{\partial S}{\partial t} &\equiv \frac{S^{t+1} S^{t1}}{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 highresolution 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 SuccessiveOverRelaxation 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 fivepoint 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_{i1,j} + \; c_{i,j} ^{EW} D_{i,j1} & \\
  \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 blocktridiagonal with
tridiagonal or diagonal blocks. \textbf{A} is a positivedefinite 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 nonlinear 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_{i1,j} +\,a_{i,j}^{W} D_{i,j1}  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}^{t1}
\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_{i1,j} ^{n+1}+\,a_{i,j}^{W} D_{i,j1} ^{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_{i1,j} ^{n}+\,_{i,j}^{W} D_{i,j1} ^{n}  D_{i,j}^n  \tilde{b}_{i,j} \\
R_{i,j}^n = &R_{i,j}^n  \omega \;a_{i,j}^{S}\; R_{i,j1}^n \\
R_{i,j}^n = &R_{i,j}^n  \omega \;a_{i,j}^{W}\; R_{i1,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 socalled redblack 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.$ tripolar 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}^{n1}$. $\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}^{n1} \rangle = 0$ the value of $\beta ^n$ is found:
 $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{r}^{n1}, \textbf{r}^{n1}} \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}^{t1} \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.

% ================================================================












+% ================================================================
+
+
+
+
+
Index: trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex (revision 6289)
@@ 257,5 +257,5 @@
%\newpage
%$\ $\newline % force a new ligne
+%$\ $\newline % force a new line
% ================================================================
@@ 988,8 +988,8 @@
\label{PE_zco_tilde}
The $\tilde{z}$coordinate has been developed by \citet{Leclair_Madec_OM10s}.
+The $\tilde{z}$coordinate has been developed by \citet{Leclair_Madec_OM11}.
It is available in \NEMO since the version 3.4. Nevertheless, it is currently not robust enough
to be used in all possible configurations. Its use is therefore not recommended.
We
+
\newpage
@@ 1113,18 +1113,14 @@
All these parameterisations of subgrid scale physics have advantages and
drawbacks. There are not all available in \NEMO. In the $z$coordinate
formulation, five options are offered for active tracers (temperature and
salinity): second order geopotential operator, second order isoneutral
operator, \citet{Gent1990} parameterisation, fourth order
geopotential operator, and various slightly diffusive advection schemes.
The same options are available for momentum, except
\citet{Gent1990} parameterisation which only involves tracers. In the
$s$coordinate formulation, additional options are offered for tracers: second
order operator acting along $s$surfaces, and for momentum: fourth order
operator acting along $s$surfaces (see \S\ref{LDF}).

\subsubsection{Lateral second order tracer diffusive operator}

The lateral second order tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}):
+drawbacks. There are not all available in \NEMO. For active tracers (temperature and
+salinity) the main ones are: Laplacian and bilaplacian operators acting along
+geopotential or isoneutral surfaces, \citet{Gent1990} parameterisation,
+and various slightly diffusive advection schemes.
+For momentum, the main ones are: Laplacian and bilaplacian operators acting along
+geopotential surfaces, and UBS advection schemes when flux form is chosen for the momentum advection.
+
+\subsubsection{Lateral Laplacian tracer diffusive operator}
+
+The lateral Laplacian tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}):
\begin{equation} \label{Eq_PE_iso_tensor}
D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
@@ 1158,11 +1154,8 @@
but have similar expressions in $z$ and $s$coordinates. In $z$coordinates:
\begin{equation} \label{Eq_PE_iso_slopes}
r_1 =\frac{e_3 }{e_1 } \left( {\frac{\partial \rho }{\partial i}} \right)
 \left( {\frac{\partial \rho }{\partial k}} \right)^{1} \ , \quad
r_1 =\frac{e_3 }{e_1 } \left( {\frac{\partial \rho }{\partial i}} \right)
 \left( {\frac{\partial \rho }{\partial k}} \right)^{1},
\end{equation}
while in $s$coordinates $\partial/\partial k$ is replaced by
$\partial/\partial s$.
+r_1 =\frac{e_3 }{e_1 } \left( \pd[\rho]{i} \right) \left( \pd[\rho]{k} \right)^{1} \, \quad
+r_2 =\frac{e_3 }{e_2 } \left( \pd[\rho]{j} \right) \left( \pd[\rho]{k} \right)^{1} \,
+\end{equation}
+while in $s$coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$.
\subsubsection{Eddy induced velocity}
@@ 1181,5 +1174,5 @@
w^\ast &= \frac{1}{e_1 e_2 }\left[
\frac{\partial }{\partial i}\left( {A^{eiv}\;e_2\,\tilde{r}_1 } \right)
 +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right) \right]
+ +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right) \right]
\end{split}
\end{equation}
@@ 1190,7 +1183,7 @@
\begin{align} \label{Eq_PE_slopes_eiv}
\tilde{r}_n = \begin{cases}
 r_n & \text{in $z$coordinate} \\
+ r_n & \text{in $z$coordinate} \\
r_n + \sigma_n & \text{in \textit{z*} and $s$coordinates}
 \end{cases}
+ \end{cases}
\quad \text{where } n=1,2
\end{align}
@@ 1200,18 +1193,18 @@
to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}).
\subsubsection{Lateral fourth order tracer diffusive operator}

The lateral fourth order tracer diffusive operator is defined by:
+\subsubsection{Lateral bilaplacian tracer diffusive operator}
+
+The lateral bilaplacian tracer diffusive operator is defined by:
\begin{equation} \label{Eq_PE_bilapT}
D^{lT}=\Delta \left( \;\Delta T \right)
+D^{lT}=  \Delta \left( \;\Delta T \right)
\qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right)
\end{equation}
It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with
+It is the Laplacian operator given by \eqref{Eq_PE_iso_tensor} applied twice with
the harmonic eddy diffusion coefficient set to the square root of the biharmonic one.
\subsubsection{Lateral second order momentum diffusive operator}

The second order momentum diffusive operator along $z$ or $s$surfaces is found by
+\subsubsection{Lateral Laplacian momentum diffusive operator}
+
+The Laplacian momentum diffusive operator along $z$ or $s$surfaces is found by
applying \eqref{Eq_PE_lap_vector} to the horizontal velocity vector (see Appendix~\ref{Apdx_B}):
\begin{equation} \label{Eq_PE_lapU}
@@ 1247,8 +1240,9 @@
of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}.
\subsubsection{lateral fourth order momentum diffusive operator}

As for tracers, the fourth order momentum diffusive operator along $z$ or $s$surfaces
is a reentering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU_iso}
with the harmonic eddy diffusion coefficient set to the square root of the biharmonic one.

+\subsubsection{lateral bilaplacian momentum diffusive operator}
+
+As for tracers, the bilaplacian order momentum diffusive operator is a
+reentering Laplacian operator with the harmonic eddy diffusion coefficient
+set to the square root of the biharmonic one. Nevertheless it is currently
+not available in the isoneutral case.
+
Index: trunk/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_SBC.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_SBC.tex (revision 6289)
@@ 1,4 +1,4 @@
% ================================================================
% Chapter � Surface Boundary Condition (SBC, ISF, ICB)
+% Chapter —— Surface Boundary Condition (SBC, ISF, ICB)
% ================================================================
\chapter{Surface Boundary Condition (SBC, ISF, ICB) }
@@ 17,5 +17,6 @@
\item the two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$
\item the incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$
 \item the surface freshwater budget $\left( {\textit{emp},\;\textit{emp}_S } \right)$
+ \item the surface freshwater budget $\left( {\textit{emp}} \right)$
+ \item the surface salt flux associated with freezing/melting of seawater $\left( {\textit{sfx}} \right)$
\end{itemize}
plus an optional field:
@@ 27,14 +28,13 @@
are controlled by namelist \ngn{namsbc} variables: an analytical formulation (\np{ln\_ana}~=~true),
a flux formulation (\np{ln\_flx}~=~true), a bulk formulae formulation (CORE
(\np{ln\_core}~=~true), CLIO (\np{ln\_clio}~=~true) or MFS
+(\np{ln\_blk\_core}~=~true), CLIO (\np{ln\_blk\_clio}~=~true) or MFS
\footnote { Note that MFS bulk formulae compute fluxes only for the ocean component}
(\np{ln\_mfs}~=~true) bulk formulae) and a coupled
formulation (exchanges with a atmospheric model via the OASIS coupler)
(\np{ln\_cpl}~=~true). When used, the atmospheric pressure forces both
ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true).
The frequency at which the six or seven fields have to be updated is the \np{nn\_fsbc}
namelist parameter.
+(\np{ln\_blk\_mfs}~=~true) bulk formulae) and a coupled or mixed forced/coupled formulation
+(exchanges with a atmospheric model via the OASIS coupler) (\np{ln\_cpl} or \np{ln\_mixcpl}~=~true).
+When used ($i.e.$ \np{ln\_apr\_dyn}~=~true), the atmospheric pressure forces both ocean and ice dynamics.
+
+The frequency at which the forcing fields have to be updated is given by the \np{nn\_fsbc} namelist parameter.
When the fields are supplied from data files (flux and bulk formulations), the input fields
need not be supplied on the model grid. Instead a file of coordinates and weights can
+need not be supplied on the model grid. Instead a file of coordinates and weights can
be supplied which maps the data from the supplied grid to the model points
(so called "Interpolation on the Fly", see \S\ref{SBC_iof}).
@@ 42,16 +42,20 @@
can be masked to avoid spurious results in proximity of the coasts as large sealand gradients characterize
most of the atmospheric variables.
+
In addition, the resulting fields can be further modified using several namelist options.
These options control the rotation of vector components supplied relative to an eastnorth
coordinate system onto the local grid directions in the model; the addition of a surface
restoring term to observed SST and/or SSS (\np{ln\_ssr}~=~true); the modification of fluxes
below icecovered areas (using observed icecover or a seaice model)
(\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater
fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of isf melting as lateral inflow (parameterisation)
 or as surface flux at the landice ocean interface (\np{ln\_isf}=~true);
the addition of a freshwater flux adjustment in order to avoid a mean sealevel drift (\np{nn\_fwb}~=~0,~1~or~2); the
transformation of the solar radiation (if provided as daily mean) into a diurnal
cycle (\np{ln\_dm2dc}~=~true); and a neutral drag coefficient can be read from an external wave
model (\np{ln\_cdgw}~=~true). The latter option is possible only in case core or mfs bulk formulas are selected.
+These options control
+\begin{itemize}
+\item the rotation of vector components supplied relative to an eastnorth
+coordinate system onto the local grid directions in the model ;
+\item the addition of a surface restoring term to observed SST and/or SSS (\np{ln\_ssr}~=~true) ;
+\item the modification of fluxes below icecovered areas (using observed icecover or a seaice model) (\np{nn\_ice}~=~0,1, 2 or 3) ;
+\item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln\_rnf}~=~true) ;
+\item the addition of isf melting as lateral inflow (parameterisation) (\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false)
+or as fluxes applied at the landice ocean interface (\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true) ;
+\item the addition of a freshwater flux adjustment in order to avoid a mean sealevel drift (\np{nn\_fwb}~=~0,~1~or~2) ;
+\item the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle (\np{ln\_dm2dc}~=~true) ;
+and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}~=~true).
+\end{itemize}
+The latter option is possible only in case core or mfs bulk formulas are selected.
In this chapter, we first discuss where the surface boundary condition appears in the
@@ 72,69 +76,31 @@
The surface ocean stress is the stress exerted by the wind and the seaice
on the ocean. The two components of stress are assumed to be interpolated
onto the ocean mesh, $i.e.$ resolved onto the model (\textbf{i},\textbf{j}) direction
at $u$ and $v$points They are applied as a surface boundary condition of the
computation of the momentum vertical mixing trend (\mdl{dynzdf} module) :
\begin{equation} \label{Eq_sbc_dynzdf}
\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right_{z=1}
 = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v }
\end{equation}
where $(\tau _u ,\;\tau _v )=(utau,vtau)$ are the two components of the wind
stress vector in the $(\textbf{i},\textbf{j})$ coordinate system.
+on the ocean. It is applied in \mdl{dynzdf} module as a surface boundary condition of the
+computation of the momentum vertical mixing trend (see \eqref{Eq_dynzdf_sbc} in \S\ref{DYN_zdf}).
+As such, it has to be provided as a 2D vector interpolated
+onto the horizontal velocity ocean mesh, $i.e.$ resolved onto the model
+(\textbf{i},\textbf{j}) direction at $u$ and $v$points.
The surface heat flux is decomposed into two parts, a non solar and a solar heat
flux, $Q_{ns}$ and $Q_{sr}$, respectively. The former is the non penetrative part
of the heat flux ($i.e.$ the sum of sensible, latent and long wave heat fluxes).
It is applied as a surface boundary condition trend of the first level temperature
time evolution equation (\mdl{trasbc} module).
\begin{equation} \label{Eq_sbc_trasbc_q}
\frac{\partial T}{\partial t}\equiv \cdots \;+\;\left. {\frac{Q_{ns} }{\rho
_o \;C_p \;e_{3t} }} \right_{k=1} \quad
\end{equation}
$Q_{sr}$ is the penetrative part of the heat flux. It is applied as a 3D
trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=True.

\begin{equation} \label{Eq_sbc_traqsr}
\frac{\partial T}{\partial t}\equiv \cdots \;+\frac{Q_{sr} }{\rho_o C_p \,e_{3t} }\delta _k \left[ {I_w } \right]
\end{equation}
where $I_w$ is a nondimensional function that describes the way the light
penetrates inside the water column. It is generally a sum of decreasing
exponentials (see \S\ref{TRA_qsr}).

The surface freshwater budget is provided by fields: \textit{emp} and $\textit{emp}_S$ which
may or may not be identical. Indeed, a surface freshwater flux has two effects:
it changes the volume of the ocean and it changes the surface concentration of
salt (and other tracers). Therefore it appears in the sea surface height as a volume
flux, \textit{emp} (\textit{dynspg\_xxx} modules), and in the salinity time evolution equations
as a concentration/dilution effect,
$\textit{emp}_{S}$ (\mdl{trasbc} module).
\begin{equation} \label{Eq_trasbc_emp}
\begin{aligned}
&\frac{\partial \eta }{\partial t}\equiv \cdots \;+\;\textit{emp}\quad \\
\\
 &\frac{\partial S}{\partial t}\equiv \cdots \;+\left. {\frac{\textit{emp}_S \;S}{e_{3t} }} \right_{k=1} \\
 \end{aligned}
\end{equation}

In the real ocean, $\textit{emp}=\textit{emp}_S$ and the ocean salt content is conserved,
but it exist several numerical reasons why this equality should be broken.
For example, when the ocean is coupled to a seaice model, the water exchanged between
ice and ocean is slightly salty (mean seaice salinity is $\sim $\textit{4 psu}). In this case,
$\textit{emp}_{S}$ take into account both concentration/dilution effect associated with
freezing/melting and the salt flux between ice and ocean, while \textit{emp} is
only the volume flux. In addition, in the current version of \NEMO, the seaice is
assumed to be above the ocean (the socalled levitating seaice). Freezing/melting does
not change the ocean volume (no impact on \textit{emp}) but it modifies the SSS.
%gm \colorbox{yellow}{(see {\S} on LIM seaice model)}.

Note that SST can also be modified by a freshwater flux. Precipitation (in
particular solid precipitation) may have a temperature significantly different from
the SST. Due to the lack of information about the temperature of
precipitation, we assume it is equal to the SST. Therefore, no
concentration/dilution term appears in the temperature equation. It has to
be emphasised that this absence does not mean that there is no heat flux
associated with precipitation! Precipitation can change the ocean volume and thus the
ocean heat content. It is therefore associated with a heat flux (not yet
diagnosed in the model) \citep{Roullet_Madec_JGR00}).
+of the heat flux ($i.e.$ the sum of sensible, latent and long wave heat fluxes
+plus the heat content of the mass exchange with the atmosphere and seaice).
+It is applied in \mdl{trasbc} module as a surface boundary condition trend of
+the first level temperature time evolution equation (see \eqref{Eq_tra_sbc}
+and \eqref{Eq_tra_sbc_lin} in \S\ref{TRA_sbc}).
+The latter is the penetrative part of the heat flux. It is applied as a 3D
+trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=\textit{true}.
+The way the light penetrates inside the water column is generally a sum of decreasing
+exponentials (see \S\ref{TRA_qsr}).
+
+The surface freshwater budget is provided by the \textit{emp} field.
+It represents the mass flux exchanged with the atmosphere (evaporation minus precipitation)
+and possibly with the seaice and ice shelves (freezing minus melting of ice).
+It affects both the ocean in two different ways:
+$(i)$ it changes the volume of the ocean and therefore appears in the sea surface height
+equation as a volume flux, and
+$(ii)$ it changes the surface temperature and salinity through the heat and salt contents
+of the mass exchanged with the atmosphere, the seaice and the ice shelves.
+
%\colorbox{yellow}{Miss: }
@@ 156,10 +122,15 @@
%
%Explain here all the namlist namsbc variable{\ldots}.
+%
+% explain : use or not of surface currents
%
%\colorbox{yellow}{End Miss }
The ocean model provides the surface currents, temperature and salinity
averaged over \np{nf\_sbc} timestep (\ref{Tab_ssm}).The computation of the
mean is done in \mdl{sbcmod} module.
+The ocean model provides, at each time step, to the surface module (\mdl{sbcmod})
+the surface currents, temperature and salinity.
+These variables are averaged over \np{nf\_sbc} timestep (\ref{Tab_ssm}),
+and it is these averaged fields which are used to computes the surface fluxes
+at a frequency of \np{nf\_sbc} timestep.
+
%TABLE
@@ 458,15 +429,19 @@
%
In some circumstances it may be useful to avoid calculating the 3D temperature, salinity and velocity fields and simply read them in from a previous run.
Options are defined through the \ngn{namsbc\_sas} namelist variables.
+In some circumstances it may be useful to avoid calculating the 3D temperature, salinity and velocity fields
+and simply read them in from a previous run or receive them from OASIS.
For example:
\begin{enumerate}
\item Multiple runs of the model are required in code development to see the affect of different algorithms in
+\begin{itemize}
+\item Multiple runs of the model are required in code development to see the effect of different algorithms in
the bulk formulae.
\item The effect of different parameter sets in the ice model is to be examined.
\end{enumerate}
+\item Development of seaice algorithms or parameterizations.
+\item spinup of the iceberg floats
+\item ocean/seaice simulation with both media running in parallel (\np{ln\_mixcpl}~=~\textit{true})
+\end{itemize}
The StandAlone Surface scheme provides this utility.
+Its options are defined through the \ngn{namsbc\_sas} namelist variables.
A new copy of the model has to be compiled with a configuration based on ORCA2\_SAS\_LIM.
However no namelist parameters need be changed from the settings of the previous run (except perhaps nn{\_}date0)
@@ 474,39 +449,33 @@
Routines replaced are:
\begin{enumerate}
\item \mdl{nemogcm}

 This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (step.F90)
+\begin{itemize}
+\item \mdl{nemogcm} : This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (step.F90)
Since the ocean state is not calculated all associated initialisations have been removed.
\item \mdl{step}

 The main time stepping routine now only needs to call the sbc routine (and a few utility functions).
\item \mdl{sbcmod}

 This has been cut down and now only calculates surface forcing and the ice model required. New surface modules
+\item \mdl{step} : The main time stepping routine now only needs to call the sbc routine (and a few utility functions).
+\item \mdl{sbcmod} : This has been cut down and now only calculates surface forcing and the ice model required. New surface modules
that can function when only the surface level of the ocean state is defined can also be added (e.g. icebergs).
\item \mdl{daymod}

 No ocean restarts are read or written (though the ice model restarts are retained), so calls to restart functions
+\item \mdl{daymod} : No ocean restarts are read or written (though the ice model restarts are retained), so calls to restart functions
have been removed. This also means that the calendar cannot be controlled by time in a restart file, so the user
must make sure that nn{\_}date0 in the model namelist is correct for his or her purposes.
\item \mdl{stpctl}

 Since there is no free surface solver, references to it have been removed from \rou{stp\_ctl} module.
\item \mdl{diawri}

 All 3D data have been removed from the output. The surface temperature, salinity and velocity components (which
+\item \mdl{stpctl} : Since there is no free surface solver, references to it have been removed from \rou{stp\_ctl} module.
+\item \mdl{diawri} : All 3D data have been removed from the output. The surface temperature, salinity and velocity components (which
have been read in) are written along with relevant forcing and ice data.
\end{enumerate}
+\end{itemize}
One new routine has been added:
\begin{enumerate}
\item \mdl{sbcsas}
 This module initialises the input files needed for reading temperature, salinity and velocity arrays at the surface.
+\begin{itemize}
+\item \mdl{sbcsas} : This module initialises the input files needed for reading temperature, salinity and velocity arrays at the surface.
These filenames are supplied in namelist namsbc{\_}sas. Unfortunately because of limitations with the \mdl{iom} module,
the full 3D fields from the mean files have to be read in and interpolated in time, before using just the top level.
Since fldread is used to read in the data, Interpolation on the Fly may be used to change input data resolution.
\end{enumerate}
+\end{itemize}
+
+
+% Missing the description of the 2 following variables:
+% ln_3d_uve = .true. ! specify whether we are supplying a 3D u,v and e3 field
+% ln_read_frq = .false. ! specify whether we must read frq or not
+
+
% ================================================================
@@ 719,7 +688,6 @@
are sent to the atmospheric component.
A generalised coupled interface has been developed. It is currently interfaced with OASIS 3
(\key{oasis3}) and does not support OASIS 4
\footnote{The \key{oasis4} exist. It activates portion of the code that are still under development.}.
+A generalised coupled interface has been developed.
+It is currently interfaced with OASIS3MCT (\key{oasis3}).
It has been successfully used to interface \NEMO to most of the European atmospheric
GCM (ARPEGE, ECHAM, ECMWF, HadAM, HadGAM, LMDz),
@@ 786,12 +754,12 @@
\label{SBC_tide}
A module is available to use the tidal potential forcing and is activated with with \key{tide}.


%nam_tide
+%nam_tide
\namdisplay{nam_tide}
%

Concerning the tidal potential, some parameters are available in namelist \ngn{nam\_tide}:
+%
+
+A module is available to compute the tidal potential and use it in the momentum equation.
+This option is activated when \key{tide} is defined.
+
+Some parameters are available in namelist \ngn{nam\_tide}:
 \np{ln\_tide\_pot} activate the tidal potential forcing
@@ 800,5 +768,4 @@
 \np{clname} is the name of constituent

The tide is generated by the forces of gravity ot the EarthMoon and EarthSun sytem;
@@ 960,6 +927,13 @@
\begin{description}
\item[\np{nn\_isf}~=~1]
The ice shelf cavity is represented. The fwf and heat flux are computed. 2 bulk formulations are available: the ISOMIP one (\np{nn\_isfblk = 1}) described in (\np{nn\_isfblk = 2}), the 3 equation formulation described in \citet{Jenkins1991}. In addition to this,
3 different way to compute the exchange coefficient are available. $\gamma\_{T/S}$ is constant (\np{nn\_gammablk = 0}), $\gamma\_{T/S}$ is velocity dependant \citep{Jenkins2010} (\np{nn\_gammablk = 1}) and $\gamma\_{T/S}$ is velocity dependant and stratification dependent \citep{Holland1999} (\np{nn\_gammablk = 2}). For each of them, the thermal/salt exchange coefficient (\np{rn\_gammat0} and \np{rn\_gammas0}) have to be specified (the default values are for the ISOMIP case).
+The ice shelf cavity is represented. The fwf and heat flux are computed. 2 bulk formulations are available:
+the ISOMIP one (\np{nn\_isfblk = 1}) described in (\np{nn\_isfblk = 2}),
+the 3 equation formulation described in \citet{Jenkins1991}.
+In addition to this, 3 different ways to compute the exchange coefficient are available.
+$\gamma\_{T/S}$ is constant (\np{nn\_gammablk = 0}), $\gamma\_{T/S}$ is velocity dependant
+\citep{Jenkins2010} (\np{nn\_gammablk = 1}) and $\gamma\_{T/S}$ is velocity dependant
+and stratification dependent \citep{Holland1999} (\np{nn\_gammablk = 2}).
+For each of them, the thermal/salt exchange coefficient (\np{rn\_gammat0} and \np{rn\_gammas0})
+have to be specified (the default values are for the ISOMIP case).
Full description, sensitivity and validation in preparation.
@@ 969,5 +943,6 @@
(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}) as in (\np{nn\_isf}~=~3).
Furthermore the fwf is computed using the \citet{Beckmann2003} parameterisation of isf melting.
The effective melting length (\np{sn\_Leff\_isf}) is read from a file and the exchange coefficients are set as (\np{rn\_gammat0}) and (\np{rn\_gammas0}).
+The effective melting length (\np{sn\_Leff\_isf}) is read from a file and the exchange coefficients
+are set as (\np{rn\_gammat0}) and (\np{rn\_gammas0}).
\item[\np{nn\_isf}~=~3]
@@ 1032,5 +1007,5 @@
% Handling of icebergs
% ================================================================
\section{ Handling of icebergs (ICB) }
+\section{Handling of icebergs (ICB)}
\label{ICB_icebergs}
%namberg
@@ 1038,8 +1013,9 @@
%
Icebergs are modelled as lagrangian particles in NEMO.
Their physical behaviour is controlled by equations as described in \citet{Martin_Adcroft_OM10} ).
(Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO.)
Icebergs are initially spawned into one of ten classes which have specific mass and thickness as described in the \ngn{namberg} namelist:
+Icebergs are modelled as lagrangian particles in NEMO \citep{Marsh_GMD2015}.
+Their physical behaviour is controlled by equations as described in \citet{Martin_Adcroft_OM10} ).
+(Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO).
+Icebergs are initially spawned into one of ten classes which have specific mass and thickness as described
+in the \ngn{namberg} namelist:
\np{rn\_initial\_mass} and \np{rn\_initial\_thickness}.
Each class has an associated scaling (\np{rn\_mass\_scaling}), which is an integer representing how many icebergs
@@ 1225,5 +1201,5 @@
The presence at the sea surface of an ice covered area modifies all the fluxes
transmitted to the ocean. There are several way to handle seaice in the system
depending on the value of the \np{nn{\_}ice} namelist parameter.
+depending on the value of the \np{nn\_ice} namelist parameter found in \ngn{namsbc} namelist.
\begin{description}
\item[nn{\_}ice = 0] there will never be seaice in the computational domain.
@@ 1300,27 +1276,41 @@
% 
\subsection [Neutral drag coefficient from external wave model (\textit{sbcwave})]
 {Neutral drag coefficient from external wave model (\mdl{sbcwave})}
+ {Neutral drag coefficient from external wave model (\mdl{sbcwave})}
\label{SBC_wave}
%namwave
\namdisplay{namsbc_wave}
%
\begin{description}

\item [??] In order to read a neutral drag coeff, from an external data source (i.e. a wave model), the
logical variable \np{ln\_cdgw}
 in $namsbc$ namelist must be defined ${.true.}$.
+
+In order to read a neutral drag coeff, from an external data source ($i.e.$ a wave model), the
+logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \textit{true}.
The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the
namelist \ngn{namsbc\_wave} (for external data names, locations, frequency, interpolation and all
the miscellanous options allowed by Input Data generic Interface see \S\ref{SBC_input})
and a 2D field of neutral drag coefficient. Then using the routine
TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided, the drag coefficient is computed according
to stable/unstable conditions of the airsea interface following \citet{Large_Yeager_Rep04}.

\end{description}
+and a 2D field of neutral drag coefficient.
+Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided,
+the drag coefficient is computed according to stable/unstable conditions of the airsea interface following \citet{Large_Yeager_Rep04}.
+
% Griffies doc:
% When running oceanice simulations, we are not explicitly representing land processes, such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift, it is important to balance the hydrological cycle in oceanice models. We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff. The result of the normalization should be a global integrated zero net water input to the oceanice system over a chosen time scale.
%How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step, so that there is always a zero net input of water to the oceanice system. Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used to alter the subsequent year�s water budget in an attempt to damp the annual water imbalance. Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing.
%When running oceanice coupled models, it is incorrect to include the water transport between the ocean and ice models when aiming to balance the hydrological cycle. The reason is that it is the sum of the water in the ocean plus ice that should be balanced when running oceanice models, not the water in any one subcomponent. As an extreme example to illustrate the issue, consider an oceanice model with zero initial sea ice. As the oceanice model spins up, there should be a net accumulation of water in the growing sea ice, and thus a net loss of water from the ocean. The total water contained in the ocean plus ice system is constant, but there is an exchange of water between the subcomponents. This exchange should not be part of the normalization used to balance the hydrological cycle in oceanice models.


+% When running oceanice simulations, we are not explicitly representing land processes,
+% such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift,
+% it is important to balance the hydrological cycle in oceanice models.
+% We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff.
+% The result of the normalization should be a global integrated zero net water input to the oceanice system over
+% a chosen time scale.
+%How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step,
+% so that there is always a zero net input of water to the oceanice system.
+% Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used
+% to alter the subsequent year�s water budget in an attempt to damp the annual water imbalance.
+% Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing.
+% When running oceanice coupled models, it is incorrect to include the water transport between the ocean
+% and ice models when aiming to balance the hydrological cycle.
+% The reason is that it is the sum of the water in the ocean plus ice that should be balanced when running oceanice models,
+% not the water in any one subcomponent. As an extreme example to illustrate the issue,
+% consider an oceanice model with zero initial sea ice. As the oceanice model spins up,
+% there should be a net accumulation of water in the growing sea ice, and thus a net loss of water from the ocean.
+% The total water contained in the ocean plus ice system is constant, but there is an exchange of water between
+% the subcomponents. This exchange should not be part of the normalization used to balance the hydrological cycle
+% in oceanice models.
+
+
Index: trunk/DOC/TexFiles/Chapters/Chap_STO.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_STO.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_STO.tex (revision 6289)
@@ 10,2 +10,9 @@
\newpage
$\ $\newline % force a new line
+%namsbc
+\namdisplay{namsto}
+%
+$\ $\newline % force a new ligne
+
+
+See \cite{Brankart_OM2013} and \cite{Brankart_al_GMD2015} papers for a description of the parameterization.
Index: trunk/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_TRA.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_TRA.tex (revision 6289)
@@ 36,8 +36,8 @@
(BBL) parametrisation, and an internal damping (DMP) term. The terms QSR,
BBC, BBL and DMP are optional. The external forcings and parameterisations
require complex inputs and complex calculations (e.g. bulk formulae, estimation
+require complex inputs and complex calculations ($e.g.$ bulk formulae, estimation
of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and
described in chapters \S\ref{SBC}, \S\ref{LDF} and \S\ref{ZDF}, respectively.
Note that \mdl{tranpc}, the nonpenetrative convection module, although
+Note that \mdl{tranpc}, the nonpenetrative convection module, although
located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields,
is described with the model vertical physics (ZDF) together with other available
@@ 45,5 +45,5 @@
In the present chapter we also describe the diagnostic equations used to compute
the seawater properties (density, BruntVais\"{a}l\"{a} frequency, specific heat and
+the seawater properties (density, BruntV\"{a}is\"{a}l\"{a} frequency, specific heat and
freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}).
@@ 54,6 +54,6 @@
found in the \textit{traTTT} or \textit{traTTT\_xxx} module, in the NEMO/OPA/TRA directory.
The user has the option of extracting each tendency term on the rhs of the tracer
equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{MISC}.
+The user has the option of extracting each tendency term on the RHS of the tracer
+equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{DIA}.
$\ $\newline % force a new ligne
@@ 68,6 +68,7 @@
%
The advection tendency of a tracer in flux form is the divergence of the advective
fluxes. Its discrete expression is given by :
+When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \textit{true}),
+the advection tendency of a tracer is expressed in flux form,
+$i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by :
\begin{equation} \label{Eq_tra_adv}
ADV_\tau =\frac{1}{b_t} \left(
@@ 171,6 +172,6 @@
% 2nd and 4th order centred schemes
% 
\subsection [centred schemes (CEN) (\np{ln\_traadv\_cen})]
 {centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
+\subsection [Centred schemes (CEN) (\np{ln\_traadv\_cen})]
+ {Centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
\label{TRA_adv_cen}
@@ 278,5 +279,5 @@
but on the latter, a splitexplicit time stepping is used, with a number of subtimestep equals
to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited
by vertical advection \citep{Lemarie_OM2015)}. Note that in this case, a similar splitexplicit
+by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar splitexplicit
time stepping should be used on vertical advection of momentum to insure a better stability
(see \S\ref{DYN_zad}).
@@ 405,5 +406,5 @@
Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979}
is used when \np{ln\_traadv\_qck}~=~\textit{true}.
QUICKEST implementation can be found in the \mdl{traadv\_mus} module.
+QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST
@@ 415,6 +416,7 @@
direction where the control of artificial diapycnal fluxes is of paramount importance.
Therefore the vertical flux is evaluated using the CEN2 scheme.
This no longer guarantees the positivity of the scheme. The use of TVD in the vertical
direction (as for the UBS case) should be implemented to restore this property.
+This no longer guarantees the positivity of the scheme.
+The use of FCT in the vertical direction (as for the UBS case) should be implemented
+to restore this property.
%%%gmcomment : Cross term are missing in the current implementation....
@@ 431,41 +433,107 @@
%
Options are defined through the \ngn{namtra\_ldf} namelist variables.
The options available for lateral diffusion are a laplacian (rotated or not)
or a biharmonic operator, the latter being more scaleselective (more
diffusive at small scales). The specification of eddy diffusivity
coefficients (either constant or variable in space and time) as well as the
computation of the slope along which the operators act, are performed in the
\mdl{ldftra} and \mdl{ldfslp} modules, respectively. This is described in Chap.~\ref{LDF}.
+Options are defined through the \ngn{namtra\_ldf} namelist variables.
+They are regrouped in four items, allowing to specify
+$(i)$ the type of operator used (none, laplacian, bilaplacian),
+$(ii)$ the direction along which the operator acts (isolevel, horizontal, isoneutral),
+$(iii)$ some specific options related to the rotated operators ($i.e.$ nonisolevel operator), and
+$(iv)$ the specification of eddy diffusivity coefficient (either constant or variable in space and time).
+Item $(iv)$ will be described in Chap.\ref{LDF} .
+The direction along which the operators act is defined through the slope between this direction and the isolevel surfaces.
+The slope is computed in the \mdl{ldfslp} module and will also be described in Chap.~\ref{LDF}.
+
The lateral diffusion of tracers is evaluated using a forward scheme,
$i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time,
except for the pure vertical component that appears when a rotation tensor
is used. This latter term is solved implicitly together with the
vertical diffusion term (see \S\ref{STP}).

% 
% Isolevel laplacian operator
% 
\subsection [Isolevel laplacian operator (lap) (\np{ln\_traldf\_lap})]
 {Isolevel laplacian operator (lap) (\np{ln\_traldf\_lap}=true) }
\label{TRA_ldf_lap}

A laplacian diffusion operator ($i.e.$ a harmonic operator) acting along the model
surfaces is given by:
+except for the pure vertical component that appears when a rotation tensor is used.
+This latter component is solved implicitly together with the vertical diffusion term (see \S\ref{STP}).
+When \np{ln\_traldf\_msc}~=~\textit{true}, a Method of Stabilizing Correction is used in which
+the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
+
+% 
+% Type of operator
+% 
+\subsection [Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, \np{ln\_traldf\_blp})]
+ {Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, or \np{ln\_traldf\_blp} = true) }
+\label{TRA_ldf_op}
+
+Three operator options are proposed and, one and only one of them must be selected:
+\begin{description}
+\item [\np{ln\_traldf\_NONE}] = true : no operator selected, the lateral diffusive tendency will not be
+applied to the tracer equation. This option can be used when the selected advection scheme
+is diffusive enough (MUSCL scheme for example).
+\item [ \np{ln\_traldf\_lap}] = true : a laplacian operator is selected. This harmonic operator
+takes the following expression: $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $,
+where the gradient operates along the selected direction (see \S\ref{TRA_ldf_dir}),
+and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see Chap.~\ref{LDF}).
+\item [\np{ln\_traldf\_blp}] = true : a bilaplacian operator is selected. This biharmonic operator
+takes the following expression:
+$\mathpzc{B}= \mathpzc{L}\left(\mathpzc{L}(T) \right) = \nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$
+where the gradient operats along the selected direction,
+and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see Chap.~\ref{LDF}).
+In the code, the bilaplacian operator is obtained by calling the laplacian twice.
+\end{description}
+
+Both laplacian and bilaplacian operators ensure the total tracer variance decrease.
+Their primary role is to provide strong dissipation at the smallest scale supported
+by the grid while minimizing the impact on the larger scale features.
+The main difference between the two operators is the scale selectiveness.
+The bilaplacian damping time ($i.e.$ its spin down time) scales like $\lambda^{4}$
+for disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones),
+whereas the laplacian damping time scales only like $\lambda^{2}$.
+
+
+% 
+% Direction of action
+% 
+\subsection [Direction of action (\np{ln\_traldf\_lev}, \np{ln\_traldf\_hor}, \np{ln\_traldf\_iso}, \np{ln\_traldf\_triad})]
+ {Direction of action (\np{ln\_traldf\_lev}, \textit{...\_hor}, \textit{...\_iso}, or \textit{...\_triad} = true) }
+\label{TRA_ldf_dir}
+
+The choice of a direction of action determines the form of operator used.
+The operator is a simple (reentrant) laplacian acting in the (\textbf{i},\textbf{j}) plane
+when isolevel option is used (\np{ln\_traldf\_lev}~=~\textit{true})
+or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}coordinate
+(\np{ln\_traldf\_hor} and \np{ln\_zco} equal \textit{true}).
+The associated code can be found in the \mdl{traldf\_lap\_blp} module.
+The operator is a rotated (reentrant) laplacian when the direction along which it acts
+does not coincide with the isolevel surfaces,
+that is when standard or triad isoneutral option is used (\np{ln\_traldf\_iso} or
+ \np{ln\_traldf\_triad} equals \textit{true}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.),
+or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}coordinate
+(\np{ln\_traldf\_hor} and \np{ln\_sco} equal \textit{true})
+\footnote{In this case, the standard isoneutral operator will be automatically selected}.
+In that case, a rotation is applied to the gradient(s) that appears in the operator
+so that diffusive fluxes acts on the three spatial direction.
+
+The resulting discret form of the three operators (one isolevel and two rotated one)
+is given in the next two subsections.
+
+
+% 
+% isolevel operator
+% 
+\subsection [Isolevel (bi)laplacian operator ( \np{ln\_traldf\_iso})]
+ {Isolevel (bi)laplacian operator ( \np{ln\_traldf\_iso}) }
+\label{TRA_ldf_lev}
+
+The laplacian diffusion operator acting along the model (\textit{i,j})surfaces is given by:
\begin{equation} \label{Eq_tra_ldf_lap}
D_T^{lT} =\frac{1}{b_tT} \left( \;
+D_t^{lT} =\frac{1}{b_t} \left( \;
\delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right]
+ \delta _{j}\left[ A_v^{lT} \; \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right] \;\right)
\end{equation}
where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$cells.
It is implemented in the \mdl{traadv\_lap} module.

This lateral operator is computed in \mdl{traldf\_lap}. It is a \emph{horizontal}
operator ($i.e.$ acting along geopotential surfaces) in the $z$coordinate with
or without partial steps, but is simply an isolevel operator in the $s$coordinate.
It is thus used when, in addition to \np{ln\_traldf\_lap}=true, we have
\np{ln\_traldf\_level}=true or \np{ln\_traldf\_hor}=\np{ln\_zco}=true.
+where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$cells
+and where zero diffusive fluxes is assumed across solid boundaries,
+first (and third in bilaplacian case) horizontal tracer derivative are masked.
+It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module.
+The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap}
+in order to compute the isolevel bilaplacian operator.
+
+It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in the $z$coordinate
+with or without partial steps, but is simply an isolevel operator in the $s$coordinate.
+It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~=~\textit{true},
+we have \np{ln\_traldf\_lev}~=~\textit{true} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~=~\textit{true}.
In both cases, it significantly contributes to diapycnal mixing.
It is therefore not recommended.
+It is therefore never recommended, even when using it in the bilaplacian case.
Note that in the partial step $z$coordinate (\np{ln\_zps}=true), tracers in horizontally
@@ 475,15 +543,19 @@
described in \S\ref{TRA_zpshde}.
% 
% Rotated laplacian operator
% 
\subsection [Rotated laplacian operator (iso) (\np{ln\_traldf\_lap})]
 {Rotated laplacian operator (iso) (\np{ln\_traldf\_lap}=true)}
+
+% 
+% Rotated laplacian operator
+% 
+\subsection [Standard and triad rotated (bi)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad})]
+ {Standard and triad (bi)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad}))}
+\label{TRA_ldf_iso_triad}
+
+%&& Standard rotated (bi)laplacian operator
+%&& 
+\subsubsection [Standard rotated (bi)laplacian operator (\mdl{traldf\_iso})]
+ {Standard rotated (bi)laplacian operator (\mdl{traldf\_iso})}
\label{TRA_ldf_iso}

If the Griffies trad scheme is not employed
(\np{ln\_traldf\_grif}=true; see App.\ref{sec:triad}) the general form of the second order lateral tracer subgrid scale physics
(\ref{Eq_PE_zdf}) takes the following semidiscrete space form in $z$ and
$s$coordinates:
+The general form of the second order lateral tracer subgrid scale physics
+(\ref{Eq_PE_zdf}) takes the following semidiscrete space form in $z$ and $s$coordinates:
\begin{equation} \label{Eq_tra_ldf_iso}
\begin{split}
@@ 527,59 +599,47 @@
of the tracer variance. Nevertheless the treatment performed on the slopes
(see \S\ref{LDF}) allows the model to run safely without any additional
background horizontal diffusion \citep{Guilyardi_al_CD01}. An alternative scheme
developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
+background horizontal diffusion \citep{Guilyardi_al_CD01}.
+
+Note that in the partial step $z$coordinate (\np{ln\_zps}=true), the horizontal derivatives
+at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific treatment.
+They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
+
+%&& Triad rotated (bi)laplacian operator
+%&& 
+\subsubsection [Triad rotated (bi)laplacian operator (\np{ln\_traldf\_triad})]
+ {Triad rotated (bi)laplacian operator (\np{ln\_traldf\_triad})}
+\label{TRA_ldf_triad}
+
+If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}=true ; see App.\ref{sec:triad})
+
+An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of
the algorithm is given in App.\ref{sec:triad}.

Note that in the partial step $z$coordinate (\np{ln\_zps}=true), the horizontal
derivatives at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific
treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.

% 
% Isolevel bilaplacian operator
% 
\subsection [Isolevel bilaplacian operator (bilap) (\np{ln\_traldf\_bilap})]
 {Isolevel bilaplacian operator (bilap) (\np{ln\_traldf\_bilap}=true)}
\label{TRA_ldf_bilap}
The lateral fourth order bilaplacian operator on tracers is obtained by
applying (\ref{Eq_tra_ldf_lap}) twice. The operator requires an additional assumption
on boundary conditions: both first and third derivative terms normal to the
coast are set to zero. It is used when, in addition to \np{ln\_traldf\_bilap}=true,
we have \np{ln\_traldf\_level}=true, or both \np{ln\_traldf\_hor}=true and
\np{ln\_zco}=false. In both cases, it can contribute diapycnal mixing,
although less than in the laplacian case. It is therefore not recommended.

Note that in the code, the bilaplacian routine does not call the laplacian
routine twice but is rather a separate routine that can be found in the
\mdl{traldf\_bilap} module. This is due to the fact that we introduce the
eddy diffusivity coefficient, A, in the operator as:
$\nabla \cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$,
instead of
$\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$
where $a=\sqrt{A}$ and $A<0$. This was a mistake: both formulations
ensure the total variance decrease, but the former requires a larger
number of codelines.

% 
% Rotated bilaplacian operator
% 
\subsection [Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap})]
 {Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap}=true)}
\label{TRA_ldf_bilapg}
+coast are set to zero.
The lateral fourth order operator formulation on tracers is obtained by
applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption
on boundary conditions: first and third derivative terms normal to the
coast, normal to the bottom and normal to the surface are set to zero. It can be found in the
\mdl{traldf\_bilapg}.

It is used when, in addition to \np{ln\_traldf\_bilap}=true, we have
\np{ln\_traldf\_iso}= .true, or both \np{ln\_traldf\_hor}=true and \np{ln\_zco}=true.
This rotated bilaplacian operator has never been seriously
tested. There are no guarantees that it is either free of bugs or correctly formulated.
Moreover, the stability range of such an operator will be probably quite
narrow, requiring a significantly smaller timestep than the one used with an
unrotated operator.
+coast, normal to the bottom and normal to the surface are set to zero.
+
+%&& Option for the rotated operators
+%&& 
+\subsubsection [Option for the rotated operators]
+ {Option for the rotated operators}
+\label{TRA_ldf_options}
+
+\np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
+
+\np{rn\_slpmax} = slope limit (both operators)
+
+\np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only)
+
+\np{rn\_sw\_triad} =1 switching triad ; =0 all 4 triads used (triad only)
+
+\np{ln\_botmix\_triad} = lateral mixing on bottom (triad only)
% ================================================================
@@ 593,5 +653,5 @@
%
Options are defined through the \ngn{namzdf} namelist variables.
+Options are defined through the \ngn{namzdf} namelist variables.
The formulation of the vertical subgrid scale tracer physics is the same
for all the vertical coordinates, and is based on a laplacian operator.
@@ 651,16 +711,10 @@
the thickness of the top model layer.
Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components ($i.e.$ atmosphere, seaice, land),
the change in the heat and salt content of the surface layer of the ocean is due both
to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)
 and to the heat and salt content of the mass exchange.
\sgacomment{ the following does not apply to the release to which this documentation is
attached and so should not be included ....
In a forthcoming release, these two parts, computed in the surface module (SBC), will be included directly
in $Q_{ns}$, the surface heat flux and $F_{salt}$, the surface salt flux.
The specification of these fluxes is further detailed in the SBC chapter (see \S\ref{SBC}).
This change will provide a forcing formulation which is the same for any tracer (including temperature and salinity).

In the current version, the situation is a little bit more complicated. }
+Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components
+($i.e.$ atmosphere, seaice, land), the change in the heat and salt content of the surface layer
+of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)
+and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$,
+the surface heat flux, and $F_{salt}$, the surface salt flux (see \S\ref{SBC} for further details).
+By doing this, the forcing formulation is the same for any tracer (including temperature and salinity).
The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following
@@ 669,22 +723,33 @@
$\bullet$ $Q_{ns}$, the nonsolar part of the net surface heat flux that crosses the sea surface
(i.e. the difference between the total surface heat flux and the fraction of the short wave flux that
penetrates into the water column, see \S\ref{TRA_qsr})

$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)

$\bullet$ $\textit{emp}_S$, an equivalent mass flux taking into account the effect of iceocean mass exchange

$\bullet$ \textit{rnf}, the mass flux associated with runoff (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)

The $\textit{emp}_S$ field is not simply the budget of evaporationprecipitation+freezingmelting because
the seaice is not currently embedded in the ocean but levitates above it. There is no mass
exchanged between the seaice and the ocean. Instead we only take into account the salt
flux associated with the nonzero salinity of seaice, and the concentration/dilution effect
due to the freezing/melting (F/M) process. These two parts of the forcing are then converted into
an equivalent mass flux given by $\textit{emp}_S  \textit{emp}$. As a result of this mess,
the surface boundary condition on temperature and salinity is applied as follows:

In the nonlinear free surface case (\key{vvl} is defined):
+penetrates into the water column, see \S\ref{TRA_qsr}) plus the heat content associated with
+of the mass exchange with the atmosphere and lands.
+
+$\bullet$ $\textit{sfx}$, the salt flux resulting from iceocean mass exchange (freezing, melting, ridging...)
+
+$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)
+ and possibly with the seaice and iceshelves.
+
+$\bullet$ \textit{rnf}, the mass flux associated with runoff
+(see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
+
+The surface boundary condition on temperature and salinity is applied as follows:
\begin{equation} \label{Eq_tra_sbc}
+\begin{aligned}
+ &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right_{k=1} } &\overline{ Q_{ns} }^t & \\
+& F^S =\frac{ 1 }{\rho _o \, \left. e_{3t} \right_{k=1} } &\overline{ \textit{sfx} }^t & \\
+ \end{aligned}
+\end{equation}
+where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps
+($t\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the
+divergence of odd and even time step (see \S\ref{STP}).
+
+In the linear free surface case (\np{ln\_linssh}~=~\textit{true}),
+an additional term has to be added on both temperature and salinity.
+On temperature, this term remove the heat content associated with mass exchange
+that has been added to $Q_{ns}$. On salinity, this term mimics the concentration/dilution effect that
+would have resulted from a change in the volume of the first level.
+The resulting surface boundary condition is applied as follows:
+\begin{equation} \label{Eq_tra_sbc_lin}
\begin{aligned}
&F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right_{k=1} }
@@ 692,50 +757,11 @@
%
& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right_{k=1} }
 &\overline{ \left( (\textit{emp}_S  \textit{emp})\;\left. S \right_{k=1} \right) }^t & \\
+ &\overline{ \left( \;\textit{sfx}  \textit{emp} \;\left. S \right_{k=1} \right) }^t & \\
\end{aligned}
\end{equation}

In the linear free surface case (\key{vvl} not defined):
\begin{equation} \label{Eq_tra_sbc_lin}
\begin{aligned}
 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right_{k=1} } &\overline{ Q_{ns} }^t & \\
%
& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right_{k=1} }
 &\overline{ \left( \textit{emp}_S\;\left. S \right_{k=1} \right) }^t & \\
 \end{aligned}
\end{equation}
where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps
($t\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the
divergence of odd and even time step (see \S\ref{STP}).

The two set of equations, \eqref{Eq_tra_sbc} and \eqref{Eq_tra_sbc_lin}, are obtained
by assuming that the temperature of precipitation and evaporation are equal to
the ocean surface temperature and that their salinity is zero. Therefore, the heat content
of the \textit{emp} budget must be added to the temperature equation in the variable volume case,
while it does not appear in the constant volume case. Similarly, the \textit{emp} budget affects
the ocean surface salinity in the constant volume case (through the concentration dilution effect)
while it does not appears explicitly in the variable volume case since salinity change will be
induced by volume change. In both constant and variable volume cases, surface salinity
will change with iceocean salt flux and F/M flux (both contained in $\textit{emp}_S  \textit{emp}$) without mass exchanges.

Note that the concentration/dilution effect due to F/M is computed using
a constant ice salinity as well as a constant ocean salinity.
This approximation suppresses the correlation between \textit{SSS}
and F/M flux, allowing the iceocean salt exchanges to be conservative.
Indeed, if this approximation is not made, even if the F/M budget is zero
on average over the whole ocean domain and over the seasonal cycle,
the associated salt flux is not zero, since seasurface salinity and F/M flux are
intrinsically correlated (high \textit{SSS} are found where freezing is
strong whilst low \textit{SSS} is usually associated with high melting areas).

Even using this approximation, an exact conservation of heat and salt content
is only achieved in the variable volume case. In the constant volume case,
there is a small imbalance associated with the product $(\partial_t\eta  \textit{emp}) * \textit{SSS}$.
Nevertheless, the salt content variation is quite small and will not induce
a long term drift as there is no physical reason for $(\partial_t\eta  \textit{emp})$
and \textit{SSS} to be correlated \citep{Roullet_Madec_JGR00}.
Note that, while quite small, the imbalance in the constant volume case is larger
+Note that an exact conservation of heat and salt content is only achieved with nonlinear free surface.
+In the linear free surface case, there is a small imbalance. The imbalance is larger
than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.
This is the reason why the modified filter is not applied in the constant volume case.
+This is the reason why the modified filter is not applied in the linear free surface case (see \S\ref{STP}).
% 
@@ 849,5 +875,5 @@
\label{TRA_bbc}
%nambbc
\namdisplay{namtra_bbc}
+\namdisplay{nambbc}
%
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ 1093,15 +1119,36 @@
\subsection[DMP\_TOOLS]{Generating resto.nc using DMP\_TOOLS}
DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input. This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
+DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$.
+Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled
+and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input.
+This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1.
+The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work.
+The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
%nam_dmp_create
\namdisplay{nam_dmp_create}
+\namtools{namelist_dmp}
%
\np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in the \nl{nam\_zoom\_dmp} name list.
The remaining switch namelist variables determine the spatial variation of the restoration coefficient in nonzoom configurations. \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for the ORCA4, ORCA2 and ORCA05 configurations. If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference configurations with previous model versions. \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. This option only has an effect if \np{ln\_full\_field} is true. \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. Finally \np{ln\_custom} specifies that the custom module will be called. This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.

The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}. Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to the full values of a 10$^{\circ}$ latitud band. This is often used because of the short adjustment time scale in the equatorial region \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.
+The remaining switch namelist variables determine the spatial variation of the restoration coefficient in nonzoom configurations.
+\np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain.
+\np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea
+for the ORCA4, ORCA2 and ORCA05 configurations.
+If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as
+a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference
+configurations with previous model versions.
+\np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines.
+This option only has an effect if \np{ln\_full\_field} is true.
+\np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer.
+Finally \np{ln\_custom} specifies that the custom module will be called.
+This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.
+
+The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}.
+Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to
+the full values of a 10$^{\circ}$ latitud band.
+This is often used because of the short adjustment time scale in the equatorial region
+\citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a
+hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.
% ================================================================
@@ 1271,10 +1318,10 @@
% 
% BruntVais\"{a}l\"{a} Frequency
% 
\subsection{BruntVais\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
+% BruntV\"{a}is\"{a}l\"{a} Frequency
+% 
+\subsection{BruntV\"{a}is\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
\label{TRA_bn2}
An accurate computation of the ocean stability (i.e. of $N$, the bruntVais\"{a}l\"{a}
+An accurate computation of the ocean stability (i.e. of $N$, the bruntV\"{a}is\"{a}l\"{a}
frequency) is of paramount importance as determine the ocean stratification and
is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent
@@ 1332,5 +1379,6 @@
\label{TRA_zpshde}
\gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, I've changed "derivative" to "difference" and "mean" to "average"}
+\gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators,
+ I've changed "derivative" to "difference" and "mean" to "average"}
With partial bottom cells (\np{ln\_zps}=true), in general, tracers in horizontally
Index: trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex (revision 6289)
@@ 34,7 +34,7 @@
coefficients can be assumed to be either constant, or a function of the local
Richardson number, or computed from a turbulent closure model (either
TKE or KPP formulation). The computation of these coefficients is initialized
+TKE or GLS formulation). The computation of these coefficients is initialized
in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or
\mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer
+\mdl{zdfgls} modules. The trends due to the vertical momentum and tracer
diffusion, including the surface forcing, are computed and added to the
general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.
@@ 355,9 +355,46 @@
%%
To be add here a description of "penetration of TKE" and the associated namelist parameters
 \np{nn\_etau}, \np{rn\_efr} and \np{nn\_htau}.
+Vertical mixing parameterizations commonly used in ocean general circulation models
+tend to produce mixedlayer depths that are too shallow during summer months and windy conditions.
+This bias is particularly acute over the Southern Ocean.
+To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{Rodgers_2014}.
+The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations,
+but rather is meant to account for observed processes that affect the density structure of
+the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme
+($i.e.$ nearinertial oscillations and ocean swells and waves).
+
+When using this parameterization ($i.e.$ when \np{nn\_etau}~=~1), the TKE input to the ocean ($S$)
+imposed by the winds in the form of nearinertial oscillations, swell and waves is parameterized
+by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by:
+\begin{equation} \label{ZDF_Ehtau}
+S = (1f_i) \; f_r \; e_s \; e^{z / h_\tau}
+\end{equation}
+where
+$z$ is the depth,
+$e_s$ is TKE surface boundary condition,
+$f_r$ is the fraction of the surface TKE that penetrate in the ocean,
+$h_\tau$ is a vertical mixing length scale that controls exponential shape of the penetration,
+and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely
+covered by seaice).
+The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter.
+The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}~=~0)
+or a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m
+at high latitudes (\np{nn\_etau}~=~1).
+
+Note that two other option existe, \np{nn\_etau}~=~2, or 3. They correspond to applying
+\eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part
+of the stress to evaluate the fraction of TKE that penetrate the ocean.
+Those two options are obsolescent features introduced for test purposes.
+They will be removed in the next release.
+
+
% from Burchard et al OM 2008 :
% the most critical process not reproduced by statistical turbulence models is the activity of internal waves and their interaction with turbulence. After the Reynolds decomposition, internal waves are in principle included in the RANS equations, but later partially excluded by the hydrostatic assumption and the model resolution. Thus far, the representation of internal wave mixing in ocean models has been relatively crude (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
+% the most critical process not reproduced by statistical turbulence models is the activity of
+% internal waves and their interaction with turbulence. After the Reynolds decomposition,
+% internal waves are in principle included in the RANS equations, but later partially
+% excluded by the hydrostatic assumption and the model resolution.
+% Thus far, the representation of internal wave mixing in ocean models has been relatively crude
+% (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
@@ 573,19 +610,4 @@
Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}.
% 
% K Profile Parametrisation (KPP)
% 
\subsection{K Profile Parametrisation (KPP) (\key{zdfkpp}) }
\label{ZDF_kpp}

%namkpp
\namdisplay{namzdf_kpp}
%

The KKP scheme has been implemented by J. Chanut ...
Options are defined through the \ngn{namzdf\_kpp} namelist variables.

\colorbox{yellow}{Add a description of KPP here.}

% ================================================================
@@ 636,5 +658,5 @@
Options are defined through the \ngn{namzdf} namelist variables.
The nonpenetrative convective adjustment is used when \np{ln\_zdfnpc}=true.
+The nonpenetrative convective adjustment is used when \np{ln\_zdfnpc}~=~\textit{true}.
It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously
the statically unstable portion of the water column, but only until the density
@@ 644,5 +666,5 @@
(Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is
found. Assume in the following that the instability is located between levels
$k$ and $k+1$. The potential temperature and salinity in the two levels are
+$k$ and $k+1$. The temperature and salinity in the two levels are
vertically mixed, conserving the heat and salt contents of the water column.
The new density is then computed by a linear approximation. If the new
@@ 664,24 +686,14 @@
\citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.
Note that in the current implementation of this algorithm presents several
limitations. First, potential density referenced to the sea surface is used to
check whether the density profile is stable or not. This is a strong
simplification which leads to large errors for realistic ocean simulations.
Indeed, many water masses of the world ocean, especially Antarctic Bottom
Water, are unstable when represented in surfacereferenced potential density.
The scheme will erroneously mix them up. Second, the mixing of potential
density is assumed to be linear. This assures the convergence of the algorithm
even when the equation of state is nonlinear. Small static instabilities can thus
persist due to cabbeling: they will be treated at the next time step.
Third, temperature and salinity, and thus density, are mixed, but the
corresponding velocity fields remain unchanged. When using a Richardson
Number dependent eddy viscosity, the mixing of momentum is done through
the vertical diffusion: after a static adjustment, the Richardson Number is zero
and thus the eddy viscosity coefficient is at a maximum. When this convective
adjustment algorithm is used with constant vertical eddy viscosity, spurious
solutions can occur since the vertical momentum diffusion remains small even
after a static adjustment. In that case, we recommend the addition of momentum
mixing in a manner that mimics the mixing in temperature and salinity
\citep{Speich_PhD92, Speich_al_JPO96}.
+The current implementation has been modified in order to deal with any non linear
+equation of seawater (L. Brodeau, personnal communication).
+Two main differences have been introduced compared to the original algorithm:
+$(i)$ the stability is now checked using the BruntV\"{a}is\"{a}l\"{a} frequency
+(not the the difference in potential density) ;
+$(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients
+are vertically mixed in the same way their temperature and salinity has been mixed.
+These two modifications allow the algorithm to perform properly and accurately
+with TEOS10 or EOS80 without having to recompute the expansion coefficients at each
+mixing iteration.
% 
@@ 689,5 +701,5 @@
% 
\subsection [Enhanced Vertical Diffusion (\np{ln\_zdfevd})]
 {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=true)}
+ {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=true)}
\label{ZDF_evd}
Index: trunk/DOC/TexFiles/Chapters/Introduction.tex
===================================================================
 trunk/DOC/TexFiles/Chapters/Introduction.tex (revision 6140)
+++ trunk/DOC/TexFiles/Chapters/Introduction.tex (revision 6289)
@@ 32,5 +32,5 @@
The equations are written in a curvilinear coordinate system, with a choice of vertical
coordinates ($z$, $s$, \textit{z*}, \textit{s*}, $\tilde{z}$, $\tilde{s}$, and a mixture of them).
Momentum equations are formulated in the vector invariant form or in the flux form.
+Momentum equations are formulated in vector invariant or flux form.
Dimensional units in the meter, kilogram, second (MKS) international system
are used throughout.
Index: trunk/DOC/TexFiles/Namelist/nam_asminc
===================================================================
 trunk/DOC/TexFiles/Namelist/nam_asminc (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/nam_asminc (revision 6289)
@@ 0,0 +1,18 @@
+!
+&nam_asminc ! assimilation increments ('key_asminc')
+!
+ ln_bkgwri = .false. ! Logical switch for writing out background state
+ ln_trainc = .false. ! Logical switch for applying tracer increments
+ ln_dyninc = .false. ! Logical switch for applying velocity increments
+ ln_sshinc = .false. ! Logical switch for applying SSH increments
+ ln_asmdin = .false. ! Logical switch for Direct Initialization (DI)
+ ln_asmiau = .false. ! Logical switch for Incremental Analysis Updating (IAU)
+ nitbkg = 0 ! Timestep of background in [0,nitendnit0001]
+ nitdin = 0 ! Timestep of background for DI in [0,nitendnit0001]
+ nitiaustr = 1 ! Timestep of start of IAU interval in [0,nitendnit0001]
+ nitiaufin = 15 ! Timestep of end of IAU interval in [0,nitendnit0001]
+ niaufn = 0 ! Type of IAU weighting function
+ ln_salfix = .false. ! Logical switch for ensuring that the sa > salfixmin
+ salfixmin = 9999 ! Minimum salinity after applying the increments
+ nn_divdmp = 0 ! Number of iterations of divergence damping operator
+/
Index: trunk/DOC/TexFiles/Namelist/nam_dia25h
===================================================================
 trunk/DOC/TexFiles/Namelist/nam_dia25h (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/nam_dia25h (revision 6289)
@@ 1,5 +1,5 @@
!
&nam_dia25h ! 25h Mean Output
+&nam_dia25h ! 25h Mean Output
!
 ln_dia25h = .true. ! Choose 25 mean output or not
+ ln_dia25h = .false. ! Choose 25h mean output or not
/
Index: trunk/DOC/TexFiles/Namelist/nam_diaharm
===================================================================
 trunk/DOC/TexFiles/Namelist/nam_diaharm (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/nam_diaharm (revision 6289)
@@ 0,0 +1,9 @@
+!
+&nam_diaharm ! Harmonic analysis of tidal constituents ('key_diaharm')
+!
+ nit000_han = 1 ! First time step used for harmonic analysis
+ nitend_han = 75 ! Last time step used for harmonic analysis
+ nstep_han = 15 ! Time step frequency for harmonic analysis
+ tname(1) = 'M2' ! Name of tidal constituents
+ tname(2) = 'K1'
+/
Index: trunk/DOC/TexFiles/Namelist/nam_diatmb
===================================================================
 trunk/DOC/TexFiles/Namelist/nam_diatmb (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/nam_diatmb (revision 6289)
@@ 1,5 +1,5 @@
!
&nam_diatmb ! 25h Mean Output
+&nam_diatmb ! Top Middle Bottom Output
!
 ln_diatmb = .true. ! Choose Top Middle and Bottom output or not
+ ln_diatmb = .false. ! Choose Top Middle and Bottom output or not
/
Index: trunk/DOC/TexFiles/Namelist/nam_dmp_create
===================================================================
 trunk/DOC/TexFiles/Namelist/nam_dmp_create (revision 6140)
+++ (revision )
@@ 1,23 +1,0 @@
&nam_dmp_create
 cp_cfg = 'orca' ! Name of model grid (orca and C1D have special options
 ! otherwise ignored)
 cp_cfz = 'antarctic' ! Name of zoom configuration (arctic and antarctic have
 ! some special treatment if lzoom=.true.)
 jp_cfg = 2 ! Resolution of the model (used for med_red_seas damping)
 lzoom = .false. ! Zoom configuration or not
 ln_full_field = .false. ! Calculate coefficient over whole of domain
 ln_med_red_seas = .true. ! Damping in Med/Red Seas
 ! (or local modifications here if ln_full_field=.true.)
 ln_old_31_lev_code = .true. ! Replicate behaviour of old online code for 31 level model
 ! (Med/Red seas damping based on level number instead of depth)
 ln_coast = .true. ! Reduce near to coastlines
 ln_zero_top_layer = .true. ! No damping in top layer
 ln_custom = .false. ! Call "custom" module to apply user modifications to the
 ! damping coefficient field
 pn_surf = 0.25 ! Surface Relaxation timescale (days)
 pn_bot = 0.25 ! Bottom relaxation timescale (days)
 pn_dep = 1000 ! Transition depth from upper to deep ocean
 nn_hdmp = 10 ! Damp poleward of this latitude (smooth transition up to maximum damping)
 jperio = 2 ! Lateral boundary condition (as specified in namelist_cfg for model run).
/

Index: trunk/DOC/TexFiles/Namelist/nam_vvl
===================================================================
 trunk/DOC/TexFiles/Namelist/nam_vvl (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/nam_vvl (revision 6289)
@@ 0,0 +1,14 @@
+!
+&nam_vvl ! vertical coordinate options (default: zstar)
+!
+ ln_vvl_zstar = .true. ! zstar vertical coordinate
+ ln_vvl_ztilde = .false. ! ztilde vertical coordinate: only high frequency variations
+ ln_vvl_layer = .false. ! full layer vertical coordinate
+ ln_vvl_ztilde_as_zstar = .false. ! ztilde vertical coordinate emulating zstar
+ ln_vvl_zstar_at_eqtor = .false. ! ztilde near the equator
+ rn_ahe3 = 0.0e0 ! thickness diffusion coefficient
+ rn_rst_e3t = 30.e0 ! ztilde to zstar restoration timescale [days]
+ rn_lf_cutoff = 5.0e0 ! cutoff frequency for lowpass filter [days]
+ rn_zdef_max = 0.9e0 ! maximum fractional e3t deformation
+ ln_vvl_dbg = .true. ! debug prints (T/F)
+/
Index: trunk/DOC/TexFiles/Namelist/namasm
===================================================================
 trunk/DOC/TexFiles/Namelist/namasm (revision 6140)
+++ (revision )
@@ 1,18 +1,0 @@
!
&nam_asminc ! assimilation increments ('key_asminc')
!
 ln_bkgwri = .false. ! Logical switch for writing out background state
 ln_trainc = .false. ! Logical switch for applying tracer increments
 ln_dyninc = .false. ! Logical switch for applying velocity increments
 ln_sshinc = .false. ! Logical switch for applying SSH increments
 ln_asmdin = .false. ! Logical switch for Direct Initialization (DI)
 ln_asmiau = .false. ! Logical switch for Incremental Analysis Updating (IAU)
 nitbkg = 0 ! Timestep of background in [0,nitendnit0001]
 nitdin = 0 ! Timestep of background for DI in [0,nitendnit0001]
 nitiaustr = 1 ! Timestep of start of IAU interval in [0,nitendnit0001]
 nitiaufin = 15 ! Timestep of end of IAU interval in [0,nitendnit0001]
 niaufn = 0 ! Type of IAU weighting function
 ln_salfix = .false. ! Logical switch for ensuring that the sa > salfixmin
 salfixmin = 9999 ! Minimum salinity after applying the increments
 nn_divdmp = 0 ! Number of iterations of divergence damping operator
/
Index: trunk/DOC/TexFiles/Namelist/nambbc
===================================================================
 trunk/DOC/TexFiles/Namelist/nambbc (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/nambbc (revision 6289)
@@ 0,0 +1,14 @@
+!
+&nambbc ! bottom temperature boundary condition (default: NO)
+!
+! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask !
+! ! ! (if <0 months) ! name ! (logical) ! (T/F ) ! 'monthly' ! filename ! pairing ! filename !
+ sn_qgh ='geothermal_heating.nc', 12. , 'heatflow', .false. , .true. , 'yearly' , '' , '' , ''
+ !
+ ln_trabbc = .false. ! Apply a geothermal heating at the ocean bottom
+ nn_geoflx = 2 ! geothermal heat flux: = 0 no flux
+ ! = 1 constant flux
+ ! = 2 variable flux (read in geothermal_heating.nc in mW/m2)
+ rn_geoflx_cst = 86.4e3 ! Constant value of geothermal heat flux [W/m2]
+ cn_dir = './' ! root directory for the location of the runoff files
+/
Index: trunk/DOC/TexFiles/Namelist/nambdy_dta2
===================================================================
 trunk/DOC/TexFiles/Namelist/nambdy_dta2 (revision 6140)
+++ (revision )
@@ 1,10 +1,0 @@
!
&nambdy_dta ! open boundaries  external data ("key_bdy")
!
! ! file name ! frequency (hours) ! variable ! time interpol. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask !
! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename !
 bn_tem = 'amm12_bdyT_clim' , 1 , 'votemper' , .true. , .true. , 'daily' , '' , '' , ''
 bn_sal = 'amm12_bdyT_clim' , 1 , 'vosaline' , .true. , .true. , 'daily' , '' , '' , ''
 cn_dir = 'bdydta/'
 ln_full_vel = .false.
/
Index: trunk/DOC/TexFiles/Namelist/nambdy_index
===================================================================
 trunk/DOC/TexFiles/Namelist/nambdy_index (revision 6140)
+++ (revision )
@@ 1,17 +1,0 @@
!
&nambdy_index ! open boundaries  definition ("key_bdy")
!
 nbdysege = 0
 nbdysegw = 1
 jpiwob = 2
 jpjwdt = 2
 jpjwft = 191
 nbdysegn = 1
 jpjnob = 191
 jpindt = 2
 jpinft = 146
 nbdysegs = 1
 jpjsob = 2
 jpisdt = 2
 jpisft = 43
/
Index: trunk/DOC/TexFiles/Namelist/nambdy_tide
===================================================================
 trunk/DOC/TexFiles/Namelist/nambdy_tide (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/nambdy_tide (revision 6289)
@@ 3,5 +3,5 @@
!
filtide = 'bdydta/amm12_bdytide_' ! file name root of tidal forcing files
 ln_bdytide_2ddta = .false. !
 ln_bdytide_conj = .false. !
+ ln_bdytide_2ddta = .false. !
+ ln_bdytide_conj = .false. !
/
Index: trunk/DOC/TexFiles/Namelist/namberg
===================================================================
 trunk/DOC/TexFiles/Namelist/namberg (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namberg (revision 6289)
@@ 27,7 +27,7 @@
rn_speed_limit = 0. ! CFL speed limit for a berg
! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask !
! ! ! (if <0 months) ! name ! (logical) ! (T/F ) ! 'monthly' ! filename ! pairing ! filename !
 sn_icb = 'calving' , 1 , 'calvingmask', .true. , .true. , 'yearly' , '' , '' , ''
+! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask !
+! ! ! (if <0 months) ! name ! (logical) ! (T/F ) ! 'monthly' ! filename ! pairing ! filename !
+ sn_icb = 'calving', 1 , 'calvingmask', .true. , .true. , 'yearly' , '' , '' , ''
cn_dir = './'
Index: trunk/DOC/TexFiles/Namelist/namc1d
===================================================================
 trunk/DOC/TexFiles/Namelist/namc1d (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namc1d (revision 6289)
@@ 0,0 +1,7 @@
+!
+&namc1d ! 1D configuration options ("key_c1d")
+!
+ rn_lat1d = 50 ! Column latitude (default at PAPA station)
+ rn_lon1d = 145 ! Column longitude (default at PAPA station)
+ ln_c1d_locpt= .true. ! Localization of 1D config in a grid (T) or independant point (F)
+/
Index: trunk/DOC/TexFiles/Namelist/namc1d_dyndmp
===================================================================
 trunk/DOC/TexFiles/Namelist/namc1d_dyndmp (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namc1d_dyndmp (revision 6289)
@@ 0,0 +1,5 @@
+!
+&namc1d_dyndmp ! U & V newtonian damping ("key_c1d")
+!
+ ln_dyndmp = .false. ! add a damping term (T) or not (F)
+/
Index: trunk/DOC/TexFiles/Namelist/namc1d_uvd
===================================================================
 trunk/DOC/TexFiles/Namelist/namc1d_uvd (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namc1d_uvd (revision 6289)
@@ 0,0 +1,12 @@
+!
+&namc1d_uvd ! data: U & V currents ("key_c1d")
+!
+! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask !
+! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename !
+ sn_ucur = 'ucurrent' , 1 ,'u_current', .false. , .true. , 'monthly' , '' , 'Ume' , ''
+ sn_vcur = 'vcurrent' , 1 ,'v_current', .false. , .true. , 'monthly' , '' , 'Vme' , ''
+!
+ cn_dir = './' ! root directory for the location of the files
+ ln_uvd_init = .false. ! Initialisation of ocean U & V with U & V input data (T) or not (F)
+ ln_uvd_dyndmp = .false. ! damping of ocean U & V toward U & V input data (T) or not (F)
+/
Index: trunk/DOC/TexFiles/Namelist/namcfg
===================================================================
 trunk/DOC/TexFiles/Namelist/namcfg (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namcfg (revision 6289)
@@ 0,0 +1,22 @@
+!
+&namcfg ! parameters of the configuration
+!
+ cp_cfg = "default" ! name of the configuration
+ cp_cfz = "no zoom" ! name of the zoom of configuration
+ jp_cfg = 0 ! resolution of the configuration
+ jpidta = 10 ! 1st lateral dimension ( >= jpi )
+ jpjdta = 12 ! 2nd " " ( >= jpj )
+ jpkdta = 31 ! number of levels ( >= jpk )
+ jpiglo = 10 ! 1st dimension of global domain > i =jpidta
+ jpjglo = 12 ! 2nd   > j =jpjdta
+ jpizoom = 1 ! left bottom (i,j) indices of the zoom
+ jpjzoom = 1 ! in data domain indices
+ jperio = 0 ! lateral cond. type (between 0 and 6)
+ ! = 0 closed ; = 1 cyclic EastWest
+ ! = 2 equatorial symmetric ; = 3 North fold Tpoint pivot
+ ! = 4 cyclic EastWest AND North fold Tpoint pivot
+ ! = 5 North fold Fpoint pivot
+ ! = 6 cyclic EastWest AND North fold Fpoint pivot
+ ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present
+ ! in netcdf input files, as the start jrow for reading
+/
Index: trunk/DOC/TexFiles/Namelist/namcfg_orca1
===================================================================
 trunk/DOC/TexFiles/Namelist/namcfg_orca1 (revision 6140)
+++ (revision )
@@ 1,12 +1,0 @@
!
&namcfg ! parameters of the configuration
!
 cp_cfg = "orca" ! name of the configuration
 jp_cfg = 1 ! resolution of the configuration
 jpidta = 362 ! 1st lateral dimension ( >= jpi )
 jpjdta = 292 ! 2nd " " ( >= jpj )
 jpkdta = 75 ! number of levels ( >= jpk )
 jpiglo = 362 ! 1st dimension of global domain > i =jpidta
 jpjglo = 292 ! 2nd   > j =jpjdta
 jperio = 6 ! lateral cond. type (between 0 and 6)
 ln_use_jattr = .true. ! use (T) the file attribute: open_ocean_jstart if present
Index: trunk/DOC/TexFiles/Namelist/namcla
===================================================================
 trunk/DOC/TexFiles/Namelist/namcla (revision 6140)
+++ (revision )
@@ 1,5 +1,0 @@
!
&namcla ! cross land advection
!
 nn_cla = 0 ! advection between 2 ocean pts separates by land
/
Index: trunk/DOC/TexFiles/Namelist/namcrs
===================================================================
 trunk/DOC/TexFiles/Namelist/namcrs (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namcrs (revision 6289)
@@ 0,0 +1,15 @@
+!
+&namcrs ! coarsened grid (for outputs and/or TOP) ("key_crs")
+!
+ nn_factx = 3 ! Reduction factor of xdirection
+ nn_facty = 3 ! Reduction factor of ydirection
+ nn_binref = 0 ! Bin centering preference: NORTH or EQUAT
+ ! 0, coarse grid is binned with preferential treatment of the north fold
+ ! 1, coarse grid is binned with centering at the equator
+ ! Symmetry with nn_facty being oddnumbered. Asymmetry with evennumbered nn_facty.
+ nn_msh_crs = 1 ! create (=1) a mesh file or not (=0)
+ nn_crs_kz = 0 ! 0, MEAN of volume boxes
+ ! 1, MAX of boxes
+ ! 2, MIN of boxes
+ ln_crs_wn = .true. ! wn coarsened (T) or computed using horizontal divergence ( F )
+/
Index: trunk/DOC/TexFiles/Namelist/namctl
===================================================================
 trunk/DOC/TexFiles/Namelist/namctl (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namctl (revision 6289)
@@ 13,4 +13,4 @@
! (no physical validity of the results)
nn_timing = 0 ! timing by routine activated (=1) creates timing.output file, or not (=0)
 nn_diacfl = 0 ! Write out cfl diagnostics (=1) in cfl_diagnostics.ascii, or not (=0)
+ nn_diacfl = 0 ! Write out CFL diagnostics (=1) in cfl_diagnostics.ascii, or not (=0)
/
Index: trunk/DOC/TexFiles/Namelist/namdia_harm
===================================================================
 trunk/DOC/TexFiles/Namelist/namdia_harm (revision 6140)
+++ (revision )
@@ 1,9 +1,0 @@
!
&nam_diaharm ! Harmonic analysis of tidal constituents ('key_diaharm')
!
 nit000_han = 1 ! First time step used for harmonic analysis
 nitend_han = 75 ! Last time step used for harmonic analysis
 nstep_han = 15 ! Time step frequency for harmonic analysis
 tname(1) = 'M2' ! Name of tidal constituents
 tname(2) = 'K1'
/
Index: trunk/DOC/TexFiles/Namelist/namdiu
===================================================================
 trunk/DOC/TexFiles/Namelist/namdiu (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namdiu (revision 6289)
@@ 1,6 +1,6 @@
!
&namdiu ! Control of diurnal SST models
+&namdiu ! Cool skin and warm layer models
!
 ln_diurnal = .TRUE. ! Turn diurnal models on/off
 ln_diurnal_only = .FALSE. ! Only run diurnal model, dynamical models switched off
+ ln_diurnal = .false. !
+ ln_diurnal_only = .false. !
/
Index: trunk/DOC/TexFiles/Namelist/namdom
===================================================================
 trunk/DOC/TexFiles/Namelist/namdom (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namdom (revision 6289)
@@ 13,9 +13,4 @@
rn_rdt = 5760. ! time step for the dynamics (and tracer if nn_acc=0)
rn_atfp = 0.1 ! asselin time filter parameter
 nn_acc = 0 ! acceleration of convergence : =1 used, rdt < rdttra(k)
 ! =0, not used, rdt = rdttra
 rn_rdtmin = 28800. ! minimum time step on tracers (used if nn_acc=1)
 rn_rdtmax = 28800. ! maximum time step on tracers (used if nn_acc=1)
 rn_rdth = 800. ! depth variation of tracer time step (used if nn_acc=1)
ln_crs = .false. ! Logical switch for coarsening module
jphgr_msh = 0 ! type of horizontal mesh
Index: trunk/DOC/TexFiles/Namelist/namdyn_nept
===================================================================
 trunk/DOC/TexFiles/Namelist/namdyn_nept (revision 6140)
+++ (revision )
@@ 1,15 +1,0 @@
!
&namdyn_nept ! Neptune effect (simplified: lateral and vertical diffusions removed)
!
 ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model
 ln_neptsimp = .false. ! yes/no use simplified neptune

 ln_smooth_neptvel = .false. ! yes/no smooth zunep, zvnep
 rn_tslse = 1.2e4 ! value of lengthscale L at the equator
 rn_tslsp = 3.0e3 ! value of lengthscale L at the pole
 ! Specify whether to ramp down the Neptune velocity in shallow
 ! water, and if so the depth range controlling such ramping down
 ln_neptramp = .true. ! ramp down Neptune velocity in shallow water
 rn_htrmin = 100.0 ! min. depth of transition range
 rn_htrmax = 200.0 ! max. depth of transition range
/
Index: trunk/DOC/TexFiles/Namelist/namdyn_vor
===================================================================
 trunk/DOC/TexFiles/Namelist/namdyn_vor (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namdyn_vor (revision 6289)
@@ 7,4 +7,4 @@
ln_dynvor_een = .false. ! energy & enstrophy scheme
nn_een_e3f = 1 ! e3f = masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
 ln_dynvor_msk = .false. ! vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) ! PLEASE DO NOT USE
+ ln_dynvor_msk = .false. ! vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) ! PLEASE DO NOT ACTIVATE
/
Index: trunk/DOC/TexFiles/Namelist/nammpp_dyndist
===================================================================
 trunk/DOC/TexFiles/Namelist/nammpp_dyndist (revision 6140)
+++ (revision )
@@ 1,7 +1,0 @@
!
&nammpp_dyndist ! Massively Parallel Distribution for AGRIF zoom ("key_agrif" && "key_mpp_dyndist")
!
 jpni = 1 ! jpni number of processors following i
 jpnj = 1 ! jpnj number of processors following j
 jpnij = 1 ! jpnij number of local domains
/
Index: trunk/DOC/TexFiles/Namelist/namobc
===================================================================
 trunk/DOC/TexFiles/Namelist/namobc (revision 6140)
+++ (revision )
@@ 1,21 +1,0 @@
!
&namobc ! open boundaries parameters ("key_obc")
!
 ln_obc_clim = .false. ! climatological obc data files (T) or not (F)
 ln_vol_cst = .true. ! impose the total volume conservation (T) or not (F)
 ln_obc_fla = .false. ! Flather open boundary condition
 nn_obcdta = 1 ! = 0 the obc data are equal to the initial state
 ! = 1 the obc data are read in 'obc.dta' files
 cn_obcdta = 'annual' ! set to annual if obc datafile hold 1 year of data
 ! set to monthly if obc datafile hold 1 month of data
 rn_dpein = 1. ! damping time scale for inflow at east open boundary
 rn_dpwin = 1. !    west  
 rn_dpnin = 1. !    north  
 rn_dpsin = 1. !    south  
 rn_dpeob = 3000. ! time relaxation (days) for the east open boundary
 rn_dpwob = 15. !    west  
 rn_dpnob = 3000. !    north  
 rn_dpsob = 15. !    south  
 rn_volemp = 1. ! = 0 the total volume change with the surface flux (EPR)
 ! = 1 the total volume remains constant
/
Index: trunk/DOC/TexFiles/Namelist/namobs
===================================================================
 trunk/DOC/TexFiles/Namelist/namobs (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namobs (revision 6289)
@@ 2,18 +2,19 @@
&namobs ! observation usage switch
!
 ln_diaobs = .true. ! Logical switch for the observation operator
 ln_t3d = .true. ! Logical switch for T profile observations
 ln_s3d = .true. ! Logical switch for S profile observations
 ln_sla = .true. ! Logical switch for SLA observations
 ln_sst = .true. ! Logical switch for SST observations
 ln_sic = .false. ! Logical switch for Sea Ice observations
 ln_vel3d = .false. ! Logical switch for velocity observations
 ln_altbias = .false. ! Logical switch for altimeter bias correction
 ln_nea = .false. ! Logical switch for rejection of observations near land
 ln_grid_global = .true. ! Logical switch for global distribution of observations
 ln_grid_search_lookup = .false. ! Logical switch for obs grid search w/lookup table
 ln_ignmis = .true. ! Logical switch for ignoring missing files
 ln_s_at_t = .false. ! Logical switch for computing model S at T obs if not there
 ln_sstnight = .false. ! Logical switch for calculating nighttime average for SST obs
+ ln_diaobs = .false. ! Logical switch for the observation operator
+ ln_t3d = .false. ! Logical switch for T profile observations
+ ln_s3d = .false. ! Logical switch for S profile observations
+ ln_sla = .false. ! Logical switch for SLA observations
+ ln_sst = .false. ! Logical switch for SST observations
+ ln_sic = .false. ! Logical switch for Sea Ice observations
+ ln_vel3d = .false. ! Logical switch for velocity observations
+ ln_altbias = .false. ! Logical switch for altimeter bias correction
+ ln_nea = .false. ! Logical switch for rejection of observations near land
+ ln_grid_global = .true. ! Logical switch for global distribution of observations
+ ln_grid_search_lookup = .false. ! Logical switch for obs grid search w/lookup table
+ ln_ignmis = .true. ! Logical switch for ignoring missing files
+ ln_s_at_t = .false. ! Logical switch for computing model S at T obs if not there
+ ln_sstnight = .false. ! Logical switch for calculating nighttime average for SST obs
+! All of the *files* variables below are arrays. Use namelist_cfg to add more files
cn_profbfiles = 'profiles_01.nc' ! Profile feedback input observation file names
cn_slafbfiles = 'sla_01.nc' ! SLA feedback input observation file names
@@ 24,6 +25,6 @@
cn_gridsearchfile = 'gridsearch.nc' ! Grid search file name
rn_gridsearchres = 0.5 ! Grid search resolution
 rn_dobsini = 20150101.000000 ! Initial date in window YYYYMMDD.HHMMSS
 rn_dobsend = 20150102.000000 ! Final date in window YYYYMMDD.HHMMSS
+ rn_dobsini = 00010101.000000 ! Initial date in window YYYYMMDD.HHMMSS
+ rn_dobsend = 00010102.000000 ! Final date in window YYYYMMDD.HHMMSS
nn_1dint = 0 ! Type of vertical interpolation method
nn_2dint = 0 ! Type of horizontal interpolation method
@@ 32,3 +33,5 @@
rn_mdtcutoff = 65.0 ! MDT cutoff for computed correction
nn_profdavtypes = 1 ! Profile daily average types  array
+ ln_sstbias = .false.
+ cn_sstbias_files = 'sstbias.nc'
/
Index: trunk/DOC/TexFiles/Namelist/namobs_example
===================================================================
 trunk/DOC/TexFiles/Namelist/namobs_example (revision 6140)
+++ (revision )
@@ 1,34 +1,0 @@
!
&namobs ! observation usage switch
!
 ln_diaobs = .true. ! Logical switch for the observation operator
 ln_t3d = .true. ! Logical switch for T profile observations
 ln_s3d = .true. ! Logical switch for S profile observations
 ln_sla = .false. ! Logical switch for SLA observations
 ln_sst = .false. ! Logical switch for SST observations
 ln_sic = .false. ! Logical switch for Sea Ice observations
 ln_vel3d = .false. ! Logical switch for velocity observations
 ln_altbias = .false. ! Logical switch for altimeter bias correction
 ln_nea = .false. ! Logical switch for rejection of observations near land
 ln_grid_global = .true. ! Logical switch for global distribution of observations
 ln_grid_search_lookup = .false. ! Logical switch for obs grid search w/lookup table
 ln_ignmis = .true. ! Logical switch for ignoring missing files
 ln_s_at_t = .false. ! Logical switch for computing model S at T obs if not there
 ln_sstnight = .false. ! Logical switch for calculating nighttime average for SST obs
 cn_profbfiles = 'profiles_01.nc' ! Profile feedback input observation file names
 cn_slafbfiles = 'sla_01.nc' ! SLA feedback input observation file names
 cn_sstfbfiles = 'sst_01.nc' ! SST feedback input observation file names
 cn_sicfbfiles = 'sic_01.nc' ! SIC feedback input observation file names
 cn_velfbfiles = 'vel_01.nc' ! Velocity feedback input observation file names
 cn_altbiasfile = 'altbias.nc' ! Altimeter bias input file name
 cn_gridsearchfile = 'gridsearch.nc' ! Grid search file name
 rn_gridsearchres = 0.5 ! Grid search resolution
 rn_dobsini = 00010101.000000 ! Initial date in window YYYYMMDD.HHMMSS
 rn_dobsend = 00010102.000000 ! Final date in window YYYYMMDD.HHMMSS
 nn_1dint = 0 ! Type of vertical interpolation method
 nn_2dint = 0 ! Type of horizontal interpolation method
 nn_msshc = 0 ! MSSH correction scheme
 rn_mdtcorr = 1.61 ! MDT correction
 rn_mdtcutoff = 65.0 ! MDT cutoff for computed correction
 nn_profdavtypes = 1 ! Profile daily average types  array
/
Index: trunk/DOC/TexFiles/Namelist/namrun
===================================================================
 trunk/DOC/TexFiles/Namelist/namrun (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namrun (revision 6289)
@@ 7,16 +7,17 @@
nn_itend = 5475 ! last time step (std 5475)
nn_date0 = 010101 ! date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
+ nn_time0 = 0 ! initial time of day in hhmm
nn_leapy = 0 ! Leap year calendar (1) or not (0)
ln_rstart = .false. ! start from rest (F) or from a restart file (T)
 nn_euler = 1 ! = 0 : start with forward time step if ln_rstart=T
 nn_rstctl = 0 ! restart control ==> activated only if ln_rstart=T
 ! = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
 ! = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart
 ! = 2 nn_date0 read in restart ; nn_it000 : check consistancy between namelist and restart
 cn_ocerst_in = "restart" ! suffix of ocean restart name (input)
 cn_ocerst_indir = "." ! directory from which to read input ocean restarts
 cn_ocerst_out = "restart" ! suffix of ocean restart name (output)
 cn_ocerst_outdir = "." ! directory in which to write output ocean restarts
 ln_iscpl = .false. ! cavity evolution forcing or coupling to ice sheet model (ln_iscpl = T => namsbc_iscpl)
+ nn_euler = 1 ! = 0 : start with forward time step if ln_rstart=T
+ nn_rstctl = 0 ! restart control ==> activated only if ln_rstart=T
+ ! ! = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
+ ! ! = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart
+ ! ! = 2 nn_date0 read in restart ; nn_it000 : check consistancy between namelist and restart
+ cn_ocerst_in = "restart" ! suffix of ocean restart name (input)
+ cn_ocerst_indir = "." ! directory from which to read input ocean restarts
+ cn_ocerst_out = "restart" ! suffix of ocean restart name (output)
+ cn_ocerst_outdir= "." ! directory in which to write output ocean restarts
+ ln_iscpl = .false. ! cavity evolution forcing or coupling to ice sheet model
nn_istate = 0 ! output the initial state (1) or not (0)
ln_rst_list = .false. ! output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
@@ 24,5 +25,4 @@
nn_stocklist = 0,0,0,0,0,0,0,0,0,0 ! List of timesteps when a restart file is to be written
nn_write = 5475 ! frequency of write in the output file (modulo referenced to nn_it000)
 ln_dimgnnn = .false. ! DIMG file format: 1 file for all processors (F) or by processor (T)
ln_mskland = .false. ! mask land points in NetCDF outputs (costly: + ~15%)
ln_cfmeta = .false. ! output additional data to netCDF files required for compliance with the CF metadata standard
Index: trunk/DOC/TexFiles/Namelist/namsbc
===================================================================
 trunk/DOC/TexFiles/Namelist/namsbc (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namsbc (revision 6289)
@@ 31,5 +31,5 @@
! Misc. options of sbc :
ln_traqsr = .true. ! Light penetration in the ocean (T => fill namtra_qsr )
 ln_dm2dc = .false. ! daily mean to diurnal cycle on short wave
+ ln_dm2dc = .false. ! daily mean to diurnal cycle on short wave
ln_rnf = .true. ! runoffs (T => fill namsbc_rnf)
ln_ssr = .true. ! Sea Surface Restoring on T and/or S (T => fill namsbc_ssr)
@@ 38,6 +38,6 @@
! =2 annual global mean of epr set to zero
ln_apr_dyn = .false. ! Patm gradient added in ocean & ice Eqs. (T => fill namsbc_apr )
 ln_isf = .false. ! ice shelf melting/freezing (T => fill namsbc_isf)
 ln_wave = .false. ! coupling with surface wave (T => fill namsbc_wave)
+ ln_isf = .false. ! ice shelf (T => fill namsbc_isf)
+ ln_wave = .false. ! coupling with surface wave (T => fill namsbc_wave)
nn_lsm = 0 ! =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
! =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field)
Index: trunk/DOC/TexFiles/Namelist/namsbc_iscpl
===================================================================
 trunk/DOC/TexFiles/Namelist/namsbc_iscpl (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namsbc_iscpl (revision 6289)
@@ 1,6 +1,6 @@
!
&namsbc_iscpl ! land ice / ocean coupling option
+&namsbc_iscpl ! land ice / ocean coupling option
!
 nn_drown = 10 ! number of iterations of the extrapolation loop (fill the new wet cells)
+ nn_drown = 10 ! number of iteration of the extrapolation loop (fill the new wet cells)
ln_hsb = .false. ! activate conservation module (conservation exact after a time of rn_fiscpl)
nn_fiscpl = 43800 ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency)
Index: trunk/DOC/TexFiles/Namelist/namsbc_isf
===================================================================
 trunk/DOC/TexFiles/Namelist/namsbc_isf (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namsbc_isf (revision 6289)
@@ 2,33 +2,32 @@
&namsbc_isf ! Top boundary layer (ISF) (nn_isf >0)
!
 nn_isf = 99 ! ice shelf melting/freezing (/=0 => fill namsbc_isf)
 ! 1 = ISF explicit (cavity open)
 ! 2 = bg03 parametrisation (cavity closed)
 ! 3 = ISF specified in depth along the calving front (cavity closed)
 ! 4 = ISF melting specified (cavity open)
+! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask !
+! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename !
+! nn_isf == 4
+ sn_fwfisf = 'rnfisf' , 12 ,'sowflisf', .false. , .true. , 'yearly' , '' , '' , ''
+! nn_isf == 3
+ sn_rnfisf = 'rnfisf' , 12 ,'sofwfisf', .false. , .true. , 'yearly' , '' , '' , ''
+! nn_isf == 2 and 3
+ sn_depmax_isf='rnfisf' , 12 ,'sozisfmax', .false. , .true. , 'yearly' , '' , '' , ''
+ sn_depmin_isf='rnfisf' , 12 ,'sozisfmin', .false. , .true. , 'yearly' , '' , '' , ''
+! nn_isf == 2
+ sn_Leff_isf = 'rnfisf' , 12 ,'Leff' , .false. , .true. , 'yearly' , '' , '' , ''
+!
+! for all case
+ nn_isf = 1 ! ice shelf melting/freezing
+ ! 1 = presence of ISF 2 = bg03 parametrisation
+ ! 3 = rnf file for isf 4 = ISF fwf specified
! option 1 and 4 need ln_isfcav = .true. (domzgr)
! ! file name ! frequency (hours) ! variable ! time interpol. ! clim ! 'yearly'/ ! weights ! rotation !
! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing !
! nn_isf == 4
 sn_fwfisf = 'rnfisf', 12 ,'sowflisf', .false. , .true. , 'yearly' , '' , ''
! nn_isf == 3
 sn_rnfisf = 'rnfisf', 12 ,'sofwfisf', .false. , .true. , 'yearly' , '' , ''
! nn_isf == 2 or 3
 sn_depmax_isf = 'rnfisf', 12 ,'sozisfmax', .false. , .true. , 'yearly' , '' , ''
 sn_depmin_isf = 'rnfisf', 12 ,'sozisfmin', .false. , .true. , 'yearly' , '' , ''
! nn_isf == 2
 sn_Leff_isf = 'rnfisf' , 12 ,'Leff' , .false. , .true. , 'yearly' , '' , ''

! only for nn_isf = 1 or 2
rn_gammat0 = 1.e4 ! gammat coefficient used in blk formula
rn_gammas0 = 1.e4 ! gammas coefficient used in blk formula
! only for nn_isf = 1 or 4
 nn_isfblk = 1 ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??)
+ rn_hisf_tbl = 30. ! thickness of the top boundary layer (Losh et al. 2008)
+ ! 0 => thickness of the tbl = thickness of the first wet cell
! only for nn_isf = 1
 rn_hisf_tbl = 30. ! thickness of the top boundary layer (Losh et al. 2008)
 ! 0 => thickness of the tbl = thickness of the first wet cell
+ nn_isfblk = 1 ! 1 ISOMIP like: 2 equations formulation (Hunter et al., 2006)
+ ! 2 ISOMIP+ like: 3 equations formulation (AsayDavis et al., 2015)
nn_gammablk = 1 ! 0 = cst Gammat (= gammat/s)
! 1 = velocity dependend Gamma (u* * gammat/s) (Jenkins et al. 2010)
 ! if you want to keep the cd as in global config, adjust rn_gammat0 to compensate
 ! 2 = velocity and stability dependent Gamma Holland et al. 1999
+ ! 2 = velocity and stability dependent Gamma (Holland et al. 1999)
/
Index: trunk/DOC/TexFiles/Namelist/namsol
===================================================================
 trunk/DOC/TexFiles/Namelist/namsol (revision 6140)
+++ (revision )
@@ 1,13 +1,0 @@
!
&namsol ! elliptic solver / island / free surface
!
 nn_solv = 1 ! elliptic solver: =1 preconditioned conjugate gradient (pcg)
 ! =2 successiveoverrelaxation (sor)
 nn_sol_arp = 0 ! absolute/relative (0/1) precision convergence test
 rn_eps = 1.e6 ! absolute precision of the solver
 nn_nmin = 300 ! minimum of iterations for the SOR solver
 nn_nmax = 800 ! maximum of iterations for the SOR solver
 nn_nmod = 10 ! frequency of test for the SOR solver
 rn_resmax = 1.e10 ! absolute precision for the SOR solver
 rn_sor = 1.92 ! optimal coefficient for SOR solver (to be adjusted with the domain)
/
Index: trunk/DOC/TexFiles/Namelist/namsplit
===================================================================
 trunk/DOC/TexFiles/Namelist/namsplit (revision 6140)
+++ (revision )
@@ 1,15 +1,0 @@
!
&namsplit ! time splitting parameters ("key_dynspg_ts")
!
 ln_bt_fw = .TRUE. ! Forward integration of barotropic equations
 ln_bt_av = .TRUE. ! Time filtering of barotropic variables
 ln_bt_nn_auto = .TRUE. ! Set nn_baro automatically to be just below
 ! a user defined maximum courant number (rn_bt_cmax)
 nn_baro = 30 ! Number of iterations of barotropic mode
 ! during rn_rdt seconds. Only used if ln_bt_nn_auto=F
 rn_bt_cmax = 0.8 ! Maximum courant number allowed if ln_bt_nn_auto=T
 nn_bt_flt = 1 ! Time filter choice
 ! = 0 None
 ! = 1 Boxcar over nn_baro barotropic steps
 ! = 2 Boxcar over 2*nn_baro " "
/
Index: trunk/DOC/TexFiles/Namelist/namsto
===================================================================
 trunk/DOC/TexFiles/Namelist/namsto (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namsto (revision 6289)
@@ 0,0 +1,16 @@
+!
+&namsto ! Stochastic parametrization of EOS (default: NO)
+!
+ ln_sto_eos = .false. ! stochastic equation of state
+ nn_sto_eos = 1 ! number of independent random walks
+ rn_eos_stdxy = 1.4 ! random walk horz. standard deviation (in grid points)
+ rn_eos_stdz = 0.7 ! random walk vert. standard deviation (in grid points)
+ rn_eos_tcor = 1440. ! random walk time correlation (in timesteps)
+ nn_eos_ord = 1 ! order of autoregressive processes
+ nn_eos_flt = 0 ! passes of Laplacian filter
+ rn_eos_lim = 2.0 ! limitation factor (default = 3.0)
+ ln_rststo = .false. ! start from mean parameter (F) or from restart file (T)
+ ln_rstseed = .true. ! read seed of RNG from restart file
+ cn_storst_in = "restart_sto" ! suffix of stochastic parameter restart file (input)
+ cn_storst_out = "restart_sto" ! suffix of stochastic parameter restart file (output)
+/
Index: trunk/DOC/TexFiles/Namelist/namtra_adv_mle
===================================================================
 trunk/DOC/TexFiles/Namelist/namtra_adv_mle (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namtra_adv_mle (revision 6289)
@@ 0,0 +1,13 @@
+!
+&namtra_adv_mle ! mixed layer eddy parametrisation (FoxKemper param) (default: NO)
+!
+ ln_mle = .false. ! (T) use the Mixed Layer Eddy (MLE) parameterisation
+ rn_ce = 0.06 ! magnitude of the MLE (typical value: 0.06 to 0.08)
+ nn_mle = 1 ! MLE type: =0 standard FoxKemper ; =1 new formulation
+ rn_lf = 5.e+3 ! typical scale of mixed layer front (meters) (case rn_mle=0)
+ rn_time = 172800. ! time scale for mixing momentum across the mixed layer (seconds) (case rn_mle=0)
+ rn_lat = 20. ! reference latitude (degrees) of MLE coef. (case rn_mle=1)
+ nn_mld_uv = 0 ! space interpolation of MLD at u & vpts (0=min,1=averaged,2=max)
+ nn_conv = 0 ! =1 no MLE in case of convection ; =0 always MLE
+ rn_rho_c_mle = 0.01 ! delta rho criterion used to calculate MLD for FK
+/
Index: trunk/DOC/TexFiles/Namelist/namtra_bbc
===================================================================
 trunk/DOC/TexFiles/Namelist/namtra_bbc (revision 6140)
+++ (revision )
@@ 1,9 +1,0 @@
!
&nambbc ! bottom temperature boundary condition
!
 ln_trabbc = .true. ! Apply a geothermal heating at the ocean bottom
 nn_geoflx = 2 ! geothermal heat flux: = 0 no flux
 ! = 1 constant flux
 ! = 2 variable flux (read in geothermal_heating.nc in mW/m2)
 rn_geoflx_cst = 86.4e3 ! Constant value of geothermal heat flux [W/m2]
/
Index: trunk/DOC/TexFiles/Namelist/namtra_ldf
===================================================================
 trunk/DOC/TexFiles/Namelist/namtra_ldf (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namtra_ldf (revision 6289)
@@ 1,7 +1,7 @@
!
&namtra_ldf ! lateral diffusion scheme for tracers (default: NO diffusion)
+&namtra_ldf ! lateral diffusion scheme for tracers (choice required)
!
! ! Operator type:
 ! ! no diffusion: set ln_traldf_lap=..._blp=F
+ ln_traldf_NONE = .false. ! no operator: no lateral diffusion applied
ln_traldf_lap = .false. ! laplacian operator
ln_traldf_blp = .false. ! bilaplacian operator
Index: trunk/DOC/TexFiles/Namelist/namtra_ldfeiv
===================================================================
 trunk/DOC/TexFiles/Namelist/namtra_ldfeiv (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namtra_ldfeiv (revision 6289)
@@ 0,0 +1,14 @@
+!
+&namtra_ldfeiv ! eddy induced velocity param. (default: NO)
+!
+ ln_ldfeiv =.false. ! use eddy induced velocity parameterization
+ ln_ldfeiv_dia =.false. ! diagnose eiv stream function and velocities
+ rn_aeiv_0 = 2000. ! eddy induced velocity coefficient [m2/s]
+ nn_aei_ijk_t = 21 ! space/time variation of the eiv coeficient
+ ! ! =20 (=30) read in eddy_induced_velocity_2D.nc (..._3D.nc) file
+ ! ! = 0 constant
+ ! ! = 10 F(k) =ldf_c1d
+ ! ! = 20 F(i,j) =ldf_c2d
+ ! ! = 21 F(i,j,t) =Treguier et al. JPO 1997 formulation
+ ! ! = 30 F(i,j,k) =ldf_c2d + ldf_c1d
+/
Index: trunk/DOC/TexFiles/Namelist/namwad
===================================================================
 trunk/DOC/TexFiles/Namelist/namwad (revision 6289)
+++ trunk/DOC/TexFiles/Namelist/namwad (revision 6289)
@@ 0,0 +1,9 @@
+!
+&namwad ! Wetting and drying
+!
+ ln_wd = .false. ! T/F activation of wetting and drying
+ rn_wdmin1 = 0.1 ! Minimum wet depth on dried cells
+ rn_wdmin2 = 0.01 ! Tolerance of min wet depth on dried cells
+ rn_wdld = 20.0 ! Land elevation below which wetting/drying is allowed
+ nn_wdit = 10 ! Max iterations for W/D limiter
+/
Index: trunk/DOC/TexFiles/Namelist/namzdf_kpp
===================================================================
 trunk/DOC/TexFiles/Namelist/namzdf_kpp (revision 6140)
+++ (revision )
@@ 1,13 +1,0 @@
!
&namzdf_kpp ! KProfile Parameterization dependent vertical mixing ("key_zdfkpp", and optionally:
! "key_kppcustom" or "key_kpplktb")
 ln_kpprimix = .true. ! shear instability mixing
 rn_difmiw = 1.0e04 ! constant internal wave viscosity [m2/s]
 rn_difsiw = 0.1e04 ! constant internal wave diffusivity [m2/s]
 rn_riinfty = 0.8 ! local Richardson Number limit for shear instability
 rn_difri = 0.0050 ! maximum shear mixing at Rig = 0 [m2/s]
 rn_bvsqcon = 0.01e07 ! BruntVaisala squared for maximum convection [1/s2]
 rn_difcon = 1. ! maximum mixing in interior convection [m2/s]
 nn_avb = 0 ! horizontal averaged (=1) or not (=0) on avt and amv
 nn_ave = 1 ! constant (=0) or profile (=1) background on avt
/
Index: trunk/DOC/TexFiles/Namelist/namzdf_ric
===================================================================
 trunk/DOC/TexFiles/Namelist/namzdf_ric (revision 6140)
+++ trunk/DOC/TexFiles/Namelist/namzdf_ric (revision 6289)
@@ 10,4 +10,4 @@
rn_wtmix = 10.0 ! vertical eddy viscosity coeff [m2/s] in the mixedlayer
rn_wvmix = 10.0 ! vertical eddy diffusion coeff [m2/s] in the mixedlayer
 ln_mldw = .true. ! Flag to use or not the mized layer depth param.
+ ln_mldw = .true. ! Flag to use or not the mixed layer depth param.
/
Index: trunk/DOC/TexFiles/Namelist_tools/namcfg_orca1
===================================================================
 trunk/DOC/TexFiles/Namelist_tools/namcfg_orca1 (revision 6289)
+++ trunk/DOC/TexFiles/Namelist_tools/namcfg_orca1 (revision 6289)
@@ 0,0 +1,12 @@
+!
+&namcfg ! parameters of the configuration
+!
+ cp_cfg = "orca" ! name of the configuration
+ jp_cfg = 1 ! resolution of the configuration
+ jpidta = 362 ! 1st lateral dimension ( >= jpi )
+ jpjdta = 292 ! 2nd " " ( >= jpj )
+ jpkdta = 75 ! number of levels ( >= jpk )
+ jpiglo = 362 ! 1st dimension of global domain > i =jpidta
+ jpjglo = 292 ! 2nd   > j =jpjdta
+ jperio = 6 ! lateral cond. type (between 0 and 6)
+ ln_use_jattr = .true. ! use (T) the file attribute: open_ocean_jstart if present
Index: trunk/DOC/TexFiles/Namelist_tools/namelist_dmp
===================================================================
 trunk/DOC/TexFiles/Namelist_tools/namelist_dmp (revision 6289)
+++ trunk/DOC/TexFiles/Namelist_tools/namelist_dmp (revision 6289)
@@ 0,0 +1,24 @@
+&nam_dmp_create
+ cp_cfg = 'orca' ! Name of model grid (orca and C1D have special options  otherwise ignored)
+ cp_cfz = 'antarctic' ! Name of zoom configuration (arctic and antarctic have some special treatment if lzoom=.true.)
+ jp_cfg = 2 ! Resolution of the model (used for med_red_seas damping)
+ lzoom = .false. ! Zoom configuration or not
+ ln_full_field = .false. ! Calculate coefficient over whole of domain
+ ln_med_red_seas = .true. ! Damping in Med/Red Seas (or local modifications here if ln_full_field=.true.)
+ ln_old_31_lev_code = .true. ! Replicate behaviour of old online code for 31 level model (Med/Red seas damping based on level number instead of depth)
+ ln_coast = .true. ! Reduce near to coastlines
+ ln_zero_top_layer = .true. ! No damping in top layer
+ ln_custom = .false. ! Call "custom" module to apply user modifications to the damping coefficient field
+ nn_hdmp = 10 ! Damp poleward of this latitude (smooth transition up to maximum damping)
+ pn_surf = 0.25 ! Surface Relaxation timescale (days)
+ pn_bot = 0.25 ! Bottom relaxation timescale (days)
+ pn_dep = 1000 ! Transition depth from upper to deep ocean
+ jperio = 2 ! Lateral boundary condition (as specified in namelist_cfg for model run).
+/
+
+&nam_zoom_dmp
+ lzoom_n = .false. ! Open boundary had northern edge?
+ lzoom_e = .false. ! Open boundary at eastern edge?
+ lzoom_w = .false. ! Open boundary at western edge?
+ lzoom_s = .false. ! Open boundary at southern edge?
+/
Index: trunk/DOC/TexFiles/Namelist_tools/namobs_example
===================================================================
 trunk/DOC/TexFiles/Namelist_tools/namobs_example (revision 6289)
+++ trunk/DOC/TexFiles/Namelist_tools/namobs_example (revision 6289)
@@ 0,0 +1,34 @@
+!
+&namobs ! observation usage switch
+!
+ ln_diaobs = .true. ! Logical switch for the observation operator
+ ln_t3d = .true. ! Logical switch for T profile observations
+ ln_s3d = .true. ! Logical switch for S profile observations
+ ln_sla = .false. ! Logical switch for SLA observations
+ ln_sst = .false. ! Logical switch for SST observations
+ ln_sic = .false. ! Logical switch for Sea Ice observations
+ ln_vel3d = .false. ! Logical switch for velocity observations
+ ln_altbias = .false. ! Logical switch for altimeter bias correction
+ ln_nea = .false. ! Logical switch for rejection of observations near land
+ ln_grid_global = .true. ! Logical switch for global distribution of observations
+ ln_grid_search_lookup = .false. ! Logical switch for obs grid search w/lookup table
+ ln_ignmis = .true. ! Logical switch for ignoring missing files
+ ln_s_at_t = .false. ! Logical switch for computing model S at T obs if not there
+ ln_sstnight = .false. ! Logical switch for calculating nighttime average for SST obs
+ cn_profbfiles = 'profiles_01.nc' ! Profile feedback input observation file names
+ cn_slafbfiles = 'sla_01.nc' ! SLA feedback input observation file names
+ cn_sstfbfiles = 'sst_01.nc' ! SST feedback input observation file names
+ cn_sicfbfiles = 'sic_01.nc' ! SIC feedback input observation file names
+ cn_velfbfiles = 'vel_01.nc' ! Velocity feedback input observation file names
+ cn_altbiasfile = 'altbias.nc' ! Altimeter bias input file name
+ cn_gridsearchfile = 'gridsearch.nc' ! Grid search file name
+ rn_gridsearchres = 0.5 ! Grid search resolution
+ rn_dobsini = 00010101.000000 ! Initial date in window YYYYMMDD.HHMMSS
+ rn_dobsend = 00010102.000000 ! Final date in window YYYYMMDD.HHMMSS
+ nn_1dint = 0 ! Type of vertical interpolation method
+ nn_2dint = 0 ! Type of horizontal interpolation method
+ nn_msshc = 0 ! MSSH correction scheme
+ rn_mdtcorr = 1.61 ! MDT correction
+ rn_mdtcutoff = 65.0 ! MDT cutoff for computed correction
+ nn_profdavtypes = 1 ! Profile daily average types  array
+/
Index: trunk/DOC/TexFiles/namelist_split.sh
===================================================================
 trunk/DOC/TexFiles/namelist_split.sh (revision 6289)
+++ trunk/DOC/TexFiles/namelist_split.sh (revision 6289)
@@ 0,0 +1,247 @@
+#! /bin/sh
+#
+# usage for NEMO doc to create an update of the Namelist directory :
+# 1 delete the existing directory (You can also choose to save it somewhere)
+# rm rf Namelist
+# 2 create the updated Namelist directory from the SHARED/namelist_ref :
+# ./namelist_split.sh i ../../NEMOGCM/CONFIG/SHARED/namelist_ref o Namelist
+#
+# .. _namelist_split.sh:
+#
+# =================
+# namelist_split.sh
+# =================
+#
+# 
+# split a namelist
+# 
+#
+# SYNOPSIS
+# ========
+#
+# .. codeblock:: bash
+#
+# namelist_split.sh i namelist o dirout [n]
+#
+# DESCRIPTION
+# ===========
+#
+# Split a namelist file (NEMO convention syntax) given in parameter
+# into files in a given output directory.
+#
+# .. option:: i
+# .. option:: o