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