- Timestamp:
- 2010-10-15T16:42:00+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_ZDF.tex
r1225 r2282 7 7 8 8 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. 9 \gmcomment{Steven remark (not taken into account : problem here with turbulent vs turbulence. I've changed "turbulent closure" to "turbulence closure", "turbulent mixing length" to "turbulence mixing length", but I've left "turbulent kinetic energy" alone - though I think it is an historical abberation! 10 Gurvan : I kept "turbulent closure etc "...} 9 10 11 \newpage 12 $\ $\newline % force a new ligne 11 13 12 14 … … 38 40 general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 39 41 These trends can be computed using either a forward time stepping scheme 40 (namelist parameter \np{ np\_zdfexp}=true) or a backward time stepping41 scheme (\np{ np\_zdfexp}=false) depending on the magnitude of the mixing42 coefficients, and thus of the formulation used (see \S\ref{ DOM_nxt}).42 (namelist parameter \np{ln\_zdfexp}=true) or a backward time stepping 43 scheme (\np{ln\_zdfexp}=false) depending on the magnitude of the mixing 44 coefficients, and thus of the formulation used (see \S\ref{STP}). 43 45 44 46 % ------------------------------------------------------------------------------------------------------------- … … 57 59 \begin{align*} 58 60 A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1} \\ 59 \\60 61 A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1} 61 62 \end{align*} 62 63 63 These values are set through the \np{ avm0} and \np{avt0} namelist parameters.64 These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. 64 65 In all cases, do not use values smaller that those associated with the molecular 65 66 viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, … … 74 75 75 76 %--------------------------------------------namric--------------------------------------------------------- 76 \namdisplay{nam ric}77 \namdisplay{namzdf_ric} 77 78 %-------------------------------------------------------------------------------------------------------------- 78 79 … … 84 85 growth of Kelvin-Helmholtz like instabilities leads to a dependency between the 85 86 vertical eddy coefficients and the local Richardson number ($i.e.$ the 86 ratio of stratification to vertical shear). Following \citet{Pac Phil1981}, the following87 ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following 87 88 formulation has been implemented: 88 89 \begin{equation} \label{Eq_zdfric} 89 90 \left\{ \begin{aligned} 90 91 A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT} \\ 91 \\92 92 A^{vm} &= \frac{A^{vT} }{\left( 1+ a \;Ri \right) } + A_b^{vm} 93 93 \end{aligned} \right. … … 99 99 is the maximum value that can be reached by the coefficient when $Ri\leq 0$, 100 100 $a=5$ and $n=2$. The last three values can be modified by setting the 101 \np{ avmri}, \np{alp} and \np{nric} namelist parameters, respectively.101 \np{rn\_avmri}, \np{rn\_alp} and \np{nn\_ric} namelist parameters, respectively. 102 102 103 103 % ------------------------------------------------------------------------------------------------------------- … … 107 107 \label{ZDF_tke} 108 108 109 %--------------------------------------------nam tke---------------------------------------------------------110 \namdisplay{nam tke}109 %--------------------------------------------namzdf_tke-------------------------------------------------- 110 \namdisplay{namzdf_tke} 111 111 %-------------------------------------------------------------------------------------------------------------- 112 112 113 113 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE 114 turbulent closure model based on a prognostic equation for $\bar 114 turbulent closure model based on a prognostic equation for $\bar{e}$, the turbulent 115 115 kinetic energy, and a closure assumption for the turbulent length scales. This 116 116 turbulent closure model has been developed by \citet{Bougeault1989} in the 117 117 atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and 118 embedded in OPA by \citet{Blanke1993} for equatorial Atlantic simulations. Since119 then, significant modifications have been introduced by \citet{Madec1998} in both120 the implementation and the formulation of the mixing length scale. The time121 evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical122 shear, its destruction through stratification, its vertical diffusion, and its dissipation123 of \citet{Kolmogorov1942} type:118 embedded in OPA, the ancestor of NEMO, by \citet{Blanke1993} for equatorial Atlantic 119 simulations. Since then, significant modifications have been introduced by 120 \citet{Madec1998} in both the implementation and the formulation of the mixing 121 length scale. The time evolution of $\bar{e}$ is the result of the production of 122 $\bar{e}$ through vertical shear, its destruction through stratification, its vertical 123 diffusion, and its dissipation of \citet{Kolmogorov1942} type: 124 124 \begin{equation} \label{Eq_zdftke_e} 125 125 \frac{\partial \bar{e}}{\partial t} = 126 \frac{ A^{vm}}{e_3}\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2127 128 - A^{vT}\,N^2126 \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2 127 +\left( {\frac{\partial v}{\partial k}} \right)^2} \right] 128 -K_\rho\,N^2 129 129 +\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 } 130 130 \;\frac{\partial \bar{e}}{\partial k}} \right] … … 133 133 \begin{equation} \label{Eq_zdftke_kz} 134 134 \begin{split} 135 A^{vm} &= C_k\ l_k\ \sqrt {\bar{e}} \\136 A^{vT}&= A^{vm} / P_{rt}135 K_m &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ 136 K_\rho &= A^{vm} / P_{rt} 137 137 \end{split} 138 138 \end{equation} 139 139 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 140 140 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 141 $P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and 142 $C_\epsilon = 0.1$ are designed to deal with vertical mixing at any depth 143 \citep{Gaspar1990}. They are set through namelist parameters \np{ediff} 144 and \np{ediss}. $P_{rt} $ can be set to unity or, following \citet{Blanke1993}, 145 be a function of the local Richardson number, $R_i $: 141 $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity 142 and diffusivity coefficients. The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ 143 $\approx 0.7$ are designed to deal with vertical mixing at any depth \citep{Gaspar1990}. 144 They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}. 145 $P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function 146 of the local Richardson number, $R_i$: 146 147 \begin{align*} \label{Eq_prt} 147 148 P_{rt} = \begin{cases} … … 151 152 \end{cases} 152 153 \end{align*} 153 Note that a horizontal Shapiro filter can optionally be applied to $R_i$. 154 However it is an obsolescent option that is not recommended. 155 The choice of $P_{rt} $ is controlled by the \np{npdl} namelist parameter. 154 The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist parameter. 156 155 157 156 For computational efficiency, the original formulation of the turbulent length 158 157 scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 159 are proposed, the choice of which is controlled by the \np{n mxl} namelist158 are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist 160 159 parameter. The first two are based on the following first order approximation 161 160 \citep{Blanke1993}: 162 161 \begin{equation} \label{Eq_tke_mxl0_1} 163 l_k = l_\epsilon = \sqrt {2 \bar e} / N164 \end{equation} 165 which is valid in a stable stratified region with constant values of the brunt-162 l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N 163 \end{equation} 164 which is valid in a stable stratified region with constant values of the Brunt- 166 165 Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance 167 to the surface or to the bottom (\np{nmxl}=0) or by the local vertical scale factor (\np{nmxl}=1). \citet{Blanke1993} notice that this simplification has two major 166 to the surface or to the bottom (\np{nn\_mxl} = 0) or by the local vertical scale factor 167 (\np{nn\_mxl} = 1). \citet{Blanke1993} notice that this simplification has two major 168 168 drawbacks: it makes no sense for locally unstable stratification and the 169 169 computation no longer uses all the information contained in the vertical density 170 170 profile. To overcome these drawbacks, \citet{Madec1998} introduces the 171 \np{n mxl}=2 or 3 cases, which add an extra assumption concerning the vertical171 \np{nn\_mxl} = 2 or 3 cases, which add an extra assumption concerning the vertical 172 172 gradient of the computed length scale. So, the length scales are first evaluated 173 173 as in \eqref{Eq_tke_mxl0_1} and then bounded such that: … … 185 185 constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, 186 186 the upward and downward length scales, and evaluate the dissipation and 187 mixing turbulent length scales as (and note that here we use numerical 188 indexing): 187 mixing length scales as (and note that here we use numerical indexing): 189 188 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 190 189 \begin{figure}[!t] \label{Fig_mixing_length} \begin{center} … … 196 195 \begin{equation} \label{Eq_tke_mxl2} 197 196 \begin{aligned} 198 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3 T}^{(k)}\ \ \ \; \right)197 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) 199 198 \quad &\text{ from $k=1$ to $jpk$ }\ \\ 200 l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3 T}^{(k-1)} \right)199 l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)} \right) 201 200 \quad &\text{ from $k=jpk$ to $1$ }\ \\ 202 201 \end{aligned} 203 202 \end{equation} 204 203 where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, 205 $i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$.206 207 In the \np{n mxl}=2 case, the dissipation and mixing length scales take the same204 $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. 205 206 In the \np{nn\_mxl}=2 case, the dissipation and mixing length scales take the same 208 207 value: $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the 209 \np{n mxl}=2 case, the dissipation and mixinglength scales are give208 \np{nn\_mxl}=2 case, the dissipation and mixing turbulent length scales are give 210 209 as in \citet{Gaspar1990}: 211 210 \begin{equation} \label{Eq_tke_mxl_gaspar} … … 217 216 218 217 At the sea surface the value of $\bar{e}$ is prescribed from the wind 219 stress field: $\bar{e}= ebb\;\left| \tau \right|$ ($ebb=60$by default)220 with a minimal threshold of $emin0=10^{-4}~m^2.s^{-2}$ (namelist218 stress field: $\bar{e}=rn\_ebb\;\left| \tau \right|$ (\np{rn\_ebb}=60 by default) 219 with a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist 221 220 parameters). Its value at the bottom of the ocean is assumed to be 222 221 equal to the value of the level just above. The time integration of the 223 222 $\bar{e}$ equation may formally lead to negative values because the 224 223 numerical scheme does not ensure its positivity. To overcome this 225 problem, a cut-off in the minimum value of $\bar{e}$ is used. Following 226 \citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 227 This allows the subsequent formulations to match that of\citet{Gargett1984} 228 for the diffusion in the thermocline and deep ocean : $(A^{vT} = 10^{-3} / N)$. 229 In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical 224 problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} 225 namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set 226 to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations 227 to match that of \citet{Gargett1984} for the diffusion in the thermocline and 228 deep ocean : $K_\rho = 10^{-3} / N$. 229 In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical 230 230 instabilities associated with too weak vertical diffusion. They must be 231 231 specified at least larger than the molecular values, and are set through 232 \textit{avm0} and \textit{avt0} (\textbf{namelist} parameters). 232 \np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}). 233 234 % ------------------------------------------------------------------------------------------------------------- 235 % TKE Turbulent Closure Scheme : new organization to energetic considerations 236 % ------------------------------------------------------------------------------------------------------------- 237 \subsection{TKE discretization considerations (\key{zdftke})} 238 \label{ZDF_tke_ene} 239 240 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 241 \begin{figure}[!t] \label{Fig_TKE_time_scheme} \begin{center} 242 \includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_ZDF_TKE_time_scheme.pdf} 243 \caption {Illustration of the TKE time integration and its links to the momentum and tracer time integration. } 244 \end{center} 245 \end{figure} 246 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 247 248 The production of turbulence by vertical shear (the first term of the right hand side 249 of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with 250 the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}). To do so a special care 251 have to be taken for both the time and space discretization of the TKE equation \citep{Burchard_OM02}. 252 253 Let us first address the time stepping issue. Fig.~\ref{Fig_TKE_time_scheme} shows 254 how the two-level Leap-Frog time stepping of the momentum and tracer equations interplays 255 with the one-level forward time stepping of TKE equation. With this framework, the total loss 256 of kinetic energy (in 1D for the demonstration) due to the vertical momentum diffusion is 257 obtained by multiplying this quantity by $u^t$ and summing the result vertically: 258 \begin{equation} \label{Eq_energ1} 259 \begin{split} 260 \int_{-H}^{\eta} u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt} \right) \,dz \\ 261 &= \Bigl[ u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta} 262 - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz } 263 \end{split} 264 \end{equation} 265 Here, the vertical diffusion of momentum is discretized backward in time 266 with a coefficient, $K_m$, known at time $t$ (Fig.~\ref{Fig_TKE_time_scheme}), 267 as it is required when using the TKE scheme (see \S\ref{STP_forward_imp}). 268 The first term of the right hand side of \eqref{Eq_energ1} represents the kinetic energy 269 transfer at the surface (atmospheric forcing) and at the bottom (friction effect). 270 The second term is always negative. It is the dissipation rate of kinetic energy, 271 and thus minus the shear production rate of $\bar{e}$. \eqref{Eq_energ1} 272 implies that, to be energetically consistent, the production rate of $\bar{e}$ 273 used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 274 ${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ (and not by the more straightforward 275 $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$). 276 277 A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 278 (second term of the right hand side of \eqref{Eq_zdftke_e}). This term 279 must balance the input of potential energy resulting from vertical mixing. 280 The rate of change of potential energy (in 1D for the demonstration) due vertical 281 mixing is obtained by multiplying vertical density diffusion 282 tendency by $g\,z$ and and summing the result vertically: 283 \begin{equation} \label{Eq_energ2} 284 \begin{split} 285 \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt} \right) \,dz \\ 286 &= \Bigl[ g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta} 287 - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz \\ 288 &= - \Bigl[ z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta} 289 + \int_{-H}^{\eta}{ \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz } 290 \end{split} 291 \end{equation} 292 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 293 The first term of the right hand side of \eqref{Eq_energ2} is always zero 294 because there is no diffusive flux through the ocean surface and bottom). 295 The second term is minus the destruction rate of $\bar{e}$ due to stratification. 296 Therefore \eqref{Eq_energ1} implies that, to be energetically consistent, the product 297 ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \eqref{Eq_zdftke_e}, the TKE equation. 298 299 Let us now address the space discretization issue. 300 The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity 301 components are in the centre of the side faces of a $t$-box in staggered C-grid 302 (Fig.\ref{Fig_cell}). A space averaging is thus required to obtain the shear TKE production term. 303 By redoing the \eqref{Eq_energ1} in the 3D case, it can be shown that the product of 304 eddy coefficient by the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. 305 Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into 306 account. 307 308 The above energetic considerations leads to 309 the following final discrete form for the TKE equation: 310 \begin{equation} \label{Eq_zdftke_ene} 311 \begin{split} 312 \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt} \equiv 313 \Biggl\{ \Biggr. 314 &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} } 315 \ \frac{\delta_{k+1/2}[u^ t ]}{{e_3u}^ t } \right) }^{\,i} \\ 316 +&\overline{ \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} } 317 \ \frac{\delta_{k+1/2}[v^ t ]}{{e_3v}^ t } \right) }^{\,j} 318 \Biggr. \Biggr\} \\ 319 % 320 - &{K_\rho}^{t-\rdt}\,{(N^2)^t} \\ 321 % 322 +&\frac{1}{{e_3w}^{t+\rdt}} \;\delta_{k+1/2} \left[ {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}} \right] \\ 323 % 324 - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt} 325 \end{split} 326 \end{equation} 327 where the last two terms in \eqref{Eq_zdftke_ene} (vertical diffusion and Kolmogorov dissipation) 328 are time stepped using a backward scheme (see\S\ref{STP_forward_imp}). 329 Note that the Kolmogorov term has been linearized in time in order to render 330 the implicit computation possible. The restart of the TKE scheme 331 requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as they all appear in 332 the right hand side of \eqref{Eq_zdftke_ene}. For the latter, it is in fact 333 the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 334 335 % ------------------------------------------------------------------------------------------------------------- 336 % GLS Generic Length Scale Scheme 337 % ------------------------------------------------------------------------------------------------------------- 338 \subsection{GLS Generic Length Scale (\key{zdfgls})} 339 \label{ZDF_gls} 340 341 %--------------------------------------------namgls--------------------------------------------------------- 342 \namdisplay{namgls} 343 %-------------------------------------------------------------------------------------------------------------- 344 345 The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on 346 two prognostic equations: one for the turbulent kinetic energy $\bar {e}$, and another 347 for the generic length scale, $\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}. 348 This later variable is defined as : $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 349 where the triplet $(p, m, n)$ value given in Tab.\ref{Tab_GLS} allows to recover 350 a number of well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982}, 351 $k$-$\epsilon$ \citep{Rodi_1987}, $k$-$\omega$ \citep{Wilcox_1988} 352 among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}). 353 The GLS scheme is given by the following set of equations: 354 \begin{equation} \label{Eq_zdfgls_e} 355 \frac{\partial \bar{e}}{\partial t} = 356 \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 357 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 358 -K_\rho \,N^2 359 +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right] 360 - \epsilon 361 \end{equation} 362 363 \begin{equation} \label{Eq_zdfgls_psi} 364 \begin{split} 365 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ 366 \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 367 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 368 - C_3 \,K_\rho\,N^2 - C_2 \,\epsilon \,Fw \right\} \\ 369 &+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 } 370 \;\frac{\partial \psi}{\partial k}} \right]\; 371 \end{split} 372 \end{equation} 373 374 \begin{equation} \label{Eq_zdfgls_kz} 375 \begin{split} 376 K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ 377 K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l 378 \end{split} 379 \end{equation} 380 381 \begin{equation} \label{Eq_zdfgls_eps} 382 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; 383 \end{equation} 384 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}) 385 and $\epsilon$ the dissipation rate. 386 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) 387 depends of the choice of the turbulence model. Four different turbulent models are pre-defined 388 (Tab.\ref{Tab_GLS}). They are made available through th \np{gls} namelist parameter. 389 390 %--------------------------------------------------TABLE-------------------------------------------------- 391 \begin{table}[htbp] \label{Tab_GLS} 392 \begin{center} 393 %\begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} 394 \begin{tabular}{ccccc} 395 & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ 396 % & \citep{Mellor_Yamada_1982} & \citep{Rodi_1987} & \citep{Wilcox_1988} & \\ 397 \hline \hline 398 \np{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ 399 \hline 400 $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ 401 $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ 402 $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ 403 $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ 404 $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ 405 $C_3$ & 1. & 1. & 1. & 1. \\ 406 $F_{wall}$ & Yes & -- & -- & -- \\ 407 \hline 408 \hline 409 \end{tabular} 410 \caption {Set of predefined GLS parameters, or equivalently predefined turbulence models available with \key{gls} and controlled by the \np{nn\_clos} namelist parameter.} 411 \end{center} 412 \end{table} 413 %-------------------------------------------------------------------------------------------------------------- 414 415 In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force 416 the convergence of the mixing length towards $K\,z_b$ ($K$: Kappa and $z_b$: rugosity length) 417 value near physical boundaries (logarithmic boundary layer law). $C_{\mu}$ and $C_{\mu'}$ 418 are calculated from stability function proposed by \citet{Galperin_al_JAS88}, or by \citet{Kantha_Clayson_1994} 419 or one of the two functions suggested by \citet{Canuto_2001} (\np{nn\_stab\_func} = 0, 1, 2 or 3, resp.}). The value of $C_{0\mu}$ depends of the choice of the stability function. 420 421 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated 422 thanks to Dirichlet or Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp. 423 The wave effect on the mixing could be also being considered \citep{Craig_Banner_1994}. 424 425 The $\psi$ equation is known to fail in stably stratified flows, and for this reason 426 almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy. 427 With this clipping, the maximum permissible length scale is determined by 428 $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. A value of $c_{lim} = 0.53$ is often used 429 \citep{Galperin_al_JAS88}. \cite{Umlauf_Burchard_CSR05} show that the value of 430 the clipping factor is of crucial importance for the entrainment depth predicted in 431 stably stratified situations, and that its value has to be chosen in accordance 432 with the algebraic model for the turbulent ßuxes. The clipping is only activated 433 if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{clim\_galp} value. 233 434 234 435 % ------------------------------------------------------------------------------------------------------------- … … 239 440 240 441 %--------------------------------------------namkpp-------------------------------------------------------- 241 \namdisplay{namkpp} 242 %-------------------------------------------------------------------------------------------------------------- 243 244 The K-Profile Parametrization (KKP) developed by \cite{Large_al_RG94} has been 245 implemented in \NEMO by J. Chanut (PhD reference to be added here!). 442 \namdisplay{namzdf_kpp} 443 %-------------------------------------------------------------------------------------------------------------- 444 445 The KKP scheme has been implemented by J. Chanut ... 246 446 247 447 \colorbox{yellow}{Add a description of KPP here.} … … 274 474 \label{ZDF_npc} 275 475 276 %--------------------------------------------namnpc-------------------------------------------------------- 277 \namdisplay{namnpc} 278 %-------------------------------------------------------------------------------------------------------------- 279 476 %--------------------------------------------namzdf-------------------------------------------------------- 477 \namdisplay{namzdf} 478 %-------------------------------------------------------------------------------------------------------------- 280 479 281 480 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 295 494 296 495 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true. 297 It is applied at each \np{nn pc1} time step and mixes downwards instantaneously496 It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously 298 497 the statically unstable portion of the water column, but only until the density 299 498 structure becomes neutrally stable ($i.e.$ until the mixed portion of the water 300 column has \textit{exactly} the density of the water just below) \citep{Madec 1991a}.499 column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}. 301 500 The associated algorithm is an iterative process used in the following way 302 501 (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is … … 320 519 convective algorithm has been proved successful in studies of the deep 321 520 water formation in the north-western Mediterranean Sea 322 \citep{Madec 1991a, Madec1991b, Madec1991c}.521 \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}. 323 522 324 523 Note that in the current implementation of this algorithm presents several … … 341 540 after a static adjustment. In that case, we recommend the addition of momentum 342 541 mixing in a manner that mimics the mixing in temperature and salinity 343 \citep{Speich 1992, Speich1996}.542 \citep{Speich_PhD92, Speich_al_JPO96}. 344 543 345 544 % ------------------------------------------------------------------------------------------------------------- … … 347 546 % ------------------------------------------------------------------------------------------------------------- 348 547 \subsection [Enhanced Vertical Diffusion (\np{ln\_zdfevd})] 349 {Enhanced Vertical Diffusion (\np{ln\_zdfevd}= .true.)}548 {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=true)} 350 549 \label{ZDF_evd} 351 550 … … 357 556 In this case, the vertical eddy mixing coefficients are assigned very large values 358 557 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 359 ($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar1997, Lazar1999}.360 This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and361 tracers (\np{n\_evdm}=1).558 ($i.e.$ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) 559 \citep{Lazar_PhD97, Lazar_al_JPO99}. This is done either on tracers only 560 (\np{nn\_evdm}=0) or on both momentum and tracers (\np{nn\_evdm}=1). 362 561 363 562 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and 364 if \np{n \_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$365 values also, are set equal to the namelist parameter \np{ avevd}. A typical value366 for $ avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of563 if \np{nn\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 564 values also, are set equal to the namelist parameter \np{rn\_avevd}. A typical value 565 for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of 367 566 convective processes is less time consuming than the convective adjustment 368 567 algorithm presented above when mixing both tracers and momentum in the 369 568 case of static instabilities. It requires the use of an implicit time stepping on 370 569 vertical diffusion terms (i.e. \np{ln\_zdfexp}=false). 570 571 Note that the stability test is performed on both \textit{before} and \textit{now} 572 values of $N^2$. This removes a potential source of divergence of odd and 573 even time step in a leapfrog environment \citep{Leclair_PhD2010} (see \S\ref{STP_mLF}). 371 574 372 575 % ------------------------------------------------------------------------------------------------------------- … … 394 597 The KPP turbulent closure scheme already includes enhanced vertical diffusion 395 598 in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ 396 found in \mdl{zdfkpp}, therefore \np{ np\_zdfevd} should notbe used with the KPP599 found in \mdl{zdfkpp}, therefore \np{ln\_zdfevd}=false should be used with the KPP 397 600 scheme. %gm% + one word on non local flux with KPP scheme trakpp.F90 module... 398 601 … … 404 607 \label{ZDF_ddm} 405 608 406 %-------------------------------------------nam ddm--------------------------------------------------------407 \namdisplay{nam ddm}609 %-------------------------------------------namzdf_ddm------------------------------------------------- 610 \namdisplay{namzdf_ddm} 408 611 %-------------------------------------------------------------------------------------------------------------- 409 612 … … 422 625 \end{align*} 423 626 where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection, 424 and $o$ by processes other than double diffusion. The rates of double-diffusive mixing 425 depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 627 and $o$ by processes other than double diffusion. The rates of double-diffusive mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 426 628 where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline 427 629 contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt … … 454 656 we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 455 657 456 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested 457 by Federov (1988) is used: 658 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by Federov (1988) is used: 458 659 \begin{align} \label{Eq_zdfddm_d} 459 660 A_d^{vT} &= \begin{cases} … … 481 682 % Bottom Friction 482 683 % ================================================================ 483 \section [Bottom Friction (\textit{zdfbfr})] 484 {Bottom Friction (\mdl{zdfbfr} module)} 684 \section [Bottom Friction (\textit{zdfbfr})] {Bottom Friction (\mdl{zdfbfr} module)} 485 685 \label{ZDF_bfr} 486 686 … … 493 693 diffusive flux. For the bottom boundary layer, one has: 494 694 \begin{equation} \label{Eq_zdfbfr_flux} 495 A^{vm} \left( \partial \textbf{U}_h / \partial z \right) = \textbf{F}_h496 \end{equation} 497 where $ \textbf{F}_h$ is supposed to represent the horizontal momentum flux695 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 696 \end{equation} 697 where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum 498 698 outside the logarithmic turbulent boundary layer (thickness of the order of 499 1~m in the ocean). How $ \textbf{F}_h$ influences the interior depends on the699 1~m in the ocean). How ${\cal F}_h^{\textbf U}$ influences the interior depends on the 500 700 vertical resolution of the model near the bottom relative to the Ekman layer 501 701 depth. For example, in order to obtain an Ekman layer depth … … 509 709 this, consider the equation for $u$ at $k$, the last ocean level: 510 710 \begin{equation} \label{Eq_zdfbfr_flux2} 511 \frac{\partial u \; (k)}{\partial t} = \frac{1}{e_{3u}} \left[ A^{vm} \; (k) \frac{U \; (k-1) - U \; (k)}{e_{3uw} \; (k-1)} - F_u \right] \approx - \frac{F_u}{e_{3u}}512 \end{equation} 513 For example, if the bottom layer thickness is 200~m, the Ekman transport will711 \frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}} 712 \end{equation} 713 If the bottom layer thickness is 200~m, the Ekman transport will 514 714 be distributed over that depth. On the other hand, if the vertical resolution 515 715 is high (1~m or less) and a turbulent closure model is used, the turbulent 516 716 Ekman layer will be represented explicitly by the model. However, the 517 717 logarithmic layer is never represented in current primitive equation model 518 applications: it is \emph{necessary} to parameterize the flux $ \textbf{F}_h $.718 applications: it is \emph{necessary} to parameterize the flux ${\cal F}^u_h $. 519 719 Two choices are available in \NEMO: a linear and a quadratic bottom friction. 520 720 Note that in both cases, the rotation between the interior velocity and the 521 721 bottom friction is neglected in the present release of \NEMO. 522 722 723 In the code, the bottom friction is imposed by adding the trend due to the bottom 724 friction to the general momentum trend in \mdl{dynbfr}. For the time-split surface 725 pressure gradient algorithm, the momentum trend due to the barotropic component 726 needs to be handled separately. For this purpose it is convenient to compute and 727 store coefficients which can be simply combined with bottom velocities and geometric 728 values to provide the momentum trend due to bottom friction. 729 These coefficients are computed in \mdl{zdfbfr} and generally take the form 730 $c_b^{\textbf U}$ where: 731 \begin{equation} \label{Eq_zdfbfr_bdef} 732 \frac{\partial {\textbf U_h}}{\partial t} = 733 - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b 734 \end{equation} 735 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 736 523 737 % ------------------------------------------------------------------------------------------------------------- 524 738 % Linear Bottom Friction 525 739 % ------------------------------------------------------------------------------------------------------------- 526 \subsection{Linear Bottom Friction (\np{n botfr} =1) }740 \subsection{Linear Bottom Friction (\np{nn\_botfr} = 0 or 1) } 527 741 \label{ZDF_bfr_linear} 528 742 529 The linear bottom friction parameterisation assumes that the bottom friction 530 is proportional to the interior velocity (i.e. the velocity of the last model level): 743 The linear bottom friction parameterisation (including the special case 744 of a free-slip condition) assumes that the bottom friction 745 is proportional to the interior velocity (i.e. the velocity of the last 746 model level): 531 747 \begin{equation} \label{Eq_zdfbfr_linear} 532 \textbf{F}_h= \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b533 \end{equation} 534 where $ \textbf{U}_h^b$ is the horizontal velocity vector of the bottom ocean535 layer and $r$ is a friction coefficient expressed in m.s$^{-1}$. This coefficient536 is generally estimated by setting a typical decay time$\tau$ in the deep ocean,748 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 749 \end{equation} 750 where $r$ is a friction coefficient expressed in ms$^{-1}$. 751 This coefficient is generally estimated by setting a typical decay time 752 $\tau$ in the deep ocean, 537 753 and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted 538 values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly 1984}.754 values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}. 539 755 A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used 540 756 in quasi-geostrophic models. One may consider the linear friction as an 541 757 approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982}, 542 758 Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed 543 of tidal currents of $U_{av} =0.1$~m .s$^{-1}$, and assuming an ocean depth544 $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m .s$^{-1}$.759 of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$, and assuming an ocean depth 760 $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 545 761 This is the default value used in \NEMO. It corresponds to a decay time scale 546 of 115~days. It can be changed by specifying \np{bfric1} (namelist parameter). 547 548 In the code, the bottom friction is imposed by updating the value of the 549 vertical eddy coefficient at the bottom level. Indeed, the discrete formulation 550 of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that 551 $\textbf {U}_h =0$ below the ocean floor, leads to 552 \begin{equation} \label{Eq_zdfbfr_linKz} 762 of 115~days. It can be changed by specifying \np{rn\_bfric1} (namelist parameter). 763 764 For the linear friction case the coefficients defined in the general 765 expression \eqref{Eq_zdfbfr_bdef} are: 766 \begin{equation} \label{Eq_zdfbfr_linbfr_b} 553 767 \begin{split} 554 A_u^{vm} &= r\;e_{3uw}\\555 A_v^{vm} &= r\;e_{3vw}\\768 c_b^u &= - r\\ 769 c_b^v &= - r\\ 556 770 \end{split} 557 771 \end{equation} 558 559 This update is done in \mdl{zdfbfr} when \np{nbotfr}=1. The value of $r$ 560 used is \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to setting $r=0$ and 561 leads to a free-slip bottom boundary condition. Setting \np{nbotfr}=0 sets 562 $r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the background 563 vertical eddy coefficient, and a no-slip boundary condition is imposed. 564 Note that this latter choice generally leads to an underestimation of the 565 bottom friction: for example with a deepest level thickness of $200~m$ 566 and $A_{vb}^{\rm {\bf U}} =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient 567 is only $r=10^{-6}$m.s$^{-1}$. 772 When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfric1}. 773 Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip 774 bottom boundary condition. These values are assigned in \mdl{zdfbfr}. 775 From v3.2 onwards there is support for local enhancement of these values 776 via an externally defined 2D mask array (\np{ln\_bfr2d}=true) given 777 in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1. 778 Locations with a non-zero mask value will have the friction coefficient increased 779 by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfric1}. 568 780 569 781 % ------------------------------------------------------------------------------------------------------------- 570 782 % Non-Linear Bottom Friction 571 783 % ------------------------------------------------------------------------------------------------------------- 572 \subsection{Non-Linear Bottom Friction (\np{n botfr} = 2)}784 \subsection{Non-Linear Bottom Friction (\np{nn\_botfr} = 2)} 573 785 \label{ZDF_bfr_nonlinear} 574 786 … … 576 788 friction is quadratic: 577 789 \begin{equation} \label{Eq_zdfbfr_nonlinear} 578 \textbf {F}_h= \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h790 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 579 791 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b 580 792 \end{equation} 581 582 with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity ($i.e.$ 583 the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient, 584 and $e_b $ a bottom turbulent kinetic energy due to tides, internal waves 585 breaking and other short time scale currents. A typical value of the drag 586 coefficient is $C_D = 10^{-3} $. As an example, the CME experiment 587 \citep{Treguier1992} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$, 588 while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ 589 and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been 590 set as default values (\np{bfric2} and \np{bfeb2} namelist parameters). 793 where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy 794 due to tides, internal waves breaking and other short time scale currents. 795 A typical value of the drag coefficient is $C_D = 10^{-3} $. As an example, 796 the CME experiment \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and 797 $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} 798 uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. 799 The CME choices have been set as default values (\np{rn\_bfric2} and \np{rn\_bfeb2} 800 namelist parameters). 591 801 592 802 As for the linear case, the bottom friction is imposed in the code by 593 updating the value of the vertical eddy coefficient at the bottom level: 594 \begin{equation} \label{Eq_zdfbfr_nonlinKz} 803 adding the trend due to the bottom friction to the general momentum trend 804 in \mdl{dynbfr}. 805 For the non-linear friction case the terms 806 computed in \mdl{zdfbfr} are: 807 \begin{equation} \label{Eq_zdfbfr_nonlinbfr} 595 808 \begin{split} 596 A_u^{vm} &=C_D\; e_{3uw} \left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^ 597 {1/2}\\ 598 A_v^{vm} &=C_D\; e_{3uw} \left[ \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ 809 c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\ 810 c_b^v &= - \; C_D\;\left[ \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ 599 811 \end{split} 600 812 \end{equation} 601 813 602 This update is done in \mdl{zdfbfr}. The coefficients that control the strength of the 603 non-linear bottom friction are initialized as namelist parameters: $C_D$= \np{bfri2}, 604 and $e_b$ =\np{bfeb2}. 605 606 % ================================================================ 814 The coefficients that control the strength of the non-linear bottom friction are 815 initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. 816 Note for applications which treat tides explicitly a low or even zero value of 817 \np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ 818 is possible via an externally defined 2D mask array (\np{ln\_bfr2d}=true). 819 See previous section for details. 820 821 % ------------------------------------------------------------------------------------------------------------- 822 % Bottom Friction stability 823 % ------------------------------------------------------------------------------------------------------------- 824 \subsection{Bottom Friction stability considerations} 825 \label{ZDF_bfr_stability} 826 827 Some care needs to exercised over the choice of parameters to ensure that the 828 implementation of bottom friction does not induce numerical instability. For 829 the purposes of stability analysis, an approximation to \eqref{Eq_zdfbfr_flux2} 830 is: 831 \begin{equation} \label{Eqn_bfrstab} 832 \begin{split} 833 \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt \\ 834 &= -\frac{ru}{e_{3u}}\;2\rdt\\ 835 \end{split} 836 \end{equation} 837 \noindent where linear bottom friction and a leapfrog timestep have been assumed. 838 To ensure that the bottom friction cannot reverse the direction of flow it is necessary to have: 839 \begin{equation} 840 |\Delta u| < \;|u| 841 \end{equation} 842 \noindent which, using \eqref{Eqn_bfrstab}, gives: 843 \begin{equation} 844 r\frac{2\rdt}{e_{3u}} < 1 \qquad \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ 845 \end{equation} 846 This same inequality can also be derived in the non-linear bottom friction case 847 if a velocity of 1 m.s$^{-1}$ is assumed. Alternatively, this criterion can be 848 rearranged to suggest a minimum bottom box thickness to ensure stability: 849 \begin{equation} 850 e_{3u} > 2\;r\;\rdt 851 \end{equation} 852 \noindent which it may be necessary to impose if partial steps are being used. 853 For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then 854 $e_{3u}$ should be greater than 3.6 m. For most applications, with physically 855 sensible parameters these restrictions should not be of concern. But 856 caution may be necessary if attempts are made to locally enhance the bottom 857 friction parameters. 858 To ensure stability limits are imposed on the bottom friction coefficients both during 859 initialisation and at each time step. Checks at initialisation are made in \mdl{zdfbfr} 860 (assuming a 1 m.s$^{-1}$ velocity in the non-linear case). 861 The number of breaches of the stability criterion are reported as well as the minimum 862 and maximum values that have been set. The criterion is also checked at each time step, 863 using the actual velocity, in \mdl{dynbfr}. Values of the bottom friction coefficient are 864 reduced as necessary to ensure stability; these changes are not reported. 865 866 % ------------------------------------------------------------------------------------------------------------- 867 % Bottom Friction with split-explicit time splitting 868 % ------------------------------------------------------------------------------------------------------------- 869 \subsection{Bottom Friction with split-explicit time splitting} 870 \label{ZDF_bfr_ts} 871 872 When calculating the momentum trend due to bottom friction in \mdl{dynbfr}, the 873 bottom velocity at the before time step is used. This velocity includes both the 874 baroclinic and barotropic components which is appropriate when using either the 875 explicit or filtered surface pressure gradient algorithms (\key{dynspg\_exp} or 876 {\key{dynspg\_flt}). Extra attention is required, however, when using 877 split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface 878 equation is solved with a small time step \np{nn\_baro}*\np{rn\_rdt}, while the three 879 dimensional prognostic variables are solved with a longer time step that is a 880 multiple of \np{rn\_rdt}. The trend in the barotropic momentum due to bottom 881 friction appropriate to this method is that given by the selected parameterisation 882 ($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities 883 at each barotropic timestep. 884 885 In the case of non-linear bottom friction, we have elected to partially linearise 886 the problem by keeping the coefficients fixed throughout the barotropic 887 time-stepping to those computed in \mdl{zdfbfr} using the now timestep. 888 This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to: 889 890 \begin{enumerate} 891 \item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before 892 barotropic velocity to the bottom friction component of the vertically 893 integrated momentum trend. Note the same stability check that is carried out 894 on the bottom friction coefficient in \mdl{dynbfr} has to be applied here to 895 ensure that the trend removed matches that which was added in \mdl{dynbfr}. 896 \item At each barotropic step, compute the contribution of the current barotropic 897 velocity to the trend due to bottom friction. Add this contribution to the 898 vertically integrated momentum trend. This contribution is handled implicitly which 899 eliminates the need to impose a stability criteria on the values of the bottom friction 900 coefficient within the barotropic loop. 901 \end{enumerate} 902 903 Note that the use of an implicit formulation 904 for the bottom friction trend means that any limiting of the bottom friction coefficient 905 in \mdl{dynbfr} does not adversely affect the solution when using split-explicit time 906 splitting. This is because the major contribution to bottom friction is likely to come from 907 the barotropic component which uses the unrestricted value of the coefficient. 908 909 The implicit formulation takes the form: 910 \begin{equation} \label{Eq_zdfbfr_implicitts} 911 \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ] 912 \end{equation} 913 where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height), 914 $c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and $RHS$ represents 915 all the components to the vertically integrated momentum trend except for that due to bottom friction. 916 917 918 919 920 % ================================================================ 921 % Tidal Mixing 922 % ================================================================ 923 \section{Tidal Mixing} 924 \label{ZDF_tmx} 925 926 %--------------------------------------------namzdf_tmx-------------------------------------------------- 927 \namdisplay{namzdf_tmx} 928 %-------------------------------------------------------------------------------------------------------------- 929 930 931 % ------------------------------------------------------------------------------------------------------------- 932 % Bottom intensified tidal mixing 933 % ------------------------------------------------------------------------------------------------------------- 934 \subsection{Bottom intensified tidal mixing} 935 \label{ZDF_tmx_bottom} 936 937 The parameterization of tidal mixing follows the general formulation for 938 the vertical eddy diffusivity proposed by \citet{St_Laurent_al_GRL02} and 939 first introduced in an OGCM by \citep{Simmons_al_OM04}. 940 In this formulation an additional vertical diffusivity resulting from internal tide breaking, 941 $A^{vT}_{tides}$ is expressed as a function of $E(x,y)$, the energy transfer from barotropic 942 tides to baroclinic tides : 943 \begin{equation} \label{Eq_Ktides} 944 A^{vT}_{tides} = q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 } 945 \end{equation} 946 where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency 947 (see \S\ref{TRA_bn2}), $\rho$ the density, $q$ the tidal dissipation efficiency, 948 and $F(z)$ the vertical structure function. 949 950 The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter) 951 and is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980). 952 The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter) 953 represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally, 954 with the remaining $1-q$ radiating away as low mode internal waves and 955 contributing to the background internal wave field. A value of $q=1/3$ is typically used 956 \citet{St_Laurent_al_GRL02}. 957 The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical. 958 It is implemented as a simple exponential decaying upward away from the bottom, 959 with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter, with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04}, 960 \begin{equation} \label{Eq_Fz} 961 F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) } 962 \end{equation} 963 and is normalized so that vertical integral over the water column is unity. 964 965 The associated vertical viscosity is calculated from the vertical 966 diffusivity assuming a Prandtl number of 1, $i.e.$ $A^{vm}_{tides}=A^{vT}_{tides}$. 967 In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity 968 is capped at $300\,cm^2/s$ and impose a lower limit on $N^2$ of \np{rn\_n2min} 969 usually set to $10^{-8} s^{-2}$. These bounds are usually rarely encountered. 970 971 The internal wave energy map, $E(x, y)$ in \eqref{Eq_Ktides}, is derived 972 from a barotropic model of the tides utilizing a parameterization of the 973 conversion of barotropic tidal energy into internal waves. 974 The essential goal of the parameterization is to represent the momentum 975 exchange between the barotropic tides and the unrepresented internal waves 976 induced by the tidal ßow over rough topography in a stratified ocean. 977 In the current version of \NEMO, the map is built from the output of 978 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}. 979 This model provides the dissipation associated with internal wave energy for the M2 and K1 980 tides component (Fig.~\ref{Fig_ZDF_M2_K1_tmx}). The S2 dissipation is simply approximated 981 as being $1/4$ of the M2 one. The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$. 982 Its global mean value is $1.1$ TW, in agreement with independent estimates 983 \citep{Egbert_Ray_Nat00, Egbert_Ray_JGR01}. 984 985 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 986 \begin{figure}[!t] \label{Fig_ZDF_M2_K1_tmx} \begin{center} 987 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_ZDF_M2_K1_tmx.pdf} 988 \caption {(a) M2 and (b) K2 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$). } 989 \end{center} 990 \end{figure} 991 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 992 993 % ------------------------------------------------------------------------------------------------------------- 994 % Indonesian area specific treatment 995 % ------------------------------------------------------------------------------------------------------------- 996 \subsection{Indonesian area specific treatment} 997 \label{ZDF_tmx_itf} 998 999 When the Indonesian Through Flow (ITF) area is included in the model domain, 1000 a specific treatment of tidal induced mixing in this area can be used. 1001 It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide 1002 an input NetCDF file, \ifile{mask\_itf}, which contains a mask array defining the ITF area 1003 where the specific treatment is applied. 1004 1005 When \np{ln\_tmx\_itf}=true, the two key parameters $q$ and $F(z)$ are adjusted following 1006 the parameterisation developed by \ref{Koch-Larrouy_al_GRL07}: 1007 1008 First, the Indonesian archipelago is a complex geographic region 1009 with a series of large, deep, semi-enclosed basins connected via 1010 numerous narrow straits. Once generated, internal tides remain 1011 confined within this semi-enclosed area and hardly radiate away. 1012 Therefore all the internal tides energy is consumed within this area. 1013 So it is assumed that $q = 1$, $i.e.$ all the energy generated is available for mixing. 1014 Note that for test purposed, the ITF tidal dissipation efficiency is a 1015 namelist parameter (\np{rn\_tfe\_itf}). A value of $1$ or close to is 1016 this recommended for this parameter. 1017 1018 Second, the vertical structure function, $F(z)$, is no more associated 1019 with a bottom intensification of the mixing, but with a maximum of 1020 energy available within the thermocline. \ref{Koch-Larrouy_al_GRL07} 1021 have suggested that the vertical distribution of the energy dissipation 1022 proportional to $N^2$ below the core of the thermocline and to $N$ above. 1023 The resulting $F(z)$ is: 1024 \begin{equation} \label{Eq_Fz_itf} 1025 F(i,j,k) \sim \left\{ \begin{aligned} 1026 \frac{q\,\Gamma E(i,j) } {\rho N \, \int N dz} \qquad \text{when $\partial_z N < 0$} \\ 1027 \frac{q\,\Gamma E(i,j) } {\rho \, \int N^2 dz} \qquad \text{when $\partial_z N > 0$} 1028 \end{aligned} \right. 1029 \end{equation} 1030 1031 Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$, 1032 which agrees with the independent estimates inferred from observations. 1033 Introduced in a regional OGCM, the parameterization improves the water mass 1034 characteristics in the different Indonesian seas, suggesting that the horizontal 1035 and vertical distributions of the mixing are adequately prescribed 1036 \citep{Koch-Larrouy_al_GRL07, Koch-Larrouy_al_OD08a, Koch-Larrouy_al_OD08b}. 1037 Note also that such a parameterisation has a sugnificant impact on the behaviour 1038 of global coupled GCMs \citep{Koch-Larrouy_al_CD10}. 1039 1040 1041 % ================================================================
Note: See TracChangeset
for help on using the changeset viewer.
