%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 表題  2 次元非静力学モデル -- 離散モデル  時間方向の離散化
%  
% 履歴  2004/08/14  杉山耕一朗: 作成開始
%       2004/09/21  小高正嗣, 北守太一
%       2005/04/13  小高正嗣
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{時間方向の離散化}

\section{運動方程式と圧力方程式}

空間離散化された運動方程式\Deqref{空間離散化された x 方向運動方程式},
\Deqref{空間離散化された z 方向運動方程式}と圧力方程式
\Deqref{空間離散化された圧力方程式}を時間方向に離散化する. 
音波に関連する項は短いタイムステップ $\Delta \tau$ で離散化し, その他
の項は長いタイムステップ $\Delta t$ で離散化する. 音波に関連する項の離
散化には HE-VI 法を採用し, $u$ の式は前進差分, $w, \pi$ の式は後退差分
(クランク・ニコルソン法)で離散化する. その他の項の離散化にはリープフロッ
グ法を用いる. 離散化した式の計算はまず $u$ の式から行う. 得られた 
$\tau +\Delta \tau$ の $u$ を用いて $\pi$ を計算し, $u, \pi$ を用いて 
$w$ を計算する.


解く方程式を以下のように書き直す. 
%
\begin{eqnarray}
&& \DP{u_{i(u),k}}{t} = - \left[\bar{c_{p}} \bar{\theta}_{v} 
                       \DP{(\pi - \alpha Div )}{x}\right]_{i(u),k} +
  F_{u,i(u),k},
  \Deqlab{uwpi:u}\\
&& \DP{w_{i,k(w)}}{t} = - \left[\bar{c_{p}} \bar{\theta}_{v} 
                       \DP{(\pi - \alpha Div )}{z}\right]_{i,k(w)} +
  F_{w,i,k(w)},
  \Deqlab{uwpi:w}\\
&& \DP{\pi_{i,k}}{t}  
  + \left[\frac{\bar{c}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \DP{(\bar{\rho} \bar{\theta}_{v} w)}{z}\right]_{i,k}
  = 
  -  \left[\frac{\bar{c}^{2}}{\bar{c_{p}} \bar{\theta}_{v}} \DP{u}{x}\right]_{i,k}.
  \Deqlab{uwpi:pi}
\end{eqnarray}
%


ただし $u, w$ の式には音波減衰項(Skamarock and Klemp, 1992)を加えてある. 
ここで $F_{u}, F_{w}$ は音波に関連しない項, 
%
 \begin{eqnarray}
  F_{u,i(u),k}^{t} &=&
   - \left[Adv_{u}\right]_{i(u),k}^{t}
   + \left[D_{u}\right]_{i(u),k}^{t - \Delta t}
  \Deqlab{uwpi:Fu}
  \\
  F_{w,i,k(w)}^{t} &=&
   - \left[Adv_{w}\right]_{i(u),k}^{t}
   + g\frac{\theta_{i,k(w)}^{t}}{\overline{\theta}_{i,k(w)}}
   + \left[D_{w}\right]_{i,k(w)}^{t - \Delta t}.
   \Deqlab{uwpi:Fw}
 \end{eqnarray}
%
であり, 時刻 $t$ で評価する.

\subsection{水平方向の運動方程式}

\Deqref{uwpi:u}を時間方向に離散化すると以下のようになる. 
%
\begin{eqnarray}
 u^{\tau + \Delta \tau}_{i(u),k}
  =  u^{\tau}_{i(u), k}
  -   \left[
  \bar{c_{p}} \bar{\theta}_{v} \Delta \tau
  \left\{
   \DP{\pi^{\tau}}{x} 
   - \DP{(\alpha Div)^{\tau}}{x} 
  \right\}
 \right]_{i(u),k}
  + 
     F_{u,i(u),k}^{t} \Delta \tau
  \Deqlab{uwpi:u_sabun}
\end{eqnarray}
%

\subsection{鉛直方向の運動方程式と圧力方程式}

HE-VI 法を用いるので, $w$ と $\pi$ の式を連立して解く.  $w$ の式におい
て音波減衰項は前進差分, 圧力項は後退差分で離散化する.  $\pi$ の式にお
いて水平微分項は\Deqref{uwpi:u_sabun}で求めた $u^{\tau +\Delta \tau}$ 
を用いて離散化し, 鉛直微分項は後退差分で離散化する.
%
\begin{eqnarray}
 w^{\tau + \Delta \tau} 
  =  w^{\tau}
  - \bar{c_{p}} \bar{\theta}_{v} \Delta \tau
  \left\{
   \beta \DP{\pi^{\tau + \Delta \tau}}{z} 
   + (1 - \beta) \DP{\pi^{\tau}}{z} 
   - \DP{(\alpha Div)^{\tau}}{z} 
  \right\}
 + F_{w}^{t} \Delta \tau.
 \Deqlab{uwpi:w_sabun}
\end{eqnarray}
%
\begin{eqnarray}
% \frac{\pi^{\tau + \Delta \tau} - \pi^{\tau}}{\Delta \tau} 
%  + \frac{\bar{c}^{2}}{\bar{c_{p}} \bar{\rho} {\bar{\theta}_{v}}^{2}}
%  \left\{
%   \beta \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau + \Delta \tau})}{z}
%   + (1 - \beta) \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
%  \right\}
%  = 
%  -  \frac{\bar{c}^{2}}{\bar{c_{p}} \bar{\theta}_{v}} \DP{u^{\tau +
%  \Delta \tau}}{x}, 
%  \nonumber \\
%
 \pi^{\tau + \Delta \tau} 
 + \beta \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
   \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau + \Delta \tau})}{z}
  = 
  \pi^{\tau}
  - (1 - \beta) \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  -  \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} \DP{u^{\tau + \Delta \tau}}{x}.
 \nonumber \\ 
 \Deqlab{uwpi:pi_sabun}
\end{eqnarray}
%
ここでは簡単のため格子点位置を表す添字は省略した. 
%
\Deqref{uwpi:pi_sabun} 式に \Deqref{uwpi:w_sabun} を代入して $w^{\tau +
\Delta \tau}$ を消去する. 
%
\begin{eqnarray}
 \pi^{\tau + \Delta \tau} 
  &-& 
  \beta^{2} 
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} 
  \bar{\theta}_{v}^{2}}
  \DP{}{z}
  \left\{
   \left(\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}\right)
   \left(
     \DP{\pi^{\tau + \Delta \tau}}{z} 
   \right)
\right\}
  \nonumber \\
 &=& 
  \pi^{\tau}
  -(1 - \beta) 
 \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
   \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  -  \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} 
  \DP{u^{\tau + \Delta \tau}}{x}
  \nonumber \\
 &&- \beta  \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
   \left[
    \DP{}{z}
    \bar{\rho} \bar{\theta}_{v}
   \left\{
   w^{\tau}
  - \bar{c_{p}} \bar{\theta}_{v} \Delta \tau
  \left\{
    (1 - \beta) \DP{\pi^{\tau}}{z} 
   - \DP{(\alpha Div)^{\tau}}{z} 
  \right\}
 + F_{w}^{t} \Delta \tau   
\right\}
\right].
  \nonumber \\
\Deqlab{uwpi:sabun}
\end{eqnarray}
%
%この方程式を解いて $\pi^{\tau + \Delta \tau}$ を求めるわけだが, 左辺の空
%間微分を陽に書き下すことで, 簡単な行列式を導くことができる. 
\Deqref{uwpi:sabun} 式右辺を空間方向に離散化し, 
格子点位置を表す添字を付けて表すと以下のようになる
(計算の詳細は \Dchapref{appendix-a} 参照). 
%
\begin{eqnarray}
&&  \left\{
   - \beta^{2}
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k}
  \Dinv{\Delta z^{2}} 
  \left(
    \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
  \right)_{k(w)}
  \right\}
   \pi^{\tau + \Delta \tau}_{i,k+1}
\nonumber \\   
 && \hspace{10mm}+ \left[
  1 + \beta^{2}
  \left(
   \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k}
  \Dinv{\Delta z^{2}} 
   \left\{
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{k(w)}
   + 
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{k-1(w)}
   \right\} 
    \right]
 \pi^{\tau + \Delta \tau}_{i,k}
   \nonumber \\   
 && \hspace{10mm}+ \left\{
   - \beta^{2}
  \left(
   \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k}
  \Dinv{\Delta z^{2}} 
   \left(
    \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
       \right)_{k-1(w)}
   \right\}
   \pi^{\tau + \Delta \tau}_{i,k-1}
   \nonumber \\
 &&= 
  \pi^{\tau}_{i,k}
  - (1 - \beta)  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k}
  \left\{
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,k}
   -\left(
     \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} 
    \right)_{k}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x} 
  \right)_{i,k}
  \nonumber \\
 &&\hspace{2mm} - \beta  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k}
   \left[
    \DP{}{z}
    \left( \bar{\rho} \bar{\theta}_{v} \right)_{i,k(w)} 
   \left\{
   w^{\tau}_{i,k(w)}
\right. \right.
\nonumber \\
&& \hspace{10mm} 
\left. \left.
  - \left( \bar{c_{p}} \bar{\theta}_{v} \right)_{i,k(w)} \Delta \tau
  \left\{
    (1 - \beta)  \DP{\pi^{\tau}}{z}
   -  \DP{(\alpha Div)^{\tau}}{z} 
  \right\}_{i,k(w)}
 + \left( F_{w}^{t} \right)_{i,k(w)} \Delta \tau   
\right\}
\right]_{i,k}.
\Deqlab{uwpi:sabun_ik}
\nonumber \\
\end{eqnarray}
%
但し平均場の量は鉛直方向にしか依存しないので $z$ 方向の添字のみ
付けてある. 


\subsubsection{境界条件}

上下境界を固定壁とする場合, 境界条件は上部下部境界で, 
%
\begin{eqnarray}
&& w(i,0(w)) = 0, \\
&& w(i,km(w)) = 0 
\end{eqnarray}
%
である.

{\bf 下部境界}: \\

下部境界($k(w) = 0(w)$)について考える. この時 \Deqref{uwpi:w_sabun} 式に
添字を付けて書き下すと, 
%
\begin{eqnarray}
 \beta \left(
	\DP{\pi^{\tau + \Delta \tau}}{z}  
       \right)_{i,0(w)}
 &=&
\left( \DP{(\alpha Div)^{\tau}}{z} \right)_{i,0(w)} 
  - (1 - \beta) \left( \DP{\pi^{\tau}}{z} \right)_{i,0(w)} 
  + \left(\Dinv{\bar{c_{p}} \bar{\theta}_{v}} F_{w}^{t}\right)_{i,0(w)} 
\nonumber \\
  &\equiv& E_{i,0(w)}
  \Deqlab{uwpi:w_sabun_kabu}
\end{eqnarray}
%
となる. したがって \Deqref{uwpi:sabun_ik} 式は以下のようになる. 
%
\begin{eqnarray}
&&  \left\{
 - \beta^{2} 
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho}
  \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z^{2}}
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{1(w)}
  \right\}
  \pi^{\tau + \Delta \tau}_{i,2}
+ 
  \left\{
  1 + \beta^{2} 
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho}
  \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z^{2}}
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{1(w)}
  \right\} \pi^{\tau + \Delta \tau}_{i,1}
 \nonumber \\
 &=& 
  \pi^{\tau}_{i,1}
 -(1 - \beta)  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{1}
  \left\{
  \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,1}
  -  
  \left(
   \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}}
  \right)_{1}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,1}
  \nonumber \\
 &&- \beta  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{1}
   \left[
    \DP{}{z}
    \bar{\rho} \bar{\theta}_{v}
   \left\{
   w^{\tau}
  - \bar{c_{p}} \bar{\theta}_{v} \Delta \tau
  \left\{
    (1 - \beta) \DP{\pi^{\tau}}{z} 
   - \DP{(\alpha Div)^{\tau}}{z} 
  \right\}
 + F_{w}^{t} \Delta \tau   
\right\}
\right]_{i,1}
\nonumber \\
 &&- \beta
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z} 
%\left[
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{i,0(w)}
   E_{i,0(w)}.
\Deqlab{uwpi:sabun_kabu}
\end{eqnarray}
%

{\bf 上部境界}: \\

上部境界($k(w) = km(w)$)について考える. この時 \Deqref{uwpi:w_sabun} 式
を添字を付けて書き下すと, 
%
\begin{eqnarray}
 \beta \left(
	\DP{\pi^{\tau + \Delta \tau}}{z}  
      \right)_{i,km(w)}
 &=&
 \left( \DP{(\alpha Div)^{\tau}}{z} \right)_{i,km(w)}
  - (1 - \beta) \left( \DP{\pi^{\tau}}{z} \right)_{i,km(w)}
  + \left(\Dinv{\bar{c_{p}} \bar{\theta}_{v}} F_{w}^{t}\right)_{i,km(w)}
\nonumber \\
 &\equiv& E_{i,km(w)}
  \Deqlab{uwpi:w_sabun_joubu}
\end{eqnarray}
%
となる. したがって \Deqref{uwpi:sabun_ik} 式は以下のようになる. 
%
\begin{eqnarray}
&&  \left\{
   1 + 
  \beta^{2} 
  \left(
   \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{km}
  \Dinv{\Delta z^{2}}
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{km-1(w)}
  \right\}
    \pi^{\tau + \Delta \tau}_{i,km}
  \nonumber \\
&& + \left\{
  - \beta^{2} 
  \left(
   \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{km}
  \Dinv{\Delta z^{2}}
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{km-1(w)}
  \right\}
     \pi^{\tau + \Delta \tau}_{i,km-1}
  \nonumber \\
 &=& 
  \pi^{\tau}_{i,km}
-(1 - \beta)  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{km}
  \left\{
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,km}
-  
  \left(
   \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} 
  \right)_{km}
  \left(
    \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,km}
  \nonumber \\
 &&- \beta  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{km}
   \left[
    \DP{}{z}
    \bar{\rho} \bar{\theta}_{v}
   \left\{
   w^{\tau}
  - \bar{c_{p}} \bar{\theta}_{v} \Delta \tau
  \left\{
    (1 - \beta) \DP{\pi^{\tau}}{z} 
   - \DP{(\alpha Div)^{\tau}}{z} 
  \right\}
 + F_{w}^{t} \Delta \tau   
\right\}
\right]_{i,km}
\nonumber \\
 &&+ \frac{\beta}{\Delta z}
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}  
  \right)_{km}
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{km(w)}
   E_{i,km(w)}.
\Deqlab{uwpi:sabun_joubu}
\end{eqnarray}
%

\subsubsection{圧力方程式の解き方}

\Deqref{uwpi:sabun_ik}, \Deqref{uwpi:sabun_kabu},
\Deqref{uwpi:sabun_joubu} 式を連立すると, 以下のような行列式の形式で書く
ことができる. 
%
 \begin{eqnarray}
  \left(\begin{array}{cccc}
          A_{1} & B_{2}  &        &   0      \\ 
          C_{1} & \ddots & \ddots &          \\
	 & \ddots & \ddots & B_{km}   \\
            0   &        & C_{km-1} & A_{km} \\
        \end{array}
  \right)
  \left(\begin{array}{cccc}
          \pi_{1,1} & \pi_{2,1}  & \cdots & \pi_{im,1} \\
          \pi_{1,2} & \ddots    & \ddots & \vdots        \\
          \vdots   &           & \ddots & \vdots        \\
          \pi_{1, km} & \cdots & \cdots & \pi_{im, km} \\
        \end{array}
  \right)^{\tau + \Delta \tau}
  \nonumber \\
  = 
  \left(\begin{array}{cccc}
   D_{1,1} & D_{2,1}  & \cdots & D_{im,1} \\
   D_{1,2} & \ddots    & \ddots & \vdots        \\
        \vdots   &           & \ddots & \vdots        \\
   D_{1,km} & \cdots & \cdots & D_{im,km} \\
        \end{array}
  \right)^{\tau}
  .
 \end{eqnarray}
%
この連立方程式を解くことで $\pi_{i,k}$ を求める. この連立方程式の係数は以下の
ように書ける. 
%
\begin{eqnarray}
 A_{k} 
  &=& 1 + \beta^{2}
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho}
  \bar{\theta}_{v}^{2}} 
  \right)_{k}
  \Dinv{\Delta z^{2}} 
   \left\{
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{k(w)}
   + 
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{k-1(w)}
   \right\} 
   \nonumber \\   
 && (k = 2, 3, \cdots km-1),
  \nonumber \\
 A_{1} 
  &=& 1 + \beta^{2}
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z^{2}} 
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{1(w)},
   \nonumber \\   
 A_{km} 
  &=& 1 + \beta^{2}
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{km}
  \Dinv{\Delta z^{2}} 
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{km-1(w)},
   \nonumber \\   
 B_{k} &=&
  - \beta^{2}
  \left(
   \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k-1}
  \Dinv{\Delta z^{2}} 
  \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
  \right)_{k-1(w)},
  \nonumber \\
 C_{k} &=&
  - \beta^{2}
  \left(
   \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k+1}
  \Dinv{\Delta z^{2}} 
   \left(
    \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} 
       \right)_{k(w)},
   \nonumber \\
 D_{i,k}
 &=& \pi^{\tau}_{i,k}
   -(1 - \beta)  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k}
  \left\{
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,k}
   -\left(
     \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} 
    \right)_{k}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x} 
  \right)_{i,k}
  - F_{i,k}
\nonumber \\
 && (k = 2, 3,  \cdots km-1),
\nonumber \\
 D_{i,1} &=& \pi^{\tau}_{i,1}
 -(1 - \beta)  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{1}
  \left\{
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,1}
-  
  \left(
   \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}}
  \right)_{1}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,1} - F_{i,k}
  \nonumber \\
 &&- \beta
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z} 
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{i,0(w)}
   E_{i,0(w)},
\nonumber \\
 D_{i,km} &=& 
  \pi^{\tau}_{i,km}
  -(1 - \beta)  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{km}
  \left\{
   \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,km}
-  
  \left(
   \frac{\bar{c}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} 
  \right)_{km}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,km} -F_{i,k}
  \nonumber \\
 &&+ \beta
  \left(
  \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}  
  \right)_{km}
  \Dinv{\Delta z}
   \left(
   \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{km(w)}
   E_{i,km(w)}.
\nonumber 
\end{eqnarray}
%
ただし, 
%
\begin{eqnarray}
 E_{i,k(w)} \equiv 
  \left( \DP{(\alpha Div)^{\tau}}{z} \right)_{i,k(w)}
  - (1 - \beta) \left( \DP{\pi^{\tau}}{z} \right)_{i,k(w)}
  + \left(\Dinv{\bar{c_{p}} \bar{\theta}_{v}} F_{w}^{t}\right)_{i,k(w)}
  \nonumber 
\end{eqnarray}
%
\begin{eqnarray}
F_{i,k} 
 &\equiv&\hspace{2mm} - \beta  
  \left(
   \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2}}
  \right)_{k}
   \left[
    \DP{}{z}
    \left( \bar{\rho} \bar{\theta}_{v} \right)_{i,k(w)} 
   \left\{
   w^{\tau}_{i,k(w)}
\right. \right.
\nonumber \\
&& \hspace{10mm} 
\left. \left.
  - \left( \bar{c_{p}} \bar{\theta}_{v} \right)_{i,k(w)} \Delta \tau
  \left\{
    (1 - \beta)  \DP{\pi^{\tau}}{z}
   -  \DP{(\alpha Div)^{\tau}}{z} 
  \right\}_{i,k(w)}
 + \left( F_{w}^{t} \right)_{i,k(w)} \Delta \tau   
\right\}
\right]_{i,k}.
\nonumber
\end{eqnarray}
%
である. 

\subsection{運動方程式の外力項}

運動方程式の外力項 \Deqref{uwpi:u}, \Deqref{uwpi:w} 式を
離散化する. 
%
 \begin{eqnarray}
  F_{u,i(u),k}^{t} &=&
   - \left[ {\rm Adv}.u \right]_{i(u),k}^{t} +
   \left[{\rm D}_{u} + {\rm Diff}.u \right]_{i(u),k}^{t-\Delta t}, 
  \Deqlab{uwpi:Fu}
  \\
  F_{w,i,k(w)} &=&
   - \left[ {\rm Adv}.w \right]_{i(u),k}^{t}
   + \left[ {\rm Buoy} \right]_{i, k(w)}^{t} 
   + \left[ {\rm D}_{w} + {\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t}.
  \Deqlab{uwpi:Fw}
 \end{eqnarray}
%
ここで, Adv は移流項, D は粘性拡散項, Buoy は浮力項である. 
中心差分でリープフロッグ法を用いるために, 数値粘性項 Diff を追加してある. 
それぞれの項を書き下すと, 
%
 \begin{eqnarray}
  \left[ {\rm Adv}_{u} \right]_{i(u),k}^{t} &=&
   - u_{i(u),k}^{t} \left[\DP{u}{x}\right]_{i(u),k}^{t}
   - w_{i(u),k}^{t} \left[\DP{u}{z}\right]_{i(u),k}^{t}  \\
  \left[ {\rm Adv}_{w} \right]_{i,k(w)}^{t} &=&
   - u_{i,k(w)}\left[\DP{w}{x}\right]_{i,k(w)}^{t}
   - w_{i,k(w)}\left[\DP{w}{z}\right]_{i,k(w)}^{t}
 \end{eqnarray}
%
であり, 浮力項は, 
%
 \begin{eqnarray}
  [{\rm Buoy}]^{t}_{i,k(w)} 
   = g \left[ \frac{\theta}{\overline{\theta}} \right]_{i,k(w)}^{t}
 \end{eqnarray}
%
であり, 粘性拡散項は, 
%
 \begin{eqnarray}
  \left[ {\rm D}_{u} \right]_{i(u),k}^{t - \Delta t} &=& 
  2 \left[ 
      \DP{}{x}\left\{ 
          \left( K_{m} \right)_{i,k} \left( \DP{u}{x} \right)_{i,k}
       \right\} 
    \right]_{i(u),k}^{t - \Delta t}
\nonumber \\
 &&
   +\left[ \DP{}{z}\left\{
       \left( K_{m} \right)_{i(u),k(w)}
       \left( \DP{w}{x} \right)_{i(u),k(w)}
     + \left( K_{m} \right)_{i(u),k(w)}
       \left( \DP{u}{z} \right)_{i(u),k(w)}
     \right\} \right]_{i(u),k}^{t - \Delta t}
\nonumber \\
 &&
   - \frac{2}{3 C_{m}^{2} l^{2}} 
     \left( \DP{ K_{m}^{2} }{x} \right)_{i(u),k}^{t - \Delta t}
   \\
%
%
  \left[ {\rm D}_{w} \right]_{i,k(w)}^{t - \Delta t} &=& 
  2 \left[ 
      \DP{}{z}\left\{ 
          \left( K_{m} \right)_{i,k} \left( \DP{w}{z} \right)_{i,k}
       \right\} 
    \right]_{i,k(w)}^{t - \Delta t}
\nonumber \\
 &&
   +\left[ \DP{}{x}\left\{
       \left( K_{m} \right)_{i(u),k(w)}
       \left( \DP{w}{x} \right)_{i(u),k(w)}
     + \left( K_{m} \right)_{i(u),k(w)}
       \left( \DP{u}{z} \right)_{i(u),k(w)}
     \right\} \right]_{i,k(w)}^{t - \Delta t}
\nonumber \\
 &&
   - \frac{2}{3 C_{m}^{2} l^{2}} 
     \left( \DP{ K_{m}^{2} }{z} \right)_{i,k(w)}^{t - \Delta t}
 \end{eqnarray}
%
である. 数値粘性項は,
%
 \begin{eqnarray}
  \left[ {\rm Diff}.u \right]_{i(u),k}^{t - \Delta t} &=& 
     \nu_{h} \left\{ \DP{}{x} \left(\DP{u}{x}\right)_{i,k} \right\}_{i(u),k}^{t - \Delta t}
   + \nu_{v} \left\{ \DP{}{z} \left(\DP{u}{z}\right)_{i(u),k(w)} \right\}_{i(u),k}^{t - \Delta t}
   \\
%
%
  \left[ {\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t} &=&  
     \nu_{h} \left\{ \DP{}{x} \left(\DP{w}{x}\right)_{i(u),k(w)} \right\}_{i,k(w)}^{t - \Delta t}
   + \nu_{v} \left\{ \DP{}{z} \left(\DP{w}{z}\right)_{i,k} \right\}_{i,k(w)}^{t - \Delta t}
 \end{eqnarray}
%
である. $K_{m}$ は乱流エネルギーの時間発展方程式から計算し(詳細は後述), 
$\nu_{h}, \nu_{v}$ は以下のように定める. 
%
\begin{eqnarray}
 \nu_{h} = \frac{\alpha_{h} \Delta x^{2}}{\Delta t}
  \\
 \nu_{v} = \frac{\alpha_{v} \Delta z^{2}}{\Delta t}
\end{eqnarray} 
%
ここで $\Delta x, \Delta z$ は水平・鉛直方向の格子間隔を意味し, 
$\alpha_{h}, \alpha_{v}$ はそれぞれ, 
%
\begin{eqnarray}
  \alpha_{h}  \le \Dinv{8}, \hspace{3em}
  \alpha_{v}  \le \Dinv{8} 
\Deqlab{nu}
\end{eqnarray} 
%
とする.