Changeset 11213


Ignore:
Timestamp:
2019-07-04T11:29:44+02:00 (9 months ago)
Author:
jchanut
Message:

#2216: Add Romain's update of ZDF chapter

Location:
NEMO/trunk/doc/latex/NEMO
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/doc/latex/NEMO/main/bibliography.bib

    r11170 r11213  
    23202320} 
    23212321 
     2322@article{         reffray.guillaume.ea_GMD15, 
     2323  title         = "Modelling turbulent vertical mixing sensitivity using 
     2324                   a 1-D version of NEMO", 
     2325  pages         = "69--86", 
     2326  journal       = "Geoscientific Model Development", 
     2327  volume        = "8", 
     2328  number        = "1", 
     2329  author        = "G. Reffray. and R. Bourdalle-Badie and C. Calone", 
     2330  year          = "2015", 
     2331  month         = "jan", 
     2332  publisher     = "Copernicus GmbH", 
     2333  issn          = "1991-9603", 
     2334  doi           = "10.5194/gmd-8-69-2015" 
     2335} 
     2336 
    23222337@book{            richtmyer.morton_bk67, 
    23232338  pages         = "405", 
     
    24842499} 
    24852500 
     2501@article{         shchepetkin_OM15, 
     2502  title         = "An Adaptive, Courant-number-dependent implicit 
     2503       scheme for vertical advection in oceanic modeling", 
     2504  pages         = "38--69", 
     2505  journal       = "Ocean Modelling", 
     2506  volume        = "91", 
     2507  author        = "A. F. Shchepetkin", 
     2508  year          = "2015", 
     2509  month         = "july", 
     2510  publisher     = "Elsevier BV", 
     2511  issn          = "1463-5003", 
     2512  doi           = "10.1016/j.ocemod.2015.03.006" 
     2513} 
     2514 
    24862515@article{         siddorn.furner_OM13, 
    24872516  title         = "An analytical stretching function that combines the best 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex

    r11179 r11213  
    2525At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}), 
    2626while at the bottom they are set to zero for heat and salt, 
    27 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie \key{trabbl} defined, 
     27unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie \np{ln\_trabbc} defined, 
    2828see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum 
    29 (see \autoref{sec:ZDF_bfr}). 
     29(see \autoref{sec:ZDF_drg}).  
    3030 
    3131In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and 
     
    3333respectively (see \autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}). 
    3434These coefficients can be assumed to be either constant, or a function of the local Richardson number, 
    35 or computed from a turbulent closure model (either TKE or GLS formulation). 
    36 The computation of these coefficients is initialized in the \mdl{zdfini} module and performed in 
    37 the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} modules. 
     35or computed from a turbulent closure model (either TKE or GLS or OSMOSIS formulation). 
     36The computation of these coefficients is initialized in the \mdl{zdfphy} module and performed in 
     37the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules. 
    3838The trends due to the vertical momentum and tracer diffusion, including the surface forcing, 
    3939are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.  
    40 These trends can be computed using either a forward time stepping scheme 
    41 (namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme 
    42 (\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing coefficients, 
    43 and thus of the formulation used (see \autoref{chap:STP}). 
    44  
    45 % ------------------------------------------------------------------------------------------------------------- 
    46 %        Constant  
    47 % ------------------------------------------------------------------------------------------------------------- 
    48 \subsection[Constant (\texttt{\textbf{key\_zdfcst}})] 
    49 {Constant (\protect\key{zdfcst})} 
    50 \label{subsec:ZDF_cst} 
    51 %--------------------------------------------namzdf--------------------------------------------------------- 
     40%These trends can be computed using either a forward time stepping scheme 
     41%(namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme 
     42%(\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing coefficients, 
     43%and thus of the formulation used (see \autoref{chap:STP}). 
     44 
     45%--------------------------------------------namzdf-------------------------------------------------------- 
    5246 
    5347\nlst{namzdf} 
    5448%-------------------------------------------------------------------------------------------------------------- 
    5549 
     50% ------------------------------------------------------------------------------------------------------------- 
     51%        Constant  
     52% ------------------------------------------------------------------------------------------------------------- 
     53\subsection[Constant (\forcode{ln_zdfcst = .true.})] 
     54{Constant (\protect\np{ln\_zdfcst}\forcode{ = .true.})} 
     55\label{subsec:ZDF_cst} 
     56 
    5657Options are defined through the \ngn{namzdf} namelist variables. 
    57 When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 
     58When \np{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to 
    5859constant values over the whole ocean. 
    5960This is the crudest way to define the vertical ocean physics. 
    60 It is recommended that this option is only used in process studies, not in basin scale simulations. 
     61It is recommended to use this option only in process studies, not in basin scale simulations. 
    6162Typical values used in this case are: 
    6263\begin{align*} 
     
    7374%        Richardson Number Dependent 
    7475% ------------------------------------------------------------------------------------------------------------- 
    75 \subsection[Richardson number dependent (\texttt{\textbf{key\_zdfric}})] 
    76 {Richardson number dependent (\protect\key{zdfric})} 
     76\subsection[Richardson number dependent (\forcode{ln_zdfric = .true.})] 
     77{Richardson number dependent (\protect\np{ln\_zdfric}\forcode{ = .true.})} 
    7778\label{subsec:ZDF_ric} 
    7879 
     
    8283%-------------------------------------------------------------------------------------------------------------- 
    8384 
    84 When \key{zdfric} is defined, a local Richardson number dependent formulation for the vertical momentum and 
     85When \np{ln\_zdfric}\forcode{ = .true.}, a local Richardson number dependent formulation for the vertical momentum and 
    8586tracer eddy coefficients is set through the \ngn{namzdf\_ric} namelist variables. 
    8687The vertical mixing coefficients are diagnosed from the large scale variables computed by the model.  
     
    131132%        TKE Turbulent Closure Scheme  
    132133% ------------------------------------------------------------------------------------------------------------- 
    133 \subsection[TKE turbulent closure scheme (\texttt{\textbf{key\_zdftke}})] 
    134 {TKE turbulent closure scheme (\protect\key{zdftke})} 
     134\subsection[TKE turbulent closure scheme (\forcode{ln_zdftke = .true.})] 
     135{TKE turbulent closure scheme (\protect\np{ln\_zdftke}\forcode{ = .true.})} 
    135136\label{subsec:ZDF_tke} 
    136  
    137137%--------------------------------------------namzdf_tke-------------------------------------------------- 
    138138 
     
    183183  \end{cases} 
    184184\end{align*} 
    185 Options are defined through the  \ngn{namzdfy\_tke} namelist variables. 
    186185The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable. 
    187186 
     
    200199too weak vertical diffusion. 
    201200They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and 
    202 \np{rn\_avt0} (namzdf namelist, see \autoref{subsec:ZDF_cst}). 
     201\np{rn\_avt0} (\ngn{namzdf} namelist, see \autoref{subsec:ZDF_cst}). 
    203202 
    204203\subsubsection{Turbulent length scale} 
     
    218217it makes no sense for locally unstable stratification and the computation no longer uses all 
    219218the information contained in the vertical density profile. 
    220 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{ = 2..3} cases, 
     219To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{ = 2, 3} cases, 
    221220which add an extra assumption concerning the vertical gradient of the computed length scale. 
    222221So, the length scales are first evaluated as in \autoref{eq:tke_mxl0_1} and then bounded such that: 
     
    285284The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and 
    286285air-sea drag coefficient.  
    287 The latter concerns the bulk formulea and is not discussed here.  
     286The latter concerns the bulk formulae and is not discussed here.  
    288287 
    289288Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is : 
     
    306305with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds  
    307306to $\alpha_{CB} = 100$. 
    308 Further setting  \np{ln\_mxl0} to true applies \autoref{eq:ZDF_Lsbc} as surface boundary condition on length scale, 
     307Further setting  \np{ln\_mxl0=.true.},  applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale, 
    309308with $\beta$ hard coded to the Stacey's value. 
    310 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on 
     309Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the  
    311310surface $\bar{e}$ value. 
    312311 
     
    325324\citep{axell_JGR02} for a $k-\epsilon$ turbulent closure. 
    326325The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in 
    327 an extra source terms of TKE, $P_{LC}$. 
     326an extra source term of TKE, $P_{LC}$. 
    328327The presence of $P_{LC}$ in \autoref{eq:zdftke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to 
    329 \forcode{.true.} in the namtke namelist. 
     328\forcode{.true.} in the \ngn{namzdf\_tke} namelist. 
    330329  
    331330By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}), 
     
    359358 
    360359The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations: 
    361 $H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by 
     360$H_{LC}$ is the depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by 
    362361converting its kinetic energy to potential energy, according to  
    363362\[ 
     
    386385\end{equation} 
    387386where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that 
    388 penetrate in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of 
     387penetrates in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of 
    389388the penetration, and $f_i$ is the ice concentration 
    390 (no penetration if $f_i=1$, that is 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). 
    391390The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter. 
    392391The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}\forcode{ = 0}) or 
     
    394393(\np{nn\_etau}\forcode{ = 1}).  
    395394 
    396 Note that two other option existe, \np{nn\_etau}\forcode{ = 2..3}. 
     395Note that two other option exist, \np{nn\_etau}\forcode{ = 2, 3}. 
    397396They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer, 
    398 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrate the ocean.  
     397or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean.  
    399398Those two options are obsolescent features introduced for test purposes. 
    400399They will be removed in the next release.  
     400 
     401% This should be explain better below what this rn_eice parameter is meant for: 
     402In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn\_eice} namelist parameter. 
     403This parameter varies from \forcode{0} for no effect to \forcode{4} to suppress the TKE input into the ocean when Sea Ice concentration 
     404is greater than 25\%.  
    401405 
    402406% from Burchard et al OM 2008 :  
     
    409413 
    410414% ------------------------------------------------------------------------------------------------------------- 
    411 %        TKE discretization considerations 
    412 % ------------------------------------------------------------------------------------------------------------- 
    413 \subsection[TKE discretization considerations (\texttt{\textbf{key\_zdftke}})] 
    414 {TKE discretization considerations (\protect\key{zdftke})} 
    415 \label{subsec:ZDF_tke_ene} 
    416  
    417 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    418 \begin{figure}[!t] 
    419   \begin{center} 
    420     \includegraphics[width=\textwidth]{Fig_ZDF_TKE_time_scheme} 
    421     \caption{ 
    422       \protect\label{fig:TKE_time_scheme} 
    423       Illustration of the TKE time integration and its links to the momentum and tracer time integration. 
    424     } 
    425   \end{center}   
    426 \end{figure} 
    427 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    428  
    429 The production of turbulence by vertical shear (the first term of the right hand side of 
    430 \autoref{eq:zdftke_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion 
    431 (first line in \autoref{eq:PE_zdf}). 
    432 To do so a special care have to be taken for both the time and space discretization of 
    433 the TKE equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}. 
    434  
    435 Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows how 
    436 the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with 
    437 the one-level forward time stepping of TKE equation. 
    438 With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to 
    439 the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and 
    440 summing the result vertically:    
    441 \begin{equation} 
    442   \label{eq:energ1} 
    443   \begin{split} 
    444     \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\ 
    445     &= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta} 
    446     - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz } 
    447   \end{split} 
    448 \end{equation} 
    449 Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$, 
    450 known at time $t$ (\autoref{fig:TKE_time_scheme}), as it is required when using the TKE scheme 
    451 (see \autoref{sec:STP_forward_imp}). 
    452 The first term of the right hand side of \autoref{eq:energ1} represents the kinetic energy transfer at 
    453 the surface (atmospheric forcing) and at the bottom (friction effect). 
    454 The second term is always negative. 
    455 It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$. 
    456 \autoref{eq:energ1} implies that, to be energetically consistent, 
    457 the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 
    458 ${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ 
    459 (and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$). 
    460  
    461 A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 
    462 (second term of the right hand side of \autoref{eq:zdftke_e}). 
    463 This term must balance the input of potential energy resulting from vertical mixing. 
    464 The rate of change of potential energy (in 1D for the demonstration) due vertical mixing is obtained by 
    465 multiplying vertical density diffusion tendency by $g\,z$ and and summing the result vertically: 
    466 \begin{equation} 
    467   \label{eq:energ2} 
    468   \begin{split} 
    469     \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\ 
    470     &= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta} 
    471     - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\ 
    472     &= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta} 
    473     + \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  } 
    474   \end{split} 
    475 \end{equation} 
    476 where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.  
    477 The first term of the right hand side of \autoref{eq:energ2} is always zero because 
    478 there is no diffusive flux through the ocean surface and bottom). 
    479 The second term is minus the destruction rate of  $\bar{e}$ due to stratification. 
    480 Therefore \autoref{eq:energ1} implies that, to be energetically consistent, 
    481 the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e}, the TKE equation. 
    482  
    483 Let us now address the space discretization issue. 
    484 The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in 
    485 the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:cell}). 
    486 A space averaging is thus required to obtain the shear TKE production term. 
    487 By redoing the \autoref{eq:energ1} in the 3D case, it can be shown that the product of eddy coefficient by 
    488 the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. 
    489 Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into account. 
    490  
    491 The above energetic considerations leads to the following final discrete form for the TKE equation: 
    492 \begin{equation} 
    493   \label{eq:zdftke_ene} 
    494   \begin{split} 
    495     \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv 
    496     \Biggl\{ \Biggr. 
    497     &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} } 
    498         \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\ 
    499     +&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} } 
    500         \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j} 
    501     \Biggr. \Biggr\}   \\ 
    502     % 
    503     - &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\ 
    504     % 
    505     +&\frac{1}{{e_3w}^{t+\rdt}}  \;\delta_{k+1/2} \left[   {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}}   \right]   \\ 
    506     % 
    507     - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt} 
    508   \end{split} 
    509 \end{equation} 
    510 where the last two terms in \autoref{eq:zdftke_ene} (vertical diffusion and Kolmogorov dissipation) 
    511 are time stepped using a backward scheme (see\autoref{sec:STP_forward_imp}). 
    512 Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible. 
    513 The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as 
    514 they all appear in the right hand side of \autoref{eq:zdftke_ene}. 
    515 For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.  
    516  
    517 % ------------------------------------------------------------------------------------------------------------- 
    518415%        GLS Generic Length Scale Scheme  
    519416% ------------------------------------------------------------------------------------------------------------- 
    520 \subsection[GLS: Generic Length Scale (\texttt{\textbf{key\_zdfgls}})] 
    521 {GLS: Generic Length Scale (\protect\key{zdfgls})} 
     417\subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls = .true.})] 
     418{GLS: Generic Length Scale (\protect\np{ln\_zdfgls}\forcode{ = .true.})} 
    522419\label{subsec:ZDF_gls} 
    523420 
     
    573470The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of 
    574471the choice of the turbulence model. 
    575 Four different turbulent models are pre-defined (Tab.\autoref{tab:GLS}). 
     472Four different turbulent models are pre-defined (\autoref{tab:GLS}). 
    576473They are made available through the \np{nn\_clo} namelist parameter.  
    577474 
     
    600497      \protect\label{tab:GLS} 
    601498      Set of predefined GLS parameters, or equivalently predefined turbulence models available with 
    602       \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls}. 
     499      \protect\np{ln\_zdfgls}\forcode{ = .true.} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls}. 
    603500    } 
    604501  \end{center} 
     
    607504 
    608505In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of 
    609 the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length) value near physical boundaries 
     506the mixing length towards $\kappa z_b$ ($\kappa$ is the Von Karman constant and $z_b$ the rugosity length scale) value near physical boundaries 
    610507(logarithmic boundary layer law). 
    611508$C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88}, 
    612509or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01} 
    613 (\np{nn\_stab\_func}\forcode{ = 0..3}, resp.).  
    614 The value of $C_{0\mu}$ depends of the choice of the stability function. 
     510(\np{nn\_stab\_func}\forcode{ = 0, 3}, resp.).  
     511The value of $C_{0\mu}$ depends on the choice of the stability function. 
    615512 
    616513The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or 
     
    633530The time and space discretization of the GLS equations follows the same energetic consideration as for 
    634531the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}. 
    635 Examples of performance of the 4 turbulent closure scheme can be found in \citet{warner.sherwood.ea_OM05}. 
     532Evaluation 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 
    636535 
    637536% ------------------------------------------------------------------------------------------------------------- 
    638537%        OSM OSMOSIS BL Scheme  
    639538% ------------------------------------------------------------------------------------------------------------- 
    640 \subsection[OSM: OSMosis boundary layer scheme (\texttt{\textbf{key\_zdfosm}})] 
    641 {OSM: OSMosis boundary layer scheme (\protect\key{zdfosm})} 
     539\subsection[OSM: OSMosis boundary layer scheme (\forcode{ln_zdfosm = .true.})] 
     540{OSM: OSMosis boundary layer scheme (\protect\np{ln\_zdfosm}\forcode{ = .true.})} 
    642541\label{subsec:ZDF_osm} 
    643  
    644542%--------------------------------------------namzdf_osm--------------------------------------------------------- 
    645543 
     
    649547The OSMOSIS turbulent closure scheme is based on......   TBC 
    650548 
     549% ------------------------------------------------------------------------------------------------------------- 
     550%        TKE and GLS discretization considerations 
     551% ------------------------------------------------------------------------------------------------------------- 
     552\subsection[ Discrete energy conservation for TKE and GLS schemes] 
     553{Discrete energy conservation for TKE and GLS schemes} 
     554\label{subsec:ZDF_tke_ene} 
     555 
     556%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     557\begin{figure}[!t] 
     558  \begin{center} 
     559    \includegraphics[width=\textwidth]{Fig_ZDF_TKE_time_scheme} 
     560    \caption{ 
     561      \protect\label{fig:TKE_time_scheme} 
     562      Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and its links to the momentum and tracer time integration. 
     563    } 
     564  \end{center}   
     565\end{figure} 
     566%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     567 
     568The production of turbulence by vertical shear (the first term of the right hand side of 
     569\autoref{eq:zdftke_e}) and  \autoref{eq:zdfgls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion 
     570(first line in \autoref{eq:PE_zdf}). 
     571To do so a special care has to be taken for both the time and space discretization of 
     572the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}. 
     573 
     574Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows how 
     575the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with 
     576the one-level forward time stepping of the equation for $\bar{e}$. 
     577With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to 
     578the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and 
     579summing the result vertically:    
     580\begin{equation} 
     581  \label{eq:energ1} 
     582  \begin{split} 
     583    \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\ 
     584    &= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta} 
     585    - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz } 
     586  \end{split} 
     587\end{equation} 
     588Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$, 
     589known at time $t$ (\autoref{fig:TKE_time_scheme}), as it is required when using the TKE scheme 
     590(see \autoref{sec:STP_forward_imp}). 
     591The first term of the right hand side of \autoref{eq:energ1} represents the kinetic energy transfer at 
     592the surface (atmospheric forcing) and at the bottom (friction effect). 
     593The second term is always negative. 
     594It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$. 
     595\autoref{eq:energ1} implies that, to be energetically consistent, 
     596the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as 
     597${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ 
     598(and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$). 
     599 
     600A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification 
     601(second term of the right hand side of \autoref{eq:zdftke_e} and \autoref{eq:zdfgls_e}). 
     602This term must balance the input of potential energy resulting from vertical mixing. 
     603The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by 
     604multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically: 
     605\begin{equation} 
     606  \label{eq:energ2} 
     607  \begin{split} 
     608    \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\ 
     609    &= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta} 
     610    - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\ 
     611    &= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta} 
     612    + \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  } 
     613  \end{split} 
     614\end{equation} 
     615where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.  
     616The first term of the right hand side of \autoref{eq:energ2} is always zero because 
     617there is no diffusive flux through the ocean surface and bottom). 
     618The second term is minus the destruction rate of  $\bar{e}$ due to stratification. 
     619Therefore \autoref{eq:energ1} implies that, to be energetically consistent, 
     620the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e} and  \autoref{eq:zdfgls_e}. 
     621 
     622Let us now address the space discretization issue. 
     623The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in 
     624the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:cell}). 
     625A space averaging is thus required to obtain the shear TKE production term. 
     626By redoing the \autoref{eq:energ1} in the 3D case, it can be shown that the product of eddy coefficient by 
     627the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. 
     628Furthermore, the time variation of $e_3$ has be taken into account. 
     629 
     630The above energetic considerations leads to the following final discrete form for the TKE equation: 
     631\begin{equation} 
     632  \label{eq:zdftke_ene} 
     633  \begin{split} 
     634    \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv 
     635    \Biggl\{ \Biggr. 
     636    &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} } 
     637        \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\ 
     638    +&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} } 
     639        \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j} 
     640    \Biggr. \Biggr\}   \\ 
     641    % 
     642    - &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\ 
     643    % 
     644    +&\frac{1}{{e_3w}^{t+\rdt}}  \;\delta_{k+1/2} \left[   {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}}   \right]   \\ 
     645    % 
     646    - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt} 
     647  \end{split} 
     648\end{equation} 
     649where the last two terms in \autoref{eq:zdftke_ene} (vertical diffusion and Kolmogorov dissipation) 
     650are time stepped using a backward scheme (see\autoref{sec:STP_forward_imp}). 
     651Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible. 
     652%The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as 
     653%they all appear in the right hand side of \autoref{eq:zdftke_ene}. 
     654%For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.  
     655 
    651656% ================================================================ 
    652657% Convection 
     
    654659\section{Convection} 
    655660\label{sec:ZDF_conv} 
    656  
    657 %--------------------------------------------namzdf-------------------------------------------------------- 
    658  
    659 \nlst{namzdf} 
    660 %-------------------------------------------------------------------------------------------------------------- 
    661661 
    662662Static instabilities (\ie light potential densities under heavy ones) may occur at particular ocean grid points. 
     
    673673{Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})} 
    674674\label{subsec:ZDF_npc} 
    675  
    676 %--------------------------------------------namzdf-------------------------------------------------------- 
    677  
    678 \nlst{namzdf} 
    679 %-------------------------------------------------------------------------------------------------------------- 
    680675 
    681676%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    733728Two main differences have been introduced compared to the original algorithm: 
    734729$(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency  
    735 (not the the difference in potential density);  
     730(not the difference in potential density);  
    736731$(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in 
    737732the same way their temperature and salinity has been mixed. 
     
    746741\label{subsec:ZDF_evd} 
    747742 
    748 %--------------------------------------------namzdf-------------------------------------------------------- 
    749  
    750 \nlst{namzdf} 
    751 %-------------------------------------------------------------------------------------------------------------- 
    752  
    753743Options are defined through the  \ngn{namzdf} namelist variables. 
    754744The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}\forcode{ = .true.}. 
    755 In this case, the vertical eddy mixing coefficients are assigned very large values 
    756 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 
     745In this case, the vertical eddy mixing coefficients are assigned very large values  
     746in regions where the stratification is unstable 
    757747(\ie when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}. 
    758748This is done either on tracers only (\np{nn\_evdm}\forcode{ = 0}) or 
     
    766756the convective adjustment algorithm presented above when mixing both tracers and 
    767757momentum in the case of static instabilities. 
    768 It requires the use of an implicit time stepping on vertical diffusion terms 
    769 (\ie np{ln\_zdfexp}\forcode{ = .false.}). 
    770758 
    771759Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$. 
     
    776764%       Turbulent Closure Scheme  
    777765% ------------------------------------------------------------------------------------------------------------- 
    778 \subsection[Turbulent closure scheme (\texttt{\textbf{key\_zdf}}\texttt{\textbf{\{tke,gls,osm\}}})] 
    779 {Turbulent Closure Scheme (\protect\key{zdftke}, \protect\key{zdfgls} or \protect\key{zdfosm})} 
     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.})} 
    780768\label{subsec:ZDF_tcs} 
    781769 
    782 The turbulent closure scheme presented in \autoref{subsec:ZDF_tke} and \autoref{subsec:ZDF_gls} 
    783 (\key{zdftke} or \key{zdftke} is defined) in theory solves the problem of statically unstable density profiles. 
     770 
     771The 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\_zdftke} or \np{ln\_zdfosm} defined) deal, in theory,  
     773with statically unstable density profiles. 
    784774In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in 
    785775\autoref{eq:zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative.  
    786 It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also the four neighbouring $A_u^{vm} {and}\;A_v^{vm}$ 
    787 (up to $1\;m^2s^{-1}$). 
     776It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also of the four neighboring values at  
     777velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$). 
    788778These large values restore the static stability of the water column in a way similar to that of 
    789779the enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}). 
     
    793783It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme, 
    794784\ie setting the \np{ln\_zdfnpc} namelist parameter to true and 
    795 defining the turbulent closure CPP key all together. 
    796  
    797 The KPP turbulent closure scheme already includes enhanced vertical diffusion in the case of convection, 
    798 as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp}, 
    799 therefore \np{ln\_zdfevd}\forcode{ = .false.} should be used with the KPP scheme. 
     785defining the turbulent closure (\np{ln\_zdftke} or \np{ln\_zdfgls} = \forcode{.true.}) all together. 
     786 
     787The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection, 
     788%as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp}, 
     789therefore \np{ln\_zdfevd}\forcode{ = .false.} should be used with the OSMOSIS scheme. 
    800790% gm%  + one word on non local flux with KPP scheme trakpp.F90 module... 
    801791 
     
    803793% Double Diffusion Mixing 
    804794% ================================================================ 
    805 \section[Double diffusion mixing (\texttt{\textbf{key\_zdfddm}})] 
    806 {Double diffusion mixing (\protect\key{zdfddm})} 
    807 \label{sec:ZDF_ddm} 
     795\section[Double diffusion mixing (\forcode{ln_zdfddm = .true.})] 
     796{Double diffusion mixing (\protect\np{ln\_zdfddm}\forcode{ = .true.})} 
     797\label{subsec:ZDF_ddm} 
     798 
    808799 
    809800%-------------------------------------------namzdf_ddm------------------------------------------------- 
     
    812803%-------------------------------------------------------------------------------------------------------------- 
    813804 
    814 Options are defined through the  \ngn{namzdf\_ddm} namelist variables. 
     805This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter 
     806\np{ln\_zdfddm} in \ngn{namzdf}. 
    815807Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. 
    816808The former condition leads to salt fingering and the latter to diffusive convection. 
     
    819811it leads to relatively minor changes in circulation but exerts significant regional influences on 
    820812temperature and salinity. 
    821 This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the \key{zdfddm} CPP key. 
     813 
    822814 
    823815Diapycnal mixing of S and T are described by diapycnal diffusion coefficients  
     
    896888% Bottom Friction 
    897889% ================================================================ 
    898 \section[Bottom and top friction (\textit{zdfbfr.F90})] 
    899 {Bottom and top friction (\protect\mdl{zdfbfr})} 
    900 \label{sec:ZDF_bfr} 
     890 \section[Bottom and top friction (\textit{zdfdrg.F90})] 
     891 {Bottom and top friction (\protect\mdl{zdfdrg})} 
     892 \label{sec:ZDF_drg} 
    901893 
    902894%--------------------------------------------nambfr-------------------------------------------------------- 
    903895% 
    904 %\nlst{nambfr} 
     896\nlst{namdrg} 
     897\nlst{namdrg_top} 
     898\nlst{namdrg_bot} 
     899 
    905900%-------------------------------------------------------------------------------------------------------------- 
    906901 
    907 Options to define the top and bottom friction are defined through the \ngn{nambfr} namelist variables. 
     902Options to define the top and bottom friction are defined through the \ngn{namdrg} namelist variables. 
    908903The bottom friction represents the friction generated by the bathymetry. 
    909904The top friction represents the friction generated by the ice shelf/ocean interface. 
    910 As the friction processes at the top and bottom are treated in similar way, 
    911 only the bottom friction is described in detail below. 
     905As the friction processes at the top and the bottom are treated in and identical way,  
     906the description below considers mostly the bottom friction case, if not stated otherwise. 
    912907 
    913908 
     
    915910a condition on the vertical diffusive flux. 
    916911For the bottom boundary layer, one has: 
    917 \[ 
    918   % \label{eq:zdfbfr_flux} 
    919   A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 
    920 \] 
     912 \[ 
     913   % \label{eq:zdfbfr_flux} 
     914   A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 
     915 \] 
    921916where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum outside 
    922917the logarithmic turbulent boundary layer (thickness of the order of 1~m in the ocean). 
     
    932927To illustrate this, consider the equation for $u$ at $k$, the last ocean level: 
    933928\begin{equation} 
    934   \label{eq:zdfbfr_flux2} 
     929  \label{eq:zdfdrg_flux2} 
    935930  \frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}} 
    936931\end{equation} 
     
    945940 
    946941In the code, the bottom friction is imposed by adding the trend due to the bottom friction to 
    947 the general momentum trend in \mdl{dynbfr}. 
     942 the general momentum trend in \mdl{dynzdf}. 
    948943For the time-split surface pressure gradient algorithm, the momentum trend due to 
    949944the barotropic component needs to be handled separately. 
    950945For this purpose it is convenient to compute and store coefficients which can be simply combined with 
    951946bottom velocities and geometric values to provide the momentum trend due to bottom friction. 
    952 These coefficients are computed in \mdl{zdfbfr} and generally take the form $c_b^{\textbf U}$ where: 
     947 These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where: 
    953948\begin{equation} 
    954949  \label{eq:zdfbfr_bdef} 
     
    956951  - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b 
    957952\end{equation} 
    958 where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity. 
     953where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.  
     954Note 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. 
    959955 
    960956% ------------------------------------------------------------------------------------------------------------- 
    961957%       Linear Bottom Friction 
    962958% ------------------------------------------------------------------------------------------------------------- 
    963 \subsection[Linear bottom friction (\forcode{nn_botfr = [01]})] 
    964 {Linear bottom friction (\protect\np{nn\_botfr}\forcode{ = [01])}} 
    965 \label{subsec:ZDF_bfr_linear} 
    966  
    967 The linear bottom friction parameterisation (including the special case of a free-slip condition) assumes that 
    968 the bottom friction is proportional to the interior velocity (\ie the velocity of the last model level): 
     959 \subsection[Linear top/bottom friction (\forcode{ln_lin = .true.})] 
     960 {Linear top/bottom friction (\protect\np{ln\_lin}\forcode{ = .true.)}} 
     961 \label{subsec:ZDF_drg_linear} 
     962 
     963The linear friction parameterisation (including the special case of a free-slip condition) assumes that 
     964the friction is proportional to the interior velocity (\ie the velocity of the first/last model level): 
    969965\[ 
    970966  % \label{eq:zdfbfr_linear} 
    971967  {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 
    972968\] 
    973 where $r$ is a friction coefficient expressed in ms$^{-1}$. 
     969where $r$ is a friction coefficient expressed in $m s^{-1}$. 
    974970This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean,  
    975971and setting $r = H / \tau$, where $H$ is the ocean depth. 
     
    981977and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 
    982978This is the default value used in \NEMO. It corresponds to a decay time scale of 115~days. 
    983 It can be changed by specifying \np{rn\_bfri1} (namelist parameter). 
    984  
    985 For the linear friction case the coefficients defined in the general expression \autoref{eq:zdfbfr_bdef} are:  
     979It can be changed by specifying \np{rn\_Uc0} (namelist parameter). 
     980 
     981 For the linear friction case the drag coefficient used in the general expression \autoref{eq:zdfbfr_bdef} is:  
    986982\[ 
    987983  % \label{eq:zdfbfr_linbfr_b} 
    988   \begin{split} 
    989     c_b^u &= - r\\ 
    990     c_b^v &= - r\\ 
    991   \end{split} 
    992 \] 
    993 When \np{nn\_botfr}\forcode{ = 1}, the value of $r$ used is \np{rn\_bfri1}. 
    994 Setting \np{nn\_botfr}\forcode{ = 0} is equivalent to setting $r=0$ and 
    995 leads to a free-slip bottom boundary condition. 
    996 These values are assigned in \mdl{zdfbfr}. 
    997 From v3.2 onwards there is support for local enhancement of these values via an externally defined 2D mask array 
    998 (\np{ln\_bfr2d}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file. 
     984    c_b^T = - r 
     985\] 
     986When \forcode{ln_lin = .true.}, the value of $r$ used is \np{rn\_Uc0}*\np{rn\_Cd0}. 
     987Setting \forcode{ln_OFF = .true.} (and \forcode{ln_lin = .true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition. 
     988These values are assigned in \mdl{zdfdrg}. 
     989Note that there is support for local enhancement of these values via an externally defined 2D mask array 
     990(\np{ln\_boost}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file. 
    999991The mask values should vary from 0 to 1. 
    1000992Locations with a non-zero mask value will have the friction coefficient increased by 
    1001 $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}. 
     993$mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}. 
    1002994 
    1003995% ------------------------------------------------------------------------------------------------------------- 
    1004996%       Non-Linear Bottom Friction 
    1005997% ------------------------------------------------------------------------------------------------------------- 
    1006 \subsection[Non-linear bottom friction (\forcode{nn_botfr = 2})] 
    1007 {Non-linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 2})} 
    1008 \label{subsec:ZDF_bfr_nonlinear} 
    1009  
    1010 The non-linear bottom friction parameterisation assumes that the bottom friction is quadratic:  
    1011 \[ 
    1012   % \label{eq:zdfbfr_nonlinear} 
     998 \subsection[Non-linear top/bottom friction (\forcode{ln_no_lin = .true.})] 
     999 {Non-linear top/bottom friction (\protect\np{ln\_no\_lin}\forcode{ = .true.})} 
     1000 \label{subsec:ZDF_drg_nonlinear} 
     1001 
     1002The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic:  
     1003\[ 
     1004  % \label{eq:zdfdrg_nonlinear} 
    10131005  {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 
    10141006  }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b 
    10151007\] 
    1016 where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy due to tides, 
     1008where $C_D$ is a drag coefficient, and $e_b $ a top/bottom turbulent kinetic energy due to tides, 
    10171009internal waves breaking and other short time scale currents. 
    10181010A typical value of the drag coefficient is $C_D = 10^{-3} $. 
     
    10201012$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{killworth_JPO92} uses $C_D = 1.4\;10^{-3}$ and 
    10211013$e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. 
    1022 The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} namelist parameters). 
    1023  
    1024 As for the linear case, the bottom friction is imposed in the code by adding the trend due to 
    1025 the bottom friction to the general momentum trend in \mdl{dynbfr}. 
    1026 For the non-linear friction case the terms computed in \mdl{zdfbfr} are: 
    1027 \[ 
    1028   % \label{eq:zdfbfr_nonlinbfr} 
    1029   \begin{split} 
    1030     c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\ 
    1031     c_b^v &= - \; C_D\;\left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ 
    1032   \end{split} 
    1033 \] 
    1034  
    1035 The coefficients that control the strength of the non-linear bottom friction are initialised as namelist parameters: 
    1036 $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. 
    1037 Note for applications which treat tides explicitly a low or even zero value of \np{rn\_bfeb2} is recommended. 
    1038 From v3.2 onwards a local enhancement of $C_D$ is possible via an externally defined 2D mask array 
    1039 (\np{ln\_bfr2d}\forcode{ = .true.}). 
    1040 This works in the same way as for the linear bottom friction case with non-zero masked locations increased by 
    1041 $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}. 
     1014The CME choices have been set as default values (\np{rn\_Cd0} and \np{rn\_ke0} namelist parameters). 
     1015 
     1016As for the linear case, the friction is imposed in the code by adding the trend due to 
     1017the friction to the general momentum trend in \mdl{dynzdf}. 
     1018For the non-linear friction case the term computed in \mdl{zdfdrg} is: 
     1019\[ 
     1020  % \label{eq:zdfdrg_nonlinbfr} 
     1021    c_b^T = - \; C_D\;\left[ \left(\bar{u_b}^{i}\right)^2 + \left(\bar{v_b}^{j}\right)^2 + e_b \right]^{1/2} 
     1022\] 
     1023 
     1024The coefficients that control the strength of the non-linear friction are initialised as namelist parameters: 
     1025$C_D$= \np{rn\_Cd0}, and $e_b$ =\np{rn\_bfeb2}. 
     1026Note that for applications which consider tides explicitly, a low or even zero value of \np{rn\_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array 
     1027(\np{ln\_boost}\forcode{ = .true.}). 
     1028This works in the same way as for the linear friction case with non-zero masked locations increased by 
     1029$mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}. 
    10421030 
    10431031% ------------------------------------------------------------------------------------------------------------- 
    10441032%       Bottom Friction Log-layer 
    10451033% ------------------------------------------------------------------------------------------------------------- 
    1046 \subsection[Log-layer bottom friction enhancement (\forcode{nn_botfr = 2}, \forcode{ln_loglayer = .true.})] 
    1047 {Log-layer bottom friction enhancement (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})} 
    1048 \label{subsec:ZDF_bfr_loglayer} 
    1049  
    1050 In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally enhanced using 
    1051 a "law of the wall" scaling. 
    1052 If  \np{ln\_loglayer} = .true., $C_D$ is no longer constant but is related to the thickness of 
    1053 the last wet layer in each column by: 
    1054 \[ 
    1055   C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2 
    1056 \] 
    1057  
    1058 \noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness length provided via 
    1059 the namelist. 
    1060  
    1061 For stability, the drag coefficient is bounded such that it is kept greater or equal to 
    1062 the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional namelist parameter: 
    1063 \np{rn\_bfri2\_max}, \ie 
    1064 \[ 
    1065   rn\_bfri2 \leq C_D \leq rn\_bfri2\_max 
    1066 \] 
    1067  
    1068 \noindent Note also that a log-layer enhancement can also be applied to the top boundary friction if 
    1069 under ice-shelf cavities are in use (\np{ln\_isfcav}\forcode{ = .true.}). 
    1070 In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}. 
    1071  
    1072 % ------------------------------------------------------------------------------------------------------------- 
    1073 %       Bottom Friction stability 
    1074 % ------------------------------------------------------------------------------------------------------------- 
    1075 \subsection{Bottom friction stability considerations} 
    1076 \label{subsec:ZDF_bfr_stability} 
    1077  
    1078 Some care needs to exercised over the choice of parameters to ensure that the implementation of 
    1079 bottom friction does not induce numerical instability. 
    1080 For the purposes of stability analysis, an approximation to \autoref{eq:zdfbfr_flux2} is: 
     1034 \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer = .true.})] 
     1035 {Log-layer top/bottom friction (\protect\np{ln\_loglayer}\forcode{ = .true.})} 
     1036 \label{subsec:ZDF_drg_loglayer} 
     1037 
     1038In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using 
     1039a "law of the wall" scaling. This assumes that the model vertical resolution can capture the logarithmic layer which typically occur for layers thinner than 1 m or so. 
     1040If  \np{ln\_loglayer} = .true., $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness): 
     1041\[ 
     1042  C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2 
     1043\] 
     1044 
     1045\noindent where $\kappa$ is the von-Karman constant and \np{rn\_z0} is a roughness length provided via the namelist. 
     1046 
     1047The drag coefficient is bounded such that it is kept greater or equal to 
     1048the base \np{rn\_Cd0} value which occurs where layer thicknesses become large and presumably logarithmic layers are not resolved at all. For stability reason, it is also not allowed to exceed the value of an additional namelist parameter: 
     1049\np{rn\_Cdmax}, \ie 
     1050\[ 
     1051  rn\_Cd0 \leq C_D \leq rn\_Cdmax 
     1052\] 
     1053 
     1054\noindent The log-layer enhancement can also be applied to the top boundary friction if 
     1055under ice-shelf cavities are activated (\np{ln\_isfcav}\forcode{ = .true.}). 
     1056%In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}. 
     1057 
     1058% ------------------------------------------------------------------------------------------------------------- 
     1059%       Explicit bottom Friction 
     1060% ------------------------------------------------------------------------------------------------------------- 
     1061 \subsection{Explicit top/bottom friction (\forcode{ln_drgimp = .false.})} 
     1062 \label{subsec:ZDF_drg_stability} 
     1063 
     1064Setting \forcode{ln_drgimp = .false.} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads: 
     1065 
     1066At the top (below an ice shelf cavity): 
     1067\[ 
     1068  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 
     1069  = c_{t}^{\textbf{U}}\textbf{u}^{n-1}_{t} 
     1070\] 
     1071 
     1072At the bottom (above the sea floor): 
     1073\[ 
     1074  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 
     1075  = c_{b}^{\textbf{U}}\textbf{u}^{n-1}_{b} 
     1076\] 
     1077 
     1078Since this is conditionally stable, some care needs to exercised over the choice of parameters to ensure that the implementation of explicit top/bottom friction does not induce numerical instability. 
     1079For the purposes of stability analysis, an approximation to \autoref{eq:zdfdrg_flux2} is: 
    10811080\begin{equation} 
    1082   \label{eq:Eqn_bfrstab} 
     1081  \label{eq:Eqn_drgstab} 
    10831082  \begin{split} 
    10841083    \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\ 
     
    10861085  \end{split} 
    10871086\end{equation} 
    1088 \noindent where linear bottom friction and a leapfrog timestep have been assumed. 
    1089 To ensure that the bottom friction cannot reverse the direction of flow it is necessary to have: 
     1087\noindent where linear friction and a leapfrog timestep have been assumed. 
     1088To ensure that the friction cannot reverse the direction of flow it is necessary to have: 
    10901089\[ 
    10911090  |\Delta u| < \;|u|  
    10921091\] 
    1093 \noindent which, using \autoref{eq:Eqn_bfrstab}, gives: 
     1092\noindent which, using \autoref{eq:Eqn_drgstab}, gives: 
    10941093\[ 
    10951094  r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ 
     
    11051104For most applications, with physically sensible parameters these restrictions should not be of concern. 
    11061105But caution may be necessary if attempts are made to locally enhance the bottom friction parameters.  
    1107 To ensure stability limits are imposed on the bottom friction coefficients both 
     1106To ensure stability limits are imposed on the top/bottom friction coefficients both 
    11081107during initialisation and at each time step. 
    1109 Checks at initialisation are made in \mdl{zdfbfr} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case). 
     1108Checks at initialisation are made in \mdl{zdfdrg} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case). 
    11101109The number of breaches of the stability criterion are reported as well as 
    11111110the minimum and maximum values that have been set. 
    1112 The criterion is also checked at each time step, using the actual velocity, in \mdl{dynbfr}. 
    1113 Values of the bottom friction coefficient are reduced as necessary to ensure stability; 
     1111The criterion is also checked at each time step, using the actual velocity, in \mdl{dynzdf}. 
     1112Values of the friction coefficient are reduced as necessary to ensure stability; 
    11141113these changes are not reported. 
    11151114 
    1116 Limits on the bottom friction coefficient are not imposed if the user has elected to 
    1117 handle the bottom friction implicitly (see \autoref{subsec:ZDF_bfr_imp}). 
     1115Limits on the top/bottom friction coefficient are not imposed if the user has elected to 
     1116handle the friction implicitly (see \autoref{subsec:ZDF_drg_imp}). 
    11181117The number of potential breaches of the explicit stability criterion are still reported for information purposes. 
    11191118 
     
    11211120%       Implicit Bottom Friction 
    11221121% ------------------------------------------------------------------------------------------------------------- 
    1123 \subsection[Implicit bottom friction (\forcode{ln_bfrimp = .true.})] 
    1124 {Implicit bottom friction (\protect\np{ln\_bfrimp}\forcode{ = .true.})} 
    1125 \label{subsec:ZDF_bfr_imp} 
     1122 \subsection[Implicit top/bottom friction (\forcode{ln_drgimp = .true.})] 
     1123 {Implicit top/bottom friction (\protect\np{ln\_drgimp}\forcode{ = .true.})} 
     1124 \label{subsec:ZDF_drg_imp} 
    11261125 
    11271126An optional implicit form of bottom friction has been implemented to improve model stability. 
    1128 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\_bfrimp} to \forcode{.true.} in the \textit{nambfr} namelist. 
    1130 This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \textit{namzdf} namelist.  
    1131  
    1132 This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, 
    1133 the bottom boundary condition is implemented implicitly. 
    1134  
    1135 \[ 
    1136   % \label{eq:dynzdf_bfr} 
    1137   \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk} 
    1138   = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}} 
    1139 \] 
    1140  
    1141 where $mbk$ is the layer number of the bottom wet layer. 
    1142 Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so, it is implicit. 
    1143  
    1144 If split-explicit time splitting is used, care must be taken to avoid the double counting of the bottom friction in 
    1145 the 2-D barotropic momentum equations. 
    1146 As NEMO only updates the barotropic pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, 
    1147 we need to remove the bottom friction induced by these two terms which has been included in the 3-D momentum trend  
    1148 and update it with the latest value. 
    1149 On the other hand, the bottom friction contributed by the other terms 
    1150 (\eg the advection term, viscosity term) has been included in the 3-D momentum equations and 
    1151 should not be added in the 2-D barotropic mode. 
    1152  
    1153 The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the following: 
    1154  
    1155 \[ 
    1156   % \label{eq:dynspg_ts_bfr1} 
    1157   \frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b} 
    1158   \left(\textbf{U}_{med}-\textbf{U}^{m-1}\right) 
    1159 \] 
    1160 \[ 
    1161   \frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+ 
    1162   \left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)- 
    1163   2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right) 
    1164 \] 
    1165  
    1166 where $\textbf{T}$ is the vertical integrated 3-D momentum trend. 
    1167 We assume the leap-frog time-stepping is used here. 
    1168 $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step. 
    1169 $c_{b}$ is the friction coefficient. 
    1170 $\eta$ is the sea surface level calculated in the barotropic loops while $\eta^{'}$ is the sea surface level used in 
    1171 the 3-D baroclinic mode. 
    1172 $\textbf{u}_{b}$ is the bottom layer horizontal velocity. 
    1173  
    1174 % ------------------------------------------------------------------------------------------------------------- 
    1175 %       Bottom Friction with split-explicit time splitting 
    1176 % ------------------------------------------------------------------------------------------------------------- 
    1177 \subsection[Bottom friction with split-explicit time splitting (\texttt{ln\_bfrimp})] 
    1178 {Bottom friction with split-explicit time splitting (\protect\np{ln\_bfrimp})} 
    1179 \label{subsec:ZDF_bfr_ts} 
    1180  
    1181 When calculating the momentum trend due to bottom friction in \mdl{dynbfr}, 
    1182 the bottom velocity at the before time step is used. 
    1183 This velocity includes both the baroclinic and barotropic components which is appropriate when 
    1184 using either the explicit or filtered surface pressure gradient algorithms 
    1185 (\key{dynspg\_exp} or \key{dynspg\_flt}). 
    1186 Extra attention is required, however, when using split-explicit time stepping (\key{dynspg\_ts}). 
    1187 In this case the free surface equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, 
    1188 while the three dimensional prognostic variables are solved with the longer time step of \np{rn\_rdt} seconds. 
    1189 The trend in the barotropic momentum due to bottom friction appropriate to this method is that given by 
    1190 the selected parameterisation (\ie linear or non-linear bottom friction) computed with 
    1191 the evolving velocities at each barotropic timestep.  
    1192  
    1193 In the case of non-linear bottom friction, we have elected to partially linearise the problem by 
    1194 keeping the coefficients fixed throughout the barotropic time-stepping to those computed in 
    1195 \mdl{zdfbfr} using the now timestep. 
    1196 This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to: 
    1197  
     1127We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting. 
     1128This option can be invoked by setting \np{ln\_drgimp} to \forcode{.true.} in the \textit{namdrg} namelist. 
     1129%This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \textit{namzdf} namelist.  
     1130 
     1131This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step:  
     1132 
     1133At the top (below an ice shelf cavity): 
     1134\[ 
     1135  % \label{eq:dynzdf_drg_top} 
     1136  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 
     1137  = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t} 
     1138\] 
     1139 
     1140At the bottom (above the sea floor): 
     1141\[ 
     1142  % \label{eq:dynzdf_drg_bot} 
     1143  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 
     1144  = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b} 
     1145\] 
     1146 
     1147where $t$ and $b$ refers to top and bottom layers respectively.  
     1148Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit. 
     1149 
     1150% ------------------------------------------------------------------------------------------------------------- 
     1151%       Bottom Friction with split-explicit free surface 
     1152% ------------------------------------------------------------------------------------------------------------- 
     1153 \subsection[Bottom friction with split-explicit free surface] 
     1154 {Bottom friction with split-explicit free surface} 
     1155 \label{subsec:ZDF_drg_ts} 
     1156 
     1157With 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 \forcode{ln_drgimp = .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 \forcode{ln_drgimp = .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 
     1160The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO is as follows: 
    11981161\begin{enumerate} 
    1199 \item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before barotropic velocity to 
    1200   the bottom friction component of the vertically integrated momentum trend. 
    1201   Note the same stability check that is carried out on the bottom friction coefficient in \mdl{dynbfr} has to 
    1202   be applied here to ensure that the trend removed matches that which was added in \mdl{dynbfr}. 
    1203 \item At each barotropic step, compute the contribution of the current barotropic velocity to 
    1204   the trend due to bottom friction. 
    1205   Add this contribution to the vertically integrated momentum trend. 
    1206   This contribution is handled implicitly which eliminates the need to impose a stability criteria on 
    1207   the values of the bottom friction coefficient within the barotropic loop.  
    1208 \end{enumerate} 
    1209  
    1210 Note that the use of an implicit formulation within the barotropic loop for the bottom friction trend means that 
    1211 any limiting of the bottom friction coefficient in \mdl{dynbfr} does not adversely affect the solution when 
    1212 using split-explicit time splitting. 
    1213 This is because the major contribution to bottom friction is likely to come from the barotropic component which 
    1214 uses the unrestricted value of the coefficient. 
    1215 However, if the limiting is thought to be having a major effect 
    1216 (a more likely prospect in coastal and shelf seas applications) then 
    1217 the fully implicit form of the bottom friction should be used (see \autoref{subsec:ZDF_bfr_imp}) 
    1218 which can be selected by setting \np{ln\_bfrimp} $=$ \forcode{.true.}. 
    1219  
    1220 Otherwise, the implicit formulation takes the form: 
    1221 \[ 
    1222   % \label{eq:zdfbfr_implicitts} 
    1223   \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ] 
    1224 \] 
    1225 where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height),  
    1226 $c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and 
    1227 $RHS$ represents all the components to the vertically integrated momentum trend except for 
    1228 that due to bottom friction. 
    1229  
    1230 % ================================================================ 
    1231 % Tidal Mixing 
    1232 % ================================================================ 
    1233 \section[Tidal mixing (\texttt{\textbf{key\_zdftmx}})] 
    1234 {Tidal mixing (\protect\key{zdftmx})} 
    1235 \label{sec:ZDF_tmx} 
    1236  
    1237 %--------------------------------------------namzdf_tmx-------------------------------------------------- 
     1162\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\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 
     1166Note that other strategies are possible, like considering vertical diffusion step in advance, \ie prior barotropic integration.   
     1167 
     1168 
     1169% ================================================================ 
     1170% Internal wave-driven mixing 
     1171% ================================================================ 
     1172\section[Internal wave-driven mixing (\forcode{ln_zdfiwm = .true.})] 
     1173{Internal wave-driven mixing (\protect\np{ln\_zdfiwm}\forcode{ = .true.})} 
     1174\label{subsec:ZDF_tmx_new} 
     1175 
     1176%--------------------------------------------namzdf_iwm------------------------------------------ 
    12381177% 
    1239 %\nlst{namzdf_tmx} 
    1240 %-------------------------------------------------------------------------------------------------------------- 
    1241  
    1242  
    1243 % ------------------------------------------------------------------------------------------------------------- 
    1244 %        Bottom intensified tidal mixing  
    1245 % ------------------------------------------------------------------------------------------------------------- 
    1246 \subsection{Bottom intensified tidal mixing} 
    1247 \label{subsec:ZDF_tmx_bottom} 
    1248  
    1249 Options are defined through the  \ngn{namzdf\_tmx} namelist variables. 
    1250 The parameterization of tidal mixing follows the general formulation for the vertical eddy diffusivity proposed by 
    1251 \citet{st-laurent.simmons.ea_GRL02} and first introduced in an OGCM by \citep{simmons.jayne.ea_OM04}.  
    1252 In this formulation an additional vertical diffusivity resulting from internal tide breaking, 
    1253 $A^{vT}_{tides}$ is expressed as a function of $E(x,y)$, 
    1254 the energy transfer from barotropic tides to baroclinic tides: 
    1255 \begin{equation} 
    1256   \label{eq:Ktides} 
    1257   A^{vT}_{tides} =  q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 } 
    1258 \end{equation} 
    1259 where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 
    1260 $\rho$ the density, $q$ the tidal dissipation efficiency, and $F(z)$ the vertical structure function.  
    1261  
    1262 The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter) and 
    1263 is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980).  
    1264 The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter) 
    1265 represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally, 
    1266 with the remaining $1-q$ radiating away as low mode internal waves and 
    1267 contributing to the background internal wave field. 
    1268 A value of $q=1/3$ is typically used \citet{st-laurent.simmons.ea_GRL02}. 
    1269 The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical. 
    1270 It is implemented as a simple exponential decaying upward away from the bottom, 
    1271 with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter, 
    1272 with a typical value of $500\,m$) \citep{st-laurent.nash_DSR04},  
    1273 \[ 
    1274   % \label{eq:Fz} 
    1275   F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) } 
    1276 \] 
    1277 and is normalized so that vertical integral over the water column is unity.  
    1278  
    1279 The associated vertical viscosity is calculated from the vertical diffusivity assuming a Prandtl number of 1, 
    1280 \ie $A^{vm}_{tides}=A^{vT}_{tides}$. 
    1281 In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity is capped at $300\,cm^2/s$ and 
    1282 impose a lower limit on $N^2$ of \np{rn\_n2min} usually set to $10^{-8} s^{-2}$. 
    1283 These bounds are usually rarely encountered. 
    1284  
    1285 The internal wave energy map, $E(x, y)$ in \autoref{eq:Ktides}, is derived from a barotropic model of 
    1286 the tides utilizing a parameterization of the conversion of barotropic tidal energy into internal waves. 
    1287 The essential goal of the parameterization is to represent the momentum exchange between the barotropic tides and 
    1288 the unrepresented internal waves induced by the tidal flow over rough topography in a stratified ocean. 
    1289 In the current version of \NEMO, the map is built from the output of 
    1290 the barotropic global ocean tide model MOG2D-G \citep{carrere.lyard_GRL03}. 
    1291 This model provides the dissipation associated with internal wave energy for the M2 and K1 tides component 
    1292 (\autoref{fig:ZDF_M2_K1_tmx}). 
    1293 The S2 dissipation is simply approximated as being $1/4$ of the M2 one. 
    1294 The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$. 
    1295 Its global mean value is $1.1$ TW, 
    1296 in agreement with independent estimates \citep{egbert.ray_N00, egbert.ray_JGR01}.  
    1297  
    1298 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1299 \begin{figure}[!t] 
    1300   \begin{center} 
    1301     \includegraphics[width=\textwidth]{Fig_ZDF_M2_K1_tmx} 
    1302     \caption{ 
    1303       \protect\label{fig:ZDF_M2_K1_tmx} 
    1304       (a) M2 and (b) K1 internal wave drag energy from \citet{carrere.lyard_GRL03} ($W/m^2$). 
    1305     } 
    1306   \end{center} 
    1307 \end{figure} 
    1308 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
    1309   
    1310 % ------------------------------------------------------------------------------------------------------------- 
    1311 %        Indonesian area specific treatment  
    1312 % ------------------------------------------------------------------------------------------------------------- 
    1313 \subsection[Indonesian area specific treatment (\texttt{ln\_zdftmx\_itf})] 
    1314 {Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})} 
    1315 \label{subsec:ZDF_tmx_itf} 
    1316  
    1317 When the Indonesian Through Flow (ITF) area is included in the model domain, 
    1318 a specific treatment of tidal induced mixing in this area can be used. 
    1319 It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide an input NetCDF file, 
    1320 \ifile{mask\_itf}, which contains a mask array defining the ITF area where the specific treatment is applied. 
    1321  
    1322 When \np{ln\_tmx\_itf}\forcode{ = .true.}, the two key parameters $q$ and $F(z)$ are adjusted following 
    1323 the parameterisation developed by \citet{koch-larrouy.madec.ea_GRL07}: 
    1324  
    1325 First, the Indonesian archipelago is a complex geographic region with a series of 
    1326 large, deep, semi-enclosed basins connected via numerous narrow straits. 
    1327 Once generated, internal tides remain confined within this semi-enclosed area and hardly radiate away. 
    1328 Therefore all the internal tides energy is consumed within this area. 
    1329 So it is assumed that $q = 1$, \ie all the energy generated is available for mixing. 
    1330 Note that for test purposed, the ITF tidal dissipation efficiency is a namelist parameter (\np{rn\_tfe\_itf}). 
    1331 A value of $1$ or close to is this recommended for this parameter. 
    1332  
    1333 Second, the vertical structure function, $F(z)$, is no more associated with a bottom intensification of the mixing, 
    1334 but with a maximum of energy available within the thermocline. 
    1335 \citet{koch-larrouy.madec.ea_GRL07} have suggested that the vertical distribution of 
    1336 the energy dissipation proportional to $N^2$ below the core of the thermocline and to $N$ above.  
    1337 The resulting $F(z)$ is: 
    1338 \[ 
    1339   % \label{eq:Fz_itf} 
    1340   F(i,j,k) \sim     \left\{ 
    1341     \begin{aligned} 
    1342       \frac{q\,\Gamma E(i,j) } {\rho N \, \int N     dz}    \qquad \text{when $\partial_z N < 0$} \\ 
    1343       \frac{q\,\Gamma E(i,j) } {\rho     \, \int N^2 dz}    \qquad \text{when $\partial_z N > 0$} 
    1344     \end{aligned} 
    1345   \right. 
    1346 \] 
    1347  
    1348 Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$,  
    1349 which agrees with the independent estimates inferred from observations. 
    1350 Introduced in a regional OGCM, the parameterization improves the water mass characteristics in 
    1351 the different Indonesian seas, suggesting that the horizontal and vertical distributions of 
    1352 the mixing are adequately prescribed \citep{koch-larrouy.madec.ea_GRL07, koch-larrouy.madec.ea_OD08*a, koch-larrouy.madec.ea_OD08*b}. 
    1353 Note also that such a parameterisation has a significant impact on the behaviour of 
    1354 global coupled GCMs \citep{koch-larrouy.lengaigne.ea_CD10}. 
    1355  
    1356 % ================================================================ 
    1357 % Internal wave-driven mixing 
    1358 % ================================================================ 
    1359 \section[Internal wave-driven mixing (\texttt{\textbf{key\_zdftmx\_new}})] 
    1360 {Internal wave-driven mixing (\protect\key{zdftmx\_new})} 
    1361 \label{sec:ZDF_tmx_new} 
    1362  
    1363 %--------------------------------------------namzdf_tmx_new------------------------------------------ 
    1364 % 
    1365 %\nlst{namzdf_tmx_new} 
     1178\nlst{namzdf_iwm} 
    13661179%-------------------------------------------------------------------------------------------------------------- 
    13671180 
     
    13761189where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of 
    13771190the energy available for mixing. 
    1378 If the \np{ln\_mevar} namelist parameter is set to false, the mixing efficiency is taken as constant and 
     1191If the \np{ln\_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and 
    13791192equal to 1/6 \citep{osborn_JPO80}. 
    13801193In the opposite (recommended) case, $R_f$ is instead a function of 
     
    13861199 
    13871200In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary  
    1388 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice.  
     1201as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice.  
    13891202This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14}, 
    13901203is implemented as in \cite{de-lavergne.madec.ea_JPO16}. 
     
    13921205The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, 
    13931206is constructed from three static maps of column-integrated internal wave energy dissipation, 
    1394 $E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures 
    1395 (de Lavergne et al., in prep): 
     1207$E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures: 
     1208 
    13961209\begin{align*} 
    13971210  F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\ 
    1398   F_{pyc}(i,j,k) &\propto N^{n\_p}\\ 
     1211  F_{pyc}(i,j,k) &\propto N^{n_p}\\ 
    13991212  F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } 
    14001213\end{align*}  
     
    14041217  h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; , 
    14051218\] 
    1406 The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist) 
     1219The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_iwm} namelist) 
    14071220controls the stratification-dependence of the pycnocline-intensified dissipation. 
    1408 It can take values of 1 (recommended) or 2. 
     1221It can take values of $1$ (recommended) or $2$. 
    14091222Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of 
    14101223the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps. 
     
    14121225$h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of 
    14131226the abyssal hill topography \citep{goff_JGR10} and the latitude. 
     1227% 
     1228% Jc: input files names ? 
     1229 
     1230% ================================================================ 
     1231% surface wave-induced mixing  
     1232% ================================================================ 
     1233\section[Surface wave-induced mixing (\forcode{ln_zdfswm = .true.})] 
     1234{Surface wave-induced mixing (\protect\np{ln\_zdfswm}\forcode{ = .true.})} 
     1235\label{subsec:ZDF_swm} 
     1236 
     1237TBC ... 
     1238 
     1239% ================================================================ 
     1240% Adaptive-implicit vertical advection 
     1241% ================================================================ 
     1242\section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp = .true.})] 
     1243{Adaptive-implicit vertical advection(\protect\np{ln\_zad\_Aimp}\forcode{ = .true.})} 
     1244\label{subsec:ZDF_aimp} 
     1245 
     1246This refers to \citep{shchepetkin_OM15}. 
     1247 
     1248TBC ... 
     1249 
     1250 
    14141251 
    14151252% ================================================================ 
Note: See TracChangeset for help on using the changeset viewer.