%\documentclass{article}
%\usepackage{amsmath}
%\begin{document}

\section{Boundary layer turbulence model}
{\bf \Large
\begin{tabular}{ccc}
\hline
  Corresponding author & : & Seiya Nishizawa\\
\hline
\end{tabular}
}

\def\half{\frac{1}{2}}

\subsection{Mellor-Yamada Nakanishi-Niino model}
level 2.5

\begin{align}
  \frac{\partial \rho u}{\partial t}
  &= -\frac{\partial}{\partial z} \rho \overline{u'w'}, \\
  \frac{\partial \rho v}{\partial t}
  &= -\frac{\partial}{\partial z} \rho \overline{v'w'}, \\
  \frac{\partial \rho \theta_l}{\partial t}
  &= -\frac{\partial}{\partial z} \rho \overline{\theta_l'w'}, \\
  \frac{\partial \rho q_a}{\partial t}
  &= -\frac{\partial}{\partial z} \rho \overline{q_a'w'}, \\
  \frac{\partial }{\partial t}\rho q^2
  &= -2\left(\rho\overline{u'w'}\frac{\partial u}{\partial z}+\rho\overline{v'w'}\frac{\partial v}{\partial z}\right)
  +2\frac{g}{\theta_0}\rho\overline{\theta_v' w'}
  -\frac{\partial}{\partial z}\rho\overline{q^2w'}
  -2\rho\epsilon, \label{eq: q2}
\end{align}
where:
\begin{equation}
  q_a = q_v + q_c + q_r + q_i + q_s + q_g,
\end{equation}
and $q^2$ is doubled turbulence kinetic energy:
\begin{equation}
  q^2 = u'^2 + v'^2 + w'^2.
\end{equation}

The higher order moments and the dissipation term are parameterized as follows:
\begin{align}
  \overline{u'w'} &= -LqS_M\frac{\partial u}{\partial z}, \\
  \overline{v'w'} &= -LqS_M\frac{\partial v}{\partial z}, \\
  \overline{\theta_l'w'} &= -LqS_H\frac{\partial \theta_l}{\partial z}, \\
  \overline{q_a'w'} &= -LqS_H\frac{\partial q_a}{\partial z}, \\
  \overline{q^2w'} &= -3LqS_M\frac{\partial q^2}{\partial z}, \\
  \overline{\theta_v' w'} &= \beta_\theta\overline{\theta_l'w}'+\beta_q\overline{q_a'w'}, \\
  \epsilon &= \frac{q^3}{B_1L},
\end{align}
where:
\begin{align}
  S_M &= \alpha_cA_1\frac{\Phi_3-3C_1\Phi_4}{D_{2.5}}, \\
  S_H &= \alpha_cA_2\frac{\Phi_2+3C_1\Phi_5}{D_{2.5}}, \\
  \beta_\theta &= 1 + 0.61 q_a - 1.61 Q_l - \tilde{R} abc, \\
  \beta_q &= 0.61\theta + \tilde{R} ac.
\end{align}

\begin{align}
  D_{2.5} &= \Phi_2\Phi_4 + \Phi_5\Phi_3, \\
  \Phi_1 &= 1-3\alpha_c^2A_2B_2(1-C_3)G_H, \\
  \Phi_2 &= 1-9\alpha_c^2A_1A_2(1-C_2)G_H, \\
  \Phi_3 &= \Phi_1+9\alpha_c^2A_2^2(1-C_2)(1-C_5)G_H, \\
  \Phi_4 &= \Phi_1-12\alpha_c^2A_1A_2(1-C_2)G_H, \\
  \Phi_5 &= 6\alpha_c^2A_1^2G_M, \\
  \alpha_c &= \left\{
  \begin{array}{ll}
    q/q_2, & q<q_2 \\
    1, & q \ge q_2
  \end{array}
  \right. , \\
  G_M &= \frac{L^2}{q^2}\left\{\left(\frac{\partial u}{\partial z}\right)^2 + \left(\frac{\partial v}{\partial z}\right)^2\right\}, \\
  G_H &= -\frac{L^2}{q^2}N^2, \\
  R &= \frac{1}{2}\left\{1+\mathrm{erf}\left(\frac{Q_1}{\sqrt{2}}\right)\right\}, \\
  \tilde{R} &= R - \frac{Q_l}{2\sigma_s}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{q_1^2}{2}\right), \\
  Q_l &= 2\sigma_s\left\{RQ_1+\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{Q_1^2}{2}\right)\right\}, \\
  Q_1 &= \frac{a}{2\sigma_s}(q_a-Q_{sl}), \\
  \sigma_s^2 &= \frac{1}{4}a^2L^2\alpha_cB_2S_H\left(\frac{\partial q_a}{\partial z} -b\frac{\partial \theta_l}{\partial z}\right)^2, \\
  \delta Q_{sl} &= \left.\frac{\partial Q_s}{\partial T}\right|_{T=T_l}, \\
  a &= \left(1+\frac{L}{C_p}\delta Q_{sl}\right)^{-1}, \\
  b &= \frac{T}{\theta}\delta Q_{sl}, \\
  c &= (1+0.61q_a - 1.61Q_l)\frac{\theta}{T}\frac{L_v}{C_p} - 1.61\theta,
\end{align}
and $Q_{sl}$ is the saturation-specific humidity at temperature $T_l (=\theta_l T/\theta)$.

The buoyancy flux term, which is the third term on the left hand side of eq. \ref{eq: q2} is:
\begin{align}
  2\frac{g}{\theta_0}\overline{\theta_v' w'}
  &= 2\frac{g}{\theta_0}\left(-\beta_\theta LqS_H\frac{\partial \theta_l}{\partial z} - \beta_q LqS_H\frac{\partial q_a}{\partial z}\right) \nonumber \\
  &= -2LqS_H\frac{g}{\theta_0}\left(\beta_\theta\frac{\partial \theta_l}{\partial z}+\beta_q\frac{\partial q_a}{\partial z}\right) \nonumber \\
  &= -2LqS_H\frac{g}{\theta_0}\frac{\partial \theta_v}{\partial z} \nonumber \\
  &= -2LqS_HN^2,
\end{align}
where $N^2$ is the square of the Brunt-Vaisala frequency.

\begin{align}
  \frac{\partial }{\partial t}\rho q^2
  &= 2\rho LqS_M\left\{\left(\frac{\partial u}{\partial z}\right)^2
                +\left(\frac{\partial v}{\partial z}\right)^2\right\} \nonumber\\
  &-2\rho LqS_HN^2
  +\frac{\partial}{\partial z}\left(3\rho LqS_M\frac{\partial}{\partial z}q^2\right)
  -2\rho\frac{q^3}{B_1L}
\end{align}

$S_{M2}, S_{H2}$, and $q_2$ are for level 2 schemes corresponding to $S_M, S_H$, and $q$, respectively:
\begin{align}
  S_{M2} &= \frac{A_1F_1}{A_2F_2}\frac{R_{f1}-Rf}{R_{f2}-Rf} S_{H2}, \\
  S_{H2} &= 3 A_2 (\gamma_1 + \gamma_2) \frac{Rf_c - Rf}{1-Rf}, \\
  q_2^2 &= B_1 L^2 S_{M2} (1-Rf) \left\{\left(\frac{\partial u}{\partial z}\right)^2+\left(\frac{\partial v}{\partial z}\right)^2\right\}.
\end{align}
$Rf$ and $Rf_c$ are the flux Richardson number and the critical flux Richardson number, respectively.
The gradient Richardson number, $Ri$, is:
\begin{equation}
  Ri = Rf \frac{S_{M2}}{S_{H2}}.
\end{equation}
$Rf$ is then:
\begin{align}
  Rf &= \frac{1}{2}\frac{A_2F_2}{A_1F_1}
  \left\{ Ri + \frac{A_1F_1}{A_2F_2}R_{f1}
        -\sqrt{Ri^2+2\frac{A_1F_1}{A_2F_2}(R_{f1}-2R_{f2})Ri+\left(\frac{A_1F_1}{A_2F_2}R_{f1}\right)^2} \right\}, \\
  Rf_C &= \frac{\gamma_1}{\gamma_1+\gamma_2}, \\
\end{align}
where:
\begin{align}
  R_{f1} &= B_1\frac{\gamma_1-C_1}{F_1}, \\
  R_{f2} &= B_1\frac{\gamma_1}{F_2}.
\end{align}

The turbulent length scale, $L$, is determined by the smallest length scale among three scales:
\begin{equation}
  \frac{1}{L} = \frac{1}{L_s} + \frac{1}{L_T} + \frac{1}{L_B}.
\end{equation}
the surface layer scale, $L_s$, the boundary layer scale, $L_T$, and buoyancy length scale, $T_B$:
\begin{align}
  L_S &= \left\{
  \begin{array}{ll}
    kz/3.7, & \zeta \ge 1 \\
    kz/(1+2.7\zeta), & 0 \le \zeta < 1 \\
    kz(1-100\zeta)^{0.2}, & \zeta < 0
  \end{array}
  \right. , \\
  L_T &= 0.23\frac{\int_0^\infty qz dz}{\int_0^\infty q dz}, \\
  L_B &= \left\{
  \begin{array}{ll}
    q/N,                      & \partial \theta_v/\partial z > 0 \;\mathrm{and}\; \zeta \ge 0 \\
    \{1+5(q_c/L_TN)^{1/2}\}q/N, & \partial \theta_v/\partial z > 0 \:\mathrm{and}\; \zeta < 0 \\
    \infty, & \partial \theta_v/\partial z \le 0
  \end{array}
  \right. ,
\end{align}
where $\zeta$ is the dimensionless height:
\begin{equation}
  \zeta = \frac{z}{L_M}.
\end{equation}
$L_M$ is the Monin-Obukhov length:
\begin{equation}
  L_M = -\frac{\theta_0 u_*^3}{kg\overline{\theta_v'w'}_g},
\end{equation}
where $u_*$ is the friction velocity, and the subscript $g$ denotes the ground surface.
$q_c$ is a velocity scale defined in a similar manner to convective velocity $w_*$, except that the depth $z_i$ of the convective boundary layer is replaced by $L_t$:
\begin{equation}
  q_c = \left\{\frac{g}{\theta_0}\overline{\theta_v' w'}_gL_T\right\}^{1/3}
\end{equation}


\begin{align}
  A_1 &= B_1\frac{1-3\gamma_1}{6}, \\
  A_2 &= \frac{1}{3\gamma_1 B_1^{1/3} Pr_N}, \\
  B_1 &= 24.0, \\
  B_2 &= 15.0, \\
  C_1 &= \gamma_1 - \frac{1}{3A_1B_1^{1/3}}, \\
  C_2 &= 0.75, \\
  C_3 &= 0.352, \\
  C_5 &= 0.2, \\
  \gamma_1 &= 0.235, \\
  \gamma_2 &= \frac{2A_1(3-2C_2)+B_2(1-C_3)}{B_1}, \\
  F_1 &= B_1(\gamma_1-C_1)+2A_1(3-2C_2)+3A_2(1-C_2)(1-C_5), \\
  F_2 &= B_1(\gamma_1+\gamma_2)-3A_1(1-C_2), \\
  Pr_N &= 0.74.
\end{align}


\subsubsection{Discretization}
The diffusion equations for $q^2a$ are solved implicitly:
\begin{align}
  \rho_k \frac{(q^2_k)^{n+1}-(q^2_k)^n}{\Delta t}
  &= 
  2\rho_k \left[ (LqS_M)_k\left\{\left(\frac{\partial u}{\partial z}\right)^2+\left(\frac{\partial v}{\partial z}\right)^2\right\} + (LqS_HN^2)_k \right] \nonumber \\
  &+ \frac{1}{\Delta z_k}\left\{ (3\rho LqS_M)_{k+\half} \frac{(q^2_{k+1})^{n+1}-(q^2_k)^{n+1}}{\Delta z_{k+\half}} - (3\rho LqS_M)_{k-\half} \frac{(q^2_k)^{n+1}-(q^2_{k-1})^{n+1}}{\Delta z_{k-\half}} \right\} \nonumber \\
  &-\frac{2\rho_k q_k}{B_1L_k}(q^2_k)^{n+1}.
\end{align}
\begin{equation}
  a_k (q^2_{k+1})^{n+1} + b_k (q^2_k)^{n+1} + c_k (q^2_{k-1})^{n+1} = d_k,
\end{equation}
where:
\begin{align}
  a_k &= -\frac{\Delta t}{\Delta z_{k+\half}\Delta z_k\rho_k}(3\rho LqS_M)_{k+\half}, \\
  b_k &= -a_k - c_k + 1 + \frac{2\Delta tq_k}{B_1L}, \\
  c_k &= -\frac{\Delta t}{\Delta z_k\Delta z_{k-\half}\rho_k}(3\rho LqS_M)_{k-\half}, \\
  d_k &= (q^2_k)^n + 2\Delta t \left[ LqS_M\left\{\left(\frac{\partial u}{\partial z}\right)^2+\left(\frac{\partial v}{\partial z}\right)^2\right\} - LqS_HN^2 \right]
\end{align}
\begin{equation}
  (q^2_k)^{n+1} = e_k (q^2_{k+1})^{n+1} + f_k,
\end{equation}
where:
\begin{align}
  e_k &= -\frac{a_k}{b_k+c_ke_{k-1}}, \\
  f_k &= \frac{d_k-c_kf_{k-1}}{b_k+c_ke_{k-1}}.
\end{align}

Vertical fluxes for $\rho u, \rho v, \rho\theta, \rho q_x$ are also solved implicitly.
For instance, the flux for $\rho u$, $F_u$ is calculated by:
\begin{equation}
  F_{u,k+\half} = (\rho LqSM)_{k+\half}\frac{u^{n+1}_{k+1}-u^{n+1}_k}{\Delta z_{k+\half}}.
\end{equation}
$u^{n+1}$ is calculated as the same way with $q^2$, but:
\begin{align}
  a_k &= -\frac{\Delta t}{\Delta z_{k+\half}\Delta z_k\rho_k}(\rho LqS_M)_{k+\half}, \\
  b_k &= -a_k - c_k + 1, \\
  c_k &= -\frac{\Delta t}{\Delta z_k\Delta z_{k-\half}\rho_k}(\rho LqS_M)_{k-\half}, \\
  d_k &= u_k^n.
\end{align}

%\end{document}

