This page explains in detail the formulations of two turbulence schemes in (ICOsahedral Nonhydrostatic) ICON models: Total turbulent energy (TTE) turbulence scheme in ICON-AES and 3D Smagorinsky turbulence scheme in ICON-LEM. Their locations in the model codes are also provided. If you have feedback or questions, please contact Junhong Lee (junhong.lee@mpimet.mpg.de), Andreas Chlond (andreas.chlond@mpimet.mpg.de), or Cathy Hohenegger (cathy.hohenegger@mpimet.mpg.de).

This page (document) is also available as pdf: exchange_coefficient_v1.8_junhong.pdf

Total turbulent energy (TTE) turbulence scheme was implemented into ICON-AES by Felix Pithan and Thorsten Mauritsen in 2016. This chapter will provide description of TTE scheme based on the codes and Louis (1979), Brinkop and Roeckner (1995), Mauritsen et al. (2007), Zilitinkevich et al. (2007), Angevine et al. (2010), Wan (2011), and Pithan (2014). These references will be referred as L79, BR95, M07, Z07, A10, W11, and P14 for the convenience.

Detail descriptions of exchange coefficient in atmosphere (from $k=1$ to $k= \textrm{nlev}$) are provided in this section, where nlev is the lowest model level.
All the equations in section 2.1 are in  subroutine  atm_exchange_coeff  of  mo_turbulence_diag.f90

The TTE scheme prognoses total turbulent energy ($E=E_k+E_p$) rather than turbulent kinetic energy unlike other turbulent prognostic schemes, where $E_k$ is the turbulent kinetic energy and $E_p$ is the turbulent potential energy. The prognostic equation is given by $E$ is given by
$$\frac{D E}{D t} = \tau \cdot S - \gamma - \cfrac{\partial F_E}{\partial z} +\left\{ \def\arraystretch{2.0} \begin{array}{ll} 0 & \textrm{for} \ N^2 \geq 0\\ 2 \cfrac{g}{\theta} \overline{w'\theta'} & \textrm{for} \ N^2 < 0 \end{array} \right., \tag{2.1} $$

from M07 Eq. (4), where $\tau$ is the turbulent stress, $S$ is the wind shear and $N$ is the Brunt-Väisälä frequency, $\gamma=\cfrac{C_{\epsilon}}{l}E^{3/2}$ is the dissipation rate, $C_{\epsilon}$ is empirical constant, $F_E=-|S| l^2 \cfrac{\partial E}{\partial z}$ is the turbulent flux of total turbulent energy, and $l$ is length scale (mixing length). $E$ without vertical flux is computed using operator splitting method as

$$\cfrac{\partial E}{\partial t} = \left( \cfrac{\partial E}{\partial t} \right)_{local} + \left( \cfrac{\partial E}{\partial t} \right)_{diffusion}, \tag{2.2}$$

where the first term of r.h.s is local process and computed first, and the second term of r.h.s is vertical diffusion process and will be computed using VDIFF module like momentum and dry static. Local process term is given by

$$ \left( \cfrac{\partial E}{\partial t} \right)_{local} = \mathcal{B} \sqrt{E} - \mathcal{C} (\sqrt{E})^3 \label{eq:2.3} \tag{2.3} $$
where $\mathcal{B}$ and $\mathcal{C}$ are
$i$) in stable case ($Ri \geq 0$): $$\mathcal{B} = K_m \ S^2 \tag{2.4}$$
$ii$) in unstable case ($Ri < 0$): $$\mathcal{B} = K_m \ S^2 - 2 K_h \ N^2 \tag{2.5}$$
and

$$\mathcal{C} = \cfrac{C_{\epsilon}}{l}. \tag{2.6}$$

Convert Eq. (\ref{eq:2.3}) into an equation of $\sqrt{E}$

$$\cfrac{D \sqrt{E}}{D t} = \cfrac{\mathcal{B}}{2} - \cfrac{\mathcal{C}}{2} (\sqrt{E})^2 , \tag{2.7}$$
and discretize using implicit time stepping

$$\cfrac{\sqrt{E^{(*)}}-\sqrt{E^{(t)}}} {\Delta t} = \cfrac{\mathcal{B}}{2} - \cfrac{\mathcal{C}}{2} (\sqrt{E^{(*)}})^2 . \tag{2.8}$$
Finally, $E^{(*)}$ can be obtained by solving a quadrature equation

$$\sqrt{ E^{(*)} } = \cfrac{ -1 + \sqrt{1 + \mathcal{C} \ \Delta t \ \left[ \mathcal{B} \Delta t + 2 \sqrt{E^{(t)}} \right] } } {\mathcal{C} \ \Delta t} \tag{2.9}$$
from BR95 Eqs. (A3-A10) and W11 section 5.3, where $E^{(*)}$ is $E$ at $t+\Delta t$, after local process, $K_m$ and $K_h$ exchange coefficients for momentum and heat. Then, $E$ is updated once more to get $E^{(t+\Delta t)}$ by solving vertical diffusion process in VDIFF module like momentum and dry static energy in the TTE scheme.
Once total turbulent energy is updated, $E_k$ and $E_p$ are derived from $E^{(*)}$. The $E_k$ and $E_p$ can be obtained from that $E = E_k + E_p$ as follows
$i$) in stable case ($Ri \geq 0$): $$E_k = E^{(t)} \left[1 + \cfrac{Ri}{3 Ri + Pr_0} \right]^{-1}, \tag{2.10}$$

and

$$E_p = E_k \left[\cfrac{Ri}{3 Ri + Pr_0} \right], \tag{2.11}$$

$ii$) in unstable case ($Ri < 0$): $$E_k = E^{(t)} \left[1 + \cfrac{Ri}{2 Ri - Pr_0} \right]^{-1}, \tag{2.12}$$

and

$$E_p = E_k \left[\cfrac{Ri}{2 Ri - Pr_0} \right], \tag{2.13}$$

based on the relationship

$$\frac{E_p}{E_k} =\left\{ \def\arraystretch{2.5} \begin{array}{lrl} \cfrac{Ri}{3 \cdot Ri + Pr_0} & \textrm{for} \ Ri \geq 0\\ \cfrac{Ri}{2 \cdot Ri - Pr_0} & \textrm{for} \ Ri < 0, \end{array} \right. \tag{2.14}$$

from P14 Eq. (5.9) and A10 Eq. (A10), where $Ri$ is moist Richardson number, $Pr_0$ is turbulent Prandtl number.

The $K_m$, exchange coefficients for momentum, and $K_h$, exchange coefficients for heat are derived from different formulas. Which formula will be used is determined by the height where $K$ is computed and convective boundary layer height ($h_d$). The $K_m$ and $K_h$ are calculated as
$i$) $z > h_d$: $$K_m = \cfrac{f^2_{\tau} E^2_k} {C_{\epsilon} E_k \sqrt{E}/l - \beta f_{\theta} \sqrt{E_k \sigma_{\theta}^2} }, \tag{2.15} $$

and

$$K_h = \cfrac{2 f_{\theta}^2 E_k l}{C_{\phi} \sqrt{E}}, \tag{2.16}$$

from M07 Eqs. (B1) and (B2),
$ii$) $z \leq h_d/2$): $$K_m = \frac{f^2_{\tau0}}{C_{\epsilon}} l_c \sqrt{E_k}, \tag{2.17}$$

and

$$K_h = Pr^{-1}_0 K_m, \tag{2.18}$$

from A10 Eqs. (A21) and (A22),
$iii$) $h_d/2 < z \leq h_d$):

$$K_m = \textrm{Max} \left[ \cfrac{f^2_{\tau} E^2_k} {C_{\epsilon} E_k \sqrt{E}/l - \beta f_{\theta} \sqrt{E_k \sigma_{\theta}^2} }\ , \ \frac{f^2_{\tau0}}{C_{\epsilon}} l_c \sqrt{E_k} \right], \tag{2.19}$$

and

$$K_h = \textrm{Max} \left[ \cfrac{2 f_{\theta}^2 E_k l}{C_{\phi} \sqrt{E}}\ , \ Pr^{-1}_0 K_m \right], \tag{2.20}$$
if in unstable case ($Ri < 0$), there is one more procedure as

$$K_m = K_m \left[ 1 - \cfrac{2c Ri} { 1 + 3 c^2 l^2 \left[ \left( \frac{\Delta z}{z} +1 \right)^{1/3} -1 \right]^{3/2} \left[ \frac{\sqrt{-Ri}}{(\Delta z)^{3/2} \sqrt{z}} \right] } \right], \tag{2.21}$$

and

$$K_h = K_h \left[ 1 - \cfrac{3c Ri} { 1 + 3 c^2 l^2 \left[ \left( \frac{\Delta z}{z} +1 \right)^{1/3} -1 \right]^{3/2} \left[ \frac{\sqrt{-Ri}}{(\Delta z)^{3/2} \sqrt{z}} \right] } \right], \tag{2.22}$$

from ECHAM6 report Eqs. (2.165) and (2.166), where $h_d$ is defined as the first model level whose dry static energy exceeds that of the lowest model level, $f_{\tau}$ ($f_{\theta}$) is nondimensional stress (heat flux), $C_{\phi}$ ($=C_{\epsilon}$) is empirical constant related to the turbulence dissipation, $l_c$ is convective length scale in convective condition, $\beta$ is $g/\theta_v$, $\sigma_{\theta}^2=2 E_p \frac{|N^2|}{\beta^2}$ is potential temperature variation, $Pr_0= K_m / K_h = f_{\tau}^2 / 2f_{\theta}^2$ is turbulent Prandtl number, and $c = 5.0$.

$Ri$, moist Richardson number can be defined as

$$N^2 = \cfrac{g}{\theta_v} \cfrac{A \Delta \theta_L + \theta D \Delta q_t}{\Delta z} \tag{2.23}$$

$$S^2 = \cfrac{(\Delta U)^2+(\Delta V)^2}{(\Delta z)^2} \tag{2.24}$$

$$ Ri = \cfrac{N^2}{S^2} \tag{2.25} \label{eq:2.25}$$

from BR95 Eq. (11) and M07 Eq. (5), where $\theta_L$ is cloud water potential temperature, and $q_t$ is total water content, last two variables are conserved under phase transitions and are defined as

$$\theta_L=\theta -\frac{L_v}{c_p} \frac{\theta}{T} m, \tag{2.26}$$

$$q_t=q+m, \tag{2.27}$$

from BR95 Eqs. (6) and (7), where $T$ is temperature, $\theta$ is potential temperature, $q$ is specific humidity and $m$ is cloud water specific humidity (cloud liquid water + cloud ice specific humidity). The coefficients $A$ and $D$ are defined as

$$A = C_f A_s + (1-C_f) A_u, \ \ \ D = C_f D_s + (1-C_f) D_u, \tag{2.28}$$

where $C_f$ is cloud cover faction, $A_s$, $A_u$, $D_s$, and $D_u$ are defined as

$$A_s = 1+0.61 q_t - \left[ \cfrac{L}{c_p T} (1+0.61 q_t) - 1.622 \right] \cfrac{0.622 \cfrac{L}{R_d T} \ q_s}{1 + 0.622 \cfrac{L}{R_d c_p T^2} \ q_s}, \tag{2.29}$$

$$D_s = A_s \frac{L}{c_pT} - 1, \tag{2.30}$$

$$A_u = 1+0.61 q, \tag{2.31}$$

$$D_u = 0.61, \tag{2.32}$$

from BR95 Eqs. (8-10), where $q_s$ is saturated specific himidity.

Nondimensional stress ($f_{\tau}$) and heat flux ($f_{\theta}$) are defined as
$i$) in stable case ($Ri \geq 0$): $$\begin{split} f_{\tau} & =\frac{\left| \tau \right|}{E_k} = f_{\tau 0} \left( 0.25+\cfrac{0.75}{1+4Ri} \right), \end{split} \tag{2.33}$$

and

$$f_{\theta} =\cfrac{\overline{w \theta}}{\sqrt{E_k \sigma^2_{\theta}}} = - \cfrac{f_{\theta 0}}{1+4Ri}, \tag{2.34}$$ from M07 Eq. (6) and Mauritsen and Svensson (2007),
$ii$) in unstable case ($Ri < 0$): $$f_{\tau} = f_{\tau 0}, \tag{2.35}$$

and

$$f_{\theta} = f_{\theta 0}, \tag{2.36}$$
where $f_{\tau 0}=0.17$, $f_{\theta 0}= -\sqrt{\cfrac{f_{\tau 0}^2} {2 Pr_0}} = - 0.1202081528$, and $Pr_0=1.0$. (In M07, $Pr_0=0.69$ and thus $f_{\theta 0}=-0.145$.)

The $l$, length scale (mixing length) is derived as

$$\cfrac{1}{l} = \cfrac{1}{k z} + \cfrac{f}{C_f \sqrt{f_{\tau} E_k}} + \cfrac{N}{C_N \sqrt{f_{\tau} E_k}}, \tag{2.37}$$

from M07 Eq. (8) and A10 Eq. (A11), where $k=0.4$ is von Kármán constant, $f$ is the Coriolis parameter, $C_f=0.185$ and $C_N=2.0$.

The $l_c$, convective length scale is used to define the length scale in convective boundary layer. The $l_c$ is computed as

$$\cfrac{1}{l_c} = \cfrac{1}{k z} + \cfrac{f}{C_f \sqrt{f_{\tau} E_k}} + \cfrac{3}{k (h_d - z)} \ \ \ \textrm{for} \ \ z < h_d, \tag{2.38}$$

from A10 Eq. (A20). The second term of r.h.s. is included in the code, but that is not included in A10 Eq. (A20). In personal conversation with Felix Pithan and Thorsten Mauritsen, they believe eddies in unstable condition should also feel the coriolis force (the second term of r.h.s) like in stable condition. Thus, I leave this as in the code.

Detail descriptions of surface transfer coefficient is provided in this section.
All the equations in section 2.2 are in  subroutine  sfc_exchange_coeff  of  mo_turbulence_diag.f90

To obtain $E$, first $E_k$ and $E_p$ are computed as Eqs. (2.10-13) and then $E$ is computed as
$i$) in stable case ($Ri \geq 0$): $$E^{(*)} = \left( 1 + \cfrac{E_p}{E_k} \right) \cfrac{u_*^2}{f_{\tau}} \tag{2.39}$$

$ii$) in unstable case ($Ri < 0$): $$E^{(*)} = \left( 1 + \cfrac{E_p}{E_k} \right) \cfrac{(u_*^3 + 2 l \frac{g}{\theta_v} \overline{w'\theta'} )^{2/3}}{f_{\tau}} \tag{2.40}$$
from P14 Eq. (5.17).
Again, $E$ is updated once more to get $E^{(t+\Delta t)}$ by solving vertical diffusion process in VDIFF module like momentum and dry static energy in the TTE scheme.

Surface transfer coefficient for momentum ($C_d$) and heat ($C_h$) are defined as

$$C_d = C_{N,d} \cdot f_d, \tag{2.41}$$

and

$$C_h = C_{N,h} \cdot f_h, \tag{2.42}$$

from ECHAM6 report Eq. (2.173) where $C_{N,d}$ and $C_{N,h}$ are surface transfer coefficients in neutral case for momentum and heat, and $f_d$ and $f_h$ are stability functions for momentum and heat.

Both of surface transfer coefficients are suggested as

$$C_{N,d} = \cfrac{l_{s}^2} { \left( f_{sl} \ z_1 \right)^2 \left[ \textrm{ln} \left( \frac{z_1}{z_0} \right) \right]^2 }, \tag{2.43}$$

and

$$C_{N,h} = \cfrac{l_{s}^2} { \left( f_{sl} \ z_1 \right)^2 \left[ \textrm{ln} \left( \frac{z_1}{z_0} \right) \right] \left[ \textrm{ln} \left( \frac{z_1}{z_{0h}} \right) \right] } \ \cfrac{1}{Pr_0}, \tag{2.44}$$

from M07 Eqs. (9) and (10) and P14 Eqs. (5.15) and (5.16) where $z_1$ is height of the first model full level, $l_s$ is length scale at the surface, $f_{sl}=0.4$ is the fraction of the first level height at which the surface fluxes are nominally evaluated, $z_0$ is roughness length for momentum, and $z_{0h}$ is roughness length for heat.

Both of stability functions can be obtained
$i$) in stable case ($Ri \geq 0$): $$f_d = \frac{f_{\tau}}{f_{\tau 0}} , \tag{2.45}$$

and

$$f_h = \frac{f_{\theta}}{f_{\theta 0}} \sqrt{ \frac{f_{\tau}}{f_{\tau 0}} } ,\tag{2.46}$$

from M07 Eqs. (9) and (10) and P14 Eqs. (5.15) and (5.16)
$ii$) in unstable case ($Ri < 0$): $$f_d = \left[ 1 - \cfrac{ 2c Ri} {1 + 3 c^2 C_{N,d} \ \sqrt{-Ri \left( \frac{z_1}{z_0} + 1 \right) }} \right] , \tag{2.47}$$

and

$$f_h = \left[ 1 - \cfrac{3 c Ri} {1 + 3 c^2 C_{N,h} \ \sqrt{-Ri \left( \frac{z_1}{z_{0h}} + 1 \right) }} \right] ,\tag{2.48}$$

from ECHAM6 report Eqs. (2.180) and (2.181) where $c=5$.

The $Ri$ at surface is determined same as Eq. (\ref{eq:2.25}) except $S$. In case of surface $S$ is suggested as

$$S^2 = \cfrac{(U_1)^2 +(V_1)^2 +(c_w w_*)^2} {z_1^2}, \tag{2.49}$$

from P14 Eq. (5.16) where $U_1$ and $V_1$ are wind speed at model full level and $c_w=0.5$ is the ratio of the mean absolute wind at the first level to the convective velocity scale under free convection.

The $l_s$, surface length scale is calculated as

$i$) in stable case ($Ri \geq 0$):
$$\cfrac{1}{l_s} = \cfrac{1}{k f_{sl} z_1} + \cfrac{f}{C_f \sqrt{f_{\tau} E_k}} + \cfrac{N}{C_N \sqrt{f_{\tau} E_k}} + \cfrac{3}{k (h_d - f_{sl} z_1)} , \tag{2.50}$$
$ii$) in unstable case ($Ri < 0$):

$$\cfrac{1}{l_s} = \cfrac{1}{k f_{sl} z_1} + \cfrac{f}{C_f \sqrt{f_{\tau} E_k}} + \cfrac{3}{k (h_d - f_{sl} z_1)} .\tag{2.51}$$

There is no reference for equation and it has a terms for unstable condition (the fourth term of r.h.s). In personal conversation with Felix Pithan and Thorsten Mauritsen, they think $h_d$ in stable condition are small enough so including this term does not make big differences from the case that this term is not included.

Detail description of momentum (and E) and temperature (and scalar) tendencies are provided in this section. The tendency in ICON-AES is calculated by VDIFF module which is implicit time stepping solver, since turbulent mixing is a very fast process compared to the typical time step used by global hydrostatic models (ECHAM6 report).

The tendency of prognostic variable ($\psi$), which can be momentum (and E) or heat (and scalar), due to turbulent motion is obtained from

$$\begin{split} \left( \frac{\partial \psi}{\partial t} \right)_{\textrm{turb}} & = \cfrac{\psi_{k}^{(t+\Delta t)} - \psi_{k}^{(t)} }{\Delta t} \\ & = -\frac{\partial \overline{w'\psi '}}{\partial p} = \frac{\partial}{\partial p} \left[ \rho g K_{\psi} \left( -\frac{\partial \hat{\psi}}{\partial z} \right) \right]. \end{split} \tag{2.52} \label{eq:2.52}$$
in  subroutine  vdiff_tendencies  of  mo_vdiff_solver.f90

To integrate Eq. (\ref{eq:2.52}), the model from time instance $t$ to $t+\Delta t$, the r.h.s of Eq. (\ref{eq:2.52}) has to be evaluated at intermediate time level $\hat{\psi}$ defined by

$$\label{eq:2.53} \hat{\psi}_k = \alpha \ \psi_k^{(t+\Delta t)} +(1-\alpha) \ \psi_k^{t}, \tag{2.53}$$

where $\alpha=1.5$ denotes the implicitness factor and follows its value the IFS model of the ECMWF.
After discretization of Eq. (\ref{eq:2.52}) and substitution of Eq. (\ref{eq:2.53}) in to Eq. (\ref{eq:2.52}), we can get $$- A_{k+1/2} \cfrac{\hat{\psi}_{k+1}}{\alpha} + B_{k+1/2} \cfrac{\hat{\psi_k}}{\alpha} - C_{k+1/2} \cfrac{\hat{\psi}_{k-1}}{\alpha} = \cfrac{\psi_k^{(t)}}{\alpha}, \tag{2.54}$$

where $A_{k+1/2}$, $B_{k+1/2}$, and $C_{k+1/2}$ are obtained as
$$A_{k+1/2} = \cfrac{1}{\Delta p_k} \Delta t \alpha g \cfrac{\rho_{k+1/2} \ K_{\psi, k+1/2}}{z_{k} - z_{k+1}}, \tag{2.55}$$

$$C_{k+1/2} = \cfrac{1}{\Delta p_k} \Delta t \alpha g \cfrac{\rho_{k-1/2} \ K_{\psi, k-1/2}}{z_{k-1} - z_{k}}, \tag{2.56}$$

$$B_{k+1/2} = 1 + A_{k+1/2} + C_{k+1/2}, \tag{2.57}$$

$i$) for $k=1$
$C_{k+1/2}$ should be rewritten as $$C_{k+1/2} = 0, \tag{2.58}$$
$ii$) for $k= \textrm{nlev}$
$A_{k+1/2}$ should be rewritten as $$A_{k+1/2} = C_{\psi} \Delta t \alpha \cfrac{\rho_k}{M_k} \sqrt{(U_1)^2 +(V_1)^2 +(c_w w_*)^2}, \tag{2.59}$$

$M_k$ is air mass [$kg \ m^{-2}$]

in  subroutine  matrix_setup_elim  of  mo_vdiff_solver.f90

Firstly, $E_{k+1/2}$ and $F_{k+1/2}$ are computed from $k=1$ to $k= \textrm{nlev}$ as

$$E_{k+1/2} = \cfrac{A_{k+1/2}} {B_{k+1/2} - C_{k+1/2} \ E_{k-1/2}} \tag{2.60}$$

and in  subroutine  matrix_setup_elim  of  mo_vdiff_solver.f90
for $k= \textrm{nlev}$, in  subroutine  matrix_to_richtmyer_coeff  of  mo_vdiff_solver.f90

$$F_{k+1/2} = \cfrac{ \psi_{k}^{(t)} / \alpha + C_{k+1/2} \ F_{k-1/2} } {B_{k+1/2} - C_{k+1/2} \ E_{k-1/2}}. \tag{2.61}$$

in  subroutine  rhs_elim  of  mo_vdiff_solver.f90
for $k= \textrm{nlev}$, in  subroutine  matrix_to_richtmyer_coeff  of  mo_vdiff_solver.f90

Then, $\hat{\psi}_k$ is computed from $k= \textrm{nlev}$ to $k=1$ as

$$\frac{\hat{\psi}_k}{\alpha} = E_{k+1/2} \ \frac{\hat{\psi}_{k+1}}{\alpha} + F_{k+1/2}. \tag{2.62}$$

in  subroutine  rhs_bksub  of  mo_vdiff_solver.f90

Finally, $\psi_k^{(t+\Delta t)}$ is calculated using Eq. (\ref{eq:2.53}) according to $\psi_k^{(t+\Delta t)} = \frac{1}{\alpha}\hat{\psi}_k - \frac{(1-\alpha)}{\alpha} \psi_k^{(t)}$ and the tendency is calculated as $\left( \frac{\partial \psi}{\partial t} \right)_{\textrm{turb}} = \frac{\psi_{k}^{(t+\Delta t)} - \psi_{k}^{(t)} }{\Delta t}$.

3D Smagrinsky turbulence scheme was implemented into ICON-LEM by Anurag Dipankar in 2013. This chapter will provide description of 3D Smagorinsky turbulence scheme based on the codes and Dipankar et al., 2015. I will refer Dipankar et al. (2015) as D15 for the convenience.

Detail descriptions of exchange coefficient in atmosphere are provided in this section. You can find further detail on grid system in Wan (2009) and Dipankar et al. (2015).

The $K_m$ at grid center and interface level is given by

$$K_m = \lambda^2 \rho \left| S \right| \left( 1-\cfrac{Ri}{Pr} \right)^{1/2} \ \ \textrm{for} \left( 1-\frac{Ri}{Pr} \right) > 0, \tag{3.1}$$

and

$$K_h = \cfrac{K_m}{Pr},\tag{3.2}$$

from D15 Eq. (10), where $\lambda$ is subgrid length scale, $|S|=|D|/\sqrt{2}$, $D$ is shear strain tensor, and $Pr=0.33333$. There is a comment on $K_m$ in the code: note that the factor $\sqrt{2}$ with $\lambda^2$ is considered into the Smagorinsky constant.
in  subroutine  smagorinsky_model  of  mo_sgs_turbulence.f90

Richardson number, $Ri$, is obtained as

$$Ri = \cfrac{N^2}{S^2}, \tag{3.3}$$

where $N^2$ and $S^2$ are

$$N^2 = \cfrac{g}{\theta_v} \cfrac{\Delta \theta_v}{\Delta z}, \tag{3.4}$$

in  subroutine  brunt_vaisala_freq  of  mo_les_utilities.f90

and

$$S^2 = \frac{D^2}{2}. \tag{3.5}$$
in  subroutine  smagorinsky_model  of  mo_sgs_turbulence.f90

The $\lambda$ subgrid length scale is given by

$$\cfrac{1}{\lambda^2} = \cfrac{1}{(C_s \Delta)^2} + \cfrac{1}{(k z)^2}, \tag{3.6}$$

where $C_s=0.23$ is Smagorinsky constant and $\Delta=(\Delta x \Delta y \Delta z)^{1/3}$.
in  subroutine  smagorinsky_model of mo_sgs_turbulence.f90

Shear strain tensor, $D$, will be suggested this subsection. It is better to define the coordinate system of the ICON before moving on the detail. $x_1$ is the horizontal axis normal to the triangle edge, $x_2$ is the horizontal axis to the triangle edge, $x_3$ is vertical axis (Figs. 1 and 2). Velocity components ($v_1$, $v_2$, $v_3$) increase along with these axises. $v_1$ is velocity component normal to the triangle edges at triangle edge and full level (can be referred as $v_n$), $v_2$ is tangential velocity component at triangle edge and full level (can be referred as $v_t$), and $v_3$ is vertical velocity component at triangle center and interface level. Furthermore, I would like to define notations used here. Position letters denote $e$: triangle edge, $v$: triangle vertex, $c$: triangle center, $f$: full level, and $i$: interface level (Fig. 1). The arrow on the superscript indicates interpolation from left letter to right letter. $\Delta$ operator with subscript indicate difference between the variables at right letter and left letter. $\Delta_{x_3}$ operator means difference between the variables at two interface levels. For example, $\Delta_{ab} \ v_1^{(e \rightarrow v)}$ indicates $v_1$ at $b$ minus $v_1$ at $a$ and two $v_1$ are interpolated into triangle vertex from triangle edge before subtraction.

Shear strain tensor is given by

$$\begin{split} D_{11} &=2 \left( \cfrac{\partial v_1}{\partial x_1} \right)\\ &=2 \left( \cfrac{\Delta_{ab} \ v_1^{(e \rightarrow v)}} {\Delta_{ab} \ x_1} \right), \end{split} \tag{3.7}$$

$$\begin{split} D_{12} &= \left( \cfrac{\partial v_1}{\partial x_2} + \cfrac{\partial v_2}{\partial x_1} \right) \\ &= \left( \cfrac{\Delta_{dc} \ v_1^{(e \rightarrow v)}} {\Delta_{dc} \ x_2} +\cfrac{\Delta_{ab} \ v_2^{(e \rightarrow v)}} {\Delta_{ab} \ x_1} \right), \end{split} \tag{3.8}$$

$$\begin{split} D_{13} &= \left( \cfrac{\partial v_1}{\partial x_3} + \cfrac{\partial v_3}{\partial x_1} \right) \\ &= \left( \cfrac{\Delta_{x_3} \ v_1^{(f \rightarrow i)}} {\Delta \ x_3} +\cfrac{\Delta_{qp} \ v_3^{(i \rightarrow f)}}{\Delta_{qp} \ x_1} \right), \end{split} \tag{3.9}$$

$$D_{21} = D_{12},\tag{3.10}$$

$$\begin{split} D_{22} &=2 \left( \cfrac{\partial v_2}{\partial x_2} \right)\\ &=2 \left( \cfrac{\Delta_{dc} \ v_2^{(e \rightarrow v)}} {\Delta_{dc} \ x_2} \right), \end{split} \tag{3.11}$$

$$\begin{split} D_{23} &= \left( \cfrac{\partial v_2}{\partial x_3} + \cfrac{\partial v_3}{\partial x_2} \right) \\ &= \left( \cfrac{\Delta_{x_3} \ v_2^{(f \rightarrow i)}} {\Delta \ x_3} +\cfrac{\Delta_{dc} \ v_3^{(c \rightarrow v)(i \rightarrow f)}} {\Delta_{dc} \ x_2} \right), \end{split} \tag{3.12}$$

$$D_{31} = D_{13}, \tag{3.13}$$

$$D_{32} = D_{23}, \tag{3.14}$$

$$\begin{split} D_{33} &=2 \left( \cfrac{\partial v_3}{\partial x_3} \right)\\ &=2 \left( \cfrac{\Delta_{x_3} \ v_3^{(c \rightarrow e)}} {\Delta \ x_3} \right), \end{split} \tag{3.15}$$

$$D^2 = D_{11}^2 + D_{22}^2 + D_{33}^2 + 2 \left( D_{12}^2 + D_{13}^2 + D_{23}^2 \right), \tag{3.16}$$

from D15 Eqs. (20-28).

in  subroutine  smagorinsky_model  of  mo_sgs_turbulence.f90
The subgrid-scale stress tensor is parameterized $$\tau_{ij} = K_m \left( D_{ij} - \cfrac{1}{3} \ \delta_{ij} \sum_{m=1}^{3} D_{mm} \right) \tag{3.17}$$

Figure 1: (left) Schematic showing the primal (black triangles) and dual (red hexagons) cells, and the associated local coordinate system. Unit vectors 1, 2, and 3 point in the direction of edge normal, tangent, and vertically upward, respectively. (right) Schematic of two adjacent triangles in ICON grid identifying the various locations used to discretize the turbulent diffusion term. The figure is taken from D15.

Figure 2: The primal and dual cells and edges, and the associated unit vectors and areas. The figure is taken from Wan, 2009.

There is nothing to describe. Surface sensible heat flux ($H$), latent heat flux, and friction velocity ($u_*$) are calculated in simple surface layer scheme.
in  subroutine  surface_condition  of  mo_surface_les.f90

Detail description of momentum and temperature (and scalar) tendencies are described in this section. The tendency in ICON-LEM is derived from horizontal diffusion terms and a vertical diffusion term. The horizontal diffusion is calculated explicitly whereas the vertical diffusion is calculated implicitly. Calculation processes for $v_1$, $v_3$, and temperature are different since $v_1$ is stored at triangle edge and full level, $v_3$ is stored at triangle center and interface level, and temperature is stored at triangle center and full level. Thus, I separate this section into three subsections.

Total tendency:

$$\left( \cfrac{\partial v_1}{\partial t} \right) = \left( \cfrac{\partial v_1}{\partial t} \right)_{\textrm{hori,turb}} + \left( \cfrac{\partial v_1}{\partial t} \right)_{\textrm{vert,turb}}, \tag{3.18}$$

Horizontal tendency: $$\begin{split} \left( \cfrac{\partial v_1}{\partial t} \right)_{\textrm{hori,turb}} & = \cfrac{1}{\rho} \left( \cfrac{\partial \tau_{11}}{\partial x_1} + \cfrac{\partial \tau_{12}}{\partial x_2} \right) \\ & = \cfrac{1}{\rho_k^{(c \rightarrow e)}} \left( \cfrac{\Delta_{qp} \ \tau_{11}}{\Delta_{qp} \ x_1} + \cfrac{\Delta_{gf} \ \tau_{12}}{\Delta_{gf} \ x_2} \right), \end{split} \tag{3.19}$$

$$\Delta_{qp} \ \tau_{11} = K_m^{(i \rightarrow f)} \left[ 2 \left( \cfrac{\Delta_{eb} \ v_1}{\Delta_{eb} \ x_1} \right) - \cfrac{2}{3} \textrm{Div.} \right] - K_m^{(i \rightarrow f)} \left[ 2 \left( \cfrac{\Delta_{ae} \ v_1}{\Delta_{ae} \ x_1} \right) - \cfrac{2}{3} \textrm{Div.} \right], \tag{3.20}$$

from D15 Eq. (B1) where $\textrm{Div.}=\frac{1}{2} (D_{11}+D_{22}+D_{33})$ and

$$\Delta_{gf} \ \tau_{12} = K_m^{(c \rightarrow v)(i \rightarrow f)} \left[ \cfrac{\Delta_{ec} \ v_1}{\Delta_{ec} \ x_2} + \cfrac{\Delta_{ab} \ v_2}{\Delta_{ab} \ x_1} \right] - K_m^{(c \rightarrow v)(i \rightarrow f)} \left[ \cfrac{\Delta_{de} \ v_1}{\Delta_{de} \ x_2} + \cfrac{\Delta_{ab} \ v_2}{\Delta_{ab} \ x_1} \right], \tag{3.21}$$

from D15 Eq. (B2).
in  subroutine  diffuse_hori_velocity  of  mo_sgs_turbulence.f90

Vertical tendency (implicit method): $$\begin{split} \left( \cfrac{\partial v_1}{\partial t} \right)_{\textrm{vert,turb}} & = \cfrac{v_1^{(t+\Delta t)} - v_1^{(t)}}{\Delta t} \\ & = \cfrac{1}{\rho} \left( \cfrac{\partial \tau_{13}}{\partial x_3} \right) \\ & = \cfrac{1}{\rho_k^{(c \rightarrow e)}} \left( \cfrac{\Delta_{x_3} \ \tau_{13}}{\Delta \ x_3} \right), \end{split} \tag{3.22} \label{eq:3.22} $$

$$\Delta_{x_3} \tau_{13} = K_{m,k-1/2}^{(c \rightarrow e)} \left[ \cfrac{\Delta_{x_3} v_1}{\Delta x_3} + \cfrac{\Delta_{qp} \ v_3}{\Delta_{qp} \ x_1} \right] - K_{m,k+1/2}^{(c \rightarrow e)} \left[ \cfrac{\Delta_{x_3} v_1}{\Delta x_3} + \cfrac{\Delta_{qp} \ v_3}{\Delta_{qp} \ x_1} \right], \tag{3.23}$$

where $v_1^{(t+\Delta t)}$ is $v_1$ at time $t+\Delta t$. For vertical tendency, implicit method is used. It is noticeable that the implicitness factor, $\alpha$, is different from that of ICON-AES as $\alpha=1$ is used here.

After discretization, Eq. (\ref{eq:3.22}) can be rewritten as

$$A_{k+1/2} \ v_{1,k-1}^{(t+\Delta t)} + B_{k+1/2} \ v_{1,k}^{(t+\Delta t)} + C_{k+1/2} \ v_{1,k+1}^{(t+\Delta t)} = \cfrac{v_{1,k}^{(t)}}{\Delta t} + \Omega_k, \tag{3.24}$$

where $A_{k+1/2}$, $B_{k+1/2}$, $A_{k+1/2}$, and $\Omega_k$ are given by

$$A_{k+1/2} = - \ \cfrac{1}{\rho_k^{(c \rightarrow e)} \Delta x_3} \cfrac{K_{m,k-1/2}^{(c \rightarrow e)}}{z_{k-1} - z_{k}}, \tag{3.25}$$

$$C_{k+1/2} = - \ \cfrac{1}{\rho_k^{(c \rightarrow e)} \Delta x_3} \cfrac{K_{m,k+1/2}^{(c \rightarrow e)}}{z_{k} - z_{k+1}}, \tag{3.26}$$

$$B_{k+1/2} = \cfrac{1}{\Delta t} - A_{k+1/2} - C_{k+1/2}, \tag{3.27}$$

$$\Omega_k = \cfrac{1}{\rho_k^{(c \rightarrow e)} \Delta x_3} \left[ K_{m,k-1/2}^{(c \rightarrow e)} \cfrac{\Delta_{qp} \ v_{3,k-1/2}}{\Delta_{qp} \ x_1} - K_{m,k+1/2}^{(c \rightarrow e)} \cfrac{\Delta_{qp} \ v_{3,k+1/2}}{\Delta_{qp} \ x_1} \right], \tag{3.28}$$

$i$) for $k=1$
$K_{m,k-1/2}$ should be rewritten as $$K_{m,k-1/2}^{(c \rightarrow e)} = 0, \tag{3.29}$$
$ii$) for $k= \textrm{nlev}$
$K_{m,k+1/2}$ and $\Omega$ should be rewritten as $$K_{m,k+1/2}^{(c \rightarrow e)} = 0, \tag{3.30}$$

and

$$\Omega_k = \cfrac{1}{\rho_k^{(c \rightarrow e)} \Delta x_3} \left[ K_{m,k-1/2}^{(c \rightarrow e)} \cfrac{\Delta_{qp} \ v_{3,k-1/2}}{\Delta_{qp} \ x_1} \right] - \cfrac{u_*^2}{\Delta z_{k-1/2}}, \tag{3.31}$$ where $u_*$ is friction velocity. in  subroutine  diffuse_hori_velocity  of  mo_sgs_turbulence.f90

Firstly, $E_{k+1/2}$ and $F_{k+1/2}$ are computed from $k=1$ to $k= \textrm{nlev}$ as

$$E_{k+1/2} = \cfrac{C_{k+1/2}} {B_{k+1/2} - A_{k+1/2} \ E_{k-1/2}}, \tag{3.32}$$

and

$$F_{k+1/2} = \cfrac{ v_{1,k}^{(t)} / \Delta t + \Omega_{k} - A_{k+1/2} F_{k-1/2} } {B_{k+1/2} - A_{k+1/2} \ E_{k-1/2}}. \tag{3.33}$$
Then, $v_{1,k}^{(t+\Delta t)}$ is computed from $k= \textrm{nlev}$ to $k=1$ as

$$v_{1,k}^{(t+\Delta t)} = - E_{k+1/2} \ v_{1,k+1}^{(t+\Delta t)} + F_{k+1/2}. \tag{3.34}$$
in  subroutine  tdma_solver  of  mo_math_utilities.f90  

Finally, the tendency is calculated as $\left( \frac{\partial v_1}{\partial t} \right)_{\textrm{vert,turb}} = \frac{v_1^{(t+\Delta t)} - v_1^{(t)}}{\Delta t}$.

Total tendency: $$\left( \cfrac{\partial v_3}{\partial t} \right) = \left( \cfrac{\partial v_3}{\partial t} \right)_{\textrm{hori,turb}} + \left( \cfrac{\partial v_3}{\partial t} \right)_{\textrm{vert,turb}}, \tag{3.35}$$

Horizontal tendency: $$\begin{split} \left( \cfrac{\partial v_3}{\partial t} \right)_{\textrm{hori,turb}} & = \cfrac{1}{\rho} \left( \cfrac{\partial \tau_{31}}{\partial x_1} + \cfrac{\partial \tau_{32}}{\partial x_2} \right) \\ & = \cfrac{1}{\rho_k^{(f \rightarrow i)}} \left( \cfrac{\Delta_{qp} \ \tau_{31}}{\Delta_{qp} \ x_1} + \cfrac{\Delta_{gf} \ \tau_{32}}{\Delta_{gf} \ x_2} \right), \end{split} \tag{3.36}$$

$$\Delta_{qp} \tau_{31} = K_m \left[ \cfrac{v_{1,k-1}^{(e \rightarrow c(p))} - v_{1,k}^{(e \rightarrow c(p))}} {z_{k-1}-z_k} + \cfrac{\Delta_{eb} \ v_3 }{\Delta_{eb} \ x_1} \right] - K_m \left[ \cfrac{v_{1,k-1}^{(e \rightarrow c(q))} - v_{1,k}^{(e \rightarrow c(q))}} {z_{k-1}-z_k} + \cfrac{\Delta_{ae} \ v_3 }{\Delta_{ae} \ x_1} \right], \tag{3.37}$$

$$\Delta_{gf} \tau_{32} = K_m^{(c \rightarrow v)} \left[ \cfrac{v_{2,k-1}^{(e \rightarrow v(c))} - v_{2,k}^{(e \rightarrow v(c))} } {z_{k-1} - z_k} + \cfrac{\Delta_{ec} \ v_3}{\Delta_{ec} \ x_2} \right] - K_m^{(c \rightarrow v)} \left[ \cfrac{v_{2,k-1}^{(e \rightarrow v(d))} - v_{2,k}^{(e \rightarrow v(d))} } {z_{k-1} - z_k} + \cfrac{\Delta_{de} \ v_3}{\Delta_{de} \ x_2} \right], \tag{3.38}$$
where $c(p)$ is interpolation to triangle center at $p$, $c(q)$ is interpolation to triangle center at $q$, $v(c)$ is interpolation to triangle vertex at $c$, and $v(d)$ is interpolation to triangle vertex at $d$.

Vertical tendency (implicit method): $$\begin{split} \left( \cfrac{\partial v_3}{\partial t} \right)_{\textrm{vert,turb}} & = \cfrac{v_{3,k-1/2}^{(t+\Delta t)} - v_{3,k-1/2}^{(t)}}{\Delta t} \\ & = \cfrac{1}{\rho} \left( \cfrac{\partial \tau_{33}}{\partial x_3} \right) \\ & = \cfrac{1}{\rho_k} \left( \cfrac{\tau_{33,k-1} - \tau_{33,k}} {z_{k-1} - z_{k}} \right), \end{split} \tag{3.39} \label{eq:3.39}$$

$$\tau_{33,k-1} - \tau_{33,k} = K_{m,k-1}^{(i \rightarrow f)} \left[ 2\left( \cfrac{v_{3,k-3/2} -v_{3,k-1/2}}{z_{k-3/2} - z_{k-1/2}} \right) - \cfrac{2}{3} \textrm{Div}_{k-1} \right] - K_{m,k}^{(i \rightarrow f)} \left[ 2 \left( \cfrac{v_{3,k-1/2} -v_{3,k+1/2}}{z_{k-1/2}-z_{k+1/2}} \right) - \cfrac{2}{3} \textrm{Div}_{k} \right]. \tag{3.40}$$

After discretization, Eq. (\ref{eq:3.39}) can be rewritten as

$$A_k \ v_{3,k-3/2}^{(t+\Delta t)} + B_k \ v_{3,k-1/2}^{(t+\Delta t)} + C_k \ v_{3,k+1/2}^{(t+\Delta t)} = \cfrac{v_{3,k-1/2}^{(t)}}{\Delta t} + \Omega_{k-1/2}, \tag{3.41}$$

where $A_k$, $B_k$, $C_k$, and $\Omega_{k-1/2}$ are computed as

$$A_k = -2 \cfrac{1}{\rho_k (z_{k-1} - z_k)} \cfrac{K_{m,k-1}^{(i \rightarrow f)}}{\Delta x_3}, \tag{3.42}$$

$$C_k = -2 \cfrac{1}{\rho_k (z_{k-1} - z_k)} \cfrac{K_{m,k}^{(i \rightarrow f)}}{\Delta x_3}, \tag{3.43}$$

$$B_k = \cfrac{1}{\Delta t} - A_k - C_k, \tag{3.44}$$

$$\Omega_{k-1/2} = - \cfrac{2}{3} K_{m,k-1}^{(i \rightarrow f)} \textrm{Div}_{k-1} + \cfrac{2}{3} K_{m,k}^{(i \rightarrow f)} \textrm{Div}_{k}, \tag{3.45}$$

$i$) for $k=2$
$A_k$ and $B_k$ should be rewritten as $$A_k = 0,\tag{3.46}$$

and

$$B_k = \cfrac{1}{\Delta t} + 2 \cfrac{1}{\rho_k (z_{k-1} - z_k)} \cfrac{K_{m,k-1}^{(i \rightarrow f)}}{\Delta x_3} - C_k, \tag{3.47}$$

$ii$) for $k= \textrm{nlev}$
$B_k$ and $C_k$ should be rewritten as $$C_k = 0, \tag{3.48}$$

and

$$B_k = \cfrac{1}{\Delta t} - A_k + 2 \cfrac{1}{\rho_k (z_{k-1} - z_k)} \cfrac{K_{m,k}^{(i \rightarrow f)}}{\Delta x_3}.\tag{3.49}$$
in  subroutine  diffuse_vert_velocity  of  mo_sgs_turbulence.f90  

Firstly, $E_{k}$ and $F_{k}$ are computed from $k=1$ to $k= \textrm{nlev}$ as

$$E_{k} = \cfrac{C_{k}} {B_{k} - A_{k} \ E_{k-1}},\tag{3.50}$$

and

$$F_{k} = \cfrac{ v_{3,k-1/2}^{(t)} / \Delta t + \Omega_{k-1/2} - A_{k} F_{k-1} } {B_{k+1} - A_{k+1} \ E_{k-1}}. \tag{3.51}$$
Then, $v_{1,k-1/2}^{(t+\Delta t)}$ is computed from $k= \textrm{nlev}$ to $k=1$ as

$$v_{3,k-1/2}^{(t+\Delta t)} = - E_{k} \ v_{3,k+1/2}^{(t+\Delta t)} + F_{k}.\tag{3.52}$$
in  subroutine  tdma_solver  of  mo_math_utilities.f90

Finally, the tendency is calculated by $\left( \frac{\partial v_3}{\partial t} \right)_{\textrm{vert,turb}} = \frac{v_{3,k-1/2}^{(t+\Delta t)} - v_{3,k-1/2}^{(t)}}{\Delta t}$.

Total tendency: $$\left( \cfrac{\partial T}{\partial t} \right) = \left( \cfrac{\partial T}{\partial t} \right)_{\textrm{hori,turb}} + \left( \cfrac{\partial T}{\partial t} \right)_{\textrm{vert,turb}}, \tag{3.53}$$

Horizontal tendency: $$\begin{split} \left( \cfrac{\partial T}{\partial t} \right)_{\textrm{hori,turb}} & = - \left( \cfrac{\partial \overline{v_1'T'}}{\partial x_1} + \cfrac{\partial \overline{v_2'T'}}{\partial x_2} \right) \\ & = \left[ \nabla_{\textrm{tri}} \cdot \left( k_h^{(c \rightarrow e)} \cfrac{\Delta_{qp} T_k}{\Delta_{qp} x_1} \right) \right] \phi, \end{split} \tag{3.54}$$

where $\nabla_{\textrm{tri}}$ is divergence at center of the triangular grid, and $\phi= \frac{c_p}{c_v}$ is coupler that converts the temperature tendency for the dynamics.

Vertical tendency (implicit method): $$ \begin{split} \left( \cfrac{\partial T}{\partial t} \right)_{\textrm{vert,turb}} & = \cfrac{T_{k}^{(t+\Delta t)} - T_{k}^{(t)}}{\Delta t} \\ & = - \cfrac{\partial \overline{v_3'T'}}{\partial x_3} \\ & = \left[ K_{h,k-1/2} \cfrac{T_{k-1}-T_{k}}{z_{k-1} - z_{k}} - K_{h,k+1/2} \cfrac{T_{k} - T_{k+1}}{z_{k} - z_{k+1}} \right] \phi. \end{split} \tag{3.55} \label{eq:3.55}$$

After discretization, Eq. (\ref{eq:3.55}) can be rewritten as

$$A_{k+1/2} \ T_{k-1}^{(t+\Delta t)} + B_{k+1/2} \ T_{k}^{(t+\Delta t)} + C_{k+1/2} \ T_{k+1}^{(t+\Delta t)} = \cfrac{T_{k}^{(t)}}{\Delta t} + \Omega_{k} \tag{3.56}$$

where $A_{k+1/2}$, $B_{k+1/2}$, $C_{k+1/2}$, and $\Omega_{k}$ are obtained as

$$A_{k+1/2} = - \cfrac{1}{\Delta x_3} \cfrac{K_{h,k-1/2}}{z_{k-1} - z_{k}} \ \phi, \tag{3.57}$$

$$C_{k+1/2} = - \cfrac{1}{\Delta x_3} \cfrac{K_{h,k+1/2}}{z_{k} - z_{k+1}} \ \phi, \tag{3.58}$$

$$B_{k+1/2} = \cfrac{1}{\Delta t} - A_{k+1/2} - C_{k+1/2}, \tag{3.59}$$

$$\Omega_k = 0, \tag{3.60}$$

$i$) for $k=1$
$A_{k+1/2}$ should be rewritten as $$A_{k+1/2} = 0, \tag{3.61}$$
$ii$) for $k= \textrm{nlev}$
$C_{k+1/2}$ and $\Omega_k$ should be rewritten as $$C_{k+1/2} = 0, \tag{3.62}$$

and

$$\Omega_k = \cfrac{1}{\Delta x_3} H \phi, \tag{3.63}$$

where $H$ is surface sensible heat flux [$K m s^{-1}$]. in  subroutine  diffuse_scalar  of  mo_sgs_turbulence.f90

Firstly, $E_{k+1/2}$ and $F_{k+1/2}$ are computed from $k=1$ to $k= \textrm{nlev}$ as
$$E_{k+1/2} = \cfrac{C_{k+1/2}} {B_{k+1/2} - A_{k+1/2} \ E_{k-1/2}},\tag{3.64}$$

and

$$F_{k+1/2} = \cfrac{ T_{k}^{(t)} / \Delta t + \Omega_{k} - A_{k+1/2} F_{k-1/2} } {B_{k+1/2} - A_{k+1/2} \ E_{k-1/2}}. \tag{3.65}$$

Then, $v_{1,k}^{(t+\Delta t)}$ is computed from $k= \textrm{nlev}$ to $k=1$ as

$$T_{k}^{(t+\Delta t)} = - E_{k+1/2} \ T_{k+1}^{(t+\Delta t)} + F_{k+1/2}. \tag{3.66}$$
in  subroutine  tdma_solver  of  mo_math_utilities.f90

Finally, the tendency is calculated by $\left( \frac{\partial T}{\partial t} \right)_{\textrm{vert,turb}} = \frac{T_{k}^{(t+\Delta t)} - T_{k}^{(t)}}{\Delta t}$.

Angevine, W. M., Jiang, H., & Mauritsen, T. (2010). Performance of an eddy diffusivity–mass flux scheme for shallow cumulus boundary layers. Monthly Weather Review, 138(7), 2895-2912.

Brinkop, S., & Roeckner, E. (1995). Sensitivity of a general circulation model to parameterizations of cloud–turbulence interactions in the atmospheric boundary layer. Tellus A, 47(2), 197-220.

Dipankar, A., Stevens, B., Heinze, R., Moseley, C., Zängl, G., Giorgetta, M., & Brdar, S. (2015). Large eddy simulation using the general circulation model ICON. Journal of Advances in Modeling Earth Systems, 7(3), 963-986.

Louis, J. F. (1979). A parametric model of vertical eddy fluxes in the atmosphere. Boundary-Layer Meteorology, 17(2), 187-202.

Mauritsen, T., Svensson, G., Zilitinkevich, S. S., Esau, I., Enger, L., & Grisogono, B. (2007). A total turbulent energy closure model for neutrally and stably stratified atmospheric boundary layers. Journal of the Atmospheric Sciences, 64(11), 4113-4126.

Pithan, F. (2014). Arctic boundary-layer processes and climate change (Doctoral dissertation, Universität Hamburg, Hamburg).

Wan, H. (2009). Developing and testing a hydrostatic atmospheric dynamical core on triangular grids, (Doctoral dissertation, Universität Hamburg, Hamburg).

Wan, H. (2011). Technical documentation on the turbulent mixing parameterization code VDIFF in ICOHAM, (Internal document, Max Planck Institute for Meteorology, Hamburg).

Zilitinkevich, S. S., Elperin, T., Rogachevskii, I. (2007). Energy and flux-budget (EFB) turbulence closure model for stably stratified flows. Part I: Steady-state, homogenous regimes. Boundary-Layer Meteorology, 125, 167-192.