Changeset 994 for trunk/DOC/TexFiles/Chapters
- Timestamp:
- 2008-05-28T11:01:09+02:00 (16 years ago)
- Location:
- trunk/DOC/TexFiles
- Files:
-
- 15 copied
- 1 moved
-
. (moved) (moved from trunk/DOC/BETA)
-
Chapters/Abstracts_Foreword.tex (copied) (copied from trunk/DOC/BETA/Chapters/Abstracts_Foreword.tex) (1 diff)
-
Chapters/Annex_A.tex (copied) (copied from trunk/DOC/BETA/Chapters/Annex_A.tex)
-
Chapters/Annex_B.tex (copied) (copied from trunk/DOC/BETA/Chapters/Annex_B.tex)
-
Chapters/Annex_C.tex (copied) (copied from trunk/DOC/BETA/Chapters/Annex_C.tex) (10 diffs)
-
Chapters/Annex_D.tex (copied) (copied from trunk/DOC/BETA/Chapters/Annex_D.tex)
-
Chapters/Chap_DOM.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_DOM.tex) (1 diff)
-
Chapters/Chap_DYN.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_DYN.tex) (8 diffs)
-
Chapters/Chap_LBC.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_LBC.tex) (18 diffs)
-
Chapters/Chap_LDF.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_LDF.tex)
-
Chapters/Chap_MISC.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_MISC.tex) (25 diffs)
-
Chapters/Chap_Model_Basics.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_Model_Basics.tex) (2 diffs)
-
Chapters/Chap_SBC.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_SBC.tex) (19 diffs)
-
Chapters/Chap_TRA.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_TRA.tex) (1 diff)
-
Chapters/Chap_ZDF.tex (copied) (copied from trunk/DOC/BETA/Chapters/Chap_ZDF.tex) (26 diffs)
-
Chapters/Introduction.tex (copied) (copied from trunk/DOC/BETA/Chapters/Introduction.tex)
Legend:
- Unmodified
- Added
- Removed
-
trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex
r817 r994 28 28 \vspace{0.5cm} 29 29 30 Madec, G., 2008: NEMO ocean engine. \textit{Note du P\^ole de mod\'{e}lisation}, Institut Pierre-Simon Laplace (IPSL), France, \colorbox{yellow}{NXX, 193pp}.\\ 30 Madec, G., 2008: NEMO ocean engine. \textit{Note du P\^ole de mod\'{e}lisation}, Institut Pierre-Simon Laplace (IPSL), France, No 27, ISSN No 1288-1619.\\ 31 31 32 32 33 \vspace{0.5cm} -
trunk/DOC/TexFiles/Chapters/Annex_C.tex
r817 r994 9 9 I'm writting this appendix. It will be available in a forthcoming release of the documentation 10 10 11 \gmcomment{11 %\gmcomment{ 12 12 13 13 % ================================================================ … … 160 160 and introducing the horizontal divergence $\chi $, it becomes: 161 161 \begin{align*} 162 \equiv& \sum\limits_{i,j,k} - \frac{1} {2} \left( \frac{\zeta} {e_{3f}} \right)^2 \; \overline{\overline{ e_ 1T\,e_2T\,e_3T\, \chi}}^{\,i+1/2,j+1/2} \;\;\equiv 0162 \equiv& \sum\limits_{i,j,k} - \frac{1} {2} \left( \frac{\zeta} {e_{3f}} \right)^2 \; \overline{\overline{ e_{1T}\,e_{2T}\,e_{3T}\, \chi}}^{\,i+1/2,j+1/2} \;\;\equiv 0 163 163 \qquad \qquad \qquad \qquad \qquad \qquad \qquad &&&&\\ 164 164 \end{align*} … … 169 169 the two components of the vorticity term is optionally offered (ENE scheme) : 170 170 \begin{equation*} 171 \frac{1} {e_3 }\nabla \times 172 \left( 173 \zeta \;{\textbf{k}}\times {\textbf {U}}_h 174 \right) 171 - \zeta \;{\textbf{k}}\times {\textbf {U}}_h 175 172 \equiv 176 173 \left( {{ … … 192 189 enstrophy on the relative vorticity term and energy on the Coriolis term. 193 190 \begin{flalign*} 194 &\int\limits_D \textbf{U}_h \times\left( \zeta \;\textbf{k} \times \textbf{U}_h \right) \; dv && \\191 &\int\limits_D - \textbf{U}_h \cdot \left( \zeta \;\textbf{k} \times \textbf{U}_h \right) \; dv && \\ 195 192 \equiv& \sum\limits_{i,j,k} \biggl\{ 196 193 \overline {\left( \frac{\zeta} {e_{3f}} \right) … … 348 345 &&& \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad 349 346 + \delta_k \left[ \overline {e_{1w}\,e_{2w}\,w}^{\,i+12}\;\overline u^{\,k+1/2} \right] 350 \; \biggr\} &&& \displaybreak[0] \\351 % 352 \equiv& \sum\limits_{i,j,k}347 \; \biggr\} \; u &&& \displaybreak[0] \\ 348 % 349 \equiv& - \sum\limits_{i,j,k} 353 350 \biggl\{ 354 351 \overline {e_{2u}\,e_{3u}\,u}^{\,i}\; \overline u^{\,i} \delta_i … … 359 356 + \overline {e_{1w}\,e_{2w}\,w}^{\,i+1/2}\; \overline u^{\,k+1/2} \delta_{k+1/2} \left[ u \right] \biggr\} && \displaybreak[0] \\ 360 357 % 361 \equiv& \sum\limits_{i,j,k}358 \equiv& - \sum\limits_{i,j,k} 362 359 \biggl\{ 363 360 \overline {e_{2u}\,e_{3u}\,u}^{\,i} \delta_i \left[ u^2 \right] … … 1224 1221 constraint of conservation of tracers: 1225 1222 \begin{flalign*} 1226 &\int\limits_D T\;\nabla \cdot \left( A\;\nabla T \right)\;dv &&&\\1223 &\int\limits_D \nabla \cdot \left( A\;\nabla T \right)\;dv &&&\\ 1227 1224 \\ 1228 1225 &\equiv \sum\limits_{i,j,k} … … 1248 1245 \end{flalign*} 1249 1246 1247 In fact, this property is simply resulting from the flux form of the operator. 1250 1248 1251 1249 % ------------------------------------------------------------------------------------------------------------- … … 1257 1255 constraint of dissipation of tracer variance: 1258 1256 \begin{flalign*} 1259 \int\limits_D T\;\nabla \cdot \left( A\;\nabla T \right)\;dv &&&\\ 1260 \end{flalign*} 1261 \begin{flalign*} 1262 \equiv \sum\limits_{i,j,k} T 1263 \biggl\{ \biggr. 1264 \delta_i 1265 \left[ 1266 A_u^{\,lT} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} 1267 \left[ T \right] 1268 \right] 1269 + \delta_j 1270 &\left[ 1271 A_v^{\,lT} \frac{e_{1v}\,e_{3v}} {e_{2v}} \delta_{j+1/2} 1272 \left[ T \right] 1273 \right] 1274 && \\ 1275 \biggl. 1276 + \delta_k 1277 &\left[ 1278 A_w^{\,vT} \frac{e_{1T}\,e_{2T}} {e_{3T}} \delta_{k+1/2} 1279 \left[ T \right] 1280 \right] 1281 \biggr\} 1282 &&\\ 1257 \int\limits_D T\;\nabla & \cdot \left( A\;\nabla T \right)\;dv &&&\\ 1258 &\equiv \sum\limits_{i,j,k} \; T 1259 \biggl\{ \biggr. 1260 \delta_i \left[ A_u^{\,lT} \frac{e_{2u}\,e_{3u}} {e_{1u}} \delta_{i+1/2} \left[T\right] \right] 1261 & + \delta_j \left[ A_v^{\,lT} \frac{e_{1v} \,e_{3v}} {e_{2v}} \delta_{j+1/2} \left[T\right] \right] 1262 \quad&& \\ 1263 \biggl. 1264 &&+ \delta_k \left[A_w^{\,vT}\frac{e_{1T}\,e_{2T}} {e_{3T}}\delta_{k+1/2}\left[T\right]\right] 1265 \biggr\} && 1283 1266 \end{flalign*} 1284 1267 \begin{flalign*} … … 1313 1296 1314 1297 %%%% end of appendix in gm comment 1315 }1298 %} -
trunk/DOC/TexFiles/Chapters/Chap_DOM.tex
r817 r994 819 819 \; 0& \text{ if $k\leq mbathy(i,j)$ } \end{cases} \\ 820 820 umask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 821 umask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i,j+1,k) \\822 umask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\821 vmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i,j+1,k) \\ 822 fmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 823 823 & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) 824 824 \end{align*} -
trunk/DOC/TexFiles/Chapters/Chap_DYN.tex
r817 r994 107 107 \left\{ 108 108 \begin{aligned} 109 { -\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} & {\overline{\overline {\left( {e_{1v} e_{3v} v} \right)}} }^{\,i, j+1/2} \\110 { +\frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j} & {\overline{\overline {\left( {e_{2u} e_{3u} u} \right)}} }^{\,i+1/2, j}109 {+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} & {\overline{\overline {\left( {e_{1v} e_{3v} v} \right)}} }^{\,i, j+1/2} \\ 110 {-\frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j} & {\overline{\overline {\left( {e_{2u} e_{3u} u} \right)}} }^{\,i+1/2, j} 111 111 \end{aligned} 112 112 \right. … … 124 124 \left\{ { 125 125 \begin{aligned} 126 { -\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)126 {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) 127 127 \;\overline {\left( {e_{1v} e_{3v} v} \right)} ^{\,i+1/2}} }^{\,j} } \\ 128 { +\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)128 {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) 129 129 \;\overline {\left( {e_{2u} e_{3u} u} \right)} ^{\,j+1/2}} }^{\,i} } 130 130 \end{aligned} … … 145 145 \left\{ { 146 146 \begin{aligned} 147 { -\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i}147 {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} 148 148 \; {\overline{\overline {\left( {e_{1v} \; e_{3v} \ v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} } 149 149 \; {\overline {\left( {\frac{f}{e_{3f} }} \right) 150 150 \;\overline {\left( {e_{1v} \; e_{3v} \ v} \right)} ^{\,i+1/2}} }^{\,j} } \\ 151 { +\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j151 {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j 152 152 \; {\overline{\overline {\left( {e_{2u} \; e_{3u} \ u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} } 153 153 \; {\overline {\left( {\frac{f}{e_{3f} }} \right) … … 220 220 -q\,e_3 \,u &\equiv -\frac{1}{e_{2v} } \left[ 221 221 {{\begin{array}{*{20}c} 222 {\,\ \ a_{j-1/2}^{i } \left( {e_{2u} e_{3 v} \ u} \right)_{j+1}^{i+1/2} }223 {\,+\,b_{j-1/2}^{i+1} \left( {e_{2u} e_{3 v} \ u} \right)_{j+1/2}^{i+1} } \\222 {\,\ \ a_{j-1/2}^{i } \left( {e_{2u} e_{3u} \ u} \right)_{j+1}^{i+1/2} } 223 {\,+\,b_{j-1/2}^{i+1} \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i+1} } \\ 224 224 \\ 225 { +\,c_{j+1/2}^{i+1} \left( {e_{2u} e_{3 v} \ u} \right)_{j+1/2}^{i+1} }226 {\,+\,d_{j+1/2}^{i } \left( {e_{2u} e_{3 v} \ u} \right)_{j+1/2}^{i } } \\225 { +\,c_{j+1/2}^{i+1} \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i+1} } 226 {\,+\,d_{j+1/2}^{i } \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i } } \\ 227 227 \end{array} }} \right] 228 228 \end{aligned} … … 670 670 \right]+\delta _j \left[ {e_{1v} e_{3v} v} \right]} \right)} 671 671 \end{equation} 672 673 where EMP is the surface freshwater budget, expressed in $Kg.m^{-2}.s^{-1}$,674 and $\rho _w =1,000\,Kg.m^{-3}$ is the volumic mass of pure water. If river675 runoff is expressed as a surface freshwater flux, see \S\ref{SBC}, then EMP676 can be written as the evaporation minus precipitation, minus the river runoff.677 The sea-surface height is evaluated using a leapfrog scheme in combination678 with an Asselin time filter, $i.e.$ the velocity appearing in679 \eqref{Eq_dynspg_ssh} is centred in time(\textit{now} velocity).672 where EMP is the surface freshwater budget, expressed in Kg/m$^2$/s 673 (which is equal to mm/s), and $\rho _w$=1,000~Kg/m$^3$ is the volumic 674 mass of pure water. If river runoff is expressed as a surface freshwater flux 675 (see \S\ref{SBC}) then EMP can be written as the evaporation minus 676 precipitation, minus the river runoff. The sea-surface height is evaluated 677 using a leapfrog scheme in combination with an Asselin time filter, $i.e.$ 678 the velocity appearing in \eqref{Eq_dynspg_ssh} is centred in time 679 (\textit{now} velocity). 680 680 681 681 The surface pressure gradient, also evaluated using a leap-frog scheme, is … … 683 683 \begin{equation} \label{Eq_dynspg_exp} 684 684 \left\{ \begin{aligned} 685 - \frac{1} {e_{1u}} \; \delta _{i+1/2} \left[ \,\eta\, \right] \\ 686 \\ 687 - \frac{1} {e_{2v}} \; \delta _{j+1/2} \left[ \,\eta\, \right] 685 - \frac{1}{e_{1u}} \; \delta _{i+1/2} \left[ \,\eta\, \right] \\ 686 - \frac{1}{e_{2v}} \; \delta _{j+1/2} \left[ \,\eta\, \right] 688 687 \end{aligned} \right. 689 688 \end{equation} … … 704 703 proposed by \citet{Griffies2004}. The general idea is to solve the free surface 705 704 equation with a small time step \np{rdtbt}, while the three dimensional 706 prognostic variables are solved with a longer time step that is a multiple of \np{rdtbt}707 (Fig.\ref {Fig_DYN_dynspg_ts}).705 prognostic variables are solved with a longer time step that is a multiple of 706 \np{rdtbt} (Fig.\ref {Fig_DYN_dynspg_ts}). 708 707 709 708 %> > > > > > > > > > > > > > > > > > > > > > > > > > > > … … 711 710 \begin{center} 712 711 \includegraphics[width=0.90\textwidth]{./Figures/Fig_DYN_dynspg_ts.pdf} 713 \caption{Schematic of the split-explicit time stepping scheme for the barotropic714 and baroclinic modes, after \citet{Griffies2004}. Time increases to the right.715 Baroclinic time steps are denoted by $t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$.716 The curved line represents a leap-frog time step, and the smaller barotropic time717 steps $N \Delta t=2\Delta t$ are denoted by the zig-zag line. The vertically718 integrated forcing \textbf{M}(t) computed at the baroclinic time step $t$719 represents the interaction between the barotropic and baroclinic motions.720 While keeping the total depth, tracer, and freshwater forcing fields fixed, a721 leap-frog integration carries the surface height and vertically integrated velocity722 from $t$ to $t+2 \Delta t$ using N barotropic time steps of length $\Delta t$.723 Time averaging the barotropic fields over the N+1time steps (endpoints724 included) centers the vertically integrated velocity a t the baroclinic timestep725 $t+\Delta t$. A baroclinic leap-frog time step carries the surface height to712 \caption{Schematic of the split-explicit time stepping scheme for the external 713 and internal modes. Time increases to the right. 714 Internal mode time steps (which are also the model time steps) are denoted 715 by $t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$. 716 The curved line represents a leap-frog time step, and the smaller time 717 steps $N \Delta t_e=\frac{3}{2}\Delta t$ are denoted by the zig-zag line. The vertically 718 integrated forcing \textbf{M}(t) computed at the model time step $t$ 719 represents the interaction between the external and internal motions. 720 While keeping \textbf{M} and freshwater forcing field fixed, a 721 leap-frog integration carries the external mode variables (surface height and vertically integrated velocity) from $t$ to $t+\frac{3}{2} \Delta t$ using N external time steps of length $\Delta t_e$. 722 Time averaging the external fields over the $\frac{2}{3}N+1$ time steps (endpoints 723 included) centers the vertically integrated velocity and the sea surface height at the model timestep $t+\Delta t$. These averaged values are used to update \textbf{M}(t) with both the surface pressure gradient and the Coriolis force. 724 A baroclinic leap-frog time step carries the surface height to The model time stepping scheme can then be achieved by 726 725 $t+\Delta t$ using the convergence of the time averaged vertically integrated 727 726 velocity taken from baroclinic time step t. } -
trunk/DOC/TexFiles/Chapters/Chap_LBC.tex
r817 r994 17 17 %-------------------------------------------------------------------------------------------------------------- 18 18 19 %The lateral ocean boundary condition continuous to coastlines are Neumann conditions for heat and salt (no flux across boundaries) and Dirichlet conditions for momentum (from free-slip to "strong" no-slip). They are handled automatically by the mask system (see \S\ref{DOM_msk}). 20 21 %OPA works with land and topography grid points in the computational domain due to the presence of continents or islands,and to the use of full or partial step representation of bottom topography. The computation is performed over the whole domain, i.e. we did not try to restrict the computation to ocean point only areas. This choice has two motivations. First, working on ocean only grid point overload the code and harms the code readability. Second,, and more importantly, it drastically reduce the vector portion of the computation, leading to a dramatic increase of CPU time requirement on vector computers. The process of defining which areas are to be masked is described in \S\ref{DOM_msk}, while this section describes how the masking affects the computation of the various terms of the equations with respect to the boundary condition at solid walls. 22 23 The discrete representation of a domain with complex boundaries (coastlines and bottom topography) leads to arrays that include large portions where a computation is not required as the model variables remain at zero. Nevertheless, vectorial supercomputers are far more efficient when computing over a whole array, and the readability of a code is greatly improved when boundary conditions are applied in an automatic way rather than by a specific computation before or after each do loop. An efficient way to work over the whole domain while specifying the boundary conditions is to use the multiplication by mask arrays in the computation. A mask array is a matrix 24 which elements are $1 $in the ocean domain and $0$ elsewhere. A simple multiplication of a variable by its own mask ensures that it will remain zero over land areas. Since most of the boundary conditions consist of a zero flux across the solid boundaries, they can be simply settled by 25 multiplying variables by the right mask arrays, i.e. the mask array of the grid point where the flux is evaluated. For example, the heat flux in the \textbf{i}-direction is evaluated at $u$-points. Evaluating this quantity as, 19 %The lateral ocean boundary conditions contiguous to coastlines are Neumann conditions for heat and salt (no flux across boundaries) and Dirichlet conditions for momentum (ranging from free-slip to "strong" no-slip). They are handled automatically by the mask system (see \S\ref{DOM_msk}). 20 21 %OPA allows land and topography grid points in the computational domain due to the presence of continents or islands, and includes the use of a full or partial step representation of bottom topography. The computation is performed over the whole domain, i.e. we do not try to restrict the computation to ocean-only points. This choice has two motivations. Firstly, working on ocean only grid points overloads the code and harms the code readability. Secondly, and more importantly, it drastically reduces the vector portion of the computation, leading to a dramatic increase of CPU time requirement on vector computers. The current section describes how the masking affects the computation of the various terms of the equations with respect to the boundary condition at solid walls. The process of defining which areas are to be masked is described in \S\ref{DOM_msk}. 22 23 The discrete representation of a domain with complex boundaries (coastlines and 24 bottom topography) leads to arrays that include large portions where a computation 25 is not required as the model variables remain at zero. Nevertheless, vectorial 26 supercomputers are far more efficient when computing over a whole array, and the 27 readability of a code is greatly improved when boundary conditions are applied in 28 an automatic way rather than by a specific computation before or after each 29 computational loop. An efficient way to work over the whole domain while specifying 30 the boundary conditions, is to use multiplication by mask arrays in the computation. 31 A mask array is a matrix whose elements are $1$ in the ocean domain and $0$ 32 elsewhere. A simple multiplication of a variable by its own mask ensures that it will 33 remain zero over land areas. Since most of the boundary conditions consist of a 34 zero flux across the solid boundaries, they can be simply applied by multiplying 35 variables by the correct mask arrays, $i.e.$ the mask array of the grid point where 36 the flux is evaluated. For example, the heat flux in the \textbf{i}-direction is evaluated 37 at $u$-points. Evaluating this quantity as, 26 38 27 39 \begin{equation} \label{Eq_lbc_aaaa} 28 40 \frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT} 29 }{e_{1u} } \delta _{i+1 / 2} \left[ T \right]\;\;mask_u41 }{e_{1u} } \; \delta _{i+1 / 2} \left[ T \right]\;\;mask_u 30 42 \end{equation} 31 where mask$_{u}$ is the mask array at $u$-point,ensures that the heat flux is32 zero inside land and at the boundaries as mask$_{u}$ is zero at solid33 boundaries defined at $u$-points in this case(normal velocity $u$ remains zero at43 (where mask$_{u}$ is the mask array at a $u$-point) ensures that the heat flux is 44 zero inside land and at the boundaries, since mask$_{u}$ is zero at solid boundaries 45 which in this case are defined at $u$-points (normal velocity $u$ remains zero at 34 46 the coast) (Fig.~\ref{Fig_LBC_uv}). 35 47 … … 37 49 \begin{figure}[!t] \label{Fig_LBC_uv} \begin{center} 38 50 \includegraphics[width=0.90\textwidth]{./Figures/Fig_LBC_uv.pdf} 39 \caption {Lateral boundary (thick line) at T-level. The velocity normal to the boundary is set to zero.} 51 \caption {Lateral boundary (thick line) at T-level. The velocity normal to the 52 boundary is set to zero.} 40 53 \end{center} \end{figure} 41 54 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 42 55 43 On momentum the situation is a bit more complex as two boundary conditions must be provided along the coast (one on the normal velocity and the other on the tangential velocity). The boundary of the ocean in C-grid is defined by the velocity-faces. For example, at a given $T$-level, the lateral boundary (coastline or intersection with the bottom topography) is made of segments 44 joining $f$-points, and normal velocity points are located between two $f-$points (Fig.~\ref{Fig_LBC_uv}). The boundary condition on the normal velocity (no flux through solid boundaries) can thus be easily settled by the mask system. The boundary condition on the tangential velocity requires a more specific treatment. It influences the relative vorticity and momentum diffusive trends, and is only required to compute the vorticity at the coast. Four different type of the lateral boundary condition are available, controlled by the value of \np{shlat} namelist parameter (which is equal to the value of the mask$_{f}$ array along the coastline): 56 For momentum the situation is a bit more complex as two boundary conditions 57 must be provided along the coast (one each for the normal and tangential velocities). 58 The boundary of the ocean in the C-grid is defined by the velocity-faces. 59 For example, at a given $T$-level, the lateral boundary (a coastline or an intersection 60 with the bottom topography) is made of segments joining $f$-points, and normal 61 velocity points are located between two $f-$points (Fig.~\ref{Fig_LBC_uv}). 62 The boundary condition on the normal velocity (no flux through solid boundaries) 63 can thus be easily implemented using the mask system. The boundary condition 64 on the tangential velocity requires a more specific treatment. This boundary 65 condition influences the relative vorticity and momentum diffusive trends, and is 66 required in order to compute the vorticity at the coast. Four different types of 67 lateral boundary condition are available, controlled by the value of the \np{shlat} 68 namelist parameter. (The value of the mask$_{f}$ array along the coastline is set 69 equal to this parameter.) These are: 45 70 46 71 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 53 78 \begin{description} 54 79 55 \item[free-slip boundary condition (\np{shlat}=0): ] the tangential velocity at the coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a). 56 57 \item[no-slip boundary condition (\np{shlat}=2): ] the tangential velocity vanishes at the coastline. Assuming that the tangential velocity decreases linearly from the closest ocean velocity grid point to the coastline, the normal derivative is evaluated as if the closest land velocity gridpoint were of the same magnitude as the closest ocean velocity gridpoint but in the opposite direction (Fig.~\ref{Fig_LBC_shlat}-b). Therefore, the vorticity along the coastlines is given by: 80 \item[free-slip boundary condition (\np{shlat}=0): ] the tangential velocity at the 81 coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the 82 tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set 83 to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a). 84 85 \item[no-slip boundary condition (\np{shlat}=2): ] the tangential velocity vanishes 86 at the coastline. Assuming that the tangential velocity decreases linearly from 87 the closest ocean velocity grid point to the coastline, the normal derivative is 88 evaluated as if the velocities at the closest land velocity gridpoint and the closest 89 ocean velocity gridpoint were of the same magnitude but in the opposite direction 90 (Fig.~\ref{Fig_LBC_shlat}-b). Therefore, the vorticity along the coastlines is given by: 91 58 92 \begin{equation*} 59 93 \zeta \equiv 2 \left(\delta_{i+1/2} \left[e_{2v} v \right] - \delta_{j+1/2} \left[e_{1u} u \right] \right) / \left(e_{1f} e_{2f} \right) \ , 60 94 \end{equation*} 61 where $u$ and $v$ are masked fields. Setting the mask$_{f}$ array to $2$ along the coastline allows to provide a vorticity field computed with the no-slip boundary condition simply by multiplying it by the mask$_{f}$ : 95 where $u$ and $v$ are masked fields. Setting the mask$_{f}$ array to $2$ along 96 the coastline provides a vorticity field computed with the no-slip boundary condition, 97 simply by multiplying it by the mask$_{f}$ : 62 98 \begin{equation} \label{Eq_lbc_bbbb} 63 99 \zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta _{i+1/2} … … 66 102 \end{equation} 67 103 68 \item["partial" free-slip boundary condition (0$<$\np{shlat}$<$2): ] the tangential velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral friction but not strong enough to vanish the tangential velocity at the coast (Fig.~\ref{Fig_LBC_shlat}-c). This can be settled by providing a value of mask$_{f}$ strictly inbetween $0$ and $2$. 69 70 \item["strong" no-slip boundary condition (2$<$\np{shlat}): ] the viscous boundary layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d). The friction is thus larger than in the no-slip case. 104 \item["partial" free-slip boundary condition (0$<$\np{shlat}$<$2): ] the tangential 105 velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral 106 friction but not strong enough to make the tangential velocity at the coast vanish 107 (Fig.~\ref{Fig_LBC_shlat}-c). This can be selected by providing a value of mask$_{f}$ 108 strictly inbetween $0$ and $2$. 109 110 \item["strong" no-slip boundary condition (2$<$\np{shlat}): ] the viscous boundary 111 layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d). 112 The friction is thus larger than in the no-slip case. 71 113 72 114 \end{description} 73 115 74 Note that when the bottom topography is entirely represented by the $s$-coordinates (pure $s$-coordinate), the lateral boundary condition on momentum tangential velocity is of much little importance as it is only applied next to the coast where the minimum water depth can be quite shallow. 75 76 The alternative numerical implementation of the no-slip boundary conditions for an arbitrary coast line of \citet{Shchepetkin1996} is also available through the activation of \key{noslip\_accurate}. It is based on a fourth order evaluation of the shear at the coast which, in turn, allows a true second order scheme in the interior of the domain ($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a technique considerably improves the quality of the numerical solution. In \NEMO, the improvement have not been found so spectacular in the half-degree global ocean (ORCA05), but significant reduction of numerically induced coast upwellings were found in eddy resolving simulation of the Alboran Sea \citep{OlivierPh2001}. Nevertheless, as no-slip boundary condition is not recommended in eddy permitting or resolving simulation \citep{Penduff2007}, the use of this option is not recommended. 77 78 In practice, the no-slip accurate option changes the way the curl is evaluated at the coast (see \mdl{divcur} module), and requires to qualify the nature of coastline grid point (convex or concave corners, straight north-south or east-west coast) which is performed in \mdl{domask} module, \rou{dom\_msk\_nsa} routine. 116 Note that when the bottom topography is entirely represented by the $s$-coor-dinates 117 (pure $s$-coordinate), the lateral boundary condition on tangential velocity is of much 118 less importance as it is only applied next to the coast where the minimum water depth 119 can be quite shallow. 120 121 The alternative numerical implementation of the no-slip boundary conditions for an 122 arbitrary coast line of \citet{Shchepetkin1996} is also available through the 123 \key{noslip\_accurate} CPP key. It is based on a fourth order evaluation of the shear at the 124 coast which, in turn, allows a true second order scheme in the interior of the domain 125 ($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical 126 scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a 127 technique considerably improves the quality of the numerical solution. In \NEMO, such 128 spectacular improvements have not been found in the half-degree global ocean 129 (ORCA05), but significant reductions of numerically induced coastal upwellings were 130 found in an eddy resolving simulation of the Alboran Sea \citep{OlivierPh2001}. 131 Nevertheless, since a no-slip boundary condition is not recommended in an eddy 132 permitting or resolving simulation \citep{Penduff2007}, the use of this option is also 133 not recommended. 134 135 In practice, the no-slip accurate option changes the way the curl is evaluated at the 136 coast (see \mdl{divcur} module), and requires the nature of each coastline grid point 137 (convex or concave corners, straight north-south or east-west coast) to be specified. 138 This is performed in routine \rou{dom\_msk\_nsa} in the \mdl{domask} module. 79 139 80 140 % ================================================================ … … 92 152 \label{LBC_jperio012} 93 153 94 The choice of closed, cyclic or symmetric model domain boundary condition is made by setting \jp{jperio} to 0, 1 or 2 in \mdl{par\_oce} file. Each time such a boundary condition is needed, it is set by a call to \mdl{lbclnk} routine. The computation of momentum and tracer trends proceed from $i=2$ to $i=jpi-1$ and from $j=2$ to $j=jpj-1$, $i.e.$in the inner model domain. To choose a lateral model boundary condition is to specify the first and last rows and columns of the model variables. 95 96 - For closed boundary (\textit{jperio=0}), solid walls are imposed at all model boundaries: 97 first and last rows and columns are set to zero. 98 99 - For cyclic east-west boundary (\textit{jperio=1}), first and last rows are set to zero 100 (closed) while first column is set to the value of the before last column 101 and last column to the value of the second one (Fig.~\ref{Fig_LBC_jperio}-a). Whatever flows 102 out of the eastern (western) end of the basin enters the western (eastern) 103 end. Note that there is neither option for north-south cyclic nor doubly 154 The choice of closed, cyclic or symmetric model domain boundary condition is made 155 by setting \jp{jperio} to 0, 1 or 2 in file \mdl{par\_oce}. Each time such a boundary 156 condition is needed, it is set by a call to routine \mdl{lbclnk}. The computation of 157 momentum and tracer trends proceeds from $i=2$ to $i=jpi-1$ and from $j=2$ to 158 $j=jpj-1$, $i.e.$ in the model interior. To choose a lateral model boundary condition 159 is to specify the first and last rows and columns of the model variables. 160 161 \begin{description} 162 163 \item[For closed boundary (\textit{jperio=0})], solid walls are imposed at all model 164 boundaries: first and last rows and columns are set to zero. 165 166 \item[For cyclic east-west boundary (\textit{jperio=1})], first and last rows are set 167 to zero (closed) whilst the first column is set to the value of the last-but-one column 168 and the last column to the value of the second one (Fig.~\ref{Fig_LBC_jperio}-a). 169 Whatever flows out of the eastern (western) end of the basin enters the western 170 (eastern) end. Note that there is no option for north-south cyclic or for doubly 104 171 cyclic cases. 105 172 106 - For symmetric boundary condition across the equator (\textit{jperio=2}), last rows, and 107 first and last columns are set to zero (closed). The row of symmetry is 108 chosen to be the $u$- and $T-$points equator line ($j=2$, i.e. at the southern end 109 of the domain). For arrays defined at $u-$ or $T-$points, the first row is set to 110 the value of the third row while for most of $v$- and $f$-points arrays (v, $\zeta$, 111 $j\psi$, but scalar arrays such as eddy coefficients) the first row is set to 112 minus the value of the second row (Fig.~\ref{Fig_LBC_jperio}-b). Note that this boundary 113 condition is not yet available on massively parallel computer 114 (\textbf{key{\_}mpp} defined). 173 \item[For symmetric boundary condition across the equator (\textit{jperio=2})], 174 last rows, and first and last columns are set to zero (closed). The row of symmetry 175 is chosen to be the $u$- and $T-$points equator line ($j=2$, i.e. at the southern 176 end of the domain). For arrays defined at $u-$ or $T-$points, the first row is set 177 to the value of the third row while for most of $v$- and $f$-point arrays ($v$, $\zeta$, 178 $j\psi$, but \gmcomment{not sure why this is "but"} scalar arrays such as eddy coefficients) 179 the first row is set to minus the value of the second row (Fig.~\ref{Fig_LBC_jperio}-b). 180 Note that this boundary condition is not yet available for the case of a massively 181 parallel computer (\textbf{key{\_}mpp} defined). 182 183 \end{description} 115 184 116 185 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 143 212 \label{LBC_mpp} 144 213 145 For massively parallel processing (mpp), a domain decomposition method is used. The basis of the method consists in splitting the large computation domain of a numerical experiment into several smaller domains and solving the set of equations by addressing independent local problems. Each processor has its own local memory and computes the model equation over a subdomain of the whole model domain. The subdomain boundary conditions are specified through communications between processors which are explicitly organized by specific statements (message passing method). 146 147 A big advantage is that the method does not need many modifications of the initial FORTRAN code. For the modeller's point of view, each sub domain running on a processor is identical to the "mono-domain" code. In addition, the programmer manages the communications between subdomains, and the code presents more scalability when the number of processors is increased. The porting of OPA code on a iPSC860 was achieved during Guyon's PhD [Guyon et al. 1994, 1995] in collaboration with CETIIS and ONERA. The implementation in the operational context and the studies of performances on a T3D and T3E Cray computers have been made in collaboration with IDRIS and CNRS. The present implementation is largely inspired from Guyon's work [Guyon 1995]. 148 149 The parallelization strategy is defined by the physical characteristics of the ocean model. Second order finite difference schemes leads to local discrete operators that depend at the very most on one neighbouring point. The only non-local computations concerne the vertical physics (implicit diffusion, 1.5 turbulent closure scheme, ...) (delocalization over the whole water column), and the solving of the elliptic equation associated with the surface pressure gradient computation (delocalization over the whole horizontal domain). Therefore, a pencil strategy is used for the data sub-structuration: the 3D initial domain is laid out on local processor memories following a 2D horizontal topological splitting. Each sub-domain computes its own surface and bottom boundary conditions and has a side wall overlapping interface which stocks lateral boundary conditions for computations in the inner sub-domain. The overlapping area is reduced to one row. After a computation, a communication phase starts: each processor sends to its neighbouring processors the update values of the point corresponding to the overlapping area of its neighbouring sub-domain. The communication is done through message passing. Usually the parallel virtual language, PVM, is used as it is a standard language available on nearly all MPP cumputers. More specific languages (i.e. computer dependant languages) can be easily used to speed up the communication, such as SHEM on T3E computer. The data exchanges between processors are required at the very place where lateral domain boundary conditions are set in the mono-domain computation (\S III.10-c): the lbc\_lnk routine which manages such conditions is substituted by mpplnk.F or mpplnk2.F routine when running on MPP computer (\key{mpp\_mpi} defined). It has to be noticed that when using MPP version of the model, the east-west cyclic boundary condition is implicitly done, while the south-symmetric boundary condition option is not available. 150 214 For massively parallel processing (mpp), a domain decomposition method is used. 215 The basic idea of the method is to split the large computation domain of a numerical 216 experiment into several smaller domains and solve the set of equations by addressing 217 independent local problems. Each processor has its own local memory and computes 218 the model equation over a subdomain of the whole model domain. The subdomain 219 boundary conditions are specified through communications between processors 220 which are organized by explicit statements (message passing method). 221 222 A big advantage is that the method does not need many modifications of the initial 223 FORTRAN code. From the modeller's point of view, each sub domain running on 224 a processor is identical to the "mono-domain" code. In addition, the programmer 225 manages the communications between subdomains, and the code is faster when 226 the number of processors is increased. The porting of OPA code on an iPSC860 227 was achieved during Guyon's PhD [Guyon et al. 1994, 1995] in collaboration with 228 CETIIS and ONERA. The implementation in the operational context and the studies 229 of performance on a T3D and T3E Cray computers have been made in collaboration 230 with IDRIS and CNRS. The present implementation is largely inspired by Guyon's 231 work [Guyon 1995]. 232 233 The parallelization strategy is defined by the physical characteristics of the 234 ocean model. Second order finite difference schemes lead to local discrete 235 operators that depend at the very most on one neighbouring point. The only 236 non-local computations concern the vertical physics (implicit diffusion, 1.5 237 turbulent closure scheme, ...) (delocalization over the whole water column), 238 and the solving of the elliptic equation associated with the surface pressure 239 gradient computation (delocalization over the whole horizontal domain). 240 Therefore, a pencil strategy is used for the data sub-structuration \gmcomment{no 241 idea what this means!}: the 3D initial domain is laid out on local processor 242 memories following a 2D horizontal topological splitting. Each sub-domain 243 computes its own surface and bottom boundary conditions and has a side 244 wall overlapping interface which defines the lateral boundary conditions for 245 computations in the inner sub-domain. The overlapping area consists of the 246 two rows at each edge of the sub-domain. After a computation, a communication 247 phase starts: each processor sends to its neighbouring processors the update 248 values of the points corresponding to the interior overlapping area to its 249 neighbouring sub-domain (i.e. the innermost of the two overlapping rows). The communication is done through message passing. Usually the parallel virtual 250 language, PVM, is used as it is a standard language available on nearly all 251 MPP computers. More specific languages (i.e. computer dependant languages) 252 can be easily used to speed up the communication, such as SHEM on a T3E 253 computer. The data exchanges between processors are required at the very 254 place where lateral domain boundary conditions are set in the mono-domain 255 computation (\S III.10-c): the lbc\_lnk routine which manages such conditions 256 is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP 257 computer (\key{mpp\_mpi} defined). It has to be pointed out that when using 258 the MPP version of the model, the east-west cyclic boundary condition is done 259 implicitly, whilst the south-symmetric boundary condition option is not available. 260 151 261 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 152 262 \begin{figure}[!t] \label{Fig_mpp} \begin{center} … … 156 266 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 157 267 158 In the standard version of the OPA model, the splitting is regular and arithmetic. the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in \mdl{par\_oce}). Each processor is independent and without message passing or synchronous process, programs run alone and access just at its own local memory. For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil) that are noted \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal domain and the overlapping rows. The number of overlapping rows is usually set to one (\jp{jpreci}=1, in \mdl{par\_oce}). The whole domain dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. The relationship between the whole domain and a sub-domain is: 268 In the standard version of the OPA model, the splitting is regular and arithmetic. 269 the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors 270 \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in 271 \mdl{par\_oce}). Each processor is independent and without message passing 272 or synchronous process \gmcomment{how does a synchronous process relate to this?}, 273 programs run alone and access just its own local memory. For this reason, the 274 main model dimensions are now the local dimensions of the subdomain (pencil) 275 that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal 276 domain and the overlapping rows. The number of rows to exchange (known as 277 the halo) is usually set to one (\jp{jpreci}=1, in \mdl{par\_oce}). The whole domain 278 dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. The relationship between 279 the whole domain and a sub-domain is: 159 280 \begin{eqnarray} 160 281 jpi & = & ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci \nonumber \\ … … 165 286 \colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and no east-west cyclic boundary conditions.} 166 287 167 One defines also variables nldi and nlei which correspond to the internal domain bounds, and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$, a global array (whole domain) by the relationship: 288 One also defines variables nldi and nlei which correspond to the internal 289 domain bounds, and the variables nimpp and njmpp which are the position 290 of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array 291 (subdomain) corresponds to an element of $T_{g}$, a global array 292 (whole domain) by the relationship: 168 293 \begin{equation} \label{Eq_lbc_nimpp} 169 294 T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k), … … 171 296 with $1 \leq i \leq jpi$, $1 \leq j \leq jpj $ , and $1 \leq k \leq jpk$. 172 297 173 Processors are numbered from 0 to $jpnij-1$, the number is saved in the variable nproc. In the standard version, a processor has no more than four neighbouring processors named nono (for north), noea (east), noso (south) and nowe (west) and two variables, nbondi and nbondj, indicate the situation of the processor \colorbox{yellow}{(see Fig.IV.3)}: 298 Processors are numbered from 0 to $jpnij-1$, the number is saved in the variable 299 nproc. In the standard version, a processor has no more than four neighbouring 300 processors named nono (for north), noea (east), noso (south) and nowe (west) 301 and two variables, nbondi and nbondj, indicate the relative position of the processor 302 \colorbox{yellow}{(see Fig.IV.3)}: 174 303 \begin{itemize} 175 304 \item nbondi = -1 an east neighbour, no west processor, … … 185 314 186 315 187 The OPA model computes equation terms with the help of mask arrays (0 onto land points and 1 onto sea points). It is easily readable and very efficient in the context of the vectorial architecture. Nevertheless, in the case of scalar processor, computations over the land regions becomes more expensive in term of CPU time. It is all the more so when we use a complex configuration with a realistic bathymetry like the global ocean where more than 50 \% of points are land points. For this reason, a pre-processing tool allows to search in the mpp domain decomposition strategy if a splitting can be found with a maximum number of only land points processors which could be eliminated: the mpp\_optimiz tools (available from the DRAKKAR web site). This optimisation is made with the knowledge of the specific bathymetry. The user chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with $jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ land processors. When those parameters are specified in module \mdl{par\_oce}, the algorithm in the \rou{inimpp2} routine set each processor name and neighbour parameters (nbound, nono, noea,...) so that the land processors are not taken into account. 316 The OPA model computes equation terms with the help of mask arrays (0 on land 317 points and 1 on sea points). It is easily readable and very efficient in the context of 318 a computer with vectorial architecture. However, in the case of a scalar processor, 319 computations over the land regions become more expensive in terms of CPU time. 320 It is worse when we use a complex configuration with a realistic bathymetry like the 321 global ocean where more than 50 \% of points are land points. For this reason, a 322 pre-processing tool can be used to choose the mpp domain decomposition with a 323 maximum number of only land points processors, which can then be eliminated. 324 (For example, the mpp\_optimiz tools, available from the DRAKKAR web site.) 325 This optimisation is dependent on the specific bathymetry employed. The user 326 then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with 327 $jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ 328 land processors. When those parameters are specified in module \mdl{par\_oce}, 329 the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound, 330 nono, noea,...) so that the land-only processors are not taken into account. 188 331 189 332 \colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp routine should be suppressed from the code.} 190 333 191 334 When land processors are eliminated, the value corresponding to these locations in the model output files is zero. Note that this is a problem for a mesh output file written by such a model configuration, because model users often divide by the scale factors ($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be best not to eliminate land processors when running the model especially to write the mesh files as outputs (when \np{nmsh} namelist parameter differs from 0). 335 \gmcomment{Steven : dont understand this, no land processor means no output file covering this part of globe; its only when files are stitched together into one that you can leave a hole} 192 336 193 337 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 194 338 \begin{figure}[!ht] \label{Fig_mppini2} \begin{center} 195 339 \includegraphics[width=0.90\textwidth]{./Figures/Fig_mppini2.pdf} 196 \caption {Example of Atlantic domain defined for the CLIPPER projet. Initial grid is composed by 773 x 1236 horizontal points. (a) the domain is splitting onto 9 \time 20 subdomains (jpni=9, jpnj=20). 52 subdomains are land areas. (b) 52 subdomains are eliminated (white rectangles) and the resulting number of processors really used during the computation is jpnij=128.} 340 \caption {Example of Atlantic domain defined for the CLIPPER projet. Initial grid is 341 composed of 773 x 1236 horizontal points. (a) the domain is split onto 9 \time 20 342 subdomains (jpni=9, jpnj=20). 52 subdomains are land areas. (b) 52 subdomains 343 are eliminated (white rectangles) and the resulting number of processors really 344 used during the computation is jpnij=128.} 197 345 \end{center} \end{figure} 198 346 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 222 370 \namdisplay{namobc} 223 371 224 It is often necessary to implement a model configuration limited to an oceanic region or a basin, which communicates with the rest of the global ocean through ``open boundaries''. As stated by \citet{Roed1986}, an open boundary is a computational border where the aim of the calculations is to allow the perturbations generated inside the computational domain to leave it without deterioration of the inner model solution. However, an open boundary also has to let information from the outer oceans enter the model and should support inflow and outflow conditions. 225 226 The open boundary package OBC is the first open boundary option developed in NEMO (originally in OPA8.2). It allows the user to 372 It is often necessary to implement a model configuration limited to an oceanic 373 region or a basin, which communicates with the rest of the global ocean through 374 ''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a 375 computational border where the aim of the calculations is to allow the perturbations 376 generated inside the computational domain to leave it without deterioration of the 377 inner model solution. However, an open boundary also has to let information from 378 the outer ocean enter the model and should support inflow and outflow conditions. 379 380 The open boundary package OBC is the first open boundary option developed in 381 NEMO (originally in OPA8.2). It allows the user to 227 382 \begin{itemize} 228 \item tell the model that a boundary is ``open'' and not closed by a wall, for example by modifying the calculation of the divergence of velocity there, and other necessary modifications;383 \item tell the model that a boundary is ''open'' and not closed by a wall, for example by modifying the calculation of the divergence of velocity there; 229 384 \item impose values of tracers and velocities at that boundary (values which may be taken from a climatology): this is the``fixed OBC'' option. 230 385 \item calculate boundary values by a sophisticated algorithm combining radiation and relaxation (``radiative OBC'' option) 231 386 \end{itemize} 232 387 233 The package resides in the OBC directory. It is described here in four parts: the boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at the boundaries (module \mdl{obcdta}), the radiation algorithm involving the namelist and module \mdl{obcrad}, and a brief presentation of boundary update and restart files. 388 The package resides in the OBC directory. It is described here in four parts: the 389 boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at 390 the boundaries (module \mdl{obcdta}), the radiation algorithm involving the 391 namelist and module \mdl{obcrad}, and a brief presentation of boundary update 392 and restart files. 234 393 235 394 %---------------------------------------------- … … 237 396 \label{OBC_geom} 238 397 % 239 First one has to realize that open boundaries may not necessarily be located at the extremities of the computational domain. They may exist in the middle of the domain, for example at Gibraltar Straits if one wants to avoid including the Mediterranean in an Atlantic domain. This flexibility has been found necessary for the CLIPPER project \citep{Treguier2001}. Because of the complexity of the geometry of ocean basins, it may even be necessary to have more than one ``west'' open boundary, more than one ``north'', etc. This is not possible with the OBC option: only one open boundary of each kind, west, east, south and north is allowed; those names refer to the grid geometry (not to the direction of the geographical ``west'', ``east'', etc). 240 241 The open boundary geometry is set by a series of parameters in the module \mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east} (true if an east open boundary exists), \jp{jpieob} the $i$-index along which the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ``d\'{e}but'' and $f$ for ``fin'' in French). Similar parameters exist for the west, south and north cases (Table~\ref{Tab_obc_param}). 398 First one has to realize that open boundaries may not necessarily be located 399 at the extremities of the computational domain. They may exist in the middle 400 of the domain, for example at Gibraltar Straits if one wants to avoid including 401 the Mediterranean in an Atlantic domain. This flexibility has been found necessary 402 for the CLIPPER project \citep{Treguier2001}. Because of the complexity of the 403 geometry of ocean basins, it may even be necessary to have more than one 404 ''west'' open boundary, more than one ''north'', etc. This is not possible with 405 the OBC option: only one open boundary of each kind, west, east, south and 406 north is allowed; these names refer to the grid geometry (not to the direction 407 of the geographical ''west'', ''east'', etc). 408 409 The open boundary geometry is set by a series of parameters in the module 410 \mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east} 411 (true if an east open boundary exists), \jp{jpieob} the $i$-index along which 412 the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which 413 it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ''d\'{e}but'' 414 and $f$ for ''fin'' in French). Similar parameters exist for the west, south and 415 north cases (Table~\ref{Tab_obc_param}). 242 416 243 417 … … 265 439 \end{tabular} 266 440 \end{center} 267 \caption{Names of different indices concerning the open boundaries. In the case of a completely open ocean domain with four ocean boundaries, the parameters take exactly the values indicated.} 441 \caption{Names of different indices relating to the open boundaries. In the case 442 of a completely open ocean domain with four ocean boundaries, the parameters 443 take exactly the values indicated.} 268 444 \end{table} 269 445 270 The open boundaries must be along coordinate lines. On the C-grid, the boundary itself is along a line of normal velocity points: $v$ points for a zonal open boundary (the south or north one), and $u$ points for a meridional open boundary (the west or east one). Another constraint is that there still must be a row of masked points all around the domain, as if the domain were a closed basin (unless periodic conditions are used together with open boundary conditions). Therefore, an open boundary cannot be located at the first/last index, namely, 1 or \jp{jpiglo}, \jp{jpjglo}. Also, the open boundary algorithm involves calculating the normal velocity points situated just on the boundary, as well as the tangential velocity and temperature, salinity just outside the boundary. This means that for a west/south boundary, normal velocities and temperature are calculated at the same index \jp{jpiwob} and \jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob} cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2. 271 272 273 The starting and ending indices are to be thought as $T$ point indices: in many cases they indicate the first land $T$-point, at the extremity of an open boundary (the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for the example for a northern open boundary. All indices are relative to the global domain. In the free surface case it is possible to have ``ocean corners'', that is, an open boundary starting and ending in the ocean. 446 The open boundaries must be along coordinate lines. On the C-grid, the boundary 447 itself is along a line of normal velocity points: $v$ points for a zonal open boundary 448 (the south or north one), and $u$ points for a meridional open boundary (the west 449 or east one). Another constraint is that there still must be a row of masked points 450 all around the domain, as if the domain were a closed basin (unless periodic conditions 451 are used together with open boundary conditions). Therefore, an open boundary 452 cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also, 453 the open boundary algorithm involves calculating the normal velocity points situated 454 just on the boundary, as well as the tangential velocity and temperature and salinity 455 just outside the boundary. This means that for a west/south boundary, normal 456 velocities and temperature are calculated at the same index \jp{jpiwob} and 457 \jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is 458 calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is 459 at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob} 460 cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2. 461 462 463 The starting and ending indices are to be thought of as $T$ point indices: in many 464 cases they indicate the first land $T$-point, at the extremity of an open boundary 465 (the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example 466 of a northern open boundary). All indices are relative to the global domain. In the 467 free surface case it is possible to have ``ocean corners'', that is, an open boundary 468 starting and ending in the ocean. 274 469 275 470 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 281 476 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 282 477 283 Although not absolutely compulsory, it is highly recommended that the bathymetry in the vicinity of an open boundary follows the following rule: in the direction perpendicular to the open line, the water depth should be constant for 4 grid points. This is in order to ensure that the radiation condition, which involves mdoel variables next to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we indicate by a $=$ symbol, the points which should have the same depth. It means that in the 4 points near the boundary, the bathymetry is cylindrical. The line behind the open T-line must be 0 in the bathymetry file (as shown on Fig.\ref{Fig_obc_north} for example). 478 Although not compulsory, it is highly recommended that the bathymetry in the 479 vicinity of an open boundary follows the following rule: in the direction perpendicular 480 to the open line, the water depth should be constant for 4 grid points. This is in 481 order to ensure that the radiation condition, which involves model variables next 482 to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we 483 indicate by an $=$ symbol, the points which should have the same depth. It means 484 that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure 485 why cylindrical}. The line behind the open $T$-line must be 0 in the bathymetry file 486 (as shown on Fig.\ref{Fig_obc_north} for example). 284 487 285 488 %---------------------------------------------- … … 287 490 \label{OBC_data} 288 491 289 It is necessary to provide information at the boundaries. The simple case happens when this information does not change in time and is equal to the initial conditions (namelist variable \np{nobc\_dta}=0). This is the case of the standard configuration EEL5 with open boundaries. When (\np{nobc\_dta}=1), open boundary information is read from netcdf files. For convenience the input files are supposed to be similar to the ''history'' NEMO output files, for dimension names and variable names. 290 Open boundary arrays must be dimensioned according to the parameters of table~\ref{Tab_obc_param}: for example, at the western boundary arrays have a dimension of \jp{jpwf}-\jp{jpwd}+1 on the horizontal and \jp{jpk} on the vertical. 291 292 When ocean observations are used to generate the boundary data (a hydrographic section for example, as in \citet{Treguier2001}) it happens often that only the velocity normal to the boundary is known, which is the reason why the initial OBC code assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be specified. As more and more global model solutions and ocean analysis products are available, it is possible to provide information about all the variables (including the tangential velocity) so that the specification of four variables at each boundaries will soon become standard. Regarding the sea surface height, one must distinguish the filtered free surface case from the time-splitting or explicit treatment of the free surface. 293 In the first case, it is assumed that the user does not wish to represent high frequency motions such as tides. The boundary condition is thus one of zero normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}. 294 No information other than the total velocity needs to be provided at the open boundaries in that case. In the other two cases (time stplitting or explicit free surface), the user must provides barotropic information (sea surface height and barotropic velocities) and the use of the Flather algorithm for barotropic variables is recommanded. However, this algorithm has not yet been fully tested and bugs remain in NEMO v2.3. Users should read the code carefully before using it. Finally, in the case of the rigid lid approximation the barotropic streamfunction must be provided, as documented in \citet{Treguier2001}). This option is no longer used but remains in NEMO V2.3. 295 296 One frequently encountered case is when an open boundary domain is constructed from a global or larger scale NEMO configuration. Assuming the domain corresponds to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the small domain can be created by using the following command on the global files: ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$. 297 The open boundary files can be constructed using ncks commands, following table~\ref{Tab_obc_ind}. 492 It is necessary to provide information at the boundaries. The simplest case is 493 when this information does not change in time and is equal to the initial conditions 494 (namelist variable \np{nobc\_dta}=0). This is the case for the standard configuration 495 EEL5 with open boundaries. When (\np{nobc\_dta}=1), open boundary information 496 is read from netcdf files. For convenience the input files are supposed to be similar 497 to the ''history'' NEMO output files, for dimension names and variable names. 498 Open boundary arrays must be dimensioned according to the parameters of table~ 499 \ref{Tab_obc_param}: for example, at the western boundary, arrays have a 500 dimension of \jp{jpwf}-\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical. 501 502 When ocean observations are used to generate the boundary data (a hydrographic 503 section for example, as in \citet{Treguier2001}) it happens often that only the velocity 504 normal to the boundary is known, which is the reason why the initial OBC code 505 assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be 506 specified. As more and more global model solutions and ocean analysis products 507 become available, it will be possible to provide information about all the variables 508 (including the tangential velocity) so that the specification of four variables at each 509 boundaries will become standard. For the sea surface height, one must distinguish 510 between the filtered free surface case and the time-splitting or explicit treatment of 511 the free surface. 512 In the first case, it is assumed that the user does not wish to represent high 513 frequency motions such as tides. The boundary condition is thus one of zero 514 normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}. 515 No information other than the total velocity needs to be provided at the open 516 boundaries in that case. In the other two cases (time splitting or explicit free surface), 517 the user must provide barotropic information (sea surface height and barotropic 518 velocities) and the use of the Flather algorithm for barotropic variables is 519 recommanded. However, this algorithm has not yet been fully tested and bugs 520 remain in NEMO v2.3. Users should read the code carefully before using it. Finally, 521 in the case of the rigid lid approximation the barotropic streamfunction must be 522 provided, as documented in \citet{Treguier2001}). This option is no longer 523 recommended but remains in NEMO V2.3. 524 525 One frequently encountered case is when an open boundary domain is constructed 526 from a global or larger scale NEMO configuration. Assuming the domain corresponds 527 to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the 528 small domain can be created by using the following netcdf utility on the global files: 529 ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities, see http://nco.sourceforge.net). The open boundary files can be constructed using ncks 530 commands, following table~\ref{Tab_obc_ind}. 298 531 299 532 %--------------------------------------------------TABLE-------------------------------------------------- … … 322 555 \end{tabular} 323 556 \end{center} 324 \caption{Indications for creating open boundary files from a global configuration, appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the $i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global configuration, starting and ending with the $j$ or $i$ indices indicated. For example, to generate file obcnorth\_V.nc, use the command ncks $-F$ $-d\;y,je-2$ $-d\;x,ib+1,ie-1$ } 557 \caption{Requirements for creating open boundary files from a global configuration, 558 appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the 559 $i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global 560 configuration, starting and ending with the $j$ or $i$ indices indicated. 561 For example, to generate file obcnorth\_V.nc, use the command ncks 562 $-F$ $-d\;y,je-2$ $-d\;x,ib+1,ie-1$ } 325 563 \end{table} 326 564 327 It is assumed that the open boundary files contain the variables for the period of the model integration. If the boundary files contain one time frame, the boundary data is held fixed in time. If the files contain 12 values, it is assumed that the input is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim} = .T.). The case of an arbitrary number of time frames is not yet implemented correctly; the user is supposed to write his own code in the module \mdl{obc\_dta} to deal with this situation. 565 It is assumed that the open boundary files contain the variables for the period of 566 the model integration. If the boundary files contain one time frame, the boundary 567 data is held fixed in time. If the files contain 12 values, it is assumed that the input 568 is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim} 569 = .True.). The case of an arbitrary number of time frames is not yet implemented 570 correctly; the user is required to write his own code in the module \mdl{obc\_dta} 571 to deal with this situation. 328 572 329 573 \subsection{Radiation algorithm} 330 574 \label{OBC_rad} 331 575 332 The art of open boundary management consists in applying a constraint strong enough so that the inner domain "feels" the rest of the ocean, but not too strong 333 in order to allow perturbations to leave the domain with minimum false reflexions of energy. The cosntraint to the specified boundary data is set by time scales 334 specified separately for each boundary and for ``inflow'' and `outflow'' conditions as defined below. These time scales are set (in days) by namelist parameters such as \np{rdpein}, \np{rdpeob} for the eastern open boundary, for example. When both time scales are zero for a given boundary, for example 335 \jp{lp\_obc\_west}=.T., \np{rdpwob}=0, \np{rdpwin}=0, this means that the boundary in question (western boundary here) is a ``fixed '' boundary where the solution is set exactly by the boundary data. This is not recommanded, excepted in compination with increased viscosity in a ``sponge'' layer next to the boundary in order to avois spurious reflexions. 336 337 338 The radiation\/relaxation algorithm is applied when either relxation time (for ``inflow'' or `outflow'') is nonzero. It has been developed and tested in the SPEM model and its successor ROMS \citep{Barnier1996, Marchesiello2001}, a $s$-coordinate model on an Arakawa C-grid. Although the algorithm has been numerically successful in the CLIPPER Atlantic models, the physics do not work as expected \citep{Treguier2001}. Users are invited to consider open boundary conditions (OBC hereafter) with some skepticism \citep{Durran2001, Blayo2005}. 339 340 The first part of the algorithm consists in calculating a phase 341 velocity to determine whether perturbations tend to propagate toward, 342 or away from, the boundary. 343 Let us consider a model variable $\phi$. 344 The phase velocity ($C_{\phi x}$,$C_{\phi y}$) for the variable 345 $\phi$, in the directions normal and tangential to 346 the boundary is 576 The art of open boundary management consists in applying a constraint strong 577 enough that the inner domain "feels" the rest of the ocean, but weak enough 578 that perturbations are allowed to leave the domain with minimum false reflections 579 of energy. The constraints are specified separately at each boundary as time 580 scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days) 581 by namelist parameters such as \np{rdpein}, \np{rdpeob} for the eastern open 582 boundary for example. When both time scales are zero for a given boundary 583 ($e.g.$ for the western boundary, \jp{lp\_obc\_west}=.True., \np{rdpwob}=0 and 584 \np{rdpwin}=0) this means that the boundary in question is a ''fixed '' boundary 585 where the solution is set exactly by the boundary data. This is not recommended, 586 except in combination with increased viscosity in a ''sponge'' layer next to the 587 boundary in order to avoid spurious reflections. 588 589 590 The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output} 591 algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is 592 non-zero. It has been developed and tested in the SPEM model and its 593 successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an 594 $s$-coordinate model on an Arakawa C-grid. Although the algorithm has 595 been numerically successful in the CLIPPER Atlantic models, the physics 596 do not work as expected \citep{Treguier2001}. Users are invited to consider 597 open boundary conditions (OBC hereafter) with some scepticism 598 \citep{Durran2001, Blayo2005}. 599 600 The first part of the algorithm calculates a phase velocity to determine 601 whether perturbations tend to propagate toward, or away from, the 602 boundary. Let us consider a model variable $\phi$. 603 The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$, 604 in the directions normal and tangential to the boundary are 347 605 \begin{equation} \label{Eq_obc_cphi} 348 606 C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x} … … 350 608 C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}. 351 609 \end{equation} 352 Following \citet{Treguier2001} and \citet{Marchesiello2001} retain only the normal353 projection of the total velocity, $C_{\phi x}$, setting 354 $C_{\phi y} =0$ (but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in the expression 355 for $C_{\phi x}$).610 Following \citet{Treguier2001} and \citet{Marchesiello2001} we retain only 611 the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$ 612 (but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in 613 the expression for $C_{\phi x}$). 356 614 357 615 The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998}, 358 616 takes into account the two rows of grid points situated inside the domain 359 617 next to the boundary, and the three previous time steps ($n$, $n-1$, 360 and $n-2$). 361 The same equation can then be discretized at the boundary at 362 time steps $n$ and $n+1$ in order to extrapolate the new boundary 363 value $\phi^{n+1}$. 364 365 In the present open boundary algorithm, the new boundary values are366 updated differently according to the sign of $C_{\phi x}$. 367 Let us take an East boundary as an example. The solution for variable $\phi$ at the boundary is given from a generalized wave equation368 with phase velocity $C_{\phi}$, with the addition of a relaxation term:618 and $n-2$). The same equation can then be discretized at the boundary at 619 time steps $n-1$, $n$ and $n+1$ \gmcomment{since the original was three time-level} 620 in order to extrapolate for the new boundary value $\phi^{n+1}$. 621 622 In the open boundary algorithm as implemented in NEMO v2.3, the new boundary 623 values are updated differently depending on the sign of $C_{\phi x}$. Let us take 624 an eastern boundary as an example. The solution for variable $\phi$ at the 625 boundary is given by a generalized wave equation with phase velocity $C_{\phi}$, 626 with the addition of a relaxation term, as: 369 627 \begin{eqnarray} 370 628 \phi_{t} & = & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi) … … 373 631 \;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi} 374 632 \end{eqnarray} 375 where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary data. 376 Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ 377 is bounded by the ratio $\delta x/\delta t$ for stability reasons. 378 When $C_{\phi x}$ is eastward (outward propagation), 379 the radiation condition (\ref{Eq_obc_rado}) is used. 380 When $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is used with 381 a strong relaxation to climatology (usually $\tau_{i}=\np{rdpein}=$1~day). 382 The time derivative in (\ref{Eq_obc_radi}) is calculated with a Euler time-stepping 383 scheme. It results from this choice that setting 384 $\tau_{i}$ smaller than, or equal to the time step is in fact equivalent to a fixed boundary condition; 385 a time scale of one day is usually a good compromise which guarantees that the inflow conditions remain close to 386 climatology while ensuring numerical stability. 387 388 In the case of a west boundary located in the Eastern Atlantic, \citet{Penduff2000} have been able to implement the radiation algorithm 389 without any boundary data, using persistence from the previous time step instead. This solution did not work in other cases \citep{Treguier2001} so that the use of boundary data is recommended. 390 Even in the outflow condition (\ref{Eq_obc_rado}), we have found desireable to maintain a weak relaxation to climatology. 391 The time step is usually chosen so as to be larger than typical turbulent scales (of order 1000~days). 392 393 The radiation condition is applied to the different model variables: temperature, salinity, tangential and normal velocities. 394 For normal and tangential velocities $u$ and $v$ radiation is applied with phase velocities calculated from $u$ and $v$ respectively. 395 For the radiation of tracers, we use the phase velocity calculated from the tangential velocity, to avoid calculating too many independent 396 radiation velocities and because tangential velocities and tracers have the same position along the boundary on a C grid. 633 where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary 634 data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio 635 $\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward 636 propagation), the radiation condition (\ref{Eq_obc_rado}) is used. 637 When $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is 638 used with a strong relaxation to climatology (usually $\tau_{i}=\np{rdpein}=$1~day). 639 Equation (\ref{Eq_obc_radi}) is solved with a Euler time-stepping scheme. As a 640 consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent 641 to a fixed boundary condition. A time scale of one day is usually a good compromise 642 which guarantees that the inflow conditions remain close to climatology while ensuring 643 numerical stability. 644 645 In the case of a western boundary located in the Eastern Atlantic, \citet{Penduff2000} 646 have been able to implement the radiation algorithm without any boundary data, 647 using persistence from the previous time step instead. This solution has not worked 648 in other cases \citep{Treguier2001}, so that the use of boundary data is recommended. 649 Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to 650 maintain a weak relaxation to climatology. The time step is usually chosen so as to 651 be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}). 652 653 The radiation condition is applied to the model variables: temperature, salinity, 654 tangential and normal velocities. For normal and tangential velocities, $u$ and $v$, 655 radiation is applied with phase velocities calculated from $u$ and $v$ respectively. 656 For the radiation of tracers, we use the phase velocity calculated from the tangential 657 velocity in order to avoid calculating too many independent radiation velocities and 658 because tangential velocities and tracers have the same position along the boundary 659 on a C-grid. 397 660 398 661 \subsection{Domain decomposition (\key{mpp\_mpi})} 399 662 \label{OBC_mpp} 400 When \key{mpp\_mpi} is active in the code, the computational domain is divided into rectangles that are attributed each to a different processor. The open boundary code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not work if there is a mpp subdomain boundary parallel to the open boundary at the index of the boundary, or the grid point after (outside), or three grid points before (inside). On the other hand, there is no problem if a mpp subdomain boundary cuts the open boundary perpendicularly. Those geometry constraints must be checked by the user (there is no safeguard in the code). 401 The general principle for the open boundary mpp code is that loops over the open boundaries are performed on local indices (nie0, nie1, nje0, nje1 for the eastern boundary for instance) that are initialized in module \mdl{obc\_ini}. Those indices have relevant values on the processors that contain a segment of an open boundary. For processors that do not include an open boundary segment, the indices are such that the calculations within the loops are not performed. 663 When \key{mpp\_mpi} is active in the code, the computational domain is divided 664 into rectangles that are attributed each to a different processor. The open boundary 665 code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not 666 work if there is an mpp subdomain boundary parallel to the open boundary at the 667 index of the boundary, or the grid point after (outside), or three grid points before 668 (inside). On the other hand, there is no problem if an mpp subdomain boundary 669 cuts the open boundary perpendicularly. These geometrical limitations must be 670 checked for by the user (there is no safeguard in the code). 671 The general principle for the open boundary mpp code is that loops over the open 672 boundaries {not sure what this means} are performed on local indices (nie0, 673 nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module 674 \mdl{obc\_ini}. Those indices have relevant values on the processors that contain 675 a segment of an open boundary. For processors that do not include an open 676 boundary segment, the indices are such that the calculations within the loops are 677 not performed. 678 \gmcomment{I dont understand most of the last few sentences} 402 679 403 Arrays of climatological data that are read in files are seen by all processors and have the same dimensions for all (for instance, for the eastern boundary, uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation are local to each processor (uebnd(jpj,jpk,3,3) for instance). This allowed to spare memory in the CLIPPER model, where the eastern boundary crossed 8 processors so that \jp{jpj} was must smaller than (\jp{jpjef}-\jp{jpjed}+1). 680 Arrays of climatological data that are read from files are seen by all processors 681 and have the same dimensions for all (for instance, for the eastern boundary, 682 uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation 683 are local to each processor (uebnd(jpj,jpk,3,3) for instance). This allowed the 684 CLIPPER model for example, to save on memory where the eastern boundary 685 crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}-\jp{jpjed}+1). 404 686 405 687 \subsection{Volume conservation} 406 688 \label{OBC_vol} 407 689 408 It is necessary to control the volume inside a domain with open boundaries. With fixed boundaries, it is enough to ensure that the total inflow/outflow has reasonable values (either zero or a value compatible with an observed volume balance). When using radiative boundary conditions it is necessary to have a volume constraint because each open boundary works independently from the others. The methodology used to control this volume is identical to the one coded in the ROMS model \citep{Marchesiello2001}. 690 It is necessary to control the volume inside a domain when using open boundaries. 691 With fixed boundaries, it is enough to ensure that the total inflow/outflow has 692 reasonable values (either zero or a value compatible with an observed volume 693 balance). When using radiative boundary conditions it is necessary to have a 694 volume constraint because each open boundary works independently from the 695 others. The methodology used to control this volume is identical to the one 696 coded in the ROMS model \citep{Marchesiello2001}. 409 697 410 698 -
trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
r817 r994 13 13 14 14 % In climate modeling, it is often necessary to allow water masses that are separated by land to exchange properties. This situation arises in models when the grid mesh is too coarse to resolve narrow passageways that in reality provide crucial connections between water masses. For example, coarse grid spacing typically closes off the Mediterranean from the Atlantic at the Straits of Gibraltar. In this case, it is important for climate models to include the effects of salty water entering the Atlantic from the Mediterranean. Likewise, it is important for the Mediterranean to replenish its supply of water from the Atlantic to balance the net evaporation occurring over the Mediterranean region. 15 % also same problem as soon as important straits are not properly resolve by the mesh. For example in ORCA 1/4\r{}, this occors in several straits of the Indonesian archipelagos(Ombai, Lombok...).16 % We describe here the three methods that can be used in NEMO to handle such unproperly resolved straits.17 %The communication consists of mixing tracers and mass/volume between non-adjacent water columns. Momentum is not mixed. The scheme conserves total tracer content, and total volume (the latter in $z*$- or $s*$-coordinate), and maintains compatibility between the tracer and mass/volume budgets.18 19 % Note: such changes are so specific to a given configuration that no attempt ha ve been made to set them in a generic way. Nevertheless,example of how they can be set up is given in the ORCA2 configuration (\key{ORCA_R2}).15 % Same problem occurs as soon as important straits are not properly resolved by the mesh. For example it arises in ORCA 1/4\r{} in several straits of the Indonesian archipelago (Ombai, Lombok...). 16 % We describe here the three methods that can be used in \NEMO to handle such improperly resolved straits. 17 %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. 18 19 % Note: such changes are so specific to a given configuration that no attempt has been made to set them in a generic way. Nevertheless, an example of how they can be set up is given in the ORCA2 configuration (\key{ORCA_R2}). 20 20 21 21 % ------------------------------------------------------------------------------------------------------------- … … 27 27 $\bullet$ reduced scale factor, also called partially open face (Hallberg, personnal communication 2006) 28 28 %also called partially open cell or partial closed cells 29 $\bullet$ increase of the viscous boundary layer by local increase of the fmask value at the coast 29 30 $\bullet$ increase of the viscous boundary layer thickness by local increase of the fmask value at the coast 30 31 31 32 \colorbox{yellow}{Add a short description of scale factor changes staff and fmask increase} … … 36 37 \includegraphics[width=0.80\textwidth]{./Figures/Fig_Gibraltar.pdf} 37 38 \includegraphics[width=0.80\textwidth]{./Figures/Fig_Gibraltar2.pdf} 38 \caption {Example of the Gibraltar strait defined in a 1\r{} x 1\r{} mesh. \textit{Top}: using partially open cells. The meridional scale factor at $V$-point is reduced on both sides of the strait to account for the real width of the strait (about 20 km). Note that the scale factors of the strait $T$-point remains unchanged. \textit{Bottom}: using viscous boundary layers. The four fmask along the strait coastlines are set to a value larger than 4, $i.e.$ "strong" no-slip case (see Fig.\ref{Fig_LBC_shlat}) creating a large viscous boundary layer that allows a reduced transport through the strait.} 39 \caption {Example of the Gibraltar strait defined in a 1\r{} x 1\r{} mesh. 40 \textit{Top}: using partially open cells. The meridional scale factor at $v$-point 41 is reduced on both sides of the strait to account for the real width of the strait 42 (about 20 km). Note that the scale factors of the strait $T$-point remains unchanged. \textit{Bottom}: using viscous boundary layers. The four fmask parameters 43 along the strait coastlines are set to a value larger than 4, $i.e.$ "strong" no-slip 44 case (see Fig.\ref{Fig_LBC_shlat}) creating a large viscous boundary layer 45 that allows a reduced transport through the strait.} 39 46 \end{center} \end{figure} 40 47 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 65 72 \label{MISC_zoom} 66 73 67 The sub-domain functionality, also improperly called zoom option (improperly because it is not associated with a change in model resolution) is a quite simple function that allows to 68 perform a simulation over a sub-domain of an already defined configuration 69 (i.e. without defining a new set of mesh, initial state and forcings). This 70 option can be useful for testing the user setting of surface boundary 71 conditions, or the initial ocean state of a huge ocean model configuration 72 while having a small central memory requirement. It can also be used to 73 easily test specific physics in a sub-domain (for example, test of the 74 coupling between sea-ice and ocean model over the Arctic or Antarctic ocean 75 as they are set in the global ocean version of OPA \citep{Madec1996}. 76 In standard, this option does not include any specific treatment for ocean 77 boundaries of the sub-domain: they are considered as artificial vertical 78 walls. Nevertheless, it is quite easy to add a restoring term toward a 79 climatology in the vicinity of those boundaries (see \S\ref{TRA_dmp}). 74 The sub-domain functionality, also improperly called the zoom option 75 (improperly because it is not associated with a change in model resolution) 76 is a quite simple function that allows a simulation over a sub-domain of an 77 already defined configuration ($i.e.$ without defining a new mesh, initial 78 state and forcings). This option can be useful for testing the user settings 79 of surface boundary conditions, or the initial ocean state of a huge ocean 80 model configuration while having a small computer memory requirement. 81 It can also be used to easily test specific physics in a sub-domain (for example, 82 see \citep{Madec1996} for a test of the coupling used in the global ocean 83 version of OPA between sea-ice and ocean model over the Arctic or Antarctic 84 ocean, using a sub-domain). In the standard model, this option does not 85 include any specific treatment for the ocean boundaries of the sub-domain: 86 they are considered as artificial vertical walls. Nevertheless, it is quite easy 87 to add a restoring term toward a climatology in the vicinity of such boundaries 88 (see \S\ref{TRA_dmp}). 80 89 81 90 In order to easily define a sub-domain over which the computation can be 82 91 performed, the dimension of all input arrays (ocean mesh, bathymetry, 83 forcing, initial state, ...) are defined as \jp{jpidta}, \jp{jpjdta} and \jp{jpkdta} (\mdl{par\_oce} module), while the 84 computational domain is defined through \jp{jpiglo}, \jp{jpjglo} and \jp{jpk} (\mdl{par\_oce} module). When running the model 85 over the whole configuration, the user set \jp{jpiglo}=\jp{jpidta} \jp{jpjglo}=\jp{jpjdta }and \jp{jpk}=\jp{jpkdta}. When running the model 86 over a sub-domain, the user have to provide the size of the sub-domain, 87 (\jp{jpiglo}, \jp{jpjglo}, \jp{jpkglo}), and the indices of the south western corner as \jp{jpizoom} and \jp{jpjzoom} in the \mdl{par\_oce} module (Fig.~\ref{Fig_LBC_zoom}). 88 89 Note that a third set of dimension exist, \jp{jpi}, \jp{jpj} and \jp{jpk} which is actually used to 90 perform the computation. It is set by default to \jp{jpi}=\jp{jpjglo} and \jp{jpj}=\jp{jpjglo}, except for massively 91 parallel computing where the computational domain is laid out on local 92 processor memories following a 2D horizontal splitting (see {\S}IV.2-c) 92 forcing, initial state, ...) are defined as \jp{jpidta}, \jp{jpjdta} and \jp{jpkdta} 93 (\mdl{par\_oce} module), while the computational domain is defined through 94 \jp{jpiglo}, \jp{jpjglo} and \jp{jpk} (\mdl{par\_oce} module). When running the 95 model over the whole domain, the user sets \jp{jpiglo}=\jp{jpidta} \jp{jpjglo}=\jp{jpjdta} 96 and \jp{jpk}=\jp{jpkdta}. When running the model over a sub-domain, the user 97 has to provide the size of the sub-domain, (\jp{jpiglo}, \jp{jpjglo}, \jp{jpkglo}), 98 and the indices of the south western corner as \jp{jpizoom} and \jp{jpjzoom} in 99 the \mdl{par\_oce} module (Fig.~\ref{Fig_LBC_zoom}). 100 101 Note that a third set of dimensions exist, \jp{jpi}, \jp{jpj} and \jp{jpk} which is 102 actually used to perform the computation. It is set by default to \jp{jpi}=\jp{jpjglo} 103 and \jp{jpj}=\jp{jpjglo}, except for massively parallel computing where the 104 computational domain is laid out on local processor memories following a 2D 105 horizontal splitting. % (see {\S}IV.2-c) ref to the section to be updated 93 106 94 107 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 95 108 \begin{figure}[!ht] \label{Fig_LBC_zoom} \begin{center} 96 109 \includegraphics[width=0.90\textwidth]{./Figures/Fig_LBC_zoom.pdf} 97 \caption { position of a model domain compared to the data input domain when the zoom functionality is used.}110 \caption {Position of a model domain compared to the data input domain when the zoom functionality is used.} 98 111 \end{center} \end{figure} 99 112 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 106 119 \label{MISC_1D} 107 120 108 The 1D model is a stand alone water column based on the 3D NEMO system. It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers or a biogeochemical model. It is set up by defining the \key{cfg\_1d} CPP key. This 1D model is a very useful tool \textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes ; \textit{(b)} to investigate suitable parameterizations of unresolved turbulence (wind steering, langmuir circulation, skin layers) ; \textit{(c)} to compare the behaviour of different mixing vertical scheme ; \textit{(d)} to perform sensitivity study to the vertical diffusion on a particular point of the ocean global domain ; \textit{(d)} to access to specific diagnostics, aside from the standard model variables, because of having small in core memory requirement. 109 110 The methodology is based on the use of the zoom functionality (see \S\ref{MISC_zoom}), and the addition of some specific routines. There is no need of defining a new set of mesh, bathymetry, initial state and forcing, as the 1D model will use those of the configuration it is a zoom of. 111 112 121 The 1D model option simulates a stand alone water column within the 3D \NEMO 122 system. It can be applied to the ocean alone or to the ocean-ice system and can 123 include passive tracers or a biogeochemical model. It is set up by defining the 124 \key{cfg\_1d} CPP key. This 1D model is a very useful tool \textit{(a)} to learn 125 about the physics and numerical treatment of vertical mixing processes ; \textit{(b)} 126 to investigate suitable parameterizations of unresolved turbulence (wind steering, 127 langmuir circulation, skin layers) ; \textit{(c)} to compare the behaviour of different 128 vertical mixing schemes ; \textit{(d)} to perform sensitivity studies on the vertical 129 diffusion at a particular point of the ocean global domain ; \textit{(d)} to produce 130 extra diagnostics, without the large memory requirement of the full 3D model. 131 132 The methodology is based on the use of the zoom functionality (see 133 \S\ref{MISC_zoom}), with some extra routines. There is no need to define a new 134 mesh, bathymetry, initial state or forcing, since the 1D model will use those of the 135 configuration it is a zoom of. 113 136 114 137 % ================================================================ … … 121 144 %-------------------------------------------------------------------------------------------------------------- 122 145 146 % Steven update just bellow... not sure it is what I vant to say 147 %Searching an equilibrium state with an ocean model requires repeated very long time 148 %integrations (each of a few thousand years for a global model). Due to the size of 123 149 Searching an equilibrium state with an ocean model requires very long time 124 150 integration (a few thousand years for a global model). Due to the size of 125 the time step required for numerical stability consideration(less than a126 few hours), this is usually requires a large elapse time. In orderovercome127 this problem, \citet{Bryan1984} introduces a technique that allowsto accelerate128 the spin up to the equilibrium. It consists in usinga larger time step in151 the time step required for numerical stability (less than a 152 few hours), this usually requires a large elapsed time. In order to overcome 153 this problem, \citet{Bryan1984} introduces a technique that is intended to accelerate 154 the spin up to equilibrium. It uses a larger time step in 129 155 the thermodynamic evolution equations than in the dynamic evolution 130 156 equations. It does not affect the equilibrium solution but modifies the 131 157 trajectory to reach it. 132 158 133 The acceleration of convergence is used when \np{nn\_acc}=1. In that case, $\Delta t=rdt$ is 134 the time step of dynamics while $\widetilde{\Delta t}=rdttra$ is the tracer time-step. Both are settled from \np{rdt} and \np{rdttra} namelist parameters. The set of prognostic equations to solve becomes: 159 The acceleration of convergence option is used when \np{nn\_acc}=1. In that case, 160 $\Delta t=rdt$ is the time step of dynamics while $\widetilde{\Delta t}=rdttra$ is the 161 tracer time-step. Their values are set from the \np{rdt} and \np{rdttra} namelist 162 parameters. The set of prognostic equations to solve becomes: 135 163 \begin{equation} \label{Eq_acc} 136 164 \begin{split} … … 142 170 \end{equation} 143 171 144 \citet{Bryan1984} has analysed the consequences of this distorted physics. Free145 waves have a slower phase speed, their meridional structure is slightly172 \citet{Bryan1984} has examined the consequences of this distorted physics. 173 Free waves have a slower phase speed, their meridional structure is slightly 146 174 modified, and the growth rate of baroclinically unstable waves is reduced 147 but there isa wider range of instability. This technique is efficient for148 searching an equilibrium state in coarse resolution models. However its175 but with a wider range of instability. This technique is efficient for 176 searching for an equilibrium state in coarse resolution models. However its 149 177 application is not suitable for many oceanic problems: it cannot be used for 150 178 transient or time evolving problems (in particular, it is very questionable 151 to keep this technique when usinga seasonal cycle in the forcing fields),179 to use this technique when there is a seasonal cycle in the forcing fields), 152 180 and it cannot be used in high-resolution models where baroclinically 153 181 unstable processes are important. Moreover, the vertical variation of 154 $\ Delta \tilde {t}$ implies that the heat and salt contents are no more182 $\widetilde{ \Delta t}$ implies that the heat and salt contents are no longer 155 183 conserved due to the vertical coupling of the ocean level through both 156 184 advection and diffusion. 185 % first mention of depth varying tilda-delta-t! and namelist parameter associated to that are to be described 157 186 158 187 … … 160 189 % Model optimisation, Control Print and Benchmark 161 190 % ================================================================ 162 \section{Model optimisation, Control Print and Benchmark}191 \section{Model Optimisation, Control Print and Benchmark} 163 192 \label{MISC_opt} 164 193 %--------------------------------------------namdom------------------------------------------------------- … … 166 195 %-------------------------------------------------------------------------------------------------------------- 167 196 168 Three points to be described here: 197 % \gmcomment{why not make these bullets into subsections?} 198 199 Three issues to be described here: 169 200 170 201 $\bullet$ Vector and memory optimisation: 171 202 172 \key{vectopt\_loop} enable the internal loop collapse, a very efficient way to increase the length of vector and thus speed up the model on vector computers. 203 \key{vectopt\_loop} enables internal loop collapse, a very efficient way to increase 204 the length of vector calculations and thus speed up the model on vector computers. 173 205 174 Add here also one word on NPROMA technique that has been found useless, since compiler have made significant progress during the last decade.206 % Add here also one word on NPROMA technique that has been found useless, since compiler have made significant progress during the last decade. 175 207 176 Add also one word on NEC specific optimisation (Novercheck option for example) 177 178 \key{vectopt\_memory} has been introduced in order to reduce the memory requirement of the model. This is obviously done by incresing the CPU time requirement, as it suppress intermediate computation saved in in-core memory. This possibility have not been intensively used. In fact up to now, it only concern the TKE physics, in which, when \key{vectopt\_memory} is defined, the coefficient used for horizontal smoothing of $A_v^T$ and $A_v^m$ are no more computed once for all. This reduces the memory requirement by three 2D arrays. 208 % Add also one word on NEC specific optimisation (Novercheck option for example) 209 210 \key{vectopt\_memory} has been introduced in order to reduce the memory 211 requirement of the model. This is obviously done at the cost of increasing the 212 CPU time requirement, since it suppress intermediate computations which would 213 have been saved in in-core memory. This feature has not been intensively used. 214 In fact, currently it is only implemented for the TKE physics, in which, when 215 \key{vectopt\_memory} is defined, the coefficients used for horizontal smoothing 216 of $A_v^T$ and $A_v^m$ are no longer computed once and for all. This reduces 217 the memory requirement by three 2D arrays. 218 179 219 180 $\bullet$ Control print: describe here 4 things: 181 182 1- \np{ln\_ctl} : compute and print the inner domain averaged trends in all TRA, DYN LDF and ZDF modules. Very useful to detect the origin of an undesired change in model results. 183 184 2- also \np{ln\_ctl} but using the nictl. njctl. namelist parameters to check the origin of differences between mono and multi processor 185 186 3- \key{esopa} (to be rename key\_nemo) : also a option of model management. When defined, this key force the activation of all options and CPP keys. For example, all the tracer and momentum advection scheme are called ! There is therefore no physical meaning associated with model results. In fact, this option allows both the compilator and the model run to go through all the Fortran lines of the model. This allows to check is there is obvious compilation or running bugs for CPP options, and running bugs for namelist options. 187 188 5- last digit comparison (\np{nbit\_cmp}). In MPP simulation, the computation of sum of the whole domain is performed as the sum over all processors of the sum of the inner domain of each processor. This double sum nevr give exactly the same result as a single sum over the whole domain, due to truncature differences. The "bit comparison" option has been introduced in order to be able to check that mono-processor and multi-processor give exactly the same results. 189 190 191 $\bullet$ Benchmark (\np{nbench}). This option, defines a benchmark run based on GYRE configuration in which the resolution remains the same whatever the domain size is. This allows to run very large model domain by just changing the domain size (\jp{jpiglo}, \jp{jpjglo}) without adjusting neither the time-step nor the physical parametrisations. 220 $\bullet$ Control print %: describe here 4 things: 221 222 1- \np{ln\_ctl} : compute and print the trends averaged over the interior domain 223 in all TRA, DYN, LDF and ZDF modules. This option is very helpful when 224 diagnosing the origin of an undesired change in model results. 225 226 2- also \np{ln\_ctl} but using the nictl and njctl namelist parameters to check 227 the source of differences between mono and multi processor runs. 228 229 3- \key{esopa} (to be rename key\_nemo) : which is another option for model 230 management. When defined, this key forces the activation of all options and 231 CPP keys. For example, all tracer and momentum advection schemes are called! 232 There is therefore no physical meaning associated with the model results. 233 However, this option forces both the compiler and the model to run through 234 all the \textsc{Fortran} lines of the model. This allows the user to check for obvious 235 compilation or execution errors with all CPP options, and errors in namelist options. 236 237 3- \key{esopa} (to be rename key\_nemo) : which is another option for model 238 management. When defined, this key forces the activation of all options and CPP keys. 239 For example, all tracer and momentum advection schemes are called! There is 240 therefore no physical meaning associated with the model results. However, this option 241 forces both the compiler and the model to run through all the Fortran lines of the model. 242 This allows the user to check for obvious compilation or execution errors with all CPP 243 options, and errors in namelist options. 244 245 4- last digit comparison (\np{nbit\_cmp}). In an MPP simulation, the computation of 246 a sum over the whole domain is performed as the summation over all processors of 247 each of their sums over their interior domains. This double sum never gives exactly 248 the same result as a single sum over the whole domain, due to truncation differences. 249 The "bit comparison" option has been introduced in order to be able to check that 250 mono-processor and multi-processor runs give exactly the same results. 251 252 $\bullet$ Benchmark (\np{nbench}). This option defines a benchmark run based on 253 a GYRE configuration in which the resolution remains the same whatever the domain 254 size. This allows a very large model domain to be used, just by changing the domain 255 size (\jp{jpiglo}, \jp{jpjglo}) and without adjusting either the time-step or the physical 256 parameterizations. 192 257 193 258 … … 202 267 203 268 The computation of the surface pressure gradient with a rigid lid assumption 204 requires to compute $\partial_t \psi$, the time evolution of the 205 barotropic streamfunction. $\partial_t \psi$ is solution of an elliptic 206 equation (I.2.4) for which two solvers are available, a 207 Successive-Over-Relaxation (SOR) or a preconditioned conjugate gradient 208 (PCG) \citep{Madec1988, Madec1990}. The PCG is a very efficient 209 method for solving elliptic equations on vector computers. It is a fast and 210 rather easy to use method, which is an attractive feature for a large number 211 of ocean situations (variable bottom topography, complex coastal geometry, 212 variable grid spacing, islands, open or cyclic boundaries, etc ...). It does 213 not require the search of an optimal parameter as in the SOR method. 214 Nevertheless, the SOR has been kept because it is a linear solver, a very 215 useful property when using the adjoint model of OPA. 216 217 The surface pressure gradient is computed in \mdl{dynspg}. The 218 default option is the use of PCG or SOR 219 depending on \np{nsolv} (namelist parameter). 269 requires the computation of $\partial_t \psi$, the time evolution of the 270 barotropic streamfunction. $\partial_t \psi$ is the solution of an elliptic 271 equation \eqref{Eq_PE_psi} for which three solvers are available: a 272 Successive-Over-Relaxation scheme (SOR), a preconditioned conjugate gradient 273 scheme(PCG) \citep{Madec1988, Madec1990} and a Finite Elements Tearing and 274 Interconnecting scheme (FETI) \citep{Guyon_al_EP99, Guyon_al_CalPar99}. 275 The PCG is a very efficient method for solving elliptic equations on vector computers. 276 It is a fast and rather easy method to use; which are attractive features for a large 277 number of ocean situations (variable bottom topography, complex coastal geometry, 278 variable grid spacing, islands, open or cyclic boundaries, etc ...). It does not require 279 a search for an optimal parameter as in the SOR method. However, the SOR has 280 been retained because it is a linear solver, which is a very useful property when 281 using the adjoint model of \NEMO. The FETI solver is a very efficient method on 282 massively parallel computers. However, it has never been used since OPA 8.0. 283 The current version in \NEMO should not even successfully go through the 284 compilation phase. 285 286 The surface pressure gradient is computed in \mdl{dynspg}. The solver used 287 (PCG or SOR) depending on the value of \np{nsolv} (namelist parameter). 220 288 At each time step the time derivative of the barotropic streamfunction is 221 the solution of (II.2.3). Introducing the following coefficients:289 the solution of \eqref{Eq_psi_total}. Introducing the following coefficients: 222 290 223 291 \begin{equation} … … 229 297 \end{equation} 230 298 231 the five-point finite difference equ ivalent equation (II.2.3)can be rewritten as:299 the five-point finite difference equation \eqref{Eq_psi_total} can be rewritten as: 232 300 233 301 \begin{multline} \label{Eq_solmat} … … 257 325 \textbf{A}). (VII.5.1) can be rewritten as: 258 326 259 \begin{multline} 327 \begin{multline} \label{Eq_solmat_p} 260 328 A_{i,j}^{N} { \left( \frac{\partial \psi}{\partial t } \right) }_{i+1,j } 261 329 +\,A_{i,j}^{E} { \left( \frac{\partial \psi}{\partial t } \right) }_{i ,j+1} \\ … … 266 334 267 335 The SOR method used is an iterative method. Its algorithm can be summarised 268 as follows [see \citet{Haltiner1980} for further discussion]:269 270 initialisation (evaluate a first guess from formertime step computations)336 as follows [see \citet{Haltiner1980} for a further discussion]: 337 338 initialisation (evaluate a first guess from previous time step computations) 271 339 \begin{equation} 272 340 \left( \frac{\partial \psi}{\partial t } \right)_{i,j}^0 = 2{\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^t} - {\left( \frac{\partial \psi}{\partial t } \right)_{i,j}^{t-1}} 273 341 \end{equation} 274 iteration $n ,$ from $n=0 $until convergence, do :342 iteration $n$, from $n=0$ until convergence, do : 275 343 \begin{multline} 276 344 R_{i,j}^n … … 286 354 + \omega \;R_{i,j}^n 287 355 \end{equation} 288 where \textit{$\omega $ }satisfies $1\leq \omega \leq 2$. An optimal value exists for \textit{$\omega $ }which 289 accelerates significantly the convergence, but it has to be adjusted 290 empirically for each model domain, except for an uniform grid where an 291 analytical expression for \textit{$\omega $ }can be found \citep{Richtmyer1967}. This 292 parameter is defined as \textit{sor}, a \textbf{namelist} parameter. The convergence test is of the form: 356 where \textit{$\omega $ }satisfies $1\leq \omega \leq 2$. An optimal value exists for 357 \textit{$\omega$} which significantly accelerates the convergence, but it has to be 358 adjusted empirically for each model domain (except for a uniform grid where an 359 analytical expression for \textit{$\omega$} can be found \citep{Richtmyer1967}). 360 The value of $\omega$ is set using \textit{sor}, a \textbf{namelist} parameter. 361 The convergence test is of the form: 293 362 \begin{equation} 294 363 \delta = \frac{\sum\limits_{i,j}{R_{i,j}^n}{R_{i,j}^n}}{\sum\limits_{i,j}{ \tilde B_{i,j}^n}{\tilde B_{i,j}^n}} \leq \epsilon 295 364 \end{equation} 296 where $\epsilon $ is the absolute precision that is required. It is highly recommended 297 to use a $\epsilon $ smaller or equal to $10^{-2}$ as larger values may leads 298 to numerically induced basin scale barotropic oscillations, and to use a two 299 or three order of magnitude smaller value in eddy resolving configuration. 300 The precision of the solver is not only a numerical parameter, but 301 influences the physics. Indeed, the zero change of kinetic energy due to the 302 work of surface pressure forces in rigid-lid approximation is only achieved 303 at the precision demanded on the solver ({\S}~C.1-e). The precision is 304 specified by setting \textit{eps} (\textbf{namelist parameter}). In addition, two other tests are used to stop the iterative 305 algorithm. They concern the number of iterations and the module of the right 306 hand side. If the former exceeds a specified value, \textit{nmax} (\textbf{namelist parameter}), or the latter is greater than 307 $10^{15}$, the whole model computation is stopped while the last 308 computed time step fields are saved in the standard output file. In both 309 cases, this usually indicates that there is something wrong in the model 310 configuration (error in the mesh, the initial state, the input forcing, or 311 the magnitude of the time step or of the mixing coefficients). A typical 312 value of $nmax$ is a few hundred for $\epsilon = 10^{-2}$, increasing to a few 313 thousand for $\epsilon = 10^{-12}$. 314 315 The vectorization of the SOR algorithm is not straightforward. (VII.5.4) 316 contains two linear recurrences on $i$ and $j$. This inhibits the vectorisation 317 ({\S}~IV.2-a). Therefore (VII.5.4a) has been rewritten as: 365 where $\epsilon$ is the absolute precision that is required. It is recommended 366 that a value smaller or equal to $10^{-3}$ is used for $\epsilon$ since larger 367 values may lead to numerically induced basin scale barotropic oscillations. In fact, 368 for an eddy resolving configuration or in a filtered free surface case, a value three 369 orders of magnitude smaller than this should be used. The precision is specified by 370 setting \textit{eps} (\textbf{namelist parameter}). In addition, two other tests are 371 used to halt the iterative algorithm. They involve the number of iterations and the 372 modulus of the right hand side. If the former exceeds a specified value, \textit{nmax} 373 (\textbf{namelist parameter}), or the latter is greater than $10^{15}$, the whole 374 model computation is stopped and the last computed time step fields are saved 375 in the standard output file. In both cases, this usually indicates that there is something 376 wrong in the model configuration (an error in the mesh, the initial state, the input forcing, 377 or the magnitude of the time step or of the mixing coefficients). A typical value of 378 $nmax$ is a few hundred when $\epsilon = 10^{-6}$, increasing to a few 379 thousand when $\epsilon = 10^{-12}$. 380 The vectorization of the SOR algorithm is not straightforward. The scheme 381 contains two linear recurrences on $i$ and $j$. This inhibits the vectorisation. 382 Therefore it has been rewritten as: 318 383 319 384 \begin{multline} … … 334 399 \end{equation} 335 400 336 If the three equations in (VII.5.6) are solved inside the same do-loop, 337 (VII.5.4a) and (VII.5.6) are strictly equivalent. In the model they are 338 solved successively over the whole domain. The convergence is slower but 339 (VII.5.6a) and (VII.5.6b) are vector loops on $i$-index (the inner loop) and 340 (VII.5.6c) is adapted to Cray vector computers by using a routine from the 341 Cray library (namely the FOLR routine) to solve the first-order linear 342 recurrence. The SOR method is very flexible and can be used under a wide 401 The SOR method is very flexible and can be used under a wide 343 402 range of conditions, including irregular boundaries, interior boundary 344 403 points, etc. Proofs of convergence, etc. may be found in the standard 345 numerical methods texts for partial equations.404 numerical methods texts for partial differential equations. 346 405 347 406 % ------------------------------------------------------------------------------------------------------------- … … 356 415 357 416 \textbf{A} is a definite positive symmetric matrix, thus solving the linear 358 system (VII.5.1)is equivalent to the minimisation of a quadratic417 system \eqref{Eq_solmat} is equivalent to the minimisation of a quadratic 359 418 functional: 360 419 \begin{equation*} … … 364 423 \end{equation*} 365 424 where $\langle , \rangle$ is the canonical dot product. The idea of the 366 conjugate gradient method is to search the solution in the following 367 iterative way: assuming that $\textbf{x}^n$ is obtained, $\textbf{x}^{n+1}$ is searched under the form $\textbf {x}^{n+1}={\textbf {x}}^n+\alpha^n{\textbf {d}}^n$ which satisfies: 425 conjugate gradient method is to search for the solution in the following 426 iterative way: assuming that $\textbf{x}^n$ has been obtained, $\textbf{x}^{n+1}$ 427 is found from $\textbf {x}^{n+1}={\textbf {x}}^n+\alpha^n{\textbf {d}}^n$ which satisfies: 368 428 \begin{equation*} 369 {\textbf{ x}}^{n+1}=\text{inf} _{{\textbf{ y}}={\textbf{ x}}^n+\alpha \,{\textbf{ d}}^n} \,\phi ({\textbf{ y}})\;\;\Leftrightarrow \;\;\frac{d\phi }{d\alpha}=0429 {\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 370 430 \end{equation*} 371 431 and expressing $\phi (\textbf{y})$ as a function of \textit{$\alpha $}, we obtain the value that minimises the functional: … … 373 433 \alpha ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{ A d}^n, \textbf{d}^n} \rangle 374 434 \end{equation*} 375 where $\textbf{r}^n = \textbf{b}-\textbf{A x}^n = \textbf{A} (\textbf{x}-\textbf{x}^n)$ is the error at rank $n$. The choice of the descent vector $\textbf{d}^n$ depends on the error: $\textbf{d}^n = \textbf{r}^n + \beta^n \,\textbf{d}^{n-1}$. $\beta ^n$ is searched such that the descent vectors form an orthogonal base for the dot product linked to \textbf{A}. Expressing the condition 435 where $\textbf{r}^n = \textbf{b}-\textbf{A x}^n = \textbf{A} (\textbf{x}-\textbf{x}^n)$ 436 is the error at rank $n$. The descent vector $\textbf{d}^n$ s chosen to be dependent 437 on the error: $\textbf{d}^n = \textbf{r}^n + \beta^n \,\textbf{d}^{n-1}$. $\beta ^n$ 438 is searched such that the descent vectors form an orthogonal basis for the dot 439 product linked to \textbf{A}. Expressing the condition 376 440 $\langle \textbf{A d}^n, \textbf{d}^{n-1} \rangle = 0$ the value of $\beta ^n$ is found: 377 $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{r}^{n-1}, \textbf{r}^{n-1}} \rangle$. As a result, the errors $ \textbf{r}^n$ form an orthogonal base378 for the canonic dot product while the descent vectors $\textbf{d}^n$ form441 $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{r}^{n-1}, \textbf{r}^{n-1}} \rangle$. As a result, the errors $ \textbf{r}^n$ form an orthogonal 442 base for the canonic dot product while the descent vectors $\textbf{d}^n$ form 379 443 an orthogonal base for the dot product linked to \textbf{A}. The resulting 380 444 algorithm is thus the following one: 381 445 382 446 initialisation : 383 384 447 385 448 \begin{equation*} … … 415 478 \end{equation} 416 479 where $\epsilon $ is the absolute precision that is required. As for the SOR algorithm, 417 both the PCG algorithm and the whole model computation are stopped when the 418 number of iteration, $nmax$, or the module of the right hand side exceeds a 419 specified value (see {\S}~VII.5-a for further discussion). The precision and 420 the maximum number of iteration are specified by setting $eps$ and $nmax$ (\textbf{namelist parameters}). 421 422 It can be demonstrated that the algorithm is optimal, provides the exact 480 the whole model computation is stopped when the number of iterations, $nmax$, or 481 the modulus of the right hand side of the convergence equation exceeds a 482 specified value (see \S\ref{MISC_solsor} for a further discussion). The required 483 precision and the maximum number of iterations allowed are specified by setting 484 $eps$ and $nmax$ (\textbf{namelist parameters}). 485 486 It can be demonstrated that the above algorithm is optimal, provides the exact 423 487 solution in a number of iterations equal to the size of the matrix, and that 424 the convergence rate is all the more fast as the matrix is closed to 425 identity (i.e. the eigen values are closed to 1). Therefore, it is more 426 efficient to solve a better conditioned system which has the same solution. 427 For that purpose, we introduce a preconditioning matrix \textbf{Q 428 }which is an approximation of \textbf{A} but much easier to invert 429 than \textbf{A} and solve the system: 430 \begin{equation} 488 the convergence rate is faster as the matrix is closer to the identity matrix, 489 $i.e.$ its eigenvalues are closer to 1. Therefore, it is more efficient to solve a 490 better conditioned system which has the same solution. For that purpose, we 491 introduce a preconditioning matrix \textbf{Q} which is an approximation of 492 \textbf{A} but much easier to invert than \textbf{A}, and solve the system: 493 \begin{equation} \label{Eq_pmat} 431 494 \textbf{Q}^{-1} \textbf{A x} = \textbf{Q}^{-1} \textbf{b} 432 495 \end{equation} 433 496 434 The same algorithm can be used to solve (VII.5.7) if instead of the 435 canonical dot product the following one is used: ${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, and if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ are substituted to \textbf{b} and \textbf{A} \citep{Madec1988}. In OPA, \textbf{Q} is chosen as the diagonal 436 of\textbf{ A}, i.e. the simplest form for \textbf{Q} so that it can be 437 easily inverted. In this case, the discrete formulation of (VII.5.8) is in 438 fact given by (VII.5.2) and thus the matrix and right hand side are computed 439 independently from the solver used. 497 The same algorithm can be used to solve \eqref{Eq_pmat} if instead of the 498 canonical dot product the following one is used: 499 ${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, 500 and if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ are substituted to \textbf{b} and \textbf{A} \citep{Madec1988}. 501 In \NEMO, \textbf{Q} is chosen as the diagonal of \textbf{ A}, i.e. the simplest form for \textbf{Q} so that it can be easily inverted. In this case, the discrete formulation of 502 \eqref{Eq_pmat} is in fact given by \eqref{Eq_solmat_p} and thus the matrix and 503 right hand side are computed independently from the solver used. 440 504 441 505 % ------------------------------------------------------------------------------------------------------------- … … 445 509 \label{MISC_solfeti} 446 510 447 FETI is a powerfull solver that was developed by Marc Guyon (ref !!!). It as been just converted from Fortan 77 to 90, but never successfully tested after that. Since then, nobody has found enough time to further investigate the implmenntationof FETI and debug it. 448 449 Its main advantaged is to save a lot of CPU time compared to SOR or PCG algorithm. Nevertheless, its main drawbacks is that the solution is depends on the domain decomposition chosen. Using a different number of processors, the solution is the same at the precision required, but not the same at the computer precision. This make it hard to debug. 511 FETI is a powerful solver that was developed by Marc Guyon \citep{Guyon_al_EP99, 512 Guyon_al_CalPar99}. It has been converted from Fortan 77 to 90, but never 513 successfully tested after that. 514 515 Its main advantage is to save a lot of CPU time when compared to the SOR or PCG 516 algorithms. However, its main drawback is that the solution is dependent on the 517 domain decomposition chosen. Using a different number of processors, the solution 518 is the same at the precision required, but not the same at the computer precision. 519 This make it hard to debug. 450 520 451 521 % ------------------------------------------------------------------------------------------------------------- … … 455 525 \label{MISC_solisl} 456 526 457 The boundary condition used for both solvers is that the time derivative of458 the barotropic streamfunction is zero along all the coastlines. When islands459 are present in the model domain, additional computations must be done to460 determine the barotropic streamfunction with the correct boundary527 The boundary condition used for both recommended solvers is that the time 528 derivative of the barotropic streamfunction is zero along all the coastlines. 529 When islands are present in the model domain, additional computations must 530 be performed to determine the barotropic streamfunction with the correct boundary 461 531 conditions. This is detailed below. 462 532 463 The model does not recognise itself the islands which must be defined using 464 bathymetry information, i.e. $mbathy$ array, equals $-1$ over the first island, $-2$ over 465 the second, ... , $-N$ over the $N^{th}$ island (see {\S}~VII.2-b). The model 466 determines the position of island grid-points and defines a closed contour 467 around each island which will be used to compute the circulation around each 468 island. The closed contour is formed of the ocean grid-points which are the 469 closest to the island. 470 471 First, the island barotropic streamfunctions $\psi_n$ are computed 472 using the PCG (they are solutions of (VII.5.1) with the right-hand side 473 equals to zero and with $\psi_n = 1$ along the island $n$ and $\psi_n = 0$ along the other coastlines) (Note that specifying $1$ as boundary 474 condition on an island for $\psi$ is equivalent to define a 475 specific right hand side for (VII.5.1)~). The precision of this computation 476 can be very high since it is done once. The absolute precision, $epsisl$, and the 477 maximum number of iteration, $nmisl$, are the \textbf{namelist parameters} used for that 478 computation. Their typical values are $epsisl = 10^{-10}$ and $nmisl = 4000$. Then the island matrix A is computed from (I.2.8) and reversed. At 479 each time step, $\psi_0$, the solution of (I.2.4) with $\psi_0 = 0$ along all coastlines, is computed using either SOR or PCG. It has to 480 be noted that the first guess of this computation is defined as in (VII.5.3) 481 except that $\partial_t \psi_0$ is used, not $\partial_t \psi$. Indeed we are 482 computing $\partial_t \psi_0$ which is usually far different from $\partial_t \psi$. Then, it is easy to find the time evolution of the barotropic 483 streamfunction on each island and to deduce $\partial_t \psi$ using (I.2.9) 484 in order to compute the surface pressure gradient. Note that the value of 485 the barotropic streamfunction itself is also computed as the time 486 integration of $\partial_t \psi$ for further diagnostics. 533 The model does not have specialised code for islands. They must instead be 534 identified to the solvers by the user via bathymetry information, i.e. the $mbathy$ 535 array should equal $-1$ over the first island, $-2$ over the second, ... , $-N$ over 536 the $N^{th}$ island. The model determines the position of island grid-points and 537 defines a closed contour around each island which is used to compute the circulation 538 around each one. The closed contour is formed from the ocean grid-points which 539 are the closest to the island. 540 541 First, the island barotropic streamfunctions $\psi_n$ are computed using the SOR 542 or PCG scheme (they are solutions of \eqref{Eq_solmat} with the right-hand side 543 equal to zero and with $\psi_n = 1$ along the island $n$ and $\psi_n = 0$ along 544 the other coastlines) (Note that specifying $1$ as boundary condition on an island 545 for $\psi$ is equivalent to defining a specific right hand side for \eqref{Eq_solmat}). 546 The requested precision of this computation can be very high since it is only 547 performed once. The absolute precision, $epsisl$, and the maximum number of 548 iterations allowed, $nmisl$, are the \textbf{namelist parameters} used for this 549 computation. Their typical values are $epsisl = 10^{-10}$ and $nmisl = 4000$. 550 Then the island matrix A is computed from \eqref{Eq_PE_psi_isl_matrix} and 551 inverted. At each time step, $\psi_0$, the solution of \eqref{Eq_solmat} with $\psi = 0$ 552 along all coastlines, is computed using either SOR or PCG. (It should 553 be noted that the first guess of this computation remains as usual except that 554 $\partial_t \psi_0$ is used, instead of $\partial_t \psi$. Indeed we are 555 computing $\partial_t \psi_0$ which is usually very different from $\partial_t \psi$.) 556 Then, it is easy to find the time evolution of the barotropic streamfunction on each 557 island and to deduce $\partial_t \psi$, and to use it to compute the surface pressure 558 gradient. Note that the value of the barotropic streamfunction itself is also computed 559 as the time integration of $\partial_t \psi$ for further diagnostics. 487 560 488 561 % ================================================================ … … 503 576 and the output file(s). The restart file is used internally by the code when 504 577 the user wants to start the model with initial conditions defined by a 505 previous simulation. It contains all the information that is necessary not506 to introduce changes in the model results (even at the computer precision)507 between a run performed with several restarts and the same run performed in508 one time. It has to be noticed that this requires that the restart file509 contain stwo consecutive time steps for all the prognostic variables, and510 that it is save in the same binary format as the one used by the computer to511 calculate (in particular, 32 bits binary IEEE format for this file must not512 be used). The output listing and file(s) aredefined but should be checked578 previous simulation. It contains all the information that is necessary in 579 order for there to be no changes in the model results (even at the computer 580 precision) between a run performed with several restarts and the same run 581 performed in one step. It should be noted that this requires that the restart file 582 contain two consecutive time steps for all the prognostic variables, and 583 that it is saved in the same binary format as the one used by the computer 584 that is to read it (in particular, 32 bits binary IEEE format must not be used for 585 this file). The output listing and file(s) are predefined but should be checked 513 586 and eventually adapted to the user's needs. The output listing is stored in 514 the $ocean.output$ file. The information are printed all the way through the code on the 515 logical unit $numout$. To locate these prints, use the UNIX command "$grep -i numout^*$" in the source code directory. 587 the $ocean.output$ file. The information is printed from within the code on the 588 logical unit $numout$. To locate these prints, use the UNIX command 589 "$grep -i numout^*$" in the source code directory. 516 590 517 591 In the standard configuration, the user will find the model results in two 518 output files for every time-step where output is demanded: a \textbf{VO} 519 file containing all the three dimensional fields in logical unit $numwvo$, and a 520 \textbf{SO} file containing all the two dimensional (horizontal) fields in 521 logical unit $numwso$. These outputs are defined in the $diawri.F$ routine. The standard and available on-line diagnostics are described in {\S}III-12-c. 522 523 The default output is 32 bits binary IEEE format, compatible with the 524 Vairmer software (see the Climate Modelling and global Change team WEB 525 server at CERFACS: http://www.cerfacs.fr). The model's reference directory 526 also contains a visualisation tool based on \textbf{NCAR Graphics 527 }(http://ngwww.ucar.edu). If a user has access to the NCAR software, he or 528 she can copy the \textbf{LODMODEL/UTILS/OPADRA} directory from the reference 529 and, following the \textbf{README}, create the graphic outputs from the 530 model's results. 592 output files containing values for every time-step where output is demanded: 593 a \textbf{VO} file containing all the three dimensional fields in logical unit 594 $numwvo$, and a \textbf{SO} file containing all the two dimensional (horizontal) 595 fields in logical unit $numwso$. These outputs are defined in the $diawri.F$ 596 routine. The default and user-selectable diagnostics are described in {\S}III-12-c. 597 598 The default output (for all output files apart from the restart file) is 32 bits binary 599 IEEE format, compatible with the Vairmer software (see the Climate Modelling 600 and global Change team WEB server at CERFACS: http://www.cerfacs.fr). 601 The model's reference directory also contains a visualisation tool based on 602 \textbf{NCAR Graphics} (http://ngwww.ucar.edu). If a user has access to the 603 NCAR software, he or she can copy the \textbf{LODMODEL/UTILS/OPADRA} 604 directory from the reference and, following the \textbf{README}, create 605 graphical outputs from the model's results. 531 606 } 532 607 … … 534 609 % Tracer/Dynamics Trends 535 610 % ------------------------------------------------------------------------------------------------------------- 536 \subsection{Tracer/Dynamics Trends (\key{trdlmd}, \key{diatrdtra} )}611 \subsection{Tracer/Dynamics Trends (\key{trdlmd}, \key{diatrdtra}, \key{diatrddyn})} 537 612 \label{MISC_tratrd} 538 613 539 614 %to be udated this corresponds to OPA8 540 When \key{diatrddyn} and/or \key{diatrddyn} cpp 541 variables are defined, each trends of the dynamics and/or temperature and 542 salinity time evolution equations are stored in three-dimensional arrays 543 just after their computation (i.e. at the end of each $dyn\cdots .F$ 544 and/or $tra\cdots .F$ routines). These trends are then used in diagnostic 545 routines $diadyn.F$ and $diatra.F$ respectively. In standard, these routines check the basin averaged properties of the momentum and tracer equations every \textit{ntrd } time-steps (\textbf{namelist parameter}). These routines 546 are given as an example; they must be adapted by the user to his/her 547 desiderata. 548 549 These two options imply the definition of several arrays in the in core 550 memory, increasing quite sensitively the code memory requirements (search 551 for \key{diatrddyn} or 552 \key{diatrdtra}in \textit{common.h} file) 615 When \key{diatrddyn} and/or \key{diatrddyn} cpp variables are defined, each 616 trend of the dynamics and/or temperature and salinity time evolution equations 617 is stored in three-dimensional arrays just after their computation ($i.e.$ at the end 618 of each $dyn\cdots .F90$ and/or $tra\cdots .F90$ routine). These trends are then 619 used in diagnostic routines $diadyn.F90$ and $diatra.F90$ respectively. 620 In the standard model, these routines check the basin averaged properties of 621 the momentum and tracer equations every \textit{ntrd } time-steps (\textbf{namelist 622 parameter}). These routines are supplied as an example; they must be adapted by 623 the user to his/her requirements. 624 625 These two options imply the creation of several extra arrays in the in-core 626 memory, increasing quite seriously the code memory requirements. 553 627 554 628 % ------------------------------------------------------------------------------------------------------------- … … 561 635 %-------------------------------------------------------------------------------------------------------------- 562 636 563 \colorbox{yellow}{a description is to be added here} 637 The on-line computation of floats adevected either by the three dimensional velocity 638 field or constraint to remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project. The algorithm used is based on 639 the work of \cite{Blanke_Raynaud_JPO97}. (see also the web site describing the off-line 640 use of this marvellous diagnostic tool (http://stockage.univ-brest.fr/~grima/Ariane/). 564 641 565 642 % ------------------------------------------------------------------------------------------------------------- … … 569 646 \label{MISC_diag_others} 570 647 571 %To be updated this corresponds to OPA 8 572 573 Aside from the standard model variables, some other diagnostics are computed 574 on-line or can be added in the model. The available ready-to-add diagnostics 575 can be found in DIA. Among the available diagnostics one can 576 quote: 648 %To be updated this mainly corresponds to OPA 8 649 650 Aside from the standard model variables, other diagnostics can be computed 651 on-line or can be added to the model. The available ready-to-add diagnostics 652 routines can be found in directory DIA. Among the available diagnostics are: 577 653 578 654 - the mixed layer depth (based on a density criterion) (\mdl{diamxl}) -
trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex
r817 r994 525 525 \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} 526 526 + \frac{\partial \left( {e_1 \,v\,u} \right)}{\partial j} \right) 527 - \frac{1}{e_3 } \frac{\partial \left( {w\,u} \right)}{\partial k} \\527 - \frac{1}{e_3 }\frac{\partial \left( { w\,u} \right)}{\partial k} \\ 528 528 - \frac{1}{e_1 }\frac{\partial}{\partial i}\left( \frac{p_s+p_h }{\rho _o} \right) 529 529 + D_u^{\vect{U}} + F_u^{\vect{U}} … … 536 536 \frac{1}{e_1 \; e_2} \left( 537 537 \frac{\partial \left( {e_2 \,u\,v} \right)}{\partial i} 538 + \frac{\partial \left( {e_1 \, u\,v} \right)}{\partial j} \right)539 - \frac{1}{e_3 } \frac{\partial \left( {w\,v} \right)}{\partial k} \\538 + \frac{\partial \left( {e_1 \,v\,v} \right)}{\partial j} \right) 539 - \frac{1}{e_3 } \frac{\partial \left( { w\,v} \right)}{\partial k} \\ 540 540 - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o} \right) 541 541 + D_v^{\vect{U}} + F_v^{\vect{U}} -
trunk/DOC/TexFiles/Chapters/Chap_SBC.tex
r817 r994 6 6 \minitoc 7 7 8 \ begin{verbatim}9 At the time of this writing, the new surface module 10 that is described in this chapter (SBC) is not yet part 11 of the current distribution. The current way to specify 12 the surface boundary condition is such a mess that we 13 did not attempt to describe it. Nevertheless, apart from 14 the way the surface forcing is implemented, the infor- 15 mation given here are relevant for a NEMO v2.3 user. 16 \ end{verbatim}17 18 The ocean needs 7 fields as surface boundary condition: 19 20 The two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$ 21 22 The incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$ 23 24 The surface freshwater budget $\left( {\text{EMP}\;,\;\text{EMP}_S } \right)$ 25 26 \colorbox {yellow}{ The river runoffs (RUNOFF)} 27 28 Four different ways are offered to provide those 7 fields to the ocean: an29 analytical formulation, a flux formulation, a bulk formulae formulation30 (CORE or CLIO bulk formulae) and a coupled formulation (exchanges with a31 atmospheric model via OASIS coupler). In addition, the resulting fields can32 be further modified on used demand via several namelist options. These options33 control the addition of a surface restoring term to observed SST and/or SSS,34 the modification of fluxes below ice-covered area (using observed ice-cover35 or a sea-ice model), the addition of river runoffs as surface freshwater36 fluxes, and the addition of a freshwater flux adjustment on order to avoid a37 mean sea-level drift.8 \newpage 9 $\ $\newline % force a new ligne 10 %---------------------------------------namsbc-------------------------------------------------- 11 \namdisplay{namsbc} 12 %-------------------------------------------------------------------------------------------------------------- 13 $\ $\newline % force a new ligne 14 15 The ocean needs six fields as surface boundary condition: 16 \begin{itemize} 17 \item the two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$ 18 \item the incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$ 19 \item the surface freshwater budget $\left( {\text{EMP},\;\text{EMP}_S } \right)$ 20 \end{itemize} 21 22 Four different ways to provide those six fields to the ocean are available which 23 are controlled by namelist variables: an analytical formulation (\np{ln\_ana}=true), 24 a flux formulation (\np{ln\_flx}=true), a bulk formulae formulation (CORE 25 (\np{ln\_core}=true) or CLIO (\np{ln\_clio}=true) bulk formulae) and a coupled 26 formulation (exchanges with a atmospheric model via the OASIS coupler) 27 (\np{ln\_cpl}=true). The frequency at which the six fields have to be updated is 28 the \np{nf\_sbc} namelist parameter. 29 In addition, the resulting fields can be further modified using 30 several namelist options. These options control the addition of a surface restoring 31 term to observed SST and/or SSS (\np{ln\_ssr}=true), the modification of fluxes 32 below ice-covered areas (using observed ice-cover or a sea-ice model) 33 (\np{nn\_ice}=0,1, 2 or 3), the addition of river runoffs as surface freshwater 34 fluxes (\np{ln\_rnf}=true), the addition of a freshwater flux adjustment in 35 order to avoid a mean sea-level drift (\np{nn\_fwb}= 0, 1 or 2), and the 36 transformation of the solar radiation (if provided as daily mean) into a diurnal 37 cycle (\np{ln\_dm2dc}=true). 38 38 39 39 In this chapter, we first discuss where the surface boundary condition 40 40 appears in the model equations. Then we present the four ways of providing 41 the surface boundary condition. Finally, the different options that modify 42 the fluxes inside the ocean are discussed. 43 44 45 46 47 48 49 50 51 41 the surface boundary condition. Finally, the different options that further modify 42 the fluxes applied to the ocean are discussed. 52 43 53 44 … … 60 51 61 52 The surface ocean stress is the stress exerted by the wind and the sea-ice 62 on the ocean. Their two components are assumed to be interpolated on the 63 ocean mesh, i.e. provided at U- and V-points and projected onto the 64 (\textbf{i},\textbf{j}) referential. They are applied as a surface boundary 65 condition of the computation of the momentum vertical mixing trend 66 (\textbf{dynzdf} module) : 53 on the ocean. The two components of stress are assumed to be interpolated 54 onto the ocean mesh, $i.e.$ resolved onto the model (\textbf{i},\textbf{j}) direction 55 at $u$- and $v$-points They are applied as a surface boundary condition of the 56 computation of the momentum vertical mixing trend (\mdl{dynzdf} module) : 67 57 \begin{equation} \label{Eq_sbc_dynzdf} 68 58 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} … … 72 62 stress vector in the $(\textbf{i},\textbf{j})$ coordinate system. 73 63 74 The surface heat flux is decomposed in two parts, a non solar andsolar heat75 flux es. The former is the non penetrative part of the heat flux (i.e.76 sensible plus latent plus long wave heat fluxes). It is applied as a surface77 boundary condition trend of the first level temperature time evolution78 equation (\mdl{trasbc} module).64 The surface heat flux is decomposed into two parts, a non solar and a solar heat 65 flux, $Q_{ns}$ and $Q_{sr}$, respectively. The former is the non penetrative part 66 of the heat flux ($i.e.$ the sum of sensible, latent and long wave heat fluxes). 67 It is applied as a surface boundary condition trend of the first level temperature 68 time evolution equation (\mdl{trasbc} module). 79 69 \begin{equation} \label{Eq_sbc_trasbc_q} 80 70 \frac{\partial T}{\partial t}\equiv \cdots \;+\;\left. {\frac{Q_{ns} }{\rho 81 71 _o \;C_p \;e_{3T} }} \right|_{k=1} \quad 82 72 \end{equation} 83 84 The latter is the penetrative part of the heat flux. It is applied as a 3D 85 trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=T. 73 $Q_{sr}$ is the penetrative part of the heat flux. It is applied as a 3D 74 trends of the temperature equation (\mdl{traqsr} module) when \np{ln\_traqsr}=True. 86 75 87 76 \begin{equation} \label{Eq_sbc_traqsr} … … 89 78 \,e_{3T} }\delta _k \left[ {I_w } \right] 90 79 \end{equation} 91 92 where $I_w$ is an adimensional function that describes the way the light 80 where $I_w$ is a non-dimensional function that describes the way the light 93 81 penetrates inside the water column. It is generally a sum of decreasing 94 exponential (see \S\ref{TRA_qsr}).95 96 The surface freshwater budget is provided through two non-necessary97 identical fields EMP and EMP$_S $. Indeed, a surface freshwater98 flux has two effects: it changes the volume of the ocean and it changes the99 s urface concentration of salt (an others tracers). Therefore it appears in100 the sea surface height and salinity time evolution equations as a volume101 flux, EMP (\textit{dynspg\_xxx} modules), andconcentration/dilution effect,102 EMP$_{S}$ (\mdl{trasbc} module) , respectively.82 exponentials (see \S\ref{TRA_qsr}). 83 84 The surface freshwater budget is provided by fields: EMP and EMP$_S$ which 85 may or may not be identical. Indeed, a surface freshwater flux has two effects: 86 it changes the volume of the ocean and it changes the surface concentration of 87 salt (and other tracers). Therefore it appears in the sea surface height as a volume 88 flux, EMP (\textit{dynspg\_xxx} modules), and in the salinity time evolution equations 89 as a concentration/dilution effect, 90 EMP$_{S}$ (\mdl{trasbc} module). 103 91 \begin{equation} \label{Eq_trasbc_emp} 104 92 \begin{aligned} … … 109 97 \end{equation} 110 98 111 In the real ocean, EMP =EMP$_S$ and the ocean salt content is conserved,99 In the real ocean, EMP$=$EMP$_S$ and the ocean salt content is conserved, 112 100 but it exist several numerical reasons why this equality should be broken. 113 101 For example: 114 102 115 103 When rigid-lid assumption is made, the ocean volume becomes constant and 116 thus, EMP =0, not EMP$_{S }$.117 118 When a sea-ice model is considered, the water exchanged between ice and119 ocean is verylightly salty (mean sea-ice salinity is $\sim $\textit{4 psu}). In this case,104 thus, EMP$=$0, not EMP$_{S }$. 105 106 When the ocean is coupled to a sea-ice model, the water exchanged between ice and 107 ocean is slightly salty (mean sea-ice salinity is $\sim $\textit{4 psu}). In this case, 120 108 EMP$_{S}$ take into account both concentration/dilution effect associated with 121 freezing/melting together withsalt flux between ice and ocean, while EMP is109 freezing/melting and the salt flux between ice and ocean, while EMP is 122 110 only the volume flux. In addition, in the current version of \NEMO, the 123 111 sea-ice is assumed to be above the ocean. Freezing/melting does not change 124 the ocean volume (not impact on EMP) while it modifies the SSS125 \colorbox{yellow}{(see {\S} on LIM sea-ice model)}.126 127 Note that SST can also be modified by a freshwater flux. Precipitation s(in128 particular solid one) may have a temperature significantly different from112 the ocean volume (not impact on EMP) but it modifies the SSS. 113 %gm \colorbox{yellow}{(see {\S} on LIM sea-ice model)}. 114 115 Note that SST can also be modified by a freshwater flux. Precipitation (in 116 particular solid precipitation) may have a temperature significantly different from 129 117 the SST. Due to the lack of information about the temperature of 130 precipitation s, we assume it is equal to the SST. Therefore, no118 precipitation, we assume it is equal to the SST. Therefore, no 131 119 concentration/dilution term appears in the temperature equation. It has to 132 be emphasised that this absence does not mean that there is no theat flux133 associated with precipitation! An excess of precipitation will change the134 ocean heat content and is therefore associated with a heat flux (not120 be emphasised that this absence does not mean that there is no heat flux 121 associated with precipitation! Precipitation can change the ocean volume and thus the 122 ocean heat content. It is therefore associated with a heat flux (not yet 135 123 diagnosed in the model) \citep{Roullet2000}). 136 124 137 \colorbox{yellow}{Miss: } 138 139 A extensive description of all namsbc namelist (parameter that have to be 140 created!) 141 142 Especially the \np{nf\_sbc}, the \mdl{sbc\_oce} module (fluxes + mean sst sss ssu 143 ssv) i.e. information required by flux computation or sea-ice 144 145 \colorbox{red}{Add nqsr = 0 / 1 replace key{\_}traqsr} 146 147 \mdl{sbc\_oce} containt the definition in memory of the 7 fields (6+runoff), add 148 a word on runoff: included in surface bc or add as lateral obc{\ldots}. 149 150 Sbcmod manage the ``providing'' (fourniture) to the ocean the 7 fields 151 152 Fluxes update only each nf{\_}sbc time step (namsbc) explain relation 153 between nf{\_}sbc and nf{\_}ice, do we define nf{\_}blk??? ? only one 154 nf{\_}sbc 155 156 Explain here all the namlist namsbc variable{\ldots}. 157 158 \colorbox{yellow}{End Miss } 159 160 The ocean model provides the following variables averaged over nf{\_}sbc 161 time-step: 125 %\colorbox{yellow}{Miss: } 126 % 127 %A extensive description of all namsbc namelist (parameter that have to be 128 %created!) 129 % 130 %Especially the \np{nf\_sbc}, the \mdl{sbc\_oce} module (fluxes + mean sst sss ssu 131 %ssv) i.e. information required by flux computation or sea-ice 132 % 133 %\mdl{sbc\_oce} containt the definition in memory of the 7 fields (6+runoff), add 134 %a word on runoff: included in surface bc or add as lateral obc{\ldots}. 135 % 136 %Sbcmod manage the ``providing'' (fourniture) to the ocean the 7 fields 137 % 138 %Fluxes update only each nf{\_}sbc time step (namsbc) explain relation 139 %between nf{\_}sbc and nf{\_}ice, do we define nf{\_}blk??? ? only one 140 %nf{\_}sbc 141 % 142 %Explain here all the namlist namsbc variable{\ldots}. 143 % 144 %\colorbox{yellow}{End Miss } 145 146 The ocean model provides the surface currents, temperature and salinity 147 averaged over \np{nf\_sbc} time-step (\ref{Tab_ssm}).The computation of the 148 mean is done in \mdl{sbcmod} module. 162 149 163 150 %-------------------------------------------------TABLE--------------------------------------------------- 164 \begin{table}[ htbp] \label{Tab_ssm}151 \begin{table}[tb] \label{Tab_ssm} 165 152 \begin{center} 166 153 \begin{tabular}{|l|l|l|l|} 167 154 \hline 168 Variable desc iption & Computer name & Units & point \\ \hline169 i-component of the surface current & ssu\_ u& $m.s^{-1}$ & U \\ \hline155 Variable description & Model variable & Units & point \\ \hline 156 i-component of the surface current & ssu\_m & $m.s^{-1}$ & U \\ \hline 170 157 j-component of the surface current & ssv\_m & $m.s^{-1}$ & V \\ \hline 171 158 Sea surface temperature & sst\_m & \r{}$K$ & T \\ \hline 172 159 Sea surface salinty & sss\_m & $psu$ & T \\ \hline 173 160 \end{tabular} 161 \caption{Ocean variables provided by the ocean to the surface module (SBC). 162 The variable are averaged over nf{\_}sbc time step, $i.e.$ the frequency of 163 computation of surface fluxes.} 174 164 \end{center} 175 165 \end{table} 176 166 %-------------------------------------------------------------------------------------------------------------- 177 167 178 The mean computation is done in sbcmod ( 179 180 \colorbox{yellow}{Penser a} mettre dans le restant l'info nf{\_}sbc ET nf{\_}sbc*rdt de sorte de 181 reinitialiser la moyenne si on change la frequence ou le pdt 182 183 NB: creer cn{\_}sbc{\_}ice (cn{\_} = character in the namelist) with 3 184 cases: 185 186 = `noice' no specific call 187 188 = `iceif ` ``ice-if'' sea ice, i.e. read observed ice-cover and modified sbc 189 bellow those area. 190 191 = `lim' LIM sea-ice model is called which update the sbc fields in ice 192 covered area 193 194 ? modify the nsbc{\_}ice variable depending of this parameter (from --1, 0 195 to 1) 196 \colorbox{yellow}{End Penser a} 168 169 170 %\colorbox{yellow}{Penser a} mettre dans le restant l'info nf{\_}sbc ET nf{\_}sbc*rdt de sorte de reinitialiser la moyenne si on change la frequence ou le pdt 171 197 172 198 173 % ================================================================ … … 203 178 \label{SBC_ana} 204 179 205 %---------------------------------------namtau - namflx-------------------------------------------------- 206 \namdisplay{namtau} 207 \namdisplay{namflx} 180 %---------------------------------------namsbc_ana-------------------------------------------------- 181 \namdisplay{namsbc_ana} 208 182 %-------------------------------------------------------------------------------------------------------------- 209 183 210 184 211 The analytical formulation of the surface boundary condition is set by 212 default. In this case, all the 6 fluxes needed by the ocean are assumed to 213 be uniform in space. They take constant values given in the namlist 214 namsbc{\_}ana : \textit{utau0}, \textit{vtau0}, \textit{qns0}, \textit{qsr0}, \textit{emp0} and \textit{emps0}. while the runoff is set to zero. In addition, 215 the wind is allowed to reach its nominal value within a given number of time 216 step (\textit{ntau000}). 217 218 If a user wants to applied a different analytical forcing, \mdl{sbcana} 219 module is the very place to do that. As an example, one can have a look to 220 the \mdl{sbc\_ana\_gyre} routine which provides the analytical forcing of the 185 The analytical formulation of the surface boundary condition is the default scheme. 186 In this case, all the six fluxes needed by the ocean are assumed to 187 be uniform in space. They take constant values given in the namelist 188 namsbc{\_}ana by the variables \np{rn\_utau0}, \np{rn\_vtau0}, \np{rn\_qns0}, 189 \np{rn\_qsr0}, and \np{rn\_emp0} (EMP$=$EMP$_S$). The runoff is set to zero. 190 In addition, the wind is allowed to reach its nominal value within a given number 191 of time steps (\np{nn\_tau000}). 192 193 If a user wants to apply a different analytical forcing, the \mdl{sbcana} 194 module can be modified to use another scheme. As an example, 195 the \mdl{sbc\_ana\_gyre} routine provides the analytical forcing for the 221 196 GYRE configuration (see GYRE configuration manual, in preparation). 222 197 … … 228 203 {Flux formulation (\mdl{sbcflx} module) } 229 204 \label{SBC_flx} 230 231 In the flux formulation (\key{sbcflx} defined), the surface boundary 205 %------------------------------------------namsbc_flx---------------------------------------------------- 206 \namdisplay{namsbc_flx} 207 %------------------------------------------------------------------------------------------------------------- 208 209 In the flux formulation (\np{ln\_flx}=true), the surface boundary 232 210 condition fields are directly read from input files. The user has to define 233 211 in the namelist namsbc{\_}flx the name of the file, the name of the variable 234 read in the file, the time frequency at which it is given , and a logical235 setting whether a time interpolation to the model time step is asked are not212 read in the file, the time frequency at which it is given (in hours), and a logical 213 setting whether a time interpolation to the model time step is required 236 214 for this field). (fld\_i namelist structure). 237 215 238 \colorbox{yellow}{ Describe the information given? }239 240 \colorbox{yellow}{ Add an info about on-line interpolation or not ? at with which scale{\ldots} }241 242 243 216 \textbf{Caution}: when the frequency is set to --12, the data are monthly 244 values. The re are assumed to be climatological values, so time interpolation217 values. These are assumed to be climatological values, so time interpolation 245 218 between December the 15$^{th}$ and January the 15$^{th}$ is done using 246 record 12 and 1219 records 12 and 1 247 220 248 221 When higher frequency is set and time interpolation is demanded, the model 249 222 will try to read the last (first) record of previous (next) year in a file 250 having the same name but a suffix {\_}prev{\_}year ( next{\_}year) being251 added . These file must only contenta single record. If they don't exist,252 the will assume that the previous year last recordis equal to the first253 record of the previousyear, and similarly, that the first record of the223 having the same name but a suffix {\_}prev{\_}year ({\_}next{\_}year) being 224 added (e.g. "{\_}1989"). These files must only contain a single record. If they don't exist, 225 the model assumes that the last record of the previous year is equal to the first 226 record of the current year, and similarly, that the first record of the 254 227 next year is equal to the last record of the current year. This will cause 255 the forcing to remain constant over the first and last half fld\_frequ 256 hours. 228 the forcing to remain constant over the first and last half fld\_frequ hours. 257 229 258 230 Note that in general, a flux formulation is used in associated with a 259 damping term to observed SST and/or SSS. See \S\ref{SBC_ssr} for its231 restoring term to observed SST and/or SSS. See \S\ref{SBC_ssr} for its 260 232 specification. 261 233 … … 271 243 using bulk formulae and atmospheric fields and ocean (and ice) variables. 272 244 273 The atmospheric fields used depends on the bulk formulae used. Two of them 274 are available : the CORE and CLIO bulk formulea. The choice is made by 275 activating the CPP key \key{sbcblk\_core} or 276 \key{sbcblk\_clio}, respectively. 277 278 \colorbox{yellow}{Note : if a sea-ice model is used then blah blah blah{\ldots}} 279 280 CORE bulk formulea 281 282 The CORE bulk formulae have been developed by \citet{LargeYeager2004}. They 283 have been design to handle the CORE forcing, a mixture of NCEP reanalysis 284 and satellite data. They use an inertial dissipative method to compute the 285 turbulent transfer coefficients (momentum, sensible heat and evaporation) 286 from the 10 meter wind speed, air temperature and specific humidity). 245 The atmospheric fields used depend on the bulk formulae used. Two bulk formulations 246 are available : the CORE and CLIO bulk formulea. The choice is made by setting to true 247 one of the following namelist variable : \np{ln\_core} and \np{ln\_clio}. 248 249 Note : in forced mode, when a sea-ice model is used, a bulk formulation have to be used. 250 Therefore the two bulk formulea provided include the computation of the fluxes over both 251 an ocean and an ice surface. 252 253 % ------------------------------------------------------------------------------------------------------------- 254 % CORE Bulk formulea 255 % ------------------------------------------------------------------------------------------------------------- 256 \subsection [CORE Bulk formulea (\np{ln\_core}=true)] 257 {CORE Bulk formulea (\np{ln\_core}=true, \mdl{sbcblk\_core})} 258 \label{SBC_blk_core} 259 %------------------------------------------namsbc_core---------------------------------------------------- 260 \namdisplay{namsbc_core} 261 %------------------------------------------------------------------------------------------------------------- 262 263 The CORE bulk formulae have been developed by \citet{LargeYeager2004}. 264 They have been designed to handle the CORE forcing, a mixture of NCEP 265 reanalysis and satellite data. They use an inertial dissipative method to compute 266 the turbulent transfer coefficients (momentum, sensible heat and evaporation) 267 from the 10 metre wind speed, air temperature and specific humidity. 268 269 Note that substituting ERA40 to NCEP reanalysis fields 270 does not require changes in the bulk formulea themself. 287 271 288 272 The required 8 input fields are: … … 293 277 \begin{tabular}{|l|l|l|l|} 294 278 \hline 295 Variable desciption & Computer name & Units& point \\ \hline296 i-component of the 10m air velocity & utau & $m.s^{-1}$ & T or U\\ \hline297 j-component of the 10m air velocity & vtau & $m.s^{-1}$ & T or V\\ \hline279 Variable desciption & Model variable & Units & point \\ \hline 280 i-component of the 10m air velocity & utau & $m.s^{-1}$ & T \\ \hline 281 j-component of the 10m air velocity & vtau & $m.s^{-1}$ & T \\ \hline 298 282 10m air temperature & tair & \r{}$K$ & T \\ \hline 299 283 Specific humidity & humi & \% & T \\ \hline … … 307 291 %-------------------------------------------------------------------------------------------------------------- 308 292 309 Note that the air velocity can be provided at either tracer ocean point or310 velocity ocean point. 311 312 \colorbox{yellow}{Explain low resolution, better to provide it at U-V, high resolution better} 313 314 \colorbox{yellow}{at T-point{\ldots} Explain why, scheme?} 315 316 \colorbox{yellow}{Add a namelist parameter to provide a switch from U/V or T (or I??) point}317 318 \colorbox{yellow}{ for utau/vtau} 319 320 CLIO bulk formulea 321 322 The CLIO bulk formulae have beendeveloped several years ago for the323 Louvain-la-neuve coupled ice-ocean model (CLIO, Goosse et al. 1997). It is a324 simpler bulk formulae that assumed the stress to be known and computes the325 radiative fluxes from a climatological cloud cover.293 Note that the air velocity is provided at a tracer ocean point, not at a velocity ocean point ($u$- and $v$-points). It is simpler and faster (less fields to be read), but it is not the recommended method when the ocean grid 294 size is the same or larger than the one of the input atmospheric fields. 295 296 % ------------------------------------------------------------------------------------------------------------- 297 % CLIO Bulk formulea 298 % ------------------------------------------------------------------------------------------------------------- 299 \subsection [CLIO Bulk formulea (\np{ln\_clio}=true)] 300 {CLIO Bulk formulea (\np{ln\_clio}=true, \mdl{sbcblk\_clio})} 301 \label{SBC_blk_clio} 302 %------------------------------------------namsbc_clio---------------------------------------------------- 303 \namdisplay{namsbc_clio} 304 %------------------------------------------------------------------------------------------------------------- 305 306 The CLIO bulk formulae were developed several years ago for the 307 Louvain-la-neuve coupled ice-ocean model (CLIO, \cite{Goosse_al_JGR99}). 308 They are simpler bulk formulae. They assume the stress to be known and 309 compute the radiative fluxes from a climatological cloud cover. 326 310 327 311 The required 7 input fields are: … … 332 316 \begin{tabular}{|l|l|l|l|} 333 317 \hline 334 Variable desciption & Computer name & Units& point \\ \hline318 Variable desciption & Model variable & Units & point \\ \hline 335 319 i-component of the ocean stress & utau & $N.m^{-2}$ & U \\ \hline 336 320 j-component of the ocean stress & vtau & $N.m^{-2}$ & V \\ \hline 337 321 Wind speed module & vatm & $m.s^{-1}$ & T \\ \hline 338 322 10m air temperature & tair & \r{}$K$ & T \\ \hline 339 S ecific humidity & humi & \% & T \\ \hline323 Specific humidity & humi & \% & T \\ \hline 340 324 Cloud cover & & \% & T \\ \hline 341 325 Total precipitation (liquid + solid) & precip & $Kg.m^{-2}.s^{-1}$ & T \\ \hline … … 346 330 %-------------------------------------------------------------------------------------------------------------- 347 331 348 As for the flux formulation, the input data informationrequired by the332 As for the flux formulation, information about the input data required by the 349 333 model is provided in the namsbc\_blk\_core or namsbc\_blk\_clio 350 namelist (via the structure fld\_i). The same assumption is made about the 351 value of the first and last record in each file. 352 334 namelist (via the structure fld\_i). The first and last record assumption is also made 335 (see \S\ref{SBC_flx}) 353 336 354 337 % ================================================================ … … 358 341 {Coupled formulation (\mdl{sbccpl} module)} 359 342 \label{SBC_cpl} 343 %------------------------------------------namsbc_cpl---------------------------------------------------- 344 \namdisplay{namsbc_cpl} 345 %------------------------------------------------------------------------------------------------------------- 360 346 361 347 In the coupled formulation of the surface boundary condition, the fluxes are … … 364 350 the atmospheric component. 365 351 352 The generalised coupled interface is under development. It should be available 353 in summer 2008. It will include the ocean interface for most of the European 354 atmospheric GCM (ARPEGE, ECHAM, ECMWF, HadAM, LMDz). 355 366 356 367 357 % ================================================================ 368 358 % Miscellanea options 369 359 % ================================================================ 370 \section{Miscellane aoptions}360 \section{Miscellaneous options} 371 361 \label{SBC_misc} 372 362 … … 377 367 {Surface restoring to observed SST and/or SSS (\mdl{sbcssr})} 378 368 \label{SBC_ssr} 379 380 In forced mode using flux formulation (default option or \key{flx} defined), a 381 feedback term \emph{must} be added to the specified surface heat flux $Q_{ns}^o$: 369 %------------------------------------------namsbc_ssr---------------------------------------------------- 370 \namdisplay{namsbc_ssr} 371 %------------------------------------------------------------------------------------------------------------- 372 373 In forced mode using a flux formulation (default option or \key{flx} defined), a 374 feedback term \emph{must} be added to the surface heat flux $Q_{ns}^o$: 382 375 \begin{equation} \label{Eq_sbc_dmp_q} 383 376 Q_{ns} = Q_{ns}^o + \frac{dQ}{dT} \left( \left. T \right|_{k=1} - SST_{Obs} \right) … … 385 378 where SST is a sea surface temperature field (observed or climatological), $T$ is 386 379 the model surface layer temperature and $\frac{dQ}{dT}$ is a negative feedback 387 coefficient usually taken equal to $-40~W.m^{-2}.$\r{}K$^{-1}$. For a $50~m$ mixed-layer depth, 388 this value corresponds to a relaxation time scale of two months. This term 389 ensures that if $T$ perfectly fits SST then $Q$ is equal to $Q_o$. 390 391 In the fresh water budget, a feedback term can also be added: 380 coefficient usually taken equal to $-40~W/m^2/K$. For a $50~m$ 381 mixed-layer depth, this value corresponds to a relaxation time scale of two months. 382 This term ensures that if $T$ perfectly matches the supplied SST, then $Q$ is 383 equal to $Q_o$. 384 385 In the fresh water budget, a feedback term can also be added. Converted into an 386 equivalent freshwater flux, it takes the following expression : 392 387 393 388 \begin{equation} \label{Eq_sbc_dmp_emp} 394 EMP = EMP_o +\gamma_s^{-1} \left(S-SSS_{Obs}\right)\left|S\right. 389 EMP = EMP_o + \gamma_s^{-1} e_{3t} \frac{ \left(\left.S\right|_{k=1}-SSS_{Obs}\right)} 390 {\left.S\right|_{k=1}} 395 391 \end{equation} 396 392 397 where EMP$_{o }$ is a net surface fresh water flux (observed, climatological or 398 atmospheric model product), \textit{SSS}$_{Obs}$is a sea surface salinity (usually a time 399 interpolation of the monthly mean PHC climatology \citep{Steele2001}, $S$ is the model 400 surface layer salinity and $\gamma_s$ is a negative feedback coefficient 401 which is provided as a namelist parameter. Unlike heat flux, there is no 402 physical justification for the feedback term in (III.4.4) as the atmosphere 403 does not care about ocean surface salinity \citep{Madec1997}. The 404 SSS restoring term can only be view as a flux correction on freshwater 405 fluxes to reduce the uncertainties we have on the observed freshwater 406 budget. 393 where EMP$_{o }$ is a net surface fresh water flux (observed, climatological or an 394 atmospheric model product), \textit{SSS}$_{Obs}$ is a sea surface salinity (usually a time 395 interpolation of the monthly mean Polar Hydrographic Climatology \citep{Steele2001}), 396 $\left.S\right|_{k=1}$ is the model surface layer salinity and $\gamma_s$ is a negative 397 feedback coefficient which is provided as a namelist parameter. Unlike heat flux, there is no 398 physical justification for the feedback term in \ref{Eq_sbc_dmp_emp} as the atmosphere 399 does not care about ocean surface salinity \citep{Madec1997}. The SSS restoring 400 term should be viewed as a flux correction on freshwater fluxes to reduce the 401 uncertainties we have on the observed freshwater budget. 407 402 408 403 % ------------------------------------------------------------------------------------------------------------- … … 411 406 \subsection{Handling of ice-covered area} 412 407 \label{SBC_ice-cover} 413 The presence of sea-ice at the top of the ocean 414 strongly modify the surface fluxes 415 416 The presence at the sea surface of an ice cover area modified all the fluxes 417 transmitted to the ocean. There is two cases whereas a sea-ice model is used 418 or not. 419 420 Without sea ice model, the information of ice-cover / open ocean is read in 421 a file (either the directly the ice-cover or the observed SST from which 422 ice-cover is deduced using a criteria on freezing point temperature). 408 409 The presence at the sea surface of an ice covered area modifies all the fluxes 410 transmitted to the ocean. There are several way to handle sea-ice in the system depending on the value of the \np{nn{\_}ice} namelist parameter. 411 \begin{description} 412 \item[nn{\_}ice = 0] there will never be sea-ice in the computational domain. This is a typical namelist value used for tropical ocean domain. The surface fluxes are simply specified for an ice-free ocean. No specific things are done for sea-ice. 413 \item[nn{\_}ice = 1] sea-ice can exist in the computational domain, but no sea-ice model is used. An observed ice covered area is read in a file. Below this area, the SST is restored to the freezing point and the heat fluxes are set to $-4~W/m^2$ ($-2~W/m^2$) in the northern (southern) hemisphere. The associated modification of the freshwater fluxes are done in such a way that the change in buoyancy fluxes remains zero. This prevents deep convection to occur when trying to reach the freezing point (and so ice covered area condition) while the SSS is too large. This manner of managing sea-ice area, just by using si IF case, is usually referred as the \textit{ice-if} model. It can be found in the \mdl{sbcice{\_}if} module. 414 \item[nn{\_}ice = 2 or more] A full sea ice model is used. This model computes the ice-ocean fluxes, that are combined with the air-sea fluxes using the ice fraction of each model cell to provide the surface ocean fluxes. Note that the activation of a sea-ice model is is done by defining a CPP key (\key{lim2} or \key{lim3}). The activation automatically ovewrite the read value of nn{\_}ice to its appropriate value ($i.e.$ $2$ for LIM-2 and $3$ for LIM-3). 415 \end{description} 416 417 % {Description of Ice-ocean interface to be added here or in LIM 2 and 3 doc ?} 423 418 424 419 % ------------------------------------------------------------------------------------------------------------- … … 428 423 {Addition of river runoffs (\mdl{sbcrnf})} 429 424 \label{SBC_rnf} 425 %------------------------------------------namsbc_rnf---------------------------------------------------- 426 \namdisplay{namsbc_rnf} 427 %------------------------------------------------------------------------------------------------------------- 430 428 431 429 It is convenient to introduce the river runoff in the model as a surface 432 fresh water flux es. \colorbox{yellow}{{\ldots} blah blah{\ldots}.}430 fresh water flux. 433 431 434 432 \colorbox{yellow}{Nevertheless, Pb of vertical resolution and increase of Kz in vicinity of } … … 444 442 {Freshwater budget control (\mdl{sbcfwb})} 445 443 \label{SBC_fwb} 446 %--------------------------------------------namfwb-------------------------------------------------------- 447 \namdisplay{namfwb} 448 %-------------------------------------------------------------------------------------------------------------- 449 450 To be written latter... 444 445 To be written later... 451 446 452 447 \gmcomment{The descrition of the technique used to control the freshwater budget has to be added here} -
trunk/DOC/TexFiles/Chapters/Chap_TRA.tex
r817 r994 202 202 % ------------------------------------------------------------------------------------------------------------- 203 203 \subsection [$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4})] 204 204 {$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=.true.)} 205 205 \label{TRA_adv_cen4} 206 206 -
trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex
r817 r994 7 7 8 8 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. 9 \gmcomment{Steven remark : problem here with turbulent vs turbulence. I've changed "turbulent closure" to "turbulence closure", "turbulent mixing length" to "turbulence mixing length", but I've left "turbulent kinetic energy" alone - though I think it is an historical abberation! 10 Gurvan : I kept "turbulent closure"...} 11 \gmcomment{Steven bis : parameterization is the american spelling, parameterisation is the british} 12 9 13 10 14 % ================================================================ … … 14 18 \label{ZDF_zdf} 15 19 16 The discrete form of the ocean subgrid scale physics has been presented in \S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries, the turbulent fluxes of momentum, heat and salt have to be defined. At the surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}), while at the bottom they are set to zero for heat and salt, unless a geothermal flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} defined, see \S\ref{TRA_bbc}), and specified through a bottom friction parameterization for momentum (see \S\ref{ZDF_bfr}). 20 The discrete form of the ocean subgrid scale physics has been presented in 21 \S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries, 22 the turbulent fluxes of momentum, heat and salt have to be defined. At the 23 surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}), 24 while at the bottom they are set to zero for heat and salt, unless a geothermal 25 flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} 26 defined, see \S\ref{TRA_bbc}), and specified through a bottom friction 27 parameterization for momentum (see \S\ref{ZDF_bfr}). 17 28 18 29 In this section we briefly discuss the various choices offered to compute 19 the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ , $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$-points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These coefficients can be assumed to be either constant, or a function of the local Richardson number, or computed from a turbulent closure model (either TKE or KPP formulation). The computation of these coefficients is initialized in \mdl{zdfini} module and performed in \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer diffusion, including the surface forcing, 20 are computed and added to the general trend in \mdl{dynzdf} and \mdl{trazdf} modules, respectively. These trends can be computed using either a forward time scheme (cpp variable 21 \np{np\_zdfexp} or a backward time scheme (default option) depending on the magnitude of the mixing coefficients used, and thus of the formulation used (see \S\ref{DOM_nxt}). 30 the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ , 31 $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- 32 points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These 33 coefficients can be assumed to be either constant, or a function of the local 34 Richardson number, or computed from a turbulent closure model (either 35 TKE or KPP formulation). The computation of these coefficients is initialized 36 in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or 37 \mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer 38 diffusion, including the surface forcing, are computed and added to the 39 general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 40 These trends can be computed using either a forward time stepping scheme 41 (namelist parameter \np{np\_zdfexp}=true) or a backward time stepping 42 scheme (\np{np\_zdfexp}=false) depending on the magnitude of the mixing 43 coefficients, and thus of the formulation used (see \S\ref{DOM_nxt}). 22 44 23 45 % ------------------------------------------------------------------------------------------------------------- … … 30 52 %-------------------------------------------------------------------------------------------------------------- 31 53 32 When the \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 33 constant values over the whole ocean. This is the crudest way to define the 34 vertical ocean physics. It is recommended to use this option only in process 35 studies, not in basin scale simulation. Typical values used in this case 36 are: 54 When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients 55 are set to constant values over the whole ocean. This is the crudest way to define 56 the vertical ocean physics. It is recommended that this option is only used in 57 process studies, not in basin scale simulations. Typical values used in this case are: 37 58 \begin{align*} 38 59 A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1} \\ … … 41 62 \end{align*} 42 63 43 These values are set through \np{avm0} and \np{avt0} namelist parameters. In any case, do not use 44 values smaller that those associated to the molecular viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity. 64 These values are set through the \np{avm0} and \np{avt0} namelist parameters. 65 In all cases, do not use values smaller that those associated with the molecular 66 viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, 67 $\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity. 45 68 46 69 … … 55 78 %-------------------------------------------------------------------------------------------------------------- 56 79 57 When \key{zdfric} is defined, a local Richardson number dependent formulation of the vertical momentum and tracer eddy coefficients is set. The vertical mixing coefficients are diagnosed from the large scale variables computed by the model (order 0.5 closure scheme). \textit{In situ} measurements allow to link vertical turbulent activity to large scale ocean structures. The hypothesis of a mixing mainly maintained by the growth of Kelvin-Helmholtz like instabilities leads to a dependency between the vertical turbulent eddy coefficients and the local Richardson number ($i.e.$ ratio of stratification over vertical shear). Following \citet{PacPhil1981}, the following formulation has been implemented: 80 When \key{zdfric} is defined, a local Richardson number dependent formulation 81 for the vertical momentum and tracer eddy coefficients is set. The vertical mixing 82 coefficients are diagnosed from the large scale variables computed by the model. 83 \textit{In situ} measurements have been used to link vertical turbulent activity to 84 large scale ocean structures. The hypothesis of a mixing mainly maintained by the 85 growth of Kelvin-Helmholtz like instabilities leads to a dependency between the 86 vertical turbulence eddy coefficients and the local Richardson number ($i.e.$ the 87 ratio of stratification to vertical shear). Following \citet{PacPhil1981}, the following 88 formulation has been implemented: 58 89 \begin{equation} \label{Eq_zdfric} 59 90 \left\{ \begin{aligned} … … 63 94 \end{aligned} \right. 64 95 \end{equation} 65 where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number, $N$ is the local brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), $A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1} $ is the maximum value that can be reached by the coefficient when $Ri\leq 0$, $a=5$ and $n=2$. The last three coefficients can be modified by setting \np{avmri}, \np{alp} and \np{nric} namelist parameter, respectively. 96 where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson 97 number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 98 $A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the 99 constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ 100 is the maximum value that can be reached by the coefficient when $Ri\leq 0$, 101 $a=5$ and $n=2$. The last three values can be modified by setting the 102 \np{avmri}, \np{alp} and \np{nric} namelist parameters, respectively. 66 103 67 104 % ------------------------------------------------------------------------------------------------------------- … … 75 112 %-------------------------------------------------------------------------------------------------------------- 76 113 77 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on a prognostic equation for $\bar {e}$, the turbulent kinetic energy, and a closure assumption for the turbulent length scales. This turbulent closure model has been developed by \citet{Bougeault1989} in atmospheric cases, adapted by \citet{Gaspar1990} for oceanic cases, and embedded in OPA by \citet{Blanke1993} for equatorial Atlantic simulations. Since then, significant modifications have been introduced by \citet{Madec1998} in both the implementation and the formulation of the mixing length scale. The time evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical shear, its destruction through stratification, its vertical diffusion and its dissipation of \citet{Kolmogorov1942} type: 114 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE 115 turbulent closure model based on a prognostic equation for $\bar {e}$, the turbulent 116 kinetic energy, and a closure assumption for the turbulence length scales. This 117 turbulent closure model has been developed by \citet{Bougeault1989} in the 118 atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and 119 embedded in OPA by \citet{Blanke1993} for equatorial Atlantic simulations. Since 120 then, significant modifications have been introduced by \citet{Madec1998} in both 121 the implementation and the formulation of the mixing length scale. The time 122 evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical 123 shear, its destruction through stratification, its vertical diffusion, and its dissipation 124 of \citet{Kolmogorov1942} type: 78 125 \begin{equation} \label{Eq_zdftke_e} 79 126 \frac{\partial \bar{e}}{\partial t} = … … 87 134 \begin{equation} \label{Eq_zdftke_kz} 88 135 \begin{split} 89 A^{vm} &= C_k\ l_k\ \sqrt {\bar{e}} 90 \\ 136 A^{vm} &= C_k\ l_k\ \sqrt {\bar{e}} \\ 91 137 A^{vT} &= A^{vm} / P_{rt} 92 138 \end{split} 93 139 \end{equation} 94 where $N$ designates the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 95 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing turbulent 96 length scales, $P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and 97 $C_\epsilon = 0.1$ are designed to deal with vertical mixing at any depth \citep{Gaspar1990}. They are set through namelist parameter \np{ediff} and \np{ediss}. $P_{rt} $ can be 98 set to unity or, following \citet{Blanke1993}, be a function of the 99 local Richardson number, $R_i $: 140 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 141 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 142 $P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and 143 $C_\epsilon = 0.1$ are designed to deal with vertical mixing at any depth 144 \citep{Gaspar1990}. They are set through namelist parameters \np{ediff} 145 and \np{ediss}. $P_{rt} $ can be set to unity or, following \citet{Blanke1993}, 146 be a function of the local Richardson number, $R_i $: 100 147 \begin{align*} \label{Eq_prt} 101 148 P_{rt} = \begin{cases} … … 105 152 \end{cases} 106 153 \end{align*} 107 Note that a horizontal Shapiro filter can be optionally applied to $R_i$. Nevertheless it is an obsolescent option that is notrecommanded. The choice of $P_{rt} $ is controlled by \np{npdl} namelist parameter. 108 109 For computational efficiency, the original formulation of the turbulent length scales proposed by \citet{Gaspar1990} has been simplified. Four formulations are proposed, the choice of which is controlled by \np{nmxl} namelist parameter. The first two are based on the following first order approximation \citep{Blanke1993}: 154 Note that a horizontal Shapiro filter can optionally be applied to $R_i$. 155 However it is an obsolescent option that is not recommended. 156 The choice of $P_{rt} $ is controlled by the \np{npdl} namelist parameter. 157 158 For computational efficiency, the original formulation of the turbulence length 159 scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 160 are proposed, the choice of which is controlled by the \np{nmxl} namelist 161 parameter. The first two are based on the following first order approximation 162 \citep{Blanke1993}: 110 163 \begin{equation} \label{Eq_tke_mxl0_1} 111 164 l_k = l_\epsilon = \sqrt {2 \bar e} / N 112 165 \end{equation} 113 which is obtained in a stable stratified region with constant values of the brunt-Vais\"{a}l\"{a} frequency. The resulting turbulent length scale is bounded by the distance to the surface or to the bottom (\np{nmxl}=0) or by the local vertical scale factor (\np{nmxl}=1). \citet{Blanke1993} notice that this simplification has two major drawbacks: it has no sense for local unstable 114 stratification and the computation no longer uses the whole information contained in the vertical density profile. To overcome this drawbacks, \citet{Madec1998} introduces the \np{nmxl}=2 or 3 cases, which add of an hypothesis on the vertical gradient of the computed length scale. So, the length scales are first evaluated as in \eqref{Eq_tke_mxl0_1} and then bounded such that: 166 which is valid in a stable stratified region with constant values of the brunt- 167 Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance 168 to the surface or to the bottom (\np{nmxl}=0) or by the local vertical scale factor (\np{nmxl}=1). \citet{Blanke1993} notice that this simplification has two major 169 drawbacks: it makes no sense for locally unstable stratification and the 170 computation no longer uses all the information contained in the vertical density 171 profile. To overcome these drawbacks, \citet{Madec1998} introduces the 172 \np{nmxl}=2 or 3 cases, which add an extra assumption concerning the vertical 173 gradient of the computed length scale. So, the length scales are first evaluated 174 as in \eqref{Eq_tke_mxl0_1} and then bounded such that: 115 175 \begin{equation} \label{Eq_tke_mxl_constraint} 116 176 \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1 117 177 \qquad \text{with }\ l = l_k = l_\epsilon 118 178 \end{equation} 119 120 \eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length scale cannot be 121 larger than the variations of depth. It provides a better approximation of the \citet{Gaspar1990} formulation while being much less time consuming. In particular, it allows the length scale to be limited not only by the distance to the surface or to the ocean bottom but also by the distance to a 122 strongly stratified portion of the water column such as the thermocline (Fig.~\ref{Fig_mixing_length}). In order to imposed \eqref{Eq_tke_mxl_constraint} constraint, we introduce two additonnal length scal: $l_{up}$ and $l_{dwn}$, the upward and downward length scale, and evaluate the dissipation and mixing turbulent length scales as (caution here we use the numerical indexation): 179 \eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length 180 scale cannot be larger than the variations of depth. It provides a better 181 approximation of the \citet{Gaspar1990} formulation while being much less 182 time consuming. In particular, it allows the length scale to be limited not only 183 by the distance to the surface or to the ocean bottom but also by the distance 184 to a strongly stratified portion of the water column such as the thermocline 185 (Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint} 186 constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, 187 the upward and downward length scales, and evaluate the dissipation and 188 mixing turbulence length scales as (and note that here we use numerical 189 indexing): 123 190 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 124 191 \begin{figure}[!t] \label{Fig_mixing_length} \begin{center} … … 130 197 \begin{equation} \label{Eq_tke_mxl2} 131 198 \begin{aligned} 132 l_{up }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3T}^{(k)}\right)199 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3T}^{(k)}\ \ \ \; \right) 133 200 \quad &\text{ from $k=1$ to $jpk$ }\ \\ 134 201 l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3T}^{(k-1)} \right) … … 136 203 \end{aligned} 137 204 \end{equation} 138 where $l^{(k)}$ is compute using \eqref{Eq_tke_mxl0_1}, $i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$. 139 140 In the \np{nmxl}=2 case, the dissipation and mixing turbulent length scales take a same value: $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np{nmxl}=2 case, the dissipation and mixing turbulent length scales are give as in \citet{Gaspar1990}: 205 where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, 206 $i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$. 207 208 In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same 209 value: $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the 210 \np{nmxl}=2 case, the dissipation and mixing turbulence length scales are give 211 as in \citet{Gaspar1990}: 141 212 \begin{equation} \label{Eq_tke_mxl_gaspar} 142 213 \begin{aligned} … … 149 220 stress field: $\bar{e}=ebb\;\left| \tau \right|$ ($ebb=60$ by default) 150 221 with a minimal threshold of $emin0=10^{-4}~m^2.s^{-2}$ (namelist 151 parameters). Its bottom value is assumed to be equal to the value of the level just above. The time 152 integration of the $\bar{e}$ equation may formally lead to negative values 153 because the numerical scheme does not ensure the positivity. To overcome 154 this problem, a cut-off in the minimum value of $\bar{e}$ is used. Following 155 \citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations to match 156 \citet{Gargett1984} one for the diffusion in the thermocline and deep ocean 157 $(A^{vT} = 10^{-3} / N)$. In addition, a 158 cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical 222 parameters). Its value at the bottom of the ocean is assumed to be 223 equal to the value of the level just above. The time integration of the 224 $\bar{e}$ equation may formally lead to negative values because the 225 numerical scheme does not ensure its positivity. To overcome this 226 problem, a cut-off in the minimum value of $\bar{e}$ is used. Following 227 \citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 228 This allows the subsequent formulations to match that of\citet{Gargett1984} 229 for the diffusion in the thermocline and deep ocean : $(A^{vT} = 10^{-3} / N)$. 230 In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical 159 231 instabilities associated with too weak vertical diffusion. They must be 160 232 specified at least larger than the molecular values, and are set through … … 172 244 173 245 The KKP scheme has been implemented by J. Chanut ... 246 174 247 \colorbox{yellow}{Add a description of KPP here.} 175 248 … … 188 261 occur at particular ocean grid points. In nature, convective processes 189 262 quickly re-establish the static stability of the water column. These 190 processes have been removed from the model via the hydrostatic assumption:191 they must be parameterized. Three parameterisations are available to deal192 with convective processes: either a non-penetrative convective adjustment or193 an enhanced vertical diffusion, or/and the use of a turbulent closure194 scheme.263 processes have been removed from the model via the hydrostatic 264 assumption so they must be parameterized. Three parameterizations 265 are available to deal with convective processes: a non-penetrative 266 convective adjustment or an enhanced vertical diffusion, or/and the 267 use of a turbulent closure scheme. 195 268 196 269 % ------------------------------------------------------------------------------------------------------------- … … 207 280 208 281 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 209 \begin{figure}[!ht ] \label{Fig_npc} \begin{center}282 \begin{figure}[!htb] \label{Fig_npc} \begin{center} 210 283 \includegraphics[width=0.90\textwidth]{./Figures/Fig_npc.pdf} 211 \caption {Example of an unstable density profile treated by the non penetrative convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from the surface to the bottom. It is found to be unstable between levels 3 and 4. They are mixed. The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed. The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed. The $1^{st}$ step ends since the density profile is then stable below the level 3. $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step: levels 2 to 5 are mixed. The new density profile is checked. It is found stable: end of algorithm.} 284 \caption {Example of an unstable density profile treated by the non penetrative 285 convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from 286 the surface to the bottom. It is found to be unstable between levels 3 and 4. 287 They are mixed. The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 288 are mixed. The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are 289 mixed. The $1^{st}$ step ends since the density profile is then stable below 290 the level 3. $2^{nd}$ step: the new $\rho$ profile is checked following the same 291 procedure as in $1^{st}$ step: levels 2 to 5 are mixed. The new density profile 292 is checked. It is found stable: end of algorithm.} 212 293 \end{center} \end{figure} 213 294 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 214 295 215 The non-penetrative convective adjustment algorithm is used when \np{ln\_zdfnpc}=T. It is applied at each \np{nnpc1} time step and mixes downwards instantaneously the statically unstable portion of the water column, but only until the density structure becomes neutrally stable ($i.e.$ until the mixed portion of the water column has \textit{exactly} the density of the water just below) \citep{Madec1991a}. This algorithm is an iterative process used in the following way 216 (Fig. \ref{Fig_npc}): going from the top of the ocean towards the bottom, the first 217 instability is searched. Assume in the following that the instability is 218 located between levels $k$ and $k+1$. The two levels are vertically mixed, for 219 potential temperature and salinity, conserving the heat and salt contents of 220 the water column. The new density is then computed by a linear 221 approximation. If the new density profile is still unstable between levels 222 $k+1$ and $k+2$, levels $k$, $k+1$ and $k+2$ are then mixed. This process is repeated until 223 stability is established below the level $k$ (the mixing process can go down to 224 the ocean bottom). The algorithm is repeated to check if the density profile 296 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true. 297 It is applied at each \np{nnpc1} time step and mixes downwards instantaneously 298 the statically unstable portion of the water column, but only until the density 299 structure becomes neutrally stable ($i.e.$ until the mixed portion of the water 300 column has \textit{exactly} the density of the water just below) \citep{Madec1991a}. 301 The associated algorithm is an iterative process used in the following way 302 (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is 303 found. Assume in the following that the instability is located between levels 304 $k$ and $k+1$. The potential temperature and salinity in the two levels are 305 vertically mixed, conserving the heat and salt contents of the water column. 306 The new density is then computed by a linear approximation. If the new 307 density profile is still unstable between levels $k+1$ and $k+2$, levels $k$, 308 $k+1$ and $k+2$ are then mixed. This process is repeated until stability is 309 established below the level $k$ (the mixing process can go down to the 310 ocean bottom). The algorithm is repeated to check if the density profile 225 311 between level $k-1$ and $k$ is unstable and/or if there is no deeper instability. 226 312 227 This algorithm is significantly different from mixing two by two statically unstable levels. The latter procedure cannot converge with a finite number 228 of iterations for some vertical profiles while the algorithm used in OPA 229 converges for any profile in a number of iterations less than the number of 230 vertical levels. This property is of paramount importance as pointed out by 231 \citet{Killworth1989}: it avoids the existence of permanent and unrealistic 232 static instabilities at the sea surface. This non-penetrative convective 233 algorithm has been proved successful in studying the deep water formation in 234 the north-western Mediterranean Sea \citep{Madec1991a, Madec1991b, Madec1991c}. 235 236 Note that in this algorithm the potential density referenced to the sea 237 surface is used to check whether the density profile is stable or not. 238 Moreover, the mixing in potential density is assumed to be linear. This 239 assures the convergence of the algorithm even when the equation of state is 240 non-linear. Small static instabilities can thus persist due to cabbeling: 241 they will be treated at the next time step. Moreover, temperature and 242 salinity, and thus density, are mixed, but the corresponding velocity fields 243 remain unchanged. When using a Richardson dependent eddy viscosity, the 244 mixing of momentum is done through the vertical diffusion: after a static 245 adjustment, the Richardson number is zero and thus the eddy viscosity 246 coefficient is at a maximum. When this algorithm is used with constant 247 vertical eddy viscosity, spurious solution can occur as the vertical 248 momentum diffusion remains small even after a static adjustment. In that 249 latter case, we recommend to add momentum mixing in a manner that mimics the 250 mixing in temperature and salinity \citep{Speich1992, Speich1996}. 251 252 %AMT This presentation fails to mention the big drawback of this scheme: many water masses of the world ocean, especially AABW, are unstable when represented in surface-referenced potential density sigma_0. The scheme erroneously mixes them up. 313 This algorithm is significantly different from mixing statically unstable levels 314 two by two. The latter procedure cannot converge with a finite number 315 of iterations for some vertical profiles while the algorithm used in \NEMO 316 converges for any profile in a number of iterations which is less than the 317 number of vertical levels. This property is of paramount importance as 318 pointed out by \citet{Killworth1989}: it avoids the existence of permanent 319 and unrealistic static instabilities at the sea surface. This non-penetrative 320 convective algorithm has been proved successful in studies of the deep 321 water formation in the north-western Mediterranean Sea 322 \citep{Madec1991a, Madec1991b, Madec1991c}. 323 324 Note that in the current implementation of this algorithm presents several 325 limitations. First, potential density referenced to the sea surface is used to 326 check whether the density profile is stable or not. This is a strong 327 simplification which leads to large errors for realistic ocean simulations. 328 Indeed, many water masses of the world ocean, especially Antarctic Bottom 329 Water, are unstable when represented in surface-referenced potential density. 330 The scheme will erroneously mix them up. Second, the mixing of potential 331 density is assumed to be linear. This assures the convergence of the algorithm 332 even when the equation of state is non-linear. Small static instabilities can thus 333 persist due to cabbeling: they will be treated at the next time step. 334 Third, temperature and salinity, and thus density, are mixed, but the 335 corresponding velocity fields remain unchanged. When using a Richardson 336 Number dependent eddy viscosity, the mixing of momentum is done through 337 the vertical diffusion: after a static adjustment, the Richardson Number is zero 338 and thus the eddy viscosity coefficient is at a maximum. When this convective 339 adjustment algorithm is used with constant vertical eddy viscosity, spurious 340 solutions can occur since the vertical momentum diffusion remains small even 341 after a static adjustment. In that case, we recommend the addition of momentum 342 mixing in a manner that mimics the mixing in temperature and salinity 343 \citep{Speich1992, Speich1996}. 253 344 254 345 % ------------------------------------------------------------------------------------------------------------- … … 263 354 %-------------------------------------------------------------------------------------------------------------- 264 355 265 The enhanced vertical diffusion parameterization is used when \np{ln\_zdfevd} is 266 defined. In this case, the vertical eddy mixing coefficients are assigned to be very large (a typical value is $1\;m^2s^{-1})$ in regions where the stratification is unstable (i.e. when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar1997, Lazar1999}. This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and 267 tracers (\np{n\_evdm}=1) mixing coefficients. 268 269 In practice, when $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$ are set to a large value, 270 \np{avevd}, and if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm} $ . Typical value for $avevd$ is in-between 1 and $100~m^2.s^{-1}$. This parameterisation of convective processes is less time consuming than the convective adjustment algorithm presented above when mixing both tracers and momentum in case of static instabilities. It requires the use of an implicit time stepping on vertical diffusion terms (i.e. \np{ln\_zdfexp}=F). 356 The enhanced vertical diffusion parameterization is used when \np{ln\_zdfevd}=true. 357 In this case, the vertical eddy mixing coefficients are assigned very large values 358 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 359 ($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar1997, Lazar1999}. 360 This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and 361 tracers (\np{n\_evdm}=1). 362 363 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and 364 if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 365 values also, are set equal to the namelist parameter \np{avevd}. A typical value 366 for $avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterization of 367 convective processes is less time consuming than the convective adjustment 368 algorithm presented above when mixing both tracers and momentum in the 369 case of static instabilities. It requires the use of an implicit time stepping on 370 vertical diffusion terms (i.e. \np{ln\_zdfexp}=false). 271 371 272 372 % ------------------------------------------------------------------------------------------------------------- … … 276 376 \label{ZDF_tcs} 277 377 278 The turbulent closure scheme presented in \S\ref{ZDF_tke} and used when the 279 \key{zdftke} is defined allows, in theory, to deal with 280 statically unstable density profiles. In such a 281 case, the term of destruction of turbulent kinetic energy through 282 stratification in \eqref{Eq_zdftke_e} becomes a source term as $N^2$ is negative. It 283 results large values of both $A_T^{vT}$ and the four neighbouring$A_u^{vm} 284 {and}\;A_v^{vm} $ (up to $1\;m^2s^{-1})$ that are able to restore the 285 static stability of the water column in a way similar to that of the 286 enhanced vertical diffusion parameterization (\S\ref{ZDF_evd}). Nevertheless, the 287 eddy coefficients computed by the turbulent scheme do usually not exceed 288 $10^{-2}m.s^{-1}$ in the vicinity of the sea surface (first ocean layer) due to 289 the bound of the turbulent length scale by the distance to the sea surface 290 (see {\S}VI.7-c). It can thus be useful to combine the enhanced vertical 291 diffusion with the turbulent closure, i.e. defining \np{np\_zdfevd} and 292 \key{zdftke} CPP variables all together. 293 294 The KPP scheme includes enhanced vertical diffusion in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp}, therefore \np{np\_zdfevd} should not be used with the KPP scheme. %gm% + one word on non local flux with KPP scheme 295 378 The TKE turbulent closure scheme presented in \S\ref{ZDF_tke} and used 379 when the \key{zdftke} is defined, in theory solves the problem of statically 380 unstable density profiles. In such a case, the term corresponding to the 381 destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_e} 382 becomes a source term, since $N^2$ is negative. It results in large values of 383 $A_T^{vT}$ and $A_T^{vT}$, and also the four neighbouring 384 $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values 385 restore the static stability of the water column in a way similar to that of the 386 enhanced vertical diffusion parameterization (\S\ref{ZDF_evd}). However, 387 in the vicinity of the sea surface (first ocean layer), the eddy coefficients 388 computed by the turbulence scheme do not usually exceed $10^{-2}m.s^{-1}$, 389 because the mixing length scale is bounded by the distance to the sea surface 390 (see \S\ref{ZDF_tke}). It can thus be useful to combine the enhanced vertical 391 diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc} 392 namelist parameter to true and defining the \key{zdftke} CPP key all together. 393 394 The KPP turbulent closure scheme already includes enhanced vertical diffusion 395 in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ 396 found in \mdl{zdfkpp}, therefore \np{np\_zdfevd} should not be used with the KPP 397 scheme. %gm% + one word on non local flux with KPP scheme trakpp.F90 module... 296 398 297 399 % ================================================================ … … 306 408 %-------------------------------------------------------------------------------------------------------------- 307 409 308 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. The former condition leads to salt fingering and the latter to diffusive convection. Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the oceans. \citet{Merryfield1999} include a parameterization of such phenomena in a global ocean model and show that it leads to relatively minor changes in circulation but exerts significant regional influences on temperature and salinity. 410 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher 411 water, or vice versa. The former condition leads to salt fingering and the latter 412 to diffusive convection. Double-diffusive phenomena contribute to diapycnal 413 mixing in extensive regions of the ocean. \citet{Merryfield1999} include a 414 parameterization of such phenomena in a global ocean model and show that 415 it leads to relatively minor changes in circulation but exerts significant regional 416 influences on temperature and salinity. 309 417 310 418 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients … … 313 421 &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS} 314 422 \end{align*} 315 where subscript $f$ represents mixing by salt fingering, 316 $d$ by diffusive convection, and $o$ by processes other than 317 double diffusion. The rates of double-diffusive mixing depend on buoyancy ratio $R_\rho = \alpha \partial_z T / \partial_z S$, where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline contraction (see \S\ref{TRA_eos}. To represent mixing of $S$ and $T$ by salt fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981): 423 where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection, 424 and $o$ by processes other than double diffusion. The rates of double-diffusive mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 425 where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline 426 contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt 427 fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981): 318 428 \begin{align} \label{Eq_zdfddm_f} 319 429 A_f^{vS} &= \begin{cases} … … 321 431 0 &\text{otherwise} 322 432 \end{cases} 323 \\ 433 \\ \label{Eq_zdfddm_f_T} 324 434 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 325 435 \end{align} … … 328 438 \begin{figure}[!t] \label{Fig_zdfddm} \begin{center} 329 439 \includegraphics[width=0.99\textwidth]{./Figures/Fig_zdfddm.pdf} 330 \caption {From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vT}$ for temperature and salt in regions of diffusive convection. Heavy curves denote the Federov parameterization and thin curves the Kelley parameterization. The latter is not implemented in \NEMO } 440 \caption {From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ 441 and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy 442 curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves 443 $A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and 444 $A_d^{vS}$ for temperature and salt in regions of diffusive convection. Heavy 445 curves denote the Federov parameterization and thin curves the Kelley 446 parameterization. The latter is not implemented in \NEMO. } 331 447 \end{center} \end{figure} 332 448 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 333 449 334 The factor 0.7 in \eqref{Eq_zdfddm_f} reflects the measured ratio $\alpha F_T /\beta F_S \approx 0.7$ of buoyancy fluxes due to transport of heat and salt (e.g., McDougall and Taylor 1984). Following \citet{Merryfield1999}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 450 The factor 0.7 in \eqref{Eq_zdfddm_f_T} reflects the measured ratio 451 $\alpha F_T /\beta F_S \approx 0.7$ of buoyancy flux of heat to buoyancy 452 flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following \citet{Merryfield1999}, 453 we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 335 454 336 455 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by Federov (1988) is used: 337 \begin{align} \label{Eq_zdfddm_d}456 \begin{align} \label{Eq_zdfddm_d} 338 457 A_d^{vT} &= \begin{cases} 339 458 1.3635 \, \exp{\left( 4.6\, \exp{ \left[ -0.54\,( R_{\rho}^{-1} - 1 ) \right] } \right)} … … 341 460 0 &\text{otherwise} 342 461 \end{cases} 343 \\ 462 \\ \label{Eq_zdfddm_d_S} 344 463 A_d^{vS} &= \begin{cases} 345 464 A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) … … 351 470 \end{align} 352 471 353 The dependences of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d} on $R_\rho$ are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing $R_\rho$ at each 354 grid point and time step. This is done in \mdl{eosbn2} at the same time as $N^2$ is computed. This avoids duplication in the computation of $\alpha$ and $\beta$ (which is usually quite expensive). 355 472 The dependencies of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$ 473 are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing 474 $R_\rho$ at each grid point on every time step. This is done in \mdl{eosbn2} at the 475 same time as $N^2$ is computed. This avoids duplication in the computation of 476 $\alpha$ and $\beta$ (which is usually quite expensive). 356 477 357 478 % ================================================================ … … 366 487 %-------------------------------------------------------------------------------------------------------------- 367 488 368 Both surface momentum flux (wind stress) and the bottom momentum flux369 (bottom friction) enter the equations as a condition on the vertical489 Both the surface momentum flux (wind stress) and the bottom momentum 490 flux (bottom friction) enter the equations as a condition on the vertical 370 491 diffusive flux. For the bottom boundary layer, one has: 371 492 \begin{equation} \label{Eq_zdfbfr_flux} … … 376 497 1~m in the ocean). How $\textbf{F}_h$ influences the interior depends on the 377 498 vertical resolution of the model near the bottom relative to the Ekman layer 378 depth. For example, in order to obtain an Ekman layer depth $d = \sqrt{2\;A^{vm}} / f = 50$~m, one needs a vertical diffusion coefficient $A^{vm} = 0.125$~m$^2$s$^{-1}$ (for a Coriolis 379 frequency $f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient 380 $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. When the vertical mixing coefficient is this small, using a flux condition is equivalent to entering the 381 viscous forces (either wind stress or bottom friction) as a body force over 382 the depth of the top or bottom model layer. To illustrate this, consider the 383 equation for $u$ at $k$, the last ocean level: 499 depth. For example, in order to obtain an Ekman layer depth 500 $d = \sqrt{2\;A^{vm}} / f = 50$~m, one needs a vertical diffusion coefficient 501 $A^{vm} = 0.125$~m$^2$s$^{-1}$ (for a Coriolis frequency 502 $f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient 503 $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. 504 When the vertical mixing coefficient is this small, using a flux condition is 505 equivalent to entering the viscous forces (either wind stress or bottom friction) 506 as a body force over the depth of the top or bottom model layer. To illustrate 507 this, consider the equation for $u$ at $k$, the last ocean level: 384 508 \begin{equation} \label{Eq_zdfbfr_flux2} 385 509 \frac{\partial u \; (k)}{\partial t} = \frac{1}{e_{3u}} \left[ A^{vm} \; (k) \frac{U \; (k-1) - U \; (k)}{e_{3uw} \; (k-1)} - F_u \right] \approx - \frac{F_u}{e_{3u}} 386 510 \end{equation} 387 For example, if the bottom layer thickness is 200~m, the Ekman transport will be388 distributed over that depth. On the other hand, if the vertical resolution511 For example, if the bottom layer thickness is 200~m, the Ekman transport will 512 be distributed over that depth. On the other hand, if the vertical resolution 389 513 is high (1~m or less) and a turbulent closure model is used, the turbulent 390 514 Ekman layer will be represented explicitly by the model. However, the 391 515 logarithmic layer is never represented in current primitive equation model 392 applications: it is \emph{necessary} to parameterize the flux $\textbf{F}_h $. Two393 choices are available in OPA: a linear and a quadratic bottom friction. Note394 that in both cases, the rotation between the interior velocity and the395 bottom friction is neglected in the present release of OPA.516 applications: it is \emph{necessary} to parameterize the flux $\textbf{F}_h $. 517 Two choices are available in \NEMO: a linear and a quadratic bottom friction. 518 Note that in both cases, the rotation between the interior velocity and the 519 bottom friction is neglected in the present release of \NEMO. 396 520 397 521 % ------------------------------------------------------------------------------------------------------------- 398 522 % Linear Bottom Friction 399 523 % ------------------------------------------------------------------------------------------------------------- 400 \subsection{Linear Bottom Friction }524 \subsection{Linear Bottom Friction (\np{nbotfr} = 1) } 401 525 \label{ZDF_bfr_linear} 402 526 403 %-------------------------------------------nambfr-------------------------------------------------------- 404 \begin{flushright} 405 (\textbf{namelist} !nbotfr :\textit{ nbotfr = 0, = 1 or = 3}) 406 \end{flushright} 407 408 The linear bottom friction parameterization assumes that the bottom friction is proportional to the interior velocity (i.e. the velocity of the last model level): 527 The linear bottom friction parameterization assumes that the bottom friction 528 is proportional to the interior velocity (i.e. the velocity of the last model level): 409 529 \begin{equation} \label{Eq_zdfbfr_linear} 410 \textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \textbf{U}_h^b 411 \end{equation} 412 where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom 413 ocean layer and $r$ a friction coefficient expressed in m.s$^{-1}$. This 414 coefficient is generally estimated by setting a typical decay time $\tau $ in the 415 deep ocean, $r = H / \tau$. Commonly accepted values of $\tau$ are of the 416 order of 100 to 200 days \citep{Weatherly1984}. A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ 417 corresponding to 115 days is usually used in quasi-geostrophic models. One may 418 consider the linear friction as an approximation of quadratic friction, 419 $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982}, Eq. 9.6.6). With a drag coefficient $C_D = 0.002$, a typical value of tidal currents $U_{av} =0.1$~m.s$^{-1}$, 420 and assuming an ocean depth $H = 4000$~m, the resulting friction 421 coefficient is $r = 4\;10^{-4}$~m.s$^{-1}$. This is the default value used in 422 OPA. It corresponds to a decay time scale of 115~days. It can be changed by 423 specifying \np{bfric1} (namelist parameter). 424 425 In the code, the bottom friction is specified by updating the value of the 426 vertical eddy coefficient at the bottom level. Indeed, the discrete 427 formulation of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that 428 $\textbf {U}_h =0$ inside the bottom, leads to 530 \textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 531 \end{equation} 532 where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom ocean 533 layer and $r$ is a friction coefficient expressed in m.s$^{-1}$. This coefficient 534 is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 535 and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted 536 values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly1984}. 537 A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used 538 in quasi-geostrophic models. One may consider the linear friction as an 539 approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982}, 540 Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed 541 of tidal currents of $U_{av} =0.1$~m.s$^{-1}$, and assuming an ocean depth 542 $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m.s$^{-1}$. 543 This is the default value used in \NEMO. It corresponds to a decay time scale 544 of 115~days. It can be changed by specifying \np{bfric1} (namelist parameter). 545 546 In the code, the bottom friction is imposed by updating the value of the 547 vertical eddy coefficient at the bottom level. Indeed, the discrete formulation 548 of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that 549 $\textbf {U}_h =0$ below the ocean floor, leads to 429 550 \begin{equation} \label{Eq_zdfbfr_linKz} 430 551 \begin{split} 431 552 A_u^{vm} &= r\;e_{3uw}\\ 432 A_v^{vm} &= r\;e_{3 uw}\\553 A_v^{vm} &= r\;e_{3vw}\\ 433 554 \end{split} 434 555 \end{equation} 435 556 436 Such an update is done in \mdl{zdfbfr} when \np{nbotfr}=1 and the value of $r$ used is 437 \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to set $r=0$ and leads to a free-slip bottom boundary condition, 438 while setting \np{nbotfr}=0 imposes $r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the 439 background vertical eddy coefficient: a no-slip boundary condition is used. 557 This update is done in \mdl{zdfbfr} when \np{nbotfr}=1. The value of $r$ 558 used is \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to setting $r=0$ and 559 leads to a free-slip bottom boundary condition. Setting \np{nbotfr}=0 sets 560 $r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the background 561 vertical eddy coefficient, and a no-slip boundary condition is imposed. 440 562 Note that this latter choice generally leads to an underestimation of the 441 bottom friction: for a deepest level thickness of $200~m$ and $A_{vb}^{\rm {\bf U}}442 =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient is only $r=10^{-6}$m.s$^{-1}$. 443 563 bottom friction: for example with a deepest level thickness of $200~m$ 564 and $A_{vb}^{\rm {\bf U}} =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient 565 is only $r=10^{-6}$m.s$^{-1}$. 444 566 445 567 % ------------------------------------------------------------------------------------------------------------- 446 568 % Non-Linear Bottom Friction 447 569 % ------------------------------------------------------------------------------------------------------------- 448 \subsection{Non-Linear Bottom Friction }570 \subsection{Non-Linear Bottom Friction (\np{nbotfr} = 2)} 449 571 \label{ZDF_bfr_nonlinear} 450 451 \begin{center}452 (\textbf{namelist} !nbotfr : \textit{nbotfr = 2})453 \end{center}454 572 455 573 The non-linear bottom friction parameterization assumes that the bottom … … 460 578 \end{equation} 461 579 462 with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity (i.e. the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient, and $e_b $ a bottom turbulent kinetic energy due to tides, internal waves breaking and other short time scale currents. A typical value of the drag coefficient is $C_D = 10^{-3} $. As an example, the CME experiment \citep{Treguier1992} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$, while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ 463 and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been set as 464 default value (\np{bfric2} and \np{bfeb2} namelist parameters). 465 466 As for the linear case, the bottom friction is specified in the code by 580 with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity ($i.e.$ 581 the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient, 582 and $e_b $ a bottom turbulent kinetic energy due to tides, internal waves 583 breaking and other short time scale currents. A typical value of the drag 584 coefficient is $C_D = 10^{-3} $. As an example, the CME experiment 585 \citep{Treguier1992} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$, 586 while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ 587 and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been 588 set as default values (\np{bfric2} and \np{bfeb2} namelist parameters). 589 590 As for the linear case, the bottom friction is imposed in the code by 467 591 updating the value of the vertical eddy coefficient at the bottom level: 468 469 592 \begin{equation} \label{Eq_zdfbfr_nonlinKz} 470 593 \begin{split} 471 594 A_u^{vm} &=C_D\; e_{3uw} \left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^ 472 595 {1/2}\\ 473 A_v^{vm} &=C_D\; e_{3uw} \left[ \left(\bar{\bar{u}}^{i +1,j}\right)^2 + v^2 + e_b \right]^{1/2}\\596 A_v^{vm} &=C_D\; e_{3uw} \left[ \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ 474 597 \end{split} 475 598 \end{equation} 476 599 477 600 This update is done in \mdl{zdfbfr}. The coefficients that control the strength of the 478 non-linear bottom friction are initialized as namelist parameters: ($C_D$= \np{bfri2}, and $e_b$ =\np{bfeb2}).479 480 % ================================================================ 601 non-linear bottom friction are initialized as namelist parameters: $C_D$= \np{bfri2}, and $e_b$ =\np{bfeb2}. 602 603 % ================================================================
Note: See TracChangeset
for help on using the changeset viewer.
