Changeset 11831 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/NEMO/subfiles/chap_ZDF.tex
- Timestamp:
- 2019-10-29T18:14:49+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc
- Property svn:ignore deleted
-
Property
svn:externals
set to
^/utils/badges badges
^/utils/logos logos
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex
- Property svn:ignore deleted
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/NEMO
- Property svn:ignore deleted
-
Property
svn:externals
set to
^/utils/figures/NEMO figures
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/NEMO/subfiles
- Property svn:ignore
-
old new 1 *.aux 2 *.bbl 3 *.blg 4 *.dvi 5 *.fdb* 6 *.fls 7 *.idx 1 *.ind 8 2 *.ilg 9 *.ind10 *.log11 *.maf12 *.mtc*13 *.out14 *.pdf15 *.toc16 _minted-*
-
- Property svn:ignore
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/NEMO/subfiles/chap_ZDF.tex
r10442 r11831 1 1 \documentclass[../main/NEMO_manual]{subfiles} 2 2 3 %% Custom aliases 4 \newcommand{\cf}{\ensuremath{C\kern-0.14em f}} 5 3 6 \begin{document} 4 % ================================================================ 5 % Chapter Vertical Ocean Physics (ZDF) 6 % ================================================================ 7 7 8 \chapter{Vertical Ocean Physics (ZDF)} 8 9 \label{chap:ZDF} 9 10 10 \minitoc 11 12 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. 13 14 \newpage 15 16 % ================================================================ 17 % Vertical Mixing 18 % ================================================================ 11 \thispagestyle{plain} 12 13 \chaptertoc 14 15 \paragraph{Changes record} ~\\ 16 17 {\footnotesize 18 \begin{tabularx}{\textwidth}{l||X|X} 19 Release & Author(s) & Modifications \\ 20 \hline 21 {\em 4.0} & {\em ...} & {\em ...} \\ 22 {\em 3.6} & {\em ...} & {\em ...} \\ 23 {\em 3.4} & {\em ...} & {\em ...} \\ 24 {\em <=3.4} & {\em ...} & {\em ...} 25 \end{tabularx} 26 } 27 28 \clearpage 29 30 \cmtgm{ Add here a small introduction to ZDF and naming of the different physics 31 (similar to what have been written for TRA and DYN).} 32 33 %% ================================================================================================= 19 34 \section{Vertical mixing} 20 \label{sec:ZDF _zdf}35 \label{sec:ZDF} 21 36 22 37 The discrete form of the ocean subgrid scale physics has been presented in … … 25 40 At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}), 26 41 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,42 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie\ \np{ln_trabbc}{ln\_trabbc} defined, 28 43 see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum 29 (see \autoref{sec:ZDF_ bfr}).44 (see \autoref{sec:ZDF_drg}). 30 45 31 46 In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and … … 33 48 respectively (see \autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}). 34 49 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.50 or computed from a turbulent closure model (either TKE or GLS or OSMOSIS formulation). 51 The computation of these coefficients is initialized in the \mdl{zdfphy} module and performed in 52 the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules. 38 53 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})} 54 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 55 %These trends can be computed using either a forward time stepping scheme 56 %(namelist parameter \np[=.true.]{ln_zdfexp}{ln\_zdfexp}) or a backward time stepping scheme 57 %(\np[=.false.]{ln_zdfexp}{ln\_zdfexp}) depending on the magnitude of the mixing coefficients, 58 %and thus of the formulation used (see \autoref{chap:TD}). 59 60 \begin{listing} 61 \nlst{namzdf} 62 \caption{\forcode{&namzdf}} 63 \label{lst:namzdf} 64 \end{listing} 65 66 %% ================================================================================================= 67 \subsection[Constant (\forcode{ln_zdfcst})]{Constant (\protect\np{ln_zdfcst}{ln\_zdfcst})} 49 68 \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 69 70 Options are defined through the \nam{zdf}{zdf} namelist variables. 71 When \np{ln_zdfcst}{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 57 72 constant values over the whole ocean. 58 73 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.74 It is recommended to use this option only in process studies, not in basin scale simulations. 60 75 Typical values used in this case are: 61 76 \begin{align*} … … 64 79 \end{align*} 65 80 66 These values are set through the \np{rn \_avm0} and \np{rn\_avt0} namelist parameters.81 These values are set through the \np{rn_avm0}{rn\_avm0} and \np{rn_avt0}{rn\_avt0} namelist parameters. 67 82 In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity, 68 83 that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and 69 84 $\sim10^{-9}~m^2.s^{-1}$ for salinity. 70 85 71 % ------------------------------------------------------------------------------------------------------------- 72 % Richardson Number Dependent 73 % ------------------------------------------------------------------------------------------------------------- 74 \subsection{Richardson number dependent (\protect\key{zdfric})} 86 %% ================================================================================================= 87 \subsection[Richardson number dependent (\forcode{ln_zdfric})]{Richardson number dependent (\protect\np{ln_zdfric}{ln\_zdfric})} 75 88 \label{subsec:ZDF_ric} 76 89 77 %--------------------------------------------namric--------------------------------------------------------- 78 79 \nlst{namzdf_ric} 80 %-------------------------------------------------------------------------------------------------------------- 81 82 When \key{zdfric} is defined, a local Richardson number dependent formulation for the vertical momentum and 83 tracer eddy coefficients is set through the \ngn{namzdf\_ric} namelist variables. 84 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. 90 \begin{listing} 91 \nlst{namzdf_ric} 92 \caption{\forcode{&namzdf_ric}} 93 \label{lst:namzdf_ric} 94 \end{listing} 95 96 When \np[=.true.]{ln_zdfric}{ln\_zdfric}, a local Richardson number dependent formulation for the vertical momentum and 97 tracer eddy coefficients is set through the \nam{zdf_ric}{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 103 117 can be reached by the coefficient when $Ri\leq 0$, $a=5$ and $n=2$. 104 The last three values can be modified by setting the \np{rn \_avmri}, \np{rn\_alp} and105 \np{nn \_ric} namelist parameters, respectively.118 The last three values can be modified by setting the \np{rn_avmri}{rn\_avmri}, \np{rn_alp}{rn\_alp} and 119 \np{nn_ric}{nn\_ric} namelist parameters, respectively. 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[=.true.]{ln_mldw}{ln\_mldw} 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 … … 122 136 \] 123 137 is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$. 124 The final $h_{e}$ is further constrained by the adjustable bounds \np{rn \_mldmin} and \np{rn\_mldmax}.138 The final $h_{e}$ is further constrained by the adjustable bounds \np{rn_mldmin}{rn\_mldmin} and \np{rn_mldmax}{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}{rn\_wtmix} and \np{rn_wvmix}{rn\_wvmix} \citep{lermusiaux_JMS01}. 141 142 %% ================================================================================================= 143 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke})]{TKE turbulent closure scheme (\protect\np{ln_zdftke}{ln\_zdftke})} 132 144 \label{subsec:ZDF_tke} 133 145 134 %--------------------------------------------namzdf_tke-------------------------------------------------- 135 136 \nlst{namzdf_tke} 137 %-------------------------------------------------------------------------------------------------------------- 146 \begin{listing} 147 \nlst{namzdf_tke} 148 \caption{\forcode{&namzdf_tke}} 149 \label{lst:namzdf_tke} 150 \end{listing} 138 151 139 152 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on 140 153 a prognostic equation for $\bar{e}$, the turbulent kinetic energy, 141 154 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 and155 This turbulent closure model has been developed by \citet{bougeault.lacarrere_MWR89} in the atmospheric case, 156 adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of \NEMO, 157 by \citet{blanke.delecluse_JPO93} for equatorial Atlantic simulations. 158 Since then, significant modifications have been introduced by \citet{madec.delecluse.ea_NPM98} in both the implementation and 146 159 the formulation of the mixing length scale. 147 160 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:161 its destruction through stratification, its vertical diffusion, and its dissipation of \citet{kolmogorov_IANS42} type: 149 162 \begin{equation} 150 \label{eq: zdftke_e}163 \label{eq:ZDF_tke_e} 151 164 \frac{\partial \bar{e}}{\partial t} = 152 165 \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2 … … 158 171 \end{equation} 159 172 \[ 160 % \label{eq: zdftke_kz}173 % \label{eq:ZDF_tke_kz} 161 174 \begin{split} 162 175 K_m &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ … … 164 177 \end{split} 165 178 \] 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, 179 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 180 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 168 181 $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients. 169 182 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}.171 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$:183 vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}. 184 They are set through namelist parameters \np{nn_ediff}{nn\_ediff} and \np{nn_ediss}{nn\_ediss}. 185 $P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$: 173 186 \begin{align*} 174 % \label{eq: prt}187 % \label{eq:ZDF_prt} 175 188 P_{rt} = 176 189 \begin{cases} … … 180 193 \end{cases} 181 194 \end{align*} 182 Options are defined through the \ngn{namzdfy\_tke} namelist variables. 183 The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable. 195 The choice of $P_{rt}$ is controlled by the \np{nn_pdl}{nn\_pdl} namelist variable. 184 196 185 197 At the sea surface, the value of $\bar{e}$ is prescribed from the wind stress field as 186 $\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 when198 $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter. 199 The default value of $e_{bb}$ is 3.75. \citep{gaspar.gregoris.ea_JGR90}), however a much larger value can be used when 188 200 taking into account the surface wave breaking (see below Eq. \autoref{eq:ZDF_Esbc}). 189 201 The bottom value of TKE is assumed to be equal to the value of the level just above. 190 202 The time integration of the $\bar{e}$ equation may formally lead to negative values because 191 203 the numerical scheme does not ensure its positivity. 192 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 in204 To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn_emin}{rn\_emin} namelist parameter). 205 Following \citet{gaspar.gregoris.ea_JGR90}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 206 This allows the subsequent formulations to match that of \citet{gargett_JMR84} for the diffusion in 195 207 the thermocline and deep ocean : $K_\rho = 10^{-3} / N$. 196 208 In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical instabilities associated with 197 209 too weak vertical diffusion. 198 They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and 199 \np{rn\_avt0} (namzdf namelist, see \autoref{subsec:ZDF_cst}). 200 210 They must be specified at least larger than the molecular values, and are set through \np{rn_avm0}{rn\_avm0} and 211 \np{rn_avt0}{rn\_avt0} (\nam{zdf}{zdf} namelist, see \autoref{subsec:ZDF_cst}). 212 213 %% ================================================================================================= 201 214 \subsubsection{Turbulent length scale} 202 215 203 216 For computational efficiency, the original formulation of the turbulent length scales proposed by 204 \citet{ Gaspar1990} has been simplified.205 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}:217 \citet{gaspar.gregoris.ea_JGR90} has been simplified. 218 Four formulations are proposed, the choice of which is controlled by the \np{nn_mxl}{nn\_mxl} namelist parameter. 219 The first two are based on the following first order approximation \citep{blanke.delecluse_JPO93}: 207 220 \begin{equation} 208 \label{eq: tke_mxl0_1}221 \label{eq:ZDF_tke_mxl0_1} 209 222 l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N 210 223 \end{equation} 211 224 which is valid in a stable stratified region with constant values of the Brunt-Vais\"{a}l\"{a} frequency. 212 225 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:226 (\np[=0]{nn_mxl}{nn\_mxl}) or by the local vertical scale factor (\np[=1]{nn_mxl}{nn\_mxl}). 227 \citet{blanke.delecluse_JPO93} notice that this simplification has two major drawbacks: 215 228 it makes no sense for locally unstable stratification and the computation no longer uses all 216 229 the information contained in the vertical density profile. 217 To overcome these drawbacks, \citet{ Madec1998} introduces the \np{nn\_mxl}\forcode{ = 2..3} cases,230 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np[=2, 3]{nn_mxl}{nn\_mxl} cases, 218 231 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:232 So, the length scales are first evaluated as in \autoref{eq:ZDF_tke_mxl0_1} and then bounded such that: 220 233 \begin{equation} 221 \label{eq: tke_mxl_constraint}234 \label{eq:ZDF_tke_mxl_constraint} 222 235 \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1 223 236 \qquad \text{with }\ l = l_k = l_\epsilon 224 237 \end{equation} 225 \autoref{eq: tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than238 \autoref{eq:ZDF_tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than 226 239 the variations of depth. 227 It provides a better approximation of the \citet{ Gaspar1990} formulation while being much less240 It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less 228 241 time consuming. 229 242 In particular, it allows the length scale to be limited not only by the distance to the surface or 230 243 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:244 the thermocline (\autoref{fig:ZDF_mixing_length}). 245 In order to impose the \autoref{eq:ZDF_tke_mxl_constraint} constraint, we introduce two additional length scales: 233 246 $l_{up}$ and $l_{dwn}$, the upward and downward length scales, and 234 247 evaluate the dissipation and mixing length scales as 235 248 (and note that here we use numerical indexing): 236 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>237 249 \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} 250 \centering 251 \includegraphics[width=0.66\textwidth]{ZDF_mixing_length} 252 \caption[Mixing length computation]{Illustration of the mixing length computation} 253 \label{fig:ZDF_mixing_length} 245 254 \end{figure} 246 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 247 \[ 248 % \label{eq:tke_mxl2} 255 \[ 256 % \label{eq:ZDF_tke_mxl2} 249 257 \begin{aligned} 250 258 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) … … 254 262 \end{aligned} 255 263 \] 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}264 where $l^{(k)}$ is computed using \autoref{eq:ZDF_tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. 265 266 In the \np[=2]{nn_mxl}{nn\_mxl} case, the dissipation and mixing length scales take the same value: 267 $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np[=3]{nn_mxl}{nn\_mxl} case, 268 the dissipation and mixing turbulent length scales are give as in \citet{gaspar.gregoris.ea_JGR90}: 269 \[ 270 % \label{eq:ZDF_tke_mxl_gaspar} 263 271 \begin{aligned} 264 272 & l_k = \sqrt{\ l_{up} \ \ l_{dwn}\ } \\ … … 267 275 \] 268 276 269 At the ocean surface, a non zero length scale is set through the \np{rn \_mxl0} namelist parameter.277 At the ocean surface, a non zero length scale is set through the \np{rn_mxl0}{rn\_mxl0} namelist parameter. 270 278 Usually the surface scale is given by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant and 271 279 $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}.280 Assuming $z_o=0.1$~m \citep{craig.banner_JPO94} leads to a 0.04~m, the default value of \np{rn_mxl0}{rn\_mxl0}. 273 281 In the ocean interior a minimum length scale is set to recover the molecular viscosity when 274 282 $\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ). 275 283 284 %% ================================================================================================= 276 285 \subsubsection{Surface wave breaking parameterization} 277 %-----------------------------------------------------------------------% 278 279 Following \citet{Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified to 286 287 Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to 280 288 include the effect of surface wave breaking energetics. 281 289 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 :290 The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and 291 air-sea drag coefficient. 292 The latter concerns the bulk formulae and is not discussed here. 293 294 Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is : 287 295 \begin{equation} 288 296 \label{eq:ZDF_Esbc} 289 297 \bar{e}_o = \frac{1}{2}\,\left( 15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o} 290 298 \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}.299 where $\alpha_{CB}$ is the \citet{craig.banner_JPO94} constant of proportionality which depends on the ''wave age'', 300 ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}. 293 301 The boundary condition on the turbulent length scale follows the Charnock's relation: 294 302 \begin{equation} … … 297 305 \end{equation} 298 306 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, and307 \citet{mellor.blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by 308 \citet{stacey_JPO99} citing observation evidence, and 301 309 $\alpha_{CB} = 100$ the Craig and Banner's value. 302 310 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} corresponds311 with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter, setting \np[=67.83]{rn_ebb}{rn\_ebb} corresponds 304 312 to $\alpha_{CB} = 100$. 305 Further setting \np {ln\_mxl0} to true applies \autoref{eq:ZDF_Lsbc} as surface boundary condition onlength scale,313 Further setting \np[=.true.]{ln_mxl0}{ln\_mxl0}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale, 306 314 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 on315 Note that a minimal threshold of \np{rn_emin0}{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the 308 316 surface $\bar{e}$ value. 309 317 318 %% ================================================================================================= 310 319 \subsubsection{Langmuir cells} 311 %--------------------------------------%312 320 313 321 Langmuir circulations (LC) can be described as ordered large-scale vertical motions in … … 315 323 Although LC have nothing to do with convection, the circulation pattern is rather similar to 316 324 so-called convective rolls in the atmospheric boundary layer. 317 The detailed physics behind LC is described in, for example, \citet{ Craik_Leibovich_JFM76}.325 The detailed physics behind LC is described in, for example, \citet{craik.leibovich_JFM76}. 318 326 The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and 319 wind drift currents. 327 wind drift currents. 320 328 321 329 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.330 \citep{axell_JGR02} for a $k-\epsilon$ turbulent closure. 323 331 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 : 332 an extra source term of TKE, $P_{LC}$. 333 The presence of $P_{LC}$ in \autoref{eq:ZDF_tke_e}, the TKE equation, is controlled by setting \np{ln_lc}{ln\_lc} to 334 \forcode{.true.} in the \nam{zdf_tke}{zdf\_tke} namelist. 335 336 By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}), 337 $P_{LC}$ is assumed to be : 330 338 \[ 331 339 P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}} 332 340 \] 333 341 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 as342 With no information about the wave field, $w_{LC}$ is assumed to be proportional to 343 the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module 344 \footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as 337 345 $u_s = 0.016 \,|U_{10m}|$. 338 346 Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of … … 341 349 For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at 342 350 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). 351 and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures). 344 352 The resulting expression for $w_{LC}$ is : 345 353 \[ … … 350 358 \end{cases} 351 359 \] 352 where $c_{LC} = 0.15$ has been chosen by \citep{ Axell_JGR02} as a good compromise to fit LES data.360 where $c_{LC} = 0.15$ has been chosen by \citep{axell_JGR02} as a good compromise to fit LES data. 353 361 The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second. 354 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}.362 The value of $c_{LC}$ is set through the \np{rn_lc}{rn\_lc} namelist parameter, 363 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}. 356 364 357 365 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 366 $H_{LC}$ is the depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by 367 converting its kinetic energy to potential energy, according to 360 368 \[ 361 369 - \int_{-H_{LC}}^0 { N^2\;z \;dz} = \frac{1}{2} u_s^2 362 370 \] 363 371 372 %% ================================================================================================= 364 373 \subsubsection{Mixing just below the mixed layer} 365 %--------------------------------------------------------------%366 374 367 375 Vertical mixing parameterizations commonly used in ocean general circulation models tend to 368 376 produce mixed-layer depths that are too shallow during summer months and windy conditions. 369 377 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}),378 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}. 379 The parameterization is an empirical one, \ie\ not derived from theoretical considerations, 380 but rather is meant to account for observed processes that affect the density structure of 381 the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme 382 (\ie\ near-inertial oscillations and ocean swells and waves). 383 384 When using this parameterization (\ie\ when \np[=1]{nn_etau}{nn\_etau}), 377 385 the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations, 378 386 swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition, … … 383 391 \end{equation} 384 392 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 of393 penetrates in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of 386 394 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).388 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}) or395 (no penetration if $f_i=1$, \ie\ if the ocean is entirely covered by sea-ice). 396 The value of $f_r$, usually a few percents, is specified through \np{rn_efr}{rn\_efr} namelist parameter. 397 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np[=0]{nn_etau}{nn\_etau}) or 390 398 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}.399 (\np[=1]{nn_etau}{nn\_etau}). 400 401 Note that two other option exist, \np[=2, 3]{nn_etau}{nn\_etau}. 394 402 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.403 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean. 396 404 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). 405 They will be removed in the next release. 406 407 % This should be explain better below what this rn_eice parameter is meant for: 408 In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn_eice}{rn\_eice} namelist parameter. 409 This parameter varies from \forcode{0} for no effect to \forcode{4} to suppress the TKE input into the ocean when Sea Ice concentration 410 is greater than 25\%. 411 412 % from Burchard et al OM 2008 : 413 % the most critical process not reproduced by statistical turbulence models is the activity of 414 % internal waves and their interaction with turbulence. After the Reynolds decomposition, 415 % internal waves are in principle included in the RANS equations, but later partially 416 % excluded by the hydrostatic assumption and the model resolution. 417 % Thus far, the representation of internal wave mixing in ocean models has been relatively crude 418 % (\eg\ Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 419 420 %% ================================================================================================= 421 \subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls})]{GLS: Generic Length Scale (\protect\np{ln_zdfgls}{ln\_zdfgls})} 422 \label{subsec:ZDF_gls} 423 424 \begin{listing} 425 \nlst{namzdf_gls} 426 \caption{\forcode{&namzdf_gls}} 427 \label{lst:namzdf_gls} 428 \end{listing} 429 430 The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations: 431 one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale, 432 $\psi$ \citep{umlauf.burchard_JMR03, umlauf.burchard_CSR05}. 433 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 434 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:ZDF_GLS} allows to recover a number of 435 well-known turbulent closures ($k$-$kl$ \citep{mellor.yamada_RG82}, $k$-$\epsilon$ \citep{rodi_JGR87}, 436 $k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}). 437 The GLS scheme is given by the following set of equations: 438 \begin{equation} 439 \label{eq:ZDF_gls_e} 440 \frac{\partial \bar{e}}{\partial t} = 441 \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 442 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 443 -K_\rho \,N^2 444 +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right] 445 - \epsilon 446 \end{equation} 447 448 \[ 449 % \label{eq:ZDF_gls_psi} 450 \begin{split} 451 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ 452 \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 453 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 454 - C_3 \,K_\rho\,N^2 - C_2 \,\epsilon \,Fw \right\} \\ 455 &+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 } 456 \;\frac{\partial \psi}{\partial k}} \right]\; 457 \end{split} 458 \] 459 460 \[ 461 % \label{eq:ZDF_gls_kz} 462 \begin{split} 463 K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ 464 K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l 465 \end{split} 466 \] 467 468 \[ 469 % \label{eq:ZDF_gls_eps} 470 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; 471 \] 472 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and 473 $\epsilon$ the dissipation rate. 474 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of 475 the choice of the turbulence model. 476 Four different turbulent models are pre-defined (\autoref{tab:ZDF_GLS}). 477 They are made available through the \np{nn_clo}{nn\_clo} namelist parameter. 478 479 \begin{table}[htbp] 480 \centering 481 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} 482 \begin{tabular}{ccccc} 483 & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ 484 % & \citep{mellor.yamada_RG82} & \citep{rodi_JGR87} & \citep{wilcox_AJ88} & \\ 485 \hline 486 \hline 487 \np{nn_clo}{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ 488 \hline 489 $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ 490 $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ 491 $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ 492 $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ 493 $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ 494 $C_3$ & 1. & 1. & 1. & 1. \\ 495 $F_{wall}$ & Yes & -- & -- & -- \\ 496 \hline 497 \hline 498 \end{tabular} 499 \caption[Set of predefined GLS parameters or equivalently predefined turbulence models available]{ 500 Set of predefined GLS parameters, or equivalently predefined turbulence models available with 501 \protect\np[=.true.]{ln_zdfgls}{ln\_zdfgls} and controlled by 502 the \protect\np{nn_clos}{nn\_clos} namelist variable in \protect\nam{zdf_gls}{zdf\_gls}.} 503 \label{tab:ZDF_GLS} 504 \end{table} 505 506 In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of 507 the mixing length towards $\kappa z_b$ ($\kappa$ is the Von Karman constant and $z_b$ the rugosity length scale) value near physical boundaries 508 (logarithmic boundary layer law). 509 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88}, 510 or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01} 511 (\np[=0, 3]{nn_stab_func}{nn\_stab\_func}, resp.). 512 The value of $C_{0\mu}$ depends on the choice of the stability function. 513 514 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or 515 Neumann condition through \np{nn_bc_surf}{nn\_bc\_surf} and \np{nn_bc_bot}{nn\_bc\_bot}, resp. 516 As for TKE closure, the wave effect on the mixing is considered when 517 \np[ > 0.]{rn_crban}{rn\_crban} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}. 518 The \np{rn_crban}{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and 519 \np{rn_charn}{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}. 520 521 The $\psi$ equation is known to fail in stably stratified flows, and for this reason 522 almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy. 523 With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. 524 A value of $c_{lim} = 0.53$ is often used \citep{galperin.kantha.ea_JAS88}. 525 \cite{umlauf.burchard_CSR05} show that the value of the clipping factor is of crucial importance for 526 the entrainment depth predicted in stably stratified situations, 527 and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes. 528 The clipping is only activated if \np[=.true.]{ln_length_lim}{ln\_length\_lim}, 529 and the $c_{lim}$ is set to the \np{rn_clim_galp}{rn\_clim\_galp} value. 530 531 The time and space discretization of the GLS equations follows the same energetic consideration as for 532 the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}. 533 Evaluation of the 4 GLS turbulent closure schemes can be found in \citet{warner.sherwood.ea_OM05} in ROMS model and 534 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO\ model. 406 535 407 536 % ------------------------------------------------------------------------------------------------------------- 408 % TKE discretization considerations537 % OSM OSMOSIS BL Scheme 409 538 % ------------------------------------------------------------------------------------------------------------- 410 \subsection{TKE discretization considerations (\protect\key{zdftke})} 411 \label{subsec:ZDF_tke_ene} 412 413 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 539 \subsection[OSM: OSMOSIS boundary layer scheme (\forcode{ln_zdfosm = .true.})] 540 {OSM: OSMOSIS boundary layer scheme (\protect\np{ln_zdfosm}{ln\_zdfosm})} 541 \label{subsec:ZDF_osm} 542 543 \begin{listing} 544 \nlst{namzdf_osm} 545 \caption{\forcode{&namzdf_osm}} 546 \label{lst:namzdf_osm} 547 \end{listing} 548 549 %-------------------------------------------------------------------------------------------------------------- 550 \paragraph{Namelist choices} 551 Most of the namelist options refer to how to specify the Stokes 552 surface drift and penetration depth. There are three options: 553 \begin{description} 554 \item \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave} Default value in \texttt{namelist\_ref}. In this case the Stokes drift is 555 assumed to be parallel to the surface wind stress, with 556 magnitude consistent with a constant turbulent Langmuir number 557 $\mathrm{La}_t=$ \protect\np{rn_m_la} {rn\_m\_la} i.e.\ 558 $u_{s0}=\tau/(\mathrm{La}_t^2\rho_0)$. Default value of 559 \protect\np{rn_m_la}{rn\_m\_la} is 0.3. The Stokes penetration 560 depth $\delta = $ \protect\np{rn_osm_dstokes}{rn\_osm\_dstokes}; this has default value 561 of 5~m. 562 563 \item \protect\np[=1]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is 564 assumed to be parallel to the surface wind stress, with 565 magnitude as in the classical Pierson-Moskowitz wind-sea 566 spectrum. Significant wave height and 567 wave-mean period taken from this spectrum are used to calculate the Stokes penetration 568 depth, following the approach set out in \citet{breivik.janssen.ea_JPO14}. 569 570 \item \protect\np[=2]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is 571 taken from ECMWF wave model output, though only the component parallel 572 to the wind stress is retained. Significant wave height and 573 wave-mean period from ECMWF wave model output are used to calculate the Stokes penetration 574 depth, again following \citet{breivik.janssen.ea_JPO14}. 575 576 \end{description} 577 578 Others refer to the treatment of diffusion and viscosity beneath 579 the surface boundary layer: 580 \begin{description} 581 \item \protect\np{ln_kpprimix} {ln\_kpprimix} Default is \np{.true.}. Switches on KPP-style Ri \#-dependent 582 mixing below the surface boundary layer. If this is set 583 \texttt{.true.} the following variable settings are honoured: 584 \item \protect\np{rn_riinfty}{rn\_riinfty} Critical value of local Ri \# below which 585 shear instability increases vertical mixing from background value. 586 \item \protect\np{rn_difri} {rn\_difri} Maximum value of Ri \#-dependent mixing at $\mathrm{Ri}=0$. 587 \item \protect\np{ln_convmix}{ln\_convmix} If \texttt{.true.} then, where water column is unstable, specify 588 diffusivity equal to \protect\np{rn_dif_conv}{rn\_dif\_conv} (default value is 1 m~s$^{-2}$). 589 \end{description} 590 Diagnostic output is controlled by: 591 \begin{description} 592 \item \protect\np{ln_dia_osm}{ln\_dia\_osm} Default is \np{.false.}; allows XIOS output of OSMOSIS internal fields. 593 \end{description} 594 Obsolete namelist parameters include: 595 \begin{description} 596 \item \protect\np{ln_use_osm_la}\np{ln\_use\_osm\_la} With \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave}, 597 \protect\np{rn_osm_dstokes} {rn\_osm\_dstokes} is always used to specify the Stokes 598 penetration depth. 599 \item \protect\np{nn_ave} {nn\_ave} Choice of averaging method for KPP-style Ri \# 600 mixing. Not taken account of. 601 \item \protect\np{rn_osm_hbl0} {rn\_osm\_hbl0} Depth of initial boundary layer is now set 602 by a density criterion similar to that used in calculating \emph{hmlp} (output as \texttt{mldr10\_1}) in \mdl{zdfmxl}. 603 \end{description} 604 605 \subsubsection{Summary} 606 Much of the time the turbulent motions in the ocean surface boundary 607 layer (OSBL) are not given by 608 classical shear turbulence. Instead they are in a regime known as 609 `Langmuir turbulence', dominated by an 610 interaction between the currents and the Stokes drift of the surface waves \citep[e.g.][]{mcwilliams.ea_JFM97}. 611 This regime is characterised by strong vertical turbulent motion, and appears when the surface Stokes drift $u_{s0}$ is much greater than the friction velocity $u_{\ast}$. More specifically Langmuir turbulence is thought to be crucial where the turbulent Langmuir number $\mathrm{La}_{t}=(u_{\ast}/u_{s0}) > 0.4$. 612 613 The OSMOSIS model is fundamentally based on results of Large Eddy 614 Simulations (LES) of Langmuir turbulence and aims to fully describe 615 this Langmuir regime. The description in this section is of necessity incomplete and further details are available in Grant. A (2019); in prep. 616 617 The OSMOSIS turbulent closure scheme is a similarity-scale scheme in 618 the same spirit as the K-profile 619 parameterization (KPP) scheme of \citet{large.ea_RG97}. 620 A specified shape of diffusivity, scaled by the (OSBL) depth 621 $h_{\mathrm{BL}}$ and a turbulent velocity scale, is imposed throughout the 622 boundary layer 623 $-h_{\mathrm{BL}}<z<\eta$. The turbulent closure model 624 also includes fluxes of tracers and momentum that are``non-local'' (independent of the local property gradient). 625 626 Rather than the OSBL 627 depth being diagnosed in terms of a bulk Richardson number criterion, 628 as in KPP, it is set by a prognostic equation that is informed by 629 energy budget considerations reminiscent of the classical mixed layer 630 models of \citet{kraus.turner_tellus67}. 631 The model also includes an explicit parametrization of the structure 632 of the pycnocline (the stratified region at the bottom of the OSBL). 633 634 Presently, mixing below the OSBL is handled by the Richardson 635 number-dependent mixing scheme used in \citet{large.ea_RG97}. 636 637 Convective parameterizations such as described in \ref{sec:ZDF_conv} 638 below should not be used with the OSMOSIS-OBL model: instabilities 639 within the OSBL are part of the model, while instabilities below the 640 ML are handled by the Ri \# dependent scheme. 641 642 \subsubsection{Depth and velocity scales} 643 The model supposes a boundary layer of thickness $h_{\mathrm{bl}}$ enclosing a well-mixed layer of thickness $h_{\mathrm{ml}}$ and a relatively thin pycnocline at the base of thickness $\Delta h$; Fig.~\ref{fig: OSBL_structure} shows typical (a) buoyancy structure and (b) turbulent buoyancy flux profile for the unstable boundary layer (losing buoyancy at the surface; e.g.\ cooling). 414 644 \begin{figure}[!t] 415 645 \begin{center} 416 \includegraphics[width=1.00\textwidth]{Fig_ZDF_TKE_time_scheme}646 %\includegraphics[width=0.7\textwidth]{ZDF_OSM_structure_of_OSBL} 417 647 \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.648 \protect\label{fig: OSBL_structure} 649 The structure of the entraining boundary layer. (a) Mean buoyancy profile. (b) Profile of the buoyancy flux. 420 650 } 421 \end{center} 651 \end{center} 422 652 \end{figure} 423 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 653 The pycnocline in the OSMOSIS scheme is assumed to have a finite thickness, and may include a number of model levels. This means that the OSMOSIS scheme must parametrize both the thickness of the pycnocline, and the turbulent fluxes within the pycnocline. 654 655 Consideration of the power input by wind acting on the Stokes drift suggests that the Langmuir turbulence has velocity scale: 656 \begin{equation}\label{eq:w_La} 657 w_{*L}= \left(u_*^2 u_{s\,0}\right)^{1/3}; 658 \end{equation} 659 but at times the Stokes drift may be weak due to e.g.\ ice cover, short fetch, misalignment with the surface stress, etc.\ so a composite velocity scale is assumed for the stable (warming) boundary layer: 660 \begin{equation}\label{eq:composite-nu} 661 \nu_{\ast}= \left\{ u_*^3 \left[1-\exp(-.5 \mathrm{La}_t^2)\right]+w_{*L}^3\right\}^{1/3}. 662 \end{equation} 663 For the unstable boundary layer this is merged with the standard convective velocity scale $w_{*C}=\left(\overline{w^\prime b^\prime}_0 \,h_\mathrm{ml}\right)^{1/3}$, where $\overline{w^\prime b^\prime}_0$ is the upwards surface buoyancy flux, to give: 664 \begin{equation}\label{eq:vel-scale-unstable} 665 \omega_* = \left(\nu_*^3 + 0.5 w_{*C}^3\right)^{1/3}. 666 \end{equation} 667 668 \subsubsection{The flux gradient model} 669 The flux-gradient relationships used in the OSMOSIS scheme take the form: 670 % 671 \begin{equation}\label{eq:flux-grad-gen} 672 \overline{w^\prime\chi^\prime}=-K\frac{\partial\overline{\chi}}{\partial z} + N_{\chi,s} +N_{\chi,b} +N_{\chi,t}, 673 \end{equation} 674 % 675 where $\chi$ is a general variable and $N_{\chi,s}, N_{\chi,b} \mathrm{and} N_{\chi,t}$ are the non-gradient terms, and represent the effects of the different terms in the turbulent flux-budget on the transport of $\chi$. $N_{\chi,s}$ represents the effects that the Stokes shear has on the transport of $\chi$, $N_{\chi,b}$ the effect of buoyancy, and $N_{\chi,t}$ the effect of the turbulent transport. The same general form for the flux-gradient relationship is used to parametrize the transports of momentum, heat and salinity. 676 677 In terms of the non-dimensionalized depth variables 678 % 679 \begin{equation}\label{eq:sigma} 680 \sigma_{\mathrm{ml}}= -z/h_{\mathrm{ml}}; \;\sigma_{\mathrm{bl}}= -z/h_{\mathrm{bl}}, 681 \end{equation} 682 % 683 in unstable conditions the eddy diffusivity ($K_d$) and eddy viscosity ($K_\nu$) profiles are parametrized as: 684 % 685 \begin{align}\label{eq:diff-unstable} 686 K_d=&0.8\, \omega_*\, h_{\mathrm{ml}} \, \sigma_{\mathrm{ml}} \left(1-\beta_d \sigma_{\mathrm{ml}}\right)^{3/2} 687 \\\label{eq:visc-unstable} 688 K_\nu =& 0.3\, \omega_* \,h_{\mathrm{ml}}\, \sigma_{\mathrm{ml}} \left(1-\beta_\nu \sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right) 689 \end{align} 690 % 691 where $\beta_d$ and $\beta_\nu$ are parameters that are determined by matching Eqs \ref{eq:diff-unstable} and \ref{eq:visc-unstable} to the eddy diffusivity and viscosity at the base of the well-mixed layer, given by 692 % 693 \begin{equation}\label{eq:diff-wml-base} 694 K_{d,\mathrm{ml}}=K_{\nu,\mathrm{ml}}=\,0.16\,\omega_* \Delta h. 695 \end{equation} 696 % 697 For stable conditions the eddy diffusivity/viscosity profiles are given by: 698 % 699 \begin{align}\label{diff-stable} 700 K_d= & 0.75\,\, \nu_*\, h_{\mathrm{ml}}\,\, \exp\left[-2.8 \left(h_{\mathrm{bl}}/L_L\right)^2\right]\sigma_{\mathrm{ml}} \left(1-\sigma_{\mathrm{ml}}\right)^{3/2} \\\label{eq:visc-stable} 701 K_\nu = & 0.375\,\, \nu_*\, h_{\mathrm{ml}} \,\, \exp\left[-2.8 \left(h_{\mathrm{bl}}/L_L\right)^2\right] \sigma_{\mathrm{ml}} \left(1-\sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right). 702 \end{align} 703 % 704 The shape of the eddy viscosity and diffusivity profiles is the same as the shape in the unstable OSBL. The eddy diffusivity/viscosity depends on the stability parameter $h_{\mathrm{bl}}/{L_L}$ where $ L_L$ is analogous to the Obukhov length, but for Langmuir turbulence: 705 \begin{equation}\label{eq:L_L} 706 L_L=-w_{*L}^3/\left<\overline{w^\prime b^\prime}\right>_L, 707 \end{equation} 708 with the mean turbulent buoyancy flux averaged over the boundary layer given in terms of its surface value $\overline{w^\prime b^\prime}_0$ and (downwards) )solar irradiance $I(z)$ by 709 \begin{equation} \label{eq:stable-av-buoy-flux} 710 \left<\overline{w^\prime b^\prime}\right>_L = \tfrac{1}{2} {\overline{w^\prime b^\prime}}_0-g\alpha_E\left[\tfrac{1}{2}(I(0)+I(-h))-\left<I\right>\right]. 711 \end{equation} 712 % 713 In unstable conditions the eddy diffusivity and viscosity depend on stability through the velocity scale $\omega_*$, which depends on the two velocity scales $\nu_*$ and $w_{*C}$. 714 715 Details of the non-gradient terms in \eqref{eq:flux-grad-gen} and of the fluxes within the pycnocline $-h_{\mathrm{bl}}<z<h_{\mathrm{ml}}$ can be found in Grant (2019). 716 717 \subsubsection{Evolution of the boundary layer depth} 718 719 The prognostic equation for the depth of the neutral/unstable boundary layer is given by \citep{grant+etal18}, 720 721 \begin{equation} \label{eq:dhdt-unstable} 722 %\frac{\partial h_\mathrm{bl}}{\partial t} + \mathbf{U}_b\cdot\nabla h_\mathrm{bl}= W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}} 723 \frac{\partial h_\mathrm{bl}}{\partial t} = W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}} 724 \end{equation} 725 where $h_\mathrm{bl}$ is the horizontally-varying depth of the OSBL, 726 $\mathbf{U}_b$ and $W_b$ are the mean horizontal and vertical 727 velocities at the base of the OSBL, ${\overline{w^\prime 728 b^\prime}}_\mathrm{ent}$ is the buoyancy flux due to entrainment 729 and $\Delta B_\mathrm{bl}$ is the difference between the buoyancy 730 averaged over the depth of the OSBL (i.e.\ including the ML and 731 pycnocline) and the buoyancy just below the base of the OSBL. This 732 equation for the case when the pycnocline has a finite thickness, 733 based on the potential energy budget of the OSBL, is the leading term 734 \citep{grant+etal18} of a generalization of that used in mixed-layer 735 models e.g.\ \citet{kraus.turner_tellus67}, in which the thickness of the pycnocline is taken to be zero. 736 737 The entrainment flux for the combination of convective and Langmuir turbulence is given by 738 \begin{equation} \label{eq:entrain-flux} 739 {\overline{w^\prime b^\prime}}_\mathrm{ent} = -\alpha_{\mathrm{B}} {\overline{w^\prime b^\prime}}_0 - \alpha_{\mathrm{S}} \frac{u_*^3}{h_{\mathrm{ml}}} 740 + G\left(\delta/h_{\mathrm{ml}} \right)\left[\alpha_{\mathrm{S}}e^{-1.5\, \mathrm{La}_t}-\alpha_{\mathrm{L}} \frac{w_{\mathrm{*L}}^3}{h_{\mathrm{ml}}}\right] 741 \end{equation} 742 where the factor $G\equiv 1 - \mathrm{e}^ {-25\delta/h_{\mathrm{bl}}}(1-4\delta/h_{\mathrm{bl}})$ models the lesser efficiency of Langmuir mixing when the boundary-layer depth is much greater than the Stokes depth, and $\alpha_{\mathrm{B}}$, $\alpha_{S}$ and $\alpha_{\mathrm{L}}$ depend on the ratio of the appropriate eddy turnover time to the inertial timescale $f^{-1}$. Results from the LES suggest $\alpha_{\mathrm{B}}=0.18 F(fh_{\mathrm{bl}}/w_{*C})$, $\alpha_{S}=0.15 F(fh_{\mathrm{bl}}/u_*$ and $\alpha_{\mathrm{L}}=0.035 F(fh_{\mathrm{bl}}/u_{*L})$, where $F(x)\equiv\tanh(x^{-1})^{0.69}$. 743 744 For the stable boundary layer, the equation for the depth of the OSBL is: 745 746 \begin{equation}\label{eq:dhdt-stable} 747 \max\left(\Delta B_{bl},\frac{w_{*L}^2}{h_\mathrm{bl}}\right)\frac{\partial h_\mathrm{bl}}{\partial t} = \left(0.06 + 0.52\,\frac{ h_\mathrm{bl}}{L_L}\right) \frac{w_{*L}^3}{h_\mathrm{bl}} +\left<\overline{w^\prime b^\prime}\right>_L. 748 \end{equation} 749 750 Equation. \ref{eq:dhdt-unstable} always leads to the depth of the entraining OSBL increasing (ignoring the effect of the mean vertical motion), but the change in the thickness of the stable OSBL given by Eq. \ref{eq:dhdt-stable} can be positive or negative, depending on the magnitudes of $\left<\overline{w^\prime b^\prime}\right>_L$ and $h_\mathrm{bl}/L_L$. The rate at which the depth of the OSBL can decrease is limited by choosing an effective buoyancy $w_{*L}^2/h_\mathrm{bl}$, in place of $\Delta B_{bl}$ which will be $\approx 0$ for the collapsing OSBL. 751 752 753 %% ================================================================================================= 754 \subsection[ Discrete energy conservation for TKE and GLS schemes]{Discrete energy conservation for TKE and GLS schemes} 755 \label{subsec:ZDF_tke_ene} 756 757 \begin{figure}[!t] 758 \centering 759 \includegraphics[width=0.66\textwidth]{ZDF_TKE_time_scheme} 760 \caption[Subgrid kinetic energy integration in GLS and TKE schemes]{ 761 Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and 762 its links to the momentum and tracer time integration.} 763 \label{fig:ZDF_TKE_time_scheme} 764 \end{figure} 424 765 425 766 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 how767 \autoref{eq:ZDF_tke_e}) and \autoref{eq:ZDF_gls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion 768 (first line in \autoref{eq:MB_zdf}). 769 To do so a special care has to be taken for both the time and space discretization of 770 the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}. 771 772 Let us first address the time stepping issue. \autoref{fig:ZDF_TKE_time_scheme} shows how 432 773 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.774 the one-level forward time stepping of the equation for $\bar{e}$. 434 775 With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to 435 776 the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and 436 summing the result vertically: 777 summing the result vertically: 437 778 \begin{equation} 438 \label{eq: energ1}779 \label{eq:ZDF_energ1} 439 780 \begin{split} 440 781 \int_{-H}^{\eta} u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt} \right) \,dz \\ … … 444 785 \end{equation} 445 786 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 at787 known at time $t$ (\autoref{fig:ZDF_TKE_time_scheme}), as it is required when using the TKE scheme 788 (see \autoref{sec:TD_forward_imp}). 789 The first term of the right hand side of \autoref{eq:ZDF_energ1} represents the kinetic energy transfer at 449 790 the surface (atmospheric forcing) and at the bottom (friction effect). 450 791 The second term is always negative. 451 792 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,793 \autoref{eq:ZDF_energ1} implies that, to be energetically consistent, 453 794 the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 454 795 ${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ … … 456 797 457 798 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}).799 (second term of the right hand side of \autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}). 459 800 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:801 The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by 802 multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically: 462 803 \begin{equation} 463 \label{eq: energ2}804 \label{eq:ZDF_energ2} 464 805 \begin{split} 465 806 \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt} \right) \,dz \\ … … 470 811 \end{split} 471 812 \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 because813 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 814 The first term of the right hand side of \autoref{eq:ZDF_energ2} is always zero because 474 815 there is no diffusive flux through the ocean surface and bottom). 475 816 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.817 Therefore \autoref{eq:ZDF_energ1} implies that, to be energetically consistent, 818 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 819 479 820 Let us now address the space discretization issue. 480 821 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}).822 the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:DOM_cell}). 482 823 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 by824 By redoing the \autoref{eq:ZDF_energ1} in the 3D case, it can be shown that the product of eddy coefficient by 484 825 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.826 Furthermore, the time variation of $e_3$ has be taken into account. 486 827 487 828 The above energetic considerations leads to the following final discrete form for the TKE equation: 488 829 \begin{equation} 489 \label{eq: zdftke_ene}830 \label{eq:ZDF_tke_ene} 490 831 \begin{split} 491 832 \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt} \equiv … … 504 845 \end{split} 505 846 \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}).847 where the last two terms in \autoref{eq:ZDF_tke_ene} (vertical diffusion and Kolmogorov dissipation) 848 are time stepped using a backward scheme (see\autoref{sec:TD_forward_imp}). 508 849 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 644 645 % ================================================================ 646 % Convection 647 % ================================================================ 850 %The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as 851 %they all appear in the right hand side of \autoref{eq:ZDF_tke_ene}. 852 %For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 853 854 %% ================================================================================================= 648 855 \section{Convection} 649 856 \label{sec:ZDF_conv} 650 857 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. 858 Static instabilities (\ie\ light potential densities under heavy ones) may occur at particular ocean grid points. 657 859 In nature, convective processes quickly re-establish the static stability of the water column. 658 860 These processes have been removed from the model via the hydrostatic assumption so they must be parameterized. … … 661 863 or/and the use of a turbulent closure scheme. 662 864 663 % ------------------------------------------------------------------------------------------------------------- 664 % Non-Penetrative Convective Adjustment 665 % ------------------------------------------------------------------------------------------------------------- 666 \subsection[Non-penetrative convective adjmt (\protect\np{ln\_tranpc}\forcode{ = .true.})] 667 {Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})} 865 %% ================================================================================================= 866 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc})]{Non-penetrative convective adjustment (\protect\np{ln_tranpc}{ln\_tranpc})} 668 867 \label{subsec:ZDF_npc} 669 868 670 %--------------------------------------------namzdf--------------------------------------------------------671 672 \nlst{namzdf}673 %--------------------------------------------------------------------------------------------------------------674 675 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>676 869 \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} 870 \centering 871 \includegraphics[width=0.66\textwidth]{ZDF_npc} 872 \caption[Unstable density profile treated by the non penetrative convective adjustment algorithm]{ 873 Example of an unstable density profile treated by 874 the non penetrative convective adjustment algorithm. 875 $1^{st}$ step: the initial profile is checked from the surface to the bottom. 876 It is found to be unstable between levels 3 and 4. 877 They are mixed. 878 The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed. 879 The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed. 880 The $1^{st}$ step ends since the density profile is then stable below the level 3. 881 $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step: 882 levels 2 to 5 are mixed. 883 The new density profile is checked. 884 It is found stable: end of algorithm.} 885 \label{fig:ZDF_npc} 694 886 \end{figure} 695 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 696 697 Options are defined through the \ngn{namzdf} namelist variables. 698 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}\forcode{ = .true.}. 699 It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of 887 888 Options are defined through the \nam{zdf}{zdf} namelist variables. 889 The non-penetrative convective adjustment is used when \np[=.true.]{ln_zdfnpc}{ln\_zdfnpc}. 890 It is applied at each \np{nn_npc}{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of 700 891 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}):892 (\ie\ until the mixed portion of the water column has \textit{exactly} the density of the water just below) 893 \citep{madec.delecluse.ea_JPO91}. 894 The associated algorithm is an iterative process used in the following way (\autoref{fig:ZDF_npc}): 704 895 starting from the top of the ocean, the first instability is found. 705 896 Assume in the following that the instability is located between levels $k$ and $k+1$. … … 716 907 This algorithm is significantly different from mixing statically unstable levels two by two. 717 908 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 than909 the algorithm used in \NEMO\ converges for any profile in a number of iterations which is less than 719 910 the number of vertical levels. 720 This property is of paramount importance as pointed out by \citet{ Killworth1989}:911 This property is of paramount importance as pointed out by \citet{killworth_iprc89}: 721 912 it avoids the existence of permanent and unrealistic static instabilities at the sea surface. 722 913 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}.914 the north-western Mediterranean Sea \citep{madec.delecluse.ea_JPO91, madec.chartier.ea_DAO91, madec.crepon_iprc91}. 724 915 725 916 The current implementation has been modified in order to deal with any non linear equation of seawater 726 917 (L. Brodeau, personnal communication). 727 918 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);919 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 920 (not the difference in potential density); 730 921 $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in 731 922 the same way their temperature and salinity has been mixed. … … 733 924 having to recompute the expansion coefficients at each mixing iteration. 734 925 735 % ------------------------------------------------------------------------------------------------------------- 736 % Enhanced Vertical Diffusion 737 % ------------------------------------------------------------------------------------------------------------- 738 \subsection{Enhanced vertical diffusion (\protect\np{ln\_zdfevd}\forcode{ = .true.})} 926 %% ================================================================================================= 927 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd})]{Enhanced vertical diffusion (\protect\np{ln_zdfevd}{ln\_zdfevd})} 739 928 \label{subsec:ZDF_evd} 740 929 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.}. 930 Options are defined through the \nam{zdf}{zdf} namelist variables. 931 The enhanced vertical diffusion parameterisation is used when \np[=.true.]{ln_zdfevd}{ln\_zdfevd}. 748 932 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},933 in regions where the stratification is unstable 934 (\ie\ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}. 935 This is done either on tracers only (\np[=0]{nn_evdm}{nn\_evdm}) or 936 on both momentum and tracers (\np[=1]{nn_evdm}{nn\_evdm}). 937 938 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np[=1]{nn_evdm}{nn\_evdm}, 755 939 the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ values also, are set equal to 756 the namelist parameter \np{rn \_avevd}.940 the namelist parameter \np{rn_avevd}{rn\_avevd}. 757 941 A typical value for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$. 758 942 This parameterisation of convective processes is less time consuming than 759 943 the convective adjustment algorithm presented above when mixing both tracers and 760 944 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 945 764 946 Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$. 765 947 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})} 948 a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:TD_mLF}). 949 950 %% ================================================================================================= 951 \subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf_}\{\forcode{tke,gls,osm}\})]{Handling convection with turbulent closure schemes (\forcode{ln_zdf{tke,gls,osm}})} 772 952 \label{subsec:ZDF_tcs} 773 953 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. 954 The turbulent closure schemes presented in \autoref{subsec:ZDF_tke}, \autoref{subsec:ZDF_gls} and 955 \autoref{subsec:ZDF_osm} (\ie\ \np{ln_zdftke}{ln\_zdftke} or \np{ln_zdfgls}{ln\_zdfgls} or \np{ln_zdfosm}{ln\_zdfosm} defined) deal, in theory, 956 with statically unstable density profiles. 776 957 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}$).958 \autoref{eq:ZDF_tke_e} or \autoref{eq:ZDF_gls_e} becomes a source term, since $N^2$ is negative. 959 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also of the four neighboring values at 960 velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$). 780 961 These large values restore the static stability of the water column in a way similar to that of 781 962 the enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}). … … 784 965 because the mixing length scale is bounded by the distance to the sea surface. 785 966 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.967 \ie\ setting the \np{ln_zdfnpc}{ln\_zdfnpc} namelist parameter to true and 968 defining the turbulent closure (\np{ln_zdftke}{ln\_zdftke} or \np{ln_zdfgls}{ln\_zdfgls} = \forcode{.true.}) all together. 969 970 The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection, 971 %as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp}, 972 therefore \np[=.false.]{ln_zdfevd}{ln\_zdfevd} should be used with the OSMOSIS scheme. 792 973 % gm% + one word on non local flux with KPP scheme trakpp.F90 module... 793 974 794 % ================================================================ 795 % Double Diffusion Mixing 796 % ================================================================ 797 \section{Double diffusion mixing (\protect\key{zdfddm})} 798 \label{sec:ZDF_ddm} 799 800 %-------------------------------------------namzdf_ddm------------------------------------------------- 801 % 975 %% ================================================================================================= 976 \section[Double diffusion mixing (\forcode{ln_zdfddm})]{Double diffusion mixing (\protect\np{ln_zdfddm}{ln\_zdfddm})} 977 \label{subsec:ZDF_ddm} 978 802 979 %\nlst{namzdf_ddm} 803 %-------------------------------------------------------------------------------------------------------------- 804 805 Options are defined through the \ngn{namzdf\_ddm} namelist variables.980 981 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter 982 \np{ln_zdfddm}{ln\_zdfddm} in \nam{zdf}{zdf}. 806 983 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. 807 984 The former condition leads to salt fingering and the latter to diffusive convection. 808 985 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 that986 \citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that 810 987 it leads to relatively minor changes in circulation but exerts significant regional influences on 811 988 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 989 990 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 815 991 \begin{align*} 816 % \label{eq: zdfddm_Kz}992 % \label{eq:ZDF_ddm_Kz} 817 993 &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\ 818 994 &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS} … … 826 1002 (1981): 827 1003 \begin{align} 828 \label{eq: zdfddm_f}1004 \label{eq:ZDF_ddm_f} 829 1005 A_f^{vS} &= 830 1006 \begin{cases} … … 832 1008 0 &\text{otherwise} 833 1009 \end{cases} 834 \\ \label{eq: zdfddm_f_T}835 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 1010 \\ \label{eq:ZDF_ddm_f_T} 1011 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 836 1012 \end{align} 837 1013 838 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>839 1014 \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} 1015 \centering 1016 \includegraphics[width=0.66\textwidth]{ZDF_ddm} 1017 \caption[Diapycnal diffusivities for temperature and salt in regions of salt fingering and 1018 diffusive convection]{ 1019 From \citet{merryfield.holloway.ea_JPO99}: 1020 (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in 1021 regions of salt fingering. 1022 Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and 1023 thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$; 1024 (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in 1025 regions of diffusive convection. 1026 Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation. 1027 The latter is not implemented in \NEMO.} 1028 \label{fig:ZDF_ddm} 853 1029 \end{figure} 854 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 855 856 The factor 0.7 in \autoref{eq:zdfddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx 0.7$ of 857 buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{McDougall_Taylor_JMR84}). 858 Following \citet{Merryfield1999}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 1030 1031 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 1032 buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{mcdougall.taylor_JMR84}). 1033 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 1034 860 1035 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by 861 Federov (1988) is used: 1036 Federov (1988) is used: 862 1037 \begin{align} 863 % \label{eq: zdfddm_d}1038 % \label{eq:ZDF_ddm_d} 864 1039 A_d^{vT} &= 865 1040 \begin{cases} … … 869 1044 \end{cases} 870 1045 \nonumber \\ 871 \label{eq: zdfddm_d_S}1046 \label{eq:ZDF_ddm_d_S} 872 1047 A_d^{vS} &= 873 1048 \begin{cases} … … 878 1053 \end{align} 879 1054 880 The dependencies of \autoref{eq: zdfddm_f} to \autoref{eq:zdfddm_d_S} on $R_\rho$ are illustrated in881 \autoref{fig: zdfddm}.1055 The dependencies of \autoref{eq:ZDF_ddm_f} to \autoref{eq:ZDF_ddm_d_S} on $R_\rho$ are illustrated in 1056 \autoref{fig:ZDF_ddm}. 882 1057 Implementing this requires computing $R_\rho$ at each grid point on every time step. 883 1058 This is done in \mdl{eosbn2} at the same time as $N^2$ is computed. 884 1059 This avoids duplication in the computation of $\alpha$ and $\beta$ (which is usually quite expensive). 885 1060 886 % ================================================================ 887 % Bottom Friction 888 % ================================================================ 889 \section{Bottom and top friction (\protect\mdl{zdfbfr})} 890 \label{sec:ZDF_bfr} 891 892 %--------------------------------------------nambfr-------------------------------------------------------- 893 % 894 %\nlst{nambfr} 895 %-------------------------------------------------------------------------------------------------------------- 896 897 Options to define the top and bottom friction are defined through the \ngn{nambfr} namelist variables. 1061 %% ================================================================================================= 1062 \section[Bottom and top friction (\textit{zdfdrg.F90})]{Bottom and top friction (\protect\mdl{zdfdrg})} 1063 \label{sec:ZDF_drg} 1064 1065 \begin{listing} 1066 \nlst{namdrg} 1067 \caption{\forcode{&namdrg}} 1068 \label{lst:namdrg} 1069 \end{listing} 1070 \begin{listing} 1071 \nlst{namdrg_top} 1072 \caption{\forcode{&namdrg_top}} 1073 \label{lst:namdrg_top} 1074 \end{listing} 1075 \begin{listing} 1076 \nlst{namdrg_bot} 1077 \caption{\forcode{&namdrg_bot}} 1078 \label{lst:namdrg_bot} 1079 \end{listing} 1080 1081 Options to define the top and bottom friction are defined through the \nam{drg}{drg} namelist variables. 898 1082 The bottom friction represents the friction generated by the bathymetry. 899 1083 The top friction represents the friction generated by the ice shelf/ocean interface. 900 As the friction processes at the top and bottom are treated in similar way, 901 only the bottom friction is described in detail below. 902 1084 As the friction processes at the top and the bottom are treated in and identical way, 1085 the description below considers mostly the bottom friction case, if not stated otherwise. 903 1086 904 1087 Both the surface momentum flux (wind stress) and the bottom momentum flux (bottom friction) enter the equations as 905 1088 a condition on the vertical diffusive flux. 906 1089 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 \]1090 \[ 1091 % \label{eq:ZDF_bfr_flux} 1092 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 1093 \] 911 1094 where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum outside 912 1095 the logarithmic turbulent boundary layer (thickness of the order of 1~m in the ocean). … … 916 1099 one needs a vertical diffusion coefficient $A^{vm} = 0.125$~m$^2$s$^{-1}$ 917 1100 (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. 1101 With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. 919 1102 When the vertical mixing coefficient is this small, using a flux condition is equivalent to 920 1103 entering the viscous forces (either wind stress or bottom friction) as a body force over the depth of the top or … … 922 1105 To illustrate this, consider the equation for $u$ at $k$, the last ocean level: 923 1106 \begin{equation} 924 \label{eq: zdfbfr_flux2}1107 \label{eq:ZDF_drg_flux2} 925 1108 \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 1109 \end{equation} … … 935 1118 936 1119 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}.1120 the general momentum trend in \mdl{dynzdf}. 938 1121 For the time-split surface pressure gradient algorithm, the momentum trend due to 939 1122 the barotropic component needs to be handled separately. 940 1123 For this purpose it is convenient to compute and store coefficients which can be simply combined with 941 1124 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:1125 These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where: 943 1126 \begin{equation} 944 \label{eq: zdfbfr_bdef}1127 \label{eq:ZDF_bfr_bdef} 945 1128 \frac{\partial {\textbf U_h}}{\partial t} = 946 1129 - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b 947 1130 \end{equation} 948 1131 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 949 950 % ------------------------------------------------------------------------------------------------------------- 951 % Linear Bottom Friction 952 % ------------------------------------------------------------------------------------------------------------- 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} 1132 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. 1133 1134 %% ================================================================================================= 1135 \subsection[Linear top/bottom friction (\forcode{ln_lin})]{Linear top/bottom friction (\protect\np{ln_lin}{ln\_lin})} 1136 \label{subsec:ZDF_drg_linear} 1137 1138 The linear friction parameterisation (including the special case of a free-slip condition) assumes that 1139 the friction is proportional to the interior velocity (\ie\ the velocity of the first/last model level): 1140 \[ 1141 % \label{eq:ZDF_bfr_linear} 960 1142 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 961 1143 \] 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, 1144 where $r$ is a friction coefficient expressed in $m s^{-1}$. 1145 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 964 1146 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}.1147 Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}. 966 1148 A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used in quasi-geostrophic models. 967 1149 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).1150 (\citet{gill_bk82}, Eq. 9.6.6). 969 1151 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 1152 and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 971 1153 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. 1154 It can be changed by specifying \np{rn_Uc0}{rn\_Uc0} (namelist parameter). 1155 1156 For the linear friction case the drag coefficient used in the general expression \autoref{eq:ZDF_bfr_bdef} is: 1157 \[ 1158 % \label{eq:ZDF_bfr_linbfr_b} 1159 c_b^T = - r 1160 \] 1161 When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}. 1162 Setting \np[=.true.]{ln_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition. 1163 1164 These values are assigned in \mdl{zdfdrg}. 1165 Note that there is support for local enhancement of these values via an externally defined 2D mask array 1166 (\np[=.true.]{ln_boost}{ln\_boost}) given in the \ifile{bfr\_coef} input NetCDF file. 988 1167 The mask values should vary from 0 to 1. 989 1168 Locations with a non-zero mask value will have the friction coefficient increased by 990 $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}. 991 992 % ------------------------------------------------------------------------------------------------------------- 993 % Non-Linear Bottom Friction 994 % ------------------------------------------------------------------------------------------------------------- 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} 1169 $mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}. 1170 1171 %% ================================================================================================= 1172 \subsection[Non-linear top/bottom friction (\forcode{ln_non_lin})]{Non-linear top/bottom friction (\protect\np{ln_non_lin}{ln\_non\_lin})} 1173 \label{subsec:ZDF_drg_nonlinear} 1174 1175 The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic: 1176 \[ 1177 % \label{eq:ZDF_drg_nonlinear} 1001 1178 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 1002 1179 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b 1003 1180 \] 1004 where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy due to tides,1181 where $C_D$ is a drag coefficient, and $e_b $ a top/bottom turbulent kinetic energy due to tides, 1005 1182 internal waves breaking and other short time scale currents. 1006 1183 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}$ and1184 As an example, the CME experiment \citep{treguier_JGR92} uses $C_D = 10^{-3}$ and 1185 $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 1186 $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}. 1030 1031 % ------------------------------------------------------------------------------------------------------------- 1032 % Bottom Friction Log-layer 1033 % ------------------------------------------------------------------------------------------------------------- 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: 1187 The CME choices have been set as default values (\np{rn_Cd0}{rn\_Cd0} and \np{rn_ke0}{rn\_ke0} namelist parameters). 1188 1189 As for the linear case, the friction is imposed in the code by adding the trend due to 1190 the friction to the general momentum trend in \mdl{dynzdf}. 1191 For the non-linear friction case the term computed in \mdl{zdfdrg} is: 1192 \[ 1193 % \label{eq:ZDF_drg_nonlinbfr} 1194 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} 1195 \] 1196 1197 The coefficients that control the strength of the non-linear friction are initialised as namelist parameters: 1198 $C_D$= \np{rn_Cd0}{rn\_Cd0}, and $e_b$ =\np{rn_bfeb2}{rn\_bfeb2}. 1199 Note that for applications which consider tides explicitly, a low or even zero value of \np{rn_bfeb2}{rn\_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array 1200 (\np[=.true.]{ln_boost}{ln\_boost}). 1201 This works in the same way as for the linear friction case with non-zero masked locations increased by 1202 $mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}. 1203 1204 %% ================================================================================================= 1205 \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer})]{Log-layer top/bottom friction (\protect\np{ln_loglayer}{ln\_loglayer})} 1206 \label{subsec:ZDF_drg_loglayer} 1207 1208 In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using 1209 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. 1210 If \np[=.true.]{ln_loglayer}{ln\_loglayer}, $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): 1211 \[ 1212 C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2 1213 \] 1214 1215 \noindent where $\kappa$ is the von-Karman constant and \np{rn_z0}{rn\_z0} is a roughness length provided via the namelist. 1216 1217 The drag coefficient is bounded such that it is kept greater or equal to 1218 the base \np{rn_Cd0}{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: 1219 \np{rn_Cdmax}{rn\_Cdmax}, \ie 1220 \[ 1221 rn\_Cd0 \leq C_D \leq rn\_Cdmax 1222 \] 1223 1224 \noindent The log-layer enhancement can also be applied to the top boundary friction if 1225 under ice-shelf cavities are activated (\np[=.true.]{ln_isfcav}{ln\_isfcav}). 1226 %In this case, the relevant namelist parameters are \np{rn_tfrz0}{rn\_tfrz0}, \np{rn_tfri2}{rn\_tfri2} and \np{rn_tfri2_max}{rn\_tfri2\_max}. 1227 1228 %% ================================================================================================= 1229 \subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np[=.false.]{ln_drgimp}{ln\_drgimp})} 1230 \label{subsec:ZDF_drg_stability} 1231 1232 Setting \np[=.false.]{ln_drgimp}{ln\_drgimp} 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: 1233 1234 At the top (below an ice shelf cavity): 1235 \[ 1236 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1237 = c_{t}^{\textbf{U}}\textbf{u}^{n-1}_{t} 1238 \] 1239 1240 At the bottom (above the sea floor): 1241 \[ 1242 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1243 = c_{b}^{\textbf{U}}\textbf{u}^{n-1}_{b} 1244 \] 1245 1246 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. 1247 For the purposes of stability analysis, an approximation to \autoref{eq:ZDF_drg_flux2} is: 1069 1248 \begin{equation} 1070 \label{eq: Eqn_bfrstab}1249 \label{eq:ZDF_Eqn_drgstab} 1071 1250 \begin{split} 1072 1251 \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt \\ … … 1074 1253 \end{split} 1075 1254 \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:1255 \noindent where linear friction and a leapfrog timestep have been assumed. 1256 To ensure that the friction cannot reverse the direction of flow it is necessary to have: 1257 \[ 1258 |\Delta u| < \;|u| 1259 \] 1260 \noindent which, using \autoref{eq:ZDF_Eqn_drgstab}, gives: 1082 1261 \[ 1083 1262 r\frac{2\rdt}{e_{3u}} < 1 \qquad \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ … … 1092 1271 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 1272 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 both1273 But caution may be necessary if attempts are made to locally enhance the bottom friction parameters. 1274 To ensure stability limits are imposed on the top/bottom friction coefficients both 1096 1275 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).1276 Checks at initialisation are made in \mdl{zdfdrg} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case). 1098 1277 The number of breaches of the stability criterion are reported as well as 1099 1278 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;1279 The criterion is also checked at each time step, using the actual velocity, in \mdl{dynzdf}. 1280 Values of the friction coefficient are reduced as necessary to ensure stability; 1102 1281 these changes are not reported. 1103 1282 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}).1283 Limits on the top/bottom friction coefficient are not imposed if the user has elected to 1284 handle the friction implicitly (see \autoref{subsec:ZDF_drg_imp}). 1106 1285 The number of potential breaches of the explicit stability criterion are still reported for information purposes. 1107 1286 1108 % ------------------------------------------------------------------------------------------------------------- 1109 % Implicit Bottom Friction 1110 % ------------------------------------------------------------------------------------------------------------- 1111 \subsection{Implicit bottom friction (\protect\np{ln\_bfrimp}\forcode{ = .true.})} 1112 \label{subsec:ZDF_bfr_imp} 1287 %% ================================================================================================= 1288 \subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np[=.true.]{ln_drgimp}{ln\_drgimp})} 1289 \label{subsec:ZDF_drg_imp} 1113 1290 1114 1291 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 1292 We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting. 1293 This option can be invoked by setting \np{ln_drgimp}{ln\_drgimp} to \forcode{.true.} in the \nam{drg}{drg} namelist. 1294 %This option requires \np{ln_zdfexp}{ln\_zdfexp} to be \forcode{.false.} in the \nam{zdf}{zdf} namelist. 1295 1296 This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step: 1297 1298 At the top (below an ice shelf cavity): 1299 \[ 1300 % \label{eq:ZDF_dynZDF__drg_top} 1301 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1302 = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t} 1303 \] 1304 1305 At the bottom (above the sea floor): 1306 \[ 1307 % \label{eq:ZDF_dynZDF__drg_bot} 1308 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1309 = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b} 1310 \] 1311 1312 where $t$ and $b$ refers to top and bottom layers respectively. 1313 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit. 1314 1315 %% ================================================================================================= 1316 \subsection[Bottom friction with split-explicit free surface]{Bottom friction with split-explicit free surface} 1317 \label{subsec:ZDF_drg_ts} 1318 1319 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[=.false.]{ln_drgimp}{ln\_drgimp} 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[=.true.]{ln_drgimp}{ln\_drgimp}, 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. 1320 1321 The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO\ is as follows: 1185 1322 \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. 1323 \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. 1324 \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 1325 \end{enumerate} 1196 1326 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-------------------------------------------------- 1224 % 1225 %\nlst{namzdf_tmx} 1226 %-------------------------------------------------------------------------------------------------------------- 1227 1228 1229 % ------------------------------------------------------------------------------------------------------------- 1230 % Bottom intensified tidal mixing 1231 % ------------------------------------------------------------------------------------------------------------- 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 by 1237 \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) and 1249 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 and 1253 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$ and 1268 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 of 1272 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 and 1274 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 of 1276 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 component 1278 (\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 treatment 1298 % ------------------------------------------------------------------------------------------------------------- 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 following 1308 the parameterisation developed by \citet{Koch-Larrouy_al_GRL07}: 1309 1310 First, the Indonesian archipelago is a complex geographic region with a series of 1311 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 of 1321 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 in 1336 the different Indonesian seas, suggesting that the horizontal and vertical distributions of 1337 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 of 1339 global coupled GCMs \citep{Koch-Larrouy_al_CD10}. 1340 1341 % ================================================================ 1342 % Internal wave-driven mixing 1343 % ================================================================ 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 %-------------------------------------------------------------------------------------------------------------- 1327 Note that other strategies are possible, like considering vertical diffusion step in advance, \ie\ prior barotropic integration. 1328 1329 %% ================================================================================================= 1330 \section[Internal wave-driven mixing (\forcode{ln_zdfiwm})]{Internal wave-driven mixing (\protect\np{ln_zdfiwm}{ln\_zdfiwm})} 1331 \label{subsec:ZDF_tmx_new} 1332 1333 \begin{listing} 1334 \nlst{namzdf_iwm} 1335 \caption{\forcode{&namzdf_iwm}} 1336 \label{lst:namzdf_iwm} 1337 \end{listing} 1351 1338 1352 1339 The parameterization of mixing induced by breaking internal waves is a generalization of 1353 the approach originally proposed by \citet{ St_Laurent_al_GRL02}.1340 the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}. 1354 1341 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}1342 and the resulting diffusivity is obtained as 1343 \[ 1344 % \label{eq:ZDF_Kwave} 1358 1345 A^{vT}_{wave} = R_f \,\frac{ \epsilon }{ \rho \, N^2 } 1359 1346 \] 1360 1347 where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of 1361 1348 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}.1349 If the \np{ln_mevar}{ln\_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and 1350 equal to 1/6 \citep{osborn_JPO80}. 1364 1351 In the opposite (recommended) case, $R_f$ is instead a function of 1365 1352 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}.1353 with $\nu$ the molecular viscosity of seawater, following the model of \cite{bouffard.boegman_DAO13} and 1354 the implementation of \cite{de-lavergne.madec.ea_JPO16}. 1368 1355 Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when 1369 1356 the mixing efficiency is constant. 1370 1357 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}.1358 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 1359 as a function of $Re_b$ by setting the \np{ln_tsdiff}{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice. 1360 This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14}, 1361 is implemented as in \cite{de-lavergne.madec.ea_JPO16}. 1375 1362 1376 1363 The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, 1377 1364 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): 1365 $E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures: 1366 1380 1367 \begin{align*} 1381 1368 F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\ 1382 F_{pyc}(i,j,k) &\propto N^{n \_p}\\1369 F_{pyc}(i,j,k) &\propto N^{n_p}\\ 1383 1370 F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } 1384 \end{align*} 1371 \end{align*} 1385 1372 In the above formula, $h_{ab}$ denotes the height above bottom, 1386 1373 $h_{wkb}$ denotes the WKB-stretched height above bottom, defined by … … 1388 1375 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz' } \; , 1389 1376 \] 1390 The $n_p$ parameter (given by \np{nn \_zpyc} in \ngn{namzdf\_tmx\_new} namelist)1377 The $n_p$ parameter (given by \np{nn_zpyc}{nn\_zpyc} in \nam{zdf_iwm}{zdf\_iwm} namelist) 1391 1378 controls the stratification-dependence of the pycnocline-intensified dissipation. 1392 It can take values of 1 (recommended) or 2.1379 It can take values of $1$ (recommended) or $2$. 1393 1380 Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of 1394 1381 the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps. 1395 1382 $h_{cri}$ is related to the large-scale topography of the ocean (etopo2) and 1396 1383 $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. 1398 1399 % ================================================================ 1400 1401 \biblio 1402 1403 \pindex 1384 the abyssal hill topography \citep{goff_JGR10} and the latitude. 1385 % Jc: input files names ? 1386 1387 %% ================================================================================================= 1388 \section[Surface wave-induced mixing (\forcode{ln_zdfswm})]{Surface wave-induced mixing (\protect\np{ln_zdfswm}{ln\_zdfswm})} 1389 \label{subsec:ZDF_swm} 1390 1391 Surface waves produce an enhanced mixing through wave-turbulence interaction. 1392 In addition to breaking waves induced turbulence (\autoref{subsec:ZDF_tke}), 1393 the influence of non-breaking waves can be accounted introducing 1394 wave-induced viscosity and diffusivity as a function of the wave number spectrum. 1395 Following \citet{qiao.yuan.ea_OD10}, a formulation of wave-induced mixing coefficient 1396 is provided as a function of wave amplitude, Stokes Drift and wave-number: 1397 1398 \begin{equation} 1399 \label{eq:ZDF_Bv} 1400 B_{v} = \alpha {A} {U}_{st} {exp(3kz)} 1401 \end{equation} 1402 1403 Where $B_{v}$ is the wave-induced mixing coefficient, $A$ is the wave amplitude, 1404 ${U}_{st}$ is the Stokes Drift velocity, $k$ is the wave number and $\alpha$ 1405 is a constant which should be determined by observations or 1406 numerical experiments and is set to be 1. 1407 1408 The coefficient $B_{v}$ is then directly added to the vertical viscosity 1409 and diffusivity coefficients. 1410 1411 In order to account for this contribution set: \forcode{ln_zdfswm=.true.}, 1412 then wave interaction has to be activated through \forcode{ln_wave=.true.}, 1413 the Stokes Drift can be evaluated by setting \forcode{ln_sdw=.true.} 1414 (see \autoref{subsec:SBC_wave_sdw}) 1415 and the needed wave fields can be provided either in forcing or coupled mode 1416 (for more information on wave parameters and settings see \autoref{sec:SBC_wave}) 1417 1418 %% ================================================================================================= 1419 \section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp})]{Adaptive-implicit vertical advection(\protect\np{ln_zad_Aimp}{ln\_zad\_Aimp})} 1420 \label{subsec:ZDF_aimp} 1421 1422 The adaptive-implicit vertical advection option in NEMO is based on the work of 1423 \citep{shchepetkin_OM15}. In common with most ocean models, the timestep used with NEMO 1424 needs to satisfy multiple criteria associated with different physical processes in order 1425 to maintain numerical stability. \citep{shchepetkin_OM15} pointed out that the vertical 1426 CFL criterion is commonly the most limiting. \citep{lemarie.debreu.ea_OM15} examined the 1427 constraints for a range of time and space discretizations and provide the CFL stability 1428 criteria for a range of advection schemes. The values for the Leap-Frog with Robert 1429 asselin filter time-stepping (as used in NEMO) are reproduced in 1430 \autoref{tab:ZDF_zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these 1431 restrictions but at the cost of large dispersive errors and, possibly, large numerical 1432 viscosity. The adaptive-implicit vertical advection option provides a targetted use of the 1433 implicit scheme only when and where potential breaches of the vertical CFL condition 1434 occur. In many practical applications these events may occur remote from the main area of 1435 interest or due to short-lived conditions such that the extra numerical diffusion or 1436 viscosity does not greatly affect the overall solution. With such applications, setting: 1437 \forcode{ln_zad_Aimp=.true.} should allow much longer model timesteps to be used whilst 1438 retaining the accuracy of the high order explicit schemes over most of the domain. 1439 1440 \begin{table}[htbp] 1441 \centering 1442 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}} 1443 \begin{tabular}{r|ccc} 1444 \hline 1445 spatial discretization & 2$^nd$ order centered & 3$^rd$ order upwind & 4$^th$ order compact \\ 1446 advective CFL criterion & 0.904 & 0.472 & 0.522 \\ 1447 \hline 1448 \end{tabular} 1449 \caption[Advective CFL criteria for the leapfrog with Robert Asselin filter time-stepping]{ 1450 The advective CFL criteria for a range of spatial discretizations for 1451 the leapfrog with Robert Asselin filter time-stepping 1452 ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}.} 1453 \label{tab:ZDF_zad_Aimp_CFLcrit} 1454 \end{table} 1455 1456 In particular, the advection scheme remains explicit everywhere except where and when 1457 local vertical velocities exceed a threshold set just below the explicit stability limit. 1458 Once the threshold is reached a tapered transition towards an implicit scheme is used by 1459 partitioning the vertical velocity into a part that can be treated explicitly and any 1460 excess that must be treated implicitly. The partitioning is achieved via a Courant-number 1461 dependent weighting algorithm as described in \citep{shchepetkin_OM15}. 1462 1463 The local cell Courant number ($Cu$) used for this partitioning is: 1464 1465 \begin{equation} 1466 \label{eq:ZDF_Eqn_zad_Aimp_Courant} 1467 \begin{split} 1468 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 ] \\ 1469 &\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 ] 1470 \big / e_{{1_t}ij}e_{{2_t}ij} \\ 1471 &\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 ] 1472 \big / e_{{1_t}ij}e_{{2_t}ij} \bigg ) \\ 1473 \end{split} 1474 \end{equation} 1475 1476 \noindent and the tapering algorithm follows \citep{shchepetkin_OM15} as: 1477 1478 \begin{align} 1479 \label{eq:ZDF_Eqn_zad_Aimp_partition} 1480 Cu_{min} &= 0.15 \nonumber \\ 1481 Cu_{max} &= 0.3 \nonumber \\ 1482 Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ 1483 Fcu &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ 1484 \cf &= 1485 \begin{cases} 1486 0.0 &\text{if $Cu \leq Cu_{min}$} \\ 1487 (Cu - Cu_{min})^2 / (Fcu + (Cu - Cu_{min})^2) &\text{else if $Cu < Cu_{cut}$} \\ 1488 (Cu - Cu_{max}) / Cu &\text{else} 1489 \end{cases} 1490 \end{align} 1491 1492 \begin{figure}[!t] 1493 \centering 1494 \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_coeff} 1495 \caption[Partitioning coefficient used to partition vertical velocities into parts]{ 1496 The value of the partitioning coefficient (\cf) used to partition vertical velocities into 1497 parts to be treated implicitly and explicitly for a range of typical Courant numbers 1498 (\forcode{ln_zad_Aimp=.true.}).} 1499 \label{fig:ZDF_zad_Aimp_coeff} 1500 \end{figure} 1501 1502 \noindent The partitioning coefficient is used to determine the part of the vertical 1503 velocity that must be handled implicitly ($w_i$) and to subtract this from the total 1504 vertical velocity ($w_n$) to leave that which can continue to be handled explicitly: 1505 1506 \begin{align} 1507 \label{eq:ZDF_Eqn_zad_Aimp_partition2} 1508 w_{i_{ijk}} &= \cf_{ijk} w_{n_{ijk}} \nonumber \\ 1509 w_{n_{ijk}} &= (1-\cf_{ijk}) w_{n_{ijk}} 1510 \end{align} 1511 1512 \noindent Note that the coefficient is such that the treatment is never fully implicit; 1513 the three cases from \autoref{eq:ZDF_Eqn_zad_Aimp_partition} can be considered as: 1514 fully-explicit; mixed explicit/implicit and mostly-implicit. With the settings shown the 1515 coefficient (\cf) varies as shown in \autoref{fig:ZDF_zad_Aimp_coeff}. Note with these values 1516 the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly 1517 implicit' is 0.45 which is just below the stability limited given in 1518 \autoref{tab:ZDF_zad_Aimp_CFLcrit} for a 3rd order scheme. 1519 1520 The $w_i$ component is added to the implicit solvers for the vertical mixing in 1521 \mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}. This is 1522 sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further 1523 intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}). 1524 For these schemes the implicit upstream fluxes must be added to both the monotonic guess 1525 and to the higher order solution when calculating the antidiffusive fluxes. The implicit 1526 vertical fluxes are then removed since they are added by the implicit solver later on. 1527 1528 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 1529 used in a wide range of simulations. The following test simulation, however, does illustrate 1530 the potential benefits and will hopefully encourage further testing and feedback from users: 1531 1532 \begin{figure}[!t] 1533 \centering 1534 \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_overflow_frames} 1535 \caption[OVERFLOW: time-series of temperature vertical cross-sections]{ 1536 A time-series of temperature vertical cross-sections for the OVERFLOW test case. 1537 These results are for the default settings with \forcode{nn_rdt=10.0} and 1538 without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}).} 1539 \label{fig:ZDF_zad_Aimp_overflow_frames} 1540 \end{figure} 1541 1542 %% ================================================================================================= 1543 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} 1544 1545 The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case} 1546 provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case 1547 by only a few extra physics choices namely: 1548 1549 \begin{verbatim} 1550 ln_dynldf_OFF = .false. 1551 ln_dynldf_lap = .true. 1552 ln_dynldf_hor = .true. 1553 ln_zdfnpc = .true. 1554 ln_traadv_fct = .true. 1555 nn_fct_h = 2 1556 nn_fct_v = 2 1557 \end{verbatim} 1558 1559 \noindent which were chosen to provide a slightly more stable and less noisy solution. The 1560 result when using the default value of \forcode{nn_rdt=10.} without adaptive-implicit 1561 vertical velocity is illustrated in \autoref{fig:ZDF_zad_Aimp_overflow_frames}. The mass of 1562 cold water, initially sitting on the shelf, moves down the slope and forms a 1563 bottom-trapped, dense plume. Even with these extra physics choices the model is close to 1564 stability limits and attempts with \forcode{nn_rdt=30.} will fail after about 5.5 hours 1565 with excessively high horizontal velocities. This time-scale corresponds with the time the 1566 plume reaches the steepest part of the topography and, although detected as a horizontal 1567 CFL breach, the instability originates from a breach of the vertical CFL limit. This is a good 1568 candidate, therefore, for use of the adaptive-implicit vertical advection scheme. 1569 1570 The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps 1571 are shown in \autoref{fig:ZDF_zad_Aimp_overflow_all_rdt} (together with the equivalent 1572 frames from the base run). In this simple example the use of the adaptive-implicit 1573 vertcal advection scheme has enabled a 12x increase in the model timestep without 1574 significantly altering the solution (although at this extreme the plume is more diffuse 1575 and has not travelled so far). Notably, the solution with and without the scheme is 1576 slightly different even with \forcode{nn_rdt=10.}; suggesting that the base run was 1577 close enough to instability to trigger the scheme despite completing successfully. 1578 To assist in diagnosing how active the scheme is, in both location and time, the 3D 1579 implicit and explicit components of the vertical velocity are available via XIOS as 1580 \texttt{wimp} and \texttt{wexp} respectively. Likewise, the partitioning coefficient 1581 (\cf) is also available as \texttt{wi\_cff}. For a quick oversight of 1582 the schemes activity the global maximum values of the absolute implicit component 1583 of the vertical velocity and the partitioning coefficient are written to the netCDF 1584 version of the run statistics file (\texttt{run.stat.nc}) if this is active (see 1585 \autoref{sec:MISC_opt} for activation details). 1586 1587 \autoref{fig:ZDF_zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for 1588 the various overflow tests. Note that the adaptive-implicit vertical advection scheme is 1589 active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the 1590 test case is close to stability limits even with this value. At the larger timesteps, the 1591 vertical velocity is treated mostly implicitly at some location throughout the run. The 1592 oscillatory nature of this measure appears to be linked to the progress of the plume front 1593 as each cusp is associated with the location of the maximum shifting to the adjacent cell. 1594 This is illustrated in \autoref{fig:ZDF_zad_Aimp_maxCf_loc} where the i- and k- locations of the 1595 maximum have been overlaid for the base run case. 1596 1597 \medskip 1598 \noindent Only limited tests have been performed in more realistic configurations. In the 1599 ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes 1600 restartability and reproducibility tests but it is unable to improve the model's stability 1601 enough to allow an increase in the model time-step. A view of the time-series of maximum 1602 partitioning coefficient (not shown here) suggests that the default time-step of 5400s is 1603 already pushing at stability limits, especially in the initial start-up phase. The 1604 time-series does not, however, exhibit any of the 'cuspiness' found with the overflow 1605 tests. 1606 1607 \medskip 1608 \noindent A short test with an eORCA1 configuration promises more since a test using a 1609 time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the 1610 time-step is limited to 2700s without. 1611 1612 \begin{figure}[!t] 1613 \centering 1614 \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_overflow_all_rdt} 1615 \caption[OVERFLOW: sample temperature vertical cross-sections from mid- and end-run]{ 1616 Sample temperature vertical cross-sections from mid- and end-run using 1617 different values for \forcode{nn_rdt} and with or without adaptive implicit vertical advection. 1618 Without the adaptive implicit vertical advection 1619 only the run with the shortest timestep is able to run to completion. 1620 Note also that the colour-scale has been chosen to confirm that 1621 temperatures remain within the original range of 10$^o$ to 20$^o$.} 1622 \label{fig:ZDF_zad_Aimp_overflow_all_rdt} 1623 \end{figure} 1624 1625 \begin{figure}[!t] 1626 \centering 1627 \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_maxCf} 1628 \caption[OVERFLOW: maximum partitioning coefficient during a series of test runs]{ 1629 The maximum partitioning coefficient during a series of test runs with 1630 increasing model timestep length. 1631 At the larger timesteps, 1632 the vertical velocity is treated mostly implicitly at some location throughout the run.} 1633 \label{fig:ZDF_zad_Aimp_maxCf} 1634 \end{figure} 1635 1636 \begin{figure}[!t] 1637 \centering 1638 \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_maxCf_loc} 1639 \caption[OVERFLOW: maximum partitioning coefficient for the case overlaid]{ 1640 The maximum partitioning coefficient for the \forcode{nn_rdt=10.0} case overlaid with 1641 information on the gridcell i- and k-locations of the maximum value.} 1642 \label{fig:ZDF_zad_Aimp_maxCf_loc} 1643 \end{figure} 1644 1645 \subinc{\input{../../global/epilogue}} 1404 1646 1405 1647 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.