Changeset 6289 for trunk/DOC/TexFiles/Chapters
- Timestamp:
- 2016-02-05T00:47:05+01:00 (8 years ago)
- Location:
- trunk/DOC/TexFiles/Chapters
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex
r6140 r6289 14 14 the earth climate system over a wide range of space and time scales. 15 15 Prognostic variables are the three-dimensional velocity field, a non-linear sea surface height, 16 the \textit{Conservative} Temperature and the \textit{Absolut } Salinity.16 the \textit{Conservative} Temperature and the \textit{Absolute} Salinity. 17 17 In the horizontal direction, the model uses a curvilinear orthogonal grid and in the vertical direction, 18 18 a full or partial step $z$-coordinate, or $s$-coordinate, or a mixture of the two. … … 31 31 interactions avec les autres composantes du syst\`{e}me climatique terrestre. 32 32 Les variables pronostiques sont le champ tridimensionnel de vitesse, une hauteur de la mer 33 lin\'{e}aire, la T temperature Conservative et la Salinit\'{e} Absolue.33 lin\'{e}aire, la Temp\'{e}erature Conservative et la Salinit\'{e} Absolue. 34 34 La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une 35 35 coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore -
trunk/DOC/TexFiles/Chapters/Annex_C.tex
r6140 r6289 1103 1103 The discrete formulation of the horizontal diffusion of momentum ensures the 1104 1104 conservation of potential vorticity and the horizontal divergence, and the 1105 dissipation of the square of these quantities ( i.e.enstrophy and the1105 dissipation of the square of these quantities ($i.e.$ enstrophy and the 1106 1106 variance of the horizontal divergence) as well as the dissipation of the 1107 1107 horizontal kinetic energy. In particular, when the eddy coefficients are … … 1127 1127 &\int \limits_D \frac{1} {e_3 } \textbf{k} \cdot \nabla \times 1128 1128 \Bigl[ \nabla_h \left( A^{\,lm}\;\chi \right) 1129 - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv = 01130 \end{flalign*}1129 - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv \\ 1130 %\end{flalign*} 1131 1131 %%%%%%%%%% recheck here.... (gm) 1132 \begin{flalign*}1133 = \int \limits_D -\frac{1} {e_3 } \textbf{k} \cdot \nabla \times1134 \Bigl[ \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv &&&\\1135 \end{flalign*}1136 \begin{flalign*}1132 %\begin{flalign*} 1133 =& \int \limits_D -\frac{1} {e_3 } \textbf{k} \cdot \nabla \times 1134 \Bigl[ \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \Bigr]\;dv \\ 1135 %\end{flalign*} 1136 %\begin{flalign*} 1137 1137 \equiv& \sum\limits_{i,j} 1138 1138 \left\{ 1139 \delta_{i+1/2} 1140 \left[ 1141 \frac {e_{2v}} {e_{1v}\,e_{3v}} \delta_i 1142 \left[ A_f^{\,lm} e_{3f} \zeta \right] 1143 \right] 1144 + \delta_{j+1/2} 1145 \left[ 1146 \frac {e_{1u}} {e_{2u}\,e_{3u}} \delta_j 1147 \left[ A_f^{\,lm} e_{3f} \zeta \right] 1148 \right] 1149 \right\} 1150 && \\ 1139 \delta_{i+1/2} \left[ \frac {e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ A_f^{\,lm} e_{3f} \zeta \right] \right] 1140 + \delta_{j+1/2} \left[ \frac {e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ A_f^{\,lm} e_{3f} \zeta \right] \right] 1141 \right\} \\ 1151 1142 % 1152 1143 \intertext{Using \eqref{DOM_di_adj}, it follows:} … … 1154 1145 \equiv& \sum\limits_{i,j,k} 1155 1146 -\,\left\{ 1156 \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i 1157 \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_i \left[ 1\right] 1158 + \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j 1159 \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_j \left[ 1\right] 1147 \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_i \left[ 1\right] 1148 + \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ A_f^{\,lm} e_{3f} \zeta \right]\;\delta_j \left[ 1\right] 1160 1149 \right\} \quad \equiv 0 1161 &&\\1150 \\ 1162 1151 \end{flalign*} 1163 1152 … … 1167 1156 \subsection{Dissipation of Horizontal Kinetic Energy} 1168 1157 \label{Apdx_C.3.2} 1169 1170 1158 1171 1159 The lateral momentum diffusion term dissipates the horizontal kinetic energy: … … 1221 1209 \label{Apdx_C.3.3} 1222 1210 1223 1224 1211 The lateral momentum diffusion term dissipates the enstrophy when the eddy 1225 1212 coefficients are horizontally uniform: … … 1228 1215 \left[ \nabla_h \left( A^{\,lm}\;\chi \right) 1229 1216 - \nabla_h \times \left( A^{\,lm}\;\zeta \; \textbf{k} \right) \right]\;dv &&&\\ 1230 & = A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times1217 &\quad = A^{\,lm} \int \limits_D \zeta \textbf{k} \cdot \nabla \times 1231 1218 \left[ \nabla_h \times \left( \zeta \; \textbf{k} \right) \right]\;dv &&&\\ 1232 &\ equiv A^{\,lm} \sum\limits_{i,j,k} \zeta \;e_{3f}1219 &\quad \equiv A^{\,lm} \sum\limits_{i,j,k} \zeta \;e_{3f} 1233 1220 \left\{ \delta_{i+1/2} \left[ \frac{e_{2v}} {e_{1v}\,e_{3v}} \delta_i \left[ e_{3f} \zeta \right] \right] 1234 1221 + \delta_{j+1/2} \left[ \frac{e_{1u}} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta \right] \right] \right\} &&&\\ … … 1236 1223 \intertext{Using \eqref{DOM_di_adj}, it follows:} 1237 1224 % 1238 &\ equiv - A^{\,lm} \sum\limits_{i,j,k}1225 &\quad \equiv - A^{\,lm} \sum\limits_{i,j,k} 1239 1226 \left\{ \left( \frac{1} {e_{1v}\,e_{3v}} \delta_i \left[ e_{3f} \zeta \right] \right)^2 b_v 1240 + \left( \frac{1} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta \right] \right)^2 b_u \right\} &&&\\ 1241 & \leq \;0 &&&\\ 1227 + \left( \frac{1} {e_{2u}\,e_{3u}} \delta_j \left[ e_{3f} \zeta \right] \right)^2 b_u \right\} \quad \leq \;0 &&&\\ 1242 1228 \end{flalign*} 1243 1229 … … 1250 1236 When the horizontal divergence of the horizontal diffusion of momentum 1251 1237 (discrete sense) is taken, the term associated with the vertical curl of the 1252 vorticity is zero locally, due to (!!! II.1.8 !!!!!). The resulting term conserves the1253 $\chi$ and dissipates $\chi^2$ when the eddy coefficients are1254 horizontally uniform.1238 vorticity is zero locally, due to \eqref{Eq_DOM_div_curl}. 1239 The resulting term conserves the $\chi$ and dissipates $\chi^2$ 1240 when the eddy coefficients are horizontally uniform. 1255 1241 \begin{flalign*} 1256 1242 & \int\limits_D \nabla_h \cdot 1257 1243 \Bigl[ \nabla_h \left( A^{\,lm}\;\chi \right) 1258 1244 - \nabla_h \times \left( A^{\,lm}\;\zeta \;\textbf{k} \right) \Bigr] dv 1259 = \int\limits_D \nabla_h \cdot \nabla_h \left( A^{\,lm}\;\chi \right) dv &&&\\1245 = \int\limits_D \nabla_h \cdot \nabla_h \left( A^{\,lm}\;\chi \right) dv \\ 1260 1246 % 1261 1247 &\equiv \sum\limits_{i,j,k} 1262 1248 \left\{ \delta_i \left[ A_u^{\,lm} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} \left[ \chi \right] \right] 1263 + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right] \right\} &&&\\1249 + \delta_j \left[ A_v^{\,lm} \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right] \right\} \\ 1264 1250 % 1265 1251 \intertext{Using \eqref{DOM_di_adj}, it follows:} … … 1267 1253 &\equiv \sum\limits_{i,j,k} 1268 1254 - \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] 1269 + \frac{e_{1v}\,e_{3v}} 1270 \q quad \equiv 0 &&&\\1255 + \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\} 1256 \quad \equiv 0 \\ 1271 1257 \end{flalign*} 1272 1258 … … 1281 1267 \left[ \nabla_h \left( A^{\,lm}\;\chi \right) 1282 1268 - \nabla_h \times \left( A^{\,lm}\;\zeta \;\textbf{k} \right) \right]\; dv 1283 = A^{\,lm} \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\; dv &&&\\1269 = A^{\,lm} \int\limits_D \chi \;\nabla_h \cdot \nabla_h \left( \chi \right)\; dv \\ 1284 1270 % 1285 1271 &\equiv A^{\,lm} \sum\limits_{i,j,k} \frac{1} {e_{1t}\,e_{2t}\,e_{3t}} \chi … … 1287 1273 \delta_i \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} \left[ \chi \right] \right] 1288 1274 + \delta_j \left[ \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right] 1289 \right\} \; e_{1t}\,e_{2t}\,e_{3t} &&&\\1275 \right\} \; e_{1t}\,e_{2t}\,e_{3t} \\ 1290 1276 % 1291 1277 \intertext{Using \eqref{DOM_di_adj}, it turns out to be:} … … 1293 1279 &\equiv - A^{\,lm} \sum\limits_{i,j,k} 1294 1280 \left\{ \left( \frac{1} {e_{1u}} \delta_{i+1/2} \left[ \chi \right] \right)^2 b_u 1295 + \left( \frac{1} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right)^2 b_v \right\} \; &&&\\ 1296 % 1297 &\leq 0 &&&\\ 1281 + \left( \frac{1} {e_{2v}} \delta_{j+1/2} \left[ \chi \right] \right)^2 b_v \right\} 1282 \quad \leq 0 \\ 1298 1283 \end{flalign*} 1299 1284 … … 1303 1288 \section{Conservation Properties on Vertical Momentum Physics} 1304 1289 \label{Apdx_C_4} 1305 1306 1290 1307 1291 As for the lateral momentum physics, the continuous form of the vertical diffusion … … 1319 1303 \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right)\; dv \quad &\leq 0 \\ 1320 1304 \end{align*} 1305 1321 1306 The first property is obvious. The second results from: 1322 1323 1307 \begin{flalign*} 1324 1308 \int\limits_D … … 1359 1343 e_{1f}\,e_{2f}\,e_{3f} \; \equiv 0 && \\ 1360 1344 \end{flalign*} 1345 1361 1346 If the vertical diffusion coefficient is uniform over the whole domain, the 1362 1347 enstrophy is dissipated, $i.e.$ … … 1366 1351 \left( \frac{A^{\,vm}} {e_3 }\; \frac{\partial \textbf{U}_h } {\partial k} \right) \right)\; dv = 0 &&&\\ 1367 1352 \end{flalign*} 1353 1368 1354 This property is only satisfied in $z$-coordinates: 1369 1370 1355 \begin{flalign*} 1371 1356 \int\limits_D \zeta \, \textbf{k} \cdot \nabla \times … … 1477 1462 1478 1463 The numerical schemes used for tracer subgridscale physics are written such 1479 that the heat and salt contents are conserved (equations in flux form, second 1480 order centered finite differences). Since a flux form is used to compute the 1481 temperature and salinity, the quadratic form of these quantities (i.e. their variance) 1482 globally tends to diminish. As for the advection term, there is generally no strict 1483 conservation of mass, even if in practice the mass is conserved to a very high 1484 accuracy. 1464 that the heat and salt contents are conserved (equations in flux form). 1465 Since a flux form is used to compute the temperature and salinity, 1466 the quadratic form of these quantities ($i.e.$ their variance) globally tends to diminish. 1467 As for the advection term, there is conservation of mass only if the Equation Of Seawater is linear. 1485 1468 1486 1469 % ------------------------------------------------------------------------------------------------------------- -
trunk/DOC/TexFiles/Chapters/Annex_D.tex
r3294 r6289 120 120 \hline 121 121 public \par or \par module variable& 122 \textbf{m n} \par \textit{but not} \par \textbf{nn\_ }&122 \textbf{m n} \par \textit{but not} \par \textbf{nn\_ np\_}& 123 123 \textbf{a b e f g h o q r} \par \textbf{t} \textit{to} \textbf{x} \par but not \par \textbf{fs rn\_}& 124 124 \textbf{l} \par \textit{but not} \par \textbf{lp ld} \par \textbf{ ll ln\_}& … … 156 156 \hline 157 157 parameter& 158 \textbf{jp }&158 \textbf{jp np\_}& 159 159 \textbf{pp}& 160 160 \textbf{lp}& … … 190 190 %-------------------------------------------------------------------------------------------------------------- 191 191 192 N.B. Parameter here, in not only parameter in the \textsc{Fortran} acceptation, it is also used for code variables 193 that are read in namelist and should never been modified during a simulation. 194 It is the case, for example, for the size of a domain (jpi,jpj,jpk). 195 192 196 \newpage 193 197 % ================================================================ -
trunk/DOC/TexFiles/Chapters/Annex_ISO.tex
r4147 r6289 11 11 \namdisplay{namtra_ldf} 12 12 %--------------------------------------------------------------------------------------------------------- 13 If the namelist variable \np{ln\_traldf\_grif} is set true (and 14 \key{ldfslp} is set), \NEMO updates both active and passive tracers 15 using the Griffies triad representation of iso-neutral diffusion and 16 the eddy-induced advective skew (GM) fluxes. Otherwise (by default) the 17 filtered version of Cox's original scheme is employed 18 (\S\ref{LDF_slp}). In the present implementation of the Griffies 19 scheme, the advective skew fluxes are implemented even if 20 \key{traldf\_eiv} is not set. 13 14 Two scheme are available to perform the iso-neutral diffusion. 15 If the namelist logical \np{ln\_traldf\_triad} is set true, 16 \NEMO updates both active and passive tracers using the Griffies triad representation 17 of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes. 18 If the namelist logical \np{ln\_traldf\_iso} is set true, 19 the filtered version of Cox's original scheme (the Standard scheme) is employed (\S\ref{LDF_slp}). 20 In the present implementation of the Griffies scheme, 21 the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false. 21 22 22 23 Values of iso-neutral diffusivity and GM coefficient are set as 23 described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd}, 24 N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and 25 GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and 26 \np{rn\_aeiv\_0}. If 2D-varying coefficients are set with 27 \key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal 28 scale factor according to \eqref{Eq_title} \footnote{Except in global ORCA 29 $0.5^{\circ}$ runs with \key{traldf\_eiv}, where 30 $A_l$ is set like $A_e$ but with a minimum vale of 31 $100\;\mathrm{m}^2\;\mathrm{s}^{-1}$}. In idealised setups with 32 \key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv} 33 is set in the global configurations with \key{traldf\_c2d}, a horizontally varying $A_e$ is 34 instead set from the Held-Larichev parameterisation\footnote{In this 35 case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further 36 reduced by a factor $|f/f_{20}|$, where $f_{20}$ is the value of $f$ 37 at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored 38 unless it is zero. 24 described in \S\ref{LDF_coef}. Note that when GM fluxes are used, 25 the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS, 26 even though the eddy advection is accomplished by means of the skew fluxes. 27 39 28 40 29 The options specific to the Griffies scheme include: 41 30 \begin{description}[font=\normalfont] 42 \item[\np{ln\_traldf\_gdia}] Default value is false. See \S\ref{sec:triad:sfdiag}. If this is set true, time-mean 43 eddy-advective (GM) velocities are output for diagnostic purposes, even 44 though the eddy advection is accomplished by means of the skew 45 fluxes. 46 \item[\np{ln\_traldf\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then 31 \item[\np{ln\_triad\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then 47 32 `iso-neutral' mixing is accomplished within the surface mixed-layer 48 33 along slopes linearly decreasing with depth from the value immediately below 49 the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}). This is the same 50 treatment as used in the default implementation 51 \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}. Where 52 \np{ln\_traldf\_iso} is set true, the vertical skew flux is further 53 reduced to ensure no vertical buoyancy flux, giving an almost pure 34 the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}). 35 This is the same treatment as used in the default implementation \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}. 36 Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced 37 to ensure no vertical buoyancy flux, giving an almost pure 54 38 horizontal diffusive tracer flux within the mixed layer. This is similar to 55 39 the tapering suggested by \citet{Gerdes1991}. See \S\ref{sec:triad:Gerdes-taper} 56 \item[\np{ln\_traldf\_botmix}] See \S\ref{sec:triad:iso_bdry}. If this 57 is set false (the default) then the lateral diffusive fluxes 58 associated with triads partly masked by topography are neglected. If 59 it is set true, however, then these lateral diffusive fluxes are 60 applied, giving smoother bottom tracer fields at the cost of 61 introducing diapycnal mixing. 40 \item[\np{ln\_botmix\_triad}] See \S\ref{sec:triad:iso_bdry}. 41 If this is set false (the default) then the lateral diffusive fluxes 42 associated with triads partly masked by topography are neglected. 43 If it is set true, however, then these lateral diffusive fluxes are applied, 44 giving smoother bottom tracer fields at the cost of introducing diapycnal mixing. 45 \item[\np{rn\_sw\_triad}] blah blah to be added.... 46 \end{description} 47 The options shared with the Standard scheme include: 48 \begin{description}[font=\normalfont] 49 \item[\np{ln\_traldf\_msc}] blah blah to be added 50 \item[\np{rn\_slpmax}] blah blah to be added 62 51 \end{description} 63 52 \section{Triad formulation of iso-neutral diffusion} 64 53 \label{sec:triad:iso} 65 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO66 framework, using scale factors rather than grid-sizes.54 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, 55 but formulated within the \NEMO framework, using scale factors rather than grid-sizes. 67 56 68 57 \subsection{The iso-neutral diffusion operator} … … 84 73 \mbox{with}\quad \;\;\Re = 85 74 \begin{pmatrix} 86 1&0&-r_1\mystrut \\87 0&1&-r_2\mystrut \\88 -r_1 &-r_2&r_1 ^2+r_2 ^2\mystrut75 1 & 0 & -r_1 \mystrut \\ 76 0 & 1 & -r_2 \mystrut \\ 77 -r_1 & -r_2 & r_1 ^2+r_2 ^2 \mystrut 89 78 \end{pmatrix} 90 79 \quad \text{and} \quad\grad T= 91 80 \begin{pmatrix} 92 \frac{1}{e_1} \pd[T]{i}\mystrut \\93 \frac{1}{e_2} \pd[T]{j}\mystrut \\94 \frac{1}{e_3} \pd[T]{k}\mystrut81 \frac{1}{e_1} \pd[T]{i} \mystrut \\ 82 \frac{1}{e_2} \pd[T]{j} \mystrut \\ 83 \frac{1}{e_3} \pd[T]{k} \mystrut 95 84 \end{pmatrix}. 96 85 \end{equation} … … 101 90 % {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 102 91 % \end{array} }} \right) 103 Here \eqref{Eq_PE_iso_slopes} 92 Here \eqref{Eq_PE_iso_slopes} 104 93 \begin{align*} 105 94 r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} … … 200 189 % the mean vertical gradient at the $u$-point, 201 190 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 202 \begin{figure}[ h] \begin{center}191 \begin{figure}[tb] \begin{center} 203 192 \includegraphics[width=1.05\textwidth]{./TexFiles/Figures/Fig_GRIFF_triad_fluxes} 204 193 \caption{ \label{fig:triad:ISO_triad} … … 256 245 \ 257 246 \frac 258 {\left(\alpha / \beta \right)_i^k \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] } 259 {\left(\alpha / \beta \right)_i^k \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }. 260 \end{equation} 261 In calculating the slopes of the local neutral 262 surfaces, the expansion coefficients $\alpha$ and $\beta$ are 263 evaluated at the anchor points of the triad \footnote{Note that in \eqref{eq:triad:R} we use the ratio $\alpha / \beta$ 264 instead of multiplying the temperature derivative by $\alpha$ and the 265 salinity derivative by $\beta$. This is more efficient as the ratio 266 $\alpha / \beta$ can to be evaluated directly}, while the metrics are 267 calculated at the $u$- and $w$-points on the arms. 247 { \alpha_i^k \ \delta_{i+i_p}[T^k] - \beta_i^k \ \delta_{i+i_p}[S^k] } 248 { \alpha_i^k \ \delta_{k+k_p}[T^i] - \beta_i^k \ \delta_{k+k_p}[S^i] }. 249 \end{equation} 250 In calculating the slopes of the local neutral surfaces, 251 the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad, 252 while the metrics are calculated at the $u$- and $w$-points on the arms. 268 253 269 254 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 270 \begin{figure}[ h] \begin{center}255 \begin{figure}[tb] \begin{center} 271 256 \includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_GRIFF_qcells} 272 257 \caption{ \label{fig:triad:qcells} … … 277 262 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 278 263 279 Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter 280 cell that is the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ 281 $u$-cell and the $i,k+k_p$ $w$-cell. Expressing the slopes $s_i$ and 282 $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation, we have 283 e.g.\ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. Each triad slope $_i^k 284 \mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to calculate the 285 lateral flux along its $u$-arm, at $(i+i_p,k)$, and then again as an 286 $s'$ to calculate the vertical flux along its $w$-arm at 287 $(i,k+k_p)$. Each vertical area $a_i$ used to calculate the lateral 288 flux and horizontal area $a'_i$ used to calculate the vertical flux 289 can also be identified as the area across the $u$- and $w$-arms of a 290 unique triad, and we notate these areas, similarly to the triad 291 slopes, as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, 292 $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq:triad:i13} 293 $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq:triad:i31} 294 $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 264 Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter 265 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. 266 Expressing the slopes $s_i$ and $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation, 267 we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. 268 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) 269 to calculate the lateral flux along its $u$-arm, at $(i+i_p,k)$, 270 and then again as an $s'$ to calculate the vertical flux along its $w$-arm at $(i,k+k_p)$. 271 Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used 272 to calculate the vertical flux can also be identified as the area across the $u$- and $w$-arms 273 of a unique triad, and we notate these areas, similarly to the triad slopes, 274 as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, 275 where $e.g.$ in \eqref{eq:triad:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, 276 and in \eqref{eq:triad:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 295 277 296 278 \subsection{The full triad fluxes} … … 667 649 or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is 668 650 masked. The associated lateral fluxes (grey-black dashed line) are 669 masked if \np{ln\_botmix\_ grif}=false, but left unmasked,670 giving bottom mixing, if \np{ln\_botmix\_ grif}=true.671 672 The default option \np{ln\_botmix\_ grif}=false is suitable when the651 masked if \np{ln\_botmix\_triad}=false, but left unmasked, 652 giving bottom mixing, if \np{ln\_botmix\_triad}=true. 653 654 The default option \np{ln\_botmix\_triad}=false is suitable when the 673 655 bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}=1), 674 656 or for simple idealized problems. For setups with topography without 675 bbl mixing, \np{ln\_botmix\_ grif}=true may be necessary.657 bbl mixing, \np{ln\_botmix\_triad}=true may be necessary. 676 658 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 677 659 \begin{figure}[h] \begin{center} … … 690 672 or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point 691 673 is masked. The associated lateral fluxes (grey-black dashed 692 line) are masked if \np{botmix\_ grif}=.false., but left693 unmasked, giving bottom mixing, if \np{botmix\_ grif}=.true.}674 line) are masked if \np{botmix\_triad}=.false., but left 675 unmasked, giving bottom mixing, if \np{botmix\_triad}=.true.} 694 676 \end{center} \end{figure} 695 677 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 931 913 it to the Eulerian velocity prior to computing the tracer 932 914 advection. This is implemented if \key{traldf\_eiv} is set in the 933 default implementation, where \np{ln\_traldf\_ grif} is set915 default implementation, where \np{ln\_traldf\_triad} is set 934 916 false. This allows us to take advantage of all the advection schemes 935 917 offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$ … … 938 920 paramount importance. 939 921 940 However, when \np{ln\_traldf\_ grif} is set true, \NEMO instead922 However, when \np{ln\_traldf\_triad} is set true, \NEMO instead 941 923 implements eddy induced advection according to the so-called skew form 942 924 \citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes … … 1137 1119 and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the 1138 1120 $i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ 1139 $u$-point is masked. The namelist parameter \np{ln\_botmix\_ grif} has1121 $u$-point is masked. The namelist parameter \np{ln\_botmix\_triad} has 1140 1122 no effect on the eddy-induced skew-fluxes. 1141 1123 -
trunk/DOC/TexFiles/Chapters/Chap_DIA.tex
r6140 r6289 12 12 % Old Model Output 13 13 % ================================================================ 14 \section{Old Model Output (default or \key{dimgout})}14 \section{Old Model Output (default)} 15 15 \label{DIA_io_old} 16 16 … … 33 33 "\textit{grep -i numout}" in the source code directory. 34 34 35 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}. 36 37 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. 38 39 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 IOIPSL-based 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 time-steps (namelist parameter). 35 By default, diagnostic output files are written in NetCDF format. 36 Since version 3.2, when defining \key{iomput}, an I/O server has been added 37 which provides more flexibility in the choice of the fields to be written 38 as well as how the writing work is distributed over the processors in massively parallel computing. 39 A complete description of the use of this I/O server is presented in the next section. 40 41 By default, \key{iomput} is not defined, NEMO produces NetCDF with the old IOIPSL library 42 which has been kept for compatibility and its easy installation. 43 However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2, 44 many diagnostic options have been added presuming the use of \key{iomput}. 45 The usefulness of the default IOIPSL-based option is expected to reduce with each new release. 46 If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module 47 and contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period 48 of nn\_write time-steps (namelist parameter). 40 49 41 50 %\gmcomment{ % start of gmcomment … … 48 57 49 58 50 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: 59 Since version 3.2, iomput is the NEMO output interface of choice. 60 It has been designed to be simple to use, flexible and efficient. 61 The two main purposes of iomput are: 51 62 \begin{enumerate} 52 \item The complete and flexible control of the output files through external XML files adapted by the user from standard templates. 53 \item To achieve high performance and scalable output through the optional distribution of all diagnostic output related tasks to dedicated processes. 63 \item The complete and flexible control of the output files through external XML files 64 adapted by the user from standard templates. 65 \item To achieve high performance and scalable output through the optional distribution 66 of all diagnostic output related tasks to dedicated processes. 54 67 \end{enumerate} 55 The first functionality allows the user to specify, without code changes or recompilation, aspects of the diagnostic output stream, such as: 68 The first functionality allows the user to specify, without code changes or recompilation, 69 aspects of the diagnostic output stream, such as: 56 70 \begin{itemize} 57 71 \item The choice of output frequencies that can be different for each file (including real months and years). 58 \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). 72 \item The choice of file contents; includes complete flexibility over which data are written in which files 73 (the same data can be written in different files). 59 74 \item The possibility to split output files at a chosen frequency. 60 75 \item The possibility to extract a vertical or an horizontal subdomain. … … 62 77 \item Control over metadata via a large XML "database" of possible output fields. 63 78 \end{itemize} 64 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: 79 In addition, iomput allows the user to add in the code the output of any new variable (scalar, 2D or 3D) 80 in a very easy way. All details of iomput functionalities are listed in the following subsections. 81 Examples of the XML files that control the outputs can be found in: 65 82 \begin{alltt} 66 83 \begin{verbatim} … … 72 89 \end{alltt} 73 90 74 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. 75 76 In version 3.6, the iom\_put interface depends on an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.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 $--$enable-parallel). Note that the files created by iomput through XIOS are incompatible with NetCDF3. All post-processsing and visualization tools must therefore be compatible with NetCDF4 and not only NetCDF3. 77 78 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 post-processing burden usually associated with using large numbers of NEMO processors. Note that for smaller configurations, the rebuilding phase can be avoided, even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server. 91 The second functionality targets output performance when running in parallel (\key{mpp\_mpi}). 92 Iomput provides the possibility to specify N dedicated I/O processes (in addition to the NEMO processes) 93 to collect and write the outputs. With an appropriate choice of N by the user, 94 the bottleneck associated with the writing of the output files can be greatly reduced. 95 96 In version 3.6, the iom\_put interface depends on an external code called 97 \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.0} 98 (use of revision 618 or higher is required). This new IO server can take advantage of 99 the parallel I/O functionality of NetCDF4 to create a single output file and therefore 100 to bypass the rebuilding phase. Note that writing in parallel into the same NetCDF files 101 requires that your NetCDF4 library is linked to an HDF5 library that has been correctly compiled 102 ($i.e.$ with the configure option $--$enable-parallel). 103 Note that the files created by iomput through XIOS are incompatible with NetCDF3. 104 All post-processsing and visualization tools must therefore be compatible with NetCDF4 and not only NetCDF3. 105 106 Even if not using the parallel I/O functionality of NetCDF4, using N dedicated I/O servers, 107 where N is typically much less than the number of NEMO processors, will reduce the number of output files created. 108 This can greatly reduce the post-processing burden usually associated with using large numbers of NEMO processors. 109 Note that for smaller configurations, the rebuilding phase can be avoided, 110 even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server. 79 111 80 112 \subsection{XIOS: the IO\_SERVER} … … 82 114 \subsubsection{Attached or detached mode?} 83 115 84 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''. 85 86 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. 87 88 In detached mode and if the type of file is ''multiple\_file'', then each stand-alone 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. 116 Iomput is based on \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS}, 117 the io\_server developed by Yann Meurdesoif from IPSL. 118 The behaviour of the I/O subsystem is controlled by settings in the external XML files listed above. 119 Key settings in the iodef.xml file are {\tt using\_server} and the {\tt type} tag associated with each defined file. 120 The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode} (as a library) 121 [{\tt false}] or in \textit{detached mode} (as an external executable on N additional, dedicated cpus) [{\tt true}]. 122 The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications. 123 The type of each file can be either ''multiple\_file'' or ''one\_file''. 124 125 In \textit{attached mode} and if the type of file is ''multiple\_file'', 126 then each NEMO process will also act as an IO server and produce its own set of output files. 127 Superficially, this emulates the standard behaviour in previous versions. 128 However, the subdomain written out by each process does not correspond to the {\tt jpi x jpj x jpk} 129 domain actually computed by the process (although it may if {\tt jpni=1}). 130 Instead each process will have collected and written out a number of complete longitudinal strips. 131 If the ''one\_file'' option is chosen then all processes will collect their longitudinal strips 132 and write (in parallel) to a single output file. 133 134 In \textit{detached mode} and if the type of file is ''multiple\_file'', 135 then each stand-alone XIOS process will collect data for a range of complete longitudinal strips 136 and write to its own set of output files. If the ''one\_file'' option is chosen then 137 all XIOS processes will collect their longitudinal strips and write (in parallel) to a single output file. 138 Note running in detached mode requires launching a Multiple Process Multiple Data (MPMD) parallel job. 139 The following subsection provides a typical example but the syntax will vary in different MPP environments. 89 140 90 141 \subsubsection{Number of cpu used by XIOS in detached mode} 91 142 92 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: 143 The number of cores used by the XIOS is specified when launching the model. 144 The number of cores dedicated to XIOS should be from ~1/10 to ~1/50 of the number or cores dedicated to NEMO. 145 Some manufacturers suggest using O($\sqrt{N}$) dedicated IO processors for N processors 146 but this is a general recommendation and not specific to NEMO. 147 It is difficult to provide precise recommendations because the optimal choice will depend on 148 the particular hardware properties of the target system (parallel filesystem performance, available memory, 149 memory bandwidth etc.) and the volume and frequency of data to be created. 150 Here is an example of 2 cpus for the io\_server and 62 cpu for nemo using mpirun: 93 151 94 152 \texttt{ mpirun -np 62 ./nemo.exe : -np 2 ./xios\_server.exe } … … 96 154 \subsubsection{Control of XIOS: the XIOS context in iodef.xml} 97 155 98 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. 156 As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml. 157 See the XML basics section below for more details on XML syntax and rules. 99 158 100 159 \begin{tabular}{|p{4cm}|p{6.0cm}|p{2.0cm}|} … … 106 165 \hline 107 166 buffer\_size & 108 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 & 167 buffer size used by XIOS to send data from NEMO to XIOS. Larger is more efficient. 168 Note that needed/used buffer sizes are summarized at the end of the job & 109 169 25000000 \\ 110 170 \hline … … 136 196 \subsubsection{Installation} 137 197 138 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 139 \href{http://www.nemo-ocean.eu/Using-NEMO/User-Guides/Basics/XIOS-IO-server-installation-and-use}{XIOS with NEMO} guide provides an example illustration of how this can be achieved. 198 As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with NEMO. 199 See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance. 200 NEMO will need to link to the compiled XIOS library. 201 The \href{http://www.nemo-ocean.eu/Using-NEMO/User-Guides/Basics/XIOS-IO-server-installation-and-use}{XIOS with NEMO} 202 guide provides an example illustration of how this can be achieved. 140 203 141 204 \subsubsection{Add your own outputs} 142 205 143 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. 206 It is very easy to add your own outputs with iomput. 207 Many standard fields and diagnostics are already prepared ($i.e.$, steps 1 to 3 below have been done) 208 and simply need to be activated by including the required output in a file definition in iodef.xml (step 4). 209 To add new output variables, all 4 of the following steps must be taken. 144 210 \begin{description} 145 211 \item[1.] in NEMO code, add a \\ … … 151 217 to the list of used modules in the upper part of your module. 152 218 153 \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: 219 \item[3.] in the field\_def.xml file, add the definition of your variable using the same identifier 220 you used in the f90 code (see subsequent sections for a details of the XML syntax and rules). 221 For example: 154 222 \vspace{-20pt} 155 223 \begin{alltt} {{\scriptsize … … 165 233 \end{verbatim} 166 234 }}\end{alltt} 167 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.: 235 Note your definition must be added to the field\_group whose reference grid is consistent 236 with the size of the array passed to iomput. 237 The grid\_ref attribute refers to definitions set in iodef.xml which, in turn, reference grids 238 and axes either defined in the code (iom\_set\_domain\_attr and iom\_set\_axis\_attr in iom.F90) 239 or defined in the domain\_def.xml file. $e.g.$: 168 240 \vspace{-20pt} 169 241 \begin{alltt} {{\scriptsize … … 173 245 }}\end{alltt} 174 246 Note, if your array is computed within the surface module each nn\_fsbc time\_step, 175 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) 247 add the field definition within the field\_group defined with the id ''SBC'': $<$field\_group id=''SBC''...$>$ 248 which has been defined with the correct frequency of operations (iom\_set\_field\_attr in iom.F90) 176 249 177 250 \item[4.] add your field in one of the output files defined in iodef.xml (again see subsequent sections for syntax and rules) \\ … … 201 274 \subsubsection{Structure of the xml file used in NEMO} 202 275 203 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): 276 The XML file used in XIOS is structured by 7 families of tags: context, axis, domain, grid, field, file and variable. 277 Each tag family has hierarchy of three flavors (except for context): 204 278 \\ 205 279 \begin{tabular}{|p{3.0cm}|p{4.5cm}|p{4.5cm}|} … … 225 299 \\ 226 300 227 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. 228 229 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):\\ 301 Each element may have several attributes. 302 Some attributes are mandatory, other are optional but have a default value and other are are completely optional. 303 Id is a special attribute used to identify an element or a group of elements. 304 It must be unique for a kind of element. 305 It is optional, but no reference to the corresponding element can be done if it is not defined. 306 307 The XML file is split into context tags that are used to isolate IO definition from different codes 308 or different parts of a code. No interference is possible between 2 different contexts. 309 Each context has its own calendar and an associated timestep. 310 In \NEMO, we used the following contexts (that can be defined in any order):\\ 230 311 \\ 231 312 \begin{tabular}{|p{3.0cm}|p{4.5cm}|p{4.5cm}|} … … 271 352 \\ 272 353 273 \noindent Each context tag related to NEMO (mother or child grids) is divided into 5 parts (that can be defined in any order):\\ 354 \noindent Each context tag related to NEMO (mother or child grids) is divided into 5 parts 355 (that can be defined in any order):\\ 274 356 \\ 275 357 \begin{tabular}{|p{3.0cm}|p{4.5cm}|p{4.5cm}|} … … 305 387 \subsubsection{Nesting XML files} 306 388 307 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: \\ 389 The XML file can be split in different parts to improve its readability and facilitate its use. 390 The inclusion of XML files into the main XML file can be done through the attribute src: \\ 308 391 {\scriptsize \verb? <context src="./nemo_def.xml" /> ?}\\ 309 392 … … 323 406 \subsubsection{Use of inheritance} 324 407 325 XML extensively uses the concept of inheritance. XML has a tree based structure with a parent-child 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. \\ 408 XML extensively uses the concept of inheritance. 409 XML has a tree based structure with a parent-child oriented relation: 410 all children inherit attributes from parent, but an attribute defined in a child replace the inherited attribute value. 411 Note that the special attribute ''id'' is never inherited. \\ 326 412 \\ 327 413 example 1: Direct inheritance. … … 362 448 \subsubsection{Use of Groups} 363 449 364 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''. 450 Groups can be used for 2 purposes. 451 Firstly, the group can be used to define common attributes to be shared by the elements of the group through inheritance. 452 In the following example, we define a group of field that will share a common grid ''grid\_T\_2D''. 453 Note that for the field ''toce'', we overwrite the grid definition inherited from the group by ''grid\_T\_3D''. 365 454 \vspace{-20pt} 366 455 \begin{alltt} {{\scriptsize … … 375 464 }}\end{alltt} 376 465 377 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: 466 Secondly, the group can be used to replace a list of elements. 467 Several examples of groups of fields are proposed at the end of the file {\tt CONFIG/SHARED/field\_def.xml}. 468 For example, a short list of the usual variables related to the U grid: 378 469 \vspace{-20pt} 379 470 \begin{alltt} {{\scriptsize … … 399 490 \subsection{Detailed functionalities } 400 491 401 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. 492 The file {\tt NEMOGCM/CONFIG/ORCA2\_LIM/iodef\_demo.xml} provides several examples of the use 493 of the new functionalities offered by the XML interface of XIOS. 402 494 403 495 \subsubsection{Define horizontal subdomains} 404 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). 496 Horizontal subdomains are defined through the attributs zoom\_ibegin, zoom\_jbegin, zoom\_ni, zoom\_nj of the tag family domain. 497 It must therefore be done in the domain part of the XML file. 498 For example, in {\tt CONFIG/SHARED/domain\_def.xml}, we provide the following example of a definition 499 of a 5 by 5 box with the bottom left corner at point (10,10). 405 500 \vspace{-20pt} 406 501 \begin{alltt} {{\scriptsize … … 419 514 \end{verbatim} 420 515 }}\end{alltt} 421 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'' ...) 516 Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain. 517 The Equatorial section, the TAO, RAMA and PIRATA moorings are alredy registered in the code 518 and can therefore be outputted without taking care of their (i,j) position in the grid. 519 These predefined domains can be activated by the use of specific domain\_ref: ''EqT'', ''EqU'' or ''EqW'' 520 for the equatorial sections and the mooring position for TAO, RAMA and PIRATA followed 521 by ''T'' (for example: ''8s137eT'', ''1.5s80.5eT'' ...) 422 522 \vspace{-20pt} 423 523 \begin{alltt} {{\scriptsize … … 1116 1216 % ------------------------------------------------------------------------------------------------------------- 1117 1217 \section[Tracer/Dynamics Trends (TRD)] 1118 {Tracer/Dynamics Trends (\key{trdtra}, \key{trddyn}, \\ 1119 \key{trddvor}, \key{trdmld})} 1218 {Tracer/Dynamics Trends (\ngn{namtrd})} 1120 1219 \label{DIA_trd} 1121 1220 … … 1124 1223 %------------------------------------------------------------------------------------------------------------- 1125 1224 1126 When \key{trddyn} and/or \key{trddyn} CPP variables are defined, each 1127 trend of the dynamics and/or temperature and salinity time evolution equations 1128 is stored in three-dimensional arrays just after their computation ($i.e.$ at the end 1129 of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). Options are defined by 1130 \ngn{namtrd} namelist variables. These trends are then 1131 used in \mdl{trdmod} (see TRD directory) every \textit{nn\_trd } time-steps. 1132 1133 What is done depends on the CPP keys defined: 1225 Each trend of the dynamics and/or temperature and salinity time evolution equations 1226 can be send to \mdl{trddyn} and/or \mdl{trdtra} modules (see TRD directory) just after their computation 1227 ($i.e.$ at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). 1228 This capability is controlled by options offered in \ngn{namtrd} namelist. 1229 Note that the output are done with xIOS, and therefore the \key{IOM} is required. 1230 1231 What is done depends on the \ngn{namtrd} logical set to \textit{true}: 1134 1232 \begin{description} 1135 \item[\key{trddyn}, \key{trdtra}] : a check of the basin averaged properties of the momentum 1136 and/or tracer equations is performed ; 1137 \item[\key{trdvor}] : a vertical summation of the moment tendencies is performed, 1138 then the curl is computed to obtain the barotropic vorticity tendencies which are output ; 1139 \item[\key{trdmld}] : output of the tracer tendencies averaged vertically 1140 either over the mixed layer (\np{nn\_ctls}=0), 1141 or over a fixed number of model levels (\np{nn\_ctls}$>$1 provides the number of level), 1142 or over a spatially varying but temporally fixed number of levels (typically the base 1143 of the winter mixed layer) read in \ifile{ctlsurf\_idx} (\np{nn\_ctls}=1) ; 1233 \item[\np{ln\_glo\_trd}] : at each \np{nn\_trd} time-step a check of the basin averaged properties 1234 of the momentum and tracer equations is performed. This also includes a check of $T^2$, $S^2$, 1235 $\tfrac{1}{2} (u^2+v2)$, and potential energy time evolution equations properties ; 1236 \item[\np{ln\_dyn\_trd}] : each 3D trend of the evolution of the two momentum components is output ; 1237 \item[\np{ln\_dyn\_mxl}] : each 3D trend of the evolution of the two momentum components averaged 1238 over the mixed layer is output ; 1239 \item[\np{ln\_vor\_trd}] : a vertical summation of the moment tendencies is performed, 1240 then the curl is computed to obtain the barotropic vorticity tendencies which are output ; 1241 \item[\np{ln\_KE\_trd}] : each 3D trend of the Kinetic Energy equation is output ; 1242 \item[\np{ln\_tra\_trd}] : each 3D trend of the evolution of temperature and salinity is output ; 1243 \item[\np{ln\_tra\_mxl}] : each 2D trend of the evolution of temperature and salinity averaged 1244 over the mixed layer is output ; 1144 1245 \end{description} 1145 1146 The units in the output file can be changed using the \np{nn\_ucf} namelist parameter.1147 For example, in case of salinity tendency the units are given by PSU/s/\np{nn\_ucf}.1148 Setting \np{nn\_ucf}=86400 ($i.e.$ the number of second in a day) provides the tendencies in PSU/d.1149 1150 When \key{trdmld} is defined, two time averaging procedure are proposed.1151 Setting \np{ln\_trdmld\_instant} to \textit{true}, a simple time averaging is performed,1152 so that the resulting tendency is the contribution to the change of a quantity between1153 the two instantaneous values taken at the extremities of the time averaging period.1154 Setting \np{ln\_trdmld\_instant} to \textit{false}, a double time averaging is performed,1155 so that the resulting tendency is the contribution to the change of a quantity between1156 two \textit{time mean} values. The later option requires the use of an extra file, \ifile{restart\_mld}1157 (\np{ln\_trdmld\_restart}=true), to restart a run.1158 1159 1246 1160 1247 Note that the mixed layer tendency diagnostic can also be used on biogeochemical models 1161 1248 via the \key{trdtrc} and \key{trdmld\_trc} CPP keys. 1249 1250 \textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested. 1251 In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} 1252 are not working, and none of the option have been tested with variable volume ($i.e.$ \key{vvl} defined). 1253 1162 1254 1163 1255 % ------------------------------------------------------------------------------------------------------------- … … 1280 1372 \label{DIA_diag_harm} 1281 1373 1282 A module is available to compute the amplitude and phase for tidal waves.1283 This diagnostic is actived with \key{diaharm}.1284 1285 1374 %------------------------------------------namdia_harm---------------------------------------------------- 1286 1375 \namdisplay{namdia_harm} 1287 1376 %---------------------------------------------------------------------------------------------------------- 1288 1377 1289 Concerning the on-line Harmonic analysis, some parameters are available in namelist 1290 \ngn{namdia\_harm} : 1291 1292 - \texttt{nit000\_han} is the first time step used for harmonic analysis 1293 1294 - \texttt{nitend\_han} is the last time step used for harmonic analysis 1295 1296 - \texttt{nstep\_han} is the time step frequency for harmonic analysis 1297 1298 - \texttt{nb\_ana} is the number of harmonics to analyse 1299 1300 - \texttt{tname} is an array with names of tidal constituents to analyse 1301 1302 \texttt{nit000\_han} and \texttt{nitend\_han} must be between \texttt{nit000} and \texttt{nitend} of the simulation. 1378 A module is available to compute the amplitude and phase of tidal waves. 1379 This on-line Harmonic analysis is actived with \key{diaharm}. 1380 Some parameters are available in namelist \ngn{namdia\_harm} : 1381 1382 - \np{nit000\_han} is the first time step used for harmonic analysis 1383 1384 - \np{nitend\_han} is the last time step used for harmonic analysis 1385 1386 - \np{nstep\_han} is the time step frequency for harmonic analysis 1387 1388 - \np{nb\_ana} is the number of harmonics to analyse 1389 1390 - \np{tname} is an array with names of tidal constituents to analyse 1391 1392 \np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation. 1303 1393 The restart capability is not implemented. 1304 1394 1305 The Harmonic analysis solve th isequation:1395 The Harmonic analysis solve the following equation: 1306 1396 \begin{equation} 1307 1397 h_{i} - A_{0} + \sum^{nb\_ana}_{j=1}[A_{j}cos(\nu_{j}t_{j}-\phi_{j})] = e_{i} … … 1350 1440 \label{DIA_diag_dct} 1351 1441 1352 A module is available to compute the transport of volume, heat and salt through sections. This diagnostic1353 is actived with \key{diadct}.1442 A module is available to compute the transport of volume, heat and salt through sections. 1443 This diagnostic is actived with \key{diadct}. 1354 1444 1355 1445 Each section is defined by the coordinates of its 2 extremities. The pathways between them are contructed … … 1373 1463 %------------------------------------------------------------------------------------------------------------- 1374 1464 1375 \ texttt{nn\_dct}: frequency of instantaneous transports computing1376 1377 \ texttt{nn\_dctwri}: frequency of writing ( mean of instantaneous transports )1378 1379 \ texttt{nn\_debug}: debugging of the section1465 \np{nn\_dct}: frequency of instantaneous transports computing 1466 1467 \np{nn\_dctwri}: frequency of writing ( mean of instantaneous transports ) 1468 1469 \np{nn\_debug}: debugging of the section 1380 1470 1381 1471 \subsubsection{ To create a binary file containing the pathway of each section } … … 1508 1598 the \key{diahth} CPP key: 1509 1599 1510 - the mixed layer depth (based on a density criterion , \citet{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})1600 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) 1511 1601 1512 1602 - the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) -
trunk/DOC/TexFiles/Chapters/Chap_DIU.tex
r6140 r6289 13 13 14 14 Code to produce an estimate of the diurnal warming and cooling of the sea surface skin 15 temperature (skin SST) is found in the DIU directory. The skin16 temperature can be split into three parts:15 temperature (skin SST) is found in the DIU directory. 16 The skin temperature can be split into three parts: 17 17 \begin{itemize} 18 18 \item A foundation SST which is free from diurnal warming. … … 25 25 Models are provided for both the warm layer, diurnal\_bulk.F90, and the cool skin, 26 26 cool\_skin.F90. Foundation SST is not considered as it can be obtained 27 either from the main NEMO model (i.e. from the temperature of the top few model levels) 28 or from 29 some other source. It must be noted that both the cool skin and 30 warm layer models produce estimates of the change in temperature ($\Delta T_{\rm{cs}}$ 31 and $\Delta T_{\rm{wl}}$) and both must 32 be added to a foundation SST to obtain the true skin temperature. 27 either from the main NEMO model ($i.e.$ from the temperature of the top few model levels) 28 or from some other source. 29 It must be noted that both the cool skin and warm layer models produce estimates of 30 the change in temperature ($\Delta T_{\rm{cs}}$ and $\Delta T_{\rm{wl}}$) 31 and both must be added to a foundation SST to obtain the true skin temperature. 33 32 34 Both the cool skin and warm layer models are controlled through the namelist `namdiu':33 Both the cool skin and warm layer models are controlled through the namelist \ngn{namdiu}: 35 34 \namdisplay{namdiu} 36 35 This namelist contains only two variables: 37 36 \begin{description} 38 \item[ ln\_diurnal] A logical switch for turning on/off both the cool skin and warm layer.39 \item[ ln\_diurnal\_only] A logical switch which if .TRUE. will run the diurnal model40 without the other dynamical parts of NEMO. ln\_diurnal\_only must be41 .FALSE. if ln\_diurnalis .FALSE.37 \item[\np{ln\_diurnal}] A logical switch for turning on/off both the cool skin and warm layer. 38 \item[\np{ln\_diurnal\_only}] A logical switch which if .TRUE. will run the diurnal model 39 without the other dynamical parts of NEMO. 40 \np{ln\_diurnal\_only} must be .FALSE. if \np{ln\_diurnal} is .FALSE. 42 41 \end{description} 43 42 … … 46 45 output if they are specified in the iodef.xml file. 47 46 48 Initialisation is through the restart file. Specifically the code will expect the presence of the 2-D variable ``Dsst'' to initialise the warm layer. The cool skin model, which is determined purely by the instantaneous fluxes, has no initialisation variable. 47 Initialisation is through the restart file. Specifically the code will expect 48 the presence of the 2-D variable ``Dsst'' to initialise the warm layer. 49 The cool skin model, which is determined purely by the instantaneous fluxes, 50 has no initialisation variable. 49 51 50 52 %=============================================================== 51 52 53 \section{Warm Layer model} 53 54 \label{warm_layer_sec} 54 55 55 %=============================================================== 56 56 … … 81 81 (\ref{ecmwf1}) is the instantaneous total thermal energy 82 82 flux into 83 the diurnal layer, i.e.83 the diurnal layer, $i.e.$ 84 84 \begin{equation} 85 85 Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,} \label{e_flux_eqn} … … 105 105 where $\zeta=\frac{D_T}{L}$. It is clear that the first derivative of 106 106 (\ref{stab_func_eqn}), and thus of (\ref{ecmwf1}), 107 is discontinuous at $\zeta=0$ ( i.e.$Q\rightarrow0$ in equation (\ref{ecmwf2})).107 is discontinuous at $\zeta=0$ ($i.e.$ $Q\rightarrow0$ in equation (\ref{ecmwf2})). 108 108 109 109 The two terms on the right hand side of (\ref{ecmwf1}) represent different processes. -
trunk/DOC/TexFiles/Chapters/Chap_DOM.tex
r6140 r6289 138 138 and $f$-points, and its divergence defined at $t$-points: 139 139 \begin{eqnarray} \label{Eq_DOM_curl} 140 \nabla \times {\rm 140 \nabla \times {\rm{\bf A}}\equiv & 141 141 \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} \\ 142 142 +& \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} \\ 143 143 +& \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} 144 144 \end{eqnarray} 145 \begin{equation} \label{Eq_DOM_div} 146 \nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}\,e_{3t}}\left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] 147 +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] 148 \end{equation} 149 150 In the special case of a pure $z$-coordinate system, \eqref{Eq_DOM_lap} and 151 \eqref{Eq_DOM_div} can be simplified. In this case, the vertical scale factor 152 becomes a function of the single variable $k$ and thus does not depend on the 153 horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to: 154 \begin{equation*} 155 \nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}} \left( \delta_i \left[e_{2u}\,a_1 \right] 156 +\delta_j \left[e_{1v}\, a_2 \right] \right) 157 +\frac{1}{e_{3t}} \delta_k \left[ a_3 \right] 158 \end{equation*} 145 \begin{eqnarray} \label{Eq_DOM_div} 146 \nabla \cdot \rm{\bf A} \equiv 147 \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] 148 +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] 149 \end{eqnarray} 159 150 160 151 The vertical average over the whole water column denoted by an overbar becomes … … 183 174 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside 184 175 continental area. Using integration by parts it can be shown that the differencing 185 operators ($\delta_i$, $\delta_j$ and $\delta_k$) are anti-symmetric linear186 operators,and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$,176 operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, 177 and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$, 187 178 $\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear 188 179 operators, $i.e.$ … … 425 416 426 417 The choice of the grid must be consistent with the boundary conditions specified 427 by the parameter \np{jperio}(see {\S\ref{LBC}).418 by \np{jperio}, a parameter found in \ngn{namcfg} namelist (see {\S\ref{LBC}). 428 419 429 420 % ------------------------------------------------------------------------------------------------------------- … … 505 496 If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft 506 497 (in meters) (\ifile{isf\_draft\_meter}) is needed and all the location where the isf cavity thinnest 507 than \np{rn\_isfhmin} meters are grounded ( iemasked).498 than \np{rn\_isfhmin} meters are grounded ($i.e.$ masked). 508 499 509 500 After reading the bathymetry, the algorithm for vertical grid definition differs … … 540 531 541 532 Three options are possible for defining the bathymetry, according to the 542 namelist variable \np{nn\_bathy} :533 namelist variable \np{nn\_bathy} (found in \ngn{namdom} namelist): 543 534 \begin{description} 544 535 \item[\np{nn\_bathy} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$ … … 721 712 usually 10\%, of the default thickness $e_{3t}(jk)$). 722 713 723 \colorbox{yellow}{Add a figure here of pstep especially at last ocean level}714 \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } } 724 715 725 716 % ------------------------------------------------------------------------------------------------------------- … … 749 740 depth, since a mixed step-like and bottom-following representation of the 750 741 topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f). 751 The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 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 terrain-following vertical coordinate is calculated. 752 753 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. 754 755 The original default NEMO s-coordinate 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}: 742 The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects 743 the sea bed and becomes a pseudo z-coordinate. 744 The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} 745 as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 746 747 Options for stretching the coordinate are provided as examples, but care must be taken to ensure 748 that the vertical stretch used is appropriate for the application. 749 750 The original default NEMO s-coordinate stretching is available if neither of the other options 751 are specified as true (\np{ln\_s\_SH94}~=~false and \np{ln\_s\_SF12}~=~false). 752 This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: 756 753 757 754 \begin{equation} … … 760 757 \end{equation} 761 758 762 where $s_{min}$ is the depth at which the s-coordinate stretching starts and allows a z-coordinate to placed on top of the stretched coordinate, and z is the depth (negative down from the asea surface). 759 where $s_{min}$ is the depth at which the s-coordinate stretching starts and 760 allows a z-coordinate to placed on top of the stretched coordinate, 761 and z is the depth (negative down from the asea surface). 763 762 764 763 \begin{equation} … … 775 774 \end{equation} 776 775 777 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: 776 A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} 777 stretching (\np{ln\_s\_SH94}~=~true), is also available and is more commonly used for shelf seas modelling: 778 778 779 779 \begin{equation} … … 792 792 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 793 793 794 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 795 bottom control parameters such that $0\leqslant \theta \leqslant 20$, and 794 where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from 795 pure $\sigma$ to the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) 796 are the surface and bottom control parameters such that $0\leqslant \theta \leqslant 20$, and 796 797 $0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom 797 798 increase of the vertical resolution (Fig.~\ref{Fig_sco_function}). 798 799 799 Another example has been provided at version 3.5 (\np{ln\_sco\_SF12}) that allows a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. In this case the a stretching function $\gamma$ is defined such that: 800 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows 801 a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. 802 In this case the a stretching function $\gamma$ is defined such that: 800 803 801 804 \begin{equation} … … 815 818 \end{equation} 816 819 817 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: 820 This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of 821 the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards 822 the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs}) 823 and bottom depths. The bottom cell depth in this example is given as a function of water depth: 818 824 819 825 \begin{equation} \label{DOM_zb} … … 843 849 \label{DOM_zgr_star} 844 850 845 This option is described in the Report by Levier \textit{et al.} (2007), available on 846 the \NEMO web site. 851 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 847 852 848 853 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances -
trunk/DOC/TexFiles/Chapters/Chap_DYN.tex
r6140 r6289 301 301 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 302 302 303 Note that a key point in \eqref{Eq_een_e3f} is that the averaging in the \textbf{i}- and 304 \textbf{j}- directions uses the masked vertical scale factor but is always divided by 305 $4$, not by the sum of the masks at the four $T$-points. This preserves the continuity of 306 $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and 307 extends by continuity the value of $e_{3f}$ into the land areas. This feature is essential for 308 the $z$-coordinate with partial steps. 303 A key point in \eqref{Eq_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. 304 It uses the sum of masked t-point vertical scale factor divided either 305 by the sum of the four t-point masks (\np{nn\_een\_e3f}~=~1), 306 or just by $4$ (\np{nn\_een\_e3f}~=~true). 307 The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ 308 tends to zero and extends by continuity the value of $e_{3f}$ into the land areas. 309 This case introduces a sub-grid-scale topography at f-points (with a systematic reduction of $e_{3f}$ 310 when a model level intercept the bathymetry) that tends to reinforce the topostrophy of the flow 311 ($i.e.$ the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}. 309 312 310 313 Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as … … 372 375 \end{aligned} \right. 373 376 \end{equation} 377 When \np{ln\_dynzad\_zts}~=~\textit{true}, a split-explicit time stepping with 5 sub-timesteps is used 378 on the vertical advection term. 379 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}. 380 Note that in this case, a similar split-explicit time stepping should be used on 381 vertical advection of tracer to ensure a better stability, 382 an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \S\ref{TRA_adv_tvd}). 383 374 384 375 385 % ================================================================ … … 716 726 $\ $\newline %force an empty line 717 727 718 %%%719 728 Options are defined through the \ngn{namdyn\_spg} namelist variables. 720 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 time-dependent 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 three-dimensional equations: the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}), and the split-explicit 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}. 721 722 %%% 729 The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}). 730 The main distinction is between the fixed volume case (linear free surface) and the variable volume case 731 (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface}) 732 the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case 733 (\S\ref{PE_free_surface}). 734 With both linear and nonlinear free surface, external gravity waves are allowed in the equations, 735 which imposes a very small time step when an explicit time stepping is used. 736 Two methods are proposed to allow a longer time step for the three-dimensional equations: 737 the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}), 738 and the split-explicit free surface described below. 739 The extra term introduced in the filtered method is calculated implicitly, 740 so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 723 741 724 742 … … 734 752 implicitly, so that a solver is used to compute it. As a consequence the update of the $next$ 735 753 velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 736 737 754 738 755 … … 777 794 $\rdt_e = \rdt / nn\_baro$. This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}=true) 778 795 considering that the stability of the barotropic system is essentially controled by external waves propagation. 779 Maximum allowed Courant number is in that case time independent, and easily computed online from the input bathymetry. 796 Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry. 797 Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn\_bt\_cmax}. 780 798 781 799 %%% … … 800 818 Schematic of the split-explicit time stepping scheme for the external 801 819 and internal modes. Time increases to the right. In this particular exemple, 802 a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_f ilt=1$) and $nn\_baro=5$.820 a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_flt=1$) and $nn\_baro=5$. 803 821 Internal mode time steps (which are also the model time steps) are denoted 804 822 by $t-\rdt$, $t$ and $t+\rdt$. Variables with $k$ superscript refer to instantaneous barotropic variables, … … 806 824 The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged 807 825 transports to advect tracers. 808 a) Forward time integration: \np{ln\_bt\_fw}=true, \np{ln\_bt\_ave}=true.809 b) Centred time integration: \np{ln\_bt\_fw}=false, \np{ln\_bt\_av e}=true.810 c) Forward time integration with no time filtering (POM-like scheme): \np{ln\_bt\_fw}=true, \np{ln\_bt\_av e}=false. }826 a) Forward time integration: \np{ln\_bt\_fw}=true, \np{ln\_bt\_av}=true. 827 b) Centred time integration: \np{ln\_bt\_fw}=false, \np{ln\_bt\_av}=true. 828 c) Forward time integration with no time filtering (POM-like scheme): \np{ln\_bt\_fw}=true, \np{ln\_bt\_av}=false. } 811 829 \end{center} \end{figure} 812 830 %> > > > > > > > > > > > > > > > > > > > > > > > > > > > … … 814 832 In the default case (\np{ln\_bt\_fw}=true), the external mode is integrated 815 833 between \textit{now} and \textit{after} baroclinic time-steps (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 816 quantities (\np{ln\_bt\_av e}=true). In that case, the integration is extended slightly beyond \textit{after} time step to provide time filtered quantities.834 quantities (\np{ln\_bt\_av}=true). In that case, the integration is extended slightly beyond \textit{after} time step to provide time filtered quantities. 817 835 These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. 818 836 Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, … … 835 853 %%% 836 854 837 One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_av e}=false).855 One can eventually choose to feedback instantaneous values by not using any time filter (\np{ln\_bt\_av}=false). 838 856 In that case, external mode equations are continuous in time, ie they are not re-initialized when starting a new 839 857 sub-stepping sequence. This is the method used so far in the POM model, the stability being maintained by refreshing at (almost) … … 987 1005 At the lateral boundaries either free slip, no slip or partial slip boundary 988 1006 conditions are applied according to the user's choice (see Chap.\ref{LBC}). 1007 1008 \gmcomment{ 1009 Hyperviscous operators are frequently used in the simulation of turbulent flows to control 1010 the dissipation of unresolved small scale features. 1011 Their primary role is to provide strong dissipation at the smallest scale supported by the grid 1012 while minimizing the impact on the larger scale features. 1013 Hyperviscous operators are thus designed to be more scale selective than the traditional, 1014 physically motivated Laplace operator. 1015 In finite difference methods, the biharmonic operator is frequently the method of choice to achieve 1016 this scale selective dissipation since its damping time ($i.e.$ its spin down time) 1017 scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$ 1018 (so that short waves damped more rapidelly than long ones), 1019 whereas the Laplace operator damping time scales only like $\lambda^{-2}$. 1020 } 989 1021 990 1022 % ================================================================ … … 1156 1188 1157 1189 Besides the surface and bottom stresses (see the above section) which are 1158 introduced as boundary conditions on the vertical mixing, two other forcings 1159 enter the dynamical equations. 1160 1161 One is the effect of atmospheric pressure on the ocean dynamics. 1162 Another forcing term is the tidal potential. 1163 Both of which will be introduced into the reference version soon. 1164 1165 \gmcomment{atmospheric pressure is there!!!! include its description } 1190 introduced as boundary conditions on the vertical mixing, three other forcings 1191 may enter the dynamical equations by affecting the surface pressure gradient. 1192 1193 (1) When \np{ln\_apr\_dyn}~=~true (see \S\ref{SBC_apr}), the atmospheric pressure is taken 1194 into account when computing the surface pressure gradient. 1195 1196 (2) When \np{ln\_tide\_pot}~=~true and \key{tide} is defined (see \S\ref{SBC_tide}), 1197 the tidal potential is taken into account when computing the surface pressure gradient. 1198 1199 (3) When \np{nn\_ice\_embd}~=~2 and LIM or CICE is used ($i.e.$ when the sea-ice is embedded in the ocean), 1200 the snow-ice mass is taken into account when computing the surface pressure gradient. 1201 1202 1203 \gmcomment{ missing : the lateral boundary condition !!! another external forcing 1204 } 1166 1205 1167 1206 % ================================================================ … … 1211 1250 1212 1251 % ================================================================ 1213 % Neptune effect1214 % ================================================================1215 \section [Neptune effect (\textit{dynnept})]1216 {Neptune effect (\mdl{dynnept})}1217 \label{DYN_nept}1218 1219 The "Neptune effect" (thus named in \citep{HollowayOM86}) is a1220 parameterisation of the potentially large effect of topographic form stress1221 (caused by eddies) in driving the ocean circulation. Originally developed for1222 low-resolution models, in which it was applied via a Laplacian (second-order)1223 diffusion-like term in the momentum equation, it can also be applied in eddy1224 permitting or resolving models, in which a more scale-selective bilaplacian1225 (fourth-order) implementation is preferred. This mechanism has a1226 significant effect on boundary currents (including undercurrents), and the1227 upwelling of deep water near continental shelves.1228 1229 The theoretical basis for the method can be found in1230 \citep{HollowayJPO92}, including the explanation of why form stress is not1231 necessarily a drag force, but may actually drive the flow.1232 \citep{HollowayJPO94} demonstrate the effects of the parameterisation in1233 the GFDL-MOM model, at a horizontal resolution of about 1.8 degrees.1234 \citep{HollowayOM08} demonstrate the biharmonic version of the1235 parameterisation in a global run of the POP model, with an average horizontal1236 grid spacing of about 32km.1237 1238 The NEMO implementation is a simplified form of that supplied by1239 Greg Holloway, the testing of which was described in \citep{HollowayJGR09}.1240 The major simplification is that a time invariant Neptune velocity1241 field is assumed. This is computed only once, during start-up, and1242 made available to the rest of the code via a module. Vertical1243 diffusive terms are also ignored, and the model topography itself1244 is used, rather than a separate topographic dataset as in1245 \citep{HollowayOM08}. This implementation is only in the iso-level1246 formulation, as is the case anyway for the bilaplacian operator.1247 1248 The velocity field is derived from a transport stream function given by:1249 1250 \begin{equation} \label{Eq_dynnept_sf}1251 \psi = -fL^2H1252 \end{equation}1253 1254 where $L$ is a latitude-dependant length scale given by:1255 1256 \begin{equation} \label{Eq_dynnept_ls}1257 L = l_1 + (l_2 -l_1)\left ( {1 + \cos 2\phi \over 2 } \right )1258 \end{equation}1259 1260 where $\phi$ is latitude and $l_1$ and $l_2$ are polar and equatorial length scales respectively.1261 Neptune velocity components, $u^*$, $v^*$ are derived from the stremfunction as:1262 1263 \begin{equation} \label{Eq_dynnept_vel}1264 u^* = -{1\over H} {\partial \psi \over \partial y}\ \ \ ,\ \ \ v^* = {1\over H} {\partial \psi \over \partial x}1265 \end{equation}1266 1267 \smallskip1268 %----------------------------------------------namdom----------------------------------------------------1269 \namdisplay{namdyn_nept}1270 %--------------------------------------------------------------------------------------------------------1271 \smallskip1272 1273 The Neptune effect is enabled when \np{ln\_neptsimp}=true (default=false).1274 \np{ln\_smooth\_neptvel} controls whether a scale-selective smoothing is applied1275 to the Neptune effect flow field (default=false) (this smoothing method is as1276 used by Holloway). \np{rn\_tslse} and \np{rn\_tslsp} are the equatorial and1277 polar values respectively of the length-scale parameter $L$ used in determining1278 the Neptune stream function \eqref{Eq_dynnept_sf} and \eqref{Eq_dynnept_ls}.1279 Values at intermediate latitudes are given by a cosine fit, mimicking the1280 variation of the deformation radius with latitude. The default values of 12km1281 and 3km are those given in \citep{HollowayJPO94}, appropriate for a coarse1282 resolution model. The finer resolution study of \citep{HollowayOM08} increased1283 the values of L by a factor of $\sqrt 2$ to 17km and 4.2km, thus doubling the1284 stream function for a given topography.1285 1286 The simple formulation for ($u^*$, $v^*$) can give unacceptably large velocities1287 in shallow water, and \citep{HollowayOM08} add an offset to the depth in the1288 denominator to control this problem. In this implementation we offer instead (at1289 the suggestion of G. Madec) the option of ramping down the Neptune flow field to1290 zero over a finite depth range. The switch \np{ln\_neptramp} activates this1291 option (default=false), in which case velocities at depths greater than1292 \np{rn\_htrmax} are unaltered, but ramp down linearly with depth to zero at a1293 depth of \np{rn\_htrmin} (and shallower).1294 1295 % ================================================================ -
trunk/DOC/TexFiles/Chapters/Chap_LBC.tex
r4147 r6289 1 1 % ================================================================ 2 % Chapter �Lateral Boundary Condition (LBC)2 % Chapter — Lateral Boundary Condition (LBC) 3 3 % ================================================================ 4 4 \chapter{Lateral Boundary Condition (LBC) } … … 127 127 can be quite shallow. 128 128 129 The alternative numerical implementation of the no-slip boundary conditions for an130 arbitrary coast line of \citet{Shchepetkin1996} is also available through the131 \key{noslip\_accurate} CPP key. It is based on a fourth order evaluation of the shear at the132 coast which, in turn, allows a true second order scheme in the interior of the domain133 ($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical134 scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a135 technique considerably improves the quality of the numerical solution. In \NEMO, such136 spectacular improvements have not been found in the half-degree global ocean137 (ORCA05), but significant reductions of numerically induced coastal upwellings were138 found in an eddy resolving simulation of the Alboran Sea \citep{Olivier_PhD01}.139 Nevertheless, since a no-slip boundary condition is not recommended in an eddy140 permitting or resolving simulation \citep{Penduff_al_OS07}, the use of this option is also141 not recommended.142 143 In practice, the no-slip accurate option changes the way the curl is evaluated at the144 coast (see \mdl{divcur} module), and requires the nature of each coastline grid point145 (convex or concave corners, straight north-south or east-west coast) to be specified.146 This is performed in routine \rou{dom\_msk\_nsa} in the \mdl{domask} module.147 129 148 130 % ================================================================ … … 204 186 % North fold (\textit{jperio = 3 }to $6)$ 205 187 % ------------------------------------------------------------------------------------------------------------- 206 \subsection{North-fold (\textit{jperio = 3 }to $6 )$}188 \subsection{North-fold (\textit{jperio = 3 }to $6$) } 207 189 \label{LBC_north_fold} 208 190 209 191 The north fold boundary condition has been introduced in order to handle the north 210 boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere. 211 \colorbox{yellow}{to be completed...} 192 boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere 193 (Fig.\ref{Fig_MISC_ORCA_msh}, and thus requires a specific treatment illustrated in Fig.\ref{Fig_North_Fold_T}. 194 Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition. 212 195 213 196 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 250 233 ocean model. Second order finite difference schemes lead to local discrete 251 234 operators that depend at the very most on one neighbouring point. The only 252 non-local computations concern the vertical physics (implicit diffusion, 1.5235 non-local computations concern the vertical physics (implicit diffusion, 253 236 turbulent closure scheme, ...) (delocalization over the whole water column), 254 237 and the solving of the elliptic equation associated with the surface pressure 255 238 gradient computation (delocalization over the whole horizontal domain). 256 239 Therefore, a pencil strategy is used for the data sub-structuration 257 \gmcomment{no idea what this means!}258 240 : the 3D initial domain is laid out on local processor 259 241 memories following a 2D horizontal topological splitting. Each sub-domain … … 264 246 phase starts: each processor sends to its neighbouring processors the update 265 247 values of the points corresponding to the interior overlapping area to its 266 neighbouring sub-domain (i.e. the innermost of the two overlapping rows). 267 The communication is done through message passing. Usually the parallel virtual 268 language, PVM, is used as it is a standard language available on nearly all 269 MPP computers. More specific languages (i.e. computer dependant languages) 270 can be easily used to speed up the communication, such as SHEM on a T3E 271 computer. The data exchanges between processors are required at the very 248 neighbouring sub-domain ($i.e.$ the innermost of the two overlapping rows). 249 The communication is done through the Message Passing Interface (MPI). 250 The data exchanges between processors are required at the very 272 251 place where lateral domain boundary conditions are set in the mono-domain 273 computation (\S III.10-c): the lbc\_lnk routine which manages such conditions 274 is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP 275 computer (\key{mpp\_mpi} defined). It has to be pointed out that when using 276 the MPP version of the model, the east-west cyclic boundary condition is done 277 implicitly, whilst the south-symmetric boundary condition option is not available. 252 computation : the \rou{lbc\_lnk} routine (found in \mdl{lbclnk} module) 253 which manages such conditions is interfaced with routines found in \mdl{lib\_mpp} module 254 when running on an MPP computer ($i.e.$ when \key{mpp\_mpi} defined). 255 It has to be pointed out that when using the MPP version of the model, 256 the east-west cyclic boundary condition is done implicitly, 257 whilst the south-symmetric boundary condition option is not available. 278 258 279 259 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 285 265 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 286 266 287 In the standard version of the OPA model, the splitting is regular and arithmetic. 288 the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors 289 \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in 290 \mdl{par\_oce}). Each processor is independent and without message passing 291 or synchronous process 292 \gmcomment{how does a synchronous process relate to this?}, 293 programs run alone and access just its own local memory. For this reason, the 294 main model dimensions are now the local dimensions of the subdomain (pencil) 267 In the standard version of \NEMO, the splitting is regular and arithmetic. 268 The i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors 269 \jp{jpnij} most often equal to $jpni \times jpnj$ (parameters set in 270 \ngn{nammpp} namelist). Each processor is independent and without message passing 271 or synchronous process, programs run alone and access just its own local memory. 272 For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil) 295 273 that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal 296 274 domain and the overlapping rows. The number of rows to exchange (known as … … 304 282 where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis. 305 283 306 \colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and 307 no east-west cyclic boundary conditions.} 308 309 One also defines variables nldi and nlei which correspond to the internal 310 domain bounds, and the variables nimpp and njmpp which are the position 311 of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array 312 (subdomain) corresponds to an element of $T_{g}$, a global array 313 (whole domain) by the relationship: 284 One also defines variables nldi and nlei which correspond to the internal domain bounds, 285 and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain. 286 An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$, 287 a global array (whole domain) by the relationship: 314 288 \begin{equation} \label{Eq_lbc_nimpp} 315 289 T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k), … … 320 294 nproc. In the standard version, a processor has no more than four neighbouring 321 295 processors named nono (for north), noea (east), noso (south) and nowe (west) 322 and two variables, nbondi and nbondj, indicate the relative position of the processor 323 \colorbox{yellow}{(see Fig.IV.3)}: 296 and two variables, nbondi and nbondj, indicate the relative position of the processor : 324 297 \begin{itemize} 325 298 \item nbondi = -1 an east neighbour, no west processor, … … 332 305 processor on its overlapping row, and sends the data issued from internal 333 306 domain corresponding to the overlapping row of the other processor. 334 335 \colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos }336 307 337 308 … … 343 314 global ocean where more than 50 \% of points are land points. For this reason, a 344 315 pre-processing tool can be used to choose the mpp domain decomposition with a 345 maximum number of only land points processors, which can then be eliminated .346 (For example, the mpp\_optimiz tools, available from the DRAKKAR web site .)316 maximum number of only land points processors, which can then be eliminated (Fig. \ref{Fig_mppini2}) 317 (For example, the mpp\_optimiz tools, available from the DRAKKAR web site). 347 318 This optimisation is dependent on the specific bathymetry employed. The user 348 319 then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with 349 320 $jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ 350 land processors. When those parameters are specified in module \mdl{par\_oce},321 land processors. When those parameters are specified in \ngn{nammpp} namelist, 351 322 the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound, 352 323 nono, noea,...) so that the land-only processors are not taken into account. 353 324 354 \ colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp325 \gmcomment{Note that the inimpp2 routine is general so that the original inimpp 355 326 routine should be suppressed from the code.} 356 327 357 328 When land processors are eliminated, the value corresponding to these locations in 358 the model output files is zero. Note that this is a problem for a mesh output file written 359 by such a model configuration, because model users often divide by the scale factors 360 ($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be 361 best not to eliminate land processors when running the model especially to write the 362 mesh files as outputs (when \np{nn\_msh} namelist parameter differs from 0). 363 %% 364 \gmcomment{Steven : dont understand this, no land processor means no output file 365 covering this part of globe; its only when files are stitched together into one that you 366 can leave a hole} 367 %% 329 the model output files is undefined. Note that this is a problem for the meshmask file 330 which requires to be defined over the whole domain. Therefore, user should not eliminate 331 land processors when creating a meshmask file ($i.e.$ when setting a non-zero value to \np{nn\_msh}). 368 332 369 333 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 380 344 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 381 345 382 383 % ================================================================384 % Open Boundary Conditions385 % ================================================================386 \section{Open Boundary Conditions (\key{obc}) (OBC)}387 \label{LBC_obc}388 %-----------------------------------------nam_obc -------------------------------------------389 %- nobc_dta = 0 ! = 0 the obc data are equal to the initial state390 %- ! = 1 the obc data are read in 'obc .dta' files391 %- rn_dpein = 1. ! ???392 %- rn_dpwin = 1. ! ???393 %- rn_dpnin = 30. ! ???394 %- rn_dpsin = 1. ! ???395 %- rn_dpeob = 1500. ! time relaxation (days) for the east open boundary396 %- rn_dpwob = 15. ! " " for the west open boundary397 %- rn_dpnob = 150. ! " " for the north open boundary398 %- rn_dpsob = 15. ! " " for the south open boundary399 %- ln_obc_clim = .true. ! climatological obc data files (default T)400 %- ln_vol_cst = .true. ! total volume conserved401 \namdisplay{namobc}402 403 It is often necessary to implement a model configuration limited to an oceanic404 region or a basin, which communicates with the rest of the global ocean through405 ''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a406 computational border where the aim of the calculations is to allow the perturbations407 generated inside the computational domain to leave it without deterioration of the408 inner model solution. However, an open boundary also has to let information from409 the outer ocean enter the model and should support inflow and outflow conditions.410 411 The open boundary package OBC is the first open boundary option developed in412 NEMO (originally in OPA8.2). It allows the user to413 \begin{itemize}414 \item tell the model that a boundary is ''open'' and not closed by a wall, for example415 by modifying the calculation of the divergence of velocity there;416 \item impose values of tracers and velocities at that boundary (values which may417 be taken from a climatology): this is the``fixed OBC'' option.418 \item calculate boundary values by a sophisticated algorithm combining radiation419 and relaxation (``radiative OBC'' option)420 \end{itemize}421 422 Options are defined through the \ngn{namobc} namelist variables.423 The package resides in the OBC directory. It is described here in four parts: the424 boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at425 the boundaries (module \mdl{obcdta}), the radiation algorithm involving the426 namelist and module \mdl{obcrad}, and a brief presentation of boundary update427 and restart files.428 429 %----------------------------------------------430 \subsection{Boundary geometry}431 \label{OBC_geom}432 %433 First one has to realize that open boundaries may not necessarily be located434 at the extremities of the computational domain. They may exist in the middle435 of the domain, for example at Gibraltar Straits if one wants to avoid including436 the Mediterranean in an Atlantic domain. This flexibility has been found necessary437 for the CLIPPER project \citep{Treguier_al_JGR01}. Because of the complexity of the438 geometry of ocean basins, it may even be necessary to have more than one439 ''west'' open boundary, more than one ''north'', etc. This is not possible with440 the OBC option: only one open boundary of each kind, west, east, south and441 north is allowed; these names refer to the grid geometry (not to the direction442 of the geographical ''west'', ''east'', etc).443 444 The open boundary geometry is set by a series of parameters in the module445 \mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east}446 (true if an east open boundary exists), \jp{jpieob} the $i$-index along which447 the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which448 it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ''d\'{e}but''449 and $f$ for ''fin'' in French). Similar parameters exist for the west, south and450 north cases (Table~\ref{Tab_obc_param}).451 452 453 %--------------------------------------------------TABLE--------------------------------------------------454 \begin{table}[htbp] \begin{center} \begin{tabular}{|l|c|c|c|}455 \hline456 Boundary and & Constant index & Starting index (d\'{e}but) & Ending index (fin) \\457 Logical flag & & & \\458 \hline459 West & \jp{jpiwob} $>= 2$ & \jp{jpjwd}$>= 2$ & \jp{jpjwf}<= \np{jpjglo}-1 \\460 lp\_obc\_west & $i$-index of a $u$ point & $j$ of a $T$ point &$j$ of a $T$ point \\461 \hline462 East & \jp{jpieob}$<=$\np{jpiglo}-2&\jp{jpjed} $>= 2$ & \jp{jpjef}$<=$ \np{jpjglo}-1 \\463 lp\_obc\_east & $i$-index of a $u$ point & $j$ of a $T$ point & $j$ of a $T$ point \\464 \hline465 South & \jp{jpjsob} $>= 2$ & \jp{jpisd} $>= 2$ & \jp{jpisf}$<=$\np{jpiglo}-1 \\466 lp\_obc\_south & $j$-index of a $v$ point & $i$ of a $T$ point & $i$ of a $T$ point \\467 \hline468 North & \jp{jpjnob} $<=$ \np{jpjglo}-2& \jp{jpind} $>= 2$ & \jp{jpinf}$<=$\np{jpiglo}-1 \\469 lp\_obc\_north & $j$-index of a $v$ point & $i$ of a $T$ point & $i$ of a $T$ point \\470 \hline471 \end{tabular} \end{center}472 \caption{ \label{Tab_obc_param}473 Names of different indices relating to the open boundaries. In the case474 of a completely open ocean domain with four ocean boundaries, the parameters475 take exactly the values indicated.}476 \end{table}477 %------------------------------------------------------------------------------------------------------------478 479 The open boundaries must be along coordinate lines. On the C-grid, the boundary480 itself is along a line of normal velocity points: $v$ points for a zonal open boundary481 (the south or north one), and $u$ points for a meridional open boundary (the west482 or east one). Another constraint is that there still must be a row of masked points483 all around the domain, as if the domain were a closed basin (unless periodic conditions484 are used together with open boundary conditions). Therefore, an open boundary485 cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also,486 the open boundary algorithm involves calculating the normal velocity points situated487 just on the boundary, as well as the tangential velocity and temperature and salinity488 just outside the boundary. This means that for a west/south boundary, normal489 velocities and temperature are calculated at the same index \jp{jpiwob} and490 \jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is491 calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is492 at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob}493 cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2.494 495 496 The starting and ending indices are to be thought of as $T$ point indices: in many497 cases they indicate the first land $T$-point, at the extremity of an open boundary498 (the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example499 of a northern open boundary). All indices are relative to the global domain. In the500 free surface case it is possible to have ``ocean corners'', that is, an open boundary501 starting and ending in the ocean.502 503 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>504 \begin{figure}[!t] \begin{center}505 \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf}506 \caption{ \label{Fig_obc_north}507 Localization of the North open boundary points.}508 \end{center} \end{figure}509 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>510 511 Although not compulsory, it is highly recommended that the bathymetry in the512 vicinity of an open boundary follows the following rule: in the direction perpendicular513 to the open line, the water depth should be constant for 4 grid points. This is in514 order to ensure that the radiation condition, which involves model variables next515 to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we516 indicate by an $=$ symbol, the points which should have the same depth. It means517 that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure518 why cylindrical}. The line behind the open $T$-line must be 0 in the bathymetry file519 (as shown on Fig.\ref{Fig_obc_north} for example).520 521 %----------------------------------------------522 \subsection{Boundary data}523 \label{OBC_data}524 525 It is necessary to provide information at the boundaries. The simplest case is526 when this information does not change in time and is equal to the initial conditions527 (namelist variable \np{nn\_obcdta}=0). This is the case for the standard configuration528 EEL5 with open boundaries. When (\np{nn\_obcdta}=1), open boundary information529 is read from netcdf files. For convenience the input files are supposed to be similar530 to the ''history'' NEMO output files, for dimension names and variable names.531 Open boundary arrays must be dimensioned according to the parameters of table~532 \ref{Tab_obc_param}: for example, at the western boundary, arrays have a533 dimension of \jp{jpwf}-\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical.534 535 When ocean observations are used to generate the boundary data (a hydrographic536 section for example, as in \citet{Treguier_al_JGR01}) it happens often that only the velocity537 normal to the boundary is known, which is the reason why the initial OBC code538 assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be539 specified. As more and more global model solutions and ocean analysis products540 become available, it will be possible to provide information about all the variables541 (including the tangential velocity) so that the specification of four variables at each542 boundaries will become standard. For the sea surface height, one must distinguish543 between the filtered free surface case and the time-splitting or explicit treatment of544 the free surface.545 In the first case, it is assumed that the user does not wish to represent high546 frequency motions such as tides. The boundary condition is thus one of zero547 normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}.548 No information other than the total velocity needs to be provided at the open549 boundaries in that case. In the other two cases (time splitting or explicit free surface),550 the user must provide barotropic information (sea surface height and barotropic551 velocities) and the use of the Flather algorithm for barotropic variables is552 recommanded. However, this algorithm has not yet been fully tested and bugs553 remain in NEMO v2.3. Users should read the code carefully before using it. Finally,554 in the case of the rigid lid approximation the barotropic streamfunction must be555 provided, as documented in \citet{Treguier_al_JGR01}). This option is no longer556 recommended but remains in NEMO V2.3.557 558 One frequently encountered case is when an open boundary domain is constructed559 from a global or larger scale NEMO configuration. Assuming the domain corresponds560 to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the561 small domain can be created by using the following netcdf utility on the global files:562 ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities,563 see their \href{http://nco.sourceforge.net}{website}).564 The open boundary files can be constructed using ncks565 commands, following table~\ref{Tab_obc_ind}.566 567 %--------------------------------------------------TABLE--------------------------------------------------568 \begin{table}[htbp] \begin{center} \begin{tabular}{|l|c|c|c|c|c|}569 \hline570 OBC & Variable & file name & Index & Start & end \\571 West & T,S & obcwest\_TS.nc & $ib$+1 & $jb$+1 & $je-1$ \\572 & U & obcwest\_U.nc & $ib$+1 & $jb$+1 & $je-1$ \\573 & V & obcwest\_V.nc & $ib$+1 & $jb$+1 & $je-1$ \\574 \hline575 East & T,S & obceast\_TS.nc & $ie$-1 & $jb$+1 & $je-1$ \\576 & U & obceast\_U.nc & $ie$-2 & $jb$+1 & $je-1$ \\577 & V & obceast\_V.nc & $ie$-1 & $jb$+1 & $je-1$ \\578 \hline579 South & T,S & obcsouth\_TS.nc & $jb$+1 & $ib$+1 & $ie-1$ \\580 & U & obcsouth\_U.nc & $jb$+1 & $ib$+1 & $ie-1$ \\581 & V & obcsouth\_V.nc & $jb$+1 & $ib$+1 & $ie-1$ \\582 \hline583 North & T,S & obcnorth\_TS.nc & $je$-1 & $ib$+1 & $ie-1$ \\584 & U & obcnorth\_U.nc & $je$-1 & $ib$+1 & $ie-1$ \\585 & V & obcnorth\_V.nc & $je$-2 & $ib$+1 & $ie-1$ \\586 \hline587 \end{tabular} \end{center}588 \caption{ \label{Tab_obc_ind}589 Requirements for creating open boundary files from a global configuration,590 appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the591 $i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global592 configuration, starting and ending with the $j$ or $i$ indices indicated.593 For example, to generate file obcnorth\_V.nc, use the command ncks594 $-F$ $-d\;y,je-2$ $-d\;x,ib+1,ie-1$ }595 \end{table}596 %-----------------------------------------------------------------------------------------------------------597 598 It is assumed that the open boundary files contain the variables for the period of599 the model integration. If the boundary files contain one time frame, the boundary600 data is held fixed in time. If the files contain 12 values, it is assumed that the input601 is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim}602 =true). The case of an arbitrary number of time frames is not yet implemented603 correctly; the user is required to write his own code in the module \mdl{obc\_dta}604 to deal with this situation.605 606 \subsection{Radiation algorithm}607 \label{OBC_rad}608 609 The art of open boundary management consists in applying a constraint strong610 enough that the inner domain "feels" the rest of the ocean, but weak enough611 that perturbations are allowed to leave the domain with minimum false reflections612 of energy. The constraints are specified separately at each boundary as time613 scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days)614 by namelist parameters such as \np{rn\_dpein}, \np{rn\_dpeob} for the eastern open615 boundary for example. When both time scales are zero for a given boundary616 ($e.g.$ for the western boundary, \jp{lp\_obc\_west}=true, \np{rn\_dpwob}=0 and617 \np{rn\_dpwin}=0) this means that the boundary in question is a ''fixed '' boundary618 where the solution is set exactly by the boundary data. This is not recommended,619 except in combination with increased viscosity in a ''sponge'' layer next to the620 boundary in order to avoid spurious reflections.621 622 623 The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output}624 algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is625 non-zero. It has been developed and tested in the SPEM model and its626 successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an627 $s$-coordinate model on an Arakawa C-grid. Although the algorithm has628 been numerically successful in the CLIPPER Atlantic models, the physics629 do not work as expected \citep{Treguier_al_JGR01}. Users are invited to consider630 open boundary conditions (OBC hereafter) with some scepticism631 \citep{Durran2001, Blayo2005}.632 633 The first part of the algorithm calculates a phase velocity to determine634 whether perturbations tend to propagate toward, or away from, the635 boundary. Let us consider a model variable $\phi$.636 The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$,637 in the directions normal and tangential to the boundary are638 \begin{equation} \label{Eq_obc_cphi}639 C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x}640 \;\;\;\;\; \;\;\;641 C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}.642 \end{equation}643 Following \citet{Treguier_al_JGR01} and \citet{Marchesiello2001} we retain only644 the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$645 (but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in646 the expression for $C_{\phi x}$).647 648 The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998},649 takes into account the two rows of grid points situated inside the domain650 next to the boundary, and the three previous time steps ($n$, $n-1$,651 and $n-2$). The same equation can then be discretized at the boundary at652 time steps $n-1$, $n$ and $n+1$ \gmcomment{since the original was three time-level}653 in order to extrapolate for the new boundary value $\phi^{n+1}$.654 655 In the open boundary algorithm as implemented in NEMO v2.3, the new boundary656 values are updated differently depending on the sign of $C_{\phi x}$. Let us take657 an eastern boundary as an example. The solution for variable $\phi$ at the658 boundary is given by a generalized wave equation with phase velocity $C_{\phi}$,659 with the addition of a relaxation term, as:660 \begin{eqnarray}661 \phi_{t} & = & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi)662 \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\663 \phi_{t} & = & \frac{1}{\tau_{i}} (\phi_{c}-\phi)664 \;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi}665 \end{eqnarray}666 where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary667 data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio668 $\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward669 propagation), the radiation condition (\ref{Eq_obc_rado}) is used.670 When $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is671 used with a strong relaxation to climatology (usually $\tau_{i}=\np{rn\_dpein}=$1~day).672 Equation (\ref{Eq_obc_radi}) is solved with a Euler time-stepping scheme. As a673 consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent674 to a fixed boundary condition. A time scale of one day is usually a good compromise675 which guarantees that the inflow conditions remain close to climatology while ensuring676 numerical stability.677 678 In the case of a western boundary located in the Eastern Atlantic, \citet{Penduff_al_JGR00}679 have been able to implement the radiation algorithm without any boundary data,680 using persistence from the previous time step instead. This solution has not worked681 in other cases \citep{Treguier_al_JGR01}, so that the use of boundary data is recommended.682 Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to683 maintain a weak relaxation to climatology. The time step is usually chosen so as to684 be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}).685 686 The radiation condition is applied to the model variables: temperature, salinity,687 tangential and normal velocities. For normal and tangential velocities, $u$ and $v$,688 radiation is applied with phase velocities calculated from $u$ and $v$ respectively.689 For the radiation of tracers, we use the phase velocity calculated from the tangential690 velocity in order to avoid calculating too many independent radiation velocities and691 because tangential velocities and tracers have the same position along the boundary692 on a C-grid.693 694 \subsection{Domain decomposition (\key{mpp\_mpi})}695 \label{OBC_mpp}696 When \key{mpp\_mpi} is active in the code, the computational domain is divided697 into rectangles that are attributed each to a different processor. The open boundary698 code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not699 work if there is an mpp subdomain boundary parallel to the open boundary at the700 index of the boundary, or the grid point after (outside), or three grid points before701 (inside). On the other hand, there is no problem if an mpp subdomain boundary702 cuts the open boundary perpendicularly. These geometrical limitations must be703 checked for by the user (there is no safeguard in the code).704 The general principle for the open boundary mpp code is that loops over the open705 boundaries {not sure what this means} are performed on local indices (nie0,706 nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module707 \mdl{obc\_ini}. Those indices have relevant values on the processors that contain708 a segment of an open boundary. For processors that do not include an open709 boundary segment, the indices are such that the calculations within the loops are710 not performed.711 \gmcomment{I dont understand most of the last few sentences}712 713 Arrays of climatological data that are read from files are seen by all processors714 and have the same dimensions for all (for instance, for the eastern boundary,715 uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation716 are local to each processor (uebnd(jpj,jpk,3,3) for instance). This allowed the717 CLIPPER model for example, to save on memory where the eastern boundary718 crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}-\jp{jpjed}+1).719 720 \subsection{Volume conservation}721 \label{OBC_vol}722 723 It is necessary to control the volume inside a domain when using open boundaries.724 With fixed boundaries, it is enough to ensure that the total inflow/outflow has725 reasonable values (either zero or a value compatible with an observed volume726 balance). When using radiative boundary conditions it is necessary to have a727 volume constraint because each open boundary works independently from the728 others. The methodology used to control this volume is identical to the one729 coded in the ROMS model \citep{Marchesiello2001}.730 731 732 %---------------------------------------- EXTRAS733 \colorbox{yellow}{Explain obc\_vol{\ldots}}734 735 \colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}}736 737 \colorbox{yellow}{OBC rigid lid? {\ldots}}738 346 739 347 % ==================================================================== -
trunk/DOC/TexFiles/Chapters/Chap_LDF.tex
r6140 r6289 69 69 r_{1u} &= \frac{e_{3u}}{ \left( e_{1u}\;\overline{\overline{e_{3w}}}^{\,i+1/2,\,k} \right)} 70 70 \;\delta_{i+1/2}[z_t] 71 &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] 71 &\approx \frac{1}{e_{1u}}\; \delta_{i+1/2}[z_t] \ \ \ 72 72 \\ 73 73 r_{2v} &= \frac{e_{3v}}{\left( e_{2v}\;\overline{\overline{e_{3w}}}^{\,j+1/2,\,k} \right)} 74 74 \;\delta_{j+1/2} [z_t] 75 &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] 75 &\approx \frac{1}{e_{2v}}\; \delta_{j+1/2}[z_t] \ \ \ 76 76 \\ 77 77 r_{1w} &= \frac{1}{e_{1w}}\;\overline{\overline{\delta_{i+1/2}[z_t]}}^{\,i,\,k+1/2} … … 448 448 \label{LDF_eiv} 449 449 450 %%gm from Triad appendix : to be incorporated.... 451 \gmcomment{ 452 Values of iso-neutral diffusivity and GM coefficient are set as 453 described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd}, 454 N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and 455 GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and 456 \np{rn\_aeiv\_0}. If 2D-varying coefficients are set with 457 \key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal 458 scale factor according to \eqref{Eq_title} \footnote{Except in global ORCA 459 $0.5^{\circ}$ runs with \key{traldf\_eiv}, where 460 $A_l$ is set like $A_e$ but with a minimum vale of 461 $100\;\mathrm{m}^2\;\mathrm{s}^{-1}$}. In idealised setups with 462 \key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv} 463 is set in the global configurations with \key{traldf\_c2d}, a horizontally varying $A_e$ is 464 instead set from the Held-Larichev parameterisation\footnote{In this 465 case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further 466 reduced by a factor $|f/f_{20}|$, where $f_{20}$ is the value of $f$ 467 at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored 468 unless it is zero. 469 } 470 450 471 When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), 451 472 an eddy induced tracer advection term is added, the formulation of which -
trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
r5118 r6289 1 1 % ================================================================ 2 % Chapter �Miscellaneous Topics2 % Chapter ——— Miscellaneous Topics 3 3 % ================================================================ 4 4 \chapter{Miscellaneous Topics} … … 34 34 has been made to set them in a generic way. However, examples of how 35 35 they can be set up is given in the ORCA 2\deg and 0.5\deg configurations. For example, 36 for details of implementation in ORCA2, search: 37 \vspace{-10pt} 38 \begin{alltt} 39 \tiny 40 \begin{verbatim} 41 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) 42 \end{verbatim} 43 \end{alltt} 36 for details of implementation in ORCA2, search: 37 \texttt{ IF( cp\_cfg == "orca" .AND. jp\_cfg == 2 ) } 44 38 45 39 % ------------------------------------------------------------------------------------------------------------- … … 80 74 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 81 75 82 % -------------------------------------------------------------------------------------------------------------83 % Cross Land Advection84 % -------------------------------------------------------------------------------------------------------------85 \subsection{Cross Land Advection (\mdl{tracla})}86 \label{MISC_strait_cla}87 %--------------------------------------------namcla--------------------------------------------------------88 \namdisplay{namcla}89 %--------------------------------------------------------------------------------------------------------------90 91 \colorbox{yellow}{Add a short description of CLA staff here or in lateral boundary condition chapter?}92 Options are defined through the \ngn{namcla} namelist variables.93 94 %The problem is resolved here by allowing the mixing of tracers and mass/volume between non-adjacent water columns at nominated regions within the model. Momentum is not mixed. The scheme conserves total tracer content, and total volume (the latter in $z*$- or $s*$-coordinate), and maintains compatibility between the tracer and mass/volume budgets.95 76 96 77 % ================================================================ … … 208 189 209 190 % ================================================================ 210 % Accelerating the Convergence211 % ================================================================212 \section{Accelerating the Convergence (\np{nn\_acc} = 1)}213 \label{MISC_acc}214 %--------------------------------------------namdom-------------------------------------------------------215 \namdisplay{namdom}216 %--------------------------------------------------------------------------------------------------------------217 218 Searching an equilibrium state with an global ocean model requires a very long time219 integration period (a few thousand years for a global model). Due to the size of220 the time step required for numerical stability (less than a few hours),221 this usually requires a large elapsed time. In order to overcome this problem,222 \citet{Bryan1984} introduces a technique that is intended to accelerate223 the spin up to equilibrium. It uses a larger time step in224 the tracer evolution equations than in the momentum evolution225 equations. It does not affect the equilibrium solution but modifies the226 trajectory to reach it.227 228 Options are defined through the \ngn{namdom} namelist variables.229 The acceleration of convergence option is used when \np{nn\_acc}=1. In that case,230 $\rdt=rn\_rdt$ is the time step of dynamics while $\widetilde{\rdt}=rdttra$ is the231 tracer time-step. the former is set from the \np{rn\_rdt} namelist parameter while the latter232 is computed using a hyperbolic tangent profile and the following namelist parameters :233 \np{rn\_rdtmin}, \np{rn\_rdtmax} and \np{rn\_rdth}. Those three parameters correspond234 to the surface value the deep ocean value and the depth at which the transition occurs, respectively.235 The set of prognostic equations to solve becomes:236 \begin{equation} \label{Eq_acc}237 \begin{split}238 \frac{\partial \textbf{U}_h }{\partial t}239 &\equiv \frac{\textbf{U}_h ^{t+1}-\textbf{U}_h^{t-1} }{2\rdt} = \ldots \\240 \frac{\partial T}{\partial t} &\equiv \frac{T^{t+1}-T^{t-1}}{2 \widetilde{\rdt}} = \ldots \\241 \frac{\partial S}{\partial t} &\equiv \frac{S^{t+1} -S^{t-1}}{2 \widetilde{\rdt}} = \ldots \\242 \end{split}243 \end{equation}244 245 \citet{Bryan1984} has examined the consequences of this distorted physics.246 Free waves have a slower phase speed, their meridional structure is slightly247 modified, and the growth rate of baroclinically unstable waves is reduced248 but with a wider range of instability. This technique is efficient for249 searching for an equilibrium state in coarse resolution models. However its250 application is not suitable for many oceanic problems: it cannot be used for251 transient or time evolving problems (in particular, it is very questionable252 to use this technique when there is a seasonal cycle in the forcing fields),253 and it cannot be used in high-resolution models where baroclinically254 unstable processes are important. Moreover, the vertical variation of255 $\widetilde{ \rdt}$ implies that the heat and salt contents are no longer256 conserved due to the vertical coupling of the ocean level through both257 advection and diffusion. Therefore \np{rn\_rdtmin} = \np{rn\_rdtmax} should be258 a more clever choice.259 260 261 % ================================================================262 191 % Accuracy and Reproducibility 263 192 % ================================================================ … … 370 299 the source of differences between mono and multi processor runs. 371 300 372 3- \key{esopa} (to be rename key\_nemo) : which is another option for model 373 management. When defined, this key forces the activation of all options and 374 CPP keys. For example, all tracer and momentum advection schemes are called! 375 Therefore the model results have no physical meaning. 376 However, this option forces both the compiler and the model to run through 377 all the \textsc{Fortran} lines of the model. This allows the user to check for obvious 378 compilation or execution errors with all CPP options, and errors in namelist options. 379 380 4- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of 301 %%gm to be removed both here and in the code 302 3- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of 381 303 a sum over the whole domain is performed as the summation over all processors of 382 304 each of their sums over their interior domains. This double sum never gives exactly … … 386 308 %THIS is to be updated with the mpp_sum_glo introduced in v3.3 387 309 % nn_bit_cmp today only check that the nn_cla = 0 (no cross land advection) 310 %%gm end 388 311 389 312 $\bullet$ Benchmark (\np{nn\_bench}). This option defines a benchmark run based on … … 393 316 or the physical parameterisations. 394 317 395 396 % ================================================================ 397 % Elliptic solvers (SOL) 398 % ================================================================ 399 \section{Elliptic solvers (SOL)} 400 \label{MISC_sol} 401 %--------------------------------------------namdom------------------------------------------------------- 402 \namdisplay{namsol} 403 %-------------------------------------------------------------------------------------------------------------- 404 405 When the filtered sea surface height option is used, the surface pressure gradient is 406 computed in \mdl{dynspg\_flt}. The force added in the momentum equation is solved implicitely. 407 It is thus solution of an elliptic equation \eqref{Eq_PE_flt} for which two solvers are available: 408 a Successive-Over-Relaxation scheme (SOR) and a preconditioned conjugate gradient 409 scheme(PCG) \citep{Madec_al_OM88, Madec_PhD90}. The solver is selected trough the 410 the value of \np{nn\_solv} \ngn{namsol} namelist variable. 411 412 The PCG is a very efficient method for solving elliptic equations on vector computers. 413 It is a fast and rather easy method to use; which are attractive features for a large 414 number of ocean situations (variable bottom topography, complex coastal geometry, 415 variable grid spacing, open or cyclic boundaries, etc ...). It does not require 416 a search for an optimal parameter as in the SOR method. However, the SOR has 417 been retained because it is a linear solver, which is a very useful property when 418 using the adjoint model of \NEMO. 419 420 At each time step, the time derivative of the sea surface height at time step $t+1$ 421 (or equivalently the divergence of the \textit{after} barotropic transport) that appears 422 in the filtering forced is the solution of the elliptic equation obtained from the horizontal 423 divergence of the vertical summation of \eqref{Eq_PE_flt}. 424 Introducing the following coefficients: 425 \begin{equation} \label{Eq_sol_matrix} 426 \begin{aligned} 427 &c_{i,j}^{NS} &&= {2 \rdt }^2 \; \frac{H_v (i,j) \; e_{1v} (i,j)}{e_{2v}(i,j)} \\ 428 &c_{i,j}^{EW} &&= {2 \rdt }^2 \; \frac{H_u (i,j) \; e_{2u} (i,j)}{e_{1u}(i,j)} \\ 429 &b_{i,j} &&= \delta_i \left[ e_{2u}M_u \right] - \delta_j \left[ e_{1v}M_v \right]\ , \\ 430 \end{aligned} 431 \end{equation} 432 the resulting five-point finite difference equation is given by: 433 \begin{equation} \label{Eq_solmat} 434 \begin{split} 435 c_{i+1,j}^{NS} D_{i+1,j} + \; c_{i,j+1}^{EW} D_{i,j+1} 436 + c_{i,j} ^{NS} D_{i-1,j} + \; c_{i,j} ^{EW} D_{i,j-1} & \\ 437 - \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} 438 \end{split} 439 \end{equation} 440 \eqref{Eq_solmat} is a linear symmetric system of equations. All the elements of 441 the corresponding matrix \textbf{A} vanish except those of five diagonals. With 442 the natural ordering of the grid points (i.e. from west to east and from 443 south to north), the structure of \textbf{A} is block-tridiagonal with 444 tridiagonal or diagonal blocks. \textbf{A} is a positive-definite symmetric 445 matrix of size $(jpi \cdot jpj)^2$, and \textbf{B}, the right hand side of 446 \eqref{Eq_solmat}, is a vector. 447 448 Note that in the linear free surface case, the depth that appears in \eqref{Eq_sol_matrix} 449 does not vary with time, and thus the matrix can be computed once for all. In non-linear free surface 450 (\key{vvl} defined) the matrix have to be updated at each time step. 451 452 % ------------------------------------------------------------------------------------------------------------- 453 % Successive Over Relaxation 454 % ------------------------------------------------------------------------------------------------------------- 455 \subsection{Successive Over Relaxation (\np{nn\_solv}=2, \mdl{solsor})} 456 \label{MISC_solsor} 457 458 Let us introduce the four cardinal coefficients: 459 \begin{align*} 460 a_{i,j}^S &= c_{i,j }^{NS}/d_{i,j} &\qquad a_{i,j}^W &= c_{i,j}^{EW}/d_{i,j} \\ 461 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} 462 \end{align*} 463 where $d_{i,j} = c_{i,j}^{NS}+ c_{i+1,j}^{NS} + c_{i,j}^{EW} + c_{i,j+1}^{EW}$ 464 (i.e. the diagonal of the matrix). \eqref{Eq_solmat} can be rewritten as: 465 \begin{equation} \label{Eq_solmat_p} 466 \begin{split} 467 a_{i,j}^{N} D_{i+1,j} +\,a_{i,j}^{E} D_{i,j+1} +\, a_{i,j}^{S} D_{i-1,j} +\,a_{i,j}^{W} D_{i,j-1} - D_{i,j} = \tilde{b}_{i,j} 468 \end{split} 469 \end{equation} 470 with $\tilde b_{i,j} = b_{i,j}/d_{i,j}$. \eqref{Eq_solmat_p} is the equation actually solved 471 with the SOR method. This method used is an iterative one. Its algorithm can be 472 summarised as follows (see \citet{Haltiner1980} for a further discussion): 473 474 initialisation (evaluate a first guess from previous time step computations) 475 \begin{equation} 476 D_{i,j}^0 = 2 \, D_{i,j}^t - D_{i,j}^{t-1} 477 \end{equation} 478 iteration $n$, from $n=0$ until convergence, do : 479 \begin{equation} \label{Eq_sor_algo} 480 \begin{split} 481 R_{i,j}^n = &a_{i,j}^{N} D_{i+1,j}^n +\,a_{i,j}^{E} D_{i,j+1} ^n 482 +\, a_{i,j}^{S} D_{i-1,j} ^{n+1}+\,a_{i,j}^{W} D_{i,j-1} ^{n+1} 483 - D_{i,j}^n - \tilde{b}_{i,j} \\ 484 D_{i,j} ^{n+1} = &D_{i,j} ^{n} + \omega \;R_{i,j}^n 485 \end{split} 486 \end{equation} 487 where \textit{$\omega $ }satisfies $1\leq \omega \leq 2$. An optimal value exists for 488 \textit{$\omega$} which significantly accelerates the convergence, but it has to be 489 adjusted empirically for each model domain (except for a uniform grid where an 490 analytical expression for \textit{$\omega$} can be found \citep{Richtmyer1967}). 491 The value of $\omega$ is set using \np{rn\_sor}, a \textbf{namelist} parameter. 492 The convergence test is of the form: 493 \begin{equation} 494 \delta = \frac{\sum\limits_{i,j}{R_{i,j}^n}{R_{i,j}^n}} 495 {\sum\limits_{i,j}{ \tilde{b}_{i,j}^n}{\tilde{b}_{i,j}^n}} \leq \epsilon 496 \end{equation} 497 where $\epsilon$ is the absolute precision that is required. It is recommended 498 that a value smaller or equal to $10^{-6}$ is used for $\epsilon$ since larger 499 values may lead to numerically induced basin scale barotropic oscillations. 500 The precision is specified by setting \np{rn\_eps} (\textbf{namelist} parameter). 501 In addition, two other tests are used to halt the iterative algorithm. They involve 502 the number of iterations and the modulus of the right hand side. If the former 503 exceeds a specified value, \np{nn\_max} (\textbf{namelist} parameter), 504 or the latter is greater than $10^{15}$, the whole model computation is stopped 505 and the last computed time step fields are saved in a abort.nc NetCDF file. 506 In both cases, this usually indicates that there is something wrong in the model 507 configuration (an error in the mesh, the initial state, the input forcing, 508 or the magnitude of the time step or of the mixing coefficients). A typical value of 509 $nn\_max$ is a few hundred when $\epsilon = 10^{-6}$, increasing to a few 510 thousand when $\epsilon = 10^{-12}$. 511 The vectorization of the SOR algorithm is not straightforward. The scheme 512 contains two linear recurrences on $i$ and $j$. This inhibits the vectorisation. 513 \eqref{Eq_sor_algo} can be been rewritten as: 514 \begin{equation} 515 \begin{split} 516 R_{i,j}^n 517 = &a_{i,j}^{N} D_{i+1,j}^n +\,a_{i,j}^{E} D_{i,j+1} ^n 518 +\,a_{i,j}^{S} D_{i-1,j} ^{n}+\,_{i,j}^{W} D_{i,j-1} ^{n} - D_{i,j}^n - \tilde{b}_{i,j} \\ 519 R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{S}\; R_{i,j-1}^n \\ 520 R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{W}\; R_{i-1,j}^n 521 \end{split} 522 \end{equation} 523 This technique slightly increases the number of iteration required to reach the convergence, 524 but this is largely compensated by the gain obtained by the suppression of the recurrences. 525 526 Another technique have been chosen, the so-called red-black SOR. It consist in solving successively 527 \eqref{Eq_sor_algo} for odd and even grid points. It also slightly reduced the convergence rate 528 but allows the vectorisation. In addition, and this is the reason why it has been chosen, it is able to handle the north fold boundary condition used in ORCA configuration ($i.e.$ tri-polar global ocean mesh). 529 530 The SOR method is very flexible and can be used under a wide range of conditions, 531 including irregular boundaries, interior boundary points, etc. Proofs of convergence, etc. 532 may be found in the standard numerical methods texts for partial differential equations. 533 534 % ------------------------------------------------------------------------------------------------------------- 535 % Preconditioned Conjugate Gradient 536 % ------------------------------------------------------------------------------------------------------------- 537 \subsection{Preconditioned Conjugate Gradient (\np{nn\_solv}=1, \mdl{solpcg}) } 538 \label{MISC_solpcg} 539 540 \textbf{A} is a definite positive symmetric matrix, thus solving the linear 541 system \eqref{Eq_solmat} is equivalent to the minimisation of a quadratic 542 functional: 543 \begin{equation*} 544 \textbf{Ax} = \textbf{b} \leftrightarrow \textbf{x} =\text{inf}_{y} \,\phi (\textbf{y}) 545 \quad , \qquad 546 \phi (\textbf{y}) = 1/2 \langle \textbf{Ay},\textbf{y}\rangle - \langle \textbf{b},\textbf{y} \rangle 547 \end{equation*} 548 where $\langle , \rangle$ is the canonical dot product. The idea of the 549 conjugate gradient method is to search for the solution in the following 550 iterative way: assuming that $\textbf{x}^n$ has been obtained, $\textbf{x}^{n+1}$ 551 is found from $\textbf {x}^{n+1}={\textbf {x}}^n+\alpha^n{\textbf {d}}^n$ which satisfies: 552 \begin{equation*} 553 {\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 554 \end{equation*} 555 and expressing $\phi (\textbf{y})$ as a function of \textit{$\alpha $}, we obtain the 556 value that minimises the functional: 557 \begin{equation*} 558 \alpha ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{ A d}^n, \textbf{d}^n} \rangle 559 \end{equation*} 560 where $\textbf{r}^n = \textbf{b}-\textbf{A x}^n = \textbf{A} (\textbf{x}-\textbf{x}^n)$ 561 is the error at rank $n$. The descent vector $\textbf{d}^n$ s chosen to be dependent 562 on the error: $\textbf{d}^n = \textbf{r}^n + \beta^n \,\textbf{d}^{n-1}$. $\beta ^n$ 563 is searched such that the descent vectors form an orthogonal basis for the dot 564 product linked to \textbf{A}. Expressing the condition 565 $\langle \textbf{A d}^n, \textbf{d}^{n-1} \rangle = 0$ the value of $\beta ^n$ is found: 566 $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{r}^{n-1}, \textbf{r}^{n-1}} \rangle$. 567 As a result, the errors $ \textbf{r}^n$ form an orthogonal 568 base for the canonic dot product while the descent vectors $\textbf{d}^n$ form 569 an orthogonal base for the dot product linked to \textbf{A}. The resulting 570 algorithm is thus the following one: 571 572 initialisation : 573 \begin{equation*} 574 \begin{split} 575 \textbf{x}^0 &= D_{i,j}^0 = 2 D_{i,j}^t - D_{i,j}^{t-1} \quad, \text{the initial guess } \\ 576 \textbf{r}^0 &= \textbf{d}^0 = \textbf{b} - \textbf{A x}^0 \\ 577 \gamma_0 &= \langle{ \textbf{r}^0 , \textbf{r}^0} \rangle 578 \end{split} 579 \end{equation*} 580 581 iteration $n,$ from $n=0$ until convergence, do : 582 \begin{equation} 583 \begin{split} 584 \text{z}^n& = \textbf{A d}^n \\ 585 \alpha_n &= \gamma_n / \langle{ \textbf{z}^n , \textbf{d}^n} \rangle \\ 586 \textbf{x}^{n+1} &= \textbf{x}^n + \alpha_n \,\textbf{d}^n \\ 587 \textbf{r}^{n+1} &= \textbf{r}^n - \alpha_n \,\textbf{z}^n \\ 588 \gamma_{n+1} &= \langle{ \textbf{r}^{n+1} , \textbf{r}^{n+1}} \rangle \\ 589 \beta_{n+1} &= \gamma_{n+1}/\gamma_{n} \\ 590 \textbf{d}^{n+1} &= \textbf{r}^{n+1} + \beta_{n+1}\; \textbf{d}^{n}\\ 591 \end{split} 592 \end{equation} 593 594 595 The convergence test is: 596 \begin{equation} 597 \delta = \gamma_{n}\; / \langle{ \textbf{b} , \textbf{b}} \rangle \leq \epsilon 598 \end{equation} 599 where $\epsilon $ is the absolute precision that is required. As for the SOR algorithm, 600 the whole model computation is stopped when the number of iterations, \np{nn\_max}, or 601 the modulus of the right hand side of the convergence equation exceeds a 602 specified value (see \S\ref{MISC_solsor} for a further discussion). The required 603 precision and the maximum number of iterations allowed are specified by setting 604 \np{rn\_eps} and \np{nn\_max} (\textbf{namelist} parameters). 605 606 It can be demonstrated that the above algorithm is optimal, provides the exact 607 solution in a number of iterations equal to the size of the matrix, and that 608 the convergence rate is faster as the matrix is closer to the identity matrix, 609 $i.e.$ its eigenvalues are closer to 1. Therefore, it is more efficient to solve 610 a better conditioned system which has the same solution. For that purpose, 611 we introduce a preconditioning matrix \textbf{Q} which is an approximation 612 of \textbf{A} but much easier to invert than \textbf{A}, and solve the system: 613 \begin{equation} \label{Eq_pmat} 614 \textbf{Q}^{-1} \textbf{A x} = \textbf{Q}^{-1} \textbf{b} 615 \end{equation} 616 617 The same algorithm can be used to solve \eqref{Eq_pmat} if instead of the 618 canonical dot product the following one is used: 619 ${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, and 620 if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ 621 are substituted to \textbf{b} and \textbf{A} \citep{Madec_al_OM88}. 622 In \NEMO, \textbf{Q} is chosen as the diagonal of \textbf{ A}, i.e. the simplest form for 623 \textbf{Q} so that it can be easily inverted. In this case, the discrete formulation of 624 \eqref{Eq_pmat} is in fact given by \eqref{Eq_solmat_p} and thus the matrix and 625 right hand side are computed independently from the solver used. 626 627 % ================================================================ 628 629 630 631 632 633 634 635 636 637 638 639 318 % ================================================================ 319 320 321 322 323 -
trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex
r6140 r6289 257 257 258 258 %\newpage 259 %$\ $\newline % force a new li gne259 %$\ $\newline % force a new line 260 260 261 261 % ================================================================ … … 988 988 \label{PE_zco_tilde} 989 989 990 The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM1 0s}.990 The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM11}. 991 991 It is available in \NEMO since the version 3.4. Nevertheless, it is currently not robust enough 992 992 to be used in all possible configurations. Its use is therefore not recommended. 993 We 993 994 994 995 995 \newpage … … 1113 1113 1114 1114 All these parameterisations of subgrid scale physics have advantages and 1115 drawbacks. There are not all available in \NEMO. In the $z$-coordinate 1116 formulation, five options are offered for active tracers (temperature and 1117 salinity): second order geopotential operator, second order isoneutral 1118 operator, \citet{Gent1990} parameterisation, fourth order 1119 geopotential operator, and various slightly diffusive advection schemes. 1120 The same options are available for momentum, except 1121 \citet{Gent1990} parameterisation which only involves tracers. In the 1122 $s$-coordinate formulation, additional options are offered for tracers: second 1123 order operator acting along $s-$surfaces, and for momentum: fourth order 1124 operator acting along $s-$surfaces (see \S\ref{LDF}). 1125 1126 \subsubsection{Lateral second order tracer diffusive operator} 1127 1128 The lateral second order tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}): 1115 drawbacks. There are not all available in \NEMO. For active tracers (temperature and 1116 salinity) the main ones are: Laplacian and bilaplacian operators acting along 1117 geopotential or iso-neutral surfaces, \citet{Gent1990} parameterisation, 1118 and various slightly diffusive advection schemes. 1119 For momentum, the main ones are: Laplacian and bilaplacian operators acting along 1120 geopotential surfaces, and UBS advection schemes when flux form is chosen for the momentum advection. 1121 1122 \subsubsection{Lateral Laplacian tracer diffusive operator} 1123 1124 The lateral Laplacian tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}): 1129 1125 \begin{equation} \label{Eq_PE_iso_tensor} 1130 1126 D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad … … 1158 1154 but have similar expressions in $z$- and $s$-coordinates. In $z$-coordinates: 1159 1155 \begin{equation} \label{Eq_PE_iso_slopes} 1160 r_1 =\frac{e_3 }{e_1 } \left( {\frac{\partial \rho }{\partial i}} \right) 1161 \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \ , \quad 1162 r_1 =\frac{e_3 }{e_1 } \left( {\frac{\partial \rho }{\partial i}} \right) 1163 \left( {\frac{\partial \rho }{\partial k}} \right)^{-1}, 1164 \end{equation} 1165 while in $s$-coordinates $\partial/\partial k$ is replaced by 1166 $\partial/\partial s$. 1156 r_1 =\frac{e_3 }{e_1 } \left( \pd[\rho]{i} \right) \left( \pd[\rho]{k} \right)^{-1} \, \quad 1157 r_2 =\frac{e_3 }{e_2 } \left( \pd[\rho]{j} \right) \left( \pd[\rho]{k} \right)^{-1} \, 1158 \end{equation} 1159 while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. 1167 1160 1168 1161 \subsubsection{Eddy induced velocity} … … 1181 1174 w^\ast &= -\frac{1}{e_1 e_2 }\left[ 1182 1175 \frac{\partial }{\partial i}\left( {A^{eiv}\;e_2\,\tilde{r}_1 } \right) 1183 +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right) \right]1176 +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right) \right] 1184 1177 \end{split} 1185 1178 \end{equation} … … 1190 1183 \begin{align} \label{Eq_PE_slopes_eiv} 1191 1184 \tilde{r}_n = \begin{cases} 1192 r_n 1185 r_n & \text{in $z$-coordinate} \\ 1193 1186 r_n + \sigma_n & \text{in \textit{z*} and $s$-coordinates} 1194 1187 \end{cases} 1195 1188 \quad \text{where } n=1,2 1196 1189 \end{align} … … 1200 1193 to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}). 1201 1194 1202 \subsubsection{Lateral fourth ordertracer diffusive operator}1203 1204 The lateral fourth ordertracer diffusive operator is defined by:1195 \subsubsection{Lateral bilaplacian tracer diffusive operator} 1196 1197 The lateral bilaplacian tracer diffusive operator is defined by: 1205 1198 \begin{equation} \label{Eq_PE_bilapT} 1206 D^{lT}= \Delta \left( \;\Delta T \right)1199 D^{lT}= - \Delta \left( \;\Delta T \right) 1207 1200 \qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right) 1208 1201 \end{equation} 1209 It is the second orderoperator given by \eqref{Eq_PE_iso_tensor} applied twice with1202 It is the Laplacian operator given by \eqref{Eq_PE_iso_tensor} applied twice with 1210 1203 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1211 1204 1212 1205 1213 \subsubsection{Lateral second ordermomentum diffusive operator}1214 1215 The second ordermomentum diffusive operator along $z$- or $s$-surfaces is found by1206 \subsubsection{Lateral Laplacian momentum diffusive operator} 1207 1208 The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by 1216 1209 applying \eqref{Eq_PE_lap_vector} to the horizontal velocity vector (see Appendix~\ref{Apdx_B}): 1217 1210 \begin{equation} \label{Eq_PE_lapU} … … 1247 1240 of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}. 1248 1241 1249 \subsubsection{lateral fourth order momentum diffusive operator} 1250 1251 As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces 1252 is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU_iso} 1253 with the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1254 1242 \subsubsection{lateral bilaplacian momentum diffusive operator} 1243 1244 As for tracers, the bilaplacian order momentum diffusive operator is a 1245 re-entering Laplacian operator with the harmonic eddy diffusion coefficient 1246 set to the square root of the biharmonic one. Nevertheless it is currently 1247 not available in the iso-neutral case. 1248 -
trunk/DOC/TexFiles/Chapters/Chap_SBC.tex
r6140 r6289 1 1 % ================================================================ 2 % Chapter �Surface Boundary Condition (SBC, ISF, ICB)2 % Chapter —— Surface Boundary Condition (SBC, ISF, ICB) 3 3 % ================================================================ 4 4 \chapter{Surface Boundary Condition (SBC, ISF, ICB) } … … 17 17 \item the two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$ 18 18 \item the incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$ 19 \item the surface freshwater budget $\left( {\textit{emp},\;\textit{emp}_S } \right)$ 19 \item the surface freshwater budget $\left( {\textit{emp}} \right)$ 20 \item the surface salt flux associated with freezing/melting of seawater $\left( {\textit{sfx}} \right)$ 20 21 \end{itemize} 21 22 plus an optional field: … … 27 28 are controlled by namelist \ngn{namsbc} variables: an analytical formulation (\np{ln\_ana}~=~true), 28 29 a flux formulation (\np{ln\_flx}~=~true), a bulk formulae formulation (CORE 29 (\np{ln\_ core}~=~true), CLIO (\np{ln\_clio}~=~true) or MFS30 (\np{ln\_blk\_core}~=~true), CLIO (\np{ln\_blk\_clio}~=~true) or MFS 30 31 \footnote { Note that MFS bulk formulae compute fluxes only for the ocean component} 31 (\np{ln\_mfs}~=~true) bulk formulae) and a coupled 32 formulation (exchanges with a atmospheric model via the OASIS coupler) 33 (\np{ln\_cpl}~=~true). When used, the atmospheric pressure forces both 34 ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true). 35 The frequency at which the six or seven fields have to be updated is the \np{nn\_fsbc} 36 namelist parameter. 32 (\np{ln\_blk\_mfs}~=~true) bulk formulae) and a coupled or mixed forced/coupled formulation 33 (exchanges with a atmospheric model via the OASIS coupler) (\np{ln\_cpl} or \np{ln\_mixcpl}~=~true). 34 When used ($i.e.$ \np{ln\_apr\_dyn}~=~true), the atmospheric pressure forces both ocean and ice dynamics. 35 36 The frequency at which the forcing fields have to be updated is given by the \np{nn\_fsbc} namelist parameter. 37 37 When the fields are supplied from data files (flux and bulk formulations), the input fields 38 need not be supplied on the model grid. 38 need not be supplied on the model grid. Instead a file of coordinates and weights can 39 39 be supplied which maps the data from the supplied grid to the model points 40 40 (so called "Interpolation on the Fly", see \S\ref{SBC_iof}). … … 42 42 can be masked to avoid spurious results in proximity of the coasts as large sea-land gradients characterize 43 43 most of the atmospheric variables. 44 44 45 In addition, the resulting fields can be further modified using several namelist options. 45 These options control the rotation of vector components supplied relative to an east-north 46 coordinate system onto the local grid directions in the model; the addition of a surface 47 restoring term to observed SST and/or SSS (\np{ln\_ssr}~=~true); the modification of fluxes 48 below ice-covered areas (using observed ice-cover or a sea-ice model) 49 (\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater 50 fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of isf melting as lateral inflow (parameterisation) 51 or as surface flux at the land-ice ocean interface (\np{ln\_isf}=~true); 52 the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the 53 transformation of the solar radiation (if provided as daily mean) into a diurnal 54 cycle (\np{ln\_dm2dc}~=~true); and a neutral drag coefficient can be read from an external wave 55 model (\np{ln\_cdgw}~=~true). The latter option is possible only in case core or mfs bulk formulas are selected. 46 These options control 47 \begin{itemize} 48 \item the rotation of vector components supplied relative to an east-north 49 coordinate system onto the local grid directions in the model ; 50 \item the addition of a surface restoring term to observed SST and/or SSS (\np{ln\_ssr}~=~true) ; 51 \item the modification of fluxes below ice-covered areas (using observed ice-cover or a sea-ice model) (\np{nn\_ice}~=~0,1, 2 or 3) ; 52 \item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln\_rnf}~=~true) ; 53 \item the addition of isf melting as lateral inflow (parameterisation) (\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false) 54 or as fluxes applied at the land-ice ocean interface (\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true) ; 55 \item the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2) ; 56 \item the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle (\np{ln\_dm2dc}~=~true) ; 57 and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}~=~true). 58 \end{itemize} 59 The latter option is possible only in case core or mfs bulk formulas are selected. 56 60 57 61 In this chapter, we first discuss where the surface boundary condition appears in the … … 72 76 73 77 The surface ocean stress is the stress exerted by the wind and the sea-ice 74 on the ocean. The two components of stress are assumed to be interpolated 75 onto the ocean mesh, $i.e.$ resolved onto the model (\textbf{i},\textbf{j}) direction 76 at $u$- and $v$-points They are applied as a surface boundary condition of the 77 computation of the momentum vertical mixing trend (\mdl{dynzdf} module) : 78 \begin{equation} \label{Eq_sbc_dynzdf} 79 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 80 = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v } 81 \end{equation} 82 where $(\tau _u ,\;\tau _v )=(utau,vtau)$ are the two components of the wind 83 stress vector in the $(\textbf{i},\textbf{j})$ coordinate system. 78 on the ocean. It is applied in \mdl{dynzdf} module as a surface boundary condition of the 79 computation of the momentum vertical mixing trend (see \eqref{Eq_dynzdf_sbc} in \S\ref{DYN_zdf}). 80 As such, it has to be provided as a 2D vector interpolated 81 onto the horizontal velocity ocean mesh, $i.e.$ resolved onto the model 82 (\textbf{i},\textbf{j}) direction at $u$- and $v$-points. 84 83 85 84 The surface heat flux is decomposed into two parts, a non solar and a solar heat 86 85 flux, $Q_{ns}$ and $Q_{sr}$, respectively. The former is the non penetrative part 87 of the heat flux ($i.e.$ the sum of sensible, latent and long wave heat fluxes). 88 It is applied as a surface boundary condition trend of the first level temperature 89 time evolution equation (\mdl{trasbc} module). 90 \begin{equation} \label{Eq_sbc_trasbc_q} 91 \frac{\partial T}{\partial t}\equiv \cdots \;+\;\left. {\frac{Q_{ns} }{\rho 92 _o \;C_p \;e_{3t} }} \right|_{k=1} \quad 93 \end{equation} 94 $Q_{sr}$ is the penetrative part of the heat flux. It is applied as a 3D 95 trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=True. 96 97 \begin{equation} \label{Eq_sbc_traqsr} 98 \frac{\partial T}{\partial t}\equiv \cdots \;+\frac{Q_{sr} }{\rho_o C_p \,e_{3t} }\delta _k \left[ {I_w } \right] 99 \end{equation} 100 where $I_w$ is a non-dimensional function that describes the way the light 101 penetrates inside the water column. It is generally a sum of decreasing 102 exponentials (see \S\ref{TRA_qsr}). 103 104 The surface freshwater budget is provided by fields: \textit{emp} and $\textit{emp}_S$ which 105 may or may not be identical. Indeed, a surface freshwater flux has two effects: 106 it changes the volume of the ocean and it changes the surface concentration of 107 salt (and other tracers). Therefore it appears in the sea surface height as a volume 108 flux, \textit{emp} (\textit{dynspg\_xxx} modules), and in the salinity time evolution equations 109 as a concentration/dilution effect, 110 $\textit{emp}_{S}$ (\mdl{trasbc} module). 111 \begin{equation} \label{Eq_trasbc_emp} 112 \begin{aligned} 113 &\frac{\partial \eta }{\partial t}\equiv \cdots \;+\;\textit{emp}\quad \\ 114 \\ 115 &\frac{\partial S}{\partial t}\equiv \cdots \;+\left. {\frac{\textit{emp}_S \;S}{e_{3t} }} \right|_{k=1} \\ 116 \end{aligned} 117 \end{equation} 118 119 In the real ocean, $\textit{emp}=\textit{emp}_S$ and the ocean salt content is conserved, 120 but it exist several numerical reasons why this equality should be broken. 121 For example, when the ocean is coupled to a sea-ice model, the water exchanged between 122 ice and ocean is slightly salty (mean sea-ice salinity is $\sim $\textit{4 psu}). In this case, 123 $\textit{emp}_{S}$ take into account both concentration/dilution effect associated with 124 freezing/melting and the salt flux between ice and ocean, while \textit{emp} is 125 only the volume flux. In addition, in the current version of \NEMO, the sea-ice is 126 assumed to be above the ocean (the so-called levitating sea-ice). Freezing/melting does 127 not change the ocean volume (no impact on \textit{emp}) but it modifies the SSS. 128 %gm \colorbox{yellow}{(see {\S} on LIM sea-ice model)}. 129 130 Note that SST can also be modified by a freshwater flux. Precipitation (in 131 particular solid precipitation) may have a temperature significantly different from 132 the SST. Due to the lack of information about the temperature of 133 precipitation, we assume it is equal to the SST. Therefore, no 134 concentration/dilution term appears in the temperature equation. It has to 135 be emphasised that this absence does not mean that there is no heat flux 136 associated with precipitation! Precipitation can change the ocean volume and thus the 137 ocean heat content. It is therefore associated with a heat flux (not yet 138 diagnosed in the model) \citep{Roullet_Madec_JGR00}). 86 of the heat flux ($i.e.$ the sum of sensible, latent and long wave heat fluxes 87 plus the heat content of the mass exchange with the atmosphere and sea-ice). 88 It is applied in \mdl{trasbc} module as a surface boundary condition trend of 89 the first level temperature time evolution equation (see \eqref{Eq_tra_sbc} 90 and \eqref{Eq_tra_sbc_lin} in \S\ref{TRA_sbc}). 91 The latter is the penetrative part of the heat flux. It is applied as a 3D 92 trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=\textit{true}. 93 The way the light penetrates inside the water column is generally a sum of decreasing 94 exponentials (see \S\ref{TRA_qsr}). 95 96 The surface freshwater budget is provided by the \textit{emp} field. 97 It represents the mass flux exchanged with the atmosphere (evaporation minus precipitation) 98 and possibly with the sea-ice and ice shelves (freezing minus melting of ice). 99 It affects both the ocean in two different ways: 100 $(i)$ it changes the volume of the ocean and therefore appears in the sea surface height 101 equation as a volume flux, and 102 $(ii)$ it changes the surface temperature and salinity through the heat and salt contents 103 of the mass exchanged with the atmosphere, the sea-ice and the ice shelves. 104 139 105 140 106 %\colorbox{yellow}{Miss: } … … 156 122 % 157 123 %Explain here all the namlist namsbc variable{\ldots}. 124 % 125 % explain : use or not of surface currents 158 126 % 159 127 %\colorbox{yellow}{End Miss } 160 128 161 The ocean model provides the surface currents, temperature and salinity 162 averaged over \np{nf\_sbc} time-step (\ref{Tab_ssm}).The computation of the 163 mean is done in \mdl{sbcmod} module. 129 The ocean model provides, at each time step, to the surface module (\mdl{sbcmod}) 130 the surface currents, temperature and salinity. 131 These variables are averaged over \np{nf\_sbc} time-step (\ref{Tab_ssm}), 132 and it is these averaged fields which are used to computes the surface fluxes 133 at a frequency of \np{nf\_sbc} time-step. 134 164 135 165 136 %-------------------------------------------------TABLE--------------------------------------------------- … … 458 429 %-------------------------------------------------------------------------------------------------------------- 459 430 460 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.461 Options are defined through the \ngn{namsbc\_sas} namelist variables. 431 In some circumstances it may be useful to avoid calculating the 3D temperature, salinity and velocity fields 432 and simply read them in from a previous run or receive them from OASIS. 462 433 For example: 463 434 464 \begin{ enumerate}465 \item Multiple runs of the model are required in code development to see the affect of different algorithms in435 \begin{itemize} 436 \item Multiple runs of the model are required in code development to see the effect of different algorithms in 466 437 the bulk formulae. 467 438 \item The effect of different parameter sets in the ice model is to be examined. 468 \end{enumerate} 439 \item Development of sea-ice algorithms or parameterizations. 440 \item spinup of the iceberg floats 441 \item ocean/sea-ice simulation with both media running in parallel (\np{ln\_mixcpl}~=~\textit{true}) 442 \end{itemize} 469 443 470 444 The StandAlone Surface scheme provides this utility. 445 Its options are defined through the \ngn{namsbc\_sas} namelist variables. 471 446 A new copy of the model has to be compiled with a configuration based on ORCA2\_SAS\_LIM. 472 447 However no namelist parameters need be changed from the settings of the previous run (except perhaps nn{\_}date0) … … 474 449 Routines replaced are: 475 450 476 \begin{enumerate} 477 \item \mdl{nemogcm} 478 479 This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (step.F90) 451 \begin{itemize} 452 \item \mdl{nemogcm} : This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (step.F90) 480 453 Since the ocean state is not calculated all associated initialisations have been removed. 481 \item \mdl{step} 482 483 The main time stepping routine now only needs to call the sbc routine (and a few utility functions). 484 \item \mdl{sbcmod} 485 486 This has been cut down and now only calculates surface forcing and the ice model required. New surface modules 454 \item \mdl{step} : The main time stepping routine now only needs to call the sbc routine (and a few utility functions). 455 \item \mdl{sbcmod} : This has been cut down and now only calculates surface forcing and the ice model required. New surface modules 487 456 that can function when only the surface level of the ocean state is defined can also be added (e.g. icebergs). 488 \item \mdl{daymod} 489 490 No ocean restarts are read or written (though the ice model restarts are retained), so calls to restart functions 457 \item \mdl{daymod} : No ocean restarts are read or written (though the ice model restarts are retained), so calls to restart functions 491 458 have been removed. This also means that the calendar cannot be controlled by time in a restart file, so the user 492 459 must make sure that nn{\_}date0 in the model namelist is correct for his or her purposes. 493 \item \mdl{stpctl} 494 495 Since there is no free surface solver, references to it have been removed from \rou{stp\_ctl} module. 496 \item \mdl{diawri} 497 498 All 3D data have been removed from the output. The surface temperature, salinity and velocity components (which 460 \item \mdl{stpctl} : Since there is no free surface solver, references to it have been removed from \rou{stp\_ctl} module. 461 \item \mdl{diawri} : All 3D data have been removed from the output. The surface temperature, salinity and velocity components (which 499 462 have been read in) are written along with relevant forcing and ice data. 500 \end{ enumerate}463 \end{itemize} 501 464 502 465 One new routine has been added: 503 466 504 \begin{enumerate} 505 \item \mdl{sbcsas} 506 This module initialises the input files needed for reading temperature, salinity and velocity arrays at the surface. 467 \begin{itemize} 468 \item \mdl{sbcsas} : This module initialises the input files needed for reading temperature, salinity and velocity arrays at the surface. 507 469 These filenames are supplied in namelist namsbc{\_}sas. Unfortunately because of limitations with the \mdl{iom} module, 508 470 the full 3D fields from the mean files have to be read in and interpolated in time, before using just the top level. 509 471 Since fldread is used to read in the data, Interpolation on the Fly may be used to change input data resolution. 510 \end{enumerate} 472 \end{itemize} 473 474 475 % Missing the description of the 2 following variables: 476 % ln_3d_uve = .true. ! specify whether we are supplying a 3D u,v and e3 field 477 % ln_read_frq = .false. ! specify whether we must read frq or not 478 479 511 480 512 481 % ================================================================ … … 719 688 are sent to the atmospheric component. 720 689 721 A generalised coupled interface has been developed. It is currently interfaced with OASIS 3 722 (\key{oasis3}) and does not support OASIS 4 723 \footnote{The \key{oasis4} exist. It activates portion of the code that are still under development.}. 690 A generalised coupled interface has been developed. 691 It is currently interfaced with OASIS-3-MCT (\key{oasis3}). 724 692 It has been successfully used to interface \NEMO to most of the European atmospheric 725 693 GCM (ARPEGE, ECHAM, ECMWF, HadAM, HadGAM, LMDz), … … 786 754 \label{SBC_tide} 787 755 788 A module is available to use the tidal potential forcing and is activated with with \key{tide}. 789 790 791 %------------------------------------------nam_tide---------------------------------------------------- 756 %------------------------------------------nam_tide--------------------------------------- 792 757 \namdisplay{nam_tide} 793 %------------------------------------------------------------------------------------------------------------- 794 795 Concerning the tidal potential, some parameters are available in namelist \ngn{nam\_tide}: 758 %----------------------------------------------------------------------------------------- 759 760 A module is available to compute the tidal potential and use it in the momentum equation. 761 This option is activated when \key{tide} is defined. 762 763 Some parameters are available in namelist \ngn{nam\_tide}: 796 764 797 765 - \np{ln\_tide\_pot} activate the tidal potential forcing … … 800 768 801 769 - \np{clname} is the name of constituent 802 803 770 804 771 The tide is generated by the forces of gravity ot the Earth-Moon and Earth-Sun sytem; … … 960 927 \begin{description} 961 928 \item[\np{nn\_isf}~=~1] 962 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, 963 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). 929 The ice shelf cavity is represented. The fwf and heat flux are computed. 2 bulk formulations are available: 930 the ISOMIP one (\np{nn\_isfblk = 1}) described in (\np{nn\_isfblk = 2}), 931 the 3 equation formulation described in \citet{Jenkins1991}. 932 In addition to this, 3 different ways to compute the exchange coefficient are available. 933 $\gamma\_{T/S}$ is constant (\np{nn\_gammablk = 0}), $\gamma\_{T/S}$ is velocity dependant 934 \citep{Jenkins2010} (\np{nn\_gammablk = 1}) and $\gamma\_{T/S}$ is velocity dependant 935 and stratification dependent \citep{Holland1999} (\np{nn\_gammablk = 2}). 936 For each of them, the thermal/salt exchange coefficient (\np{rn\_gammat0} and \np{rn\_gammas0}) 937 have to be specified (the default values are for the ISOMIP case). 964 938 Full description, sensitivity and validation in preparation. 965 939 … … 969 943 (\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). 970 944 Furthermore the fwf is computed using the \citet{Beckmann2003} parameterisation of isf melting. 971 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}). 945 The effective melting length (\np{sn\_Leff\_isf}) is read from a file and the exchange coefficients 946 are set as (\np{rn\_gammat0}) and (\np{rn\_gammas0}). 972 947 973 948 \item[\np{nn\_isf}~=~3] … … 1032 1007 % Handling of icebergs 1033 1008 % ================================================================ 1034 \section{ Handling of icebergs (ICB)}1009 \section{Handling of icebergs (ICB)} 1035 1010 \label{ICB_icebergs} 1036 1011 %------------------------------------------namberg---------------------------------------------------- … … 1038 1013 %------------------------------------------------------------------------------------------------------------- 1039 1014 1040 Icebergs are modelled as lagrangian particles in NEMO. 1041 Their physical behaviour is controlled by equations as described in \citet{Martin_Adcroft_OM10} ). 1042 (Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO.) 1043 Icebergs are initially spawned into one of ten classes which have specific mass and thickness as described in the \ngn{namberg} namelist: 1015 Icebergs are modelled as lagrangian particles in NEMO \citep{Marsh_GMD2015}. 1016 Their physical behaviour is controlled by equations as described in \citet{Martin_Adcroft_OM10} ). 1017 (Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO). 1018 Icebergs are initially spawned into one of ten classes which have specific mass and thickness as described 1019 in the \ngn{namberg} namelist: 1044 1020 \np{rn\_initial\_mass} and \np{rn\_initial\_thickness}. 1045 1021 Each class has an associated scaling (\np{rn\_mass\_scaling}), which is an integer representing how many icebergs … … 1225 1201 The presence at the sea surface of an ice covered area modifies all the fluxes 1226 1202 transmitted to the ocean. There are several way to handle sea-ice in the system 1227 depending on the value of the \np{nn {\_}ice} namelist parameter.1203 depending on the value of the \np{nn\_ice} namelist parameter found in \ngn{namsbc} namelist. 1228 1204 \begin{description} 1229 1205 \item[nn{\_}ice = 0] there will never be sea-ice in the computational domain. … … 1300 1276 % ------------------------------------------------------------------------------------------------------------- 1301 1277 \subsection [Neutral drag coefficient from external wave model (\textit{sbcwave})] 1302 1278 {Neutral drag coefficient from external wave model (\mdl{sbcwave})} 1303 1279 \label{SBC_wave} 1304 1280 %------------------------------------------namwave---------------------------------------------------- 1305 1281 \namdisplay{namsbc_wave} 1306 1282 %------------------------------------------------------------------------------------------------------------- 1307 \begin{description} 1308 1309 \item [??] In order to read a neutral drag coeff, from an external data source (i.e. a wave model), the 1310 logical variable \np{ln\_cdgw} 1311 in $namsbc$ namelist must be defined ${.true.}$. 1283 1284 In order to read a neutral drag coeff, from an external data source ($i.e.$ a wave model), the 1285 logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \textit{true}. 1312 1286 The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the 1313 1287 namelist \ngn{namsbc\_wave} (for external data names, locations, frequency, interpolation and all 1314 1288 the miscellanous options allowed by Input Data generic Interface see \S\ref{SBC_input}) 1315 and a 2D field of neutral drag coefficient. Then using the routine 1316 TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided, the drag coefficient is computed according 1317 to stable/unstable conditions of the air-sea interface following \citet{Large_Yeager_Rep04}. 1318 1319 \end{description} 1289 and a 2D field of neutral drag coefficient. 1290 Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided, 1291 the drag coefficient is computed according to stable/unstable conditions of the air-sea interface following \citet{Large_Yeager_Rep04}. 1292 1320 1293 1321 1294 % Griffies doc: 1322 % When running ocean-ice 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 ocean-ice 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 ocean-ice system over a chosen time scale. 1323 %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 ocean-ice 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. 1324 %When running ocean-ice 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 ocean-ice models, not the water in any one sub-component. As an extreme example to illustrate the issue, consider an ocean-ice model with zero initial sea ice. As the ocean-ice 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 ocean-ice models. 1325 1326 1295 % When running ocean-ice simulations, we are not explicitly representing land processes, 1296 % such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift, 1297 % it is important to balance the hydrological cycle in ocean-ice models. 1298 % We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff. 1299 % The result of the normalization should be a global integrated zero net water input to the ocean-ice system over 1300 % a chosen time scale. 1301 %How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step, 1302 % so that there is always a zero net input of water to the ocean-ice system. 1303 % Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used 1304 % to alter the subsequent year�s water budget in an attempt to damp the annual water imbalance. 1305 % Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing. 1306 % When running ocean-ice coupled models, it is incorrect to include the water transport between the ocean 1307 % and ice models when aiming to balance the hydrological cycle. 1308 % The reason is that it is the sum of the water in the ocean plus ice that should be balanced when running ocean-ice models, 1309 % not the water in any one sub-component. As an extreme example to illustrate the issue, 1310 % consider an ocean-ice model with zero initial sea ice. As the ocean-ice model spins up, 1311 % there should be a net accumulation of water in the growing sea ice, and thus a net loss of water from the ocean. 1312 % The total water contained in the ocean plus ice system is constant, but there is an exchange of water between 1313 % the subcomponents. This exchange should not be part of the normalization used to balance the hydrological cycle 1314 % in ocean-ice models. 1315 1316 -
trunk/DOC/TexFiles/Chapters/Chap_STO.tex
r5404 r6289 10 10 \newpage 11 11 $\ $\newline % force a new line 12 %---------------------------------------namsbc-------------------------------------------------- 13 \namdisplay{namsto} 14 %-------------------------------------------------------------------------------------------------------------- 15 $\ $\newline % force a new ligne 16 17 18 See \cite{Brankart_OM2013} and \cite{Brankart_al_GMD2015} papers for a description of the parameterization. -
trunk/DOC/TexFiles/Chapters/Chap_TRA.tex
r6140 r6289 36 36 (BBL) parametrisation, and an internal damping (DMP) term. The terms QSR, 37 37 BBC, BBL and DMP are optional. The external forcings and parameterisations 38 require complex inputs and complex calculations ( e.g.bulk formulae, estimation38 require complex inputs and complex calculations ($e.g.$ bulk formulae, estimation 39 39 of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and 40 40 described in chapters \S\ref{SBC}, \S\ref{LDF} and \S\ref{ZDF}, respectively. 41 Note that \mdl{tranpc}, the non-penetrative convection module, 41 Note that \mdl{tranpc}, the non-penetrative convection module, although 42 42 located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields, 43 43 is described with the model vertical physics (ZDF) together with other available … … 45 45 46 46 In the present chapter we also describe the diagnostic equations used to compute 47 the sea-water properties (density, Brunt-V ais\"{a}l\"{a} frequency, specific heat and47 the sea-water properties (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and 48 48 freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}). 49 49 … … 54 54 found in the \textit{traTTT} or \textit{traTTT\_xxx} module, in the NEMO/OPA/TRA directory. 55 55 56 The user has the option of extracting each tendency term on the rhsof the tracer57 equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{ MISC}.56 The user has the option of extracting each tendency term on the RHS of the tracer 57 equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{DIA}. 58 58 59 59 $\ $\newline % force a new ligne … … 68 68 %------------------------------------------------------------------------------------------------------------- 69 69 70 The advection tendency of a tracer in flux form is the divergence of the advective 71 fluxes. Its discrete expression is given by : 70 When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \textit{true}), 71 the advection tendency of a tracer is expressed in flux form, 72 $i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by : 72 73 \begin{equation} \label{Eq_tra_adv} 73 74 ADV_\tau =-\frac{1}{b_t} \left( … … 171 172 % 2nd and 4th order centred schemes 172 173 % ------------------------------------------------------------------------------------------------------------- 173 \subsection [ centred schemes (CEN) (\np{ln\_traadv\_cen})]174 { centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}174 \subsection [Centred schemes (CEN) (\np{ln\_traadv\_cen})] 175 {Centred schemes (CEN) (\np{ln\_traadv\_cen}=true)} 175 176 \label{TRA_adv_cen} 176 177 … … 278 279 but on the latter, a split-explicit time stepping is used, with a number of sub-timestep equals 279 280 to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited 280 by vertical advection \citep{Lemarie_OM2015 )}. Note that in this case, a similar split-explicit281 by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar split-explicit 281 282 time stepping should be used on vertical advection of momentum to insure a better stability 282 283 (see \S\ref{DYN_zad}). … … 405 406 Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} 406 407 is used when \np{ln\_traadv\_qck}~=~\textit{true}. 407 QUICKEST implementation can be found in the \mdl{traadv\_ mus} module.408 QUICKEST implementation can be found in the \mdl{traadv\_qck} module. 408 409 409 410 QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST … … 415 416 direction where the control of artificial diapycnal fluxes is of paramount importance. 416 417 Therefore the vertical flux is evaluated using the CEN2 scheme. 417 This no longer guarantees the positivity of the scheme. The use of TVD in the vertical 418 direction (as for the UBS case) should be implemented to restore this property. 418 This no longer guarantees the positivity of the scheme. 419 The use of FCT in the vertical direction (as for the UBS case) should be implemented 420 to restore this property. 419 421 420 422 %%%gmcomment : Cross term are missing in the current implementation.... … … 431 433 %------------------------------------------------------------------------------------------------------------- 432 434 433 Options are defined through the \ngn{namtra\_ldf} namelist variables. 434 The options available for lateral diffusion are a laplacian (rotated or not) 435 or a biharmonic operator, the latter being more scale-selective (more 436 diffusive at small scales). The specification of eddy diffusivity 437 coefficients (either constant or variable in space and time) as well as the 438 computation of the slope along which the operators act, are performed in the 439 \mdl{ldftra} and \mdl{ldfslp} modules, respectively. This is described in Chap.~\ref{LDF}. 435 Options are defined through the \ngn{namtra\_ldf} namelist variables. 436 They are regrouped in four items, allowing to specify 437 $(i)$ the type of operator used (none, laplacian, bilaplacian), 438 $(ii)$ the direction along which the operator acts (iso-level, horizontal, iso-neutral), 439 $(iii)$ some specific options related to the rotated operators ($i.e.$ non-iso-level operator), and 440 $(iv)$ the specification of eddy diffusivity coefficient (either constant or variable in space and time). 441 Item $(iv)$ will be described in Chap.\ref{LDF} . 442 The direction along which the operators act is defined through the slope between this direction and the iso-level surfaces. 443 The slope is computed in the \mdl{ldfslp} module and will also be described in Chap.~\ref{LDF}. 444 440 445 The lateral diffusion of tracers is evaluated using a forward scheme, 441 446 $i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time, 442 except for the pure vertical component that appears when a rotation tensor 443 is used. This latter term is solved implicitly together with the 444 vertical diffusion term (see \S\ref{STP}). 445 446 % ------------------------------------------------------------------------------------------------------------- 447 % Iso-level laplacian operator 448 % ------------------------------------------------------------------------------------------------------------- 449 \subsection [Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap})] 450 {Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap}=true) } 451 \label{TRA_ldf_lap} 452 453 A laplacian diffusion operator ($i.e.$ a harmonic operator) acting along the model 454 surfaces is given by: 447 except for the pure vertical component that appears when a rotation tensor is used. 448 This latter component is solved implicitly together with the vertical diffusion term (see \S\ref{STP}). 449 When \np{ln\_traldf\_msc}~=~\textit{true}, a Method of Stabilizing Correction is used in which 450 the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}. 451 452 % ------------------------------------------------------------------------------------------------------------- 453 % Type of operator 454 % ------------------------------------------------------------------------------------------------------------- 455 \subsection [Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, \np{ln\_traldf\_blp})] 456 {Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, or \np{ln\_traldf\_blp} = true) } 457 \label{TRA_ldf_op} 458 459 Three operator options are proposed and, one and only one of them must be selected: 460 \begin{description} 461 \item [\np{ln\_traldf\_NONE}] = true : no operator selected, the lateral diffusive tendency will not be 462 applied to the tracer equation. This option can be used when the selected advection scheme 463 is diffusive enough (MUSCL scheme for example). 464 \item [ \np{ln\_traldf\_lap}] = true : a laplacian operator is selected. This harmonic operator 465 takes the following expression: $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $, 466 where the gradient operates along the selected direction (see \S\ref{TRA_ldf_dir}), 467 and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see Chap.~\ref{LDF}). 468 \item [\np{ln\_traldf\_blp}] = true : a bilaplacian operator is selected. This biharmonic operator 469 takes the following expression: 470 $\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 471 where the gradient operats along the selected direction, 472 and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see Chap.~\ref{LDF}). 473 In the code, the bilaplacian operator is obtained by calling the laplacian twice. 474 \end{description} 475 476 Both laplacian and bilaplacian operators ensure the total tracer variance decrease. 477 Their primary role is to provide strong dissipation at the smallest scale supported 478 by the grid while minimizing the impact on the larger scale features. 479 The main difference between the two operators is the scale selectiveness. 480 The bilaplacian damping time ($i.e.$ its spin down time) scales like $\lambda^{-4}$ 481 for disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones), 482 whereas the laplacian damping time scales only like $\lambda^{-2}$. 483 484 485 % ------------------------------------------------------------------------------------------------------------- 486 % Direction of action 487 % ------------------------------------------------------------------------------------------------------------- 488 \subsection [Direction of action (\np{ln\_traldf\_lev}, \np{ln\_traldf\_hor}, \np{ln\_traldf\_iso}, \np{ln\_traldf\_triad})] 489 {Direction of action (\np{ln\_traldf\_lev}, \textit{...\_hor}, \textit{...\_iso}, or \textit{...\_triad} = true) } 490 \label{TRA_ldf_dir} 491 492 The choice of a direction of action determines the form of operator used. 493 The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane 494 when iso-level option is used (\np{ln\_traldf\_lev}~=~\textit{true}) 495 or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}-coordinate 496 (\np{ln\_traldf\_hor} and \np{ln\_zco} equal \textit{true}). 497 The associated code can be found in the \mdl{traldf\_lap\_blp} module. 498 The operator is a rotated (re-entrant) laplacian when the direction along which it acts 499 does not coincide with the iso-level surfaces, 500 that is when standard or triad iso-neutral option is used (\np{ln\_traldf\_iso} or 501 \np{ln\_traldf\_triad} equals \textit{true}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.), 502 or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}-coordinate 503 (\np{ln\_traldf\_hor} and \np{ln\_sco} equal \textit{true}) 504 \footnote{In this case, the standard iso-neutral operator will be automatically selected}. 505 In that case, a rotation is applied to the gradient(s) that appears in the operator 506 so that diffusive fluxes acts on the three spatial direction. 507 508 The resulting discret form of the three operators (one iso-level and two rotated one) 509 is given in the next two sub-sections. 510 511 512 % ------------------------------------------------------------------------------------------------------------- 513 % iso-level operator 514 % ------------------------------------------------------------------------------------------------------------- 515 \subsection [Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso})] 516 {Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso}) } 517 \label{TRA_ldf_lev} 518 519 The laplacian diffusion operator acting along the model (\textit{i,j})-surfaces is given by: 455 520 \begin{equation} \label{Eq_tra_ldf_lap} 456 D_ T^{lT} =\frac{1}{b_tT} \left( \;521 D_t^{lT} =\frac{1}{b_t} \left( \; 457 522 \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right] 458 523 + \delta _{j}\left[ A_v^{lT} \; \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right] \;\right) 459 524 \end{equation} 460 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells. 461 It is implemented in the \mdl{traadv\_lap} module. 462 463 This lateral operator is computed in \mdl{traldf\_lap}. It is a \emph{horizontal} 464 operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate with 465 or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 466 It is thus used when, in addition to \np{ln\_traldf\_lap}=true, we have 467 \np{ln\_traldf\_level}=true or \np{ln\_traldf\_hor}=\np{ln\_zco}=true. 525 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells 526 and where zero diffusive fluxes is assumed across solid boundaries, 527 first (and third in bilaplacian case) horizontal tracer derivative are masked. 528 It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module. 529 The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} 530 in order to compute the iso-level bilaplacian operator. 531 532 It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate 533 with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 534 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~=~\textit{true}, 535 we have \np{ln\_traldf\_lev}~=~\textit{true} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~=~\textit{true}. 468 536 In both cases, it significantly contributes to diapycnal mixing. 469 It is therefore n ot recommended.537 It is therefore never recommended, even when using it in the bilaplacian case. 470 538 471 539 Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), tracers in horizontally … … 475 543 described in \S\ref{TRA_zpshde}. 476 544 477 % ------------------------------------------------------------------------------------------------------------- 478 % Rotated laplacian operator 479 % ------------------------------------------------------------------------------------------------------------- 480 \subsection [Rotated laplacian operator (iso) (\np{ln\_traldf\_lap})] 481 {Rotated laplacian operator (iso) (\np{ln\_traldf\_lap}=true)} 545 546 % ------------------------------------------------------------------------------------------------------------- 547 % Rotated laplacian operator 548 % ------------------------------------------------------------------------------------------------------------- 549 \subsection [Standard and triad rotated (bi-)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad})] 550 {Standard and triad (bi-)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad}))} 551 \label{TRA_ldf_iso_triad} 552 553 %&& Standard rotated (bi-)laplacian operator 554 %&& ---------------------------------------------- 555 \subsubsection [Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})] 556 {Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})} 482 557 \label{TRA_ldf_iso} 483 484 If the Griffies trad scheme is not employed 485 (\np{ln\_traldf\_grif}=true; see App.\ref{sec:triad}) the general form of the second order lateral tracer subgrid scale physics 486 (\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and 487 $s$-coordinates: 558 The general form of the second order lateral tracer subgrid scale physics 559 (\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and $s$-coordinates: 488 560 \begin{equation} \label{Eq_tra_ldf_iso} 489 561 \begin{split} … … 527 599 of the tracer variance. Nevertheless the treatment performed on the slopes 528 600 (see \S\ref{LDF}) allows the model to run safely without any additional 529 background horizontal diffusion \citep{Guilyardi_al_CD01}. An alternative scheme 530 developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 601 background horizontal diffusion \citep{Guilyardi_al_CD01}. 602 603 Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal derivatives 604 at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific treatment. 605 They are calculated in module zpshde, described in \S\ref{TRA_zpshde}. 606 607 %&& Triad rotated (bi-)laplacian operator 608 %&& ------------------------------------------- 609 \subsubsection [Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})] 610 {Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})} 611 \label{TRA_ldf_triad} 612 613 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}=true ; see App.\ref{sec:triad}) 614 615 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 531 616 is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of 532 617 the algorithm is given in App.\ref{sec:triad}. 533 534 Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal535 derivatives at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific536 treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.537 538 % -------------------------------------------------------------------------------------------------------------539 % Iso-level bilaplacian operator540 % -------------------------------------------------------------------------------------------------------------541 \subsection [Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap})]542 {Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap}=true)}543 \label{TRA_ldf_bilap}544 618 545 619 The lateral fourth order bilaplacian operator on tracers is obtained by 546 620 applying (\ref{Eq_tra_ldf_lap}) twice. The operator requires an additional assumption 547 621 on boundary conditions: both first and third derivative terms normal to the 548 coast are set to zero. It is used when, in addition to \np{ln\_traldf\_bilap}=true, 549 we have \np{ln\_traldf\_level}=true, or both \np{ln\_traldf\_hor}=true and 550 \np{ln\_zco}=false. In both cases, it can contribute diapycnal mixing, 551 although less than in the laplacian case. It is therefore not recommended. 552 553 Note that in the code, the bilaplacian routine does not call the laplacian 554 routine twice but is rather a separate routine that can be found in the 555 \mdl{traldf\_bilap} module. This is due to the fact that we introduce the 556 eddy diffusivity coefficient, A, in the operator as: 557 $\nabla \cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$, 558 instead of 559 $-\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$ 560 where $a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations 561 ensure the total variance decrease, but the former requires a larger 562 number of code-lines. 563 564 % ------------------------------------------------------------------------------------------------------------- 565 % Rotated bilaplacian operator 566 % ------------------------------------------------------------------------------------------------------------- 567 \subsection [Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap})] 568 {Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap}=true)} 569 \label{TRA_ldf_bilapg} 622 coast are set to zero. 570 623 571 624 The lateral fourth order operator formulation on tracers is obtained by 572 625 applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption 573 626 on boundary conditions: first and third derivative terms normal to the 574 coast, normal to the bottom and normal to the surface are set to zero. It can be found in the 575 \mdl{traldf\_bilapg}. 576 577 It is used when, in addition to \np{ln\_traldf\_bilap}=true, we have 578 \np{ln\_traldf\_iso}= .true, or both \np{ln\_traldf\_hor}=true and \np{ln\_zco}=true. 579 This rotated bilaplacian operator has never been seriously 580 tested. There are no guarantees that it is either free of bugs or correctly formulated. 581 Moreover, the stability range of such an operator will be probably quite 582 narrow, requiring a significantly smaller time-step than the one used with an 583 unrotated operator. 627 coast, normal to the bottom and normal to the surface are set to zero. 628 629 %&& Option for the rotated operators 630 %&& ---------------------------------------------- 631 \subsubsection [Option for the rotated operators] 632 {Option for the rotated operators} 633 \label{TRA_ldf_options} 634 635 \np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators) 636 637 \np{rn\_slpmax} = slope limit (both operators) 638 639 \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 640 641 \np{rn\_sw\_triad} =1 switching triad ; =0 all 4 triads used (triad only) 642 643 \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only) 584 644 585 645 % ================================================================ … … 593 653 %-------------------------------------------------------------------------------------------------------------- 594 654 595 Options are defined through the 655 Options are defined through the \ngn{namzdf} namelist variables. 596 656 The formulation of the vertical subgrid scale tracer physics is the same 597 657 for all the vertical coordinates, and is based on a laplacian operator. … … 651 711 the thickness of the top model layer. 652 712 653 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components ($i.e.$ atmosphere, sea-ice, land), 654 the change in the heat and salt content of the surface layer of the ocean is due both 655 to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) 656 and to the heat and salt content of the mass exchange. 657 \sgacomment{ the following does not apply to the release to which this documentation is 658 attached and so should not be included .... 659 In a forthcoming release, these two parts, computed in the surface module (SBC), will be included directly 660 in $Q_{ns}$, the surface heat flux and $F_{salt}$, the surface salt flux. 661 The specification of these fluxes is further detailed in the SBC chapter (see \S\ref{SBC}). 662 This change will provide a forcing formulation which is the same for any tracer (including temperature and salinity). 663 664 In the current version, the situation is a little bit more complicated. } 713 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components 714 ($i.e.$ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer 715 of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) 716 and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$, 717 the surface heat flux, and $F_{salt}$, the surface salt flux (see \S\ref{SBC} for further details). 718 By doing this, the forcing formulation is the same for any tracer (including temperature and salinity). 665 719 666 720 The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following … … 669 723 $\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 670 724 (i.e. the difference between the total surface heat flux and the fraction of the short wave flux that 671 penetrates into the water column, see \S\ref{TRA_qsr}) 672 673 $\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) 674 675 $\bullet$ $\textit{emp}_S$, an equivalent mass flux taking into account the effect of ice-ocean mass exchange 676 677 $\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) 678 679 The $\textit{emp}_S$ field is not simply the budget of evaporation-precipitation+freezing-melting because 680 the sea-ice is not currently embedded in the ocean but levitates above it. There is no mass 681 exchanged between the sea-ice and the ocean. Instead we only take into account the salt 682 flux associated with the non-zero salinity of sea-ice, and the concentration/dilution effect 683 due to the freezing/melting (F/M) process. These two parts of the forcing are then converted into 684 an equivalent mass flux given by $\textit{emp}_S - \textit{emp}$. As a result of this mess, 685 the surface boundary condition on temperature and salinity is applied as follows: 686 687 In the nonlinear free surface case (\key{vvl} is defined): 725 penetrates into the water column, see \S\ref{TRA_qsr}) plus the heat content associated with 726 of the mass exchange with the atmosphere and lands. 727 728 $\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 729 730 $\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) 731 and possibly with the sea-ice and ice-shelves. 732 733 $\bullet$ \textit{rnf}, the mass flux associated with runoff 734 (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 735 736 The surface boundary condition on temperature and salinity is applied as follows: 688 737 \begin{equation} \label{Eq_tra_sbc} 738 \begin{aligned} 739 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} } &\overline{ Q_{ns} }^t & \\ 740 & F^S =\frac{ 1 }{\rho _o \, \left. e_{3t} \right|_{k=1} } &\overline{ \textit{sfx} }^t & \\ 741 \end{aligned} 742 \end{equation} 743 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps 744 ($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the 745 divergence of odd and even time step (see \S\ref{STP}). 746 747 In the linear free surface case (\np{ln\_linssh}~=~\textit{true}), 748 an additional term has to be added on both temperature and salinity. 749 On temperature, this term remove the heat content associated with mass exchange 750 that has been added to $Q_{ns}$. On salinity, this term mimics the concentration/dilution effect that 751 would have resulted from a change in the volume of the first level. 752 The resulting surface boundary condition is applied as follows: 753 \begin{equation} \label{Eq_tra_sbc_lin} 689 754 \begin{aligned} 690 755 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} } … … 692 757 % 693 758 & F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} } 694 &\overline{ \left( (\textit{emp}_S - \textit{emp})\;\left. S \right|_{k=1} \right) }^t & \\759 &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1} \right) }^t & \\ 695 760 \end{aligned} 696 761 \end{equation} 697 698 In the linear free surface case (\key{vvl} not defined): 699 \begin{equation} \label{Eq_tra_sbc_lin} 700 \begin{aligned} 701 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} } &\overline{ Q_{ns} }^t & \\ 702 % 703 & F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} } 704 &\overline{ \left( \textit{emp}_S\;\left. S \right|_{k=1} \right) }^t & \\ 705 \end{aligned} 706 \end{equation} 707 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps 708 ($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the 709 divergence of odd and even time step (see \S\ref{STP}). 710 711 The two set of equations, \eqref{Eq_tra_sbc} and \eqref{Eq_tra_sbc_lin}, are obtained 712 by assuming that the temperature of precipitation and evaporation are equal to 713 the ocean surface temperature and that their salinity is zero. Therefore, the heat content 714 of the \textit{emp} budget must be added to the temperature equation in the variable volume case, 715 while it does not appear in the constant volume case. Similarly, the \textit{emp} budget affects 716 the ocean surface salinity in the constant volume case (through the concentration dilution effect) 717 while it does not appears explicitly in the variable volume case since salinity change will be 718 induced by volume change. In both constant and variable volume cases, surface salinity 719 will change with ice-ocean salt flux and F/M flux (both contained in $\textit{emp}_S - \textit{emp}$) without mass exchanges. 720 721 Note that the concentration/dilution effect due to F/M is computed using 722 a constant ice salinity as well as a constant ocean salinity. 723 This approximation suppresses the correlation between \textit{SSS} 724 and F/M flux, allowing the ice-ocean salt exchanges to be conservative. 725 Indeed, if this approximation is not made, even if the F/M budget is zero 726 on average over the whole ocean domain and over the seasonal cycle, 727 the associated salt flux is not zero, since sea-surface salinity and F/M flux are 728 intrinsically correlated (high \textit{SSS} are found where freezing is 729 strong whilst low \textit{SSS} is usually associated with high melting areas). 730 731 Even using this approximation, an exact conservation of heat and salt content 732 is only achieved in the variable volume case. In the constant volume case, 733 there is a small imbalance associated with the product $(\partial_t\eta - \textit{emp}) * \textit{SSS}$. 734 Nevertheless, the salt content variation is quite small and will not induce 735 a long term drift as there is no physical reason for $(\partial_t\eta - \textit{emp})$ 736 and \textit{SSS} to be correlated \citep{Roullet_Madec_JGR00}. 737 Note that, while quite small, the imbalance in the constant volume case is larger 762 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface. 763 In the linear free surface case, there is a small imbalance. The imbalance is larger 738 764 than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}. 739 This is the reason why the modified filter is not applied in the constant volume case.765 This is the reason why the modified filter is not applied in the linear free surface case (see \S\ref{STP}). 740 766 741 767 % ------------------------------------------------------------------------------------------------------------- … … 849 875 \label{TRA_bbc} 850 876 %--------------------------------------------nambbc-------------------------------------------------------- 851 \namdisplay{nam tra_bbc}877 \namdisplay{nambbc} 852 878 %-------------------------------------------------------------------------------------------------------------- 853 879 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 1093 1119 \subsection[DMP\_TOOLS]{Generating resto.nc using DMP\_TOOLS} 1094 1120 1095 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. 1121 DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. 1122 Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled 1123 and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input. 1124 This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. 1125 The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. 1126 The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient. 1096 1127 1097 1128 %--------------------------------------------nam_dmp_create------------------------------------------------- 1098 \nam display{nam_dmp_create}1129 \namtools{namelist_dmp} 1099 1130 %------------------------------------------------------------------------------------------------------- 1100 1131 1101 1132 \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. 1102 1133 1103 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom 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. 1104 1105 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}. 1134 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations. 1135 \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. 1136 \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea 1137 for the ORCA4, ORCA2 and ORCA05 configurations. 1138 If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as 1139 a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference 1140 configurations with previous model versions. 1141 \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. 1142 This option only has an effect if \np{ln\_full\_field} is true. 1143 \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. 1144 Finally \np{ln\_custom} specifies that the custom module will be called. 1145 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. 1146 1147 The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}. 1148 Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to 1149 the full values of a 10$^{\circ}$ latitud band. 1150 This is often used because of the short adjustment time scale in the equatorial region 1151 \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a 1152 hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}. 1106 1153 1107 1154 % ================================================================ … … 1271 1318 1272 1319 % ------------------------------------------------------------------------------------------------------------- 1273 % Brunt-V ais\"{a}l\"{a} Frequency1274 % ------------------------------------------------------------------------------------------------------------- 1275 \subsection{Brunt-V ais\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}1320 % Brunt-V\"{a}is\"{a}l\"{a} Frequency 1321 % ------------------------------------------------------------------------------------------------------------- 1322 \subsection{Brunt-V\"{a}is\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)} 1276 1323 \label{TRA_bn2} 1277 1324 1278 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V ais\"{a}l\"{a}1325 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} 1279 1326 frequency) is of paramount importance as determine the ocean stratification and 1280 1327 is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent … … 1332 1379 \label{TRA_zpshde} 1333 1380 1334 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, I've changed "derivative" to "difference" and "mean" to "average"} 1381 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, 1382 I've changed "derivative" to "difference" and "mean" to "average"} 1335 1383 1336 1384 With partial bottom cells (\np{ln\_zps}=true), in general, tracers in horizontally -
trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex
r5120 r6289 34 34 coefficients can be assumed to be either constant, or a function of the local 35 35 Richardson number, or computed from a turbulent closure model (either 36 TKE or KPPformulation). The computation of these coefficients is initialized36 TKE or GLS formulation). The computation of these coefficients is initialized 37 37 in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or 38 \mdl{zdf kpp} modules. The trends due to the vertical momentum and tracer38 \mdl{zdfgls} modules. The trends due to the vertical momentum and tracer 39 39 diffusion, including the surface forcing, are computed and added to the 40 40 general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. … … 355 355 %--------------------------------------------------------------% 356 356 357 To be add here a description of "penetration of TKE" and the associated namelist parameters 358 \np{nn\_etau}, \np{rn\_efr} and \np{nn\_htau}. 357 Vertical mixing parameterizations commonly used in ocean general circulation models 358 tend to produce mixed-layer depths that are too shallow during summer months and windy conditions. 359 This bias is particularly acute over the Southern Ocean. 360 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{Rodgers_2014}. 361 The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations, 362 but rather is meant to account for observed processes that affect the density structure of 363 the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme 364 ($i.e.$ near-inertial oscillations and ocean swells and waves). 365 366 When using this parameterization ($i.e.$ when \np{nn\_etau}~=~1), the TKE input to the ocean ($S$) 367 imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized 368 by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by: 369 \begin{equation} \label{ZDF_Ehtau} 370 S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau} 371 \end{equation} 372 where 373 $z$ is the depth, 374 $e_s$ is TKE surface boundary condition, 375 $f_r$ is the fraction of the surface TKE that penetrate in the ocean, 376 $h_\tau$ is a vertical mixing length scale that controls exponential shape of the penetration, 377 and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely 378 covered by sea-ice). 379 The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. 380 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}~=~0) 381 or a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m 382 at high latitudes (\np{nn\_etau}~=~1). 383 384 Note that two other option existe, \np{nn\_etau}~=~2, or 3. They correspond to applying 385 \eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part 386 of the stress to evaluate the fraction of TKE that penetrate the ocean. 387 Those two options are obsolescent features introduced for test purposes. 388 They will be removed in the next release. 389 390 359 391 360 392 % from Burchard et al OM 2008 : 361 % 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). 393 % the most critical process not reproduced by statistical turbulence models is the activity of 394 % internal waves and their interaction with turbulence. After the Reynolds decomposition, 395 % internal waves are in principle included in the RANS equations, but later partially 396 % excluded by the hydrostatic assumption and the model resolution. 397 % Thus far, the representation of internal wave mixing in ocean models has been relatively crude 398 % (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 362 399 363 400 … … 573 610 Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}. 574 611 575 % -------------------------------------------------------------------------------------------------------------576 % K Profile Parametrisation (KPP)577 % -------------------------------------------------------------------------------------------------------------578 \subsection{K Profile Parametrisation (KPP) (\key{zdfkpp}) }579 \label{ZDF_kpp}580 581 %--------------------------------------------namkpp--------------------------------------------------------582 \namdisplay{namzdf_kpp}583 %--------------------------------------------------------------------------------------------------------------584 585 The KKP scheme has been implemented by J. Chanut ...586 Options are defined through the \ngn{namzdf\_kpp} namelist variables.587 588 \colorbox{yellow}{Add a description of KPP here.}589 590 612 591 613 % ================================================================ … … 636 658 637 659 Options are defined through the \ngn{namzdf} namelist variables. 638 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc} =true.660 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}~=~\textit{true}. 639 661 It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously 640 662 the statically unstable portion of the water column, but only until the density … … 644 666 (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is 645 667 found. Assume in the following that the instability is located between levels 646 $k$ and $k+1$. The potentialtemperature and salinity in the two levels are668 $k$ and $k+1$. The temperature and salinity in the two levels are 647 669 vertically mixed, conserving the heat and salt contents of the water column. 648 670 The new density is then computed by a linear approximation. If the new … … 664 686 \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}. 665 687 666 Note that in the current implementation of this algorithm presents several 667 limitations. First, potential density referenced to the sea surface is used to 668 check whether the density profile is stable or not. This is a strong 669 simplification which leads to large errors for realistic ocean simulations. 670 Indeed, many water masses of the world ocean, especially Antarctic Bottom 671 Water, are unstable when represented in surface-referenced potential density. 672 The scheme will erroneously mix them up. Second, the mixing of potential 673 density is assumed to be linear. This assures the convergence of the algorithm 674 even when the equation of state is non-linear. Small static instabilities can thus 675 persist due to cabbeling: they will be treated at the next time step. 676 Third, temperature and salinity, and thus density, are mixed, but the 677 corresponding velocity fields remain unchanged. When using a Richardson 678 Number dependent eddy viscosity, the mixing of momentum is done through 679 the vertical diffusion: after a static adjustment, the Richardson Number is zero 680 and thus the eddy viscosity coefficient is at a maximum. When this convective 681 adjustment algorithm is used with constant vertical eddy viscosity, spurious 682 solutions can occur since the vertical momentum diffusion remains small even 683 after a static adjustment. In that case, we recommend the addition of momentum 684 mixing in a manner that mimics the mixing in temperature and salinity 685 \citep{Speich_PhD92, Speich_al_JPO96}. 688 The current implementation has been modified in order to deal with any non linear 689 equation of seawater (L. Brodeau, personnal communication). 690 Two main differences have been introduced compared to the original algorithm: 691 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 692 (not the the difference in potential density) ; 693 $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients 694 are vertically mixed in the same way their temperature and salinity has been mixed. 695 These two modifications allow the algorithm to perform properly and accurately 696 with TEOS10 or EOS-80 without having to recompute the expansion coefficients at each 697 mixing iteration. 686 698 687 699 % ------------------------------------------------------------------------------------------------------------- … … 689 701 % ------------------------------------------------------------------------------------------------------------- 690 702 \subsection [Enhanced Vertical Diffusion (\np{ln\_zdfevd})] 691 703 {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=true)} 692 704 \label{ZDF_evd} 693 705 -
trunk/DOC/TexFiles/Chapters/Introduction.tex
r6140 r6289 32 32 The equations are written in a curvilinear coordinate system, with a choice of vertical 33 33 coordinates ($z$, $s$, \textit{z*}, \textit{s*}, $\tilde{z}$, $\tilde{s}$, and a mixture of them). 34 Momentum equations are formulated in the vector invariant form or in theflux form.34 Momentum equations are formulated in vector invariant or flux form. 35 35 Dimensional units in the meter, kilogram, second (MKS) international system 36 36 are used throughout.
Note: See TracChangeset
for help on using the changeset viewer.