Changeset 11573 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/doc/latex/NEMO/subfiles/chap_ZDF.tex
- Timestamp:
- 2019-09-19T11:18:03+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/doc/latex/NEMO/subfiles/chap_ZDF.tex
r11225 r11573 1 1 \documentclass[../main/NEMO_manual]{subfiles} 2 3 %% Custom aliases 4 \newcommand{\cf}{\ensuremath{C\kern-0.14em f}} 2 5 3 6 \begin{document} … … 8 11 \label{chap:ZDF} 9 12 10 \ minitoc13 \chaptertoc 11 14 12 15 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. … … 18 21 % ================================================================ 19 22 \section{Vertical mixing} 20 \label{sec:ZDF _zdf}23 \label{sec:ZDF} 21 24 22 25 The discrete form of the ocean subgrid scale physics has been presented in … … 25 28 At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}), 26 29 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 \np{ln\_trabbc} defined,30 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie\ \np{ln\_trabbc} defined, 28 31 see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum 29 (see \autoref{sec:ZDF_drg}). 32 (see \autoref{sec:ZDF_drg}). 30 33 31 34 In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and … … 37 40 the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules. 38 41 The trends due to the vertical momentum and tracer diffusion, including the surface forcing, 39 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 42 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 40 43 %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 scheme42 %(\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 %(namelist parameter \np{ln\_zdfexp}\forcode{=.true.}) or a backward time stepping scheme 45 %(\np{ln\_zdfexp}\forcode{=.false.}) depending on the magnitude of the mixing coefficients, 46 %and thus of the formulation used (see \autoref{chap:TD}). 44 47 45 48 %--------------------------------------------namzdf-------------------------------------------------------- 46 49 47 \nlst{namzdf} 50 \begin{listing} 51 \nlst{namzdf} 52 \caption{\forcode{&namzdf}} 53 \label{lst:namzdf} 54 \end{listing} 48 55 %-------------------------------------------------------------------------------------------------------------- 49 56 50 57 % ------------------------------------------------------------------------------------------------------------- 51 % Constant 52 % ------------------------------------------------------------------------------------------------------------- 53 \subsection[Constant (\forcode{ln_zdfcst = .true.})] 54 {Constant (\protect\np{ln\_zdfcst}\forcode{ = .true.})} 58 % Constant 59 % ------------------------------------------------------------------------------------------------------------- 60 \subsection[Constant (\forcode{ln_zdfcst})]{Constant (\protect\np{ln\_zdfcst})} 55 61 \label{subsec:ZDF_cst} 56 62 57 Options are defined through the \n gn{namzdf} namelist variables.63 Options are defined through the \nam{zdf} namelist variables. 58 64 When \np{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 59 65 constant values over the whole ocean. … … 66 72 \end{align*} 67 73 68 These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. 74 These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. 69 75 In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity, 70 76 that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and … … 74 80 % Richardson Number Dependent 75 81 % ------------------------------------------------------------------------------------------------------------- 76 \subsection[Richardson number dependent (\forcode{ln_zdfric = .true.})] 77 {Richardson number dependent (\protect\np{ln\_zdfric}\forcode{ = .true.})} 82 \subsection[Richardson number dependent (\forcode{ln_zdfric})]{Richardson number dependent (\protect\np{ln\_zdfric})} 78 83 \label{subsec:ZDF_ric} 79 84 80 85 %--------------------------------------------namric--------------------------------------------------------- 81 86 82 \nlst{namzdf_ric} 87 \begin{listing} 88 \nlst{namzdf_ric} 89 \caption{\forcode{&namzdf_ric}} 90 \label{lst:namzdf_ric} 91 \end{listing} 83 92 %-------------------------------------------------------------------------------------------------------------- 84 93 85 When \np{ln\_zdfric}\forcode{ =.true.}, a local Richardson number dependent formulation for the vertical momentum and86 tracer eddy coefficients is set through the \n gn{namzdf\_ric} namelist variables.87 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. 94 When \np{ln\_zdfric}\forcode{=.true.}, a local Richardson number dependent formulation for the vertical momentum and 95 tracer eddy coefficients is set through the \nam{zdf\_ric} namelist variables. 96 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. 88 97 \textit{In situ} measurements have been used to link vertical turbulent activity to large scale ocean structures. 89 98 The hypothesis of a mixing mainly maintained by the growth of Kelvin-Helmholtz like instabilities leads to 90 99 a dependency between the vertical eddy coefficients and the local Richardson number 91 (\ie the ratio of stratification to vertical shear).100 (\ie\ the ratio of stratification to vertical shear). 92 101 Following \citet{pacanowski.philander_JPO81}, the following formulation has been implemented: 93 102 \[ 94 % \label{eq: zdfric}103 % \label{eq:ZDF_ric} 95 104 \left\{ 96 105 \begin{aligned} … … 101 110 \] 102 111 where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number, 103 $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 112 $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 104 113 $A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the constant case 105 114 (see \autoref{subsec:ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ is the maximum value that … … 109 118 110 119 A simple mixing-layer model to transfer and dissipate the atmospheric forcings 111 (wind-stress and buoyancy fluxes) can be activated setting the \np{ln\_mldw}\forcode{ =.true.} in the namelist.120 (wind-stress and buoyancy fluxes) can be activated setting the \np{ln\_mldw}\forcode{=.true.} in the namelist. 112 121 113 122 In this case, the local depth of turbulent wind-mixing or "Ekman depth" $h_{e}(x,y,t)$ is evaluated and … … 130 139 131 140 % ------------------------------------------------------------------------------------------------------------- 132 % TKE Turbulent Closure Scheme 133 % ------------------------------------------------------------------------------------------------------------- 134 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke = .true.})] 135 {TKE turbulent closure scheme (\protect\np{ln\_zdftke}\forcode{ = .true.})} 141 % TKE Turbulent Closure Scheme 142 % ------------------------------------------------------------------------------------------------------------- 143 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke})]{TKE turbulent closure scheme (\protect\np{ln\_zdftke})} 136 144 \label{subsec:ZDF_tke} 137 145 %--------------------------------------------namzdf_tke-------------------------------------------------- 138 146 139 \nlst{namzdf_tke} 147 \begin{listing} 148 \nlst{namzdf_tke} 149 \caption{\forcode{&namzdf_tke}} 150 \label{lst:namzdf_tke} 151 \end{listing} 140 152 %-------------------------------------------------------------------------------------------------------------- 141 153 … … 144 156 and a closure assumption for the turbulent length scales. 145 157 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,158 adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of \NEMO, 147 159 by \citet{blanke.delecluse_JPO93} for equatorial Atlantic simulations. 148 160 Since then, significant modifications have been introduced by \citet{madec.delecluse.ea_NPM98} in both the implementation and … … 151 163 its destruction through stratification, its vertical diffusion, and its dissipation of \citet{kolmogorov_IANS42} type: 152 164 \begin{equation} 153 \label{eq: zdftke_e}165 \label{eq:ZDF_tke_e} 154 166 \frac{\partial \bar{e}}{\partial t} = 155 167 \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2 … … 161 173 \end{equation} 162 174 \[ 163 % \label{eq: zdftke_kz}175 % \label{eq:ZDF_tke_kz} 164 176 \begin{split} 165 177 K_m &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ … … 167 179 \end{split} 168 180 \] 169 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 170 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 181 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 182 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 171 183 $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients. 172 184 The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with 173 vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}. 185 vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}. 174 186 They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}. 175 187 $P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$: 176 188 \begin{align*} 177 % \label{eq: prt}189 % \label{eq:ZDF_prt} 178 190 P_{rt} = 179 191 \begin{cases} … … 199 211 too weak vertical diffusion. 200 212 They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and 201 \np{rn\_avt0} (\n gn{namzdf} namelist, see \autoref{subsec:ZDF_cst}).213 \np{rn\_avt0} (\nam{zdf} namelist, see \autoref{subsec:ZDF_cst}). 202 214 203 215 \subsubsection{Turbulent length scale} … … 208 220 The first two are based on the following first order approximation \citep{blanke.delecluse_JPO93}: 209 221 \begin{equation} 210 \label{eq: tke_mxl0_1}222 \label{eq:ZDF_tke_mxl0_1} 211 223 l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N 212 224 \end{equation} 213 225 which is valid in a stable stratified region with constant values of the Brunt-Vais\"{a}l\"{a} frequency. 214 226 The resulting length scale is bounded by the distance to the surface or to the bottom 215 (\np{nn\_mxl}\forcode{ = 0}) or by the local vertical scale factor (\np{nn\_mxl}\forcode{ =1}).227 (\np{nn\_mxl}\forcode{=0}) or by the local vertical scale factor (\np{nn\_mxl}\forcode{=1}). 216 228 \citet{blanke.delecluse_JPO93} notice that this simplification has two major drawbacks: 217 229 it makes no sense for locally unstable stratification and the computation no longer uses all 218 230 the information contained in the vertical density profile. 219 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{ =2, 3} cases,231 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{=2, 3} cases, 220 232 which add an extra assumption concerning the vertical gradient of the computed length scale. 221 So, the length scales are first evaluated as in \autoref{eq: tke_mxl0_1} and then bounded such that:233 So, the length scales are first evaluated as in \autoref{eq:ZDF_tke_mxl0_1} and then bounded such that: 222 234 \begin{equation} 223 \label{eq: tke_mxl_constraint}235 \label{eq:ZDF_tke_mxl_constraint} 224 236 \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1 225 237 \qquad \text{with }\ l = l_k = l_\epsilon 226 238 \end{equation} 227 \autoref{eq: tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than239 \autoref{eq:ZDF_tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than 228 240 the variations of depth. 229 It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less 241 It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less 230 242 time consuming. 231 243 In particular, it allows the length scale to be limited not only by the distance to the surface or 232 244 to the ocean bottom but also by the distance to a strongly stratified portion of the water column such as 233 the thermocline (\autoref{fig: mixing_length}).234 In order to impose the \autoref{eq: tke_mxl_constraint} constraint, we introduce two additional length scales:245 the thermocline (\autoref{fig:ZDF_mixing_length}). 246 In order to impose the \autoref{eq:ZDF_tke_mxl_constraint} constraint, we introduce two additional length scales: 235 247 $l_{up}$ and $l_{dwn}$, the upward and downward length scales, and 236 248 evaluate the dissipation and mixing length scales as … … 238 250 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 239 251 \begin{figure}[!t] 240 \begin{center} 241 \includegraphics[width=\textwidth]{Fig_mixing_length} 242 \caption{ 243 \protect\label{fig:mixing_length} 244 Illustration of the mixing length computation. 245 } 246 \end{center} 252 \centering 253 \includegraphics[width=0.66\textwidth]{Fig_mixing_length} 254 \caption[Mixing length computation]{Illustration of the mixing length computation} 255 \label{fig:ZDF_mixing_length} 247 256 \end{figure} 248 257 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 249 258 \[ 250 % \label{eq: tke_mxl2}259 % \label{eq:ZDF_tke_mxl2} 251 260 \begin{aligned} 252 261 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) … … 256 265 \end{aligned} 257 266 \] 258 where $l^{(k)}$ is computed using \autoref{eq: tke_mxl0_1}, \ie$l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.259 260 In the \np{nn\_mxl}\forcode{ =2} case, the dissipation and mixing length scales take the same value:261 $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np{nn\_mxl}\forcode{ =3} case,267 where $l^{(k)}$ is computed using \autoref{eq:ZDF_tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. 268 269 In the \np{nn\_mxl}\forcode{=2} case, the dissipation and mixing length scales take the same value: 270 $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np{nn\_mxl}\forcode{=3} case, 262 271 the dissipation and mixing turbulent length scales are give as in \citet{gaspar.gregoris.ea_JGR90}: 263 272 \[ 264 % \label{eq: tke_mxl_gaspar}273 % \label{eq:ZDF_tke_mxl_gaspar} 265 274 \begin{aligned} 266 275 & l_k = \sqrt{\ l_{up} \ \ l_{dwn}\ } \\ … … 283 292 This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow. 284 293 The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and 285 air-sea drag coefficient. 286 The latter concerns the bulk formulae and is not discussed here. 294 air-sea drag coefficient. 295 The latter concerns the bulk formulae and is not discussed here. 287 296 288 297 Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is : … … 292 301 \end{equation} 293 302 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}. 303 ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}. 295 304 The boundary condition on the turbulent length scale follows the Charnock's relation: 296 305 \begin{equation} … … 303 312 $\alpha_{CB} = 100$ the Craig and Banner's value. 304 313 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 305 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds 314 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds 306 315 to $\alpha_{CB} = 100$. 307 Further setting \np{ln\_mxl0 =.true.}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,316 Further setting \np{ln\_mxl0}\forcode{ =.true.}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale, 308 317 with $\beta$ hard coded to the Stacey's value. 309 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the 318 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the 310 319 surface $\bar{e}$ value. 311 320 … … 319 328 The detailed physics behind LC is described in, for example, \citet{craik.leibovich_JFM76}. 320 329 The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and 321 wind drift currents. 330 wind drift currents. 322 331 323 332 Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by … … 325 334 The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in 326 335 an extra source term of TKE, $P_{LC}$. 327 The presence of $P_{LC}$ in \autoref{eq: zdftke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to328 \forcode{.true.} in the \n gn{namzdf\_tke} namelist.329 336 The presence of $P_{LC}$ in \autoref{eq:ZDF_tke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to 337 \forcode{.true.} in the \nam{zdf\_tke} namelist. 338 330 339 By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}), 331 $P_{LC}$ is assumed to be : 340 $P_{LC}$ is assumed to be : 332 341 \[ 333 342 P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}} 334 343 \] 335 344 where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth. 336 With no information about the wave field, $w_{LC}$ is assumed to be proportional to 337 the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module 345 With no information about the wave field, $w_{LC}$ is assumed to be proportional to 346 the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module 338 347 \footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as 339 348 $u_s = 0.016 \,|U_{10m}|$. … … 343 352 For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at 344 353 a finite depth $H_{LC}$ (which is often close to the mixed layer depth), 345 and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures). 354 and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures). 346 355 The resulting expression for $w_{LC}$ is : 347 356 \[ … … 355 364 The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second. 356 365 The value of $c_{LC}$ is set through the \np{rn\_lc} namelist parameter, 357 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}. 366 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}. 358 367 359 368 The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations: 360 369 $H_{LC}$ is the depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by 361 converting its kinetic energy to potential energy, according to 370 converting its kinetic energy to potential energy, according to 362 371 \[ 363 372 - \int_{-H_{LC}}^0 { N^2\;z \;dz} = \frac{1}{2} u_s^2 … … 370 379 produce mixed-layer depths that are too shallow during summer months and windy conditions. 371 380 This bias is particularly acute over the Southern Ocean. 372 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}. 373 The parameterization is an empirical one, \ie not derived from theoretical considerations,374 but rather is meant to account for observed processes that affect the density structure of 375 the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme 376 (\ie near-inertial oscillations and ocean swells and waves).377 378 When using this parameterization (\ie when \np{nn\_etau}\forcode{ =1}),381 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}. 382 The parameterization is an empirical one, \ie\ not derived from theoretical considerations, 383 but rather is meant to account for observed processes that affect the density structure of 384 the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme 385 (\ie\ near-inertial oscillations and ocean swells and waves). 386 387 When using this parameterization (\ie\ when \np{nn\_etau}\forcode{=1}), 379 388 the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations, 380 389 swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition, … … 387 396 penetrates in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of 388 397 the penetration, and $f_i$ is the ice concentration 389 (no penetration if $f_i=1$, \ie if the ocean is entirely covered by sea-ice).398 (no penetration if $f_i=1$, \ie\ if the ocean is entirely covered by sea-ice). 390 399 The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. 391 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}\forcode{ =0}) or400 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}\forcode{=0}) or 392 401 a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m at high latitudes 393 (\np{nn\_etau}\forcode{ = 1}).394 395 Note that two other option exist, \np{nn\_etau}\forcode{ =2, 3}.402 (\np{nn\_etau}\forcode{=1}). 403 404 Note that two other option exist, \np{nn\_etau}\forcode{=2, 3}. 396 405 They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer, 397 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean. 406 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean. 398 407 Those two options are obsolescent features introduced for test purposes. 399 They will be removed in the next release. 408 They will be removed in the next release. 400 409 401 410 % This should be explain better below what this rn_eice parameter is meant for: 402 411 In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn\_eice} namelist parameter. 403 412 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\%. 405 406 % from Burchard et al OM 2008 : 407 % the most critical process not reproduced by statistical turbulence models is the activity of 408 % internal waves and their interaction with turbulence. After the Reynolds decomposition, 409 % internal waves are in principle included in the RANS equations, but later partially 410 % excluded by the hydrostatic assumption and the model resolution. 411 % Thus far, the representation of internal wave mixing in ocean models has been relatively crude 412 % (\eg Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 413 414 % ------------------------------------------------------------------------------------------------------------- 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.})} 413 is greater than 25\%. 414 415 % from Burchard et al OM 2008 : 416 % the most critical process not reproduced by statistical turbulence models is the activity of 417 % internal waves and their interaction with turbulence. After the Reynolds decomposition, 418 % internal waves are in principle included in the RANS equations, but later partially 419 % excluded by the hydrostatic assumption and the model resolution. 420 % Thus far, the representation of internal wave mixing in ocean models has been relatively crude 421 % (\eg\ Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 422 423 % ------------------------------------------------------------------------------------------------------------- 424 % GLS Generic Length Scale Scheme 425 % ------------------------------------------------------------------------------------------------------------- 426 \subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls})]{GLS: Generic Length Scale (\protect\np{ln\_zdfgls})} 419 427 \label{subsec:ZDF_gls} 420 428 421 429 %--------------------------------------------namzdf_gls--------------------------------------------------------- 422 430 423 \nlst{namzdf_gls} 431 \begin{listing} 432 \nlst{namzdf_gls} 433 \caption{\forcode{&namzdf_gls}} 434 \label{lst:namzdf_gls} 435 \end{listing} 424 436 %-------------------------------------------------------------------------------------------------------------- 425 437 … … 427 439 one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale, 428 440 $\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 of441 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 442 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:ZDF_GLS} allows to recover a number of 431 443 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}). 444 $k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}). 433 445 The GLS scheme is given by the following set of equations: 434 446 \begin{equation} 435 \label{eq: zdfgls_e}447 \label{eq:ZDF_gls_e} 436 448 \frac{\partial \bar{e}}{\partial t} = 437 449 \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 … … 443 455 444 456 \[ 445 % \label{eq: zdfgls_psi}457 % \label{eq:ZDF_gls_psi} 446 458 \begin{split} 447 459 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ … … 455 467 456 468 \[ 457 % \label{eq: zdfgls_kz}469 % \label{eq:ZDF_gls_kz} 458 470 \begin{split} 459 471 K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ … … 463 475 464 476 \[ 465 % \label{eq: zdfgls_eps}477 % \label{eq:ZDF_gls_eps} 466 478 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; 467 479 \] 468 480 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and 469 $\epsilon$ the dissipation rate. 481 $\epsilon$ the dissipation rate. 470 482 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of 471 483 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. 484 Four different turbulent models are pre-defined (\autoref{tab:ZDF_GLS}). 485 They are made available through the \np{nn\_clo} namelist parameter. 474 486 475 487 %--------------------------------------------------TABLE-------------------------------------------------- 476 488 \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} 489 \centering 490 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} 491 \begin{tabular}{ccccc} 492 & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ 493 % & \citep{mellor.yamada_RG82} & \citep{rodi_JGR87} & \citep{wilcox_AJ88} & \\ 494 \hline 495 \hline 496 \np{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ 497 \hline 498 $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ 499 $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ 500 $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ 501 $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ 502 $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ 503 $C_3$ & 1. & 1. & 1. & 1. \\ 504 $F_{wall}$ & Yes & -- & -- & -- \\ 505 \hline 506 \hline 507 \end{tabular} 508 \caption[Set of predefined GLS parameters or equivalently predefined turbulence models available]{ 509 Set of predefined GLS parameters, or equivalently predefined turbulence models available with 510 \protect\np{ln\_zdfgls}\forcode{=.true.} and controlled by 511 the \protect\np{nn\_clos} namelist variable in \protect\nam{zdf\_gls}.} 512 \label{tab:ZDF_GLS} 502 513 \end{table} 503 514 %-------------------------------------------------------------------------------------------------------------- … … 508 519 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88}, 509 520 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.).521 (\np{nn\_stab\_func}\forcode{=0, 3}, resp.). 511 522 The value of $C_{0\mu}$ depends on the choice of the stability function. 512 523 … … 516 527 \np{rn\_crban}\forcode{ > 0.} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}. 517 528 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}. 529 \np{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}. 519 530 520 531 The $\psi$ equation is known to fail in stably stratified flows, and for this reason … … 525 536 the entrainment depth predicted in stably stratified situations, 526 537 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.},538 The clipping is only activated if \np{ln\_length\_lim}\forcode{=.true.}, 528 539 and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value. 529 540 … … 531 542 the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}. 532 543 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.})} 544 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO\ model. 545 546 547 % ------------------------------------------------------------------------------------------------------------- 548 % OSM OSMOSIS BL Scheme 549 % ------------------------------------------------------------------------------------------------------------- 550 \subsection[OSM: OSMosis boundary layer scheme (\forcode{ln_zdfosm})]{OSM: OSMosis boundary layer scheme (\protect\np{ln\_zdfosm})} 541 551 \label{subsec:ZDF_osm} 542 552 %--------------------------------------------namzdf_osm--------------------------------------------------------- 543 553 544 \nlst{namzdf_osm} 554 \begin{listing} 555 \nlst{namzdf_osm} 556 \caption{\forcode{&namzdf_osm}} 557 \label{lst:namzdf_osm} 558 \end{listing} 545 559 %-------------------------------------------------------------------------------------------------------------- 546 560 … … 550 564 % TKE and GLS discretization considerations 551 565 % ------------------------------------------------------------------------------------------------------------- 552 \subsection[ Discrete energy conservation for TKE and GLS schemes] 553 {Discrete energy conservation for TKE and GLS schemes} 566 \subsection[ Discrete energy conservation for TKE and GLS schemes]{Discrete energy conservation for TKE and GLS schemes} 554 567 \label{subsec:ZDF_tke_ene} 555 568 556 569 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 557 570 \begin{figure}[!t] 558 \begin{center} 559 \includegraphics[width=\textwidth]{Fig_ZDF_TKE_time_scheme} 560 \caption{ 561 \protect\label{fig:TKE_time_scheme} 562 Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and its links to the momentum and tracer time integration. 563 } 564 \end{center} 571 \centering 572 \includegraphics[width=0.66\textwidth]{Fig_ZDF_TKE_time_scheme} 573 \caption[Subgrid kinetic energy integration in GLS and TKE schemes]{ 574 Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and 575 its links to the momentum and tracer time integration.} 576 \label{fig:ZDF_TKE_time_scheme} 565 577 \end{figure} 566 578 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 567 579 568 580 The production of turbulence by vertical shear (the first term of the right hand side of 569 \autoref{eq: zdftke_e}) and \autoref{eq:zdfgls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion570 (first line in \autoref{eq: PE_zdf}).581 \autoref{eq:ZDF_tke_e}) and \autoref{eq:ZDF_gls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion 582 (first line in \autoref{eq:MB_zdf}). 571 583 To do so a special care has to be taken for both the time and space discretization of 572 584 the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}. 573 585 574 Let us first address the time stepping issue. \autoref{fig: TKE_time_scheme} shows how586 Let us first address the time stepping issue. \autoref{fig:ZDF_TKE_time_scheme} shows how 575 587 the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with 576 588 the one-level forward time stepping of the equation for $\bar{e}$. 577 589 With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to 578 590 the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and 579 summing the result vertically: 591 summing the result vertically: 580 592 \begin{equation} 581 \label{eq: energ1}593 \label{eq:ZDF_energ1} 582 594 \begin{split} 583 595 \int_{-H}^{\eta} u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt} \right) \,dz \\ … … 587 599 \end{equation} 588 600 Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$, 589 known at time $t$ (\autoref{fig: TKE_time_scheme}), as it is required when using the TKE scheme590 (see \autoref{sec: STP_forward_imp}).591 The first term of the right hand side of \autoref{eq: energ1} represents the kinetic energy transfer at601 known at time $t$ (\autoref{fig:ZDF_TKE_time_scheme}), as it is required when using the TKE scheme 602 (see \autoref{sec:TD_forward_imp}). 603 The first term of the right hand side of \autoref{eq:ZDF_energ1} represents the kinetic energy transfer at 592 604 the surface (atmospheric forcing) and at the bottom (friction effect). 593 605 The second term is always negative. 594 606 It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$. 595 \autoref{eq: energ1} implies that, to be energetically consistent,607 \autoref{eq:ZDF_energ1} implies that, to be energetically consistent, 596 608 the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 597 609 ${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ … … 599 611 600 612 A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 601 (second term of the right hand side of \autoref{eq: zdftke_e} and \autoref{eq:zdfgls_e}).613 (second term of the right hand side of \autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}). 602 614 This term must balance the input of potential energy resulting from vertical mixing. 603 615 The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by 604 616 multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically: 605 617 \begin{equation} 606 \label{eq: energ2}618 \label{eq:ZDF_energ2} 607 619 \begin{split} 608 620 \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt} \right) \,dz \\ … … 613 625 \end{split} 614 626 \end{equation} 615 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 616 The first term of the right hand side of \autoref{eq: energ2} is always zero because627 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 628 The first term of the right hand side of \autoref{eq:ZDF_energ2} is always zero because 617 629 there is no diffusive flux through the ocean surface and bottom). 618 630 The second term is minus the destruction rate of $\bar{e}$ due to stratification. 619 Therefore \autoref{eq: energ1} implies that, to be energetically consistent,620 the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq: zdftke_e} and \autoref{eq:zdfgls_e}.631 Therefore \autoref{eq:ZDF_energ1} implies that, to be energetically consistent, 632 the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}. 621 633 622 634 Let us now address the space discretization issue. 623 635 The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in 624 the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig: cell}).636 the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:DOM_cell}). 625 637 A space averaging is thus required to obtain the shear TKE production term. 626 By redoing the \autoref{eq: energ1} in the 3D case, it can be shown that the product of eddy coefficient by638 By redoing the \autoref{eq:ZDF_energ1} in the 3D case, it can be shown that the product of eddy coefficient by 627 639 the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. 628 640 Furthermore, the time variation of $e_3$ has be taken into account. … … 630 642 The above energetic considerations leads to the following final discrete form for the TKE equation: 631 643 \begin{equation} 632 \label{eq: zdftke_ene}644 \label{eq:ZDF_tke_ene} 633 645 \begin{split} 634 646 \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt} \equiv … … 647 659 \end{split} 648 660 \end{equation} 649 where the last two terms in \autoref{eq: zdftke_ene} (vertical diffusion and Kolmogorov dissipation)650 are time stepped using a backward scheme (see\autoref{sec: STP_forward_imp}).661 where the last two terms in \autoref{eq:ZDF_tke_ene} (vertical diffusion and Kolmogorov dissipation) 662 are time stepped using a backward scheme (see\autoref{sec:TD_forward_imp}). 651 663 Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible. 652 664 %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. 665 %they all appear in the right hand side of \autoref{eq:ZDF_tke_ene}. 666 %For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 655 667 656 668 % ================================================================ … … 660 672 \label{sec:ZDF_conv} 661 673 662 Static instabilities (\ie light potential densities under heavy ones) may occur at particular ocean grid points.674 Static instabilities (\ie\ light potential densities under heavy ones) may occur at particular ocean grid points. 663 675 In nature, convective processes quickly re-establish the static stability of the water column. 664 676 These processes have been removed from the model via the hydrostatic assumption so they must be parameterized. … … 668 680 669 681 % ------------------------------------------------------------------------------------------------------------- 670 % Non-Penetrative Convective Adjustment 671 % ------------------------------------------------------------------------------------------------------------- 672 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc = .true.})] 673 {Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})} 682 % Non-Penetrative Convective Adjustment 683 % ------------------------------------------------------------------------------------------------------------- 684 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc})]{Non-penetrative convective adjustment (\protect\np{ln\_tranpc})} 674 685 \label{subsec:ZDF_npc} 675 686 676 687 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 677 688 \begin{figure}[!htb] 678 \begin{center} 679 \includegraphics[width=\textwidth]{Fig_npc} 680 \caption{ 681 \protect\label{fig:npc} 682 Example of an unstable density profile treated by the non penetrative convective adjustment algorithm. 683 $1^{st}$ step: the initial profile is checked from the surface to the bottom. 684 It is found to be unstable between levels 3 and 4. 685 They are mixed. 686 The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed. 687 The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed. 688 The $1^{st}$ step ends since the density profile is then stable below the level 3. 689 $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step: 690 levels 2 to 5 are mixed. 691 The new density profile is checked. 692 It is found stable: end of algorithm. 693 } 694 \end{center} 689 \centering 690 \includegraphics[width=0.66\textwidth]{Fig_npc} 691 \caption[Unstable density profile treated by the non penetrative convective adjustment algorithm]{ 692 Example of an unstable density profile treated by 693 the non penetrative convective adjustment algorithm. 694 $1^{st}$ step: the initial profile is checked from the surface to the bottom. 695 It is found to be unstable between levels 3 and 4. 696 They are mixed. 697 The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed. 698 The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed. 699 The $1^{st}$ step ends since the density profile is then stable below the level 3. 700 $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step: 701 levels 2 to 5 are mixed. 702 The new density profile is checked. 703 It is found stable: end of algorithm.} 704 \label{fig:ZDF_npc} 695 705 \end{figure} 696 706 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 697 707 698 Options are defined through the \n gn{namzdf} namelist variables.699 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}\forcode{ =.true.}.708 Options are defined through the \nam{zdf} namelist variables. 709 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}\forcode{=.true.}. 700 710 It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of 701 711 the water column, but only until the density structure becomes neutrally stable 702 (\ie until the mixed portion of the water column has \textit{exactly} the density of the water just below)712 (\ie\ until the mixed portion of the water column has \textit{exactly} the density of the water just below) 703 713 \citep{madec.delecluse.ea_JPO91}. 704 The associated algorithm is an iterative process used in the following way (\autoref{fig: npc}):714 The associated algorithm is an iterative process used in the following way (\autoref{fig:ZDF_npc}): 705 715 starting from the top of the ocean, the first instability is found. 706 716 Assume in the following that the instability is located between levels $k$ and $k+1$. … … 717 727 This algorithm is significantly different from mixing statically unstable levels two by two. 718 728 The latter procedure cannot converge with a finite number of iterations for some vertical profiles while 719 the algorithm used in \NEMO converges for any profile in a number of iterations which is less than729 the algorithm used in \NEMO\ converges for any profile in a number of iterations which is less than 720 730 the number of vertical levels. 721 731 This property is of paramount importance as pointed out by \citet{killworth_iprc89}: … … 727 737 (L. Brodeau, personnal communication). 728 738 Two main differences have been introduced compared to the original algorithm: 729 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 730 (not the difference in potential density); 739 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 740 (not the difference in potential density); 731 741 $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in 732 742 the same way their temperature and salinity has been mixed. … … 735 745 736 746 % ------------------------------------------------------------------------------------------------------------- 737 % Enhanced Vertical Diffusion 738 % ------------------------------------------------------------------------------------------------------------- 739 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd = .true.})] 740 {Enhanced vertical diffusion (\protect\np{ln\_zdfevd}\forcode{ = .true.})} 747 % Enhanced Vertical Diffusion 748 % ------------------------------------------------------------------------------------------------------------- 749 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd})]{Enhanced vertical diffusion (\protect\np{ln\_zdfevd})} 741 750 \label{subsec:ZDF_evd} 742 751 743 Options are defined through the \n gn{namzdf} namelist variables.744 The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}\forcode{ =.true.}.745 In this case, the vertical eddy mixing coefficients are assigned very large values 752 Options are defined through the \nam{zdf} namelist variables. 753 The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}\forcode{=.true.}. 754 In this case, the vertical eddy mixing coefficients are assigned very large values 746 755 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}.748 This is done either on tracers only (\np{nn\_evdm}\forcode{ =0}) or749 on both momentum and tracers (\np{nn\_evdm}\forcode{ =1}).750 751 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np{nn\_evdm}\forcode{ =1},756 (\ie\ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}. 757 This is done either on tracers only (\np{nn\_evdm}\forcode{=0}) or 758 on both momentum and tracers (\np{nn\_evdm}\forcode{=1}). 759 760 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np{nn\_evdm}\forcode{=1}, 752 761 the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ values also, are set equal to 753 762 the namelist parameter \np{rn\_avevd}. … … 759 768 Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$. 760 769 This removes a potential source of divergence of odd and even time step in 761 a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:STP_mLF}). 762 763 % ------------------------------------------------------------------------------------------------------------- 764 % Turbulent Closure Scheme 765 % ------------------------------------------------------------------------------------------------------------- 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.})} 770 a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:TD_mLF}). 771 772 % ------------------------------------------------------------------------------------------------------------- 773 % Turbulent Closure Scheme 774 % ------------------------------------------------------------------------------------------------------------- 775 \subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf_}\{\forcode{tke,gls,osm}\})]{Handling convection with turbulent closure schemes (\forcode{ln_zdf{tke,gls,osm}})} 768 776 \label{subsec:ZDF_tcs} 769 777 770 778 771 779 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,780 \autoref{subsec:ZDF_osm} (\ie\ \np{ln\_zdftke} or \np{ln\_zdfgls} or \np{ln\_zdfosm} defined) deal, in theory, 773 781 with statically unstable density profiles. 774 782 In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in 775 \autoref{eq: zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative.776 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also of the four neighboring values at 783 \autoref{eq:ZDF_tke_e} or \autoref{eq:ZDF_gls_e} becomes a source term, since $N^2$ is negative. 784 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also of the four neighboring values at 777 785 velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$). 778 786 These large values restore the static stability of the water column in a way similar to that of … … 782 790 because the mixing length scale is bounded by the distance to the sea surface. 783 791 It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme, 784 \ie setting the \np{ln\_zdfnpc} namelist parameter to true and792 \ie\ setting the \np{ln\_zdfnpc} namelist parameter to true and 785 793 defining the turbulent closure (\np{ln\_zdftke} or \np{ln\_zdfgls} = \forcode{.true.}) all together. 786 794 787 795 The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection, 788 796 %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.797 therefore \np{ln\_zdfevd}\forcode{=.false.} should be used with the OSMOSIS scheme. 790 798 % gm% + one word on non local flux with KPP scheme trakpp.F90 module... 791 799 … … 793 801 % Double Diffusion Mixing 794 802 % ================================================================ 795 \section[Double diffusion mixing (\forcode{ln_zdfddm = .true.})] 796 {Double diffusion mixing (\protect\np{ln\_zdfddm}\forcode{ = .true.})} 803 \section[Double diffusion mixing (\forcode{ln_zdfddm})]{Double diffusion mixing (\protect\np{ln\_zdfddm})} 797 804 \label{subsec:ZDF_ddm} 798 805 … … 804 811 805 812 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter 806 \np{ln\_zdfddm} in \n gn{namzdf}.813 \np{ln\_zdfddm} in \nam{zdf}. 807 814 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. 808 815 The former condition leads to salt fingering and the latter to diffusive convection. 809 816 Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the ocean. 810 \citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that 817 \citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that 811 818 it leads to relatively minor changes in circulation but exerts significant regional influences on 812 819 temperature and salinity. 813 820 814 821 815 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 822 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 816 823 \begin{align*} 817 % \label{eq: zdfddm_Kz}824 % \label{eq:ZDF_ddm_Kz} 818 825 &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\ 819 826 &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS} … … 827 834 (1981): 828 835 \begin{align} 829 \label{eq: zdfddm_f}836 \label{eq:ZDF_ddm_f} 830 837 A_f^{vS} &= 831 838 \begin{cases} … … 833 840 0 &\text{otherwise} 834 841 \end{cases} 835 \\ \label{eq: zdfddm_f_T}836 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 842 \\ \label{eq:ZDF_ddm_f_T} 843 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 837 844 \end{align} 838 845 839 846 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 840 847 \begin{figure}[!t] 841 \begin{center} 842 \includegraphics[width=\textwidth]{Fig_zdfddm} 843 \caption{ 844 \protect\label{fig:zdfddm} 845 From \citet{merryfield.holloway.ea_JPO99} : 846 (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in regions of salt fingering. 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}$; 848 (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in regions of 849 diffusive convection. 850 Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation. 851 The latter is not implemented in \NEMO. 852 } 853 \end{center} 848 \centering 849 \includegraphics[width=0.66\textwidth]{Fig_zdfddm} 850 \caption[Diapycnal diffusivities for temperature and salt in regions of salt fingering and 851 diffusive convection]{ 852 From \citet{merryfield.holloway.ea_JPO99}: 853 (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in 854 regions of salt fingering. 855 Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and 856 thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$; 857 (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in 858 regions of diffusive convection. 859 Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation. 860 The latter is not implemented in \NEMO.} 861 \label{fig:ZDF_ddm} 854 862 \end{figure} 855 863 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 856 864 857 The factor 0.7 in \autoref{eq: zdfddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx 0.7$ of865 The factor 0.7 in \autoref{eq:ZDF_ddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx 0.7$ of 858 866 buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{mcdougall.taylor_JMR84}). 859 867 Following \citet{merryfield.holloway.ea_JPO99}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 860 868 861 869 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by 862 Federov (1988) is used: 870 Federov (1988) is used: 863 871 \begin{align} 864 % \label{eq: zdfddm_d}872 % \label{eq:ZDF_ddm_d} 865 873 A_d^{vT} &= 866 874 \begin{cases} … … 870 878 \end{cases} 871 879 \nonumber \\ 872 \label{eq: zdfddm_d_S}880 \label{eq:ZDF_ddm_d_S} 873 881 A_d^{vS} &= 874 882 \begin{cases} … … 879 887 \end{align} 880 888 881 The dependencies of \autoref{eq: zdfddm_f} to \autoref{eq:zdfddm_d_S} on $R_\rho$ are illustrated in882 \autoref{fig: zdfddm}.889 The dependencies of \autoref{eq:ZDF_ddm_f} to \autoref{eq:ZDF_ddm_d_S} on $R_\rho$ are illustrated in 890 \autoref{fig:ZDF_ddm}. 883 891 Implementing this requires computing $R_\rho$ at each grid point on every time step. 884 892 This is done in \mdl{eosbn2} at the same time as $N^2$ is computed. … … 888 896 % Bottom Friction 889 897 % ================================================================ 890 \section[Bottom and top friction (\textit{zdfdrg.F90})] 891 {Bottom and top friction (\protect\mdl{zdfdrg})} 892 \label{sec:ZDF_drg} 893 894 %--------------------------------------------nambfr-------------------------------------------------------- 898 \section[Bottom and top friction (\textit{zdfdrg.F90})]{Bottom and top friction (\protect\mdl{zdfdrg})} 899 \label{sec:ZDF_drg} 900 901 %--------------------------------------------namdrg-------------------------------------------------------- 895 902 % 896 \nlst{namdrg} 897 \nlst{namdrg_top} 898 \nlst{namdrg_bot} 903 \begin{listing} 904 \nlst{namdrg} 905 \caption{\forcode{&namdrg}} 906 \label{lst:namdrg} 907 \end{listing} 908 \begin{listing} 909 \nlst{namdrg_top} 910 \caption{\forcode{&namdrg_top}} 911 \label{lst:namdrg_top} 912 \end{listing} 913 \begin{listing} 914 \nlst{namdrg_bot} 915 \caption{\forcode{&namdrg_bot}} 916 \label{lst:namdrg_bot} 917 \end{listing} 899 918 900 919 %-------------------------------------------------------------------------------------------------------------- 901 920 902 Options to define the top and bottom friction are defined through the \n gn{namdrg} namelist variables.921 Options to define the top and bottom friction are defined through the \nam{drg} namelist variables. 903 922 The bottom friction represents the friction generated by the bathymetry. 904 923 The top friction represents the friction generated by the ice shelf/ocean interface. 905 As the friction processes at the top and the bottom are treated in and identical way, 924 As the friction processes at the top and the bottom are treated in and identical way, 906 925 the description below considers mostly the bottom friction case, if not stated otherwise. 907 926 … … 911 930 For the bottom boundary layer, one has: 912 931 \[ 913 % \label{eq: zdfbfr_flux}932 % \label{eq:ZDF_bfr_flux} 914 933 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 915 934 \] … … 921 940 one needs a vertical diffusion coefficient $A^{vm} = 0.125$~m$^2$s$^{-1}$ 922 941 (for a Coriolis frequency $f = 10^{-4}$~m$^2$s$^{-1}$). 923 With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. 942 With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. 924 943 When the vertical mixing coefficient is this small, using a flux condition is equivalent to 925 944 entering the viscous forces (either wind stress or bottom friction) as a body force over the depth of the top or … … 927 946 To illustrate this, consider the equation for $u$ at $k$, the last ocean level: 928 947 \begin{equation} 929 \label{eq: zdfdrg_flux2}948 \label{eq:ZDF_drg_flux2} 930 949 \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}} 931 950 \end{equation} … … 947 966 These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where: 948 967 \begin{equation} 949 \label{eq: zdfbfr_bdef}968 \label{eq:ZDF_bfr_bdef} 950 969 \frac{\partial {\textbf U_h}}{\partial t} = 951 970 - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b 952 971 \end{equation} 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 (\ieat 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.972 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 973 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. 955 974 956 975 % ------------------------------------------------------------------------------------------------------------- 957 976 % Linear Bottom Friction 958 977 % ------------------------------------------------------------------------------------------------------------- 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} 978 \subsection[Linear top/bottom friction (\forcode{ln_lin})]{Linear top/bottom friction (\protect\np{ln\_lin})} 979 \label{subsec:ZDF_drg_linear} 962 980 963 981 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):965 \[ 966 % \label{eq: zdfbfr_linear}982 the friction is proportional to the interior velocity (\ie\ the velocity of the first/last model level): 983 \[ 984 % \label{eq:ZDF_bfr_linear} 967 985 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 968 986 \] 969 987 where $r$ is a friction coefficient expressed in $m s^{-1}$. 970 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 988 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 971 989 and setting $r = H / \tau$, where $H$ is the ocean depth. 972 990 Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}. … … 979 997 It can be changed by specifying \np{rn\_Uc0} (namelist parameter). 980 998 981 For the linear friction case the drag coefficient used in the general expression \autoref{eq: zdfbfr_bdef} is:982 \[ 983 % \label{eq: zdfbfr_linbfr_b}999 For the linear friction case the drag coefficient used in the general expression \autoref{eq:ZDF_bfr_bdef} is: 1000 \[ 1001 % \label{eq:ZDF_bfr_linbfr_b} 984 1002 c_b^T = - r 985 1003 \] 986 1004 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.1005 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 1006 989 1007 These values are assigned in \mdl{zdfdrg}. 990 1008 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.1009 (\np{ln\_boost}\forcode{=.true.}) given in the \ifile{bfr\_coef} input NetCDF file. 992 1010 The mask values should vary from 0 to 1. 993 1011 Locations with a non-zero mask value will have the friction coefficient increased by … … 997 1015 % Non-Linear Bottom Friction 998 1016 % ------------------------------------------------------------------------------------------------------------- 999 \subsection[Non-linear top/bottom friction (\forcode{ln_no_lin = .true.})] 1000 {Non-linear top/bottom friction (\protect\np{ln\_no\_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} 1017 \subsection[Non-linear top/bottom friction (\forcode{ln_non_lin})]{Non-linear top/bottom friction (\protect\np{ln\_non\_lin})} 1018 \label{subsec:ZDF_drg_nonlinear} 1019 1020 The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic: 1021 \[ 1022 % \label{eq:ZDF_drg_nonlinear} 1006 1023 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 1007 1024 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b … … 1019 1036 For the non-linear friction case the term computed in \mdl{zdfdrg} is: 1020 1037 \[ 1021 % \label{eq: zdfdrg_nonlinbfr}1038 % \label{eq:ZDF_drg_nonlinbfr} 1022 1039 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 1040 \] … … 1026 1043 $C_D$= \np{rn\_Cd0}, and $e_b$ =\np{rn\_bfeb2}. 1027 1044 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.}).1045 (\np{ln\_boost}\forcode{=.true.}). 1029 1046 This works in the same way as for the linear friction case with non-zero masked locations increased by 1030 1047 $mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}. … … 1033 1050 % Bottom Friction Log-layer 1034 1051 % ------------------------------------------------------------------------------------------------------------- 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} 1052 \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer})]{Log-layer top/bottom friction (\protect\np{ln\_loglayer})} 1053 \label{subsec:ZDF_drg_loglayer} 1038 1054 1039 1055 In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using … … 1054 1070 1055 1071 \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.}).1072 under ice-shelf cavities are activated (\np{ln\_isfcav}\forcode{=.true.}). 1057 1073 %In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}. 1058 1074 … … 1060 1076 % Explicit bottom Friction 1061 1077 % ------------------------------------------------------------------------------------------------------------- 1062 \subsection{Explicit top/bottom friction (\forcode{ln_drgimp =.false.})}1063 1078 \subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np{ln\_drgimp}\forcode{=.false.})} 1079 \label{subsec:ZDF_drg_stability} 1064 1080 1065 1081 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: … … 1078 1094 1079 1095 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:1096 For the purposes of stability analysis, an approximation to \autoref{eq:ZDF_drg_flux2} is: 1081 1097 \begin{equation} 1082 \label{eq: Eqn_drgstab}1098 \label{eq:ZDF_Eqn_drgstab} 1083 1099 \begin{split} 1084 1100 \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt \\ … … 1089 1105 To ensure that the friction cannot reverse the direction of flow it is necessary to have: 1090 1106 \[ 1091 |\Delta u| < \;|u| 1092 \] 1093 \noindent which, using \autoref{eq: Eqn_drgstab}, gives:1107 |\Delta u| < \;|u| 1108 \] 1109 \noindent which, using \autoref{eq:ZDF_Eqn_drgstab}, gives: 1094 1110 \[ 1095 1111 r\frac{2\rdt}{e_{3u}} < 1 \qquad \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ … … 1104 1120 For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then $e_{3u}$ should be greater than 3.6 m. 1105 1121 For most applications, with physically sensible parameters these restrictions should not be of concern. 1106 But caution may be necessary if attempts are made to locally enhance the bottom friction parameters. 1122 But caution may be necessary if attempts are made to locally enhance the bottom friction parameters. 1107 1123 To ensure stability limits are imposed on the top/bottom friction coefficients both 1108 1124 during initialisation and at each time step. … … 1121 1137 % Implicit Bottom Friction 1122 1138 % ------------------------------------------------------------------------------------------------------------- 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} 1139 \subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np{ln\_drgimp}\forcode{=.true.})} 1140 \label{subsec:ZDF_drg_imp} 1126 1141 1127 1142 An optional implicit form of bottom friction has been implemented to improve model stability. 1128 1143 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: 1144 This option can be invoked by setting \np{ln\_drgimp} to \forcode{.true.} in the \nam{drg} namelist. 1145 %This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \nam{zdf} namelist. 1146 1147 This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step: 1133 1148 1134 1149 At the top (below an ice shelf cavity): 1135 1150 \[ 1136 % \label{eq: dynzdf_drg_top}1151 % \label{eq:ZDF_dynZDF__drg_top} 1137 1152 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1138 1153 = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t} … … 1141 1156 At the bottom (above the sea floor): 1142 1157 \[ 1143 % \label{eq: dynzdf_drg_bot}1158 % \label{eq:ZDF_dynZDF__drg_bot} 1144 1159 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1145 1160 = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b} 1146 1161 \] 1147 1162 1148 where $t$ and $b$ refers to top and bottom layers respectively. 1163 where $t$ and $b$ refers to top and bottom layers respectively. 1149 1164 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit. 1150 1165 … … 1152 1167 % Bottom Friction with split-explicit free surface 1153 1168 % ------------------------------------------------------------------------------------------------------------- 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: 1169 \subsection[Bottom friction with split-explicit free surface]{Bottom friction with split-explicit free surface} 1170 \label{subsec:ZDF_drg_ts} 1171 1172 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. 1173 1174 The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO\ is as follows: 1161 1175 \begin{enumerate} 1162 1176 \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 1177 \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.1178 \end{enumerate} 1179 1180 Note that other strategies are possible, like considering vertical diffusion step in advance, \ie\ prior barotropic integration. 1167 1181 1168 1182 … … 1170 1184 % Internal wave-driven mixing 1171 1185 % ================================================================ 1172 \section[Internal wave-driven mixing (\forcode{ln_zdfiwm = .true.})] 1173 {Internal wave-driven mixing (\protect\np{ln\_zdfiwm}\forcode{ = .true.})} 1186 \section[Internal wave-driven mixing (\forcode{ln_zdfiwm})]{Internal wave-driven mixing (\protect\np{ln\_zdfiwm})} 1174 1187 \label{subsec:ZDF_tmx_new} 1175 1188 1176 1189 %--------------------------------------------namzdf_iwm------------------------------------------ 1177 1190 % 1178 \nlst{namzdf_iwm} 1191 \begin{listing} 1192 \nlst{namzdf_iwm} 1193 \caption{\forcode{&namzdf_iwm}} 1194 \label{lst:namzdf_iwm} 1195 \end{listing} 1179 1196 %-------------------------------------------------------------------------------------------------------------- 1180 1197 … … 1182 1199 the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}. 1183 1200 A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed, 1184 and the resulting diffusivity is obtained as 1185 \[ 1186 % \label{eq: Kwave}1201 and the resulting diffusivity is obtained as 1202 \[ 1203 % \label{eq:ZDF_Kwave} 1187 1204 A^{vT}_{wave} = R_f \,\frac{ \epsilon }{ \rho \, N^2 } 1188 1205 \] … … 1198 1215 the mixing efficiency is constant. 1199 1216 1200 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 1201 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice. 1217 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 1218 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice. 1202 1219 This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14}, 1203 1220 is implemented as in \cite{de-lavergne.madec.ea_JPO16}. … … 1211 1228 F_{pyc}(i,j,k) &\propto N^{n_p}\\ 1212 1229 F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } 1213 \end{align*} 1230 \end{align*} 1214 1231 In the above formula, $h_{ab}$ denotes the height above bottom, 1215 1232 $h_{wkb}$ denotes the WKB-stretched height above bottom, defined by … … 1217 1234 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz' } \; , 1218 1235 \] 1219 The $n_p$ parameter (given by \np{nn\_zpyc} in \n gn{namzdf\_iwm} namelist)1236 The $n_p$ parameter (given by \np{nn\_zpyc} in \nam{zdf\_iwm} namelist) 1220 1237 controls the stratification-dependence of the pycnocline-intensified dissipation. 1221 1238 It can take values of $1$ (recommended) or $2$. … … 1229 1246 1230 1247 % ================================================================ 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.})} 1248 % surface wave-induced mixing 1249 % ================================================================ 1250 \section[Surface wave-induced mixing (\forcode{ln_zdfswm})]{Surface wave-induced mixing (\protect\np{ln\_zdfswm})} 1235 1251 \label{subsec:ZDF_swm} 1236 1252 1237 TBC ... 1253 Surface waves produce an enhanced mixing through wave-turbulence interaction. 1254 In addition to breaking waves induced turbulence (\autoref{subsec:ZDF_tke}), 1255 the influence of non-breaking waves can be accounted introducing 1256 wave-induced viscosity and diffusivity as a function of the wave number spectrum. 1257 Following \citet{qiao.yuan.ea_OD10}, a formulation of wave-induced mixing coefficient 1258 is provided as a function of wave amplitude, Stokes Drift and wave-number: 1259 1260 \begin{equation} 1261 \label{eq:ZDF_Bv} 1262 B_{v} = \alpha {A} {U}_{st} {exp(3kz)} 1263 \end{equation} 1264 1265 Where $B_{v}$ is the wave-induced mixing coefficient, $A$ is the wave amplitude, 1266 ${U}_{st}$ is the Stokes Drift velocity, $k$ is the wave number and $\alpha$ 1267 is a constant which should be determined by observations or 1268 numerical experiments and is set to be 1. 1269 1270 The coefficient $B_{v}$ is then directly added to the vertical viscosity 1271 and diffusivity coefficients. 1272 1273 In order to account for this contribution set: \forcode{ln_zdfswm=.true.}, 1274 then wave interaction has to be activated through \forcode{ln_wave=.true.}, 1275 the Stokes Drift can be evaluated by setting \forcode{ln_sdw=.true.} 1276 (see \autoref{subsec:SBC_wave_sdw}) 1277 and the needed wave fields can be provided either in forcing or coupled mode 1278 (for more information on wave parameters and settings see \autoref{sec:SBC_wave}) 1238 1279 1239 1280 % ================================================================ 1240 1281 % Adaptive-implicit vertical advection 1241 1282 % ================================================================ 1242 \section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp = .true.})] 1243 {Adaptive-implicit vertical advection(\protect\np{ln\_zad\_Aimp}\forcode{ = .true.})} 1283 \section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp})]{Adaptive-implicit vertical advection(\protect\np{ln\_zad\_Aimp})} 1244 1284 \label{subsec:ZDF_aimp} 1245 1285 1246 This refers to \citep{shchepetkin_OM15}. 1247 1248 TBC ... 1249 1250 1286 The adaptive-implicit vertical advection option in NEMO is based on the work of 1287 \citep{shchepetkin_OM15}. In common with most ocean models, the timestep used with NEMO 1288 needs to satisfy multiple criteria associated with different physical processes in order 1289 to maintain numerical stability. \citep{shchepetkin_OM15} pointed out that the vertical 1290 CFL criterion is commonly the most limiting. \citep{lemarie.debreu.ea_OM15} examined the 1291 constraints for a range of time and space discretizations and provide the CFL stability 1292 criteria for a range of advection schemes. The values for the Leap-Frog with Robert 1293 asselin filter time-stepping (as used in NEMO) are reproduced in 1294 \autoref{tab:ZDF_zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these 1295 restrictions but at the cost of large dispersive errors and, possibly, large numerical 1296 viscosity. The adaptive-implicit vertical advection option provides a targetted use of the 1297 implicit scheme only when and where potential breaches of the vertical CFL condition 1298 occur. In many practical applications these events may occur remote from the main area of 1299 interest or due to short-lived conditions such that the extra numerical diffusion or 1300 viscosity does not greatly affect the overall solution. With such applications, setting: 1301 \forcode{ln_zad_Aimp=.true.} should allow much longer model timesteps to be used whilst 1302 retaining the accuracy of the high order explicit schemes over most of the domain. 1303 1304 \begin{table}[htbp] 1305 \centering 1306 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}} 1307 \begin{tabular}{r|ccc} 1308 \hline 1309 spatial discretization & 2$^nd$ order centered & 3$^rd$ order upwind & 4$^th$ order compact \\ 1310 advective CFL criterion & 0.904 & 0.472 & 0.522 \\ 1311 \hline 1312 \end{tabular} 1313 \caption[Advective CFL criteria for the leapfrog with Robert Asselin filter time-stepping]{ 1314 The advective CFL criteria for a range of spatial discretizations for 1315 the leapfrog with Robert Asselin filter time-stepping 1316 ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}.} 1317 \label{tab:ZDF_zad_Aimp_CFLcrit} 1318 \end{table} 1319 1320 In particular, the advection scheme remains explicit everywhere except where and when 1321 local vertical velocities exceed a threshold set just below the explicit stability limit. 1322 Once the threshold is reached a tapered transition towards an implicit scheme is used by 1323 partitioning the vertical velocity into a part that can be treated explicitly and any 1324 excess that must be treated implicitly. The partitioning is achieved via a Courant-number 1325 dependent weighting algorithm as described in \citep{shchepetkin_OM15}. 1326 1327 The local cell Courant number ($Cu$) used for this partitioning is: 1328 1329 \begin{equation} 1330 \label{eq:ZDF_Eqn_zad_Aimp_Courant} 1331 \begin{split} 1332 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 ] \\ 1333 &\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 ] 1334 \big / e_{{1_t}ij}e_{{2_t}ij} \\ 1335 &\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 ] 1336 \big / e_{{1_t}ij}e_{{2_t}ij} \bigg ) \\ 1337 \end{split} 1338 \end{equation} 1339 1340 \noindent and the tapering algorithm follows \citep{shchepetkin_OM15} as: 1341 1342 \begin{align} 1343 \label{eq:ZDF_Eqn_zad_Aimp_partition} 1344 Cu_{min} &= 0.15 \nonumber \\ 1345 Cu_{max} &= 0.3 \nonumber \\ 1346 Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ 1347 Fcu &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ 1348 \cf &= 1349 \begin{cases} 1350 0.0 &\text{if $Cu \leq Cu_{min}$} \\ 1351 (Cu - Cu_{min})^2 / (Fcu + (Cu - Cu_{min})^2) &\text{else if $Cu < Cu_{cut}$} \\ 1352 (Cu - Cu_{max}) / Cu &\text{else} 1353 \end{cases} 1354 \end{align} 1355 1356 \begin{figure}[!t] 1357 \centering 1358 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_coeff} 1359 \caption[Partitioning coefficient used to partition vertical velocities into parts]{ 1360 The value of the partitioning coefficient (\cf) used to partition vertical velocities into 1361 parts to be treated implicitly and explicitly for a range of typical Courant numbers 1362 (\forcode{ln_zad_Aimp=.true.}).} 1363 \label{fig:ZDF_zad_Aimp_coeff} 1364 \end{figure} 1365 1366 \noindent The partitioning coefficient is used to determine the part of the vertical 1367 velocity that must be handled implicitly ($w_i$) and to subtract this from the total 1368 vertical velocity ($w_n$) to leave that which can continue to be handled explicitly: 1369 1370 \begin{align} 1371 \label{eq:ZDF_Eqn_zad_Aimp_partition2} 1372 w_{i_{ijk}} &= \cf_{ijk} w_{n_{ijk}} \nonumber \\ 1373 w_{n_{ijk}} &= (1-\cf_{ijk}) w_{n_{ijk}} 1374 \end{align} 1375 1376 \noindent Note that the coefficient is such that the treatment is never fully implicit; 1377 the three cases from \autoref{eq:ZDF_Eqn_zad_Aimp_partition} can be considered as: 1378 fully-explicit; mixed explicit/implicit and mostly-implicit. With the settings shown the 1379 coefficient (\cf) varies as shown in \autoref{fig:ZDF_zad_Aimp_coeff}. Note with these values 1380 the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly 1381 implicit' is 0.45 which is just below the stability limited given in 1382 \autoref{tab:ZDF_zad_Aimp_CFLcrit} for a 3rd order scheme. 1383 1384 The $w_i$ component is added to the implicit solvers for the vertical mixing in 1385 \mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}. This is 1386 sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further 1387 intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}). 1388 For these schemes the implicit upstream fluxes must be added to both the monotonic guess 1389 and to the higher order solution when calculating the antidiffusive fluxes. The implicit 1390 vertical fluxes are then removed since they are added by the implicit solver later on. 1391 1392 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 1393 used in a wide range of simulations. The following test simulation, however, does illustrate 1394 the potential benefits and will hopefully encourage further testing and feedback from users: 1395 1396 \begin{figure}[!t] 1397 \centering 1398 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} 1399 \caption[OVERFLOW: time-series of temperature vertical cross-sections]{ 1400 A time-series of temperature vertical cross-sections for the OVERFLOW test case. 1401 These results are for the default settings with \forcode{nn_rdt=10.0} and 1402 without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}).} 1403 \label{fig:ZDF_zad_Aimp_overflow_frames} 1404 \end{figure} 1405 1406 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} 1407 1408 The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case} 1409 provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case 1410 by only a few extra physics choices namely: 1411 1412 \begin{verbatim} 1413 ln_dynldf_OFF = .false. 1414 ln_dynldf_lap = .true. 1415 ln_dynldf_hor = .true. 1416 ln_zdfnpc = .true. 1417 ln_traadv_fct = .true. 1418 nn_fct_h = 2 1419 nn_fct_v = 2 1420 \end{verbatim} 1421 1422 \noindent which were chosen to provide a slightly more stable and less noisy solution. The 1423 result when using the default value of \forcode{nn_rdt=10.} without adaptive-implicit 1424 vertical velocity is illustrated in \autoref{fig:ZDF_zad_Aimp_overflow_frames}. The mass of 1425 cold water, initially sitting on the shelf, moves down the slope and forms a 1426 bottom-trapped, dense plume. Even with these extra physics choices the model is close to 1427 stability limits and attempts with \forcode{nn_rdt=30.} will fail after about 5.5 hours 1428 with excessively high horizontal velocities. This time-scale corresponds with the time the 1429 plume reaches the steepest part of the topography and, although detected as a horizontal 1430 CFL breach, the instability originates from a breach of the vertical CFL limit. This is a good 1431 candidate, therefore, for use of the adaptive-implicit vertical advection scheme. 1432 1433 The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps 1434 are shown in \autoref{fig:ZDF_zad_Aimp_overflow_all_rdt} (together with the equivalent 1435 frames from the base run). In this simple example the use of the adaptive-implicit 1436 vertcal advection scheme has enabled a 12x increase in the model timestep without 1437 significantly altering the solution (although at this extreme the plume is more diffuse 1438 and has not travelled so far). Notably, the solution with and without the scheme is 1439 slightly different even with \forcode{nn_rdt=10.}; suggesting that the base run was 1440 close enough to instability to trigger the scheme despite completing successfully. 1441 To assist in diagnosing how active the scheme is, in both location and time, the 3D 1442 implicit and explicit components of the vertical velocity are available via XIOS as 1443 \texttt{wimp} and \texttt{wexp} respectively. Likewise, the partitioning coefficient 1444 (\cf) is also available as \texttt{wi\_cff}. For a quick oversight of 1445 the schemes activity the global maximum values of the absolute implicit component 1446 of the vertical velocity and the partitioning coefficient are written to the netCDF 1447 version of the run statistics file (\texttt{run.stat.nc}) if this is active (see 1448 \autoref{sec:MISC_opt} for activation details). 1449 1450 \autoref{fig:ZDF_zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for 1451 the various overflow tests. Note that the adaptive-implicit vertical advection scheme is 1452 active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the 1453 test case is close to stability limits even with this value. At the larger timesteps, the 1454 vertical velocity is treated mostly implicitly at some location throughout the run. The 1455 oscillatory nature of this measure appears to be linked to the progress of the plume front 1456 as each cusp is associated with the location of the maximum shifting to the adjacent cell. 1457 This is illustrated in \autoref{fig:ZDF_zad_Aimp_maxCf_loc} where the i- and k- locations of the 1458 maximum have been overlaid for the base run case. 1459 1460 \medskip 1461 \noindent Only limited tests have been performed in more realistic configurations. In the 1462 ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes 1463 restartability and reproducibility tests but it is unable to improve the model's stability 1464 enough to allow an increase in the model time-step. A view of the time-series of maximum 1465 partitioning coefficient (not shown here) suggests that the default time-step of 5400s is 1466 already pushing at stability limits, especially in the initial start-up phase. The 1467 time-series does not, however, exhibit any of the 'cuspiness' found with the overflow 1468 tests. 1469 1470 \medskip 1471 \noindent A short test with an eORCA1 configuration promises more since a test using a 1472 time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the 1473 time-step is limited to 2700s without. 1474 1475 \begin{figure}[!t] 1476 \centering 1477 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_overflow_all_rdt} 1478 \caption[OVERFLOW: sample temperature vertical cross-sections from mid- and end-run]{ 1479 Sample temperature vertical cross-sections from mid- and end-run using 1480 different values for \forcode{nn_rdt} and with or without adaptive implicit vertical advection. 1481 Without the adaptive implicit vertical advection 1482 only the run with the shortest timestep is able to run to completion. 1483 Note also that the colour-scale has been chosen to confirm that 1484 temperatures remain within the original range of 10$^o$ to 20$^o$.} 1485 \label{fig:ZDF_zad_Aimp_overflow_all_rdt} 1486 \end{figure} 1487 1488 \begin{figure}[!t] 1489 \centering 1490 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_maxCf} 1491 \caption[OVERFLOW: maximum partitioning coefficient during a series of test runs]{ 1492 The maximum partitioning coefficient during a series of test runs with 1493 increasing model timestep length. 1494 At the larger timesteps, 1495 the vertical velocity is treated mostly implicitly at some location throughout the run.} 1496 \label{fig:ZDF_zad_Aimp_maxCf} 1497 \end{figure} 1498 1499 \begin{figure}[!t] 1500 \centering 1501 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_maxCf_loc} 1502 \caption[OVERFLOW: maximum partitioning coefficient for the case overlaid]{ 1503 The maximum partitioning coefficient for the \forcode{nn_rdt=10.0} case overlaid with 1504 information on the gridcell i- and k-locations of the maximum value.} 1505 \label{fig:ZDF_zad_Aimp_maxCf_loc} 1506 \end{figure} 1251 1507 1252 1508 % ================================================================
Note: See TracChangeset
for help on using the changeset viewer.
