- Timestamp:
- 2019-08-08T15:40:47+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/fix_vvl_ticket1791/doc
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/fix_vvl_ticket1791/doc
- Property svn:ignore deleted
-
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex
- Property svn:ignore
-
old new 1 *.aux 2 *.bbl 3 *.blg 4 *.dvi 5 *.fdb* 6 *.fls 7 *.idx 8 *.ilg 9 *.ind 10 *.log 11 *.maf 12 *.mtc* 13 *.out 14 *.pdf 15 *.toc 16 _minted-* 1 figures
-
- Property svn:ignore
-
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO
- Property svn:ignore deleted
-
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_ZDF.tex
r10442 r11422 25 25 At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}), 26 26 while at the bottom they are set to zero for heat and salt, 27 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie \ key{trabbl} defined,27 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie \np{ln\_trabbc} defined, 28 28 see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum 29 (see \autoref{sec:ZDF_ bfr}).29 (see \autoref{sec:ZDF_drg}). 30 30 31 31 In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and … … 33 33 respectively (see \autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}). 34 34 These coefficients can be assumed to be either constant, or a function of the local Richardson number, 35 or computed from a turbulent closure model (either TKE or GLS formulation).36 The computation of these coefficients is initialized in the \mdl{zdf ini} module and performed in37 the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} modules.35 or computed from a turbulent closure model (either TKE or GLS or OSMOSIS formulation). 36 The computation of these coefficients is initialized in the \mdl{zdfphy} module and performed in 37 the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules. 38 38 The trends due to the vertical momentum and tracer diffusion, including the surface forcing, 39 39 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 40 These trends can be computed using either a forward time stepping scheme 41 (namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme 42 (\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing coefficients, 43 and thus of the formulation used (see \autoref{chap:STP}). 44 45 % ------------------------------------------------------------------------------------------------------------- 46 % Constant 47 % ------------------------------------------------------------------------------------------------------------- 48 \subsection{Constant (\protect\key{zdfcst})} 49 \label{subsec:ZDF_cst} 50 %--------------------------------------------namzdf--------------------------------------------------------- 40 %These trends can be computed using either a forward time stepping scheme 41 %(namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme 42 %(\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing coefficients, 43 %and thus of the formulation used (see \autoref{chap:STP}). 44 45 %--------------------------------------------namzdf-------------------------------------------------------- 51 46 52 47 \nlst{namzdf} 53 48 %-------------------------------------------------------------------------------------------------------------- 54 49 50 % ------------------------------------------------------------------------------------------------------------- 51 % Constant 52 % ------------------------------------------------------------------------------------------------------------- 53 \subsection[Constant (\forcode{ln_zdfcst = .true.})] 54 {Constant (\protect\np{ln\_zdfcst}\forcode{ = .true.})} 55 \label{subsec:ZDF_cst} 56 55 57 Options are defined through the \ngn{namzdf} namelist variables. 56 When \ key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to58 When \np{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 57 59 constant values over the whole ocean. 58 60 This is the crudest way to define the vertical ocean physics. 59 It is recommended t hat this option is only usedin process studies, not in basin scale simulations.61 It is recommended to use this option only in process studies, not in basin scale simulations. 60 62 Typical values used in this case are: 61 63 \begin{align*} … … 72 74 % Richardson Number Dependent 73 75 % ------------------------------------------------------------------------------------------------------------- 74 \subsection{Richardson number dependent (\protect\key{zdfric})} 76 \subsection[Richardson number dependent (\forcode{ln_zdfric = .true.})] 77 {Richardson number dependent (\protect\np{ln\_zdfric}\forcode{ = .true.})} 75 78 \label{subsec:ZDF_ric} 76 79 … … 80 83 %-------------------------------------------------------------------------------------------------------------- 81 84 82 When \ key{zdfric} is defined, a local Richardson number dependent formulation for the vertical momentum and85 When \np{ln\_zdfric}\forcode{ = .true.}, a local Richardson number dependent formulation for the vertical momentum and 83 86 tracer eddy coefficients is set through the \ngn{namzdf\_ric} namelist variables. 84 87 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. … … 87 90 a dependency between the vertical eddy coefficients and the local Richardson number 88 91 (\ie the ratio of stratification to vertical shear). 89 Following \citet{ Pacanowski_Philander_JPO81}, the following formulation has been implemented:92 Following \citet{pacanowski.philander_JPO81}, the following formulation has been implemented: 90 93 \[ 91 94 % \label{eq:zdfric} … … 124 127 The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}. 125 128 Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to 126 the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{ Lermusiaux2001}.129 the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{lermusiaux_JMS01}. 127 130 128 131 % ------------------------------------------------------------------------------------------------------------- 129 132 % TKE Turbulent Closure Scheme 130 133 % ------------------------------------------------------------------------------------------------------------- 131 \subsection{TKE turbulent closure scheme (\protect\key{zdftke})} 134 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke = .true.})] 135 {TKE turbulent closure scheme (\protect\np{ln\_zdftke}\forcode{ = .true.})} 132 136 \label{subsec:ZDF_tke} 133 134 137 %--------------------------------------------namzdf_tke-------------------------------------------------- 135 138 … … 140 143 a prognostic equation for $\bar{e}$, the turbulent kinetic energy, 141 144 and a closure assumption for the turbulent length scales. 142 This turbulent closure model has been developed by \citet{ Bougeault1989} in the atmospheric case,143 adapted by \citet{ Gaspar1990} for the oceanic case, and embedded in OPA, the ancestor of NEMO,144 by \citet{ Blanke1993} for equatorial Atlantic simulations.145 Since then, significant modifications have been introduced by \citet{ Madec1998} in both the implementation and145 This turbulent closure model has been developed by \citet{bougeault.lacarrere_MWR89} in the atmospheric case, 146 adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of NEMO, 147 by \citet{blanke.delecluse_JPO93} for equatorial Atlantic simulations. 148 Since then, significant modifications have been introduced by \citet{madec.delecluse.ea_NPM98} in both the implementation and 146 149 the formulation of the mixing length scale. 147 150 The time evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical shear, 148 its destruction through stratification, its vertical diffusion, and its dissipation of \citet{ Kolmogorov1942} type:151 its destruction through stratification, its vertical diffusion, and its dissipation of \citet{kolmogorov_IANS42} type: 149 152 \begin{equation} 150 153 \label{eq:zdftke_e} … … 168 171 $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients. 169 172 The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with 170 vertical mixing at any depth \citep{ Gaspar1990}.173 vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}. 171 174 They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}. 172 $P_{rt}$ can be set to unity or, following \citet{ Blanke1993}, be a function of the local Richardson number, $R_i$:175 $P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$: 173 176 \begin{align*} 174 177 % \label{eq:prt} … … 180 183 \end{cases} 181 184 \end{align*} 182 Options are defined through the \ngn{namzdfy\_tke} namelist variables.183 185 The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable. 184 186 185 187 At the sea surface, the value of $\bar{e}$ is prescribed from the wind stress field as 186 188 $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn\_ebb} namelist parameter. 187 The default value of $e_{bb}$ is 3.75. \citep{ Gaspar1990}), however a much larger value can be used when189 The default value of $e_{bb}$ is 3.75. \citep{gaspar.gregoris.ea_JGR90}), however a much larger value can be used when 188 190 taking into account the surface wave breaking (see below Eq. \autoref{eq:ZDF_Esbc}). 189 191 The bottom value of TKE is assumed to be equal to the value of the level just above. … … 191 193 the numerical scheme does not ensure its positivity. 192 194 To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} namelist parameter). 193 Following \citet{ Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$.194 This allows the subsequent formulations to match that of \citet{ Gargett1984} for the diffusion in195 Following \citet{gaspar.gregoris.ea_JGR90}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 196 This allows the subsequent formulations to match that of \citet{gargett_JMR84} for the diffusion in 195 197 the thermocline and deep ocean : $K_\rho = 10^{-3} / N$. 196 198 In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical instabilities associated with 197 199 too weak vertical diffusion. 198 200 They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and 199 \np{rn\_avt0} ( namzdfnamelist, see \autoref{subsec:ZDF_cst}).201 \np{rn\_avt0} (\ngn{namzdf} namelist, see \autoref{subsec:ZDF_cst}). 200 202 201 203 \subsubsection{Turbulent length scale} 202 204 203 205 For computational efficiency, the original formulation of the turbulent length scales proposed by 204 \citet{ Gaspar1990} has been simplified.206 \citet{gaspar.gregoris.ea_JGR90} has been simplified. 205 207 Four formulations are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist parameter. 206 The first two are based on the following first order approximation \citep{ Blanke1993}:208 The first two are based on the following first order approximation \citep{blanke.delecluse_JPO93}: 207 209 \begin{equation} 208 210 \label{eq:tke_mxl0_1} … … 212 214 The resulting length scale is bounded by the distance to the surface or to the bottom 213 215 (\np{nn\_mxl}\forcode{ = 0}) or by the local vertical scale factor (\np{nn\_mxl}\forcode{ = 1}). 214 \citet{ Blanke1993} notice that this simplification has two major drawbacks:216 \citet{blanke.delecluse_JPO93} notice that this simplification has two major drawbacks: 215 217 it makes no sense for locally unstable stratification and the computation no longer uses all 216 218 the information contained in the vertical density profile. 217 To overcome these drawbacks, \citet{ Madec1998} introduces the \np{nn\_mxl}\forcode{ = 2..3} cases,219 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{ = 2, 3} cases, 218 220 which add an extra assumption concerning the vertical gradient of the computed length scale. 219 221 So, the length scales are first evaluated as in \autoref{eq:tke_mxl0_1} and then bounded such that: … … 225 227 \autoref{eq:tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than 226 228 the variations of depth. 227 It provides a better approximation of the \citet{ Gaspar1990} formulation while being much less229 It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less 228 230 time consuming. 229 231 In particular, it allows the length scale to be limited not only by the distance to the surface or … … 237 239 \begin{figure}[!t] 238 240 \begin{center} 239 \includegraphics[width= 1.00\textwidth]{Fig_mixing_length}241 \includegraphics[width=\textwidth]{Fig_mixing_length} 240 242 \caption{ 241 243 \protect\label{fig:mixing_length} … … 258 260 In the \np{nn\_mxl}\forcode{ = 2} case, the dissipation and mixing length scales take the same value: 259 261 $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np{nn\_mxl}\forcode{ = 3} case, 260 the dissipation and mixing turbulent length scales are give as in \citet{ Gaspar1990}:262 the dissipation and mixing turbulent length scales are give as in \citet{gaspar.gregoris.ea_JGR90}: 261 263 \[ 262 264 % \label{eq:tke_mxl_gaspar} … … 270 272 Usually the surface scale is given by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant and 271 273 $z_o$ the roughness parameter of the surface. 272 Assuming $z_o=0.1$~m \citep{ Craig_Banner_JPO94} leads to a 0.04~m, the default value of \np{rn\_mxl0}.274 Assuming $z_o=0.1$~m \citep{craig.banner_JPO94} leads to a 0.04~m, the default value of \np{rn\_mxl0}. 273 275 In the ocean interior a minimum length scale is set to recover the molecular viscosity when 274 276 $\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ). … … 277 279 %-----------------------------------------------------------------------% 278 280 279 Following \citet{ Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified to281 Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to 280 282 include the effect of surface wave breaking energetics. 281 283 This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow. 282 The \citet{ Mellor_Blumberg_JPO04} modifications acts on surface length scale and TKE values and284 The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and 283 285 air-sea drag coefficient. 284 The latter concerns the bulk formul eaand is not discussed here.285 286 Following \citet{ Craig_Banner_JPO94}, the boundary condition on surface TKE value is :286 The latter concerns the bulk formulae and is not discussed here. 287 288 Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is : 287 289 \begin{equation} 288 290 \label{eq:ZDF_Esbc} 289 291 \bar{e}_o = \frac{1}{2}\,\left( 15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o} 290 292 \end{equation} 291 where $\alpha_{CB}$ is the \citet{ Craig_Banner_JPO94} constant of proportionality which depends on the ''wave age'',292 ranging from 57 for mature waves to 146 for younger waves \citep{ Mellor_Blumberg_JPO04}.293 where $\alpha_{CB}$ is the \citet{craig.banner_JPO94} constant of proportionality which depends on the ''wave age'', 294 ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}. 293 295 The boundary condition on the turbulent length scale follows the Charnock's relation: 294 296 \begin{equation} … … 297 299 \end{equation} 298 300 where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant. 299 \citet{ Mellor_Blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by300 \citet{ Stacey_JPO99} citing observation evidence, and301 \citet{mellor.blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by 302 \citet{stacey_JPO99} citing observation evidence, and 301 303 $\alpha_{CB} = 100$ the Craig and Banner's value. 302 304 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 303 305 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds 304 306 to $\alpha_{CB} = 100$. 305 Further setting \np{ln\_mxl0 } to true applies \autoref{eq:ZDF_Lsbc} as surface boundary condition onlength scale,307 Further setting \np{ln\_mxl0=.true.}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale, 306 308 with $\beta$ hard coded to the Stacey's value. 307 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on 309 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the 308 310 surface $\bar{e}$ value. 309 311 … … 315 317 Although LC have nothing to do with convection, the circulation pattern is rather similar to 316 318 so-called convective rolls in the atmospheric boundary layer. 317 The detailed physics behind LC is described in, for example, \citet{ Craik_Leibovich_JFM76}.319 The detailed physics behind LC is described in, for example, \citet{craik.leibovich_JFM76}. 318 320 The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and 319 321 wind drift currents. 320 322 321 323 Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by 322 \citep{ Axell_JGR02} for a $k-\epsilon$ turbulent closure.324 \citep{axell_JGR02} for a $k-\epsilon$ turbulent closure. 323 325 The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in 324 an extra source term sof TKE, $P_{LC}$.326 an extra source term of TKE, $P_{LC}$. 325 327 The presence of $P_{LC}$ in \autoref{eq:zdftke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to 326 \forcode{.true.} in the namtkenamelist.328 \forcode{.true.} in the \ngn{namzdf\_tke} namelist. 327 329 328 By making an analogy with the characteristic convective velocity scale (\eg, \citet{ D'Alessio_al_JPO98}),330 By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}), 329 331 $P_{LC}$ is assumed to be : 330 332 \[ … … 334 336 With no information about the wave field, $w_{LC}$ is assumed to be proportional to 335 337 the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module 336 \footnote{Following \citet{ Li_Garrett_JMR93}, the surface Stoke drift velocity may be expressed as338 \footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as 337 339 $u_s = 0.016 \,|U_{10m}|$. 338 340 Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of … … 350 352 \end{cases} 351 353 \] 352 where $c_{LC} = 0.15$ has been chosen by \citep{ Axell_JGR02} as a good compromise to fit LES data.354 where $c_{LC} = 0.15$ has been chosen by \citep{axell_JGR02} as a good compromise to fit LES data. 353 355 The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second. 354 356 The value of $c_{LC}$ is set through the \np{rn\_lc} namelist parameter, 355 having in mind that it should stay between 0.15 and 0.54 \citep{ Axell_JGR02}.357 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}. 356 358 357 359 The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations: 358 $H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by360 $H_{LC}$ is the depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by 359 361 converting its kinetic energy to potential energy, according to 360 362 \[ … … 368 370 produce mixed-layer depths that are too shallow during summer months and windy conditions. 369 371 This bias is particularly acute over the Southern Ocean. 370 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{ Rodgers_2014}.372 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}. 371 373 The parameterization is an empirical one, \ie not derived from theoretical considerations, 372 374 but rather is meant to account for observed processes that affect the density structure of … … 383 385 \end{equation} 384 386 where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that 385 penetrate in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of387 penetrates in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of 386 388 the penetration, and $f_i$ is the ice concentration 387 (no penetration if $f_i=1$, that isif the ocean is entirely covered by sea-ice).389 (no penetration if $f_i=1$, \ie if the ocean is entirely covered by sea-ice). 388 390 The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. 389 391 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}\forcode{ = 0}) or … … 391 393 (\np{nn\_etau}\forcode{ = 1}). 392 394 393 Note that two other option exist e, \np{nn\_etau}\forcode{ = 2..3}.395 Note that two other option exist, \np{nn\_etau}\forcode{ = 2, 3}. 394 396 They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer, 395 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrate the ocean.397 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean. 396 398 Those two options are obsolescent features introduced for test purposes. 397 399 They will be removed in the next release. 400 401 % This should be explain better below what this rn_eice parameter is meant for: 402 In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn\_eice} namelist parameter. 403 This parameter varies from \forcode{0} for no effect to \forcode{4} to suppress the TKE input into the ocean when Sea Ice concentration 404 is greater than 25\%. 398 405 399 406 % from Burchard et al OM 2008 : … … 406 413 407 414 % ------------------------------------------------------------------------------------------------------------- 408 % TKE discretization considerations 409 % ------------------------------------------------------------------------------------------------------------- 410 \subsection{TKE discretization considerations (\protect\key{zdftke})} 415 % GLS Generic Length Scale Scheme 416 % ------------------------------------------------------------------------------------------------------------- 417 \subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls = .true.})] 418 {GLS: Generic Length Scale (\protect\np{ln\_zdfgls}\forcode{ = .true.})} 419 \label{subsec:ZDF_gls} 420 421 %--------------------------------------------namzdf_gls--------------------------------------------------------- 422 423 \nlst{namzdf_gls} 424 %-------------------------------------------------------------------------------------------------------------- 425 426 The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations: 427 one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale, 428 $\psi$ \citep{umlauf.burchard_JMR03, umlauf.burchard_CSR05}. 429 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 430 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:GLS} allows to recover a number of 431 well-known turbulent closures ($k$-$kl$ \citep{mellor.yamada_RG82}, $k$-$\epsilon$ \citep{rodi_JGR87}, 432 $k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}). 433 The GLS scheme is given by the following set of equations: 434 \begin{equation} 435 \label{eq:zdfgls_e} 436 \frac{\partial \bar{e}}{\partial t} = 437 \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 438 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 439 -K_\rho \,N^2 440 +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right] 441 - \epsilon 442 \end{equation} 443 444 \[ 445 % \label{eq:zdfgls_psi} 446 \begin{split} 447 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ 448 \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 449 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 450 - C_3 \,K_\rho\,N^2 - C_2 \,\epsilon \,Fw \right\} \\ 451 &+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 } 452 \;\frac{\partial \psi}{\partial k}} \right]\; 453 \end{split} 454 \] 455 456 \[ 457 % \label{eq:zdfgls_kz} 458 \begin{split} 459 K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ 460 K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l 461 \end{split} 462 \] 463 464 \[ 465 % \label{eq:zdfgls_eps} 466 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; 467 \] 468 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and 469 $\epsilon$ the dissipation rate. 470 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of 471 the choice of the turbulence model. 472 Four different turbulent models are pre-defined (\autoref{tab:GLS}). 473 They are made available through the \np{nn\_clo} namelist parameter. 474 475 %--------------------------------------------------TABLE-------------------------------------------------- 476 \begin{table}[htbp] 477 \begin{center} 478 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} 479 \begin{tabular}{ccccc} 480 & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ 481 % & \citep{mellor.yamada_RG82} & \citep{rodi_JGR87} & \citep{wilcox_AJ88} & \\ 482 \hline 483 \hline 484 \np{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ 485 \hline 486 $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ 487 $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ 488 $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ 489 $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ 490 $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ 491 $C_3$ & 1. & 1. & 1. & 1. \\ 492 $F_{wall}$ & Yes & -- & -- & -- \\ 493 \hline 494 \hline 495 \end{tabular} 496 \caption{ 497 \protect\label{tab:GLS} 498 Set of predefined GLS parameters, or equivalently predefined turbulence models available with 499 \protect\np{ln\_zdfgls}\forcode{ = .true.} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls}. 500 } 501 \end{center} 502 \end{table} 503 %-------------------------------------------------------------------------------------------------------------- 504 505 In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of 506 the mixing length towards $\kappa z_b$ ($\kappa$ is the Von Karman constant and $z_b$ the rugosity length scale) value near physical boundaries 507 (logarithmic boundary layer law). 508 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88}, 509 or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01} 510 (\np{nn\_stab\_func}\forcode{ = 0, 3}, resp.). 511 The value of $C_{0\mu}$ depends on the choice of the stability function. 512 513 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or 514 Neumann condition through \np{nn\_bc\_surf} and \np{nn\_bc\_bot}, resp. 515 As for TKE closure, the wave effect on the mixing is considered when 516 \np{rn\_crban}\forcode{ > 0.} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}. 517 The \np{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and 518 \np{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}. 519 520 The $\psi$ equation is known to fail in stably stratified flows, and for this reason 521 almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy. 522 With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. 523 A value of $c_{lim} = 0.53$ is often used \citep{galperin.kantha.ea_JAS88}. 524 \cite{umlauf.burchard_CSR05} show that the value of the clipping factor is of crucial importance for 525 the entrainment depth predicted in stably stratified situations, 526 and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes. 527 The clipping is only activated if \np{ln\_length\_lim}\forcode{ = .true.}, 528 and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value. 529 530 The time and space discretization of the GLS equations follows the same energetic consideration as for 531 the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}. 532 Evaluation of the 4 GLS turbulent closure schemes can be found in \citet{warner.sherwood.ea_OM05} in ROMS model and 533 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO model. 534 535 536 % ------------------------------------------------------------------------------------------------------------- 537 % OSM OSMOSIS BL Scheme 538 % ------------------------------------------------------------------------------------------------------------- 539 \subsection[OSM: OSMosis boundary layer scheme (\forcode{ln_zdfosm = .true.})] 540 {OSM: OSMosis boundary layer scheme (\protect\np{ln\_zdfosm}\forcode{ = .true.})} 541 \label{subsec:ZDF_osm} 542 %--------------------------------------------namzdf_osm--------------------------------------------------------- 543 544 \nlst{namzdf_osm} 545 %-------------------------------------------------------------------------------------------------------------- 546 547 The OSMOSIS turbulent closure scheme is based on...... TBC 548 549 % ------------------------------------------------------------------------------------------------------------- 550 % TKE and GLS discretization considerations 551 % ------------------------------------------------------------------------------------------------------------- 552 \subsection[ Discrete energy conservation for TKE and GLS schemes] 553 {Discrete energy conservation for TKE and GLS schemes} 411 554 \label{subsec:ZDF_tke_ene} 412 555 … … 414 557 \begin{figure}[!t] 415 558 \begin{center} 416 \includegraphics[width= 1.00\textwidth]{Fig_ZDF_TKE_time_scheme}559 \includegraphics[width=\textwidth]{Fig_ZDF_TKE_time_scheme} 417 560 \caption{ 418 561 \protect\label{fig:TKE_time_scheme} 419 Illustration of the TKE time integrationand its links to the momentum and tracer time integration.562 Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and its links to the momentum and tracer time integration. 420 563 } 421 564 \end{center} … … 424 567 425 568 The production of turbulence by vertical shear (the first term of the right hand side of 426 \autoref{eq:zdftke_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion569 \autoref{eq:zdftke_e}) and \autoref{eq:zdfgls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion 427 570 (first line in \autoref{eq:PE_zdf}). 428 To do so a special care ha veto be taken for both the time and space discretization of429 the TKE equation \citep{Burchard_OM02,Marsaleix_al_OM08}.571 To do so a special care has to be taken for both the time and space discretization of 572 the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}. 430 573 431 574 Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows how 432 575 the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with 433 the one-level forward time stepping of TKE equation.576 the one-level forward time stepping of the equation for $\bar{e}$. 434 577 With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to 435 578 the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and … … 456 599 457 600 A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 458 (second term of the right hand side of \autoref{eq:zdftke_e} ).601 (second term of the right hand side of \autoref{eq:zdftke_e} and \autoref{eq:zdfgls_e}). 459 602 This term must balance the input of potential energy resulting from vertical mixing. 460 The rate of change of potential energy (in 1D for the demonstration) due vertical mixing is obtained by461 multiplying vertical density diffusion tendency by $g\,z$ and and summing the result vertically:603 The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by 604 multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically: 462 605 \begin{equation} 463 606 \label{eq:energ2} … … 475 618 The second term is minus the destruction rate of $\bar{e}$ due to stratification. 476 619 Therefore \autoref{eq:energ1} implies that, to be energetically consistent, 477 the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e} , the TKE equation.620 the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e} and \autoref{eq:zdfgls_e}. 478 621 479 622 Let us now address the space discretization issue. … … 483 626 By redoing the \autoref{eq:energ1} in the 3D case, it can be shown that the product of eddy coefficient by 484 627 the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. 485 Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have tobe taken into account.628 Furthermore, the time variation of $e_3$ has be taken into account. 486 629 487 630 The above energetic considerations leads to the following final discrete form for the TKE equation: … … 507 650 are time stepped using a backward scheme (see\autoref{sec:STP_forward_imp}). 508 651 Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible. 509 The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as 510 they all appear in the right hand side of \autoref{eq:zdftke_ene}. 511 For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 512 513 % ------------------------------------------------------------------------------------------------------------- 514 % GLS Generic Length Scale Scheme 515 % ------------------------------------------------------------------------------------------------------------- 516 \subsection{GLS: Generic Length Scale (\protect\key{zdfgls})} 517 \label{subsec:ZDF_gls} 518 519 %--------------------------------------------namzdf_gls--------------------------------------------------------- 520 521 \nlst{namzdf_gls} 522 %-------------------------------------------------------------------------------------------------------------- 523 524 The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations: 525 one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale, 526 $\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}. 527 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 528 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:GLS} allows to recover a number of 529 well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982}, $k$-$\epsilon$ \citep{Rodi_1987}, 530 $k$-$\omega$ \citep{Wilcox_1988} among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}). 531 The GLS scheme is given by the following set of equations: 532 \begin{equation} 533 \label{eq:zdfgls_e} 534 \frac{\partial \bar{e}}{\partial t} = 535 \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 536 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 537 -K_\rho \,N^2 538 +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right] 539 - \epsilon 540 \end{equation} 541 542 \[ 543 % \label{eq:zdfgls_psi} 544 \begin{split} 545 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ 546 \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 547 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 548 - C_3 \,K_\rho\,N^2 - C_2 \,\epsilon \,Fw \right\} \\ 549 &+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 } 550 \;\frac{\partial \psi}{\partial k}} \right]\; 551 \end{split} 552 \] 553 554 \[ 555 % \label{eq:zdfgls_kz} 556 \begin{split} 557 K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ 558 K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l 559 \end{split} 560 \] 561 562 \[ 563 % \label{eq:zdfgls_eps} 564 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; 565 \] 566 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and 567 $\epsilon$ the dissipation rate. 568 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of 569 the choice of the turbulence model. 570 Four different turbulent models are pre-defined (Tab.\autoref{tab:GLS}). 571 They are made available through the \np{nn\_clo} namelist parameter. 572 573 %--------------------------------------------------TABLE-------------------------------------------------- 574 \begin{table}[htbp] 575 \begin{center} 576 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} 577 \begin{tabular}{ccccc} 578 & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ 579 % & \citep{Mellor_Yamada_1982} & \citep{Rodi_1987} & \citep{Wilcox_1988} & \\ 580 \hline 581 \hline 582 \np{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ 583 \hline 584 $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ 585 $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ 586 $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ 587 $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ 588 $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ 589 $C_3$ & 1. & 1. & 1. & 1. \\ 590 $F_{wall}$ & Yes & -- & -- & -- \\ 591 \hline 592 \hline 593 \end{tabular} 594 \caption{ 595 \protect\label{tab:GLS} 596 Set of predefined GLS parameters, or equivalently predefined turbulence models available with 597 \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls}. 598 } 599 \end{center} 600 \end{table} 601 %-------------------------------------------------------------------------------------------------------------- 602 603 In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of 604 the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length) value near physical boundaries 605 (logarithmic boundary layer law). 606 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{Galperin_al_JAS88}, 607 or by \citet{Kantha_Clayson_1994} or one of the two functions suggested by \citet{Canuto_2001} 608 (\np{nn\_stab\_func}\forcode{ = 0..3}, resp.). 609 The value of $C_{0\mu}$ depends of the choice of the stability function. 610 611 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or 612 Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp. 613 As for TKE closure, the wave effect on the mixing is considered when 614 \np{ln\_crban}\forcode{ = .true.} \citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}. 615 The \np{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and 616 \np{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}. 617 618 The $\psi$ equation is known to fail in stably stratified flows, and for this reason 619 almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy. 620 With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. 621 A value of $c_{lim} = 0.53$ is often used \citep{Galperin_al_JAS88}. 622 \cite{Umlauf_Burchard_CSR05} show that the value of the clipping factor is of crucial importance for 623 the entrainment depth predicted in stably stratified situations, 624 and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes. 625 The clipping is only activated if \np{ln\_length\_lim}\forcode{ = .true.}, 626 and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value. 627 628 The time and space discretization of the GLS equations follows the same energetic consideration as for 629 the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{Burchard_OM02}. 630 Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}. 631 632 % ------------------------------------------------------------------------------------------------------------- 633 % OSM OSMOSIS BL Scheme 634 % ------------------------------------------------------------------------------------------------------------- 635 \subsection{OSM: OSMOSIS boundary layer scheme (\protect\key{zdfosm})} 636 \label{subsec:ZDF_osm} 637 638 %--------------------------------------------namzdf_osm--------------------------------------------------------- 639 640 \nlst{namzdf_osm} 641 %-------------------------------------------------------------------------------------------------------------- 642 643 The OSMOSIS turbulent closure scheme is based on...... TBC 652 %The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as 653 %they all appear in the right hand side of \autoref{eq:zdftke_ene}. 654 %For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 644 655 645 656 % ================================================================ … … 648 659 \section{Convection} 649 660 \label{sec:ZDF_conv} 650 651 %--------------------------------------------namzdf--------------------------------------------------------652 653 \nlst{namzdf}654 %--------------------------------------------------------------------------------------------------------------655 661 656 662 Static instabilities (\ie light potential densities under heavy ones) may occur at particular ocean grid points. … … 664 670 % Non-Penetrative Convective Adjustment 665 671 % ------------------------------------------------------------------------------------------------------------- 666 \subsection[Non-penetrative convective adj mt (\protect\np{ln\_tranpc}\forcode{= .true.})]667 672 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc = .true.})] 673 {Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})} 668 674 \label{subsec:ZDF_npc} 669 670 %--------------------------------------------namzdf--------------------------------------------------------671 672 \nlst{namzdf}673 %--------------------------------------------------------------------------------------------------------------674 675 675 676 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 676 677 \begin{figure}[!htb] 677 678 \begin{center} 678 \includegraphics[width= 0.90\textwidth]{Fig_npc}679 \includegraphics[width=\textwidth]{Fig_npc} 679 680 \caption{ 680 681 \protect\label{fig:npc} … … 700 701 the water column, but only until the density structure becomes neutrally stable 701 702 (\ie until the mixed portion of the water column has \textit{exactly} the density of the water just below) 702 \citep{ Madec_al_JPO91}.703 \citep{madec.delecluse.ea_JPO91}. 703 704 The associated algorithm is an iterative process used in the following way (\autoref{fig:npc}): 704 705 starting from the top of the ocean, the first instability is found. … … 718 719 the algorithm used in \NEMO converges for any profile in a number of iterations which is less than 719 720 the number of vertical levels. 720 This property is of paramount importance as pointed out by \citet{ Killworth1989}:721 This property is of paramount importance as pointed out by \citet{killworth_iprc89}: 721 722 it avoids the existence of permanent and unrealistic static instabilities at the sea surface. 722 723 This non-penetrative convective algorithm has been proved successful in studies of the deep water formation in 723 the north-western Mediterranean Sea \citep{ Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.724 the north-western Mediterranean Sea \citep{madec.delecluse.ea_JPO91, madec.chartier.ea_DAO91, madec.crepon_iprc91}. 724 725 725 726 The current implementation has been modified in order to deal with any non linear equation of seawater … … 727 728 Two main differences have been introduced compared to the original algorithm: 728 729 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 729 (not the thedifference in potential density);730 (not the difference in potential density); 730 731 $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in 731 732 the same way their temperature and salinity has been mixed. … … 736 737 % Enhanced Vertical Diffusion 737 738 % ------------------------------------------------------------------------------------------------------------- 738 \subsection{Enhanced vertical diffusion (\protect\np{ln\_zdfevd}\forcode{ = .true.})} 739 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd = .true.})] 740 {Enhanced vertical diffusion (\protect\np{ln\_zdfevd}\forcode{ = .true.})} 739 741 \label{subsec:ZDF_evd} 740 741 %--------------------------------------------namzdf--------------------------------------------------------742 743 \nlst{namzdf}744 %--------------------------------------------------------------------------------------------------------------745 742 746 743 Options are defined through the \ngn{namzdf} namelist variables. 747 744 The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}\forcode{ = .true.}. 748 In this case, the vertical eddy mixing coefficients are assigned very large values 749 (a typical value is $10\;m^2s^{-1})$in regions where the stratification is unstable750 (\ie when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{ Lazar_PhD97, Lazar_al_JPO99}.745 In this case, the vertical eddy mixing coefficients are assigned very large values 746 in regions where the stratification is unstable 747 (\ie when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}. 751 748 This is done either on tracers only (\np{nn\_evdm}\forcode{ = 0}) or 752 749 on both momentum and tracers (\np{nn\_evdm}\forcode{ = 1}). … … 759 756 the convective adjustment algorithm presented above when mixing both tracers and 760 757 momentum in the case of static instabilities. 761 It requires the use of an implicit time stepping on vertical diffusion terms762 (\ie np{ln\_zdfexp}\forcode{ = .false.}).763 758 764 759 Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$. 765 760 This removes a potential source of divergence of odd and even time step in 766 a leapfrog environment \citep{ Leclair_PhD2010} (see \autoref{sec:STP_mLF}).761 a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:STP_mLF}). 767 762 768 763 % ------------------------------------------------------------------------------------------------------------- 769 764 % Turbulent Closure Scheme 770 765 % ------------------------------------------------------------------------------------------------------------- 771 \subsection[Turbulent closure scheme (\protect\key{zdf}\{tke,gls,osm\})]{Turbulent Closure Scheme (\protect\key{zdftke}, \protect\key{zdfgls} or \protect\key{zdfosm})} 766 \subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf/tke/gls/osm = .true.})] 767 {Handling convection with turbulent closure schemes (\protect\np{ln\_zdf/tke/gls/osm}\forcode{ = .true.})} 772 768 \label{subsec:ZDF_tcs} 773 769 774 The turbulent closure scheme presented in \autoref{subsec:ZDF_tke} and \autoref{subsec:ZDF_gls} 775 (\key{zdftke} or \key{zdftke} is defined) in theory solves the problem of statically unstable density profiles. 770 771 The turbulent closure schemes presented in \autoref{subsec:ZDF_tke}, \autoref{subsec:ZDF_gls} and 772 \autoref{subsec:ZDF_osm} (\ie \np{ln\_zdftke} or \np{ln\_zdfgls} or \np{ln\_zdfosm} defined) deal, in theory, 773 with statically unstable density profiles. 776 774 In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in 777 775 \autoref{eq:zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative. 778 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also the four neighbouring $A_u^{vm} {and}\;A_v^{vm}$779 (up to $1\;m^2s^{-1}$).776 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also of the four neighboring values at 777 velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$). 780 778 These large values restore the static stability of the water column in a way similar to that of 781 779 the enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}). … … 785 783 It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme, 786 784 \ie setting the \np{ln\_zdfnpc} namelist parameter to true and 787 defining the turbulent closure CPP keyall together.788 789 The KPPturbulent closure scheme already includes enhanced vertical diffusion in the case of convection,790 as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp},791 therefore \np{ln\_zdfevd}\forcode{ = .false.} should be used with the KPPscheme.785 defining the turbulent closure (\np{ln\_zdftke} or \np{ln\_zdfgls} = \forcode{.true.}) all together. 786 787 The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection, 788 %as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp}, 789 therefore \np{ln\_zdfevd}\forcode{ = .false.} should be used with the OSMOSIS scheme. 792 790 % gm% + one word on non local flux with KPP scheme trakpp.F90 module... 793 791 … … 795 793 % Double Diffusion Mixing 796 794 % ================================================================ 797 \section{Double diffusion mixing (\protect\key{zdfddm})} 798 \label{sec:ZDF_ddm} 795 \section[Double diffusion mixing (\forcode{ln_zdfddm = .true.})] 796 {Double diffusion mixing (\protect\np{ln\_zdfddm}\forcode{ = .true.})} 797 \label{subsec:ZDF_ddm} 798 799 799 800 800 %-------------------------------------------namzdf_ddm------------------------------------------------- … … 803 803 %-------------------------------------------------------------------------------------------------------------- 804 804 805 Options are defined through the \ngn{namzdf\_ddm} namelist variables. 805 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter 806 \np{ln\_zdfddm} in \ngn{namzdf}. 806 807 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. 807 808 The former condition leads to salt fingering and the latter to diffusive convection. 808 809 Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the ocean. 809 \citet{ Merryfield1999} include a parameterisation of such phenomena in a global ocean model and show that810 \citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that 810 811 it leads to relatively minor changes in circulation but exerts significant regional influences on 811 812 temperature and salinity. 812 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the \key{zdfddm} CPP key. 813 813 814 814 815 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients … … 839 840 \begin{figure}[!t] 840 841 \begin{center} 841 \includegraphics[width= 0.99\textwidth]{Fig_zdfddm}842 \includegraphics[width=\textwidth]{Fig_zdfddm} 842 843 \caption{ 843 844 \protect\label{fig:zdfddm} 844 From \citet{ Merryfield1999} :845 From \citet{merryfield.holloway.ea_JPO99} : 845 846 (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in regions of salt fingering. 846 847 Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$; … … 855 856 856 857 The factor 0.7 in \autoref{eq:zdfddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx 0.7$ of 857 buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{ McDougall_Taylor_JMR84}).858 Following \citet{ Merryfield1999}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.858 buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{mcdougall.taylor_JMR84}). 859 Following \citet{merryfield.holloway.ea_JPO99}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 859 860 860 861 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by … … 887 888 % Bottom Friction 888 889 % ================================================================ 889 \section{Bottom and top friction (\protect\mdl{zdfbfr})} 890 \label{sec:ZDF_bfr} 890 \section[Bottom and top friction (\textit{zdfdrg.F90})] 891 {Bottom and top friction (\protect\mdl{zdfdrg})} 892 \label{sec:ZDF_drg} 891 893 892 894 %--------------------------------------------nambfr-------------------------------------------------------- 893 895 % 894 %\nlst{nambfr} 896 \nlst{namdrg} 897 \nlst{namdrg_top} 898 \nlst{namdrg_bot} 899 895 900 %-------------------------------------------------------------------------------------------------------------- 896 901 897 Options to define the top and bottom friction are defined through the \ngn{nam bfr} namelist variables.902 Options to define the top and bottom friction are defined through the \ngn{namdrg} namelist variables. 898 903 The bottom friction represents the friction generated by the bathymetry. 899 904 The top friction represents the friction generated by the ice shelf/ocean interface. 900 As the friction processes at the top and bottom are treated in similar way,901 only the bottom friction is described in detail below.905 As the friction processes at the top and the bottom are treated in and identical way, 906 the description below considers mostly the bottom friction case, if not stated otherwise. 902 907 903 908 … … 905 910 a condition on the vertical diffusive flux. 906 911 For the bottom boundary layer, one has: 907 \[908 % \label{eq:zdfbfr_flux}909 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}910 \]912 \[ 913 % \label{eq:zdfbfr_flux} 914 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 915 \] 911 916 where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum outside 912 917 the logarithmic turbulent boundary layer (thickness of the order of 1~m in the ocean). … … 922 927 To illustrate this, consider the equation for $u$ at $k$, the last ocean level: 923 928 \begin{equation} 924 \label{eq:zdf bfr_flux2}929 \label{eq:zdfdrg_flux2} 925 930 \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}} 926 931 \end{equation} … … 935 940 936 941 In the code, the bottom friction is imposed by adding the trend due to the bottom friction to 937 the general momentum trend in \mdl{dynbfr}.942 the general momentum trend in \mdl{dynzdf}. 938 943 For the time-split surface pressure gradient algorithm, the momentum trend due to 939 944 the barotropic component needs to be handled separately. 940 945 For this purpose it is convenient to compute and store coefficients which can be simply combined with 941 946 bottom velocities and geometric values to provide the momentum trend due to bottom friction. 942 These coefficients are computed in \mdl{zdfbfr} and generally take the form $c_b^{\textbf U}$ where:947 These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where: 943 948 \begin{equation} 944 949 \label{eq:zdfbfr_bdef} … … 946 951 - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b 947 952 \end{equation} 948 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 953 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 954 Note than from \NEMO 4.0, drag coefficients are only computed at cell centers (\ie at T-points) and refer to as $c_b^T$ in the following. These are then linearly interpolated in space to get $c_b^\textbf{U}$ at velocity points. 949 955 950 956 % ------------------------------------------------------------------------------------------------------------- 951 957 % Linear Bottom Friction 952 958 % ------------------------------------------------------------------------------------------------------------- 953 \subsection{Linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 0..1})} 954 \label{subsec:ZDF_bfr_linear} 955 956 The linear bottom friction parameterisation (including the special case of a free-slip condition) assumes that 957 the bottom friction is proportional to the interior velocity (\ie the velocity of the last model level): 959 \subsection[Linear top/bottom friction (\forcode{ln_lin = .true.})] 960 {Linear top/bottom friction (\protect\np{ln\_lin}\forcode{ = .true.)}} 961 \label{subsec:ZDF_drg_linear} 962 963 The linear friction parameterisation (including the special case of a free-slip condition) assumes that 964 the friction is proportional to the interior velocity (\ie the velocity of the first/last model level): 958 965 \[ 959 966 % \label{eq:zdfbfr_linear} 960 967 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 961 968 \] 962 where $r$ is a friction coefficient expressed in ms$^{-1}$.969 where $r$ is a friction coefficient expressed in $m s^{-1}$. 963 970 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 964 971 and setting $r = H / \tau$, where $H$ is the ocean depth. 965 Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{ Weatherly_JMR84}.972 Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}. 966 973 A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used in quasi-geostrophic models. 967 974 One may consider the linear friction as an approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ 968 (\citet{ Gill1982}, Eq. 9.6.6).975 (\citet{gill_bk82}, Eq. 9.6.6). 969 976 For example, with a drag coefficient $C_D = 0.002$, a typical speed of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$, 970 977 and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 971 978 This is the default value used in \NEMO. It corresponds to a decay time scale of 115~days. 972 It can be changed by specifying \np{rn\_ bfri1} (namelist parameter).973 974 For the linear friction case the coefficients defined in the general expression \autoref{eq:zdfbfr_bdef} are:979 It can be changed by specifying \np{rn\_Uc0} (namelist parameter). 980 981 For the linear friction case the drag coefficient used in the general expression \autoref{eq:zdfbfr_bdef} is: 975 982 \[ 976 983 % \label{eq:zdfbfr_linbfr_b} 977 \begin{split} 978 c_b^u &= - r\\ 979 c_b^v &= - r\\ 980 \end{split} 981 \] 982 When \np{nn\_botfr}\forcode{ = 1}, the value of $r$ used is \np{rn\_bfri1}. 983 Setting \np{nn\_botfr}\forcode{ = 0} is equivalent to setting $r=0$ and 984 leads to a free-slip bottom boundary condition. 985 These values are assigned in \mdl{zdfbfr}. 986 From v3.2 onwards there is support for local enhancement of these values via an externally defined 2D mask array 987 (\np{ln\_bfr2d}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file. 984 c_b^T = - r 985 \] 986 When \np{ln\_lin} \forcode{= .true.}, the value of $r$ used is \np{rn\_Uc0}*\np{rn\_Cd0}. 987 Setting \np{ln\_OFF} \forcode{= .true.} (and \forcode{ln_lin = .true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition. 988 989 These values are assigned in \mdl{zdfdrg}. 990 Note that there is support for local enhancement of these values via an externally defined 2D mask array 991 (\np{ln\_boost}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file. 988 992 The mask values should vary from 0 to 1. 989 993 Locations with a non-zero mask value will have the friction coefficient increased by 990 $mask\_value$ *\np{rn\_bfrien}*\np{rn\_bfri1}.994 $mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}. 991 995 992 996 % ------------------------------------------------------------------------------------------------------------- 993 997 % Non-Linear Bottom Friction 994 998 % ------------------------------------------------------------------------------------------------------------- 995 \subsection{Non-linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 2})} 996 \label{subsec:ZDF_bfr_nonlinear} 997 998 The non-linear bottom friction parameterisation assumes that the bottom friction is quadratic: 999 \[ 1000 % \label{eq:zdfbfr_nonlinear} 999 \subsection[Non-linear top/bottom friction (\forcode{ln_non_lin = .true.})] 1000 {Non-linear top/bottom friction (\protect\np{ln\_non\_lin}\forcode{ = .true.})} 1001 \label{subsec:ZDF_drg_nonlinear} 1002 1003 The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic: 1004 \[ 1005 % \label{eq:zdfdrg_nonlinear} 1001 1006 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 1002 1007 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b 1003 1008 \] 1004 where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy due to tides,1009 where $C_D$ is a drag coefficient, and $e_b $ a top/bottom turbulent kinetic energy due to tides, 1005 1010 internal waves breaking and other short time scale currents. 1006 1011 A typical value of the drag coefficient is $C_D = 10^{-3} $. 1007 As an example, the CME experiment \citep{ Treguier_JGR92} uses $C_D = 10^{-3}$ and1008 $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{ Killworth1992} uses $C_D = 1.4\;10^{-3}$ and1012 As an example, the CME experiment \citep{treguier_JGR92} uses $C_D = 10^{-3}$ and 1013 $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{killworth_JPO92} uses $C_D = 1.4\;10^{-3}$ and 1009 1014 $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. 1010 The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} namelist parameters). 1011 1012 As for the linear case, the bottom friction is imposed in the code by adding the trend due to 1013 the bottom friction to the general momentum trend in \mdl{dynbfr}. 1014 For the non-linear friction case the terms computed in \mdl{zdfbfr} are: 1015 \[ 1016 % \label{eq:zdfbfr_nonlinbfr} 1017 \begin{split} 1018 c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\ 1019 c_b^v &= - \; C_D\;\left[ \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ 1020 \end{split} 1021 \] 1022 1023 The coefficients that control the strength of the non-linear bottom friction are initialised as namelist parameters: 1024 $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. 1025 Note for applications which treat tides explicitly a low or even zero value of \np{rn\_bfeb2} is recommended. 1026 From v3.2 onwards a local enhancement of $C_D$ is possible via an externally defined 2D mask array 1027 (\np{ln\_bfr2d}\forcode{ = .true.}). 1028 This works in the same way as for the linear bottom friction case with non-zero masked locations increased by 1029 $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}. 1015 The CME choices have been set as default values (\np{rn\_Cd0} and \np{rn\_ke0} namelist parameters). 1016 1017 As for the linear case, the friction is imposed in the code by adding the trend due to 1018 the friction to the general momentum trend in \mdl{dynzdf}. 1019 For the non-linear friction case the term computed in \mdl{zdfdrg} is: 1020 \[ 1021 % \label{eq:zdfdrg_nonlinbfr} 1022 c_b^T = - \; C_D\;\left[ \left(\bar{u_b}^{i}\right)^2 + \left(\bar{v_b}^{j}\right)^2 + e_b \right]^{1/2} 1023 \] 1024 1025 The coefficients that control the strength of the non-linear friction are initialised as namelist parameters: 1026 $C_D$= \np{rn\_Cd0}, and $e_b$ =\np{rn\_bfeb2}. 1027 Note that for applications which consider tides explicitly, a low or even zero value of \np{rn\_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array 1028 (\np{ln\_boost}\forcode{ = .true.}). 1029 This works in the same way as for the linear friction case with non-zero masked locations increased by 1030 $mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}. 1030 1031 1031 1032 % ------------------------------------------------------------------------------------------------------------- 1032 1033 % Bottom Friction Log-layer 1033 1034 % ------------------------------------------------------------------------------------------------------------- 1034 \subsection[Log-layer btm frict enhncmnt (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})] 1035 {Log-layer bottom friction enhancement (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})} 1036 \label{subsec:ZDF_bfr_loglayer} 1037 1038 In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally enhanced using 1039 a "law of the wall" scaling. 1040 If \np{ln\_loglayer} = .true., $C_D$ is no longer constant but is related to the thickness of 1041 the last wet layer in each column by: 1042 \[ 1043 C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2 1044 \] 1045 1046 \noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness length provided via 1047 the namelist. 1048 1049 For stability, the drag coefficient is bounded such that it is kept greater or equal to 1050 the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional namelist parameter: 1051 \np{rn\_bfri2\_max}, \ie 1052 \[ 1053 rn\_bfri2 \leq C_D \leq rn\_bfri2\_max 1054 \] 1055 1056 \noindent Note also that a log-layer enhancement can also be applied to the top boundary friction if 1057 under ice-shelf cavities are in use (\np{ln\_isfcav}\forcode{ = .true.}). 1058 In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}. 1059 1060 % ------------------------------------------------------------------------------------------------------------- 1061 % Bottom Friction stability 1062 % ------------------------------------------------------------------------------------------------------------- 1063 \subsection{Bottom friction stability considerations} 1064 \label{subsec:ZDF_bfr_stability} 1065 1066 Some care needs to exercised over the choice of parameters to ensure that the implementation of 1067 bottom friction does not induce numerical instability. 1068 For the purposes of stability analysis, an approximation to \autoref{eq:zdfbfr_flux2} is: 1035 \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer = .true.})] 1036 {Log-layer top/bottom friction (\protect\np{ln\_loglayer}\forcode{ = .true.})} 1037 \label{subsec:ZDF_drg_loglayer} 1038 1039 In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using 1040 a "law of the wall" scaling. This assumes that the model vertical resolution can capture the logarithmic layer which typically occur for layers thinner than 1 m or so. 1041 If \np{ln\_loglayer} \forcode{= .true.}, $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness): 1042 \[ 1043 C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2 1044 \] 1045 1046 \noindent where $\kappa$ is the von-Karman constant and \np{rn\_z0} is a roughness length provided via the namelist. 1047 1048 The drag coefficient is bounded such that it is kept greater or equal to 1049 the base \np{rn\_Cd0} value which occurs where layer thicknesses become large and presumably logarithmic layers are not resolved at all. For stability reason, it is also not allowed to exceed the value of an additional namelist parameter: 1050 \np{rn\_Cdmax}, \ie 1051 \[ 1052 rn\_Cd0 \leq C_D \leq rn\_Cdmax 1053 \] 1054 1055 \noindent The log-layer enhancement can also be applied to the top boundary friction if 1056 under ice-shelf cavities are activated (\np{ln\_isfcav}\forcode{ = .true.}). 1057 %In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}. 1058 1059 % ------------------------------------------------------------------------------------------------------------- 1060 % Explicit bottom Friction 1061 % ------------------------------------------------------------------------------------------------------------- 1062 \subsection{Explicit top/bottom friction (\forcode{ln_drgimp = .false.})} 1063 \label{subsec:ZDF_drg_stability} 1064 1065 Setting \np{ln\_drgimp} \forcode{= .false.} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads: 1066 1067 At the top (below an ice shelf cavity): 1068 \[ 1069 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1070 = c_{t}^{\textbf{U}}\textbf{u}^{n-1}_{t} 1071 \] 1072 1073 At the bottom (above the sea floor): 1074 \[ 1075 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1076 = c_{b}^{\textbf{U}}\textbf{u}^{n-1}_{b} 1077 \] 1078 1079 Since this is conditionally stable, some care needs to exercised over the choice of parameters to ensure that the implementation of explicit top/bottom friction does not induce numerical instability. 1080 For the purposes of stability analysis, an approximation to \autoref{eq:zdfdrg_flux2} is: 1069 1081 \begin{equation} 1070 \label{eq:Eqn_ bfrstab}1082 \label{eq:Eqn_drgstab} 1071 1083 \begin{split} 1072 1084 \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt \\ … … 1074 1086 \end{split} 1075 1087 \end{equation} 1076 \noindent where linear bottomfriction and a leapfrog timestep have been assumed.1077 To ensure that the bottomfriction cannot reverse the direction of flow it is necessary to have:1088 \noindent where linear friction and a leapfrog timestep have been assumed. 1089 To ensure that the friction cannot reverse the direction of flow it is necessary to have: 1078 1090 \[ 1079 1091 |\Delta u| < \;|u| 1080 1092 \] 1081 \noindent which, using \autoref{eq:Eqn_ bfrstab}, gives:1093 \noindent which, using \autoref{eq:Eqn_drgstab}, gives: 1082 1094 \[ 1083 1095 r\frac{2\rdt}{e_{3u}} < 1 \qquad \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ … … 1093 1105 For most applications, with physically sensible parameters these restrictions should not be of concern. 1094 1106 But caution may be necessary if attempts are made to locally enhance the bottom friction parameters. 1095 To ensure stability limits are imposed on the bottom friction coefficients both1107 To ensure stability limits are imposed on the top/bottom friction coefficients both 1096 1108 during initialisation and at each time step. 1097 Checks at initialisation are made in \mdl{zdf bfr} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case).1109 Checks at initialisation are made in \mdl{zdfdrg} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case). 1098 1110 The number of breaches of the stability criterion are reported as well as 1099 1111 the minimum and maximum values that have been set. 1100 The criterion is also checked at each time step, using the actual velocity, in \mdl{dyn bfr}.1101 Values of the bottomfriction coefficient are reduced as necessary to ensure stability;1112 The criterion is also checked at each time step, using the actual velocity, in \mdl{dynzdf}. 1113 Values of the friction coefficient are reduced as necessary to ensure stability; 1102 1114 these changes are not reported. 1103 1115 1104 Limits on the bottom friction coefficient are not imposed if the user has elected to1105 handle the bottom friction implicitly (see \autoref{subsec:ZDF_bfr_imp}).1116 Limits on the top/bottom friction coefficient are not imposed if the user has elected to 1117 handle the friction implicitly (see \autoref{subsec:ZDF_drg_imp}). 1106 1118 The number of potential breaches of the explicit stability criterion are still reported for information purposes. 1107 1119 … … 1109 1121 % Implicit Bottom Friction 1110 1122 % ------------------------------------------------------------------------------------------------------------- 1111 \subsection{Implicit bottom friction (\protect\np{ln\_bfrimp}\forcode{ = .true.})} 1112 \label{subsec:ZDF_bfr_imp} 1123 \subsection[Implicit top/bottom friction (\forcode{ln_drgimp = .true.})] 1124 {Implicit top/bottom friction (\protect\np{ln\_drgimp}\forcode{ = .true.})} 1125 \label{subsec:ZDF_drg_imp} 1113 1126 1114 1127 An optional implicit form of bottom friction has been implemented to improve model stability. 1115 We recommend this option for shelf sea and coastal ocean applications, especially for split-explicit time splitting. 1116 This option can be invoked by setting \np{ln\_bfrimp} to \forcode{.true.} in the \textit{nambfr} namelist. 1117 This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \textit{namzdf} namelist. 1118 1119 This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, 1120 the bottom boundary condition is implemented implicitly. 1121 1122 \[ 1123 % \label{eq:dynzdf_bfr} 1124 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk} 1125 = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}} 1126 \] 1127 1128 where $mbk$ is the layer number of the bottom wet layer. 1129 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so, it is implicit. 1130 1131 If split-explicit time splitting is used, care must be taken to avoid the double counting of the bottom friction in 1132 the 2-D barotropic momentum equations. 1133 As NEMO only updates the barotropic pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, 1134 we need to remove the bottom friction induced by these two terms which has been included in the 3-D momentum trend 1135 and update it with the latest value. 1136 On the other hand, the bottom friction contributed by the other terms 1137 (\eg the advection term, viscosity term) has been included in the 3-D momentum equations and 1138 should not be added in the 2-D barotropic mode. 1139 1140 The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the following: 1141 1142 \[ 1143 % \label{eq:dynspg_ts_bfr1} 1144 \frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b} 1145 \left(\textbf{U}_{med}-\textbf{U}^{m-1}\right) 1146 \] 1147 \[ 1148 \frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+ 1149 \left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)- 1150 2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right) 1151 \] 1152 1153 where $\textbf{T}$ is the vertical integrated 3-D momentum trend. 1154 We assume the leap-frog time-stepping is used here. 1155 $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step. 1156 $c_{b}$ is the friction coefficient. 1157 $\eta$ is the sea surface level calculated in the barotropic loops while $\eta^{'}$ is the sea surface level used in 1158 the 3-D baroclinic mode. 1159 $\textbf{u}_{b}$ is the bottom layer horizontal velocity. 1160 1161 % ------------------------------------------------------------------------------------------------------------- 1162 % Bottom Friction with split-explicit time splitting 1163 % ------------------------------------------------------------------------------------------------------------- 1164 \subsection[Bottom friction w/ split-explicit time splitting (\protect\np{ln\_bfrimp})] 1165 {Bottom friction with split-explicit time splitting (\protect\np{ln\_bfrimp})} 1166 \label{subsec:ZDF_bfr_ts} 1167 1168 When calculating the momentum trend due to bottom friction in \mdl{dynbfr}, 1169 the bottom velocity at the before time step is used. 1170 This velocity includes both the baroclinic and barotropic components which is appropriate when 1171 using either the explicit or filtered surface pressure gradient algorithms 1172 (\key{dynspg\_exp} or \key{dynspg\_flt}). 1173 Extra attention is required, however, when using split-explicit time stepping (\key{dynspg\_ts}). 1174 In this case the free surface equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, 1175 while the three dimensional prognostic variables are solved with the longer time step of \np{rn\_rdt} seconds. 1176 The trend in the barotropic momentum due to bottom friction appropriate to this method is that given by 1177 the selected parameterisation (\ie linear or non-linear bottom friction) computed with 1178 the evolving velocities at each barotropic timestep. 1179 1180 In the case of non-linear bottom friction, we have elected to partially linearise the problem by 1181 keeping the coefficients fixed throughout the barotropic time-stepping to those computed in 1182 \mdl{zdfbfr} using the now timestep. 1183 This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to: 1184 1128 We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting. 1129 This option can be invoked by setting \np{ln\_drgimp} to \forcode{.true.} in the \textit{namdrg} namelist. 1130 %This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \textit{namzdf} namelist. 1131 1132 This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step: 1133 1134 At the top (below an ice shelf cavity): 1135 \[ 1136 % \label{eq:dynzdf_drg_top} 1137 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1138 = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t} 1139 \] 1140 1141 At the bottom (above the sea floor): 1142 \[ 1143 % \label{eq:dynzdf_drg_bot} 1144 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1145 = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b} 1146 \] 1147 1148 where $t$ and $b$ refers to top and bottom layers respectively. 1149 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit. 1150 1151 % ------------------------------------------------------------------------------------------------------------- 1152 % Bottom Friction with split-explicit free surface 1153 % ------------------------------------------------------------------------------------------------------------- 1154 \subsection[Bottom friction with split-explicit free surface] 1155 {Bottom friction with split-explicit free surface} 1156 \label{subsec:ZDF_drg_ts} 1157 1158 With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \np{ln\_drgimp}\forcode{= .false.} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie \forcode{ln_bt_fw = .false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \np{ln\_drgimp}\forcode{= .true.}, stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions. 1159 1160 The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO is as follows: 1185 1161 \begin{enumerate} 1186 \item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before barotropic velocity to 1187 the bottom friction component of the vertically integrated momentum trend. 1188 Note the same stability check that is carried out on the bottom friction coefficient in \mdl{dynbfr} has to 1189 be applied here to ensure that the trend removed matches that which was added in \mdl{dynbfr}. 1190 \item At each barotropic step, compute the contribution of the current barotropic velocity to 1191 the trend due to bottom friction. 1192 Add this contribution to the vertically integrated momentum trend. 1193 This contribution is handled implicitly which eliminates the need to impose a stability criteria on 1194 the values of the bottom friction coefficient within the barotropic loop. 1195 \end{enumerate} 1196 1197 Note that the use of an implicit formulation within the barotropic loop for the bottom friction trend means that 1198 any limiting of the bottom friction coefficient in \mdl{dynbfr} does not adversely affect the solution when 1199 using split-explicit time splitting. 1200 This is because the major contribution to bottom friction is likely to come from the barotropic component which 1201 uses the unrestricted value of the coefficient. 1202 However, if the limiting is thought to be having a major effect 1203 (a more likely prospect in coastal and shelf seas applications) then 1204 the fully implicit form of the bottom friction should be used (see \autoref{subsec:ZDF_bfr_imp}) 1205 which can be selected by setting \np{ln\_bfrimp} $=$ \forcode{.true.}. 1206 1207 Otherwise, the implicit formulation takes the form: 1208 \[ 1209 % \label{eq:zdfbfr_implicitts} 1210 \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ] 1211 \] 1212 where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height), 1213 $c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and 1214 $RHS$ represents all the components to the vertically integrated momentum trend except for 1215 that due to bottom friction. 1216 1217 % ================================================================ 1218 % Tidal Mixing 1219 % ================================================================ 1220 \section{Tidal mixing (\protect\key{zdftmx})} 1221 \label{sec:ZDF_tmx} 1222 1223 %--------------------------------------------namzdf_tmx-------------------------------------------------- 1162 \item To extend the stability of the barotropic sub-stepping, bottom stresses are refreshed at each sub-iteration. The baroclinic part of the flow entering the stresses is frozen at the initial time of the barotropic iteration. In case of non-linear friction, the drag coefficient is also constant. 1163 \item In case of an implicit drag, specific computations are performed in \mdl{dynzdf} which renders the overall scheme mixed explicit/implicit: the barotropic components of 3d velocities are removed before seeking for the implicit vertical diffusion result. Top/bottom stresses due to the barotropic components are explicitly accounted for thanks to the updated values of barotropic velocities. Then the implicit solution of 3d velocities is obtained. Lastly, the residual barotropic component is replaced by the time split estimate. 1164 \end{enumerate} 1165 1166 Note that other strategies are possible, like considering vertical diffusion step in advance, \ie prior barotropic integration. 1167 1168 1169 % ================================================================ 1170 % Internal wave-driven mixing 1171 % ================================================================ 1172 \section[Internal wave-driven mixing (\forcode{ln_zdfiwm = .true.})] 1173 {Internal wave-driven mixing (\protect\np{ln\_zdfiwm}\forcode{ = .true.})} 1174 \label{subsec:ZDF_tmx_new} 1175 1176 %--------------------------------------------namzdf_iwm------------------------------------------ 1224 1177 % 1225 %\nlst{namzdf_tmx}1178 \nlst{namzdf_iwm} 1226 1179 %-------------------------------------------------------------------------------------------------------------- 1227 1180 1228 1229 % -------------------------------------------------------------------------------------------------------------1230 % Bottom intensified tidal mixing1231 % -------------------------------------------------------------------------------------------------------------1232 \subsection{Bottom intensified tidal mixing}1233 \label{subsec:ZDF_tmx_bottom}1234 1235 Options are defined through the \ngn{namzdf\_tmx} namelist variables.1236 The parameterization of tidal mixing follows the general formulation for the vertical eddy diffusivity proposed by1237 \citet{St_Laurent_al_GRL02} and first introduced in an OGCM by \citep{Simmons_al_OM04}.1238 In this formulation an additional vertical diffusivity resulting from internal tide breaking,1239 $A^{vT}_{tides}$ is expressed as a function of $E(x,y)$,1240 the energy transfer from barotropic tides to baroclinic tides:1241 \begin{equation}1242 \label{eq:Ktides}1243 A^{vT}_{tides} = q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 }1244 \end{equation}1245 where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),1246 $\rho$ the density, $q$ the tidal dissipation efficiency, and $F(z)$ the vertical structure function.1247 1248 The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter) and1249 is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980).1250 The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter)1251 represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally,1252 with the remaining $1-q$ radiating away as low mode internal waves and1253 contributing to the background internal wave field.1254 A value of $q=1/3$ is typically used \citet{St_Laurent_al_GRL02}.1255 The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical.1256 It is implemented as a simple exponential decaying upward away from the bottom,1257 with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter,1258 with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04},1259 \[1260 % \label{eq:Fz}1261 F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }1262 \]1263 and is normalized so that vertical integral over the water column is unity.1264 1265 The associated vertical viscosity is calculated from the vertical diffusivity assuming a Prandtl number of 1,1266 \ie $A^{vm}_{tides}=A^{vT}_{tides}$.1267 In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity is capped at $300\,cm^2/s$ and1268 impose a lower limit on $N^2$ of \np{rn\_n2min} usually set to $10^{-8} s^{-2}$.1269 These bounds are usually rarely encountered.1270 1271 The internal wave energy map, $E(x, y)$ in \autoref{eq:Ktides}, is derived from a barotropic model of1272 the tides utilizing a parameterization of the conversion of barotropic tidal energy into internal waves.1273 The essential goal of the parameterization is to represent the momentum exchange between the barotropic tides and1274 the unrepresented internal waves induced by the tidal flow over rough topography in a stratified ocean.1275 In the current version of \NEMO, the map is built from the output of1276 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.1277 This model provides the dissipation associated with internal wave energy for the M2 and K1 tides component1278 (\autoref{fig:ZDF_M2_K1_tmx}).1279 The S2 dissipation is simply approximated as being $1/4$ of the M2 one.1280 The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$.1281 Its global mean value is $1.1$ TW,1282 in agreement with independent estimates \citep{Egbert_Ray_Nat00, Egbert_Ray_JGR01}.1283 1284 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>1285 \begin{figure}[!t]1286 \begin{center}1287 \includegraphics[width=0.90\textwidth]{Fig_ZDF_M2_K1_tmx}1288 \caption{1289 \protect\label{fig:ZDF_M2_K1_tmx}1290 (a) M2 and (b) K1 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$).1291 }1292 \end{center}1293 \end{figure}1294 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>1295 1296 % -------------------------------------------------------------------------------------------------------------1297 % Indonesian area specific treatment1298 % -------------------------------------------------------------------------------------------------------------1299 \subsection{Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})}1300 \label{subsec:ZDF_tmx_itf}1301 1302 When the Indonesian Through Flow (ITF) area is included in the model domain,1303 a specific treatment of tidal induced mixing in this area can be used.1304 It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide an input NetCDF file,1305 \ifile{mask\_itf}, which contains a mask array defining the ITF area where the specific treatment is applied.1306 1307 When \np{ln\_tmx\_itf}\forcode{ = .true.}, the two key parameters $q$ and $F(z)$ are adjusted following1308 the parameterisation developed by \citet{Koch-Larrouy_al_GRL07}:1309 1310 First, the Indonesian archipelago is a complex geographic region with a series of1311 large, deep, semi-enclosed basins connected via numerous narrow straits.1312 Once generated, internal tides remain confined within this semi-enclosed area and hardly radiate away.1313 Therefore all the internal tides energy is consumed within this area.1314 So it is assumed that $q = 1$, \ie all the energy generated is available for mixing.1315 Note that for test purposed, the ITF tidal dissipation efficiency is a namelist parameter (\np{rn\_tfe\_itf}).1316 A value of $1$ or close to is this recommended for this parameter.1317 1318 Second, the vertical structure function, $F(z)$, is no more associated with a bottom intensification of the mixing,1319 but with a maximum of energy available within the thermocline.1320 \citet{Koch-Larrouy_al_GRL07} have suggested that the vertical distribution of1321 the energy dissipation proportional to $N^2$ below the core of the thermocline and to $N$ above.1322 The resulting $F(z)$ is:1323 \[1324 % \label{eq:Fz_itf}1325 F(i,j,k) \sim \left\{1326 \begin{aligned}1327 \frac{q\,\Gamma E(i,j) } {\rho N \, \int N dz} \qquad \text{when $\partial_z N < 0$} \\1328 \frac{q\,\Gamma E(i,j) } {\rho \, \int N^2 dz} \qquad \text{when $\partial_z N > 0$}1329 \end{aligned}1330 \right.1331 \]1332 1333 Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$,1334 which agrees with the independent estimates inferred from observations.1335 Introduced in a regional OGCM, the parameterization improves the water mass characteristics in1336 the different Indonesian seas, suggesting that the horizontal and vertical distributions of1337 the mixing are adequately prescribed \citep{Koch-Larrouy_al_GRL07, Koch-Larrouy_al_OD08a, Koch-Larrouy_al_OD08b}.1338 Note also that such a parameterisation has a significant impact on the behaviour of1339 global coupled GCMs \citep{Koch-Larrouy_al_CD10}.1340 1341 % ================================================================1342 % Internal wave-driven mixing1343 % ================================================================1344 \section{Internal wave-driven mixing (\protect\key{zdftmx\_new})}1345 \label{sec:ZDF_tmx_new}1346 1347 %--------------------------------------------namzdf_tmx_new------------------------------------------1348 %1349 %\nlst{namzdf_tmx_new}1350 %--------------------------------------------------------------------------------------------------------------1351 1352 1181 The parameterization of mixing induced by breaking internal waves is a generalization of 1353 the approach originally proposed by \citet{ St_Laurent_al_GRL02}.1182 the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}. 1354 1183 A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed, 1355 1184 and the resulting diffusivity is obtained as … … 1360 1189 where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of 1361 1190 the energy available for mixing. 1362 If the \np{ln\_mevar} namelist parameter is set to false, the mixing efficiency is taken as constant and1363 equal to 1/6 \citep{ Osborn_JPO80}.1191 If the \np{ln\_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and 1192 equal to 1/6 \citep{osborn_JPO80}. 1364 1193 In the opposite (recommended) case, $R_f$ is instead a function of 1365 1194 the turbulence intensity parameter $Re_b = \frac{ \epsilon}{\nu \, N^2}$, 1366 with $\nu$ the molecular viscosity of seawater, following the model of \cite{ Bouffard_Boegman_DAO2013} and1367 the implementation of \cite{de _lavergne_JPO2016_efficiency}.1195 with $\nu$ the molecular viscosity of seawater, following the model of \cite{bouffard.boegman_DAO13} and 1196 the implementation of \cite{de-lavergne.madec.ea_JPO16}. 1368 1197 Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when 1369 1198 the mixing efficiency is constant. 1370 1199 1371 1200 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 1372 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice.1373 This parameterization of differential mixing, due to \cite{ Jackson_Rehmann_JPO2014},1374 is implemented as in \cite{de _lavergne_JPO2016_efficiency}.1201 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice. 1202 This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14}, 1203 is implemented as in \cite{de-lavergne.madec.ea_JPO16}. 1375 1204 1376 1205 The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, 1377 1206 is constructed from three static maps of column-integrated internal wave energy dissipation, 1378 $E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures 1379 (de Lavergne et al., in prep): 1207 $E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures: 1208 1380 1209 \begin{align*} 1381 1210 F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\ 1382 F_{pyc}(i,j,k) &\propto N^{n \_p}\\1211 F_{pyc}(i,j,k) &\propto N^{n_p}\\ 1383 1212 F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } 1384 1213 \end{align*} … … 1388 1217 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz' } \; , 1389 1218 \] 1390 The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_ tmx\_new} namelist)1219 The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_iwm} namelist) 1391 1220 controls the stratification-dependence of the pycnocline-intensified dissipation. 1392 It can take values of 1 (recommended) or 2.1221 It can take values of $1$ (recommended) or $2$. 1393 1222 Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of 1394 1223 the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps. 1395 1224 $h_{cri}$ is related to the large-scale topography of the ocean (etopo2) and 1396 1225 $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of 1397 the abyssal hill topography \citep{Goff_JGR2010} and the latitude. 1226 the abyssal hill topography \citep{goff_JGR10} and the latitude. 1227 % 1228 % Jc: input files names ? 1229 1230 % ================================================================ 1231 % surface wave-induced mixing 1232 % ================================================================ 1233 \section[Surface wave-induced mixing (\forcode{ln_zdfswm = .true.})] 1234 {Surface wave-induced mixing (\protect\np{ln\_zdfswm}\forcode{ = .true.})} 1235 \label{subsec:ZDF_swm} 1236 1237 Surface waves produce an enhanced mixing through wave-turbulence interaction. 1238 In addition to breaking waves induced turbulence (\autoref{subsec:ZDF_tke}), 1239 the influence of non-breaking waves can be accounted introducing 1240 wave-induced viscosity and diffusivity as a function of the wave number spectrum. 1241 Following \citet{qiao.yuan.ea_OD10}, a formulation of wave-induced mixing coefficient 1242 is provided as a function of wave amplitude, Stokes Drift and wave-number: 1243 1244 \begin{equation} 1245 \label{eq:Bv} 1246 B_{v} = \alpha {A} {U}_{st} {exp(3kz)} 1247 \end{equation} 1248 1249 Where $B_{v}$ is the wave-induced mixing coefficient, $A$ is the wave amplitude, 1250 ${U}_{st}$ is the Stokes Drift velocity, $k$ is the wave number and $\alpha$ 1251 is a constant which should be determined by observations or 1252 numerical experiments and is set to be 1. 1253 1254 The coefficient $B_{v}$ is then directly added to the vertical viscosity 1255 and diffusivity coefficients. 1256 1257 In order to account for this contribution set: \forcode{ln_zdfswm = .true.}, 1258 then wave interaction has to be activated through \forcode{ln_wave = .true.}, 1259 the Stokes Drift can be evaluated by setting \forcode{ln_sdw = .true.} 1260 (see \autoref{subsec:SBC_wave_sdw}) 1261 and the needed wave fields can be provided either in forcing or coupled mode 1262 (for more information on wave parameters and settings see \autoref{sec:SBC_wave}) 1263 1264 % ================================================================ 1265 % Adaptive-implicit vertical advection 1266 % ================================================================ 1267 \section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp = .true.})] 1268 {Adaptive-implicit vertical advection(\protect\np{ln\_zad\_Aimp}\forcode{ = .true.})} 1269 \label{subsec:ZDF_aimp} 1270 1271 The adaptive-implicit vertical advection option in NEMO is based on the work of 1272 \citep{shchepetkin_OM15}. In common with most ocean models, the timestep used with NEMO 1273 needs to satisfy multiple criteria associated with different physical processes in order 1274 to maintain numerical stability. \citep{shchepetkin_OM15} pointed out that the vertical 1275 CFL criterion is commonly the most limiting. \citep{lemarie.debreu.ea_OM15} examined the 1276 constraints for a range of time and space discretizations and provide the CFL stability 1277 criteria for a range of advection schemes. The values for the Leap-Frog with Robert 1278 asselin filter time-stepping (as used in NEMO) are reproduced in 1279 \autoref{tab:zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these 1280 restrictions but at the cost of large dispersive errors and, possibly, large numerical 1281 viscosity. The adaptive-implicit vertical advection option provides a targetted use of the 1282 implicit scheme only when and where potential breaches of the vertical CFL condition 1283 occur. In many practical applications these events may occur remote from the main area of 1284 interest or due to short-lived conditions such that the extra numerical diffusion or 1285 viscosity does not greatly affect the overall solution. With such applications, setting: 1286 \forcode{ln_zad_Aimp = .true.} should allow much longer model timesteps to be used whilst 1287 retaining the accuracy of the high order explicit schemes over most of the domain. 1288 1289 \begin{table}[htbp] 1290 \begin{center} 1291 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}} 1292 \begin{tabular}{r|ccc} 1293 \hline 1294 spatial discretization & 2nd order centered & 3rd order upwind & 4th order compact \\ 1295 advective CFL criterion & 0.904 & 0.472 & 0.522 \\ 1296 \hline 1297 \end{tabular} 1298 \caption{ 1299 \protect\label{tab:zad_Aimp_CFLcrit} 1300 The advective CFL criteria for a range of spatial discretizations for the Leap-Frog with Robert Asselin filter time-stepping 1301 ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}. 1302 } 1303 \end{center} 1304 \end{table} 1305 1306 In particular, the advection scheme remains explicit everywhere except where and when 1307 local vertical velocities exceed a threshold set just below the explicit stability limit. 1308 Once the threshold is reached a tapered transition towards an implicit scheme is used by 1309 partitioning the vertical velocity into a part that can be treated explicitly and any 1310 excess that must be treated implicitly. The partitioning is achieved via a Courant-number 1311 dependent weighting algorithm as described in \citep{shchepetkin_OM15}. 1312 1313 The local cell Courant number ($Cu$) used for this partitioning is: 1314 1315 \begin{equation} 1316 \label{eq:Eqn_zad_Aimp_Courant} 1317 \begin{split} 1318 Cu &= {2 \rdt \over e^n_{3t_{ijk}}} \bigg (\big [ \texttt{Max}(w^n_{ijk},0.0) - \texttt{Min}(w^n_{ijk+1},0.0) \big ] \\ 1319 &\phantom{=} +\big [ \texttt{Max}(e_{{2_u}ij}e^n_{{3_{u}}ijk}u^n_{ijk},0.0) - \texttt{Min}(e_{{2_u}i-1j}e^n_{{3_{u}}i-1jk}u^n_{i-1jk},0.0) \big ] 1320 \big / e_{{1_t}ij}e_{{2_t}ij} \\ 1321 &\phantom{=} +\big [ \texttt{Max}(e_{{1_v}ij}e^n_{{3_{v}}ijk}v^n_{ijk},0.0) - \texttt{Min}(e_{{1_v}ij-1}e^n_{{3_{v}}ij-1k}v^n_{ij-1k},0.0) \big ] 1322 \big / e_{{1_t}ij}e_{{2_t}ij} \bigg ) \\ 1323 \end{split} 1324 \end{equation} 1325 1326 \noindent and the tapering algorithm follows \citep{shchepetkin_OM15} as: 1327 1328 \begin{align} 1329 \label{eq:Eqn_zad_Aimp_partition} 1330 Cu_{min} &= 0.15 \nonumber \\ 1331 Cu_{max} &= 0.3 \nonumber \\ 1332 Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ 1333 Fcu &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ 1334 C\kern-0.14em f &= 1335 \begin{cases} 1336 0.0 &\text{if $Cu \leq Cu_{min}$} \\ 1337 (Cu - Cu_{min})^2 / (Fcu + (Cu - Cu_{min})^2) &\text{else if $Cu < Cu_{cut}$} \\ 1338 (Cu - Cu_{max}) / Cu &\text{else} 1339 \end{cases} 1340 \end{align} 1341 1342 \begin{figure}[!t] 1343 \begin{center} 1344 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_coeff} 1345 \caption{ 1346 \protect\label{fig:zad_Aimp_coeff} 1347 The value of the partitioning coefficient ($C\kern-0.14em f$) used to partition vertical velocities into parts to 1348 be treated implicitly and explicitly for a range of typical Courant numbers (\forcode{ln_zad_Aimp=.true.}) 1349 } 1350 \end{center} 1351 \end{figure} 1352 1353 \noindent The partitioning coefficient is used to determine the part of the vertical 1354 velocity that must be handled implicitly ($w_i$) and to subtract this from the total 1355 vertical velocity ($w_n$) to leave that which can continue to be handled explicitly: 1356 1357 \begin{align} 1358 \label{eq:Eqn_zad_Aimp_partition2} 1359 w_{i_{ijk}} &= C\kern-0.14em f_{ijk} w_{n_{ijk}} \nonumber \\ 1360 w_{n_{ijk}} &= (1-C\kern-0.14em f_{ijk}) w_{n_{ijk}} 1361 \end{align} 1362 1363 \noindent Note that the coefficient is such that the treatment is never fully implicit; 1364 the three cases from \autoref{eq:Eqn_zad_Aimp_partition} can be considered as: 1365 fully-explicit; mixed explicit/implicit and mostly-implicit. With the settings shown the 1366 coefficient ($C\kern-0.14em f$) varies as shown in \autoref{fig:zad_Aimp_coeff}. Note with these values 1367 the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly 1368 implicit' is 0.45 which is just below the stability limited given in 1369 \autoref{tab:zad_Aimp_CFLcrit} for a 3rd order scheme. 1370 1371 The $w_i$ component is added to the implicit solvers for the vertical mixing in 1372 \mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}. This is 1373 sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further 1374 intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}). 1375 For these schemes the implicit upstream fluxes must be added to both the monotonic guess 1376 and to the higher order solution when calculating the antidiffusive fluxes. The implicit 1377 vertical fluxes are then removed since they are added by the implicit solver later on. 1378 1379 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 1380 used in a wide range of simulations. The following test simulation, however, does illustrate 1381 the potential benefits and will hopefully encourage further testing and feedback from users: 1382 1383 \begin{figure}[!t] 1384 \begin{center} 1385 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} 1386 \caption{ 1387 \protect\label{fig:zad_Aimp_overflow_frames} 1388 A time-series of temperature vertical cross-sections for the OVERFLOW test case. These results are for the default 1389 settings with \forcode{nn_rdt=10.0} and without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}). 1390 } 1391 \end{center} 1392 \end{figure} 1393 1394 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} 1395 The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case} 1396 provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case 1397 by only a few extra physics choices namely: 1398 1399 \begin{verbatim} 1400 ln_dynldf_OFF = .false. 1401 ln_dynldf_lap = .true. 1402 ln_dynldf_hor = .true. 1403 ln_zdfnpc = .true. 1404 ln_traadv_fct = .true. 1405 nn_fct_h = 2 1406 nn_fct_v = 2 1407 \end{verbatim} 1408 1409 \noindent which were chosen to provide a slightly more stable and less noisy solution. The 1410 result when using the default value of \forcode{nn_rdt = 10.} without adaptive-implicit 1411 vertical velocity is illustrated in \autoref{fig:zad_Aimp_overflow_frames}. The mass of 1412 cold water, initially sitting on the shelf, moves down the slope and forms a 1413 bottom-trapped, dense plume. Even with these extra physics choices the model is close to 1414 stability limits and attempts with \forcode{nn_rdt = 30.} will fail after about 5.5 hours 1415 with excessively high horizontal velocities. This time-scale corresponds with the time the 1416 plume reaches the steepest part of the topography and, although detected as a horizontal 1417 CFL breach, the instability originates from a breach of the vertical CFL limit. This is a good 1418 candidate, therefore, for use of the adaptive-implicit vertical advection scheme. 1419 1420 The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps 1421 are shown in \autoref{fig:zad_Aimp_overflow_all_rdt} (together with the equivalent 1422 frames from the base run). In this simple example the use of the adaptive-implicit 1423 vertcal advection scheme has enabled a 12x increase in the model timestep without 1424 significantly altering the solution (although at this extreme the plume is more diffuse 1425 and has not travelled so far). Notably, the solution with and without the scheme is 1426 slightly different even with \forcode{nn_rdt = 10.}; suggesting that the base run was 1427 close enough to instability to trigger the scheme despite completing successfully. 1428 To assist in diagnosing how active the scheme is, in both location and time, the 3D 1429 implicit and explicit components of the vertical velocity are available via XIOS as 1430 \texttt{wimp} and \texttt{wexp} respectively. Likewise, the partitioning coefficient 1431 ($C\kern-0.14em f$) is also available as \texttt{wi\_cff}. For a quick oversight of 1432 the schemes activity the global maximum values of the absolute implicit component 1433 of the vertical velocity and the partitioning coefficient are written to the netCDF 1434 version of the run statistics file (\texttt{run.stat.nc}) if this is active (see 1435 \autoref{sec:MISC_opt} for activation details). 1436 1437 \autoref{fig:zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for 1438 the various overflow tests. Note that the adaptive-implicit vertical advection scheme is 1439 active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the 1440 test case is close to stability limits even with this value. At the larger timesteps, the 1441 vertical velocity is treated mostly implicitly at some location throughout the run. The 1442 oscillatory nature of this measure appears to be linked to the progress of the plume front 1443 as each cusp is associated with the location of the maximum shifting to the adjacent cell. 1444 This is illustrated in \autoref{fig:zad_Aimp_maxCf_loc} where the i- and k- locations of the 1445 maximum have been overlaid for the base run case. 1446 1447 \medskip 1448 \noindent Only limited tests have been performed in more realistic configurations. In the 1449 ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes 1450 restartability and reproducibility tests but it is unable to improve the model's stability 1451 enough to allow an increase in the model time-step. A view of the time-series of maximum 1452 partitioning coefficient (not shown here) suggests that the default time-step of 5400s is 1453 already pushing at stability limits, especially in the initial start-up phase. The 1454 time-series does not, however, exhibit any of the 'cuspiness' found with the overflow 1455 tests. 1456 1457 \medskip 1458 \noindent A short test with an eORCA1 configuration promises more since a test using a 1459 time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the 1460 time-step is limited to 2700s without. 1461 1462 \begin{figure}[!t] 1463 \begin{center} 1464 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_all_rdt} 1465 \caption{ 1466 \protect\label{fig:zad_Aimp_overflow_all_rdt} 1467 Sample temperature vertical cross-sections from mid- and end-run using different values for \forcode{nn_rdt} 1468 and with or without adaptive implicit vertical advection. Without the adaptive implicit vertical advection only 1469 the run with the shortest timestep is able to run to completion. Note also that the colour-scale has been 1470 chosen to confirm that temperatures remain within the original range of 10$^o$ to 20$^o$. 1471 } 1472 \end{center} 1473 \end{figure} 1474 1475 \begin{figure}[!t] 1476 \begin{center} 1477 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf} 1478 \caption{ 1479 \protect\label{fig:zad_Aimp_maxCf} 1480 The maximum partitioning coefficient during a series of test runs with increasing model timestep length. 1481 At the larger timesteps, the vertical velocity is treated mostly implicitly at some location throughout 1482 the run. 1483 } 1484 \end{center} 1485 \end{figure} 1486 1487 \begin{figure}[!t] 1488 \begin{center} 1489 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf_loc} 1490 \caption{ 1491 \protect\label{fig:zad_Aimp_maxCf_loc} 1492 The maximum partitioning coefficient for the \forcode{nn_rdt=10.0s} case overlaid with information on the gridcell i- and k- 1493 locations of the maximum value. 1494 } 1495 \end{center} 1496 \end{figure} 1398 1497 1399 1498 % ================================================================
Note: See TracChangeset
for help on using the changeset viewer.