Changeset 10502 for NEMO/trunk/doc
- Timestamp:
- 2019-01-10T18:45:21+01:00 (6 years ago)
- Location:
- NEMO/trunk/doc/latex/NEMO/subfiles
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex
r10442 r10502 3 3 \begin{document} 4 4 % ================================================================ 5 % Chapter 2 Space and Time Domain (DOM)5 % Chapter 2 ——— Space and Time Domain (DOM) 6 6 % ================================================================ 7 7 \chapter{Space Domain (DOM)} … … 40 40 \begin{figure}[!tb] 41 41 \begin{center} 42 \includegraphics[ width=0.90\textwidth]{Fig_cell}42 \includegraphics[]{Fig_cell} 43 43 \caption{ 44 44 \protect\label{fig:cell} … … 46 46 $t$ indicates scalar points where temperature, salinity, density, pressure and 47 47 horizontal divergence are defined. 48 ($u$,$v$,$w$) indicates vector points,49 and $f$ indicates vorticity points where both relative and planetary vorticities are defined48 $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and 49 planetary vorticities are defined. 50 50 } 51 51 \end{center} … … 64 64 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 65 65 66 The ocean mesh (\ie the position of all the scalar and vector points) is defined by 67 the transformation that gives ($\lambda$ ,$\varphi$ ,$z$)as a function of $(i,j,k)$.66 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 67 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 68 68 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}. 69 69 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 70 70 the grid-point where the scale factors are defined. 71 71 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. 72 As a result, 73 the mesh on which partial derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, 74 and $\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity. 75 Discrete partial derivatives are formulated by the traditional, 76 centred second order finite difference approximation while 77 the scale factors are chosen equal to their local analytical value. 72 As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 73 $\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity. 74 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 75 while the scale factors are chosen equal to their local analytical value. 78 76 An important point here is that the partial derivative of the scale factors must be evaluated by 79 77 centred finite difference approximation, not from their analytical expression. 80 This preserves the symmetry of the discrete set of equations and 81 the refore satisfies many of thecontinuous properties (see \autoref{apdx:C}).78 This preserves the symmetry of the discrete set of equations and therefore satisfies many of 79 the continuous properties (see \autoref{apdx:C}). 82 80 A similar, related remark can be made about the domain size: 83 81 when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors 84 (see \autoref{eq:DOM_bar} )in the next section).82 (see \autoref{eq:DOM_bar} in the next section). 85 83 86 84 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 89 87 \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 90 88 \hline 91 T &$i$ & $j$ & $k$ \\ \hline 92 u & $i+1/2$ & $j$ & $k$ \\ \hline 93 v & $i$ & $j+1/2$ & $k$ \\ \hline 94 w & $i$ & $j$ & $k+1/2$ \\ \hline 95 f & $i+1/2$ & $j+1/2$ & $k$ \\ \hline 96 uw & $i+1/2$ & $j$ & $k+1/2$ \\ \hline 97 vw & $i$ & $j+1/2$ & $k+1/2$ \\ \hline 98 fw & $i+1/2$ & $j+1/2$ & $k+1/2$ \\ \hline 89 T & $i $ & $j $ & $k $ \\ 90 \hline 91 u & $i + 1/2$ & $j $ & $k $ \\ 92 \hline 93 v & $i $ & $j + 1/2$ & $k $ \\ 94 \hline 95 w & $i $ & $j $ & $k + 1/2$ \\ 96 \hline 97 f & $i + 1/2$ & $j + 1/2$ & $k $ \\ 98 \hline 99 uw & $i + 1/2$ & $j $ & $k + 1/2$ \\ 100 \hline 101 vw & $i $ & $j + 1/2$ & $k + 1/2$ \\ 102 \hline 103 fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 104 \hline 99 105 \end{tabular} 100 106 \caption{ 101 107 \protect\label{tab:cell} 102 108 Location of grid-points as a function of integer or integer and a half value of the column, line or level. 103 This indexing is only used for the writing of the semi -discrete equation.109 This indexing is only used for the writing of the semi -discrete equation. 104 110 In the code, the indexing uses integer values only and has a reverse direction in the vertical 105 111 (see \autoref{subsec:DOM_Num_Index}) … … 115 121 \label{subsec:DOM_operators} 116 122 117 Given the values of a variable $q$ at adjacent points, 118 the differencing and averaging operators at themidpoint between them are:119 \ [123 Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at 124 the midpoint between them are: 125 \begin{alignat*}{2} 120 126 % \label{eq:di_mi} 121 \begin{split} 122 \delta_i [q] &= \ \ q(i+1/2) - q(i-1/2) \\ 123 \overline q^{\,i} &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2 124 \end{split} 125 \] 126 127 Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and $k+1/2$. 127 \delta_i [q] &= & &q (i + 1/2) - q (i - 1/2) \\ 128 \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2 129 \end{alignat*} 130 131 Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 128 132 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at 129 133 a $t$-point has its three components defined at $u$-, $v$- and $w$-points while 130 its Laplaci en is defined at $t$-point.131 These operators have the following discrete forms in the curvilinear $s$-coordinate system:134 its Laplacian is defined at $t$-point. 135 These operators have the following discrete forms in the curvilinear $s$-coordinates system: 132 136 \[ 133 137 % \label{eq:DOM_grad} 134 \nabla q \equiv \frac{1}{e_{1u} } \delta_{i+1/2 } [q] \;\,\mathbf{i}135 + \frac{1}{e_{2v} } \delta_{j+1/2 } [q] \;\,\mathbf{j}136 + \frac{1}{e_{3w}} \delta_{k+1/2} [q] \;\,\mathbf{k}138 \nabla q \equiv \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i 139 + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j 140 + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k 137 141 \] 138 142 \begin{multline*} 139 143 % \label{eq:DOM_lap} 140 \Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 141 \;\left( \delta_i \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right] 142 + \delta_j \left[ \frac{e_{1v}\,e_{3v}} {e_{2v}} \;\delta_{j+1/2} [q] \right] \; \right) \\ 143 +\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} } \;\delta_{k+1/2} [q] \right] 144 \Delta q \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 145 \; \lt[ \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) 146 + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\ 147 + \frac{1}{e_{3t}} 148 \delta_k \lt[ \frac{1 }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] 144 149 \end{multline*} 145 150 146 Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ 147 defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, 148 and its divergence defined at $t$-points: 149 \begin{align*} 150 % \label{eq:DOM_curl} 151 \nabla \times {\rm{\bf A}}\equiv & 152 \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right) &\ \mathbf{i} \\ 153 +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1 \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right) &\ \mathbf{j} \\ 154 +& \frac{1}{e_{1f} \,e_{2f} } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2 \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right) &\ \mathbf{k} 155 \end{align*} 156 \begin{align*} 157 % \label{eq:DOM_div} 158 \nabla \cdot \rm{\bf A} \equiv 159 \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] 160 +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] 161 \end{align*} 151 Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 152 vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and 153 its divergence defined at $t$-points: 154 \begin{multline} 155 % \label{eq:DOM_curl} 156 \nabla \times \vect A \equiv \frac{1}{e_{2v} \, e_{3vw}} 157 \Big[ \delta_{j + 1/2} (e_{3w} \, a_3) 158 - \delta_{k + 1/2} (e_{2v} \, a_2) \Big] \vect i \\ 159 + \frac{1}{e_{2u} \, e_{3uw}} 160 \Big[ \delta_{k + 1/2} (e_{1u} \, a_1) 161 - \delta_{i + 1/2} (e_{3w} \, a_3) \Big] \vect j \\ 162 + \frac{1}{e_{1f} \, e_{2f}} 163 \Big[ \delta_{i + 1/2} (e_{2v} \, a_2) 164 - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k 165 \end{multline} 166 \begin{equation} 167 % \label{eq:DOM_div} 168 \nabla \cdot \vect A \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 169 \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] 170 + \frac{1}{e_{3t}} \delta_k (a_3) 171 \end{equation} 162 172 163 173 The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which 164 is a masked field ( \ieequal to zero inside solid area):174 is a masked field (i.e. equal to zero inside solid area): 165 175 \begin{equation} 166 176 \label{eq:DOM_bar} 167 \bar q = \frac{1}{H} \int_{k^b}^{k^o} {q\;e_{3q} \,dk} 168 \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } 177 \bar q = \frac{1}{H} \int_{k^b}^{k^o} q \; e_{3q} \, dk \equiv \frac{1}{H_q} \sum \limits_k q \; e_{3q} 169 178 \end{equation} 170 179 where $H_q$ is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 171 $k^b$ and $k^o$ are the bottom and surface $k$-indices, 172 and the symbol $k^o$ refers to a summation over all grid points of the same type in the direction indicated by 173 the subscript (here $k$). 180 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over 181 all grid points of the same type in the direction indicated by the subscript (here $k$). 174 182 175 183 In continuous form, the following properties are satisfied: 176 \begin{ equation}184 \begin{gather} 177 185 \label{eq:DOM_curl_grad} 178 \nabla \times \nabla q ={\rm {\bf {0}}} 179 \end{equation} 180 \begin{equation} 186 \nabla \times \nabla q = \vect 0 \\ 181 187 \label{eq:DOM_div_curl} 182 \nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0183 \end{ equation}188 \nabla \cdot (\nabla \times \vect A) = 0 189 \end{gather} 184 190 185 191 It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as 186 the scalar $q$ is taken at $t$-points and 187 the vector \textbf{A} has its components defined atvector points $(u,v,w)$.192 the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 193 vector points $(u,v,w)$. 188 194 189 195 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area. 190 Using integration by parts it can be shown that 191 the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, 192 and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$, $\overline{\,\cdot\,}^{\,k}$ and 193 $\overline{\,\cdot\,}^{\,k}$) are symmetric linear operators, \ie 194 \begin{align} 196 Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 197 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 198 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 199 \begin{alignat}{4} 195 200 \label{eq:DOM_di_adj} 196 \sum\limits_i { a_i \;\delta_i \left[ b \right]} 197 &\equiv -\sum\limits_i {\delta_{i+1/2} \left[ a \right]\;b_{i+1/2} } \\ 201 &\sum \limits_i a_i \; \delta_i [b] &\equiv &- &&\sum \limits_i \delta _{ i + 1/2} [a] &b_{i + 1/2} \\ 198 202 \label{eq:DOM_mi_adj} 199 \sum\limits_i { a_i \;\overline b^{\,i}} 200 & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } 201 \end{align} 202 203 In other words, the adjoint of the differencing and averaging operators are $\delta_i^*=\delta_{i+1/2}$ and 204 ${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively. 203 &\sum \limits_i a_i \; \overline b^{\, i} &\equiv & &&\sum \limits_i \overline a ^{\, i + 1/2} &b_{i + 1/2} 204 \end{alignat} 205 206 In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 207 $(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 205 208 These two properties will be used extensively in the \autoref{apdx:C} to 206 209 demonstrate integral conservative properties of the discrete formulation chosen. … … 215 218 \begin{figure}[!tb] 216 219 \begin{center} 217 \includegraphics[ width=0.90\textwidth]{Fig_index_hor}220 \includegraphics[]{Fig_index_hor} 218 221 \caption{ 219 222 \protect\label{fig:index_hor} … … 230 233 Therefore a specific integer indexing must be defined for points other than $t$-points 231 234 (\ie velocity and vorticity grid-points). 232 Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k =1$.235 Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$. 233 236 234 237 % ----------------------------------- … … 250 253 \label{subsec:DOM_Num_Index_vertical} 251 254 252 In the vertical, the chosen indexing requires special attention since 253 the $k$-axis is re-orientated downward in the \fortran code compared to254 the indexing used in the semi-discrete equations andgiven in \autoref{subsec:DOM_cell}.255 The sea surface corresponds to the $w$-level $k =1$ which is the same index as $t$-level just below255 In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in 256 the \fortran code compared to the indexing used in the semi -discrete equations and 257 given in \autoref{subsec:DOM_cell}. 258 The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below 256 259 (\autoref{fig:index_vert}). 257 The last $w$-level ($k =jpk$) either corresponds to the ocean floor or is inside the bathymetry while260 The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while 258 261 the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}). 259 262 Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, … … 262 265 have the same $i$ or $j$ index 263 266 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). 264 Since the scale factors are chosen to be strictly positive, a \emph{minus sign} appears in the \fortran 265 code \emph{before all the vertical derivatives} of the discrete equations given in this documentation. 267 Since the scale factors are chosen to be strictly positive, 268 a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of 269 the discrete equations given in this documentation. 266 270 267 271 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 268 272 \begin{figure}[!pt] 269 273 \begin{center} 270 \includegraphics[ width=.90\textwidth]{Fig_index_vert}274 \includegraphics[]{Fig_index_vert} 271 275 \caption{ 272 276 \protect\label{fig:index_vert} … … 287 291 The total size of the computational domain is set by the parameters \np{jpiglo}, 288 292 \np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively. 289 %%%290 %%%291 %%%292 293 Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when 293 294 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, … … 299 300 \section{Needed fields} 300 301 \label{sec:DOM_fields} 301 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 302 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 303 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 302 304 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 303 305 The associated scale factors are defined using the analytical first derivative of the transformation 304 306 \autoref{eq:scale_factors}. 305 Necessary fields for configuration definition are: \\ 306 Geographic position : 307 308 longitude: glamt, glamu, glamv and glamf (at T, U, V and F point) 309 310 latitude: gphit, gphiu, gphiv and gphif (at T, U, V and F point)\\ 311 Coriolis parameter (if domain not on the sphere): 312 313 ff\_f and ff\_t (at T and F point)\\ 314 Scale factors : 307 Necessary fields for configuration definition are: 308 309 \begin{itemize} 310 \item 311 Geographic position: 312 longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and 313 latitude with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif} 314 (all respectively at T, U, V and F point) 315 \item 316 Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t} 317 (at T and F point) 318 \item 319 Scale factors: 320 \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction), 321 \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and 322 \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\ 323 \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits), 324 \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed. 325 \end{itemize} 315 326 316 e1t, e1u, e1v and e1f (on i direction), 317 318 e2t, e2u, e2v and e2f (on j direction) and 319 320 ie1e2u\_v, e1e2u , e1e2v 321 322 e1e2u , e1e2v are u and v surfaces (if gridsize reduction in some straits)\\ 323 ie1e2u\_v is a flag to flag set u and v surfaces are neither read nor computed.\\ 324 325 These fields can be read in an domain input file which name is setted in 326 \np{cn\_domcfg} parameter specified in \ngn{namcfg}. 327 These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in 328 \ngn{namcfg}. 327 329 328 330 \nlst{namcfg} 329 or they can be defined in an analytical way in MY\_SRC directory of the configuration. 331 332 Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration. 330 333 For Reference Configurations of NEMO input domain files are supplied by NEMO System Team. 331 For analytical definition of input fields two routines are supplied: \mdl{us erdef\_hgr} and \mdl{userdef\_zgr}.332 They are an example of GYRE configuration parameters, and they are available in NEMO/OPA\_SRC/USRdirectory,333 they provide the horizontal and vertical mesh. 334 For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 335 They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory, 336 they provide the horizontal and vertical mesh. 334 337 % ------------------------------------------------------------------------------------------------------------- 335 338 % Needed fields … … 366 369 ($i$ and $j$, respectively) (geographical configuration of the mesh), 367 370 the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$, 368 and their derivatives $\lambda'(i) $ $\varphi'(j)$ in the \mdl{domhgr} module.371 and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module. 369 372 The model computes the grid-point positions and scale factors in the horizontal plane as follows: 370 \begin{flalign*} 371 \lambda_t &\equiv \text{glamt}= \lambda(i) & \varphi_t &\equiv \text{gphit} = \varphi(j)\\ 372 \lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\ 373 \lambda_v &\equiv \text{glamv}= \lambda(i) & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\ 374 \lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2) 375 \end{flalign*} 376 \begin{flalign*} 377 e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i) \; \cos\varphi(j) |& 378 e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)| \\ 379 e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2) \; \cos\varphi(j) |& 380 e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\ 381 e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i) \; \cos\varphi(j+1/2) |& 382 e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|\\ 383 e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)\; \cos\varphi(j+1/2) |& 384 e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)| 385 \end{flalign*} 373 \begin{align*} 374 \lambda_t &\equiv \text{glamt} = \lambda (i ) 375 &\varphi_t &\equiv \text{gphit} = \varphi (j ) \\ 376 \lambda_u &\equiv \text{glamu} = \lambda (i + 1/2) 377 &\varphi_u &\equiv \text{gphiu} = \varphi (j ) \\ 378 \lambda_v &\equiv \text{glamv} = \lambda (i ) 379 &\varphi_v &\equiv \text{gphiv} = \varphi (j + 1/2) \\ 380 \lambda_f &\equiv \text{glamf} = \lambda (i + 1/2) 381 &\varphi_f &\equiv \text{gphif} = \varphi (j + 1/2) \\ 382 e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i ) \; \cos\varphi(j ) | 383 &e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j ) | \\ 384 e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i + 1/2) \; \cos\varphi(j ) | 385 &e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j ) | \\ 386 e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i ) \; \cos\varphi(j + 1/2) | 387 &e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j + 1/2) | \\ 388 e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) | 389 &e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j + 1/2) | 390 \end{align*} 386 391 where the last letter of each computational name indicates the grid point considered and 387 392 $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). 388 393 Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points, 389 thus no specific arrays are defined at $w$-points. 394 thus no specific arrays are defined at $w$-points. 390 395 391 396 Note that the definition of the scale factors … … 405 410 \begin{figure}[!t] 406 411 \begin{center} 407 \includegraphics[ width=0.90\textwidth]{Fig_zgr_e3}412 \includegraphics[]{Fig_zgr_e3} 408 413 \caption{ 409 414 \protect\label{fig:zgr_e3} 410 415 Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 411 416 and (b) analytically derived grid-point position and scale factors. 412 For both grids here, 413 the same $w$-point depth has been chosen but in (a) the $t$-points are set half way between $w$-points while 414 in (b) they are defined from an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$. 417 For both grids here, the same $w$-point depth has been chosen but 418 in (a) the $t$-points are set half way between $w$-points while 419 in (b) they are defined from an analytical function: 420 $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 415 421 Note the resulting difference between the value of the grid-size $\Delta_k$ and 416 422 those of the scale factor $e_k$. … … 426 432 \label{subsec:DOM_hgr_msh_choice} 427 433 428 429 434 % ------------------------------------------------------------------------------------------------------------- 430 435 % Grid files … … 434 439 435 440 All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks) 436 can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}).441 can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}). 437 442 This can be particularly useful for plots and off-line diagnostics. 438 443 In some cases, the user may choose to make a local modification of a scale factor in the code. … … 441 446 An example is Gibraltar Strait in the ORCA2 configuration. 442 447 When such modifications are done, 443 the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid.448 the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid. 444 449 445 450 % ================================================================ … … 466 471 \begin{figure}[!tb] 467 472 \begin{center} 468 \includegraphics[ width=1.0\textwidth]{Fig_z_zps_s_sps}473 \includegraphics[]{Fig_z_zps_s_sps} 469 474 \caption{ 470 475 \protect\label{fig:z_zps_s_sps} … … 475 480 (d) hybrid $s-z$ coordinate, 476 481 (e) hybrid $s-z$ coordinate with partial step, and 477 (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh} \forcode{= .false.}).482 (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}). 478 483 Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e). 479 484 } … … 485 490 must be done once of all at the beginning of an experiment. 486 491 It is not intended as an option which can be enabled or disabled in the middle of an experiment. 487 Three main choices are offered (\autoref{fig:z_zps_s_sps} a to c):488 $z$-coordinate with full step bathymetry (\np{ln\_zco} \forcode{= .true.}),489 $z$-coordinate with partial step bathymetry (\np{ln\_zps} \forcode{= .true.}),490 or generalized, $s$-coordinate (\np{ln\_sco} \forcode{= .true.}).492 Three main choices are offered (\autoref{fig:z_zps_s_sps}): 493 $z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}), 494 $z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}), 495 or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}). 491 496 Hybridation of the three main coordinates are available: 492 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps} e).497 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}). 493 498 By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that 494 the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps} f).495 When a linear free surface is assumed (\np{ln\_linssh} \forcode{= .true.}),496 the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0surface497 (in other words, the top of the ocean in not a rigid-lid). 499 the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}). 500 When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}), 501 the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface 502 (in other words, the top of the ocean in not a rigid-lid). 498 503 The last choice in terms of vertical coordinate concerns the presence (or not) in 499 504 the model domain of ocean cavities beneath ice shelves. … … 502 507 and partial step are also applied at the ocean/ice shelf interface. 503 508 504 Contrary to the horizontal grid, the vertical grid is computed in the code and 505 no provision is made forreading it from a file.509 Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for 510 reading it from a file. 506 511 The only input file is the bathymetry (in meters) (\ifile{bathy\_meter}) 507 512 \footnote{ 508 513 N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file, 509 so that the computation of the number of wet ocean point in each water column is by-passed 510 }. 511 If \np{ln\_isfcav}\forcode{ = .true.}, 512 an extra file input file describing the ice shelf draft (in meters) (\ifile{isf\_draft\_meter}) is needed. 514 so that the computation of the number of wet ocean point in each water column is by-passed}. 515 If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing 516 the ice shelf draft (in meters) is needed. 513 517 514 518 After reading the bathymetry, the algorithm for vertical grid definition differs between the different options: 515 519 \begin{description} 516 520 \item[\textit{zco}] 517 set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0(k)$.521 set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$. 518 522 \item[\textit{zps}] 519 set a reference coordinate transformation $z_0 (k)$, 520 and calculate the thickness of the deepest level at each $(i,j)$ point using the bathymetry, 521 to obtain the final three-dimensional depth and scale factor arrays. 523 set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at 524 each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays. 522 525 \item[\textit{sco}] 523 smooth the bathymetry to fulfil the hydrostatic consistency criteria and526 smooth the bathymetry to fulfill the hydrostatic consistency criteria and 524 527 set the three-dimensional transformation. 525 528 \item[\textit{s-z} and \textit{s-zps}] 526 smooth the bathymetry to fulfil the hydrostatic consistency criteria and529 smooth the bathymetry to fulfill the hydrostatic consistency criteria and 527 530 set the three-dimensional transformation $z(i,j,k)$, 528 531 and possibly introduce masking of extra land points to better fit the original bathymetry file. … … 532 535 %%% 533 536 534 Unless a linear free surface is used (\np{ln\_linssh} \forcode{= .false.}),537 Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}), 535 538 the arrays describing the grid point depths and vertical scale factors are three set of 536 539 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 537 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively.540 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 538 541 They are updated at each model time step using a fixed reference coordinate system which 539 542 computer names have a $\_0$ suffix. 540 When the linear free surface option is used (\np{ln\_linssh}\forcode{ = .true.}), 541 \textit{before}, \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 542 543 When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before}, 544 \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 543 545 544 546 % ------------------------------------------------------------------------------------------------------------- … … 551 553 (found in \ngn{namdom} namelist): 552 554 \begin{description} 553 \item[\np{nn\_bathy} \forcode{= 0}]:555 \item[\np{nn\_bathy}~\forcode{= 0}]: 554 556 a flat-bottom domain is defined. 555 557 The total depth $z_w (jpk)$ is given by the coordinate transformation. 556 The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 557 \item[\np{nn\_bathy} \forcode{= -1}]:558 The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 559 \item[\np{nn\_bathy}~\forcode{= -1}]: 558 560 a domain with a bump of topography one third of the domain width at the central latitude. 559 This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 560 \item[\np{nn\_bathy} \forcode{= 1}]:561 This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 562 \item[\np{nn\_bathy}~\forcode{= 1}]: 561 563 read a bathymetry and ice shelf draft (if needed). 562 564 The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at … … 569 571 The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at 570 572 each grid point of the model grid. 571 This file is only needed if \np{ln\_isfcav} \forcode{= .true.}.573 This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}. 572 574 Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. 573 575 \end{description} 574 576 575 577 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 576 ( e.g, great lakes, Caspian sea...)577 even if the model resolution does not allow their communication withthe rest of the ocean.578 (\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with 579 the rest of the ocean. 578 580 This is unnecessary when the ocean is forced by fixed atmospheric conditions, 579 581 so these seas can be removed from the ocean domain. 580 582 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), 581 but the code has to be adapted to the user's configuration. 583 but the code has to be adapted to the user's configuration. 582 584 583 585 % ------------------------------------------------------------------------------------------------------------- 584 586 % z-coordinate and reference coordinate transformation 585 587 % ------------------------------------------------------------------------------------------------------------- 586 \subsection[$Z$-coordinate (\protect\np{ln\_zco} \forcode{= .true.}) and ref. coordinate]587 {$Z$-coordinate (\protect\np{ln\_zco} \forcode{= .true.}) and reference coordinate}588 \subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate] 589 {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate} 588 590 \label{subsec:DOM_zco} 589 591 … … 591 593 \begin{figure}[!tb] 592 594 \begin{center} 593 \includegraphics[ width=0.90\textwidth]{Fig_zgr}595 \includegraphics[]{Fig_zgr} 594 596 \caption{ 595 597 \protect\label{fig:zgr} … … 602 604 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 603 605 604 The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ and $gdepw_0$ for 605 $t$- and $w$-points, respectively. 606 As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the ocean surface. 606 The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points, 607 respectively. 608 As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. 609 $gdepw_0(1)$ is the ocean surface. 607 610 There are at most \jp{jpk}-1 $t$-points inside the ocean, 608 the additional $t$-point at $jk =jpk$ is below the sea floor and is not used.611 the additional $t$-point at $jk = jpk$ is below the sea floor and is not used. 609 612 The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose 610 613 analytical derivative with respect to $k$ provides the vertical scale factors. … … 613 616 using parameters provided in the \ngn{namcfg} namelist. 614 617 615 It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr =0}).616 In that case, 617 the parameters \jp{jpk} (number of $w$-levels) and \np{pphmax} (total ocean depth in meters) fully define the grid. 618 It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}). 619 In that case, the parameters \jp{jpk} (number of $w$-levels) and 620 \np{pphmax} (total ocean depth in meters) fully define the grid. 618 621 619 622 For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface. 620 623 The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps): 621 \begin{ equation}624 \begin{gather} 622 625 \label{eq:DOM_zgr_ana_1} 623 \begin{split} 624 z_0 (k) &= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ 625 e_3^0 (k) &= \left| -h_0 -h_1 \;\tanh \left( {{(k-h_{th} )} / {h_{cr} }} \right) \right| 626 \end{split} 627 \end{equation} 628 where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. 626 z_0 (k) = h_{sur} - h_0 \; k - \; h_1 \; \log \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\ 627 e_3^0(k) = \lt| - h_0 - h_1 \; \tanh \big[ (k - h_{th}) / h_{cr} \big] \rt| 628 \end{gather} 629 where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels. 629 630 Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with 630 631 a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}). 631 632 632 If the ice shelf cavities are opened (\np{ln\_isfcav} \forcode{= .true.}), the definition of $z_0$ is the same.633 If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same. 633 634 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 634 635 \begin{equation} 635 636 \label{eq:DOM_zgr_ana_2} 636 637 \begin{split} 637 e_3^T(k) &= z_W (k +1) - z_W (k)\\638 e_3^W(k) &= z_T (k ) - z_T (k-1) \\638 e_3^T(k) &= z_W (k + 1) - z_W (k ) \\ 639 e_3^W(k) &= z_T (k ) - z_T (k - 1) 639 640 \end{split} 640 641 \end{equation} 641 642 This formulation decrease the self-generated circulation into the ice shelf cavity 642 643 (which can, in extreme case, leads to blow up).\\ 643 644 644 645 The most used vertical grid for ORCA2 has $10~m$ ($500~m )$resolution in the surface (bottom) layers and645 The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and 646 646 a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. 647 647 This leads to the following conditions: 648 648 \begin{equation} 649 649 \label{eq:DOM_zgr_coef} 650 \begin{split} 651 e_3 (1+1/2) &=10. \\ 652 e_3 (jpk-1/2) &=500. \\ 653 z(1) &=0. \\ 654 z(jpk) &=-5000. \\ 655 \end{split} 650 \begin{array}{ll} 651 e_3 (1 + 1/2) = 10. & z(1 ) = 0. \\ 652 e_3 (jpk - 1/2) = 500. & z(jpk) = -5000. 653 \end{array} 656 654 \end{equation} 657 655 658 With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}=$31$,659 the four coefficients $h_{sur}$, $h_ {0}$, $h_{1}$, and $h_{th}$ in656 With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$, 657 the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in 660 658 \autoref{eq:DOM_zgr_ana_2} have been determined such that 661 659 \autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. 662 660 For the first standard ORCA2 vertical grid this led to the following values: 663 $h_{sur} = 4762.96$, $h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$.661 $h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$. 664 662 The resulting depths and scale factors as a function of the model levels are shown in 665 663 \autoref{fig:zgr} and given in \autoref{tab:orca_zgr}. 666 Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 667 668 Rather than entering parameters $h_{sur}$, $h_ {0}$, and $h_{1}$ directly, it is possible to recalculate them.669 In that case the user sets \np{ppsur} \forcode{ = }\np{ppa0}\forcode{ = }\np{ppa1}\forcode{ = 999999}.,664 Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 665 666 Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them. 667 In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$., 670 668 in \ngn{namcfg} namelist, and specifies instead the four following parameters: 671 669 \begin{itemize} 672 670 \item 673 \np{ppacr} =$h_{cr}$: stretching factor (nondimensional).671 \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional). 674 672 The larger \np{ppacr}, the smaller the stretching. 675 673 Values from $3$ to $10$ are usual. 676 674 \item 677 \np{ppkth} =$h_{th}$: is approximately the model level at which maximum stretching occurs675 \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs 678 676 (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) 679 677 \item … … 683 681 \end{itemize} 684 682 As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are: 685 \jp{jpk}\forcode{ = 46}, \np{ppacr}\forcode{ = 9}, \np{ppkth}\forcode{ = 23.563}, 686 \np{ppdzmin}\forcode{ = 6}m, \np{pphmax}\forcode{ = 5750}m. 683 \jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$. 687 684 688 685 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 691 688 \begin{tabular}{c||r|r|r|r} 692 689 \hline 693 \textbf{LEVEL}& \textbf{gdept\_1d}& \textbf{gdepw\_1d}& \textbf{e3t\_1d }& \textbf{e3w\_1d } \\ \hline 694 1 & \textbf{ 5.00} & 0.00 & \textbf{ 10.00} & 10.00 \\ \hline 695 2 & \textbf{15.00} & 10.00 & \textbf{ 10.00} & 10.00 \\ \hline 696 3 & \textbf{25.00} & 20.00 & \textbf{ 10.00} & 10.00 \\ \hline 697 4 & \textbf{35.01} & 30.00 & \textbf{ 10.01} & 10.00 \\ \hline 698 5 & \textbf{45.01} & 40.01 & \textbf{ 10.01} & 10.01 \\ \hline 699 6 & \textbf{55.03} & 50.02 & \textbf{ 10.02} & 10.02 \\ \hline 700 7 & \textbf{65.06} & 60.04 & \textbf{ 10.04} & 10.03 \\ \hline 701 8 & \textbf{75.13} & 70.09 & \textbf{ 10.09} & 10.06 \\ \hline 702 9 & \textbf{85.25} & 80.18 & \textbf{ 10.17} & 10.12 \\ \hline 703 10 & \textbf{95.49} & 90.35 & \textbf{ 10.33} & 10.24 \\ \hline 704 11 & \textbf{105.97} & 100.69 & \textbf{ 10.65} & 10.47 \\ \hline 705 12 & \textbf{116.90} & 111.36 & \textbf{ 11.27} & 10.91 \\ \hline 706 13 & \textbf{128.70} & 122.65 & \textbf{ 12.47} & 11.77 \\ \hline 707 14 & \textbf{142.20} & 135.16 & \textbf{ 14.78} & 13.43 \\ \hline 708 15 & \textbf{158.96} & 150.03 & \textbf{ 19.23} & 16.65 \\ \hline 709 16 & \textbf{181.96} & 169.42 & \textbf{ 27.66} & 22.78 \\ \hline 710 17 & \textbf{216.65} & 197.37 & \textbf{ 43.26} & 34.30 \\ \hline 711 18 & \textbf{272.48} & 241.13 & \textbf{ 70.88} & 55.21 \\ \hline 712 19 & \textbf{364.30} & 312.74 & \textbf{116.11} & 90.99 \\ \hline 713 20 & \textbf{511.53} & 429.72 & \textbf{181.55} & 146.43 \\ \hline 714 21 & \textbf{732.20} & 611.89 & \textbf{261.03} & 220.35 \\ \hline 715 22 & \textbf{1033.22}& 872.87 & \textbf{339.39} & 301.42 \\ \hline 716 23 & \textbf{1405.70}& 1211.59 & \textbf{402.26} & 373.31 \\ \hline 717 24 & \textbf{1830.89}& 1612.98 & \textbf{444.87} & 426.00 \\ \hline 718 25 & \textbf{2289.77}& 2057.13 & \textbf{470.55} & 459.47 \\ \hline 719 26 & \textbf{2768.24}& 2527.22 & \textbf{484.95} & 478.83 \\ \hline 720 27 & \textbf{3257.48}& 3011.90 & \textbf{492.70} & 489.44 \\ \hline 721 28 & \textbf{3752.44}& 3504.46 & \textbf{496.78} & 495.07 \\ \hline 722 29 & \textbf{4250.40}& 4001.16 & \textbf{498.90} & 498.02 \\ \hline 723 30 & \textbf{4749.91}& 4500.02 & \textbf{500.00} & 499.54 \\ \hline 724 31 & \textbf{5250.23}& 5000.00 & \textbf{500.56} & 500.33 \\ \hline 690 \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\ 691 \hline 692 1 & \textbf{ 5.00} & 0.00 & \textbf{ 10.00} & 10.00 \\ 693 \hline 694 2 & \textbf{ 15.00} & 10.00 & \textbf{ 10.00} & 10.00 \\ 695 \hline 696 3 & \textbf{ 25.00} & 20.00 & \textbf{ 10.00} & 10.00 \\ 697 \hline 698 4 & \textbf{ 35.01} & 30.00 & \textbf{ 10.01} & 10.00 \\ 699 \hline 700 5 & \textbf{ 45.01} & 40.01 & \textbf{ 10.01} & 10.01 \\ 701 \hline 702 6 & \textbf{ 55.03} & 50.02 & \textbf{ 10.02} & 10.02 \\ 703 \hline 704 7 & \textbf{ 65.06} & 60.04 & \textbf{ 10.04} & 10.03 \\ 705 \hline 706 8 & \textbf{ 75.13} & 70.09 & \textbf{ 10.09} & 10.06 \\ 707 \hline 708 9 & \textbf{ 85.25} & 80.18 & \textbf{ 10.17} & 10.12 \\ 709 \hline 710 10 & \textbf{ 95.49} & 90.35 & \textbf{ 10.33} & 10.24 \\ 711 \hline 712 11 & \textbf{ 105.97} & 100.69 & \textbf{ 10.65} & 10.47 \\ 713 \hline 714 12 & \textbf{ 116.90} & 111.36 & \textbf{ 11.27} & 10.91 \\ 715 \hline 716 13 & \textbf{ 128.70} & 122.65 & \textbf{ 12.47} & 11.77 \\ 717 \hline 718 14 & \textbf{ 142.20} & 135.16 & \textbf{ 14.78} & 13.43 \\ 719 \hline 720 15 & \textbf{ 158.96} & 150.03 & \textbf{ 19.23} & 16.65 \\ 721 \hline 722 16 & \textbf{ 181.96} & 169.42 & \textbf{ 27.66} & 22.78 \\ 723 \hline 724 17 & \textbf{ 216.65} & 197.37 & \textbf{ 43.26} & 34.30 \\ 725 \hline 726 18 & \textbf{ 272.48} & 241.13 & \textbf{ 70.88} & 55.21 \\ 727 \hline 728 19 & \textbf{ 364.30} & 312.74 & \textbf{ 116.11} & 90.99 \\ 729 \hline 730 20 & \textbf{ 511.53} & 429.72 & \textbf{ 181.55} & 146.43 \\ 731 \hline 732 21 & \textbf{ 732.20} & 611.89 & \textbf{ 261.03} & 220.35 \\ 733 \hline 734 22 & \textbf{ 1033.22} & 872.87 & \textbf{ 339.39} & 301.42 \\ 735 \hline 736 23 & \textbf{ 1405.70} & 1211.59 & \textbf{ 402.26} & 373.31 \\ 737 \hline 738 24 & \textbf{ 1830.89} & 1612.98 & \textbf{ 444.87} & 426.00 \\ 739 \hline 740 25 & \textbf{ 2289.77} & 2057.13 & \textbf{ 470.55} & 459.47 \\ 741 \hline 742 26 & \textbf{ 2768.24} & 2527.22 & \textbf{ 484.95} & 478.83 \\ 743 \hline 744 27 & \textbf{ 3257.48} & 3011.90 & \textbf{ 492.70} & 489.44 \\ 745 \hline 746 28 & \textbf{ 3752.44} & 3504.46 & \textbf{ 496.78} & 495.07 \\ 747 \hline 748 29 & \textbf{ 4250.40} & 4001.16 & \textbf{ 498.90} & 498.02 \\ 749 \hline 750 30 & \textbf{ 4749.91} & 4500.02 & \textbf{ 500.00} & 499.54 \\ 751 \hline 752 31 & \textbf{ 5250.23} & 5000.00 & \textbf{ 500.56} & 500.33 \\ 753 \hline 725 754 \end{tabular} 726 755 \end{center} … … 736 765 % z-coordinate with partial step 737 766 % ------------------------------------------------------------------------------------------------------------- 738 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps} \forcode{= .true.})}767 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})} 739 768 \label{subsec:DOM_zps} 740 769 %--------------------------------------------namdom------------------------------------------------------- … … 744 773 745 774 In $z$-coordinate partial step, 746 the depths of the model levels are defined by the reference analytical function $z_0 747 the previous section, \ emph{except} in the bottom layer.775 the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in 776 the previous section, \textit{except} in the bottom layer. 748 777 The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to 749 778 allow a better representation of the bathymetry, especially in the case of small slopes … … 752 781 With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. 753 782 The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: 754 the maximum thickness allowed is $2*e_{3t}(jpk -1)$.783 the maximum thickness allowed is $2*e_{3t}(jpk - 1)$. 755 784 This has to be kept in mind when specifying values in \ngn{namdom} namelist, 756 785 as the maximum depth \np{pphmax} in partial steps: 757 for example, with \np{pphmax} $=5750~m$ for the DRAKKAR 45 layer grid,758 the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk -1)$ being $250~m$).786 for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, 787 the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$). 759 788 Two variables in the namdom namelist are used to define the partial step vertical grid. 760 789 The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is … … 767 796 % s-coordinate 768 797 % ------------------------------------------------------------------------------------------------------------- 769 \subsection{$S$-coordinate (\protect\np{ln\_sco} \forcode{= .true.})}798 \subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})} 770 799 \label{subsec:DOM_sco} 771 800 %------------------------------------------nam_zgr_sco--------------------------------------------------- … … 774 803 %-------------------------------------------------------------------------------------------------------------- 775 804 Options are defined in \ngn{namzgr\_sco}. 776 In $s$-coordinate (\np{ln\_sco} \forcode{= .true.}), the depth and thickness of the model levels are defined from805 In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from 777 806 the product of a depth field and either a stretching function or its derivative, respectively: 778 807 779 \ [808 \begin{align*} 780 809 % \label{eq:DOM_sco_ana} 781 \begin{split} 782 z(k) &= h(i,j) \; z_0(k) \\ 783 e_3(k) &= h(i,j) \; z_0'(k) 784 \end{split} 785 \] 810 z(k) &= h(i,j) \; z_0 (k) \\ 811 e_3(k) &= h(i,j) \; z_0'(k) 812 \end{align*} 786 813 787 814 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and … … 789 816 The depth field $h$ is not necessary the ocean depth, 790 817 since a mixed step-like and bottom-following representation of the topography can be used 791 (\autoref{fig:z_zps_s_sps} d-e) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}f).818 (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). 792 819 The namelist parameter \np{rn\_rmax} determines the slope at which 793 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 820 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 794 821 The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as 795 822 the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. … … 799 826 800 827 The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true 801 (\np{ln\_s\_SH94} \forcode{ = .false.} and \np{ln\_s\_SF12}\forcode{ = .false.}).828 (\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}). 802 829 This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: 803 830 804 831 \[ 805 z = s_{min} +C\left(s\right)\left(H-s_{min}\right)832 z = s_{min} + C (s) (H - s_{min}) 806 833 % \label{eq:SH94_1} 807 834 \] … … 810 837 allows a $z$-coordinate to placed on top of the stretched coordinate, 811 838 and $z$ is the depth (negative down from the asea surface). 839 \begin{gather*} 840 s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1 841 % \label{eq:DOM_s} 842 \\ 843 % \label{eq:DOM_sco_function} 844 C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)} 845 \end{gather*} 846 847 A stretching function, 848 modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}), 849 is also available and is more commonly used for shelf seas modelling: 812 850 813 851 \[ 814 s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1 815 % \label{eq:DOM_s} 816 \] 817 818 \[ 819 % \label{eq:DOM_sco_function} 820 \begin{split} 821 C(s) &= \frac{ \left[ \tanh{ \left( \theta \, (s+b) \right)} 822 - \tanh{ \left( \theta \, b \right)} \right]} 823 {2\;\sinh \left( \theta \right)} 824 \end{split} 825 \] 826 827 A stretching function, 828 modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}\forcode{ = .true.}), 829 is also available and is more commonly used for shelf seas modelling: 830 831 \[ 832 C\left(s\right) = \left(1 - b \right)\frac{ \sinh\left( \theta s\right)}{\sinh\left(\theta\right)} + \\ 833 b\frac{ \tanh \left[ \theta \left(s + \frac{1}{2} \right)\right] - \tanh\left(\frac{\theta}{2}\right)}{ 2\tanh\left (\frac{\theta}{2}\right)} 852 C(s) = (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)} 853 + b \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] - \tanh \lt( \frac{\theta}{2} \rt)} 854 { 2 \tanh \lt( \frac{\theta}{2} \rt)} 834 855 % \label{eq:SH94_2} 835 856 \] … … 838 859 \begin{figure}[!ht] 839 860 \begin{center} 840 \includegraphics[ width=1.0\textwidth]{Fig_sco_function}861 \includegraphics[]{Fig_sco_function} 841 862 \caption{ 842 863 \protect\label{fig:sco_function} … … 848 869 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 849 870 850 where $H_c$ is the critical depth (\np{rn\_hc}) at which 851 the coordinate transitions from pure $\sigma$ to the stretched coordinate, 852 and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and bottom control parameters such that 853 $0\leqslant \theta \leqslant 20$, and $0\leqslant b\leqslant 1$. 871 where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to 872 the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and 873 bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. 854 874 $b$ has been designed to allow surface and/or bottom increase of the vertical resolution 855 875 (\autoref{fig:sco_function}). … … 859 879 In this case the a stretching function $\gamma$ is defined such that: 860 880 861 \ [862 z = - \gamma h \quad \text{ with} \quad 0 \leq \gamma \leq 1881 \begin{equation} 882 z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 863 883 % \label{eq:z} 864 \ ]884 \end{equation} 865 885 866 886 The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 867 887 868 \ [888 \begin{gather*} 869 889 % \label{eq:DOM_gamma_deriv} 870 \gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right) 871 \] 872 873 Where: 874 \[ 890 \gamma = A \lt( \sigma - \frac{1}{2} (\sigma^2 + f (\sigma)) \rt) 891 + B \lt( \sigma^3 - f (\sigma) \rt) + f (\sigma) \\ 892 \intertext{Where:} 875 893 % \label{eq:DOM_gamma} 876 f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} \quad \text{ and } \quad \sigma = \frac{k}{n-1} 877 \] 894 f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2} 895 \quad \text{and} \quad \sigma = \frac{k}{n - 1} 896 \end{gather*} 878 897 879 898 This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of … … 892 911 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 893 912 \begin{figure}[!ht] 894 \includegraphics[width=1.0\textwidth]{Fig_DOM_compare_coordinates_surface}895 896 897 898 the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in899 the surface 100m for a idealised bathymetry that goes from 50m to 5500mdepth.900 901 902 913 \includegraphics[]{Fig_DOM_compare_coordinates_surface} 914 \caption{ 915 A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), 916 a 50 level $Z$-coordinate (contoured surfaces) and 917 the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for 918 a idealised bathymetry that goes from $50~m$ to $5500~m$ depth. 919 For clarity every third coordinate surface is shown. 920 } 921 \label{fig:fig_compare_coordinates_surface} 903 922 \end{figure} 904 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>923 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 905 924 906 925 This gives a smooth analytical stretching in computational space that is constrained to … … 925 944 926 945 % ------------------------------------------------------------------------------------------------------------- 927 % \zstar- or \sstar-coordinate928 % ------------------------------------------------------------------------------------------------------------- 929 \subsection{ $Z^*$- or $S^*$-coordinate (\protect\np{ln\_linssh}\forcode{ = .false.})}946 % z*- or s*-coordinate 947 % ------------------------------------------------------------------------------------------------------------- 948 \subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})} 930 949 \label{subsec:DOM_zgr_star} 931 950 932 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 951 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 933 952 934 953 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances … … 940 959 \label{subsec:DOM_msk} 941 960 942 Whatever the vertical coordinate used, 943 the model offers the possibility of representing the bottom topography with steps that 944 follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. 945 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, 946 which gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 947 mbathy is computed from the meter bathymetry using the definiton of gdept as 948 the number of $t$-points which gdept $\leq$ bathy. 961 Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with 962 steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. 963 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which 964 gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 965 mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which 966 gdept $\leq$ bathy. 949 967 950 968 Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after … … 954 972 As for the representation of bathymetry, a 2D integer array, misfdep, is created. 955 973 misfdep defines the level of the first wet $t$-point. 956 All the cells between $k =1$ and $misfdep(i,j)-1$ are masked.957 By default, misfdep(:,:)=1and no cells are masked.974 All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked. 975 By default, $misfdep(:,:) = 1$ and no cells are masked. 958 976 959 977 In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into 960 978 the cavities are performed in the \textit{zgr\_isf} routine. 961 The compatibility between ice shelf draft and bathymetry is checked. 962 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 979 The compatibility between ice shelf draft and bathymetry is checked. 980 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 963 981 If only one cell on the water column is opened at $t$-, $u$- or $v$-points, 964 982 the bathymetry or the ice shelf draft is dug to fit this constrain. 965 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. \\983 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. 966 984 967 985 From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: 968 \begin{align*} 969 tmask(i,j,k) &= \begin{cases} \; 0& \text{ if $k < misfdep(i,j) $ } \\ 970 \; 1& \text{ if $misfdep(i,j) \leq k\leq mbathy(i,j)$ } \\ 971 \; 0& \text{ if $k > mbathy(i,j)$ } \end{cases} \\ 972 umask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 973 vmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i,j+1,k) \\ 974 fmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 975 & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 976 wmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1) 977 \end{align*} 986 \begin{alignat*}{2} 987 tmask(i,j,k) &= & & 988 \begin{cases} 989 0 &\text{if $ k < misfdep(i,j)$} \\ 990 1 &\text{if $misfdep(i,j) \leq k \leq mbathy(i,j)$} \\ 991 0 &\text{if $ k > mbathy(i,j)$} 992 \end{cases} 993 \\ 994 umask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ 995 vmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j + 1,k) \\ 996 fmask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ 997 & &* &tmask(i,j,k) * tmask(i + 1,j, k) \\ 998 wmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j,k - 1) \\ 999 \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) 1000 \end{alignat*} 978 1001 979 1002 Note that, without ice shelves cavities, 980 1003 masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). 981 1004 Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 982 exactly in the same way as for the bottom boundary. 1005 exactly in the same way as for the bottom boundary. 983 1006 984 1007 The specification of closed lateral boundaries requires that at least 985 1008 the first and last rows and columns of the \textit{mbathy} array are set to zero. 986 In the particular case of an east-west cyclical boundary condition, 987 \textit{mbathy} has its last column equal to the second one and its first column equal to the last but one 988 (and so too the mask arrays) (see \autoref{fig:LBC_jperio}). 989 1009 In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 1010 the second one and its first column equal to the last but one (and so too the mask arrays) 1011 (see \autoref{fig:LBC_jperio}). 990 1012 991 1013 % ================================================================ … … 1000 1022 1001 1023 Options are defined in \ngn{namtsd}. 1002 By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 1024 By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and 1025 salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 1003 1026 \begin{description} 1004 \item[\np{ln\_tsd\_init} \forcode{= .true.}]1027 \item[\np{ln\_tsd\_init}~\forcode{= .true.}] 1005 1028 use a T and S input files that can be given on the model grid itself or on their native input data grid. 1006 1029 In the latter case, … … 1009 1032 The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. 1010 1033 The computation is done in the \mdl{dtatsd} module. 1011 \item[\np{ln\_tsd\_init} \forcode{= .false.}]1012 use constant salinity value of 35.5 psu and an analytical profile of temperature (typical of the tropical ocean),1013 see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.1034 \item[\np{ln\_tsd\_init}~\forcode{= .false.}] 1035 use constant salinity value of $35.5~psu$ and an analytical profile of temperature 1036 (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. 1014 1037 \end{description} 1015 1038 -
NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex
r10442 r10502 15 15 %STEVEN : is the use of the word "positive" to describe a scheme enough, or should it be "positive definite"? I added a comment to this effect on some instances of this below 16 16 17 %\newpage 18 19 Using the representation described in \autoref{chap:DOM}, 20 several semi-discrete space forms of the tracer equations are available depending on 21 the vertical coordinate used and on the physics used. 17 Using the representation described in \autoref{chap:DOM}, several semi -discrete space forms of 18 the tracer equations are available depending on the vertical coordinate used and on the physics used. 22 19 In all the equations presented here, the masking has been omitted for simplicity. 23 One must be aware that all the quantities are masked fields and 24 that each time a mean or difference operator is used, 25 the resulting field is multiplied by a mask. 20 One must be aware that all the quantities are masked fields and that each time a mean or 21 difference operator is used, the resulting field is multiplied by a mask. 26 22 27 23 The two active tracers are potential temperature and salinity. 28 24 Their prognostic equations can be summarized as follows: 29 25 \[ 30 \text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC}31 \ (+\text{QSR})\ (+\text{BBC})\ (+\text{BBL})\ (+\text{DMP})26 \text{NXT} = \text{ADV} + \text{LDF} + \text{ZDF} + \text{SBC} 27 + \{\text{QSR}, \text{BBC}, \text{BBL}, \text{DMP}\} 32 28 \] 33 29 … … 39 35 The terms QSR, BBC, BBL and DMP are optional. 40 36 The external forcings and parameterisations require complex inputs and complex calculations 41 (\eg bulk formulae, estimation of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and 42 described in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 43 Note that \mdl{tranpc}, the non-penetrative convection module, although located in the NEMO/OPA/TRA directory as 44 it directly modifies the tracer fields, is described with the model vertical physics (ZDF) together with 37 (\eg bulk formulae, estimation of mixing coefficients) that are carried out in the SBC, 38 LDF and ZDF modules and described in \autoref{chap:SBC}, \autoref{chap:LDF} and 39 \autoref{chap:ZDF}, respectively. 40 Note that \mdl{tranpc}, the non-penetrative convection module, although located in 41 the \path{./src/OCE/TRA} directory as it directly modifies the tracer fields, 42 is described with the model vertical physics (ZDF) together with 45 43 other available parameterization of convection. 46 44 … … 50 48 51 49 The different options available to the user are managed by namelist logicals or CPP keys. 52 For each equation term 50 For each equation term \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 53 51 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 54 52 The CPP key (when it exists) is \key{traTTT}. 55 53 The equivalent code can be found in the \textit{traTTT} or \textit{traTTT\_xxx} module, 56 in the NEMO/OPA/TRAdirectory.54 in the \path{./src/OCE/TRA} directory. 57 55 58 56 The user has the option of extracting each tendency term on the RHS of the tracer equation for output 59 (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl} \forcode{= .true.}), as described in \autoref{chap:DIA}.57 (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~\forcode{= .true.}), as described in \autoref{chap:DIA}. 60 58 61 59 % ================================================================ … … 75 73 \begin{equation} 76 74 \label{eq:tra_adv} 77 ADV_\tau =-\frac{1}{b_t} \left( 78 \;\delta_i \left[ e_{2u}\,e_{3u} \; u\; \tau_u \right] 79 +\delta_j \left[ e_{1v}\,e_{3v} \; v\; \tau_v \right] \; \right) 80 -\frac{1}{e_{3t}} \;\delta_k \left[ w\; \tau_w \right] 75 ADV_\tau = - \frac{1}{b_t} \Big( \delta_i [ e_{2u} \, e_{3u} \; u \; \tau_u] 76 + \delta_j [ e_{1v} \, e_{3v} \; v \; \tau_v] \Big) 77 - \frac{1}{e_{3t}} \delta_k [w \; \tau_w] 81 78 \end{equation} 82 where $\tau$ is either T or S, and $b_t = e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells.79 where $\tau$ is either T or S, and $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells. 83 80 The flux form in \autoref{eq:tra_adv} implicitly requires the use of the continuity equation. 84 Indeed, it is obtained by using the following equality: 85 $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ which 86 results from the use of the continuity equation, $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 87 (which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, \ie \np{ln\_linssh}\forcode{ = .true.}). 81 Indeed, it is obtained by using the following equality: $\nabla \cdot (\vect U \, T) = \vect U \cdot \nabla T$ which 82 results from the use of the continuity equation, $\partial_t e_3 + e_3 \; \nabla \cdot \vect U = 0$ 83 (which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, \ie \np{ln\_linssh}~\forcode{= .true.}). 88 84 Therefore it is of paramount importance to design the discrete analogue of the advection tendency so that 89 85 it is consistent with the continuity equation in order to enforce the conservation properties of … … 94 90 \begin{figure}[!t] 95 91 \begin{center} 96 \includegraphics[ width=0.9\textwidth]{Fig_adv_scheme}92 \includegraphics[]{Fig_adv_scheme} 97 93 \caption{ 98 94 \protect\label{fig:adv_scheme} … … 120 116 since the normal velocity is zero there. 121 117 At the sea surface the boundary condition depends on the type of sea surface chosen: 118 122 119 \begin{description} 123 120 \item[linear free surface:] 124 (\np{ln\_linssh} \forcode{= .true.})121 (\np{ln\_linssh}~\forcode{= .true.}) 125 122 the first level thickness is constant in time: 126 the vertical boundary condition is applied at the fixed surface $z=0$ rather than on the moving surface $z=\eta$. 123 the vertical boundary condition is applied at the fixed surface $z = 0$ rather than on 124 the moving surface $z = \eta$. 127 125 There is a non-zero advective flux which is set for all advection schemes as 128 $\ left. {\tau_w } \right|_{k=1/2} =T_{k=1} $,129 \ie the product of surface velocity (at $z=0$) bythe first level tracer value.126 $\tau_w|_{k = 1/2} = T_{k = 1}$, \ie the product of surface velocity (at $z = 0$) by 127 the first level tracer value. 130 128 \item[non-linear free surface:] 131 (\np{ln\_linssh} \forcode{= .false.})129 (\np{ln\_linssh}~\forcode{= .false.}) 132 130 convergence/divergence in the first ocean level moves the free surface up/down. 133 131 There is no tracer advection through it so that the advective fluxes through the surface are also zero. 134 132 \end{description} 133 135 134 In all cases, this boundary condition retains local conservation of tracer. 136 135 Global conservation is obtained in non-linear free surface case, but \textit{not} in the linear free surface case. 137 136 Nevertheless, in the latter case, it is achieved to a good approximation since 138 137 the non-conservative term is the product of the time derivative of the tracer and the free surface height, 139 two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}. 140 141 The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco}) 142 is the centred (\textit{now}) \textit{effective} ocean velocity, 143 \ie the \textit{eulerian} velocity (see \autoref{chap:DYN}) plus 144 the eddy induced velocity (\textit{eiv}) and/or 145 the mixed layer eddy induced velocity (\textit{eiv}) when 146 those parameterisations are used (see \autoref{chap:LDF}). 138 two quantities that are not correlated \citep{Roullet_Madec_JGR00, Griffies_al_MWR01, Campin2004}. 139 140 The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco}) is 141 the centred (\textit{now}) \textit{effective} ocean velocity, \ie the \textit{eulerian} velocity 142 (see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) and/or 143 the mixed layer eddy induced velocity (\textit{eiv}) when those parameterisations are used 144 (see \autoref{chap:LDF}). 147 145 148 146 Several tracer advection scheme are proposed, namely a $2^{nd}$ or $4^{th}$ order centred schemes (CEN), 149 a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), 150 a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL), 151 a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), 152 and a Quadratic Upstream Interpolation for Convective Kinematics with 153 Estimated Streaming Terms scheme (QUICKEST). 154 The choice is made in the \textit{\ngn{namtra\_adv}} namelist, 155 by setting to \forcode{.true.} one of the logicals \textit{ln\_traadv\_xxx}. 156 The corresponding code can be found in the \mdl{traadv\_xxx} module, 157 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 158 By default (\ie in the reference namelist, \ngn{namelist\_ref}), all the logicals are set to \forcode{.false.}. 159 If the user does not select an advection scheme in the configuration namelist (\ngn{namelist\_cfg}), 147 a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), a Monotone Upstream Scheme for 148 Conservative Laws scheme (MUSCL), a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), 149 and a Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms scheme (QUICKEST). 150 The choice is made in the \ngn{namtra\_adv} namelist, by setting to \forcode{.true.} one of 151 the logicals \textit{ln\_traadv\_xxx}. 152 The corresponding code can be found in the \textit{traadv\_xxx.F90} module, where 153 \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 154 By default (\ie in the reference namelist, \textit{namelist\_ref}), all the logicals are set to \forcode{.false.}. 155 If the user does not select an advection scheme in the configuration namelist (\textit{namelist\_cfg}), 160 156 the tracers will \textit{not} be advected! 161 157 … … 163 159 The choosing an advection scheme is a complex matter which depends on the model physics, model resolution, 164 160 type of tracer, as well as the issue of numerical cost. In particular, we note that 165 (1) CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that 166 they do not necessarily need additional diffusion; 167 (2) CEN and UBS are not \textit{positive} schemes 168 \footnote{negative values can appear in an initially strictly positive tracer field which is advected}, 169 implying that false extrema are permitted. 170 Their use is not recommended on passive tracers; 171 (3) It is recommended that the same advection-diffusion scheme is used on both active and passive tracers. 172 Indeed, if a source or sink of a passive tracer depends on an active one, 173 the difference of treatment of active and passive tracers can create very nice-looking frontal structures that 174 are pure numerical artefacts. 161 162 \begin{enumerate} 163 \item 164 CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that 165 they do not necessarily need additional diffusion; 166 \item 167 CEN and UBS are not \textit{positive} schemes 168 \footnote{negative values can appear in an initially strictly positive tracer field which is advected}, 169 implying that false extrema are permitted. 170 Their use is not recommended on passive tracers; 171 \item 172 It is recommended that the same advection-diffusion scheme is used on both active and passive tracers. 173 \end{enumerate} 174 175 Indeed, if a source or sink of a passive tracer depends on an active one, the difference of treatment of active and 176 passive tracers can create very nice-looking frontal structures that are pure numerical artefacts. 175 177 Nevertheless, most of our users set a different treatment on passive and active tracers, 176 178 that's the reason why this possibility is offered. 177 We strongly suggest them to perform a sensitivity experiment using a same treatment to 178 assess the robustness oftheir results.179 We strongly suggest them to perform a sensitivity experiment using a same treatment to assess the robustness of 180 their results. 179 181 180 182 % ------------------------------------------------------------------------------------------------------------- 181 183 % 2nd and 4th order centred schemes 182 184 % ------------------------------------------------------------------------------------------------------------- 183 \subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen} \forcode{= .true.})}185 \subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}~\forcode{= .true.})} 184 186 \label{subsec:TRA_adv_cen} 185 187 186 188 % 2nd order centred scheme 187 189 188 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen} \forcode{= .true.}.190 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~\forcode{= .true.}. 189 191 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 190 192 setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$. … … 196 198 \begin{equation} 197 199 \label{eq:tra_adv_cen2} 198 \tau_u^{cen2} = \overline T ^{i+1/2}200 \tau_u^{cen2} = \overline T ^{i + 1/2} 199 201 \end{equation} 200 202 201 CEN2 is non diffusive (\ie it conserves the tracer variance, $\tau^2 )$but dispersive203 CEN2 is non diffusive (\ie it conserves the tracer variance, $\tau^2$) but dispersive 202 204 (\ie it may create false extrema). 203 205 It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to 204 206 produce a sensible solution. 205 207 The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 206 so $T$ in (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value. 208 so $T$ in (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value. 207 209 208 210 Note that using the CEN2, the overall tracer advection is of second order accuracy since … … 216 218 \begin{equation} 217 219 \label{eq:tra_adv_cen4} 218 \tau_u^{cen4} = \overline{ T - \frac{1}{6}\,\delta_i \left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2}220 \tau_u^{cen4} = \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \, \Big]}^{\,i + 1/2} 219 221 \end{equation} 220 In the vertical direction (\np{nn\_cen\_v} \forcode{= 4}),222 In the vertical direction (\np{nn\_cen\_v}~\forcode{= 4}), 221 223 a $4^{th}$ COMPACT interpolation has been prefered \citep{Demange_PhD2014}. 222 224 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion, 223 spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}. 224 225 spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}. 225 226 226 227 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but … … 249 250 % FCT scheme 250 251 % ------------------------------------------------------------------------------------------------------------- 251 \subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct} \forcode{= .true.})}252 \subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}~\forcode{= .true.})} 252 253 \label{subsec:TRA_adv_tvd} 253 254 254 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct} \forcode{= .true.}.255 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~\forcode{= .true.}. 255 256 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 256 257 setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$. … … 263 264 \label{eq:tra_adv_fct} 264 265 \begin{split} 265 \tau_u^{ups} &=266 \tau_u^{ups} &= 266 267 \begin{cases} 267 T_{i+1} & \text{if $\ u_{i+1/2} < 0$} \hfill\\268 T_i & \text{if $\ u_{i+1/2} \geq 0$} \hfill\\268 T_{i + 1} & \text{if~} u_{i + 1/2} < 0 \\ 269 T_i & \text{if~} u_{i + 1/2} \geq 0 \\ 269 270 \end{cases} 270 \\ \\271 \tau_u^{fct} &=\tau_u^{ups} +c_u \;\left( {\tau_u^{cen} -\tau_u^{ups} } \right)271 \\ 272 \tau_u^{fct} &= \tau_u^{ups} + c_u \, \big( \tau_u^{cen} - \tau_u^{ups} \big) 272 273 \end{split} 273 274 \end{equation} … … 278 279 The one chosen in \NEMO is described in \citet{Zalesak_JCP79}. 279 280 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field. 280 The resulting scheme is quite expensive but \ emph{positive}.281 The resulting scheme is quite expensive but \textit{positive}. 281 282 It can be used on both active and passive tracers. 282 283 A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}. … … 294 295 $\tau_u^{ups}$ is evaluated using the \textit{before} tracer. 295 296 In other words, the advective part of the scheme is time stepped with a leap-frog scheme 296 while a forward scheme is used for the diffusive part. 297 while a forward scheme is used for the diffusive part. 297 298 298 299 % ------------------------------------------------------------------------------------------------------------- 299 300 % MUSCL scheme 300 301 % ------------------------------------------------------------------------------------------------------------- 301 \subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus} \forcode{= .true.})}302 \subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}~\forcode{= .true.})} 302 303 \label{subsec:TRA_adv_mus} 303 304 304 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus} \forcode{= .true.}.305 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~\forcode{= .true.}. 305 306 MUSCL implementation can be found in the \mdl{traadv\_mus} module. 306 307 … … 309 310 two $T$-points (\autoref{fig:adv_scheme}). 310 311 For example, in the $i$-direction : 311 \ [312 \begin{equation} 312 313 % \label{eq:tra_adv_mus} 313 \tau_u^{mus} = \l eft\{314 \begin{aligned}315 &\tau_i &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)316 &\ \widetilde{\partial _i \tau} & \quad \text{if }\;u_{i+1/2} \geqslant 0\\317 &\tau_{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\rdt}{e_{1u} } \right)318 &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0319 \end{aligned}320 \right.321 \ ]322 where $\widetilde{\partial 314 \tau_u^{mus} = \lt\{ 315 \begin{split} 316 \tau_i &+ \frac{1}{2} \lt( 1 - \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 317 \widetilde{\partial_i \tau} & \text{if~} u_{i + 1/2} \geqslant 0 \\ 318 \tau_{i + 1/2} &+ \frac{1}{2} \lt( 1 + \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 319 \widetilde{\partial_{i + 1/2} \tau} & \text{if~} u_{i + 1/2} < 0 320 \end{split} 321 \rt. 322 \end{equation} 323 where $\widetilde{\partial_i \tau}$ is the slope of the tracer on which a limitation is imposed to 323 324 ensure the \textit{positive} character of the scheme. 324 325 325 The time stepping is performed using a forward scheme, 326 that is the \textit{before} tracer field is used toevaluate $\tau_u^{mus}$.326 The time stepping is performed using a forward scheme, that is the \textit{before} tracer field is used to 327 evaluate $\tau_u^{mus}$. 327 328 328 329 For an ocean grid point adjacent to land and where the ocean velocity is directed toward land, … … 330 331 This choice ensure the \textit{positive} character of the scheme. 331 332 In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using upstream fluxes 332 (\np{ln\_mus\_ups} \forcode{= .true.}).333 (\np{ln\_mus\_ups}~\forcode{= .true.}). 333 334 334 335 % ------------------------------------------------------------------------------------------------------------- 335 336 % UBS scheme 336 337 % ------------------------------------------------------------------------------------------------------------- 337 \subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs} \forcode{= .true.})}338 \subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}~\forcode{= .true.})} 338 339 \label{subsec:TRA_adv_ubs} 339 340 340 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs} \forcode{= .true.}.341 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~\forcode{= .true.}. 341 342 UBS implementation can be found in the \mdl{traadv\_mus} module. 342 343 … … 347 348 \begin{equation} 348 349 \label{eq:tra_adv_ubs} 349 \tau_u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{ 350 \begin{aligned} 351 &\tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ 352 &\tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 353 \end{aligned} 354 \right. 350 \tau_u^{ubs} = \overline T ^{i + 1/2} - \frac{1}{6} 351 \begin{cases} 352 \tau"_i & \text{if~} u_{i + 1/2} \geqslant 0 \\ 353 \tau"_{i + 1} & \text{if~} u_{i + 1/2} < 0 354 \end{cases} 355 \quad 356 \text{where~} \tau"_i = \delta_i \lt[ \delta_{i + 1/2} [\tau] \rt] 355 357 \end{equation} 356 where $\tau "_i =\delta_i \left[ {\delta_{i+1/2} \left[ \tau \right]} \right]$. 357 358 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error 358 359 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error 359 360 \citep{Shchepetkin_McWilliams_OM05}. 360 361 The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}. 361 362 It is a relatively good compromise between accuracy and smoothness. 362 Nevertheless the scheme is not \ emph{positive}, meaning that false extrema are permitted,363 Nevertheless the scheme is not \textit{positive}, meaning that false extrema are permitted, 363 364 but the amplitude of such are significantly reduced over the centred second or fourth order method. 364 365 Therefore it is not recommended that it should be applied to a passive tracer that requires positivity. … … 368 369 \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}. 369 370 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme 370 (\np{nn\_cen\_v} \forcode{= 2 or 4}).371 (\np{nn\_cen\_v}~\forcode{= 2 or 4}). 371 372 372 373 For stability reasons (see \autoref{chap:STP}), the first term in \autoref{eq:tra_adv_ubs} … … 382 383 383 384 Note that it is straightforward to rewrite \autoref{eq:tra_adv_ubs} as follows: 384 \ [385 \begin{gather} 385 386 \label{eq:traadv_ubs2} 386 \tau_u^{ubs} = \tau_u^{cen4} + \frac{1}{12} \left\{ 387 \begin{aligned} 388 & + \tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ 389 & - \tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 390 \end{aligned} 391 \right. 392 \] 393 or equivalently 394 \[ 387 \tau_u^{ubs} = \tau_u^{cen4} + \frac{1}{12} 388 \begin{cases} 389 + \tau"_i & \text{if} \ u_{i + 1/2} \geqslant 0 \\ 390 - \tau"_{i + 1} & \text{if} \ u_{i + 1/2} < 0 391 \end{cases} 392 \intertext{or equivalently} 395 393 % \label{eq:traadv_ubs2b} 396 u_{i+1/2} \ \tau_u^{ubs} 397 =u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta_i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2} 398 - \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] 399 \] 394 u_{i + 1/2} \ \tau_u^{ubs} = u_{i + 1/2} \, \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \Big]}^{\,i + 1/2} 395 - \frac{1}{2} |u|_{i + 1/2} \, \frac{1}{6} \, \delta_{i + 1/2} [\tau"_i] \nonumber 396 \end{gather} 400 397 401 398 \autoref{eq:traadv_ubs2} has several advantages. … … 403 400 an upstream-biased diffusion term is added. 404 401 Secondly, this emphasises that the $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has to 405 be evaluated at the \ emph{now} time step using \autoref{eq:tra_adv_ubs}.402 be evaluated at the \textit{now} time step using \autoref{eq:tra_adv_ubs}. 406 403 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy coefficient which 407 is simply proportional to the velocity: 408 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. 404 is simply proportional to the velocity: $A_u^{lm} = \frac{1}{12} \, {e_{1u}}^3 \, |u|$. 409 405 Note the current version of NEMO uses the computationally more efficient formulation \autoref{eq:tra_adv_ubs}. 410 406 … … 412 408 % QCK scheme 413 409 % ------------------------------------------------------------------------------------------------------------- 414 \subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck} \forcode{= .true.})}410 \subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}~\forcode{= .true.})} 415 411 \label{subsec:TRA_adv_qck} 416 412 417 413 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme 418 proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck} \forcode{= .true.}.414 proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}~\forcode{= .true.}. 419 415 QUICKEST implementation can be found in the \mdl{traadv\_qck} module. 420 416 … … 422 418 \citep{Leonard1991}. 423 419 It has been implemented in NEMO by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. 424 The resulting scheme is quite expensive but \ emph{positive}.420 The resulting scheme is quite expensive but \textit{positive}. 425 421 It can be used on both active and passive tracers. 426 422 However, the intrinsic diffusion of QCK makes its use risky in the vertical direction where … … 431 427 432 428 %%%gmcomment : Cross term are missing in the current implementation.... 433 434 429 435 430 % ================================================================ … … 458 453 except for the pure vertical component that appears when a rotation tensor is used. 459 454 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}). 460 When \np{ln\_traldf\_msc} \forcode{= .true.}, a Method of Stabilizing Correction is used in which455 When \np{ln\_traldf\_msc}~\forcode{= .true.}, a Method of Stabilizing Correction is used in which 461 456 the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}. 462 457 … … 464 459 % Type of operator 465 460 % ------------------------------------------------------------------------------------------------------------- 466 \subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})] 467 {Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 461 \subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})]{Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 468 462 \label{subsec:TRA_ldf_op} 469 463 470 464 Three operator options are proposed and, one and only one of them must be selected: 465 471 466 \begin{description} 472 \item[\np{ln\_traldf\_NONE} \forcode{= .true.}:]467 \item[\np{ln\_traldf\_NONE}~\forcode{= .true.}:] 473 468 no operator selected, the lateral diffusive tendency will not be applied to the tracer equation. 474 469 This option can be used when the selected advection scheme is diffusive enough (MUSCL scheme for example). 475 \item[\np{ln\_traldf\_lap} \forcode{= .true.}:]470 \item[\np{ln\_traldf\_lap}~\forcode{= .true.}:] 476 471 a laplacian operator is selected. 477 This harmonic operator takes the following expression: $\mathpzc{L}(T) =\nabla \cdot A_{ht}\;\nabla T $,472 This harmonic operator takes the following expression: $\mathpzc{L}(T) = \nabla \cdot A_{ht} \; \nabla T $, 478 473 where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}), 479 474 and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}). 480 \item[\np{ln\_traldf\_blp} \forcode{= .true.}]:475 \item[\np{ln\_traldf\_blp}~\forcode{= .true.}]: 481 476 a bilaplacian operator is selected. 482 477 This biharmonic operator takes the following expression: 483 $\mathpzc{B} =- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$478 $\mathpzc{B} = - \mathpzc{L}(\mathpzc{L}(T)) = - \nabla \cdot b \nabla (\nabla \cdot b \nabla T)$ 484 479 where the gradient operats along the selected direction, 485 and $b^2 =B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$(see \autoref{chap:LDF}).480 and $b^2 = B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}). 486 481 In the code, the bilaplacian operator is obtained by calling the laplacian twice. 487 482 \end{description} … … 495 490 whereas the laplacian damping time scales only like $\lambda^{-2}$. 496 491 497 498 492 % ------------------------------------------------------------------------------------------------------------- 499 493 % Direction of action 500 494 % ------------------------------------------------------------------------------------------------------------- 501 \subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})] 502 {Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 495 \subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})]{Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 503 496 \label{subsec:TRA_ldf_dir} 504 497 505 498 The choice of a direction of action determines the form of operator used. 506 499 The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane when 507 iso-level option is used (\np{ln\_traldf\_lev} \forcode{= .true.}) or508 when a horizontal (\ie geopotential) operator is demanded in \ zstar-coordinate500 iso-level option is used (\np{ln\_traldf\_lev}~\forcode{= .true.}) or 501 when a horizontal (\ie geopotential) operator is demanded in \textit{z}-coordinate 509 502 (\np{ln\_traldf\_hor} and \np{ln\_zco} equal \forcode{.true.}). 510 503 The associated code can be found in the \mdl{traldf\_lap\_blp} module. … … 521 514 522 515 The resulting discret form of the three operators (one iso-level and two rotated one) is given in 523 the next two sub-sections. 524 516 the next two sub-sections. 525 517 526 518 % ------------------------------------------------------------------------------------------------------------- 527 519 % iso-level operator 528 520 % ------------------------------------------------------------------------------------------------------------- 529 \subsection{Iso-level (bi -)laplacian operator ( \protect\np{ln\_traldf\_iso}) }521 \subsection{Iso-level (bi -)laplacian operator ( \protect\np{ln\_traldf\_iso}) } 530 522 \label{subsec:TRA_ldf_lev} 531 523 … … 533 525 \begin{equation} 534 526 \label{eq:tra_ldf_lap} 535 D_t^{lT} =\frac{1}{b_t} \left( \; 536 \delta_{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta_{i+1/2} [T] \right] 537 + \delta_{j}\left[ A_v^{lT} \; \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta_{j+1/2} [T] \right] \;\right) 527 D_t^{lT} = \frac{1}{b_t} \Bigg( \delta_{i} \lt[ A_u^{lT} \; \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [T] \rt] 528 + \delta_{j} \lt[ A_v^{lT} \; \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [T] \rt] \Bigg) 538 529 \end{equation} 539 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells and530 where $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells and 540 531 where zero diffusive fluxes is assumed across solid boundaries, 541 532 first (and third in bilaplacian case) horizontal tracer derivative are masked. 542 533 It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module. 543 534 The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} in order to 544 compute the iso-level bilaplacian operator. 545 546 It is a \ emph{horizontal} operator (\ie acting along geopotential surfaces) in535 compute the iso-level bilaplacian operator. 536 537 It is a \textit{horizontal} operator (\ie acting along geopotential surfaces) in 547 538 the $z$-coordinate with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 548 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp} \forcode{= .true.},549 we have \np{ln\_traldf\_lev} \forcode{ = .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}\forcode{= .true.}.539 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~\forcode{= .true.}, 540 we have \np{ln\_traldf\_lev}~\forcode{= .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~\forcode{= .true.}. 550 541 In both cases, it significantly contributes to diapycnal mixing. 551 542 It is therefore never recommended, even when using it in the bilaplacian case. 552 543 553 Note that in the partial step $z$-coordinate (\np{ln\_zps} \forcode{= .true.}),544 Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}), 554 545 tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom. 555 546 In this case, horizontal derivatives in (\autoref{eq:tra_ldf_lap}) at the bottom level require a specific treatment. 556 547 They are calculated in the \mdl{zpshde} module, described in \autoref{sec:TRA_zpshde}. 557 548 558 559 549 % ------------------------------------------------------------------------------------------------------------- 560 550 % Rotated laplacian operator 561 551 % ------------------------------------------------------------------------------------------------------------- 562 \subsection{Standard and triad (bi -)laplacian operator}552 \subsection{Standard and triad (bi -)laplacian operator} 563 553 \label{subsec:TRA_ldf_iso_triad} 564 554 565 %&& Standard rotated (bi -)laplacian operator555 %&& Standard rotated (bi -)laplacian operator 566 556 %&& ---------------------------------------------- 567 \subsubsection{Standard rotated (bi -)laplacian operator (\protect\mdl{traldf\_iso})}557 \subsubsection{Standard rotated (bi -)laplacian operator (\protect\mdl{traldf\_iso})} 568 558 \label{subsec:TRA_ldf_iso} 569 570 559 The general form of the second order lateral tracer subgrid scale physics (\autoref{eq:PE_zdf}) 571 takes the following semi -discrete space form in $z$- and $s$-coordinates:560 takes the following semi -discrete space form in $z$- and $s$-coordinates: 572 561 \begin{equation} 573 562 \label{eq:tra_ldf_iso} 574 563 \begin{split} 575 D_T^{lT} = \frac{1}{b_t} & \left\{ \,\;\delta_i \left[ A_u^{lT} \left( 576 \frac{e_{2u}\,e_{3u}}{e_{1u}} \,\delta_{i+1/2}[T] 577 - e_{2u}\;r_{1u} \,\overline{\overline{ \delta_{k+1/2}[T] }}^{\,i+1/2,k} 578 \right) \right] \right. \\ 579 & +\delta_j \left[ A_v^{lT} \left( 580 \frac{e_{1v}\,e_{3v}}{e_{2v}} \,\delta_{j+1/2} [T] 581 - e_{1v}\,r_{2v} \,\overline{\overline{ \delta_{k+1/2} [T] }}^{\,j+1/2,k} 582 \right) \right] \\ 583 & +\delta_k \left[ A_w^{lT} \left( 584 -\;e_{2w}\,r_{1w} \,\overline{\overline{ \delta_{i+1/2} [T] }}^{\,i,k+1/2} 585 \right. \right. \\ 586 & \qquad \qquad \quad 587 - e_{1w}\,r_{2w} \,\overline{\overline{ \delta_{j+1/2} [T] }}^{\,j,k+1/2} \\ 588 & \left. {\left. { \qquad \qquad \ \ \ \left. { 589 +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right) 590 \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} 564 D_T^{lT} = \frac{1}{b_t} \Bigg[ \quad &\delta_i A_u^{lT} \lt( \frac{e_{2u} e_{3u}}{e_{1u}} \, \delta_{i + 1/2} [T] 565 - e_{2u} r_{1u} \, \overline{\overline{\delta_{k + 1/2} [T]}}^{\,i + 1/2,k} \rt) \Bigg. \\ 566 + &\delta_j A_v^{lT} \lt( \frac{e_{1v} e_{3v}}{e_{2v}} \, \delta_{j + 1/2} [T] 567 - e_{1v} r_{2v} \, \overline{\overline{\delta_{k + 1/2} [T]}}^{\,j + 1/2,k} \rt) \\ 568 + &\delta_k A_w^{lT} \lt( \frac{e_{1w} e_{2w}}{e_{3w}} (r_{1w}^2 + r_{2w}^2) \, \delta_{k + 1/2} [T] \rt. \\ 569 & \qquad \quad \Bigg. \lt. - e_{2w} r_{1w} \, \overline{\overline{\delta_{i + 1/2} [T]}}^{\,i,k + 1/2} 570 - e_{1w} r_{2w} \, \overline{\overline{\delta_{j + 1/2} [T]}}^{\,j,k + 1/2} \rt) \Bigg] 591 571 \end{split} 592 572 \end{equation} 593 where $b_t $=$e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells,573 where $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells, 594 574 $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and 595 575 the surface along which the diffusion operator acts (\ie horizontal or iso-neutral surfaces). 596 It is thus used when, in addition to \np{ln\_traldf\_lap} \forcode{= .true.},597 we have \np{ln\_traldf\_iso} \forcode{= .true.},598 or both \np{ln\_traldf\_hor} \forcode{ = .true.} and \np{ln\_zco}\forcode{= .true.}.576 It is thus used when, in addition to \np{ln\_traldf\_lap}~\forcode{= .true.}, 577 we have \np{ln\_traldf\_iso}~\forcode{= .true.}, 578 or both \np{ln\_traldf\_hor}~\forcode{= .true.} and \np{ln\_zco}~\forcode{= .true.}. 599 579 The way these slopes are evaluated is given in \autoref{sec:LDF_slp}. 600 580 At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using 601 the mask technique (see \autoref{sec:LBC_coast}). 581 the mask technique (see \autoref{sec:LBC_coast}). 602 582 603 583 The operator in \autoref{eq:tra_ldf_iso} involves both lateral and vertical derivatives. … … 606 586 For computer efficiency reasons, this term is not computed in the \mdl{traldf\_iso} module, 607 587 but in the \mdl{trazdf} module where, if iso-neutral mixing is used, 608 the vertical mixing coefficient is simply increased by 609 $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$. 588 the vertical mixing coefficient is simply increased by $\frac{e_{1w} e_{2w}}{e_{3w}}(r_{1w}^2 + r_{2w}^2)$. 610 589 611 590 This formulation conserves the tracer but does not ensure the decrease of the tracer variance. 612 591 Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without 613 any additional background horizontal diffusion \citep{Guilyardi_al_CD01}. 614 615 Note that in the partial step $z$-coordinate (\np{ln\_zps} \forcode{= .true.}),592 any additional background horizontal diffusion \citep{Guilyardi_al_CD01}. 593 594 Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}), 616 595 the horizontal derivatives at the bottom level in \autoref{eq:tra_ldf_iso} require a specific treatment. 617 596 They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}. 618 597 619 %&& Triad rotated (bi -)laplacian operator598 %&& Triad rotated (bi -)laplacian operator 620 599 %&& ------------------------------------------- 621 \subsubsection{Triad rotated (bi -)laplacian operator (\protect\np{ln\_traldf\_triad})}600 \subsubsection{Triad rotated (bi -)laplacian operator (\protect\np{ln\_traldf\_triad})} 622 601 \label{subsec:TRA_ldf_triad} 623 602 624 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad} \forcode{= .true.}; see \autoref{apdx:triad})603 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}~\forcode{= .true.}; see \autoref{apdx:triad}) 625 604 626 605 An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 627 is also available in \NEMO (\np{ln\_traldf\_grif} \forcode{= .true.}).606 is also available in \NEMO (\np{ln\_traldf\_grif}~\forcode{= .true.}). 628 607 A complete description of the algorithm is given in \autoref{apdx:triad}. 629 608 … … 635 614 It requires an additional assumption on boundary conditions: 636 615 first and third derivative terms normal to the coast, 637 normal to the bottom and normal to the surface are set to zero. 616 normal to the bottom and normal to the surface are set to zero. 638 617 639 618 %&& Option for the rotated operators … … 642 621 \label{subsec:TRA_ldf_options} 643 622 644 \np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators) 645 646 \np{rn\_slpmax} = slope limit (both operators) 647 648 \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 649 650 \np{rn\_sw\_triad} =1 switching triad; 651 =0 all 4 triads used (triad only) 652 653 \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only) 623 \begin{itemize} 624 \item \np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators) 625 \item \np{rn\_slpmax} = slope limit (both operators) 626 \item \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 627 \item \np{rn\_sw\_triad} $= 1$ switching triad; $= 0$ all 4 triads used (triad only) 628 \item \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only) 629 \end{itemize} 654 630 655 631 % ================================================================ … … 666 642 The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates, 667 643 and is based on a laplacian operator. 668 The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the following semi -discrete space form:669 \ [644 The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the following semi -discrete space form: 645 \begin{gather*} 670 646 % \label{eq:tra_zdf} 671 \begin{split} 672 D^{vT}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}} \delta_{k+1/2}[T] \;\right] \\ 673 D^{vS}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vS}_w}{e_{3w}} \delta_{k+1/2}[S] \;\right] 674 \end{split} 675 \] 647 D^{vT}_T = \frac{1}{e_{3t}} \, \delta_k \lt[ \, \frac{A^{vT}_w}{e_{3w}} \delta_{k + 1/2}[T] \, \rt] \\ 648 D^{vS}_T = \frac{1}{e_{3t}} \; \delta_k \lt[ \, \frac{A^{vS}_w}{e_{3w}} \delta_{k + 1/2}[S] \, \rt] 649 \end{gather*} 676 650 where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity coefficients on temperature and salinity, 677 651 respectively. 678 Generally, $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is parameterised (\ie \key{zdfddm} is defined). 652 Generally, $A_w^{vT} = A_w^{vS}$ except when double diffusive mixing is parameterised 653 (\ie \key{zdfddm} is defined). 679 654 The way these coefficients are evaluated is given in \autoref{chap:ZDF} (ZDF). 680 655 Furthermore, when iso-neutral mixing is used, both mixing coefficients are increased by 681 $\frac{e_{1w} \,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ to account for682 the vertical second derivative of\autoref{eq:tra_ldf_iso}.656 $\frac{e_{1w} e_{2w}}{e_{3w} }({r_{1w}^2 + r_{2w}^2})$ to account for the vertical second derivative of 657 \autoref{eq:tra_ldf_iso}. 683 658 684 659 At the surface and bottom boundaries, the turbulent fluxes of heat and salt must be specified. 685 660 At the surface they are prescribed from the surface forcing and added in a dedicated routine 686 661 (see \autoref{subsec:TRA_sbc}), whilst at the bottom they are set to zero for heat and salt unless 687 a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}). 662 a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}). 688 663 689 664 The large eddy coefficient found in the mixed layer together with high vertical resolution implies that 690 in the case of explicit time stepping (\np{ln\_zdfexp} \forcode{= .true.})665 in the case of explicit time stepping (\np{ln\_zdfexp}~\forcode{= .true.}) 691 666 there would be too restrictive a constraint on the time step. 692 667 Therefore, the default implicit time stepping is preferred for the vertical diffusion since 693 668 it overcomes the stability constraint. 694 A forward time differencing scheme (\np{ln\_zdfexp} \forcode{= .true.}) using669 A forward time differencing scheme (\np{ln\_zdfexp}~\forcode{= .true.}) using 695 670 a time splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative. 696 Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 671 Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 697 672 698 673 % ================================================================ … … 712 687 This has been found to enhance readability of the code. 713 688 The two formulations are completely equivalent; 714 the forcing terms in trasbc are the surface fluxes divided by the thickness of the top model layer. 689 the forcing terms in trasbc are the surface fluxes divided by the thickness of the top model layer. 715 690 716 691 Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components … … 724 699 The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following forcing fields (used on tracers): 725 700 726 $\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 727 (\ie the difference between the total surface heat flux and the fraction of the short wave flux that 728 penetrates into the water column, see \autoref{subsec:TRA_qsr}) 729 plus the heat content associated with of the mass exchange with the atmosphere and lands. 730 731 $\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 732 733 $\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and 734 possibly with the sea-ice and ice-shelves. 735 736 $\bullet$ \textit{rnf}, the mass flux associated with runoff 737 (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 738 739 $\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, 740 (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied). 701 \begin{itemize} 702 \item 703 $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 704 (\ie the difference between the total surface heat flux and the fraction of the short wave flux that 705 penetrates into the water column, see \autoref{subsec:TRA_qsr}) 706 plus the heat content associated with of the mass exchange with the atmosphere and lands. 707 \item 708 $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 709 \item 710 \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and 711 possibly with the sea-ice and ice-shelves. 712 \item 713 \textit{rnf}, the mass flux associated with runoff 714 (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 715 \item 716 \textit{fwfisf}, the mass flux associated with ice shelf melt, 717 (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied). 718 \end{itemize} 741 719 742 720 The surface boundary condition on temperature and salinity is applied as follows: 743 721 \begin{equation} 744 722 \label{eq:tra_sbc} 745 \begin{aligned} 746 &F^T = \frac{ 1 }{\rho_o \;C_p \,\left. e_{3t} \right|_{k=1} } &\overline{ Q_{ns} }^t & \\ 747 & F^S =\frac{ 1 }{\rho_o \, \left. e_{3t} \right|_{k=1} } &\overline{ \textit{sfx} }^t & \\ 748 \end{aligned} 749 \end{equation} 750 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps ($t-\rdt/2$ and $t+\rdt/2$). 723 \begin{alignedat}{2} 724 F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{Q_{ns} }^t \\ 725 F^S &= &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{\textit{sfx}}^t 726 \end{alignedat} 727 \end{equation} 728 where $\overline x^t$ means that $x$ is averaged over two consecutive time steps 729 ($t - \rdt / 2$ and $t + \rdt / 2$). 751 730 Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:STP}). 752 731 753 In the linear free surface case (\np{ln\_linssh} \forcode{= .true.}), an additional term has to be added on732 In the linear free surface case (\np{ln\_linssh}~\forcode{= .true.}), an additional term has to be added on 754 733 both temperature and salinity. 755 734 On temperature, this term remove the heat content associated with mass exchange that has been added to $Q_{ns}$. … … 759 738 \begin{equation} 760 739 \label{eq:tra_sbc_lin} 761 \begin{aligned} 762 &F^T = \frac{ 1 }{\rho_o \;C_p \,\left. e_{3t} \right|_{k=1} } 763 &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^t & \\ 764 % 765 & F^S =\frac{ 1 }{\rho_o \,\left. e_{3t} \right|_{k=1} } 766 &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1} \right) }^t & \\ 767 \end{aligned} 740 \begin{alignedat}{2} 741 F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 742 &\overline{(Q_{ns} - C_p \, \textit{emp} \lt. T \rt|_{k = 1})}^t \\ 743 F^S &= &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 744 &\overline{(\textit{sfx} - \textit{emp} \lt. S \rt|_{k = 1})}^t 745 \end{alignedat} 768 746 \end{equation} 769 747 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface. … … 783 761 784 762 Options are defined through the \ngn{namtra\_qsr} namelist variables. 785 When the penetrative solar radiation option is used (\np{ln\_flxqsr} \forcode{= .true.}),763 When the penetrative solar radiation option is used (\np{ln\_flxqsr}~\forcode{= .true.}), 786 764 the solar radiation penetrates the top few tens of meters of the ocean. 787 If it is not used (\np{ln\_flxqsr} \forcode{= .false.}) all the heat flux is absorbed in the first ocean level.765 If it is not used (\np{ln\_flxqsr}~\forcode{= .false.}) all the heat flux is absorbed in the first ocean level. 788 766 Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq:PE_tra_T} and 789 767 the surface boundary condition is modified to take into account only the non-penetrative part of the surface … … 791 769 \begin{equation} 792 770 \label{eq:PE_qsr} 793 \begin{ split}794 \ frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k}\\795 Q_{ns} &= Q_\text{Total} - Q_{sr}796 \end{ split}771 \begin{gathered} 772 \pd[T]{t} = \ldots + \frac{1}{\rho_o \, C_p \, e_3} \; \pd[I]{k} \\ 773 Q_{ns} = Q_\text{Total} - Q_{sr} 774 \end{gathered} 797 775 \end{equation} 798 776 where $Q_{sr}$ is the penetrative part of the surface heat flux (\ie the shortwave radiation) and 799 $I$ is the downward irradiance ($\l eft. I \right|_{z=\eta}=Q_{sr}$).777 $I$ is the downward irradiance ($\lt. I \rt|_{z = \eta} = Q_{sr}$). 800 778 The additional term in \autoref{eq:PE_qsr} is discretized as follows: 801 779 \begin{equation} 802 780 \label{eq:tra_qsr} 803 \frac{1}{\rho_o \, C_p \,e_3} \; \frac{\partial I}{\partial k} \equiv \frac{1}{\rho_o\, C_p\, e_{3t}} \delta_k \left[ I_w \right]781 \frac{1}{\rho_o \, C_p \, e_3} \, \pd[I]{k} \equiv \frac{1}{\rho_o \, C_p \, e_{3t}} \delta_k [I_w] 804 782 \end{equation} 805 783 … … 810 788 (specified through namelist parameter \np{rn\_abs}). 811 789 It is assumed to penetrate the ocean with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 812 of a few tens of centimetres (typically $\xi_0 =0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist).790 of a few tens of centimetres (typically $\xi_0 = 0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist). 813 791 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy propagates to 814 792 larger depths where it contributes to local heating. 815 793 The way this second part of the solar energy penetrates into the ocean depends on which formulation is chosen. 816 In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd} \forcode{= .true.})794 In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}~\forcode{= .true.}) 817 795 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 818 796 leading to the following expression \citep{Paulson1977}: 819 797 \[ 820 798 % \label{eq:traqsr_iradiance} 821 I(z) = Q_{sr} \l eft[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right]799 I(z) = Q_{sr} \lt[ Re^{- z / \xi_0} + (1 - R) e^{- z / \xi_1} \rt] 822 800 \] 823 801 where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths. 824 802 It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter. 825 The set of default values ($\xi_0 $, $\xi_1$, $R$) corresponds to a Type I water in Jerlov's (1968) classification803 The set of default values ($\xi_0, \xi_1, R$) corresponds to a Type I water in Jerlov's (1968) classification 826 804 (oligotrophic waters). 827 805 … … 840 818 reproduces quite closely the light penetration profiles predicted by the full spectal model, 841 819 but with much greater computational efficiency. 842 The 2-bands formulation does not reproduce the full model very well. 843 844 The RGB formulation is used when \np{ln\_qsr\_rgb} \forcode{= .true.}.820 The 2-bands formulation does not reproduce the full model very well. 821 822 The RGB formulation is used when \np{ln\_qsr\_rgb}~\forcode{= .true.}. 845 823 The RGB attenuation coefficients (\ie the inverses of the extinction length scales) are tabulated over 846 824 61 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L 847 825 (see the routine \rou{trc\_oce\_rgb} in \mdl{trc\_oce} module). 848 826 Four types of chlorophyll can be chosen in the RGB formulation: 849 \begin{description} 850 \item[\np{nn\_chdta}\forcode{ = 0}] 827 828 \begin{description} 829 \item[\np{nn\_chdta}~\forcode{= 0}] 851 830 a constant 0.05 g.Chl/L value everywhere ; 852 \item[\np{nn\_chdta} \forcode{= 1}]831 \item[\np{nn\_chdta}~\forcode{= 1}] 853 832 an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in 854 833 the vertical direction; 855 \item[\np{nn\_chdta} \forcode{= 2}]834 \item[\np{nn\_chdta}~\forcode{= 2}] 856 835 same as previous case except that a vertical profile of chlorophyl is used. 857 836 Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value; 858 \item[\np{ln\_qsr\_bio} \forcode{= .true.}]837 \item[\np{ln\_qsr\_bio}~\forcode{= .true.}] 859 838 simulated time varying chlorophyll by TOP biogeochemical model. 860 839 In this case, the RGB formulation is used to calculate both the phytoplankton light limitation in 861 PISCES or LOBSTER and the oceanic heating rate. 840 PISCES or LOBSTER and the oceanic heating rate. 862 841 \end{description} 842 863 843 The trend in \autoref{eq:tra_qsr} associated with the penetration of the solar radiation is added to 864 844 the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. … … 871 851 Finally, note that when the ocean is shallow ($<$ 200~m), part of the solar radiation can reach the ocean floor. 872 852 In this case, we have chosen that all remaining radiation is absorbed in the last ocean level 873 (\ie $I$ is masked). 853 (\ie $I$ is masked). 874 854 875 855 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 876 856 \begin{figure}[!t] 877 857 \begin{center} 878 \includegraphics[ width=1.0\textwidth]{Fig_TRA_Irradiance}858 \includegraphics[]{Fig_TRA_Irradiance} 879 859 \caption{ 880 860 \protect\label{fig:traqsr_irradiance} … … 903 883 \begin{figure}[!t] 904 884 \begin{center} 905 \includegraphics[ width=1.0\textwidth]{Fig_TRA_geoth}885 \includegraphics[]{Fig_TRA_geoth} 906 886 \caption{ 907 887 \protect\label{fig:geothermal} … … 917 897 This is the default option in \NEMO, and it is implemented using the masking technique. 918 898 However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling. 919 This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1\;W/m^2$ \citep{Stein_Stein_Nat92}),899 This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{Stein_Stein_Nat92}), 920 900 but it warms systematically the ocean and acts on the densest water masses. 921 901 Taking this flux into account in a global ocean model increases the deepest overturning cell 922 (\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}.902 (\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}. 923 903 924 904 Options are defined through the \ngn{namtra\_bbc} namelist variables. … … 939 919 %-------------------------------------------------------------------------------------------------------------- 940 920 941 Options are defined through the 921 Options are defined through the \ngn{nambbl} namelist variables. 942 922 In a $z$-coordinate configuration, the bottom topography is represented by a series of discrete steps. 943 923 This is not adequate to represent gravity driven downslope flows. … … 951 931 sometimes over a thickness much larger than the thickness of the observed gravity plume. 952 932 A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of 953 a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved. 933 a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved. 954 934 955 935 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{Beckmann_Doscher1997}, … … 964 944 % Diffusive BBL 965 945 % ------------------------------------------------------------------------------------------------------------- 966 \subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf} \forcode{= 1})}946 \subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}~\forcode{= 1})} 967 947 \label{subsec:TRA_bbl_diff} 968 948 … … 971 951 \[ 972 952 % \label{eq:tra_bbl_diff} 973 {\rm {\bf F}}_\sigma=A_l^\sigma \;\nabla_\sigma T953 \vect F_\sigma = A_l^\sigma \, \nabla_\sigma T 974 954 \] 975 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 976 and$A_l^\sigma$ the lateral diffusivity in the BBL.955 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, and 956 $A_l^\sigma$ the lateral diffusivity in the BBL. 977 957 Following \citet{Beckmann_Doscher1997}, the latter is prescribed with a spatial dependence, 978 958 \ie in the conditional form 979 959 \begin{equation} 980 960 \label{eq:tra_bbl_coef} 981 A_l^\sigma (i,j,t) =\left\{ {982 \begin{ array}{l}983 A_{bbl} \quad \quad \mbox{if} \quad \nabla_\sigma \rho \cdot \nabla H<0 \\\\984 0\quad \quad \;\,\mbox{otherwise}\\985 \end{array}}986 \right.987 \end{equation} 961 A_l^\sigma (i,j,t) = 962 \begin{cases} 963 A_{bbl} & \text{if~} \nabla_\sigma \rho \cdot \nabla H < 0 \\ 964 \\ 965 0 & \text{otherwise} \\ 966 \end{cases} 967 \end{equation} 988 968 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist parameter \np{rn\_ahtbbl} and 989 969 usually set to a value much larger than the one used for lateral mixing in the open ocean. … … 995 975 \[ 996 976 % \label{eq:tra_bbl_Drho} 997 \nabla_\sigma \rho / \rho = \alpha \, \nabla_\sigma T + \beta \,\nabla_\sigma S977 \nabla_\sigma \rho / \rho = \alpha \, \nabla_\sigma T + \beta \, \nabla_\sigma S 998 978 \] 999 where $\rho$, $\alpha$ and $\beta$ are functions of $\overline {T}^\sigma$,1000 $\overline {S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature, salinity and depth, respectively.979 where $\rho$, $\alpha$ and $\beta$ are functions of $\overline T^\sigma$, $\overline S^\sigma$ and 980 $\overline H^\sigma$, the along bottom mean temperature, salinity and depth, respectively. 1001 981 1002 982 % ------------------------------------------------------------------------------------------------------------- 1003 983 % Advective BBL 1004 984 % ------------------------------------------------------------------------------------------------------------- 1005 \subsection{Advective bottom boundary layer (\protect\np{nn\_bbl\_adv} \forcode{= 1..2})}985 \subsection{Advective bottom boundary layer (\protect\np{nn\_bbl\_adv}~\forcode{= 1..2})} 1006 986 \label{subsec:TRA_bbl_adv} 1007 987 1008 %\sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following 1009 %if this is not what is meant then "downwards sloping flow" is also a possibility"} 988 %\sgacomment{ 989 % "downsloping flow" has been replaced by "downslope flow" in the following 990 % if this is not what is meant then "downwards sloping flow" is also a possibility" 991 %} 1010 992 1011 993 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1012 994 \begin{figure}[!t] 1013 995 \begin{center} 1014 \includegraphics[ width=0.7\textwidth]{Fig_BBL_adv}996 \includegraphics[]{Fig_BBL_adv} 1015 997 \caption{ 1016 998 \protect\label{fig:bbl} 1017 999 Advective/diffusive Bottom Boundary Layer. 1018 The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i +1}_{kdnw}$.1000 The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i + 1}_{kdnw}$. 1019 1001 Red arrows indicate the additional overturning circulation due to the advective BBL. 1020 1002 The transport of the downslope flow is defined either as the transport of the bottom ocean cell (black arrow), … … 1026 1008 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1027 1009 1028 1029 1010 %!! nn_bbl_adv = 1 use of the ocean velocity as bbl velocity 1030 1011 %!! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation 1031 %!! \ietransport proportional to the along-slope density gradient1012 %!! i.e. transport proportional to the along-slope density gradient 1032 1013 1033 1014 %%%gmcomment : this section has to be really written 1034 1015 1035 When applying an advective BBL (\np{nn\_bbl\_adv} \forcode{= 1..2}), an overturning circulation is added which1016 When applying an advective BBL (\np{nn\_bbl\_adv}~\forcode{= 1..2}), an overturning circulation is added which 1036 1017 connects two adjacent bottom grid-points only if dense water overlies less dense water on the slope. 1037 The density difference causes dense water to move down the slope. 1038 1039 \np{nn\_bbl\_adv} \forcode{= 1}:1018 The density difference causes dense water to move down the slope. 1019 1020 \np{nn\_bbl\_adv}~\forcode{= 1}: 1040 1021 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step 1041 1022 (see black arrow in \autoref{fig:bbl}) \citep{Beckmann_Doscher1997}. 1042 1023 It is a \textit{conditional advection}, that is, advection is allowed only 1043 if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho \cdot \nabla H<0$) and1044 if the velocity is directed towards greater depth (\ie $\vect {U} \cdot \nabla H>0$).1045 1046 \np{nn\_bbl\_adv} \forcode{= 2}:1024 if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$) and 1025 if the velocity is directed towards greater depth (\ie $\vect U \cdot \nabla H > 0$). 1026 1027 \np{nn\_bbl\_adv}~\forcode{= 2}: 1047 1028 the downslope velocity is chosen to be proportional to $\Delta \rho$, 1048 1029 the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}. 1049 1030 The advection is allowed only if dense water overlies less dense water on the slope 1050 (\ie $\nabla_\sigma \rho \cdot \nabla H<0$).1031 (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$). 1051 1032 For example, the resulting transport of the downslope flow, here in the $i$-direction (\autoref{fig:bbl}), 1052 1033 is simply given by the following expression: 1053 1034 \[ 1054 1035 % \label{eq:bbl_Utr} 1055 u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o} e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right)1036 u^{tr}_{bbl} = \gamma g \frac{\Delta \rho}{\rho_o} e_{1u} \, min ({e_{3u}}_{kup},{e_{3u}}_{kdwn}) 1056 1037 \] 1057 1038 where $\gamma$, expressed in seconds, is the coefficient of proportionality provided as \np{rn\_gambbl}, … … 1062 1043 The possible values for $\gamma$ range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}. 1063 1044 1064 Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl})$ using the upwind scheme.1045 Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using the upwind scheme. 1065 1046 Such a diffusive advective scheme has been chosen to mimic the entrainment between the downslope plume and 1066 1047 the surrounding water at intermediate depths. … … 1071 1052 the downslope flow \autoref{eq:bbl_dw}, the horizontal \autoref{eq:bbl_hor} and 1072 1053 the upward \autoref{eq:bbl_up} return flows as follows: 1073 \begin{align} 1054 \begin{alignat}{3} 1055 \label{eq:bbl_dw} 1074 1056 \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw} 1075 + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} \left( T^{sh}_{kup} - T^{do}_{kdw} \right) \label{eq:bbl_dw}\\1076 %1057 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} &&\lt( T^{sh}_{kup} - T^{do}_{kdw} \rt) \\ 1058 \label{eq:bbl_hor} 1077 1059 \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 1078 + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} \left( T^{do}_{kup} - T^{sh}_{kup} \right) \label{eq:bbl_hor}\\1079 1060 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} &&\lt( T^{do}_{kup} - T^{sh}_{kup} \rt) \\ 1061 % 1080 1062 \intertext{and for $k =kdw-1,\;..., \; kup$ :} 1081 1063 % 1064 \label{eq:bbl_up} 1082 1065 \partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k} 1083 + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}} \left( T^{do}_{k+1} - T^{sh}_{k} \right) \label{eq:bbl_up}1084 \end{align }1085 where $b_t$ is the $T$-cell volume. 1086 1087 Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl})$, is available in the model outputs.1066 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}} &&\lt( T^{do}_{k +1} - T^{sh}_{k} \rt) 1067 \end{alignat} 1068 where $b_t$ is the $T$-cell volume. 1069 1070 Note that the BBL transport, $(u^{tr}_{bbl},v^{tr}_{bbl})$, is available in the model outputs. 1088 1071 It has to be used to compute the effective velocity as well as the effective overturning circulation. 1089 1072 … … 1101 1084 \begin{equation} 1102 1085 \label{eq:tra_dmp} 1103 \begin{ split}1104 \ frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right)\\1105 \ frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right)1106 \end{ split}1086 \begin{gathered} 1087 \pd[T]{t} = \cdots - \gamma (T - T_o) \\ 1088 \pd[S]{t} = \cdots - \gamma (S - S_o) 1089 \end{gathered} 1107 1090 \end{equation} 1108 1091 where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ are given temperature and salinity fields … … 1111 1094 The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true. 1112 1095 It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true in 1113 \ textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set1096 \ngn{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set 1114 1097 (\ie that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread}, 1115 1098 see \autoref{subsec:SBC_fldread}). … … 1128 1111 The second case corresponds to the use of the robust diagnostic method \citep{Sarmiento1982}. 1129 1112 It allows us to find the velocity field consistent with the model dynamics whilst 1130 having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$). 1113 having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$). 1131 1114 1132 1115 The robust diagnostic method is very efficient in preventing temperature drift in intermediate waters but … … 1140 1123 \citep{Madec_al_JPO96}. 1141 1124 1142 \subsection{Generating \ifile{resto} using DMP\_TOOLS} 1143 1144 DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. 1145 Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and 1146 run on the same machine as the NEMO model. 1147 A \ifile{mesh\_mask} file for the model configuration is required as an input. 1148 This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. 1149 The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. 1150 The \ngn{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for 1151 the restoration coefficient. 1152 1153 %--------------------------------------------nam_dmp_create------------------------------------------------- 1154 %\namtools{namelist_dmp} 1155 %------------------------------------------------------------------------------------------------------- 1156 1157 \np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and 1158 should be the same as specified in \ngn{namcfg}. 1159 The variable \np{lzoom} is used to specify that the damping is being used as in case \textit{a} above to 1160 provide boundary conditions to a zoom configuration. 1161 In the case of the arctic or antarctic zoom configurations this includes some specific treatment. 1162 Otherwise damping is applied to the 6 grid points along the ocean boundaries. 1163 The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in 1164 the \ngn{nam\_zoom\_dmp} name list. 1165 1166 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in 1167 non-zoom configurations. 1168 \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. 1169 \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for 1170 the ORCA4, ORCA2 and ORCA05 configurations. 1171 If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as 1172 a function of the model number. 1173 This option is included to allow backwards compatability of the ORCA2 reference configurations with 1174 previous model versions. 1175 \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. 1176 This option only has an effect if \np{ln\_full\_field} is true. 1177 \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. 1178 Finally \np{ln\_custom} specifies that the custom module will be called. 1179 This module is contained in the file \mdl{custom} and can be edited by users. 1180 For example damping could be applied in a specific region. 1181 1182 The restoration coefficient can be set to zero in equatorial regions by 1183 specifying a positive value of \np{nn\_hdmp}. 1184 Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to 1185 the full values of a 10\deg latitud band. 1186 This is often used because of the short adjustment time scale in the equatorial region 1187 \citep{Reverdin1991, Fujio1991, Marti_PhD92}. 1188 The time scale associated with the damping depends on the depth as a hyperbolic tangent, 1189 with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}. 1125 For generating \ifile{resto}, see the documentation for the DMP tool provided with the source code under 1126 \path{./tools/DMP_TOOLS}. 1190 1127 1191 1128 % ================================================================ … … 1199 1136 %-------------------------------------------------------------------------------------------------------------- 1200 1137 1201 Options are defined through the 1138 Options are defined through the \ngn{namdom} namelist variables. 1202 1139 The general framework for tracer time stepping is a modified leap-frog scheme \citep{Leclair_Madec_OM09}, 1203 1140 \ie a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:STP_mLF}): 1204 1141 \begin{equation} 1205 1142 \label{eq:tra_nxt} 1206 \begin{aligned} 1207 (e_{3t}T)^{t+\rdt} &= (e_{3t}T)_f^{t-\rdt} &+ 2 \, \rdt \,e_{3t}^t\ \text{RHS}^t & \\ \\ 1208 (e_{3t}T)_f^t \;\ \quad &= (e_{3t}T)^t \;\quad 1209 &+\gamma \,\left[ {(e_{3t}T)_f^{t-\rdt} -2(e_{3t}T)^t+(e_{3t}T)^{t+\rdt}} \right] & \\ 1210 & &- \gamma\,\rdt \, \left[ Q^{t+\rdt/2} - Q^{t-\rdt/2} \right] & 1211 \end{aligned} 1143 \begin{alignedat}{3} 1144 &(e_{3t}T)^{t + \rdt} &&= (e_{3t}T)_f^{t - \rdt} &&+ 2 \, \rdt \,e_{3t}^t \ \text{RHS}^t \\ 1145 &(e_{3t}T)_f^t &&= (e_{3t}T)^t &&+ \, \gamma \, \lt[ (e_{3t}T)_f^{t - \rdt} - 2(e_{3t}T)^t + (e_{3t}T)^{t + \rdt} \rt] \\ 1146 & && &&- \, \gamma \, \rdt \, \lt[ Q^{t + \rdt/2} - Q^{t - \rdt/2} \rt] 1147 \end{alignedat} 1212 1148 \end{equation} 1213 1149 where RHS is the right hand side of the temperature equation, the subscript $f$ denotes filtered values, … … 1215 1151 (\ie fluxes plus content in mass exchanges). 1216 1152 $\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 1217 Its default value is \np{rn\_atfp} \forcode{= 10.e-3}.1153 Its default value is \np{rn\_atfp}~\forcode{= 10.e-3}. 1218 1154 Note that the forcing correction term in the filter is not applied in linear free surface 1219 (\jp{lk\_vvl} \forcode{ = .false.}) (see \autoref{subsec:TRA_sbc}.1155 (\jp{lk\_vvl}~\forcode{= .false.}) (see \autoref{subsec:TRA_sbc}). 1220 1156 Not also that in constant volume case, the time stepping is performed on $T$, not on its content, $e_{3t}T$. 1221 1157 1222 When the vertical mixing is solved implicitly, 1223 the update of the \textit{next} tracer fields is done in module \mdl{trazdf}.1158 When the vertical mixing is solved implicitly, the update of the \textit{next} tracer fields is done in 1159 \mdl{trazdf} module. 1224 1160 In this case only the swapping of arrays and the Asselin filtering is done in the \mdl{tranxt} module. 1225 1161 1226 1162 In order to prepare for the computation of the \textit{next} time step, a swap of tracer arrays is performed: 1227 $T^{t -\rdt} = T^t$ and $T^t = T_f$.1163 $T^{t - \rdt} = T^t$ and $T^t = T_f$. 1228 1164 1229 1165 % ================================================================ … … 1240 1176 % Equation of State 1241 1177 % ------------------------------------------------------------------------------------------------------------- 1242 \subsection{Equation of seawater (\protect\np{nn\_eos} \forcode{= -1..1})}1178 \subsection{Equation of seawater (\protect\np{nn\_eos}~\forcode{= -1..1})} 1243 1179 \label{subsec:TRA_eos} 1244 1180 … … 1264 1200 To that purposed, a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available. 1265 1201 1266 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density.1202 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density. 1267 1203 Called \textit{rau0} in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 1268 1204 This is a sensible choice for the reference density used in a Boussinesq ocean climate model, as, … … 1270 1206 density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}. 1271 1207 1272 Options are defined through the 1208 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which 1273 1209 controls the EOS used (\forcode{= -1} for TEOS10 ; \forcode{= 0} for EOS-80 ; \forcode{= 1} for S-EOS). 1210 1274 1211 \begin{description} 1275 \item[\np{nn\_eos} \forcode{= -1}]1212 \item[\np{nn\_eos}~\forcode{= -1}] 1276 1213 the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 1277 1214 The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, … … 1282 1219 the TEOS-10 rational function approximation for hydrographic data analysis \citep{TEOS10}. 1283 1220 A key point is that conservative state variables are used: 1284 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \deg {C}, notation: $\Theta$).1221 Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \degC, notation: $\Theta$). 1285 1222 The pressure in decibars is approximated by the depth in meters. 1286 1223 With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. 1287 It is set to $C_p=3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}. 1288 1224 It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}. 1289 1225 Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$. 1290 1226 In particular, the initial state deined by the user have to be given as \textit{Conservative} Temperature and … … 1293 1229 either computing the air-sea and ice-sea fluxes (forced mode) or 1294 1230 sending the SST field to the atmosphere (coupled mode). 1295 1296 \item[\np{nn\_eos}\forcode{ = 0}] 1231 \item[\np{nn\_eos}~\forcode{= 0}] 1297 1232 the polyEOS80-bsq equation of seawater is used. 1298 1233 It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized to … … 1305 1240 pressure \citep{UNESCO1983}. 1306 1241 Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which 1307 is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value. 1308 1309 \item[\np{nn\_eos}\forcode{ = 1}] 1242 is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value. 1243 \item[\np{nn\_eos}~\forcode{= 1}] 1310 1244 a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 1311 1245 the coefficients of which has been optimized to fit the behavior of TEOS10 … … 1317 1251 as well as between \textit{absolute} and \textit{practical} salinity. 1318 1252 S-EOS takes the following expression: 1319 \ [1253 \begin{gather*} 1320 1254 % \label{eq:tra_S-EOS} 1321 \begin{split} 1322 d_a(T,S,z) = ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a \\ 1323 & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a \\ 1324 & - \nu \; T_a \; S_a \; ) \; / \; \rho_o \\ 1325 with \ \ T_a = T-10 \; ; & \; S_a = S-35 \; ;\; \rho_o = 1026~Kg/m^3 1326 \end{split} 1327 \] 1255 \begin{alignedat}{2} 1256 &d_a(T,S,z) = \frac{1}{\rho_o} \big[ &- a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * &T_a \big. \\ 1257 & &+ b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * &S_a \\ 1258 & \big. &- \nu \; T_a &S_a \big] \\ 1259 \end{alignedat} 1260 \\ 1261 \text{with~} T_a = T - 10 \, ; \, S_a = S - 35 \, ; \, \rho_o = 1026~Kg/m^3 1262 \end{gather*} 1328 1263 where the computer name of the coefficients as well as their standard value are given in \autoref{tab:SEOS}. 1329 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing the associated coefficients. 1330 Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS. 1331 setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS. 1264 In fact, when choosing S-EOS, various approximation of EOS can be specified simply by 1265 changing the associated coefficients. 1266 Setting to zero the two thermobaric coefficients $(\mu_1,\mu_2)$ remove thermobaric effect from S-EOS. 1267 setting to zero the three cabbeling coefficients $(\lambda_1,\lambda_2,\nu)$ remove cabbeling effect from 1268 S-EOS. 1332 1269 Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 1333 1270 \end{description} 1334 1335 1271 1336 1272 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1337 1273 \begin{table}[!tb] 1338 1274 \begin{center} 1339 \begin{tabular}{| p{26pt}|p{72pt}|p{56pt}|p{136pt}|}1275 \begin{tabular}{|l|l|l|l|} 1340 1276 \hline 1341 coeff. & computer name & S-EOS & description \\ \hline 1342 $a_0$ & \np{rn\_a0} & 1.6550 $10^{-1}$ & linear thermal expansion coeff. \\ \hline 1343 $b_0$ & \np{rn\_b0} & 7.6554 $10^{-1}$ & linear haline expansion coeff. \\ \hline 1344 $\lambda_1$ & \np{rn\_lambda1}& 5.9520 $10^{-2}$ & cabbeling coeff. in $T^2$ \\ \hline 1345 $\lambda_2$ & \np{rn\_lambda2}& 5.4914 $10^{-4}$ & cabbeling coeff. in $S^2$ \\ \hline 1346 $\nu$ & \np{rn\_nu} & 2.4341 $10^{-3}$ & cabbeling coeff. in $T \, S$ \\ \hline 1347 $\mu_1$ & \np{rn\_mu1} & 1.4970 $10^{-4}$ & thermobaric coeff. in T \\ \hline 1348 $\mu_2$ & \np{rn\_mu2} & 1.1090 $10^{-5}$ & thermobaric coeff. in S \\ \hline 1277 coeff. & computer name & S-EOS & description \\ 1278 \hline 1279 $a_0$ & \np{rn\_a0} & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\ 1280 \hline 1281 $b_0$ & \np{rn\_b0} & $7.6554~10^{-1}$ & linear haline expansion coeff. \\ 1282 \hline 1283 $\lambda_1$ & \np{rn\_lambda1}& $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$ \\ 1284 \hline 1285 $\lambda_2$ & \np{rn\_lambda2}& $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$ \\ 1286 \hline 1287 $\nu$ & \np{rn\_nu} & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$ \\ 1288 \hline 1289 $\mu_1$ & \np{rn\_mu1} & $1.4970~10^{-4}$ & thermobaric coeff. in T \\ 1290 \hline 1291 $\mu_2$ & \np{rn\_mu2} & $1.1090~10^{-5}$ & thermobaric coeff. in S \\ 1292 \hline 1349 1293 \end{tabular} 1350 1294 \caption{ … … 1352 1296 Standard value of S-EOS coefficients. 1353 1297 } 1354 1298 \end{center} 1355 1299 \end{table} 1356 1300 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1357 1301 1358 1359 1302 % ------------------------------------------------------------------------------------------------------------- 1360 1303 % Brunt-V\"{a}is\"{a}l\"{a} Frequency 1361 1304 % ------------------------------------------------------------------------------------------------------------- 1362 \subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos} \forcode{= 0..2})}1305 \subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}~\forcode{= 0..2})} 1363 1306 \label{subsec:TRA_bn2} 1364 1307 1365 An accurate computation of the ocean stability ( \ieof $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of1308 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of 1366 1309 paramount importance as determine the ocean stratification and is used in several ocean parameterisations 1367 1310 (namely TKE, GLS, Richardson number dependent vertical diffusion, enhanced vertical diffusion, … … 1372 1315 \[ 1373 1316 % \label{eq:tra_bn2} 1374 N^2 = \frac{g}{e_{3w}} \left( \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T] \right)1317 N^2 = \frac{g}{e_{3w}} \lt( \beta \; \delta_{k + 1/2}[S] - \alpha \; \delta_{k + 1/2}[T] \rt) 1375 1318 \] 1376 where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS,1377 and,$\alpha$ and $\beta$ are the thermal and haline expansion coefficients.1378 The coefficients are a polynomial function of temperature, salinity and depth which 1379 expression depends onthe chosen EOS.1319 where $(T,S) = (\Theta,S_A)$ for TEOS10, $(\theta,S_p)$ for TEOS-80, or $(T,S)$ for S-EOS, and, 1320 $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 1321 The coefficients are a polynomial function of temperature, salinity and depth which expression depends on 1322 the chosen EOS. 1380 1323 They are computed through \textit{eos\_rab}, a \fortran function that can be found in \mdl{eosbn2}. 1381 1324 … … 1390 1333 \label{eq:tra_eos_fzp} 1391 1334 \begin{split} 1392 T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} - 2.154996 \;10^{-4} \,S \right) \ S \\ 1393 - 7.53\,10^{-3} \ \ p 1394 \end{split} 1335 &T_f (S,p) = \lt( a + b \, \sqrt{S} + c \, S \rt) \, S + d \, p \\ 1336 &\text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \\ 1337 &\text{and~} d = -7.53~10^{-3} 1338 \end{split} 1395 1339 \end{equation} 1396 1340 1397 1341 \autoref{eq:tra_eos_fzp} is only used to compute the potential freezing point of sea water 1398 (\ie referenced to the surface $p =0$),1342 (\ie referenced to the surface $p = 0$), 1399 1343 thus the pressure dependent terms in \autoref{eq:tra_eos_fzp} (last term) have been dropped. 1400 1344 The freezing point is computed through \textit{eos\_fzp}, 1401 a \fortran function that can be found in \mdl{eosbn2}. 1402 1345 a \fortran function that can be found in \mdl{eosbn2}. 1403 1346 1404 1347 % ------------------------------------------------------------------------------------------------------------- … … 1411 1354 % 1412 1355 1413 1414 1356 % ================================================================ 1415 1357 % Horizontal Derivative in zps-coordinate … … 1421 1363 I've changed "derivative" to "difference" and "mean" to "average"} 1422 1364 1423 With partial cells (\np{ln\_zps} \forcode{ = .true.}) at bottom and top (\np{ln\_isfcav}\forcode{= .true.}),1365 With partial cells (\np{ln\_zps}~\forcode{= .true.}) at bottom and top (\np{ln\_isfcav}~\forcode{= .true.}), 1424 1366 in general, tracers in horizontally adjacent cells live at different depths. 1425 1367 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and 1426 1368 the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 1427 The partial cell properties at the top (\np{ln\_isfcav} \forcode{= .true.}) are computed in the same way as1369 The partial cell properties at the top (\np{ln\_isfcav}~\forcode{= .true.}) are computed in the same way as 1428 1370 for the bottom. 1429 1371 So, only the bottom interpolation is explained below. … … 1432 1374 a linear interpolation in the vertical is used to approximate the deeper tracer as if 1433 1375 it actually lived at the depth of the shallower tracer point (\autoref{fig:Partial_step_scheme}). 1434 For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde {T}$, is:1376 For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde T$, is: 1435 1377 1436 1378 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1437 1379 \begin{figure}[!p] 1438 1380 \begin{center} 1439 \includegraphics[ width=0.9\textwidth]{Fig_partial_step_scheme}1381 \includegraphics[]{Fig_partial_step_scheme} 1440 1382 \caption{ 1441 1383 \protect\label{fig:Partial_step_scheme} 1442 1384 Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate 1443 (\protect\np{ln\_zps} \forcode{ = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i )>0$.1444 A linear interpolation is used to estimate $\widetilde {T}_k^{i+1}$,1385 (\protect\np{ln\_zps}~\forcode{= .true.}) in the case $(e3w_k^{i + 1} - e3w_k^i) > 0$. 1386 A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$, 1445 1387 the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 1446 The horizontal difference is then given by: $\delta_{i +1/2} T_k= \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and1447 the average by: $\overline {T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i}) / 2$.1388 The horizontal difference is then given by: $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and 1389 the average by: $\overline T_k^{\, i + 1/2} = (\widetilde T_k^{\, i + 1/2} - T_k^{\, i}) / 2$. 1448 1390 } 1449 1391 \end{center} … … 1451 1393 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1452 1394 \[ 1453 \widetilde {T}= \left\{1454 \begin{aligned }1455 &T^{\, i+1} -\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right)}{ e_{3w}^{i+1} }\;\delta_k T^{i+1}1456 & & \quad\text{if $\ e_{3w}^{i+1} \geq e_{3w}^i$ }\\ \\1457 &T^{\, i} \ \ \ \,+\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right) }{e_{3w}^i }\;\delta_k T^{i+1}1458 & & \quad\text{if $\ e_{3w}^{i+1} < e_{3w}^i$}1459 \end{aligned }1460 \r ight.1395 \widetilde T = \lt\{ 1396 \begin{alignedat}{2} 1397 &T^{\, i + 1} &-\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt) }{ e_{3w}^{i + 1} } \; \delta_k T^{i + 1} 1398 & \quad \text{if $e_{3w}^{i + 1} \geq e_{3w}^i$} \\ \\ 1399 &T^{\, i} &+\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt )}{e_{3w}^i } \; \delta_k T^{i + 1} 1400 & \quad \text{if $e_{3w}^{i + 1} < e_{3w}^i$} 1401 \end{alignedat} 1402 \rt. 1461 1403 \] 1462 1404 and the resulting forms for the horizontal difference and the horizontal average value of $T$ at a $U$-point are: 1463 1405 \begin{equation} 1464 1406 \label{eq:zps_hde} 1465 \begin{ aligned}1466 \delta_{i +1/2} T=1407 \begin{split} 1408 \delta_{i + 1/2} T &= 1467 1409 \begin{cases} 1468 \ \ \ \widetilde {T}\quad\ -T^i & \ \ \quad\quad\text{if $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ \\ 1469 \ \ \ T^{\,i+1}-\widetilde{T} & \ \ \quad\quad\text{if $\ e_{3w}^{i+1} < e_{3w}^i$ } 1410 \widetilde T - T^i & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 1411 \\ 1412 T^{\, i + 1} - \widetilde T & \text{if~} e_{3w}^{i + 1} < e_{3w}^i 1470 1413 \end{cases} 1471 \\ \\1472 \overline {T}^{\,i+1/2} \=1414 \\ 1415 \overline T^{\, i + 1/2} &= 1473 1416 \begin{cases} 1474 ( \widetilde {T}\ \ \;\,-T^{\,i}) / 2 & \;\ \ \quad\text{if $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ \\ 1475 ( T^{\,i+1}-\widetilde{T} ) / 2 & \;\ \ \quad\text{if $\ e_{3w}^{i+1} < e_{3w}^i$ } 1417 (\widetilde T - T^{\, i} ) / 2 & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 1418 \\ 1419 (T^{\, i + 1} - \widetilde T) / 2 & \text{if~} e_{3w}^{i + 1} < e_{3w}^i 1476 1420 \end{cases} 1477 \end{ aligned}1421 \end{split} 1478 1422 \end{equation} 1479 1423 1480 1424 The computation of horizontal derivative of tracers as well as of density is performed once for all at 1481 1425 each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed. 1482 It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde {\rho}$,1426 It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde \rho$, 1483 1427 is not the same as that used for $T$ and $S$. 1484 Instead of forming a linear approximation of density, we compute $\widetilde {\rho }$ from the interpolated values of1428 Instead of forming a linear approximation of density, we compute $\widetilde \rho$ from the interpolated values of 1485 1429 $T$ and $S$, and the pressure at a $u$-point 1486 (in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos} 1430 (in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos}): 1487 1431 \[ 1488 1432 % \label{eq:zps_hde_rho} 1489 \widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u }) 1490 \quad \text{where }\ z_u = \min \left( {z_T^{i+1} ,z_T^i } \right) 1433 \widetilde \rho = \rho (\widetilde T,\widetilde S,z_u) \quad \text{where~} z_u = \min \lt( z_T^{i + 1},z_T^i \rt) 1491 1434 \] 1492 1435
Note: See TracChangeset
for help on using the changeset viewer.