- Timestamp:
- 2019-10-25T16:27:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc
-
Property
svn:externals
set to
^/utils/badges badges
^/utils/logos logos
-
Property
svn:externals
set to
-
NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex
- Property svn:ignore deleted
-
NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO
-
Property
svn:externals
set to
^/utils/figures/NEMO figures
-
Property
svn:externals
set to
-
NEMO/branches/2019/dev_r11470_HPC_12_mpi3/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_r11470_HPC_12_mpi3/doc/latex/NEMO/subfiles/chap_ZDF.tex
r11435 r11799 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 11 \thispagestyle{plain} 12 10 13 \chaptertoc 11 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 12 30 %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 31 14 \newpage 15 16 % ================================================================ 17 % Vertical Mixing 18 % ================================================================ 32 %% ================================================================================================= 19 33 \section{Vertical mixing} 20 \label{sec:ZDF _zdf}34 \label{sec:ZDF} 21 35 22 36 The discrete form of the ocean subgrid scale physics has been presented in … … 25 39 At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}), 26 40 while at the bottom they are set to zero for heat and salt, 27 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie\ \np{ln \_trabbc} defined,41 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie\ \np{ln_trabbc}{ln\_trabbc} defined, 28 42 see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum 29 43 (see \autoref{sec:ZDF_drg}). … … 39 53 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 40 54 %These trends can be computed using either a forward time stepping scheme 41 %(namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme 42 %(\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing coefficients, 43 %and thus of the formulation used (see \autoref{chap:STP}). 44 45 %--------------------------------------------namzdf-------------------------------------------------------- 46 47 \nlst{namzdf} 48 %-------------------------------------------------------------------------------------------------------------- 49 50 % ------------------------------------------------------------------------------------------------------------- 51 % Constant 52 % ------------------------------------------------------------------------------------------------------------- 53 \subsection[Constant (\forcode{ln_zdfcst = .true.})] 54 {Constant (\protect\np{ln\_zdfcst}\forcode{ = .true.})} 55 %(namelist parameter \np[=.true.]{ln_zdfexp}{ln\_zdfexp}) or a backward time stepping scheme 56 %(\np[=.false.]{ln_zdfexp}{ln\_zdfexp}) depending on the magnitude of the mixing coefficients, 57 %and thus of the formulation used (see \autoref{chap:TD}). 58 59 \begin{listing} 60 \nlst{namzdf} 61 \caption{\forcode{&namzdf}} 62 \label{lst:namzdf} 63 \end{listing} 64 65 %% ================================================================================================= 66 \subsection[Constant (\forcode{ln_zdfcst})]{Constant (\protect\np{ln_zdfcst}{ln\_zdfcst})} 55 67 \label{subsec:ZDF_cst} 56 68 57 Options are defined through the \nam{zdf} namelist variables.58 When \np{ln \_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to69 Options are defined through the \nam{zdf}{zdf} namelist variables. 70 When \np{ln_zdfcst}{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 59 71 constant values over the whole ocean. 60 72 This is the crudest way to define the vertical ocean physics. … … 66 78 \end{align*} 67 79 68 These values are set through the \np{rn \_avm0} and \np{rn\_avt0} namelist parameters.80 These values are set through the \np{rn_avm0}{rn\_avm0} and \np{rn_avt0}{rn\_avt0} namelist parameters. 69 81 In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity, 70 82 that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and 71 83 $\sim10^{-9}~m^2.s^{-1}$ for salinity. 72 84 73 % ------------------------------------------------------------------------------------------------------------- 74 % Richardson Number Dependent 75 % ------------------------------------------------------------------------------------------------------------- 76 \subsection[Richardson number dependent (\forcode{ln_zdfric = .true.})] 77 {Richardson number dependent (\protect\np{ln\_zdfric}\forcode{ = .true.})} 85 %% ================================================================================================= 86 \subsection[Richardson number dependent (\forcode{ln_zdfric})]{Richardson number dependent (\protect\np{ln_zdfric}{ln\_zdfric})} 78 87 \label{subsec:ZDF_ric} 79 88 80 %--------------------------------------------namric--------------------------------------------------------- 81 82 \nlst{namzdf_ric} 83 %-------------------------------------------------------------------------------------------------------------- 84 85 When \np{ln\_zdfric}\forcode{ = .true.}, a local Richardson number dependent formulation for the vertical momentum and 86 tracer eddy coefficients is set through the \nam{zdf\_ric} namelist variables. 89 \begin{listing} 90 \nlst{namzdf_ric} 91 \caption{\forcode{&namzdf_ric}} 92 \label{lst:namzdf_ric} 93 \end{listing} 94 95 When \np[=.true.]{ln_zdfric}{ln\_zdfric}, a local Richardson number dependent formulation for the vertical momentum and 96 tracer eddy coefficients is set through the \nam{zdf_ric}{zdf\_ric} namelist variables. 87 97 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. 88 98 \textit{In situ} measurements have been used to link vertical turbulent activity to large scale ocean structures. … … 92 102 Following \citet{pacanowski.philander_JPO81}, the following formulation has been implemented: 93 103 \[ 94 % \label{eq: zdfric}104 % \label{eq:ZDF_ric} 95 105 \left\{ 96 106 \begin{aligned} … … 105 115 (see \autoref{subsec:ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ is the maximum value that 106 116 can be reached by the coefficient when $Ri\leq 0$, $a=5$ and $n=2$. 107 The last three values can be modified by setting the \np{rn \_avmri}, \np{rn\_alp} and108 \np{nn \_ric} namelist parameters, respectively.117 The last three values can be modified by setting the \np{rn_avmri}{rn\_avmri}, \np{rn_alp}{rn\_alp} and 118 \np{nn_ric}{nn\_ric} namelist parameters, respectively. 109 119 110 120 A simple mixing-layer model to transfer and dissipate the atmospheric forcings 111 (wind-stress and buoyancy fluxes) can be activated setting the \np {ln\_mldw}\forcode{ = .true.} in the namelist.121 (wind-stress and buoyancy fluxes) can be activated setting the \np[=.true.]{ln_mldw}{ln\_mldw} in the namelist. 112 122 113 123 In this case, the local depth of turbulent wind-mixing or "Ekman depth" $h_{e}(x,y,t)$ is evaluated and … … 125 135 \] 126 136 is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$. 127 The final $h_{e}$ is further constrained by the adjustable bounds \np{rn \_mldmin} and \np{rn\_mldmax}.137 The final $h_{e}$ is further constrained by the adjustable bounds \np{rn_mldmin}{rn\_mldmin} and \np{rn_mldmax}{rn\_mldmax}. 128 138 Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to 129 the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{lermusiaux_JMS01}. 130 131 % ------------------------------------------------------------------------------------------------------------- 132 % TKE Turbulent Closure Scheme 133 % ------------------------------------------------------------------------------------------------------------- 134 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke = .true.})] 135 {TKE turbulent closure scheme (\protect\np{ln\_zdftke}\forcode{ = .true.})} 139 the empirical values \np{rn_wtmix}{rn\_wtmix} and \np{rn_wvmix}{rn\_wvmix} \citep{lermusiaux_JMS01}. 140 141 %% ================================================================================================= 142 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke})]{TKE turbulent closure scheme (\protect\np{ln_zdftke}{ln\_zdftke})} 136 143 \label{subsec:ZDF_tke} 137 %--------------------------------------------namzdf_tke-------------------------------------------------- 138 139 \nlst{namzdf_tke} 140 %-------------------------------------------------------------------------------------------------------------- 144 145 \begin{listing} 146 \nlst{namzdf_tke} 147 \caption{\forcode{&namzdf_tke}} 148 \label{lst:namzdf_tke} 149 \end{listing} 141 150 142 151 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on … … 151 160 its destruction through stratification, its vertical diffusion, and its dissipation of \citet{kolmogorov_IANS42} type: 152 161 \begin{equation} 153 \label{eq: zdftke_e}162 \label{eq:ZDF_tke_e} 154 163 \frac{\partial \bar{e}}{\partial t} = 155 164 \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2 … … 161 170 \end{equation} 162 171 \[ 163 % \label{eq: zdftke_kz}172 % \label{eq:ZDF_tke_kz} 164 173 \begin{split} 165 174 K_m &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ … … 172 181 The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with 173 182 vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}. 174 They are set through namelist parameters \np{nn \_ediff} and \np{nn\_ediss}.183 They are set through namelist parameters \np{nn_ediff}{nn\_ediff} and \np{nn_ediss}{nn\_ediss}. 175 184 $P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$: 176 185 \begin{align*} 177 % \label{eq: prt}186 % \label{eq:ZDF_prt} 178 187 P_{rt} = 179 188 \begin{cases} … … 183 192 \end{cases} 184 193 \end{align*} 185 The choice of $P_{rt}$ is controlled by the \np{nn \_pdl} namelist variable.194 The choice of $P_{rt}$ is controlled by the \np{nn_pdl}{nn\_pdl} namelist variable. 186 195 187 196 At the sea surface, the value of $\bar{e}$ is prescribed from the wind stress field as 188 $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn \_ebb} namelist parameter.197 $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter. 189 198 The default value of $e_{bb}$ is 3.75. \citep{gaspar.gregoris.ea_JGR90}), however a much larger value can be used when 190 199 taking into account the surface wave breaking (see below Eq. \autoref{eq:ZDF_Esbc}). … … 192 201 The time integration of the $\bar{e}$ equation may formally lead to negative values because 193 202 the numerical scheme does not ensure its positivity. 194 To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn \_emin} namelist parameter).203 To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn_emin}{rn\_emin} namelist parameter). 195 204 Following \citet{gaspar.gregoris.ea_JGR90}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 196 205 This allows the subsequent formulations to match that of \citet{gargett_JMR84} for the diffusion in … … 198 207 In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical instabilities associated with 199 208 too weak vertical diffusion. 200 They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and 201 \np{rn\_avt0} (\nam{zdf} namelist, see \autoref{subsec:ZDF_cst}). 202 209 They must be specified at least larger than the molecular values, and are set through \np{rn_avm0}{rn\_avm0} and 210 \np{rn_avt0}{rn\_avt0} (\nam{zdf}{zdf} namelist, see \autoref{subsec:ZDF_cst}). 211 212 %% ================================================================================================= 203 213 \subsubsection{Turbulent length scale} 204 214 205 215 For computational efficiency, the original formulation of the turbulent length scales proposed by 206 216 \citet{gaspar.gregoris.ea_JGR90} has been simplified. 207 Four formulations are proposed, the choice of which is controlled by the \np{nn \_mxl} namelist parameter.217 Four formulations are proposed, the choice of which is controlled by the \np{nn_mxl}{nn\_mxl} namelist parameter. 208 218 The first two are based on the following first order approximation \citep{blanke.delecluse_JPO93}: 209 219 \begin{equation} 210 \label{eq: tke_mxl0_1}220 \label{eq:ZDF_tke_mxl0_1} 211 221 l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N 212 222 \end{equation} 213 223 which is valid in a stable stratified region with constant values of the Brunt-Vais\"{a}l\"{a} frequency. 214 224 The resulting length scale is bounded by the distance to the surface or to the bottom 215 (\np {nn\_mxl}\forcode{ = 0}) or by the local vertical scale factor (\np{nn\_mxl}\forcode{ = 1}).225 (\np[=0]{nn_mxl}{nn\_mxl}) or by the local vertical scale factor (\np[=1]{nn_mxl}{nn\_mxl}). 216 226 \citet{blanke.delecluse_JPO93} notice that this simplification has two major drawbacks: 217 227 it makes no sense for locally unstable stratification and the computation no longer uses all 218 228 the information contained in the vertical density profile. 219 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np {nn\_mxl}\forcode{ = 2, 3} cases,229 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np[=2, 3]{nn_mxl}{nn\_mxl} cases, 220 230 which add an extra assumption concerning the vertical gradient of the computed length scale. 221 So, the length scales are first evaluated as in \autoref{eq: tke_mxl0_1} and then bounded such that:231 So, the length scales are first evaluated as in \autoref{eq:ZDF_tke_mxl0_1} and then bounded such that: 222 232 \begin{equation} 223 \label{eq: tke_mxl_constraint}233 \label{eq:ZDF_tke_mxl_constraint} 224 234 \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1 225 235 \qquad \text{with }\ l = l_k = l_\epsilon 226 236 \end{equation} 227 \autoref{eq: tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than237 \autoref{eq:ZDF_tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than 228 238 the variations of depth. 229 239 It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less … … 231 241 In particular, it allows the length scale to be limited not only by the distance to the surface or 232 242 to the ocean bottom but also by the distance to a strongly stratified portion of the water column such as 233 the thermocline (\autoref{fig: mixing_length}).234 In order to impose the \autoref{eq: tke_mxl_constraint} constraint, we introduce two additional length scales:243 the thermocline (\autoref{fig:ZDF_mixing_length}). 244 In order to impose the \autoref{eq:ZDF_tke_mxl_constraint} constraint, we introduce two additional length scales: 235 245 $l_{up}$ and $l_{dwn}$, the upward and downward length scales, and 236 246 evaluate the dissipation and mixing length scales as 237 247 (and note that here we use numerical indexing): 238 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>239 248 \begin{figure}[!t] 240 \begin{center} 241 \includegraphics[width=\textwidth]{Fig_mixing_length} 242 \caption{ 243 \protect\label{fig:mixing_length} 244 Illustration of the mixing length computation. 245 } 246 \end{center} 249 \centering 250 \includegraphics[width=0.66\textwidth]{Fig_mixing_length} 251 \caption[Mixing length computation]{Illustration of the mixing length computation} 252 \label{fig:ZDF_mixing_length} 247 253 \end{figure} 248 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 249 \[ 250 % \label{eq:tke_mxl2} 254 \[ 255 % \label{eq:ZDF_tke_mxl2} 251 256 \begin{aligned} 252 257 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) … … 256 261 \end{aligned} 257 262 \] 258 where $l^{(k)}$ is computed using \autoref{eq: tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.259 260 In the \np {nn\_mxl}\forcode{ = 2} case, the dissipation and mixing length scales take the same value:261 $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np {nn\_mxl}\forcode{ = 3} case,263 where $l^{(k)}$ is computed using \autoref{eq:ZDF_tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. 264 265 In the \np[=2]{nn_mxl}{nn\_mxl} case, the dissipation and mixing length scales take the same value: 266 $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np[=3]{nn_mxl}{nn\_mxl} case, 262 267 the dissipation and mixing turbulent length scales are give as in \citet{gaspar.gregoris.ea_JGR90}: 263 268 \[ 264 % \label{eq: tke_mxl_gaspar}269 % \label{eq:ZDF_tke_mxl_gaspar} 265 270 \begin{aligned} 266 271 & l_k = \sqrt{\ l_{up} \ \ l_{dwn}\ } \\ … … 269 274 \] 270 275 271 At the ocean surface, a non zero length scale is set through the \np{rn \_mxl0} namelist parameter.276 At the ocean surface, a non zero length scale is set through the \np{rn_mxl0}{rn\_mxl0} namelist parameter. 272 277 Usually the surface scale is given by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant and 273 278 $z_o$ the roughness parameter of the surface. 274 Assuming $z_o=0.1$~m \citep{craig.banner_JPO94} leads to a 0.04~m, the default value of \np{rn \_mxl0}.279 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}. 275 280 In the ocean interior a minimum length scale is set to recover the molecular viscosity when 276 281 $\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ). 277 282 283 %% ================================================================================================= 278 284 \subsubsection{Surface wave breaking parameterization} 279 %-----------------------------------------------------------------------%280 285 281 286 Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to … … 303 308 $\alpha_{CB} = 100$ the Craig and Banner's value. 304 309 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 305 with $e_{bb}$ the \np{rn \_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds310 with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter, setting \np[=67.83]{rn_ebb}{rn\_ebb} corresponds 306 311 to $\alpha_{CB} = 100$. 307 Further setting \np {ln\_mxl0}\forcode{ =.true.}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,312 Further setting \np[=.true.]{ln_mxl0}{ln\_mxl0}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale, 308 313 with $\beta$ hard coded to the Stacey's value. 309 Note that a minimal threshold of \np{rn \_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the314 Note that a minimal threshold of \np{rn_emin0}{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the 310 315 surface $\bar{e}$ value. 311 316 317 %% ================================================================================================= 312 318 \subsubsection{Langmuir cells} 313 %--------------------------------------%314 319 315 320 Langmuir circulations (LC) can be described as ordered large-scale vertical motions in … … 325 330 The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in 326 331 an extra source term of TKE, $P_{LC}$. 327 The presence of $P_{LC}$ in \autoref{eq: zdftke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to328 \forcode{.true.} in the \nam{zdf \_tke} namelist.332 The presence of $P_{LC}$ in \autoref{eq:ZDF_tke_e}, the TKE equation, is controlled by setting \np{ln_lc}{ln\_lc} to 333 \forcode{.true.} in the \nam{zdf_tke}{zdf\_tke} namelist. 329 334 330 335 By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}), … … 354 359 where $c_{LC} = 0.15$ has been chosen by \citep{axell_JGR02} as a good compromise to fit LES data. 355 360 The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second. 356 The value of $c_{LC}$ is set through the \np{rn \_lc} namelist parameter,361 The value of $c_{LC}$ is set through the \np{rn_lc}{rn\_lc} namelist parameter, 357 362 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}. 358 363 … … 364 369 \] 365 370 371 %% ================================================================================================= 366 372 \subsubsection{Mixing just below the mixed layer} 367 %--------------------------------------------------------------%368 373 369 374 Vertical mixing parameterizations commonly used in ocean general circulation models tend to … … 376 381 (\ie\ near-inertial oscillations and ocean swells and waves). 377 382 378 When using this parameterization (\ie\ when \np {nn\_etau}\forcode{ = 1}),383 When using this parameterization (\ie\ when \np[=1]{nn_etau}{nn\_etau}), 379 384 the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations, 380 385 swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition, … … 388 393 the penetration, and $f_i$ is the ice concentration 389 394 (no penetration if $f_i=1$, \ie\ if the ocean is entirely covered by sea-ice). 390 The value of $f_r$, usually a few percents, is specified through \np{rn \_efr} namelist parameter.391 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np {nn\_etau}\forcode{ = 0}) or395 The value of $f_r$, usually a few percents, is specified through \np{rn_efr}{rn\_efr} namelist parameter. 396 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np[=0]{nn_etau}{nn\_etau}) or 392 397 a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m at high latitudes 393 (\np {nn\_etau}\forcode{ = 1}).394 395 Note that two other option exist, \np {nn\_etau}\forcode{ = 2, 3}.398 (\np[=1]{nn_etau}{nn\_etau}). 399 400 Note that two other option exist, \np[=2, 3]{nn_etau}{nn\_etau}. 396 401 They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer, 397 402 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean. … … 400 405 401 406 % This should be explain better below what this rn_eice parameter is meant for: 402 In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn \_eice} namelist parameter.407 In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn_eice}{rn\_eice} namelist parameter. 403 408 This parameter varies from \forcode{0} for no effect to \forcode{4} to suppress the TKE input into the ocean when Sea Ice concentration 404 409 is greater than 25\%. … … 412 417 % (\eg\ Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 413 418 414 % ------------------------------------------------------------------------------------------------------------- 415 % GLS Generic Length Scale Scheme 416 % ------------------------------------------------------------------------------------------------------------- 417 \subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls = .true.})] 418 {GLS: Generic Length Scale (\protect\np{ln\_zdfgls}\forcode{ = .true.})} 419 %% ================================================================================================= 420 \subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls})]{GLS: Generic Length Scale (\protect\np{ln_zdfgls}{ln\_zdfgls})} 419 421 \label{subsec:ZDF_gls} 420 422 421 %--------------------------------------------namzdf_gls--------------------------------------------------------- 422 423 \nlst{namzdf_gls} 424 %-------------------------------------------------------------------------------------------------------------- 423 \begin{listing} 424 \nlst{namzdf_gls} 425 \caption{\forcode{&namzdf_gls}} 426 \label{lst:namzdf_gls} 427 \end{listing} 425 428 426 429 The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations: … … 428 431 $\psi$ \citep{umlauf.burchard_JMR03, umlauf.burchard_CSR05}. 429 432 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 430 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab: GLS} allows to recover a number of433 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:ZDF_GLS} allows to recover a number of 431 434 well-known turbulent closures ($k$-$kl$ \citep{mellor.yamada_RG82}, $k$-$\epsilon$ \citep{rodi_JGR87}, 432 435 $k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}). 433 436 The GLS scheme is given by the following set of equations: 434 437 \begin{equation} 435 \label{eq: zdfgls_e}438 \label{eq:ZDF_gls_e} 436 439 \frac{\partial \bar{e}}{\partial t} = 437 440 \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 … … 443 446 444 447 \[ 445 % \label{eq: zdfgls_psi}448 % \label{eq:ZDF_gls_psi} 446 449 \begin{split} 447 450 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ … … 455 458 456 459 \[ 457 % \label{eq: zdfgls_kz}460 % \label{eq:ZDF_gls_kz} 458 461 \begin{split} 459 462 K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ … … 463 466 464 467 \[ 465 % \label{eq: zdfgls_eps}468 % \label{eq:ZDF_gls_eps} 466 469 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; 467 470 \] … … 470 473 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of 471 474 the choice of the turbulence model. 472 Four different turbulent models are pre-defined (\autoref{tab:GLS}). 473 They are made available through the \np{nn\_clo} namelist parameter. 474 475 %--------------------------------------------------TABLE-------------------------------------------------- 475 Four different turbulent models are pre-defined (\autoref{tab:ZDF_GLS}). 476 They are made available through the \np{nn_clo}{nn\_clo} namelist parameter. 477 476 478 \begin{table}[htbp] 477 \begin{center} 478 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} 479 \begin{tabular}{ccccc} 480 & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ 481 % & \citep{mellor.yamada_RG82} & \citep{rodi_JGR87} & \citep{wilcox_AJ88} & \\ 482 \hline 483 \hline 484 \np{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ 485 \hline 486 $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ 487 $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ 488 $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ 489 $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ 490 $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ 491 $C_3$ & 1. & 1. & 1. & 1. \\ 492 $F_{wall}$ & Yes & -- & -- & -- \\ 493 \hline 494 \hline 495 \end{tabular} 496 \caption{ 497 \protect\label{tab:GLS} 498 Set of predefined GLS parameters, or equivalently predefined turbulence models available with 499 \protect\np{ln\_zdfgls}\forcode{ = .true.} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\nam{zdf\_gls}. 500 } 501 \end{center} 479 \centering 480 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c} 481 \begin{tabular}{ccccc} 482 & $k-kl$ & $k-\epsilon$ & $k-\omega$ & generic \\ 483 % & \citep{mellor.yamada_RG82} & \citep{rodi_JGR87} & \citep{wilcox_AJ88} & \\ 484 \hline 485 \hline 486 \np{nn_clo}{nn\_clo} & \textbf{0} & \textbf{1} & \textbf{2} & \textbf{3} \\ 487 \hline 488 $( p , n , m )$ & ( 0 , 1 , 1 ) & ( 3 , 1.5 , -1 ) & ( -1 , 0.5 , -1 ) & ( 2 , 1 , -0.67 ) \\ 489 $\sigma_k$ & 2.44 & 1. & 2. & 0.8 \\ 490 $\sigma_\psi$ & 2.44 & 1.3 & 2. & 1.07 \\ 491 $C_1$ & 0.9 & 1.44 & 0.555 & 1. \\ 492 $C_2$ & 0.5 & 1.92 & 0.833 & 1.22 \\ 493 $C_3$ & 1. & 1. & 1. & 1. \\ 494 $F_{wall}$ & Yes & -- & -- & -- \\ 495 \hline 496 \hline 497 \end{tabular} 498 \caption[Set of predefined GLS parameters or equivalently predefined turbulence models available]{ 499 Set of predefined GLS parameters, or equivalently predefined turbulence models available with 500 \protect\np[=.true.]{ln_zdfgls}{ln\_zdfgls} and controlled by 501 the \protect\np{nn_clos}{nn\_clos} namelist variable in \protect\nam{zdf_gls}{zdf\_gls}.} 502 \label{tab:ZDF_GLS} 502 503 \end{table} 503 %--------------------------------------------------------------------------------------------------------------504 504 505 505 In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of … … 508 508 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88}, 509 509 or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01} 510 (\np {nn\_stab\_func}\forcode{ = 0, 3}, resp.).510 (\np[=0, 3]{nn_stab_func}{nn\_stab\_func}, resp.). 511 511 The value of $C_{0\mu}$ depends on the choice of the stability function. 512 512 513 513 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or 514 Neumann condition through \np{nn \_bc\_surf} and \np{nn\_bc\_bot}, resp.514 Neumann condition through \np{nn_bc_surf}{nn\_bc\_surf} and \np{nn_bc_bot}{nn\_bc\_bot}, resp. 515 515 As for TKE closure, the wave effect on the mixing is considered when 516 \np {rn\_crban}\forcode{ > 0.} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}.517 The \np{rn \_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and518 \np{rn \_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}.516 \np[ > 0.]{rn_crban}{rn\_crban} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}. 517 The \np{rn_crban}{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and 518 \np{rn_charn}{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}. 519 519 520 520 The $\psi$ equation is known to fail in stably stratified flows, and for this reason … … 525 525 the entrainment depth predicted in stably stratified situations, 526 526 and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes. 527 The clipping is only activated if \np {ln\_length\_lim}\forcode{ = .true.},528 and the $c_{lim}$ is set to the \np{rn \_clim\_galp} value.527 The clipping is only activated if \np[=.true.]{ln_length_lim}{ln\_length\_lim}, 528 and the $c_{lim}$ is set to the \np{rn_clim_galp}{rn\_clim\_galp} value. 529 529 530 530 The time and space discretization of the GLS equations follows the same energetic consideration as for … … 533 533 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO\ model. 534 534 535 536 535 % ------------------------------------------------------------------------------------------------------------- 537 % OSM OSMOSIS BL Scheme 536 % OSM OSMOSIS BL Scheme 538 537 % ------------------------------------------------------------------------------------------------------------- 539 \subsection[OSM: OSM osisboundary layer scheme (\forcode{ln_zdfosm = .true.})]540 {OSM: OSM osis boundary layer scheme (\protect\np{ln\_zdfosm}\forcode{ = .true.})}538 \subsection[OSM: OSMOSIS boundary layer scheme (\forcode{ln_zdfosm = .true.})] 539 {OSM: OSMOSIS boundary layer scheme (\protect\np{ln_zdfosm}{ln\_zdfosm})} 541 540 \label{subsec:ZDF_osm} 542 %--------------------------------------------namzdf_osm--------------------------------------------------------- 543 544 \nlst{namzdf_osm} 541 542 \begin{listing} 543 \nlst{namzdf_osm} 544 \caption{\forcode{&namzdf_osm}} 545 \label{lst:namzdf_osm} 546 \end{listing} 547 545 548 %-------------------------------------------------------------------------------------------------------------- 546 547 The OSMOSIS turbulent closure scheme is based on...... TBC 548 549 % ------------------------------------------------------------------------------------------------------------- 550 % TKE and GLS discretization considerations 551 % ------------------------------------------------------------------------------------------------------------- 552 \subsection[ Discrete energy conservation for TKE and GLS schemes] 553 {Discrete energy conservation for TKE and GLS schemes} 554 \label{subsec:ZDF_tke_ene} 555 556 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 549 \paragraph{Namelist choices} 550 Most of the namelist options refer to how to specify the Stokes 551 surface drift and penetration depth. There are three options: 552 \begin{description} 553 \item \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave} Default value in \texttt{namelist\_ref}. In this case the Stokes drift is 554 assumed to be parallel to the surface wind stress, with 555 magnitude consistent with a constant turbulent Langmuir number 556 $\mathrm{La}_t=$ \protect\np{rn_m_la} {rn\_m\_la} i.e.\ 557 $u_{s0}=\tau/(\mathrm{La}_t^2\rho_0)$. Default value of 558 \protect\np{rn_m_la}{rn\_m\_la} is 0.3. The Stokes penetration 559 depth $\delta = $ \protect\np{rn_osm_dstokes}{rn\_osm\_dstokes}; this has default value 560 of 5~m. 561 562 \item \protect\np[=1]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is 563 assumed to be parallel to the surface wind stress, with 564 magnitude as in the classical Pierson-Moskowitz wind-sea 565 spectrum. Significant wave height and 566 wave-mean period taken from this spectrum are used to calculate the Stokes penetration 567 depth, following the approach set out in \citet{breivik.janssen.ea_JPO14}. 568 569 \item \protect\np[=2]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is 570 taken from ECMWF wave model output, though only the component parallel 571 to the wind stress is retained. Significant wave height and 572 wave-mean period from ECMWF wave model output are used to calculate the Stokes penetration 573 depth, again following \citet{breivik.janssen.ea_JPO14}. 574 575 \end{description} 576 577 Others refer to the treatment of diffusion and viscosity beneath 578 the surface boundary layer: 579 \begin{description} 580 \item \protect\np{ln_kpprimix} {ln\_kpprimix} Default is \np{.true.}. Switches on KPP-style Ri \#-dependent 581 mixing below the surface boundary layer. If this is set 582 \texttt{.true.} the following variable settings are honoured: 583 \item \protect\np{rn_riinfty}{rn\_riinfty} Critical value of local Ri \# below which 584 shear instability increases vertical mixing from background value. 585 \item \protect\np{rn_difri} {rn\_difri} Maximum value of Ri \#-dependent mixing at $\mathrm{Ri}=0$. 586 \item \protect\np{ln_convmix}{ln\_convmix} If \texttt{.true.} then, where water column is unstable, specify 587 diffusivity equal to \protect\np{rn_dif_conv}{rn\_dif\_conv} (default value is 1 m~s$^{-2}$). 588 \end{description} 589 Diagnostic output is controlled by: 590 \begin{description} 591 \item \protect\np{ln_dia_osm}{ln\_dia\_osm} Default is \np{.false.}; allows XIOS output of OSMOSIS internal fields. 592 \end{description} 593 Obsolete namelist parameters include: 594 \begin{description} 595 \item \protect\np{ln_use_osm_la}\np{ln\_use\_osm\_la} With \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave}, 596 \protect\np{rn_osm_dstokes} {rn\_osm\_dstokes} is always used to specify the Stokes 597 penetration depth. 598 \item \protect\np{nn_ave} {nn\_ave} Choice of averaging method for KPP-style Ri \# 599 mixing. Not taken account of. 600 \item \protect\np{rn_osm_hbl0} {rn\_osm\_hbl0} Depth of initial boundary layer is now set 601 by a density criterion similar to that used in calculating \emph{hmlp} (output as \texttt{mldr10\_1}) in \mdl{zdfmxl}. 602 \end{description} 603 604 \subsubsection{Summary} 605 Much of the time the turbulent motions in the ocean surface boundary 606 layer (OSBL) are not given by 607 classical shear turbulence. Instead they are in a regime known as 608 `Langmuir turbulence', dominated by an 609 interaction between the currents and the Stokes drift of the surface waves \citep[e.g.][]{mcwilliams.ea_JFM97}. 610 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$. 611 612 The OSMOSIS model is fundamentally based on results of Large Eddy 613 Simulations (LES) of Langmuir turbulence and aims to fully describe 614 this Langmuir regime. The description in this section is of necessity incomplete and further details are available in Grant. A (2019); in prep. 615 616 The OSMOSIS turbulent closure scheme is a similarity-scale scheme in 617 the same spirit as the K-profile 618 parameterization (KPP) scheme of \citet{large.ea_RG97}. 619 A specified shape of diffusivity, scaled by the (OSBL) depth 620 $h_{\mathrm{BL}}$ and a turbulent velocity scale, is imposed throughout the 621 boundary layer 622 $-h_{\mathrm{BL}}<z<\eta$. The turbulent closure model 623 also includes fluxes of tracers and momentum that are``non-local'' (independent of the local property gradient). 624 625 Rather than the OSBL 626 depth being diagnosed in terms of a bulk Richardson number criterion, 627 as in KPP, it is set by a prognostic equation that is informed by 628 energy budget considerations reminiscent of the classical mixed layer 629 models of \citet{kraus.turner_tellus67}. 630 The model also includes an explicit parametrization of the structure 631 of the pycnocline (the stratified region at the bottom of the OSBL). 632 633 Presently, mixing below the OSBL is handled by the Richardson 634 number-dependent mixing scheme used in \citet{large.ea_RG97}. 635 636 Convective parameterizations such as described in \ref{sec:ZDF_conv} 637 below should not be used with the OSMOSIS-OBL model: instabilities 638 within the OSBL are part of the model, while instabilities below the 639 ML are handled by the Ri \# dependent scheme. 640 641 \subsubsection{Depth and velocity scales} 642 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). 557 643 \begin{figure}[!t] 558 644 \begin{center} 559 \includegraphics[width= \textwidth]{Fig_ZDF_TKE_time_scheme}645 \includegraphics[width=0.7\textwidth]{Fig_ZDF_OSM_structure_of_OSBL} 560 646 \caption{ 561 \protect\label{fig: TKE_time_scheme}562 Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and its links to the momentum and tracer time integration.647 \protect\label{fig: OSBL_structure} 648 The structure of the entraining boundary layer. (a) Mean buoyancy profile. (b) Profile of the buoyancy flux. 563 649 } 564 650 \end{center} 565 651 \end{figure} 566 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 652 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. 653 654 Consideration of the power input by wind acting on the Stokes drift suggests that the Langmuir turbulence has velocity scale: 655 \begin{equation}\label{eq:w_La} 656 w_{*L}= \left(u_*^2 u_{s\,0}\right)^{1/3}; 657 \end{equation} 658 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: 659 \begin{equation}\label{eq:composite-nu} 660 \nu_{\ast}= \left\{ u_*^3 \left[1-\exp(-.5 \mathrm{La}_t^2)\right]+w_{*L}^3\right\}^{1/3}. 661 \end{equation} 662 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: 663 \begin{equation}\label{eq:vel-scale-unstable} 664 \omega_* = \left(\nu_*^3 + 0.5 w_{*C}^3\right)^{1/3}. 665 \end{equation} 666 667 \subsubsection{The flux gradient model} 668 The flux-gradient relationships used in the OSMOSIS scheme take the form: 669 % 670 \begin{equation}\label{eq:flux-grad-gen} 671 \overline{w^\prime\chi^\prime}=-K\frac{\partial\overline{\chi}}{\partial z} + N_{\chi,s} +N_{\chi,b} +N_{\chi,t}, 672 \end{equation} 673 % 674 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. 675 676 In terms of the non-dimensionalized depth variables 677 % 678 \begin{equation}\label{eq:sigma} 679 \sigma_{\mathrm{ml}}= -z/h_{\mathrm{ml}}; \;\sigma_{\mathrm{bl}}= -z/h_{\mathrm{bl}}, 680 \end{equation} 681 % 682 in unstable conditions the eddy diffusivity ($K_d$) and eddy viscosity ($K_\nu$) profiles are parametrized as: 683 % 684 \begin{align}\label{eq:diff-unstable} 685 K_d=&0.8\, \omega_*\, h_{\mathrm{ml}} \, \sigma_{\mathrm{ml}} \left(1-\beta_d \sigma_{\mathrm{ml}}\right)^{3/2} 686 \\\label{eq:visc-unstable} 687 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) 688 \end{align} 689 % 690 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 691 % 692 \begin{equation}\label{eq:diff-wml-base} 693 K_{d,\mathrm{ml}}=K_{\nu,\mathrm{ml}}=\,0.16\,\omega_* \Delta h. 694 \end{equation} 695 % 696 For stable conditions the eddy diffusivity/viscosity profiles are given by: 697 % 698 \begin{align}\label{diff-stable} 699 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} 700 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). 701 \end{align} 702 % 703 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: 704 \begin{equation}\label{eq:L_L} 705 L_L=-w_{*L}^3/\left<\overline{w^\prime b^\prime}\right>_L, 706 \end{equation} 707 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 708 \begin{equation} \label{eq:stable-av-buoy-flux} 709 \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]. 710 \end{equation} 711 % 712 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}$. 713 714 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). 715 716 \subsubsection{Evolution of the boundary layer depth} 717 718 The prognostic equation for the depth of the neutral/unstable boundary layer is given by \citep{grant+etal18}, 719 720 \begin{equation} \label{eq:dhdt-unstable} 721 %\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}} 722 \frac{\partial h_\mathrm{bl}}{\partial t} = W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}} 723 \end{equation} 724 where $h_\mathrm{bl}$ is the horizontally-varying depth of the OSBL, 725 $\mathbf{U}_b$ and $W_b$ are the mean horizontal and vertical 726 velocities at the base of the OSBL, ${\overline{w^\prime 727 b^\prime}}_\mathrm{ent}$ is the buoyancy flux due to entrainment 728 and $\Delta B_\mathrm{bl}$ is the difference between the buoyancy 729 averaged over the depth of the OSBL (i.e.\ including the ML and 730 pycnocline) and the buoyancy just below the base of the OSBL. This 731 equation for the case when the pycnocline has a finite thickness, 732 based on the potential energy budget of the OSBL, is the leading term 733 \citep{grant+etal18} of a generalization of that used in mixed-layer 734 models e.g.\ \citet{kraus.turner_tellus67}, in which the thickness of the pycnocline is taken to be zero. 735 736 The entrainment flux for the combination of convective and Langmuir turbulence is given by 737 \begin{equation} \label{eq:entrain-flux} 738 {\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}}} 739 + 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] 740 \end{equation} 741 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}$. 742 743 For the stable boundary layer, the equation for the depth of the OSBL is: 744 745 \begin{equation}\label{eq:dhdt-stable} 746 \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. 747 \end{equation} 748 749 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. 750 751 752 %% ================================================================================================= 753 \subsection[ Discrete energy conservation for TKE and GLS schemes]{Discrete energy conservation for TKE and GLS schemes} 754 \label{subsec:ZDF_tke_ene} 755 756 \begin{figure}[!t] 757 \centering 758 \includegraphics[width=0.66\textwidth]{Fig_ZDF_TKE_time_scheme} 759 \caption[Subgrid kinetic energy integration in GLS and TKE schemes]{ 760 Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and 761 its links to the momentum and tracer time integration.} 762 \label{fig:ZDF_TKE_time_scheme} 763 \end{figure} 567 764 568 765 The production of turbulence by vertical shear (the first term of the right hand side of 569 \autoref{eq: zdftke_e}) and \autoref{eq:zdfgls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion570 (first line in \autoref{eq: PE_zdf}).766 \autoref{eq:ZDF_tke_e}) and \autoref{eq:ZDF_gls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion 767 (first line in \autoref{eq:MB_zdf}). 571 768 To do so a special care has to be taken for both the time and space discretization of 572 769 the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}. 573 770 574 Let us first address the time stepping issue. \autoref{fig: TKE_time_scheme} shows how771 Let us first address the time stepping issue. \autoref{fig:ZDF_TKE_time_scheme} shows how 575 772 the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with 576 773 the one-level forward time stepping of the equation for $\bar{e}$. … … 579 776 summing the result vertically: 580 777 \begin{equation} 581 \label{eq: energ1}778 \label{eq:ZDF_energ1} 582 779 \begin{split} 583 780 \int_{-H}^{\eta} u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt} \right) \,dz \\ … … 587 784 \end{equation} 588 785 Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$, 589 known at time $t$ (\autoref{fig: TKE_time_scheme}), as it is required when using the TKE scheme590 (see \autoref{sec: STP_forward_imp}).591 The first term of the right hand side of \autoref{eq: energ1} represents the kinetic energy transfer at786 known at time $t$ (\autoref{fig:ZDF_TKE_time_scheme}), as it is required when using the TKE scheme 787 (see \autoref{sec:TD_forward_imp}). 788 The first term of the right hand side of \autoref{eq:ZDF_energ1} represents the kinetic energy transfer at 592 789 the surface (atmospheric forcing) and at the bottom (friction effect). 593 790 The second term is always negative. 594 791 It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$. 595 \autoref{eq: energ1} implies that, to be energetically consistent,792 \autoref{eq:ZDF_energ1} implies that, to be energetically consistent, 596 793 the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 597 794 ${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ … … 599 796 600 797 A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 601 (second term of the right hand side of \autoref{eq: zdftke_e} and \autoref{eq:zdfgls_e}).798 (second term of the right hand side of \autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}). 602 799 This term must balance the input of potential energy resulting from vertical mixing. 603 800 The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by 604 801 multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically: 605 802 \begin{equation} 606 \label{eq: energ2}803 \label{eq:ZDF_energ2} 607 804 \begin{split} 608 805 \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt} \right) \,dz \\ … … 614 811 \end{equation} 615 812 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 616 The first term of the right hand side of \autoref{eq: energ2} is always zero because813 The first term of the right hand side of \autoref{eq:ZDF_energ2} is always zero because 617 814 there is no diffusive flux through the ocean surface and bottom). 618 815 The second term is minus the destruction rate of $\bar{e}$ due to stratification. 619 Therefore \autoref{eq: energ1} implies that, to be energetically consistent,620 the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq: zdftke_e} and \autoref{eq:zdfgls_e}.816 Therefore \autoref{eq:ZDF_energ1} implies that, to be energetically consistent, 817 the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}. 621 818 622 819 Let us now address the space discretization issue. 623 820 The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in 624 the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig: cell}).821 the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:DOM_cell}). 625 822 A space averaging is thus required to obtain the shear TKE production term. 626 By redoing the \autoref{eq: energ1} in the 3D case, it can be shown that the product of eddy coefficient by823 By redoing the \autoref{eq:ZDF_energ1} in the 3D case, it can be shown that the product of eddy coefficient by 627 824 the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. 628 825 Furthermore, the time variation of $e_3$ has be taken into account. … … 630 827 The above energetic considerations leads to the following final discrete form for the TKE equation: 631 828 \begin{equation} 632 \label{eq: zdftke_ene}829 \label{eq:ZDF_tke_ene} 633 830 \begin{split} 634 831 \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt} \equiv … … 647 844 \end{split} 648 845 \end{equation} 649 where the last two terms in \autoref{eq: zdftke_ene} (vertical diffusion and Kolmogorov dissipation)650 are time stepped using a backward scheme (see\autoref{sec: STP_forward_imp}).846 where the last two terms in \autoref{eq:ZDF_tke_ene} (vertical diffusion and Kolmogorov dissipation) 847 are time stepped using a backward scheme (see\autoref{sec:TD_forward_imp}). 651 848 Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible. 652 849 %The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as 653 %they all appear in the right hand side of \autoref{eq: zdftke_ene}.850 %they all appear in the right hand side of \autoref{eq:ZDF_tke_ene}. 654 851 %For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 655 852 656 % ================================================================ 657 % Convection 658 % ================================================================ 853 %% ================================================================================================= 659 854 \section{Convection} 660 855 \label{sec:ZDF_conv} … … 667 862 or/and the use of a turbulent closure scheme. 668 863 669 % ------------------------------------------------------------------------------------------------------------- 670 % Non-Penetrative Convective Adjustment 671 % ------------------------------------------------------------------------------------------------------------- 672 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc = .true.})] 673 {Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})} 864 %% ================================================================================================= 865 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc})]{Non-penetrative convective adjustment (\protect\np{ln_tranpc}{ln\_tranpc})} 674 866 \label{subsec:ZDF_npc} 675 867 676 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>677 868 \begin{figure}[!htb] 678 \begin{center} 679 \includegraphics[width=\textwidth]{Fig_npc} 680 \caption{ 681 \protect\label{fig:npc} 682 Example of an unstable density profile treated by the non penetrative convective adjustment algorithm. 683 $1^{st}$ step: the initial profile is checked from the surface to the bottom. 684 It is found to be unstable between levels 3 and 4. 685 They are mixed. 686 The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed. 687 The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed. 688 The $1^{st}$ step ends since the density profile is then stable below the level 3. 689 $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step: 690 levels 2 to 5 are mixed. 691 The new density profile is checked. 692 It is found stable: end of algorithm. 693 } 694 \end{center} 869 \centering 870 \includegraphics[width=0.66\textwidth]{Fig_npc} 871 \caption[Unstable density profile treated by the non penetrative convective adjustment algorithm]{ 872 Example of an unstable density profile treated by 873 the non penetrative convective adjustment algorithm. 874 $1^{st}$ step: the initial profile is checked from the surface to the bottom. 875 It is found to be unstable between levels 3 and 4. 876 They are mixed. 877 The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed. 878 The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed. 879 The $1^{st}$ step ends since the density profile is then stable below the level 3. 880 $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step: 881 levels 2 to 5 are mixed. 882 The new density profile is checked. 883 It is found stable: end of algorithm.} 884 \label{fig:ZDF_npc} 695 885 \end{figure} 696 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 697 698 Options are defined through the \nam{zdf} namelist variables. 699 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}\forcode{ = .true.}. 700 It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of 886 887 Options are defined through the \nam{zdf}{zdf} namelist variables. 888 The non-penetrative convective adjustment is used when \np[=.true.]{ln_zdfnpc}{ln\_zdfnpc}. 889 It is applied at each \np{nn_npc}{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of 701 890 the water column, but only until the density structure becomes neutrally stable 702 891 (\ie\ until the mixed portion of the water column has \textit{exactly} the density of the water just below) 703 892 \citep{madec.delecluse.ea_JPO91}. 704 The associated algorithm is an iterative process used in the following way (\autoref{fig: npc}):893 The associated algorithm is an iterative process used in the following way (\autoref{fig:ZDF_npc}): 705 894 starting from the top of the ocean, the first instability is found. 706 895 Assume in the following that the instability is located between levels $k$ and $k+1$. … … 734 923 having to recompute the expansion coefficients at each mixing iteration. 735 924 736 % ------------------------------------------------------------------------------------------------------------- 737 % Enhanced Vertical Diffusion 738 % ------------------------------------------------------------------------------------------------------------- 739 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd = .true.})] 740 {Enhanced vertical diffusion (\protect\np{ln\_zdfevd}\forcode{ = .true.})} 925 %% ================================================================================================= 926 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd})]{Enhanced vertical diffusion (\protect\np{ln_zdfevd}{ln\_zdfevd})} 741 927 \label{subsec:ZDF_evd} 742 928 743 Options are defined through the \nam{zdf} namelist variables.744 The enhanced vertical diffusion parameterisation is used when \np {ln\_zdfevd}\forcode{ = .true.}.929 Options are defined through the \nam{zdf}{zdf} namelist variables. 930 The enhanced vertical diffusion parameterisation is used when \np[=.true.]{ln_zdfevd}{ln\_zdfevd}. 745 931 In this case, the vertical eddy mixing coefficients are assigned very large values 746 932 in regions where the stratification is unstable 747 933 (\ie\ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}. 748 This is done either on tracers only (\np {nn\_evdm}\forcode{ = 0}) or749 on both momentum and tracers (\np {nn\_evdm}\forcode{ = 1}).750 751 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np {nn\_evdm}\forcode{ = 1},934 This is done either on tracers only (\np[=0]{nn_evdm}{nn\_evdm}) or 935 on both momentum and tracers (\np[=1]{nn_evdm}{nn\_evdm}). 936 937 In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np[=1]{nn_evdm}{nn\_evdm}, 752 938 the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ values also, are set equal to 753 the namelist parameter \np{rn \_avevd}.939 the namelist parameter \np{rn_avevd}{rn\_avevd}. 754 940 A typical value for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$. 755 941 This parameterisation of convective processes is less time consuming than … … 759 945 Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$. 760 946 This removes a potential source of divergence of odd and even time step in 761 a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:STP_mLF}). 762 763 % ------------------------------------------------------------------------------------------------------------- 764 % Turbulent Closure Scheme 765 % ------------------------------------------------------------------------------------------------------------- 766 \subsection{Handling convection with turbulent closure schemes (\forcode{ln_zdf\{tke,gls,osm\} = .true.})} 947 a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:TD_mLF}). 948 949 %% ================================================================================================= 950 \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}})} 767 951 \label{subsec:ZDF_tcs} 768 952 769 770 953 The turbulent closure schemes presented in \autoref{subsec:ZDF_tke}, \autoref{subsec:ZDF_gls} and 771 \autoref{subsec:ZDF_osm} (\ie\ \np{ln \_zdftke} or \np{ln\_zdfgls} or \np{ln\_zdfosm} defined) deal, in theory,954 \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, 772 955 with statically unstable density profiles. 773 956 In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in 774 \autoref{eq: zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative.957 \autoref{eq:ZDF_tke_e} or \autoref{eq:ZDF_gls_e} becomes a source term, since $N^2$ is negative. 775 958 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also of the four neighboring values at 776 959 velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$). … … 781 964 because the mixing length scale is bounded by the distance to the sea surface. 782 965 It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme, 783 \ie\ setting the \np{ln \_zdfnpc} namelist parameter to true and784 defining the turbulent closure (\np{ln \_zdftke} or \np{ln\_zdfgls} = \forcode{.true.}) all together.966 \ie\ setting the \np{ln_zdfnpc}{ln\_zdfnpc} namelist parameter to true and 967 defining the turbulent closure (\np{ln_zdftke}{ln\_zdftke} or \np{ln_zdfgls}{ln\_zdfgls} = \forcode{.true.}) all together. 785 968 786 969 The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection, 787 970 %as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp}, 788 therefore \np {ln\_zdfevd}\forcode{ = .false.} should be used with the OSMOSIS scheme.971 therefore \np[=.false.]{ln_zdfevd}{ln\_zdfevd} should be used with the OSMOSIS scheme. 789 972 % gm% + one word on non local flux with KPP scheme trakpp.F90 module... 790 973 791 % ================================================================ 792 % Double Diffusion Mixing 793 % ================================================================ 794 \section[Double diffusion mixing (\forcode{ln_zdfddm = .true.})] 795 {Double diffusion mixing (\protect\np{ln\_zdfddm}\forcode{ = .true.})} 974 %% ================================================================================================= 975 \section[Double diffusion mixing (\forcode{ln_zdfddm})]{Double diffusion mixing (\protect\np{ln_zdfddm}{ln\_zdfddm})} 796 976 \label{subsec:ZDF_ddm} 797 977 798 799 %-------------------------------------------namzdf_ddm-------------------------------------------------800 %801 978 %\nlst{namzdf_ddm} 802 %--------------------------------------------------------------------------------------------------------------803 979 804 980 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter 805 \np{ln \_zdfddm} in \nam{zdf}.981 \np{ln_zdfddm}{ln\_zdfddm} in \nam{zdf}{zdf}. 806 982 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. 807 983 The former condition leads to salt fingering and the latter to diffusive convection. … … 811 987 temperature and salinity. 812 988 813 814 989 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 815 990 \begin{align*} 816 % \label{eq: zdfddm_Kz}991 % \label{eq:ZDF_ddm_Kz} 817 992 &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\ 818 993 &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS} … … 826 1001 (1981): 827 1002 \begin{align} 828 \label{eq: zdfddm_f}1003 \label{eq:ZDF_ddm_f} 829 1004 A_f^{vS} &= 830 1005 \begin{cases} … … 832 1007 0 &\text{otherwise} 833 1008 \end{cases} 834 \\ \label{eq: zdfddm_f_T}1009 \\ \label{eq:ZDF_ddm_f_T} 835 1010 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 836 1011 \end{align} 837 1012 838 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>839 1013 \begin{figure}[!t] 840 \begin{center} 841 \includegraphics[width=\textwidth]{Fig_zdfddm} 842 \caption{ 843 \protect\label{fig:zdfddm} 844 From \citet{merryfield.holloway.ea_JPO99} : 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} 1014 \centering 1015 \includegraphics[width=0.66\textwidth]{Fig_zdfddm} 1016 \caption[Diapycnal diffusivities for temperature and salt in regions of salt fingering and 1017 diffusive convection]{ 1018 From \citet{merryfield.holloway.ea_JPO99}: 1019 (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in 1020 regions of salt fingering. 1021 Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and 1022 thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$; 1023 (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in 1024 regions of diffusive convection. 1025 Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation. 1026 The latter is not implemented in \NEMO.} 1027 \label{fig:ZDF_ddm} 853 1028 \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 1029 1030 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 857 1031 buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{mcdougall.taylor_JMR84}). 858 1032 Following \citet{merryfield.holloway.ea_JPO99}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. … … 861 1035 Federov (1988) is used: 862 1036 \begin{align} 863 % \label{eq: zdfddm_d}1037 % \label{eq:ZDF_ddm_d} 864 1038 A_d^{vT} &= 865 1039 \begin{cases} … … 869 1043 \end{cases} 870 1044 \nonumber \\ 871 \label{eq: zdfddm_d_S}1045 \label{eq:ZDF_ddm_d_S} 872 1046 A_d^{vS} &= 873 1047 \begin{cases} … … 878 1052 \end{align} 879 1053 880 The dependencies of \autoref{eq: zdfddm_f} to \autoref{eq:zdfddm_d_S} on $R_\rho$ are illustrated in881 \autoref{fig: zdfddm}.1054 The dependencies of \autoref{eq:ZDF_ddm_f} to \autoref{eq:ZDF_ddm_d_S} on $R_\rho$ are illustrated in 1055 \autoref{fig:ZDF_ddm}. 882 1056 Implementing this requires computing $R_\rho$ at each grid point on every time step. 883 1057 This is done in \mdl{eosbn2} at the same time as $N^2$ is computed. 884 1058 This avoids duplication in the computation of $\alpha$ and $\beta$ (which is usually quite expensive). 885 1059 886 % ================================================================ 887 % Bottom Friction 888 % ================================================================ 889 \section[Bottom and top friction (\textit{zdfdrg.F90})] 890 {Bottom and top friction (\protect\mdl{zdfdrg})} 891 \label{sec:ZDF_drg} 892 893 %--------------------------------------------nambfr-------------------------------------------------------- 894 % 895 \nlst{namdrg} 896 \nlst{namdrg_top} 897 \nlst{namdrg_bot} 898 899 %-------------------------------------------------------------------------------------------------------------- 900 901 Options to define the top and bottom friction are defined through the \nam{drg} namelist variables. 1060 %% ================================================================================================= 1061 \section[Bottom and top friction (\textit{zdfdrg.F90})]{Bottom and top friction (\protect\mdl{zdfdrg})} 1062 \label{sec:ZDF_drg} 1063 1064 \begin{listing} 1065 \nlst{namdrg} 1066 \caption{\forcode{&namdrg}} 1067 \label{lst:namdrg} 1068 \end{listing} 1069 \begin{listing} 1070 \nlst{namdrg_top} 1071 \caption{\forcode{&namdrg_top}} 1072 \label{lst:namdrg_top} 1073 \end{listing} 1074 \begin{listing} 1075 \nlst{namdrg_bot} 1076 \caption{\forcode{&namdrg_bot}} 1077 \label{lst:namdrg_bot} 1078 \end{listing} 1079 1080 Options to define the top and bottom friction are defined through the \nam{drg}{drg} namelist variables. 902 1081 The bottom friction represents the friction generated by the bathymetry. 903 1082 The top friction represents the friction generated by the ice shelf/ocean interface. … … 905 1084 the description below considers mostly the bottom friction case, if not stated otherwise. 906 1085 907 908 1086 Both the surface momentum flux (wind stress) and the bottom momentum flux (bottom friction) enter the equations as 909 1087 a condition on the vertical diffusive flux. 910 1088 For the bottom boundary layer, one has: 911 1089 \[ 912 % \label{eq: zdfbfr_flux}1090 % \label{eq:ZDF_bfr_flux} 913 1091 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 914 1092 \] … … 926 1104 To illustrate this, consider the equation for $u$ at $k$, the last ocean level: 927 1105 \begin{equation} 928 \label{eq: zdfdrg_flux2}1106 \label{eq:ZDF_drg_flux2} 929 1107 \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}} 930 1108 \end{equation} … … 946 1124 These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where: 947 1125 \begin{equation} 948 \label{eq: zdfbfr_bdef}1126 \label{eq:ZDF_bfr_bdef} 949 1127 \frac{\partial {\textbf U_h}}{\partial t} = 950 1128 - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b … … 953 1131 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. 954 1132 955 % ------------------------------------------------------------------------------------------------------------- 956 % Linear Bottom Friction 957 % ------------------------------------------------------------------------------------------------------------- 958 \subsection[Linear top/bottom friction (\forcode{ln_lin = .true.})] 959 {Linear top/bottom friction (\protect\np{ln\_lin}\forcode{ = .true.)}} 960 \label{subsec:ZDF_drg_linear} 1133 %% ================================================================================================= 1134 \subsection[Linear top/bottom friction (\forcode{ln_lin})]{Linear top/bottom friction (\protect\np{ln_lin}{ln\_lin})} 1135 \label{subsec:ZDF_drg_linear} 961 1136 962 1137 The linear friction parameterisation (including the special case of a free-slip condition) assumes that 963 1138 the friction is proportional to the interior velocity (\ie\ the velocity of the first/last model level): 964 1139 \[ 965 % \label{eq: zdfbfr_linear}1140 % \label{eq:ZDF_bfr_linear} 966 1141 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 967 1142 \] … … 976 1151 and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 977 1152 This is the default value used in \NEMO. It corresponds to a decay time scale of 115~days. 978 It can be changed by specifying \np{rn \_Uc0} (namelist parameter).979 980 For the linear friction case the drag coefficient used in the general expression \autoref{eq: zdfbfr_bdef} is:981 \[ 982 % \label{eq: zdfbfr_linbfr_b}1153 It can be changed by specifying \np{rn_Uc0}{rn\_Uc0} (namelist parameter). 1154 1155 For the linear friction case the drag coefficient used in the general expression \autoref{eq:ZDF_bfr_bdef} is: 1156 \[ 1157 % \label{eq:ZDF_bfr_linbfr_b} 983 1158 c_b^T = - r 984 1159 \] 985 When \np {ln\_lin} \forcode{= .true.}, the value of $r$ used is \np{rn\_Uc0}*\np{rn\_Cd0}.986 Setting \np {ln\_OFF} \forcode{= .true.} (and \forcode{ln_lin =.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.1160 When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}. 1161 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. 987 1162 988 1163 These values are assigned in \mdl{zdfdrg}. 989 1164 Note that there is support for local enhancement of these values via an externally defined 2D mask array 990 (\np {ln\_boost}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file.1165 (\np[=.true.]{ln_boost}{ln\_boost}) given in the \ifile{bfr\_coef} input NetCDF file. 991 1166 The mask values should vary from 0 to 1. 992 1167 Locations with a non-zero mask value will have the friction coefficient increased by 993 $mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}. 994 995 % ------------------------------------------------------------------------------------------------------------- 996 % Non-Linear Bottom Friction 997 % ------------------------------------------------------------------------------------------------------------- 998 \subsection[Non-linear top/bottom friction (\forcode{ln_non_lin = .true.})] 999 {Non-linear top/bottom friction (\protect\np{ln\_non\_lin}\forcode{ = .true.})} 1000 \label{subsec:ZDF_drg_nonlinear} 1168 $mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}. 1169 1170 %% ================================================================================================= 1171 \subsection[Non-linear top/bottom friction (\forcode{ln_non_lin})]{Non-linear top/bottom friction (\protect\np{ln_non_lin}{ln\_non\_lin})} 1172 \label{subsec:ZDF_drg_nonlinear} 1001 1173 1002 1174 The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic: 1003 1175 \[ 1004 % \label{eq: zdfdrg_nonlinear}1176 % \label{eq:ZDF_drg_nonlinear} 1005 1177 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 1006 1178 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b … … 1012 1184 $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 1013 1185 $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. 1014 The CME choices have been set as default values (\np{rn \_Cd0} and \np{rn\_ke0} namelist parameters).1186 The CME choices have been set as default values (\np{rn_Cd0}{rn\_Cd0} and \np{rn_ke0}{rn\_ke0} namelist parameters). 1015 1187 1016 1188 As for the linear case, the friction is imposed in the code by adding the trend due to … … 1018 1190 For the non-linear friction case the term computed in \mdl{zdfdrg} is: 1019 1191 \[ 1020 % \label{eq: zdfdrg_nonlinbfr}1192 % \label{eq:ZDF_drg_nonlinbfr} 1021 1193 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} 1022 1194 \] 1023 1195 1024 1196 The coefficients that control the strength of the non-linear friction are initialised as namelist parameters: 1025 $C_D$= \np{rn \_Cd0}, and $e_b$ =\np{rn\_bfeb2}.1026 Note that for applications which consider tides explicitly, a low or even zero value of \np{rn \_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array1027 (\np {ln\_boost}\forcode{ = .true.}).1197 $C_D$= \np{rn_Cd0}{rn\_Cd0}, and $e_b$ =\np{rn_bfeb2}{rn\_bfeb2}. 1198 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 1199 (\np[=.true.]{ln_boost}{ln\_boost}). 1028 1200 This works in the same way as for the linear friction case with non-zero masked locations increased by 1029 $mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}. 1030 1031 % ------------------------------------------------------------------------------------------------------------- 1032 % Bottom Friction Log-layer 1033 % ------------------------------------------------------------------------------------------------------------- 1034 \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer = .true.})] 1035 {Log-layer top/bottom friction (\protect\np{ln\_loglayer}\forcode{ = .true.})} 1036 \label{subsec:ZDF_drg_loglayer} 1201 $mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}. 1202 1203 %% ================================================================================================= 1204 \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer})]{Log-layer top/bottom friction (\protect\np{ln_loglayer}{ln\_loglayer})} 1205 \label{subsec:ZDF_drg_loglayer} 1037 1206 1038 1207 In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using 1039 1208 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. 1040 If \np {ln\_loglayer} \forcode{= .true.}, $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness):1209 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): 1041 1210 \[ 1042 1211 C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2 1043 1212 \] 1044 1213 1045 \noindent where $\kappa$ is the von-Karman constant and \np{rn \_z0} is a roughness length provided via the namelist.1214 \noindent where $\kappa$ is the von-Karman constant and \np{rn_z0}{rn\_z0} is a roughness length provided via the namelist. 1046 1215 1047 1216 The drag coefficient is bounded such that it is kept greater or equal to 1048 the base \np{rn \_Cd0} value which occurs where layer thicknesses become large and presumably logarithmic layers are not resolved at all. For stability reason, it is also not allowed to exceed the value of an additional namelist parameter:1049 \np{rn \_Cdmax}, \ie1217 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: 1218 \np{rn_Cdmax}{rn\_Cdmax}, \ie 1050 1219 \[ 1051 1220 rn\_Cd0 \leq C_D \leq rn\_Cdmax … … 1053 1222 1054 1223 \noindent The log-layer enhancement can also be applied to the top boundary friction if 1055 under ice-shelf cavities are activated (\np{ln\_isfcav}\forcode{ = .true.}). 1056 %In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}. 1057 1058 % ------------------------------------------------------------------------------------------------------------- 1059 % Explicit bottom Friction 1060 % ------------------------------------------------------------------------------------------------------------- 1061 \subsection{Explicit top/bottom friction (\forcode{ln_drgimp = .false.})} 1062 \label{subsec:ZDF_drg_stability} 1063 1064 Setting \np{ln\_drgimp} \forcode{= .false.} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads: 1224 under ice-shelf cavities are activated (\np[=.true.]{ln_isfcav}{ln\_isfcav}). 1225 %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}. 1226 1227 %% ================================================================================================= 1228 \subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np[=.false.]{ln_drgimp}{ln\_drgimp})} 1229 \label{subsec:ZDF_drg_stability} 1230 1231 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: 1065 1232 1066 1233 At the top (below an ice shelf cavity): … … 1077 1244 1078 1245 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. 1079 For the purposes of stability analysis, an approximation to \autoref{eq: zdfdrg_flux2} is:1246 For the purposes of stability analysis, an approximation to \autoref{eq:ZDF_drg_flux2} is: 1080 1247 \begin{equation} 1081 \label{eq: Eqn_drgstab}1248 \label{eq:ZDF_Eqn_drgstab} 1082 1249 \begin{split} 1083 1250 \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt \\ … … 1090 1257 |\Delta u| < \;|u| 1091 1258 \] 1092 \noindent which, using \autoref{eq: Eqn_drgstab}, gives:1259 \noindent which, using \autoref{eq:ZDF_Eqn_drgstab}, gives: 1093 1260 \[ 1094 1261 r\frac{2\rdt}{e_{3u}} < 1 \qquad \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ … … 1117 1284 The number of potential breaches of the explicit stability criterion are still reported for information purposes. 1118 1285 1119 % ------------------------------------------------------------------------------------------------------------- 1120 % Implicit Bottom Friction 1121 % ------------------------------------------------------------------------------------------------------------- 1122 \subsection[Implicit top/bottom friction (\forcode{ln_drgimp = .true.})] 1123 {Implicit top/bottom friction (\protect\np{ln\_drgimp}\forcode{ = .true.})} 1124 \label{subsec:ZDF_drg_imp} 1286 %% ================================================================================================= 1287 \subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np[=.true.]{ln_drgimp}{ln\_drgimp})} 1288 \label{subsec:ZDF_drg_imp} 1125 1289 1126 1290 An optional implicit form of bottom friction has been implemented to improve model stability. 1127 1291 We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting. 1128 This option can be invoked by setting \np{ln \_drgimp} to \forcode{.true.} in the \nam{drg} namelist.1129 %This option requires \np{ln \_zdfexp} to be \forcode{.false.} in the \nam{zdf} namelist.1292 This option can be invoked by setting \np{ln_drgimp}{ln\_drgimp} to \forcode{.true.} in the \nam{drg}{drg} namelist. 1293 %This option requires \np{ln_zdfexp}{ln\_zdfexp} to be \forcode{.false.} in the \nam{zdf}{zdf} namelist. 1130 1294 1131 1295 This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step: … … 1133 1297 At the top (below an ice shelf cavity): 1134 1298 \[ 1135 % \label{eq: dynzdf_drg_top}1299 % \label{eq:ZDF_dynZDF__drg_top} 1136 1300 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1137 1301 = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t} … … 1140 1304 At the bottom (above the sea floor): 1141 1305 \[ 1142 % \label{eq: dynzdf_drg_bot}1306 % \label{eq:ZDF_dynZDF__drg_bot} 1143 1307 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1144 1308 = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b} … … 1148 1312 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit. 1149 1313 1150 % ------------------------------------------------------------------------------------------------------------- 1151 % Bottom Friction with split-explicit free surface 1152 % ------------------------------------------------------------------------------------------------------------- 1153 \subsection[Bottom friction with split-explicit free surface] 1154 {Bottom friction with split-explicit free surface} 1155 \label{subsec:ZDF_drg_ts} 1156 1157 With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \np{ln\_drgimp}\forcode{= .false.} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie\ \forcode{ln_bt_fw = .false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \np{ln\_drgimp}\forcode{= .true.}, stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions. 1314 %% ================================================================================================= 1315 \subsection[Bottom friction with split-explicit free surface]{Bottom friction with split-explicit free surface} 1316 \label{subsec:ZDF_drg_ts} 1317 1318 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. 1158 1319 1159 1320 The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO\ is as follows: … … 1165 1326 Note that other strategies are possible, like considering vertical diffusion step in advance, \ie\ prior barotropic integration. 1166 1327 1167 1168 % ================================================================ 1169 % Internal wave-driven mixing 1170 % ================================================================ 1171 \section[Internal wave-driven mixing (\forcode{ln_zdfiwm = .true.})] 1172 {Internal wave-driven mixing (\protect\np{ln\_zdfiwm}\forcode{ = .true.})} 1328 %% ================================================================================================= 1329 \section[Internal wave-driven mixing (\forcode{ln_zdfiwm})]{Internal wave-driven mixing (\protect\np{ln_zdfiwm}{ln\_zdfiwm})} 1173 1330 \label{subsec:ZDF_tmx_new} 1174 1331 1175 %--------------------------------------------namzdf_iwm------------------------------------------ 1176 % 1177 \nlst{namzdf_iwm} 1178 %-------------------------------------------------------------------------------------------------------------- 1332 \begin{listing} 1333 \nlst{namzdf_iwm} 1334 \caption{\forcode{&namzdf_iwm}} 1335 \label{lst:namzdf_iwm} 1336 \end{listing} 1179 1337 1180 1338 The parameterization of mixing induced by breaking internal waves is a generalization of … … 1183 1341 and the resulting diffusivity is obtained as 1184 1342 \[ 1185 % \label{eq: Kwave}1343 % \label{eq:ZDF_Kwave} 1186 1344 A^{vT}_{wave} = R_f \,\frac{ \epsilon }{ \rho \, N^2 } 1187 1345 \] 1188 1346 where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of 1189 1347 the energy available for mixing. 1190 If the \np{ln \_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and1348 If the \np{ln_mevar}{ln\_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and 1191 1349 equal to 1/6 \citep{osborn_JPO80}. 1192 1350 In the opposite (recommended) case, $R_f$ is instead a function of … … 1198 1356 1199 1357 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 1200 as a function of $Re_b$ by setting the \np{ln \_tsdiff} parameter to \forcode{.true.}, a recommended choice.1358 as a function of $Re_b$ by setting the \np{ln_tsdiff}{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice. 1201 1359 This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14}, 1202 1360 is implemented as in \cite{de-lavergne.madec.ea_JPO16}. … … 1216 1374 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz' } \; , 1217 1375 \] 1218 The $n_p$ parameter (given by \np{nn \_zpyc} in \nam{zdf\_iwm} namelist)1376 The $n_p$ parameter (given by \np{nn_zpyc}{nn\_zpyc} in \nam{zdf_iwm}{zdf\_iwm} namelist) 1219 1377 controls the stratification-dependence of the pycnocline-intensified dissipation. 1220 1378 It can take values of $1$ (recommended) or $2$. … … 1224 1382 $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of 1225 1383 the abyssal hill topography \citep{goff_JGR10} and the latitude. 1226 %1227 1384 % Jc: input files names ? 1228 1385 1229 % ================================================================ 1230 % surface wave-induced mixing 1231 % ================================================================ 1232 \section[Surface wave-induced mixing (\forcode{ln_zdfswm = .true.})] 1233 {Surface wave-induced mixing (\protect\np{ln\_zdfswm}\forcode{ = .true.})} 1386 %% ================================================================================================= 1387 \section[Surface wave-induced mixing (\forcode{ln_zdfswm})]{Surface wave-induced mixing (\protect\np{ln_zdfswm}{ln\_zdfswm})} 1234 1388 \label{subsec:ZDF_swm} 1235 1389 … … 1242 1396 1243 1397 \begin{equation} 1244 \label{eq: Bv}1398 \label{eq:ZDF_Bv} 1245 1399 B_{v} = \alpha {A} {U}_{st} {exp(3kz)} 1246 1400 \end{equation} … … 1254 1408 and diffusivity coefficients. 1255 1409 1256 In order to account for this contribution set: \forcode{ln_zdfswm =.true.},1257 then wave interaction has to be activated through \forcode{ln_wave =.true.},1258 the Stokes Drift can be evaluated by setting \forcode{ln_sdw =.true.}1410 In order to account for this contribution set: \forcode{ln_zdfswm=.true.}, 1411 then wave interaction has to be activated through \forcode{ln_wave=.true.}, 1412 the Stokes Drift can be evaluated by setting \forcode{ln_sdw=.true.} 1259 1413 (see \autoref{subsec:SBC_wave_sdw}) 1260 1414 and the needed wave fields can be provided either in forcing or coupled mode 1261 1415 (for more information on wave parameters and settings see \autoref{sec:SBC_wave}) 1262 1416 1263 % ================================================================ 1264 % Adaptive-implicit vertical advection 1265 % ================================================================ 1266 \section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp = .true.})] 1267 {Adaptive-implicit vertical advection(\protect\np{ln\_zad\_Aimp}\forcode{ = .true.})} 1417 %% ================================================================================================= 1418 \section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp})]{Adaptive-implicit vertical advection(\protect\np{ln_zad_Aimp}{ln\_zad\_Aimp})} 1268 1419 \label{subsec:ZDF_aimp} 1269 1420 … … 1276 1427 criteria for a range of advection schemes. The values for the Leap-Frog with Robert 1277 1428 asselin filter time-stepping (as used in NEMO) are reproduced in 1278 \autoref{tab: zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these1429 \autoref{tab:ZDF_zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these 1279 1430 restrictions but at the cost of large dispersive errors and, possibly, large numerical 1280 1431 viscosity. The adaptive-implicit vertical advection option provides a targetted use of the … … 1283 1434 interest or due to short-lived conditions such that the extra numerical diffusion or 1284 1435 viscosity does not greatly affect the overall solution. With such applications, setting: 1285 \forcode{ln_zad_Aimp =.true.} should allow much longer model timesteps to be used whilst1436 \forcode{ln_zad_Aimp=.true.} should allow much longer model timesteps to be used whilst 1286 1437 retaining the accuracy of the high order explicit schemes over most of the domain. 1287 1438 1288 1439 \begin{table}[htbp] 1289 \begin{center} 1290 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}} 1291 \begin{tabular}{r|ccc} 1292 \hline 1293 spatial discretization & 2nd order centered & 3rd order upwind & 4th order compact \\ 1294 advective CFL criterion & 0.904 & 0.472 & 0.522 \\ 1295 \hline 1296 \end{tabular} 1297 \caption{ 1298 \protect\label{tab:zad_Aimp_CFLcrit} 1299 The advective CFL criteria for a range of spatial discretizations for the Leap-Frog with Robert Asselin filter time-stepping 1300 ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}. 1301 } 1302 \end{center} 1440 \centering 1441 % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}} 1442 \begin{tabular}{r|ccc} 1443 \hline 1444 spatial discretization & 2$^nd$ order centered & 3$^rd$ order upwind & 4$^th$ order compact \\ 1445 advective CFL criterion & 0.904 & 0.472 & 0.522 \\ 1446 \hline 1447 \end{tabular} 1448 \caption[Advective CFL criteria for the leapfrog with Robert Asselin filter time-stepping]{ 1449 The advective CFL criteria for a range of spatial discretizations for 1450 the leapfrog with Robert Asselin filter time-stepping 1451 ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}.} 1452 \label{tab:ZDF_zad_Aimp_CFLcrit} 1303 1453 \end{table} 1304 1454 … … 1313 1463 1314 1464 \begin{equation} 1315 \label{eq: Eqn_zad_Aimp_Courant}1465 \label{eq:ZDF_Eqn_zad_Aimp_Courant} 1316 1466 \begin{split} 1317 1467 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 ] \\ … … 1326 1476 1327 1477 \begin{align} 1328 \label{eq: Eqn_zad_Aimp_partition}1478 \label{eq:ZDF_Eqn_zad_Aimp_partition} 1329 1479 Cu_{min} &= 0.15 \nonumber \\ 1330 1480 Cu_{max} &= 0.3 \nonumber \\ 1331 1481 Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ 1332 1482 Fcu &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ 1333 C\kern-0.14emf &=1483 \cf &= 1334 1484 \begin{cases} 1335 1485 0.0 &\text{if $Cu \leq Cu_{min}$} \\ … … 1340 1490 1341 1491 \begin{figure}[!t] 1342 \begin{center} 1343 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_coeff} 1344 \caption{ 1345 \protect\label{fig:zad_Aimp_coeff} 1346 The value of the partitioning coefficient ($C\kern-0.14em f$) used to partition vertical velocities into parts to 1347 be treated implicitly and explicitly for a range of typical Courant numbers (\forcode{ln_zad_Aimp=.true.}) 1348 } 1349 \end{center} 1492 \centering 1493 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_coeff} 1494 \caption[Partitioning coefficient used to partition vertical velocities into parts]{ 1495 The value of the partitioning coefficient (\cf) used to partition vertical velocities into 1496 parts to be treated implicitly and explicitly for a range of typical Courant numbers 1497 (\forcode{ln_zad_Aimp=.true.}).} 1498 \label{fig:ZDF_zad_Aimp_coeff} 1350 1499 \end{figure} 1351 1500 … … 1355 1504 1356 1505 \begin{align} 1357 \label{eq: Eqn_zad_Aimp_partition2}1358 w_{i_{ijk}} &= C\kern-0.14emf_{ijk} w_{n_{ijk}} \nonumber \\1359 w_{n_{ijk}} &= (1- C\kern-0.14em f_{ijk}) w_{n_{ijk}}1506 \label{eq:ZDF_Eqn_zad_Aimp_partition2} 1507 w_{i_{ijk}} &= \cf_{ijk} w_{n_{ijk}} \nonumber \\ 1508 w_{n_{ijk}} &= (1-\cf_{ijk}) w_{n_{ijk}} 1360 1509 \end{align} 1361 1510 1362 1511 \noindent Note that the coefficient is such that the treatment is never fully implicit; 1363 the three cases from \autoref{eq: Eqn_zad_Aimp_partition} can be considered as:1512 the three cases from \autoref{eq:ZDF_Eqn_zad_Aimp_partition} can be considered as: 1364 1513 fully-explicit; mixed explicit/implicit and mostly-implicit. With the settings shown the 1365 coefficient ( $C\kern-0.14em f$) varies as shown in \autoref{fig:zad_Aimp_coeff}. Note with these values1514 coefficient (\cf) varies as shown in \autoref{fig:ZDF_zad_Aimp_coeff}. Note with these values 1366 1515 the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly 1367 1516 implicit' is 0.45 which is just below the stability limited given in 1368 \autoref{tab: zad_Aimp_CFLcrit} for a 3rd order scheme.1517 \autoref{tab:ZDF_zad_Aimp_CFLcrit} for a 3rd order scheme. 1369 1518 1370 1519 The $w_i$ component is added to the implicit solvers for the vertical mixing in … … 1376 1525 vertical fluxes are then removed since they are added by the implicit solver later on. 1377 1526 1378 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 1527 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 1379 1528 used in a wide range of simulations. The following test simulation, however, does illustrate 1380 1529 the potential benefits and will hopefully encourage further testing and feedback from users: 1381 1530 1382 1531 \begin{figure}[!t] 1383 \begin{center} 1384 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} 1385 \caption{ 1386 \protect\label{fig:zad_Aimp_overflow_frames} 1387 A time-series of temperature vertical cross-sections for the OVERFLOW test case. These results are for the default 1388 settings with \forcode{nn_rdt=10.0} and without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}). 1389 } 1390 \end{center} 1532 \centering 1533 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} 1534 \caption[OVERFLOW: time-series of temperature vertical cross-sections]{ 1535 A time-series of temperature vertical cross-sections for the OVERFLOW test case. 1536 These results are for the default settings with \forcode{nn_rdt=10.0} and 1537 without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}).} 1538 \label{fig:ZDF_zad_Aimp_overflow_frames} 1391 1539 \end{figure} 1392 1540 1541 %% ================================================================================================= 1393 1542 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} 1543 1394 1544 The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case} 1395 1545 provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case … … 1407 1557 1408 1558 \noindent which were chosen to provide a slightly more stable and less noisy solution. The 1409 result when using the default value of \forcode{nn_rdt =10.} without adaptive-implicit1410 vertical velocity is illustrated in \autoref{fig: zad_Aimp_overflow_frames}. The mass of1559 result when using the default value of \forcode{nn_rdt=10.} without adaptive-implicit 1560 vertical velocity is illustrated in \autoref{fig:ZDF_zad_Aimp_overflow_frames}. The mass of 1411 1561 cold water, initially sitting on the shelf, moves down the slope and forms a 1412 1562 bottom-trapped, dense plume. Even with these extra physics choices the model is close to 1413 stability limits and attempts with \forcode{nn_rdt =30.} will fail after about 5.5 hours1563 stability limits and attempts with \forcode{nn_rdt=30.} will fail after about 5.5 hours 1414 1564 with excessively high horizontal velocities. This time-scale corresponds with the time the 1415 1565 plume reaches the steepest part of the topography and, although detected as a horizontal … … 1418 1568 1419 1569 The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps 1420 are shown in \autoref{fig: zad_Aimp_overflow_all_rdt} (together with the equivalent1570 are shown in \autoref{fig:ZDF_zad_Aimp_overflow_all_rdt} (together with the equivalent 1421 1571 frames from the base run). In this simple example the use of the adaptive-implicit 1422 1572 vertcal advection scheme has enabled a 12x increase in the model timestep without 1423 1573 significantly altering the solution (although at this extreme the plume is more diffuse 1424 1574 and has not travelled so far). Notably, the solution with and without the scheme is 1425 slightly different even with \forcode{nn_rdt =10.}; suggesting that the base run was1575 slightly different even with \forcode{nn_rdt=10.}; suggesting that the base run was 1426 1576 close enough to instability to trigger the scheme despite completing successfully. 1427 1577 To assist in diagnosing how active the scheme is, in both location and time, the 3D 1428 1578 implicit and explicit components of the vertical velocity are available via XIOS as 1429 1579 \texttt{wimp} and \texttt{wexp} respectively. Likewise, the partitioning coefficient 1430 ( $C\kern-0.14em f$) is also available as \texttt{wi\_cff}. For a quick oversight of1580 (\cf) is also available as \texttt{wi\_cff}. For a quick oversight of 1431 1581 the schemes activity the global maximum values of the absolute implicit component 1432 1582 of the vertical velocity and the partitioning coefficient are written to the netCDF … … 1434 1584 \autoref{sec:MISC_opt} for activation details). 1435 1585 1436 \autoref{fig: zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for1586 \autoref{fig:ZDF_zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for 1437 1587 the various overflow tests. Note that the adaptive-implicit vertical advection scheme is 1438 1588 active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the … … 1441 1591 oscillatory nature of this measure appears to be linked to the progress of the plume front 1442 1592 as each cusp is associated with the location of the maximum shifting to the adjacent cell. 1443 This is illustrated in \autoref{fig: zad_Aimp_maxCf_loc} where the i- and k- locations of the1593 This is illustrated in \autoref{fig:ZDF_zad_Aimp_maxCf_loc} where the i- and k- locations of the 1444 1594 maximum have been overlaid for the base run case. 1445 1595 … … 1460 1610 1461 1611 \begin{figure}[!t] 1462 \ begin{center}1463 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_all_rdt}1464 \caption{1465 \protect\label{fig:zad_Aimp_overflow_all_rdt}1466 Sample temperature vertical cross-sections from mid- and end-run using different values for \forcode{nn_rdt}1467 and with or without adaptive implicit vertical advection. Without the adaptive implicit vertical advection only1468 the run with the shortest timestep is able to run to completion. Note also that the colour-scale has been1469 chosen to confirm that temperatures remain within the original range of 10$^o$ to 20$^o$.1470 }1471 \ end{center}1612 \centering 1613 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_overflow_all_rdt} 1614 \caption[OVERFLOW: sample temperature vertical cross-sections from mid- and end-run]{ 1615 Sample temperature vertical cross-sections from mid- and end-run using 1616 different values for \forcode{nn_rdt} and with or without adaptive implicit vertical advection. 1617 Without the adaptive implicit vertical advection 1618 only the run with the shortest timestep is able to run to completion. 1619 Note also that the colour-scale has been chosen to confirm that 1620 temperatures remain within the original range of 10$^o$ to 20$^o$.} 1621 \label{fig:ZDF_zad_Aimp_overflow_all_rdt} 1472 1622 \end{figure} 1473 1623 1474 1624 \begin{figure}[!t] 1475 \begin{center} 1476 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf} 1477 \caption{ 1478 \protect\label{fig:zad_Aimp_maxCf} 1479 The maximum partitioning coefficient during a series of test runs with increasing model timestep length. 1480 At the larger timesteps, the vertical velocity is treated mostly implicitly at some location throughout 1481 the run. 1482 } 1483 \end{center} 1625 \centering 1626 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_maxCf} 1627 \caption[OVERFLOW: maximum partitioning coefficient during a series of test runs]{ 1628 The maximum partitioning coefficient during a series of test runs with 1629 increasing model timestep length. 1630 At the larger timesteps, 1631 the vertical velocity is treated mostly implicitly at some location throughout the run.} 1632 \label{fig:ZDF_zad_Aimp_maxCf} 1484 1633 \end{figure} 1485 1634 1486 1635 \begin{figure}[!t] 1487 \begin{center} 1488 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf_loc} 1489 \caption{ 1490 \protect\label{fig:zad_Aimp_maxCf_loc} 1491 The maximum partitioning coefficient for the \forcode{nn_rdt=10.0s} case overlaid with information on the gridcell i- and k- 1492 locations of the maximum value. 1493 } 1494 \end{center} 1636 \centering 1637 \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_maxCf_loc} 1638 \caption[OVERFLOW: maximum partitioning coefficient for the case overlaid]{ 1639 The maximum partitioning coefficient for the \forcode{nn_rdt=10.0} case overlaid with 1640 information on the gridcell i- and k-locations of the maximum value.} 1641 \label{fig:ZDF_zad_Aimp_maxCf_loc} 1495 1642 \end{figure} 1496 1643 1497 % ================================================================ 1498 1499 \biblio 1500 1501 \pindex 1644 \onlyinsubfile{\input{../../global/epilogue}} 1502 1645 1503 1646 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.