Changeset 11512 for NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/doc/latex/NEMO/subfiles/chap_ZDF.tex
- Timestamp:
- 2019-09-09T12:05:20+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/doc/latex/NEMO/subfiles/chap_ZDF.tex
r11353 r11512 8 8 \label{chap:ZDF} 9 9 10 \ minitoc10 \chaptertoc 11 11 12 12 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. … … 25 25 At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}), 26 26 while at the bottom they are set to zero for heat and salt, 27 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie \np{ln\_trabbc} defined,27 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie\ \np{ln\_trabbc} defined, 28 28 see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum 29 (see \autoref{sec:ZDF_drg}). 29 (see \autoref{sec:ZDF_drg}). 30 30 31 31 In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and … … 37 37 the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules. 38 38 The trends due to the vertical momentum and tracer diffusion, including the surface forcing, 39 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 39 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 40 40 %These trends can be computed using either a forward time stepping scheme 41 41 %(namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme … … 49 49 50 50 % ------------------------------------------------------------------------------------------------------------- 51 % Constant 51 % Constant 52 52 % ------------------------------------------------------------------------------------------------------------- 53 53 \subsection[Constant (\forcode{ln_zdfcst = .true.})] … … 55 55 \label{subsec:ZDF_cst} 56 56 57 Options are defined through the \n gn{namzdf} namelist variables.57 Options are defined through the \nam{zdf} namelist variables. 58 58 When \np{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 59 59 constant values over the whole ocean. … … 66 66 \end{align*} 67 67 68 These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. 68 These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters. 69 69 In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity, 70 70 that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and … … 84 84 85 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 \n gn{namzdf\_ric} namelist variables.87 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. 86 tracer eddy coefficients is set through the \nam{zdf\_ric} namelist variables. 87 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. 88 88 \textit{In situ} measurements have been used to link vertical turbulent activity to large scale ocean structures. 89 89 The hypothesis of a mixing mainly maintained by the growth of Kelvin-Helmholtz like instabilities leads to 90 90 a dependency between the vertical eddy coefficients and the local Richardson number 91 (\ie the ratio of stratification to vertical shear).91 (\ie\ the ratio of stratification to vertical shear). 92 92 Following \citet{pacanowski.philander_JPO81}, the following formulation has been implemented: 93 93 \[ … … 101 101 \] 102 102 where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number, 103 $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 103 $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 104 104 $A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the constant case 105 105 (see \autoref{subsec:ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ is the maximum value that … … 130 130 131 131 % ------------------------------------------------------------------------------------------------------------- 132 % TKE Turbulent Closure Scheme 132 % TKE Turbulent Closure Scheme 133 133 % ------------------------------------------------------------------------------------------------------------- 134 134 \subsection[TKE turbulent closure scheme (\forcode{ln_zdftke = .true.})] … … 144 144 and a closure assumption for the turbulent length scales. 145 145 This turbulent closure model has been developed by \citet{bougeault.lacarrere_MWR89} in the atmospheric case, 146 adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of NEMO,146 adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of \NEMO, 147 147 by \citet{blanke.delecluse_JPO93} for equatorial Atlantic simulations. 148 148 Since then, significant modifications have been introduced by \citet{madec.delecluse.ea_NPM98} in both the implementation and … … 167 167 \end{split} 168 168 \] 169 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 170 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 169 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 170 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 171 171 $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients. 172 172 The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with 173 vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}. 173 vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}. 174 174 They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}. 175 175 $P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$: … … 199 199 too weak vertical diffusion. 200 200 They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and 201 \np{rn\_avt0} (\n gn{namzdf} namelist, see \autoref{subsec:ZDF_cst}).201 \np{rn\_avt0} (\nam{zdf} namelist, see \autoref{subsec:ZDF_cst}). 202 202 203 203 \subsubsection{Turbulent length scale} … … 227 227 \autoref{eq:tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than 228 228 the variations of depth. 229 It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less 229 It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less 230 230 time consuming. 231 231 In particular, it allows the length scale to be limited not only by the distance to the surface or … … 256 256 \end{aligned} 257 257 \] 258 where $l^{(k)}$ is computed using \autoref{eq:tke_mxl0_1}, \ie $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.258 where $l^{(k)}$ is computed using \autoref{eq:tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. 259 259 260 260 In the \np{nn\_mxl}\forcode{ = 2} case, the dissipation and mixing length scales take the same value: … … 283 283 This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow. 284 284 The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and 285 air-sea drag coefficient. 286 The latter concerns the bulk formulae and is not discussed here. 285 air-sea drag coefficient. 286 The latter concerns the bulk formulae and is not discussed here. 287 287 288 288 Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is : … … 292 292 \end{equation} 293 293 where $\alpha_{CB}$ is the \citet{craig.banner_JPO94} constant of proportionality which depends on the ''wave age'', 294 ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}. 294 ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}. 295 295 The boundary condition on the turbulent length scale follows the Charnock's relation: 296 296 \begin{equation} … … 303 303 $\alpha_{CB} = 100$ the Craig and Banner's value. 304 304 As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 305 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds 305 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds 306 306 to $\alpha_{CB} = 100$. 307 Further setting \np{ln\_mxl0 =.true.}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,307 Further setting \np{ln\_mxl0}\forcode{ =.true.}, applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale, 308 308 with $\beta$ hard coded to the Stacey's value. 309 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the 309 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the 310 310 surface $\bar{e}$ value. 311 311 … … 319 319 The detailed physics behind LC is described in, for example, \citet{craik.leibovich_JFM76}. 320 320 The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and 321 wind drift currents. 321 wind drift currents. 322 322 323 323 Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by … … 326 326 an extra source term of TKE, $P_{LC}$. 327 327 The presence of $P_{LC}$ in \autoref{eq:zdftke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to 328 \forcode{.true.} in the \n gn{namzdf\_tke} namelist.329 328 \forcode{.true.} in the \nam{zdf\_tke} namelist. 329 330 330 By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}), 331 $P_{LC}$ is assumed to be : 331 $P_{LC}$ is assumed to be : 332 332 \[ 333 333 P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}} 334 334 \] 335 335 where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth. 336 With no information about the wave field, $w_{LC}$ is assumed to be proportional to 337 the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module 336 With no information about the wave field, $w_{LC}$ is assumed to be proportional to 337 the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module 338 338 \footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as 339 339 $u_s = 0.016 \,|U_{10m}|$. … … 343 343 For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at 344 344 a finite depth $H_{LC}$ (which is often close to the mixed layer depth), 345 and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures). 345 and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures). 346 346 The resulting expression for $w_{LC}$ is : 347 347 \[ … … 355 355 The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second. 356 356 The value of $c_{LC}$ is set through the \np{rn\_lc} namelist parameter, 357 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}. 357 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}. 358 358 359 359 The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations: 360 360 $H_{LC}$ is the depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by 361 converting its kinetic energy to potential energy, according to 361 converting its kinetic energy to potential energy, according to 362 362 \[ 363 363 - \int_{-H_{LC}}^0 { N^2\;z \;dz} = \frac{1}{2} u_s^2 … … 370 370 produce mixed-layer depths that are too shallow during summer months and windy conditions. 371 371 This bias is particularly acute over the Southern Ocean. 372 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}. 373 The parameterization is an empirical one, \ie not derived from theoretical considerations,374 but rather is meant to account for observed processes that affect the density structure of 375 the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme 376 (\ie near-inertial oscillations and ocean swells and waves).377 378 When using this parameterization (\ie when \np{nn\_etau}\forcode{ = 1}),372 To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}. 373 The parameterization is an empirical one, \ie\ not derived from theoretical considerations, 374 but rather is meant to account for observed processes that affect the density structure of 375 the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme 376 (\ie\ near-inertial oscillations and ocean swells and waves). 377 378 When using this parameterization (\ie\ when \np{nn\_etau}\forcode{ = 1}), 379 379 the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations, 380 380 swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition, … … 387 387 penetrates in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of 388 388 the penetration, and $f_i$ is the ice concentration 389 (no penetration if $f_i=1$, \ie if the ocean is entirely covered by sea-ice).389 (no penetration if $f_i=1$, \ie\ if the ocean is entirely covered by sea-ice). 390 390 The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. 391 391 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}\forcode{ = 0}) or 392 392 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}). 393 (\np{nn\_etau}\forcode{ = 1}). 394 394 395 395 Note that two other option exist, \np{nn\_etau}\forcode{ = 2, 3}. 396 396 They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer, 397 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean. 397 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean. 398 398 Those two options are obsolescent features introduced for test purposes. 399 They will be removed in the next release. 399 They will be removed in the next release. 400 400 401 401 % This should be explain better below what this rn_eice parameter is meant for: 402 402 In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn\_eice} namelist parameter. 403 403 This parameter varies from \forcode{0} for no effect to \forcode{4} to suppress the TKE input into the ocean when Sea Ice concentration 404 is greater than 25\%. 405 406 % from Burchard et al OM 2008 : 407 % the most critical process not reproduced by statistical turbulence models is the activity of 408 % internal waves and their interaction with turbulence. After the Reynolds decomposition, 409 % internal waves are in principle included in the RANS equations, but later partially 410 % excluded by the hydrostatic assumption and the model resolution. 411 % Thus far, the representation of internal wave mixing in ocean models has been relatively crude 412 % (\eg Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).413 414 % ------------------------------------------------------------------------------------------------------------- 415 % GLS Generic Length Scale Scheme 404 is greater than 25\%. 405 406 % from Burchard et al OM 2008 : 407 % the most critical process not reproduced by statistical turbulence models is the activity of 408 % internal waves and their interaction with turbulence. After the Reynolds decomposition, 409 % internal waves are in principle included in the RANS equations, but later partially 410 % excluded by the hydrostatic assumption and the model resolution. 411 % Thus far, the representation of internal wave mixing in ocean models has been relatively crude 412 % (\eg\ Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002). 413 414 % ------------------------------------------------------------------------------------------------------------- 415 % GLS Generic Length Scale Scheme 416 416 % ------------------------------------------------------------------------------------------------------------- 417 417 \subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls = .true.})] … … 427 427 one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale, 428 428 $\psi$ \citep{umlauf.burchard_JMR03, umlauf.burchard_CSR05}. 429 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 429 This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$, 430 430 where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:GLS} allows to recover a number of 431 431 well-known turbulent closures ($k$-$kl$ \citep{mellor.yamada_RG82}, $k$-$\epsilon$ \citep{rodi_JGR87}, 432 $k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}). 432 $k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}). 433 433 The GLS scheme is given by the following set of equations: 434 434 \begin{equation} … … 467 467 \] 468 468 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and 469 $\epsilon$ the dissipation rate. 469 $\epsilon$ the dissipation rate. 470 470 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of 471 471 the choice of the turbulence model. 472 472 Four different turbulent models are pre-defined (\autoref{tab:GLS}). 473 They are made available through the \np{nn\_clo} namelist parameter. 473 They are made available through the \np{nn\_clo} namelist parameter. 474 474 475 475 %--------------------------------------------------TABLE-------------------------------------------------- … … 497 497 \protect\label{tab:GLS} 498 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\n gn{namzdf\_gls}.499 \protect\np{ln\_zdfgls}\forcode{ = .true.} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\nam{zdf\_gls}. 500 500 } 501 501 \end{center} … … 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{nn\_stab\_func}\forcode{ = 0, 3}, resp.). 511 511 The value of $C_{0\mu}$ depends on the choice of the stability function. 512 512 … … 516 516 \np{rn\_crban}\forcode{ > 0.} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}. 517 517 The \np{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and 518 \np{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}. 518 \np{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 … … 531 531 the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}. 532 532 Evaluation of the 4 GLS turbulent closure schemes can be found in \citet{warner.sherwood.ea_OM05} in ROMS model and 533 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO model.534 535 536 % ------------------------------------------------------------------------------------------------------------- 537 % OSM OSMOSIS BL Scheme 533 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO\ model. 534 535 536 % ------------------------------------------------------------------------------------------------------------- 537 % OSM OSMOSIS BL Scheme 538 538 % ------------------------------------------------------------------------------------------------------------- 539 539 \subsection[OSM: OSMosis boundary layer scheme (\forcode{ln_zdfosm = .true.})] … … 562 562 Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and its links to the momentum and tracer time integration. 563 563 } 564 \end{center} 564 \end{center} 565 565 \end{figure} 566 566 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 577 577 With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to 578 578 the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and 579 summing the result vertically: 579 summing the result vertically: 580 580 \begin{equation} 581 581 \label{eq:energ1} … … 613 613 \end{split} 614 614 \end{equation} 615 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 615 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$. 616 616 The first term of the right hand side of \autoref{eq:energ2} is always zero because 617 617 there is no diffusive flux through the ocean surface and bottom). … … 652 652 %The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as 653 653 %they all appear in the right hand side of \autoref{eq:zdftke_ene}. 654 %For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 654 %For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored. 655 655 656 656 % ================================================================ … … 660 660 \label{sec:ZDF_conv} 661 661 662 Static instabilities (\ie light potential densities under heavy ones) may occur at particular ocean grid points.662 Static instabilities (\ie\ light potential densities under heavy ones) may occur at particular ocean grid points. 663 663 In nature, convective processes quickly re-establish the static stability of the water column. 664 664 These processes have been removed from the model via the hydrostatic assumption so they must be parameterized. … … 668 668 669 669 % ------------------------------------------------------------------------------------------------------------- 670 % Non-Penetrative Convective Adjustment 670 % Non-Penetrative Convective Adjustment 671 671 % ------------------------------------------------------------------------------------------------------------- 672 672 \subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc = .true.})] … … 696 696 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 697 697 698 Options are defined through the \n gn{namzdf} namelist variables.698 Options are defined through the \nam{zdf} namelist variables. 699 699 The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}\forcode{ = .true.}. 700 700 It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of 701 701 the water column, but only until the density structure becomes neutrally stable 702 (\ie until the mixed portion of the water column has \textit{exactly} the density of the water just below)702 (\ie\ until the mixed portion of the water column has \textit{exactly} the density of the water just below) 703 703 \citep{madec.delecluse.ea_JPO91}. 704 704 The associated algorithm is an iterative process used in the following way (\autoref{fig:npc}): … … 717 717 This algorithm is significantly different from mixing statically unstable levels two by two. 718 718 The latter procedure cannot converge with a finite number of iterations for some vertical profiles while 719 the algorithm used in \NEMO converges for any profile in a number of iterations which is less than719 the algorithm used in \NEMO\ converges for any profile in a number of iterations which is less than 720 720 the number of vertical levels. 721 721 This property is of paramount importance as pointed out by \citet{killworth_iprc89}: … … 727 727 (L. Brodeau, personnal communication). 728 728 Two main differences have been introduced compared to the original algorithm: 729 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 730 (not the difference in potential density); 729 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency 730 (not the difference in potential density); 731 731 $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in 732 732 the same way their temperature and salinity has been mixed. … … 735 735 736 736 % ------------------------------------------------------------------------------------------------------------- 737 % Enhanced Vertical Diffusion 737 % Enhanced Vertical Diffusion 738 738 % ------------------------------------------------------------------------------------------------------------- 739 739 \subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd = .true.})] … … 741 741 \label{subsec:ZDF_evd} 742 742 743 Options are defined through the \n gn{namzdf} namelist variables.743 Options are defined through the \nam{zdf} namelist variables. 744 744 The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}\forcode{ = .true.}. 745 In this case, the vertical eddy mixing coefficients are assigned very large values 745 In this case, the vertical eddy mixing coefficients are assigned very large values 746 746 in regions where the stratification is unstable 747 (\ie when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}.747 (\ie\ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}. 748 748 This is done either on tracers only (\np{nn\_evdm}\forcode{ = 0}) or 749 749 on both momentum and tracers (\np{nn\_evdm}\forcode{ = 1}). … … 762 762 763 763 % ------------------------------------------------------------------------------------------------------------- 764 % Turbulent Closure Scheme 765 % ------------------------------------------------------------------------------------------------------------- 766 \subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf/tke/gls/osm = .true.})] 767 {Handling convection with turbulent closure schemes (\protect\np{ln\_zdf/tke/gls/osm}\forcode{ = .true.})} 764 % Turbulent Closure Scheme 765 % ------------------------------------------------------------------------------------------------------------- 766 \subsection{Handling convection with turbulent closure schemes (\forcode{ln_zdf\{tke,gls,osm\} = .true.})} 768 767 \label{subsec:ZDF_tcs} 769 768 770 769 771 770 The turbulent closure schemes presented in \autoref{subsec:ZDF_tke}, \autoref{subsec:ZDF_gls} and 772 \autoref{subsec:ZDF_osm} (\ie \np{ln\_zdftke} or \np{ln\_zdfgls} or \np{ln\_zdfosm} defined) deal, in theory,771 \autoref{subsec:ZDF_osm} (\ie\ \np{ln\_zdftke} or \np{ln\_zdfgls} or \np{ln\_zdfosm} defined) deal, in theory, 773 772 with statically unstable density profiles. 774 773 In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in 775 \autoref{eq:zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative. 776 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also of the four neighboring values at 774 \autoref{eq:zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative. 775 It results in large values of $A_T^{vT}$ and $A_T^{vT}$, and also of the four neighboring values at 777 776 velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$). 778 777 These large values restore the static stability of the water column in a way similar to that of … … 782 781 because the mixing length scale is bounded by the distance to the sea surface. 783 782 It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme, 784 \ie setting the \np{ln\_zdfnpc} namelist parameter to true and783 \ie\ setting the \np{ln\_zdfnpc} namelist parameter to true and 785 784 defining the turbulent closure (\np{ln\_zdftke} or \np{ln\_zdfgls} = \forcode{.true.}) all together. 786 785 … … 804 803 805 804 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter 806 \np{ln\_zdfddm} in \n gn{namzdf}.805 \np{ln\_zdfddm} in \nam{zdf}. 807 806 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. 808 807 The former condition leads to salt fingering and the latter to diffusive convection. 809 808 Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the ocean. 810 \citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that 809 \citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that 811 810 it leads to relatively minor changes in circulation but exerts significant regional influences on 812 811 temperature and salinity. 813 812 814 813 815 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 814 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 816 815 \begin{align*} 817 816 % \label{eq:zdfddm_Kz} … … 834 833 \end{cases} 835 834 \\ \label{eq:zdfddm_f_T} 836 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 835 A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 837 836 \end{align} 838 837 … … 860 859 861 860 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by 862 Federov (1988) is used: 861 Federov (1988) is used: 863 862 \begin{align} 864 863 % \label{eq:zdfddm_d} … … 900 899 %-------------------------------------------------------------------------------------------------------------- 901 900 902 Options to define the top and bottom friction are defined through the \n gn{namdrg} namelist variables.901 Options to define the top and bottom friction are defined through the \nam{drg} namelist variables. 903 902 The bottom friction represents the friction generated by the bathymetry. 904 903 The top friction represents the friction generated by the ice shelf/ocean interface. 905 As the friction processes at the top and the bottom are treated in and identical way, 904 As the friction processes at the top and the bottom are treated in and identical way, 906 905 the description below considers mostly the bottom friction case, if not stated otherwise. 907 906 … … 921 920 one needs a vertical diffusion coefficient $A^{vm} = 0.125$~m$^2$s$^{-1}$ 922 921 (for a Coriolis frequency $f = 10^{-4}$~m$^2$s$^{-1}$). 923 With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. 922 With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. 924 923 When the vertical mixing coefficient is this small, using a flux condition is equivalent to 925 924 entering the viscous forces (either wind stress or bottom friction) as a body force over the depth of the top or … … 951 950 - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b 952 951 \end{equation} 953 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 954 Note than from \NEMO 4.0, drag coefficients are only computed at cell centers (\ieat T-points) and refer to as $c_b^T$ in the following. These are then linearly interpolated in space to get $c_b^\textbf{U}$ at velocity points.952 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 953 Note than from \NEMO\ 4.0, drag coefficients are only computed at cell centers (\ie\ at T-points) and refer to as $c_b^T$ in the following. These are then linearly interpolated in space to get $c_b^\textbf{U}$ at velocity points. 955 954 956 955 % ------------------------------------------------------------------------------------------------------------- … … 962 961 963 962 The linear friction parameterisation (including the special case of a free-slip condition) assumes that 964 the friction is proportional to the interior velocity (\ie the velocity of the first/last model level):963 the friction is proportional to the interior velocity (\ie\ the velocity of the first/last model level): 965 964 \[ 966 965 % \label{eq:zdfbfr_linear} … … 968 967 \] 969 968 where $r$ is a friction coefficient expressed in $m s^{-1}$. 970 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 969 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 971 970 and setting $r = H / \tau$, where $H$ is the ocean depth. 972 971 Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}. … … 979 978 It can be changed by specifying \np{rn\_Uc0} (namelist parameter). 980 979 981 For the linear friction case the drag coefficient used in the general expression \autoref{eq:zdfbfr_bdef} is: 980 For the linear friction case the drag coefficient used in the general expression \autoref{eq:zdfbfr_bdef} is: 982 981 \[ 983 982 % \label{eq:zdfbfr_linbfr_b} … … 1001 1000 \label{subsec:ZDF_drg_nonlinear} 1002 1001 1003 The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic: 1002 The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic: 1004 1003 \[ 1005 1004 % \label{eq:zdfdrg_nonlinear} … … 1089 1088 To ensure that the friction cannot reverse the direction of flow it is necessary to have: 1090 1089 \[ 1091 |\Delta u| < \;|u| 1090 |\Delta u| < \;|u| 1092 1091 \] 1093 1092 \noindent which, using \autoref{eq:Eqn_drgstab}, gives: … … 1104 1103 For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then $e_{3u}$ should be greater than 3.6 m. 1105 1104 For most applications, with physically sensible parameters these restrictions should not be of concern. 1106 But caution may be necessary if attempts are made to locally enhance the bottom friction parameters. 1105 But caution may be necessary if attempts are made to locally enhance the bottom friction parameters. 1107 1106 To ensure stability limits are imposed on the top/bottom friction coefficients both 1108 1107 during initialisation and at each time step. … … 1127 1126 An optional implicit form of bottom friction has been implemented to improve model stability. 1128 1127 We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting. 1129 This option can be invoked by setting \np{ln\_drgimp} to \forcode{.true.} in the \ textit{namdrg} namelist.1130 %This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \ textit{namzdf} namelist.1131 1132 This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step: 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. 1130 1131 This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step: 1133 1132 1134 1133 At the top (below an ice shelf cavity): … … 1146 1145 \] 1147 1146 1148 where $t$ and $b$ refers to top and bottom layers respectively. 1147 where $t$ and $b$ refers to top and bottom layers respectively. 1149 1148 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit. 1150 1149 … … 1156 1155 \label{subsec:ZDF_drg_ts} 1157 1156 1158 With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \np{ln\_drgimp}\forcode{= .false.} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie \forcode{ln_bt_fw = .false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \np{ln\_drgimp}\forcode{= .true.}, stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions.1159 1160 The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO is as follows: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. 1158 1159 The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO\ is as follows: 1161 1160 \begin{enumerate} 1162 1161 \item To extend the stability of the barotropic sub-stepping, bottom stresses are refreshed at each sub-iteration. The baroclinic part of the flow entering the stresses is frozen at the initial time of the barotropic iteration. In case of non-linear friction, the drag coefficient is also constant. 1163 1162 \item In case of an implicit drag, specific computations are performed in \mdl{dynzdf} which renders the overall scheme mixed explicit/implicit: the barotropic components of 3d velocities are removed before seeking for the implicit vertical diffusion result. Top/bottom stresses due to the barotropic components are explicitly accounted for thanks to the updated values of barotropic velocities. Then the implicit solution of 3d velocities is obtained. Lastly, the residual barotropic component is replaced by the time split estimate. 1164 \end{enumerate} 1165 1166 Note that other strategies are possible, like considering vertical diffusion step in advance, \ie prior barotropic integration.1163 \end{enumerate} 1164 1165 Note that other strategies are possible, like considering vertical diffusion step in advance, \ie\ prior barotropic integration. 1167 1166 1168 1167 … … 1182 1181 the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}. 1183 1182 A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed, 1184 and the resulting diffusivity is obtained as 1183 and the resulting diffusivity is obtained as 1185 1184 \[ 1186 1185 % \label{eq:Kwave} … … 1198 1197 the mixing efficiency is constant. 1199 1198 1200 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 1201 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice. 1199 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. 1202 1201 This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14}, 1203 1202 is implemented as in \cite{de-lavergne.madec.ea_JPO16}. … … 1211 1210 F_{pyc}(i,j,k) &\propto N^{n_p}\\ 1212 1211 F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } 1213 \end{align*} 1212 \end{align*} 1214 1213 In the above formula, $h_{ab}$ denotes the height above bottom, 1215 1214 $h_{wkb}$ denotes the WKB-stretched height above bottom, defined by … … 1217 1216 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz' } \; , 1218 1217 \] 1219 The $n_p$ parameter (given by \np{nn\_zpyc} in \n gn{namzdf\_iwm} namelist)1218 The $n_p$ parameter (given by \np{nn\_zpyc} in \nam{zdf\_iwm} namelist) 1220 1219 controls the stratification-dependence of the pycnocline-intensified dissipation. 1221 1220 It can take values of $1$ (recommended) or $2$. … … 1229 1228 1230 1229 % ================================================================ 1231 % surface wave-induced mixing 1230 % surface wave-induced mixing 1232 1231 % ================================================================ 1233 1232 \section[Surface wave-induced mixing (\forcode{ln_zdfswm = .true.})] … … 1237 1236 Surface waves produce an enhanced mixing through wave-turbulence interaction. 1238 1237 In addition to breaking waves induced turbulence (\autoref{subsec:ZDF_tke}), 1239 the influence of non-breaking waves can be accounted introducing 1238 the influence of non-breaking waves can be accounted introducing 1240 1239 wave-induced viscosity and diffusivity as a function of the wave number spectrum. 1241 1240 Following \citet{qiao.yuan.ea_OD10}, a formulation of wave-induced mixing coefficient … … 1247 1246 \end{equation} 1248 1247 1249 Where $B_{v}$ is the wave-induced mixing coefficient, $A$ is the wave amplitude, 1250 ${U}_{st}$ is the Stokes Drift velocity, $k$ is the wave number and $\alpha$ 1251 is a constant which should be determined by observations or 1248 Where $B_{v}$ is the wave-induced mixing coefficient, $A$ is the wave amplitude, 1249 ${U}_{st}$ is the Stokes Drift velocity, $k$ is the wave number and $\alpha$ 1250 is a constant which should be determined by observations or 1252 1251 numerical experiments and is set to be 1. 1253 1252 1254 The coefficient $B_{v}$ is then directly added to the vertical viscosity 1253 The coefficient $B_{v}$ is then directly added to the vertical viscosity 1255 1254 and diffusivity coefficients. 1256 1255 1257 1256 In order to account for this contribution set: \forcode{ln_zdfswm = .true.}, 1258 1257 then wave interaction has to be activated through \forcode{ln_wave = .true.}, 1259 the Stokes Drift can be evaluated by setting \forcode{ln_sdw = .true.} 1258 the Stokes Drift can be evaluated by setting \forcode{ln_sdw = .true.} 1260 1259 (see \autoref{subsec:SBC_wave_sdw}) 1261 1260 and the needed wave fields can be provided either in forcing or coupled mode … … 1269 1268 \label{subsec:ZDF_aimp} 1270 1269 1271 This refers to \citep{shchepetkin_OM15}. 1272 1273 TBC ... 1274 1275 1270 The adaptive-implicit vertical advection option in NEMO is based on the work of 1271 \citep{shchepetkin_OM15}. In common with most ocean models, the timestep used with NEMO 1272 needs to satisfy multiple criteria associated with different physical processes in order 1273 to maintain numerical stability. \citep{shchepetkin_OM15} pointed out that the vertical 1274 CFL criterion is commonly the most limiting. \citep{lemarie.debreu.ea_OM15} examined the 1275 constraints for a range of time and space discretizations and provide the CFL stability 1276 criteria for a range of advection schemes. The values for the Leap-Frog with Robert 1277 asselin filter time-stepping (as used in NEMO) are reproduced in 1278 \autoref{tab:zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these 1279 restrictions but at the cost of large dispersive errors and, possibly, large numerical 1280 viscosity. The adaptive-implicit vertical advection option provides a targetted use of the 1281 implicit scheme only when and where potential breaches of the vertical CFL condition 1282 occur. In many practical applications these events may occur remote from the main area of 1283 interest or due to short-lived conditions such that the extra numerical diffusion or 1284 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 whilst 1286 retaining the accuracy of the high order explicit schemes over most of the domain. 1287 1288 \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} 1303 \end{table} 1304 1305 In particular, the advection scheme remains explicit everywhere except where and when 1306 local vertical velocities exceed a threshold set just below the explicit stability limit. 1307 Once the threshold is reached a tapered transition towards an implicit scheme is used by 1308 partitioning the vertical velocity into a part that can be treated explicitly and any 1309 excess that must be treated implicitly. The partitioning is achieved via a Courant-number 1310 dependent weighting algorithm as described in \citep{shchepetkin_OM15}. 1311 1312 The local cell Courant number ($Cu$) used for this partitioning is: 1313 1314 \begin{equation} 1315 \label{eq:Eqn_zad_Aimp_Courant} 1316 \begin{split} 1317 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 ] \\ 1318 &\phantom{=} +\big [ \texttt{Max}(e_{{2_u}ij}e^n_{{3_{u}}ijk}u^n_{ijk},0.0) - \texttt{Min}(e_{{2_u}i-1j}e^n_{{3_{u}}i-1jk}u^n_{i-1jk},0.0) \big ] 1319 \big / e_{{1_t}ij}e_{{2_t}ij} \\ 1320 &\phantom{=} +\big [ \texttt{Max}(e_{{1_v}ij}e^n_{{3_{v}}ijk}v^n_{ijk},0.0) - \texttt{Min}(e_{{1_v}ij-1}e^n_{{3_{v}}ij-1k}v^n_{ij-1k},0.0) \big ] 1321 \big / e_{{1_t}ij}e_{{2_t}ij} \bigg ) \\ 1322 \end{split} 1323 \end{equation} 1324 1325 \noindent and the tapering algorithm follows \citep{shchepetkin_OM15} as: 1326 1327 \begin{align} 1328 \label{eq:Eqn_zad_Aimp_partition} 1329 Cu_{min} &= 0.15 \nonumber \\ 1330 Cu_{max} &= 0.3 \nonumber \\ 1331 Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ 1332 Fcu &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ 1333 C\kern-0.14em f &= 1334 \begin{cases} 1335 0.0 &\text{if $Cu \leq Cu_{min}$} \\ 1336 (Cu - Cu_{min})^2 / (Fcu + (Cu - Cu_{min})^2) &\text{else if $Cu < Cu_{cut}$} \\ 1337 (Cu - Cu_{max}) / Cu &\text{else} 1338 \end{cases} 1339 \end{align} 1340 1341 \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} 1350 \end{figure} 1351 1352 \noindent The partitioning coefficient is used to determine the part of the vertical 1353 velocity that must be handled implicitly ($w_i$) and to subtract this from the total 1354 vertical velocity ($w_n$) to leave that which can continue to be handled explicitly: 1355 1356 \begin{align} 1357 \label{eq:Eqn_zad_Aimp_partition2} 1358 w_{i_{ijk}} &= C\kern-0.14em f_{ijk} w_{n_{ijk}} \nonumber \\ 1359 w_{n_{ijk}} &= (1-C\kern-0.14em f_{ijk}) w_{n_{ijk}} 1360 \end{align} 1361 1362 \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: 1364 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 values 1366 the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly 1367 implicit' is 0.45 which is just below the stability limited given in 1368 \autoref{tab:zad_Aimp_CFLcrit} for a 3rd order scheme. 1369 1370 The $w_i$ component is added to the implicit solvers for the vertical mixing in 1371 \mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}. This is 1372 sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further 1373 intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}). 1374 For these schemes the implicit upstream fluxes must be added to both the monotonic guess 1375 and to the higher order solution when calculating the antidiffusive fluxes. The implicit 1376 vertical fluxes are then removed since they are added by the implicit solver later on. 1377 1378 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 1379 used in a wide range of simulations. The following test simulation, however, does illustrate 1380 the potential benefits and will hopefully encourage further testing and feedback from users: 1381 1382 \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} 1391 \end{figure} 1392 1393 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} 1394 The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case} 1395 provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case 1396 by only a few extra physics choices namely: 1397 1398 \begin{verbatim} 1399 ln_dynldf_OFF = .false. 1400 ln_dynldf_lap = .true. 1401 ln_dynldf_hor = .true. 1402 ln_zdfnpc = .true. 1403 ln_traadv_fct = .true. 1404 nn_fct_h = 2 1405 nn_fct_v = 2 1406 \end{verbatim} 1407 1408 \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-implicit 1410 vertical velocity is illustrated in \autoref{fig:zad_Aimp_overflow_frames}. The mass of 1411 cold water, initially sitting on the shelf, moves down the slope and forms a 1412 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 hours 1414 with excessively high horizontal velocities. This time-scale corresponds with the time the 1415 plume reaches the steepest part of the topography and, although detected as a horizontal 1416 CFL breach, the instability originates from a breach of the vertical CFL limit. This is a good 1417 candidate, therefore, for use of the adaptive-implicit vertical advection scheme. 1418 1419 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 equivalent 1421 frames from the base run). In this simple example the use of the adaptive-implicit 1422 vertcal advection scheme has enabled a 12x increase in the model timestep without 1423 significantly altering the solution (although at this extreme the plume is more diffuse 1424 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 was 1426 close enough to instability to trigger the scheme despite completing successfully. 1427 To assist in diagnosing how active the scheme is, in both location and time, the 3D 1428 implicit and explicit components of the vertical velocity are available via XIOS as 1429 \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 of 1431 the schemes activity the global maximum values of the absolute implicit component 1432 of the vertical velocity and the partitioning coefficient are written to the netCDF 1433 version of the run statistics file (\texttt{run.stat.nc}) if this is active (see 1434 \autoref{sec:MISC_opt} for activation details). 1435 1436 \autoref{fig:zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for 1437 the various overflow tests. Note that the adaptive-implicit vertical advection scheme is 1438 active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the 1439 test case is close to stability limits even with this value. At the larger timesteps, the 1440 vertical velocity is treated mostly implicitly at some location throughout the run. The 1441 oscillatory nature of this measure appears to be linked to the progress of the plume front 1442 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 the 1444 maximum have been overlaid for the base run case. 1445 1446 \medskip 1447 \noindent Only limited tests have been performed in more realistic configurations. In the 1448 ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes 1449 restartability and reproducibility tests but it is unable to improve the model's stability 1450 enough to allow an increase in the model time-step. A view of the time-series of maximum 1451 partitioning coefficient (not shown here) suggests that the default time-step of 5400s is 1452 already pushing at stability limits, especially in the initial start-up phase. The 1453 time-series does not, however, exhibit any of the 'cuspiness' found with the overflow 1454 tests. 1455 1456 \medskip 1457 \noindent A short test with an eORCA1 configuration promises more since a test using a 1458 time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the 1459 time-step is limited to 2700s without. 1460 1461 \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 only 1468 the run with the shortest timestep is able to run to completion. Note also that the colour-scale has been 1469 chosen to confirm that temperatures remain within the original range of 10$^o$ to 20$^o$. 1470 } 1471 \end{center} 1472 \end{figure} 1473 1474 \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} 1484 \end{figure} 1485 1486 \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} 1495 \end{figure} 1276 1496 1277 1497 % ================================================================
Note: See TracChangeset
for help on using the changeset viewer.
