Changeset 11435 for NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex
- Timestamp:
- 2019-08-14T14:45:08+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex
r11407 r11435 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
Note: See TracChangeset
for help on using the changeset viewer.