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

\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$ を計算する.


運動方程式の各項のうち, 音波に関係しない項を $F_u, F_w$ として
まとめると, 運動方程式と圧力方程式は以下のように書ける. 
%
\begin{eqnarray}
&& \DP{u_{i(u),k}}{t} = - \left[\langle c_p \rangle \bar{\theta}_{v} 
                       \DP{(\Pi - \alpha Div )}{x}\right]_{i(u),k} 
  + [F_{u}]_{i(u),k}^{t},
  \Deqlab{uwpi:u}\\
&& \DP{w_{i,k(w)}}{t} = - \left[\langle c_p \rangle \bar{\theta}_{v} 
                       \DP{(\Pi - \alpha Div )}{z}\right]_{i,k(w)} 
  + [F_{w}]_{i,k(w)}^{t},
  \Deqlab{uwpi:w}\\
&& \DP{\Pi_{i,k}}{t}  
  + \left[\frac{\overline{ \langle c_s \rangle}^{2}}
   {\langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2}}
  \DP{(\bar{\rho} \bar{\theta}_{v} w)}{z}\right]_{i,k}
  = 
  -  \left[ \frac{\overline{ \langle c_s \rangle}^{2}}
   { \langle c_p \rangle \bar{\theta}_{v}}
      \DP{u}{x} \right]_{i,k}.
  \Deqlab{uwpi:pi}
\end{eqnarray}
%
ただし $u, w$ の式には音波減衰項 $\alpha Div$ を加えてある
(Skamarock and Klemp, 1992). 
音波に関連しない項 $F_{u}, F_{w}$ は, 
%
 \begin{eqnarray}
&&  [F_{u}]_{i(u),k}^{t} =
   - \left[{\rm Adv}.{u}\right]_{i(u),k}^{t}
   + \left[{\rm Turb}.{u}\right]_{i(u),k}^{t - \Delta t}
   + \left[{\rm Diff}.u\right]_{i(u),k}^{t - \Delta t}
  \Deqlab{uwpi:Fu}
  \\
&&  \left[F_{w}\right]_{i,k(w)}^{t} =
   - \left[{\rm Adv}.{w}\right]_{i,k(w)}^{t}
   + [{\rm Buoy}]_{i,k(w)}^{t}
   + \left[{\rm Turb}.{w}\right]_{i,k(w)}^{t - \Delta t}.
   + \left[{\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t}.
   \Deqlab{uwpi:Fw}
 \end{eqnarray}
%
であり, 時刻 $t$ で評価することにする.
但し, 中心差分でリープフロッグ法を用いるため, 数値粘性項 Diff を追加して
ある.  

\subsection{音波に関連する項の時間方向の離散化}

\subsubsection{水平方向の運動方程式の離散化}

\Deqref{uwpi:u} を時間方向に離散化すると以下のようになる. 
%
\begin{eqnarray}
 u^{\tau + \Delta \tau}_{i(u),k}
  =  u^{\tau}_{i(u), k}
  -   \left[
  \langle c_p \rangle \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}
%

\subsubsection{鉛直方向の運動方程式と圧力方程式の離散化}

HE-VI 法を用いるので, $w$ と $\Pi$ の式を連立して解く.  $w$ の式におい
て音波減衰項は前進差分, 圧力項は後退差分で離散化する.  $\Pi$ の式にお
いて水平微分項は\Deqref{uwpi:u_sabun}で求めた $u^{\tau +\Delta \tau}$ 
を用いて離散化し, 鉛直微分項は後退差分で離散化する.
%
\begin{eqnarray}
 w^{\tau + \Delta \tau}_{i,k(w)}
  &=&  w^{\tau}_{i,k(w)}
  - \langle c_p \rangle \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\} \nonumber \\
 & &
  + F_{w,i,k(w)}^{t} \Delta \tau. 
 \Deqlab{uwpi:w_sabun}
\end{eqnarray}
%
\begin{eqnarray}
% \frac{\Pi^{\tau + \Delta \tau} - \Pi^{\tau}}{\Delta \tau} 
%  + \frac{\bar{c_s}^{2}}{\langle c_p \rangle \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_s}^{2}}{\langle c_p \rangle \bar{\theta}_{v}} \DP{u^{\tau +
%  \Delta \tau}}{x}, 
%  \nonumber \\
%
\Pi^{\tau + \Delta \tau} 
 + \beta \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
 \bar{\theta}_{v}^{2}} 
   \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau + \Delta \tau})}{z}
&=& 
  \Pi^{\tau}
  - (1 - \beta) \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle
  \bar{\rho} \bar{\theta}_{v}^{2}} 
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
 \nonumber \\
& &
  -  \frac{\bar{c_s}^{2} \Delta \tau}{\langle c_p \rangle
  \bar{\theta}_{v}} \DP{u^{\tau + \Delta \tau}}{x}
 - \overline{u} \Delta \tau \DP{\Pi^{\tau}}{x}
 + H_{i,k}^{\tau} \Delta \tau.
  \nonumber \\ 
 \Deqlab{uwpi:pi_sabun}
\end{eqnarray}
但し
\begin{eqnarray}
 H_{i,k}^{\tau}
  &=& 
   \frac{\bar{c_s}^{2}}{\langle c_p \rangle \bar{\theta}_v}
   \Bigg[
    \Dinv{\bar{\rho}} \sum M_{fall}^{\tau}
    + \Dinv{\bar{\theta}} 
      \Bigg\{ \Dinv{\bar{\Pi}} 
       (Q_{cond}^{\tau} + Q_{rad}^{\tau} + Q_{dis}^{\tau}) 
       + D_{\theta}^{t - \Delta t}
      \Bigg\}
  \nonumber \\
  & &
   \Dinv{\bar{q}_d 
     \langle M \rangle / M_d + \sum \bar{q}_v \langle M \rangle / M_v}
   \Bigg\{
    \frac{\langle M \rangle}{M_d}
     \Bigg( - \frac{\bar{q}_d}{\bar{\rho}} \sum M_{fall}^{\tau} 
       + D_{q_d}^{t - \Delta t } \Bigg) \nonumber \\
  & &
    + \sum \frac{\langle M \rangle}{M_v}
     \Bigg( - \frac{\bar{q}_v}{\bar{\rho}} \sum M_{fall}^{\tau} 
       + M_{src}^{\tau}(q_v) D_{q_v}^{t - \Delta t }
     \Bigg)    
   \Bigg\}
   \Bigg]_{i,k}
 \Deqlab{Hdef}
\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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho} 
  \bar{\theta}_{v}^{2}}
  \DP{}{z}
  \left\{
   \left(\langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2}\right)
   \left(
     \DP{\Pi^{\tau + \Delta \tau}}{z} 
   \right)
\right\}
  \nonumber \\
 &=& 
  \Pi^{\tau}
  -(1 - \beta) 
 \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
 \bar{\theta}_{v}^{2}} 
   \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  -  \frac{\bar{c_s}^{2} \Delta \tau}{\langle c_p \rangle \bar{\theta}_{v}} 
  \DP{u^{\tau + \Delta \tau}}{x}
  \nonumber \\
 && - \beta  \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle
  \bar{\rho} \bar{\theta}_{v}^{2}} 
   \left[
    \DP{}{z}
    \bar{\rho} \bar{\theta}_{v}
   \left\{
   w^{\tau}
  - \langle c_p \rangle \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 \\
& &
 - \overline{u} \Delta \tau \DP{\Pi^{\tau}}{x}
 + H_{i,k}^{\tau} \Delta \tau.
\Deqlab{uwpi:sabun}
\end{eqnarray}
%
%この方程式を解いて $\Pi^{\tau + \Delta \tau}$ を求めるわけだが, 左辺の空
%間微分を陽に書き下すことで, 簡単な行列式を導くことができる. 
\Deqref{uwpi:sabun} 式右辺を空間方向に離散化し, 
格子点位置を表す添字を付けて表すと以下のようになる
(計算の詳細は \Dchapref{appendix-a} 参照). 
%
\begin{eqnarray}
&&  \left\{
   - \beta^{2}
  \left(
  \frac{\bar{c_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}} 
  \right)_{k}
  \Dinv{\Delta z^{2}} 
  \left(
    \langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{k}
  \Dinv{\Delta z^{2}} 
   \left\{
   \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{k(w)}
   + 
   \left(
   \langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{k}
  \Dinv{\Delta z^{2}} 
   \left(
    \langle c_p \rangle \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_s}^{2}\Delta \tau}{\langle c_p \rangle \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_s}^{2} \Delta \tau}{\langle c_p \rangle \bar{\theta}_{v}} 
    \right)_{k}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x} 
  \right)_{i,k}
  \nonumber \\
 &&\hspace{2mm} - \beta  
  \left(
   \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \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 \\
&&
\left. \left.
  - \left( \langle c_p \rangle \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 \\
& &
 - \left( \overline{u} \Delta \tau \DP{\Pi^{\tau}}{x} \right)_{i,k}
 + H_{i,k}^{\tau} \Delta \tau.
\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{\langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z^{2}}
   \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{1(w)}
  \right\}
  \Pi^{\tau + \Delta \tau}_{i,2} \nonumber \\
&&
+ 
  \left\{
  1 + \beta^{2} 
  \left(
  \frac{\bar{c_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z^{2}}
   \left(
   \langle c_p \rangle \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_s}^{2}\Delta \tau}{\langle c_p \rangle \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_s}^{2} \Delta \tau}{\langle c_p \rangle \bar{\theta}_{v}}
  \right)_{1}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,1}
  \nonumber \\
 &&- \beta  
  \left(
   \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{1} \nonumber \\
 &&
   \left[
    \DP{}{z}
    \bar{\rho} \bar{\theta}_{v}
   \left\{
   w^{\tau}
  - \langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}} 
  \right)_{1}
  \Dinv{\Delta z} 
%\left[
   \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{i,0(w)}
   E_{i,0(w)}
   - \left( \overline{u} \DP{\Pi^{\tau}}{x} \right)_{i,1} \Delta \tau
   + H_{i,1}^{\tau} \Delta \tau. \nonumber \\
\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)} \nonumber \\
 & &
  + \left(\Dinv{\langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{km}
  \Dinv{\Delta z^{2}}
   \left(
   \langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{km}
  \Dinv{\Delta z^{2}}
   \left(
   \langle c_p \rangle \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_s}^{2}\Delta \tau}{\langle c_p \rangle \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_s}^{2} \Delta \tau}{\langle c_p \rangle \bar{\theta}_{v}} 
  \right)_{km}
  \left(
    \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,km}
  \nonumber \\
 &&- \beta  
  \left(
   \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{km} 
  \nonumber \\
 &&
   \left[
    \DP{}{z}
    \bar{\rho} \bar{\theta}_{v}
   \left\{
   w^{\tau}
  - \langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}}   
  \right)_{km}
   \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2}
   \right)_{km(w)}
   E_{i,km(w)} 
   - \left( \overline{u} \DP{\Pi^{\tau}}{x} \right)_{i,km} \Delta \tau
   + H_{i,km}^{\tau} \Delta \tau. \nonumber \\
\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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}} 
  \right)_{k}
  \Dinv{\Delta z^{2}} 
   \left\{
   \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{k(w)}
   + 
   \left(
   \langle c_p \rangle \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_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}} 
  \right)_{1}
  \Dinv{\Delta z^{2}} 
   \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{1(w)},
   \nonumber \\   
 A_{km} 
  &=& 1 + \beta^{2}
  \left(
  \frac{\bar{c_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}} 
  \right)_{km}
  \Dinv{\Delta z^{2}} 
   \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2} 
   \right)_{km-1(w)},
   \nonumber \\   
 B_{k} &=&
  - \beta^{2}
  \left(
   \frac{\bar{c_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{k-1}
  \Dinv{\Delta z^{2}} 
  \left(
   \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2}
  \right)_{k-1(w)},
   \nonumber \\   
 && (k = 2, 3, \cdots km),
  \nonumber \\
 C_{k} &=&
  - \beta^{2}
  \left(
   \frac{\bar{c_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{k+1}
  \Dinv{\Delta z^{2}} 
   \left(
    \langle c_p \rangle \bar{\rho} \bar{\theta}_{v}^{2} 
       \right)_{k(w)},
   \nonumber \\   
 && (k = 1, 2, \cdots km-1),
   \nonumber \\
 D_{i,k}
 &=& \Pi^{\tau}_{i,k}
   -(1 - \beta)  
  \left(
   \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{k}
  \left\{
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,k} \nonumber \\
 &&
   -\left(
     \frac{\bar{c_s}^{2} \Delta \tau}{\langle c_p \rangle \bar{\theta}_{v}} 
    \right)_{k}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x} 
  \right)_{i,k}
%  - F_{i,k}
  + F_{i,k}
  \hspace{20mm} (k = 2, 3,  \cdots km-1),
\nonumber \\
 D_{i,1} &=& \Pi^{\tau}_{i,1}
 -(1 - \beta)  
  \left(
   \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{1}
  \left\{
    \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,1} \nonumber \\
&&
-  
  \left(
   \frac{\bar{c_s}^{2} \Delta \tau}{\langle c_p \rangle \bar{\theta}_{v}}
  \right)_{1}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,1} 
%  - F_{i,k}
  + F_{i,1}
  \nonumber \\
 &&- \beta
  \left(
  \frac{\bar{c_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}}
  \right)_{1}
  \Dinv{\Delta z} 
   \left(
   \langle c_p \rangle \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_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{km}
  \left\{
   \DP{(\bar{\rho} \bar{\theta}_{v} w^{\tau})}{z}
  \right\}_{i,km} \nonumber \\
&&
-  
  \left(
   \frac{\bar{c_s}^{2} \Delta \tau}{\langle c_p \rangle \bar{\theta}_{v}} 
  \right)_{km}
  \left(
  \DP{u^{\tau + \Delta \tau}}{x}
  \right)_{i,km} 
%  - F_{i,k}
  + F_{i,km}
  \nonumber \\
 &&+ \beta
  \left(
  \frac{\bar{c_s}^{2}{\Delta \tau}^{2}}{\langle c_p \rangle \bar{\rho}
  \bar{\theta}_{v}^{2}} 
  \right)_{km}
  \Dinv{\Delta z}
   \left(
   \langle c_p \rangle \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{\langle c_p \rangle \bar{\theta}_{v}}
     F_{w}^{t}\right)_{i,k(w)} 
  \nonumber \\
%\end{eqnarray}
%
%\begin{eqnarray}
F_{i,k} 
 &\equiv& - \beta  
  \left(
   \frac{\bar{c_s}^{2}\Delta \tau}{\langle c_p \rangle \bar{\rho}
   \bar{\theta}_{v}^{2}} 
  \right)_{k} \nonumber \\
&&
   \Bigg[
    \DP{}{z}
    \left( \bar{\rho} \bar{\theta}_{v} \right)_{i,k(w)} 
   \Bigg\{
   w^{\tau}_{i,k(w)} \nonumber \\
&&
  - \left( \langle c_p \rangle \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   
\Bigg\}
\Bigg]_{i,k}. \nonumber \\
& &
 - \left( \overline{u} \DP{\Pi^{\tau}}{x} \right)_{i,k} \Delta \tau
 + H_{i,k}^{\tau}
\nonumber \\
%\end{eqnarray}
%\begin{eqnarray}
 H_{i,k}^{\tau}
  &\equiv& 
   \frac{\bar{c_s}^{2}}{\langle c_p \rangle \bar{\theta}_v}
   \Bigg[
    \Dinv{\bar{\rho}} \sum M_{fall}^{\tau}
    + \Dinv{\bar{\theta}} 
      \Bigg\{ \Dinv{\bar{\Pi}} 
       (Q_{cond}^{\tau} + Q_{rad}^{\tau} + Q_{dis}^{\tau}) 
       + D_{\theta}^{t - \Delta t}
      \Bigg\}
  \nonumber \\
  & &
   \Dinv{\bar{q}_d 
     \langle M \rangle / M_d + \sum \bar{q}_v \langle M \rangle / M_v}
   \Bigg\{
    \frac{\langle M \rangle}{M_d}
     \Bigg( - \frac{\bar{q}_d}{\bar{\rho}} \sum M_{fall}^{\tau} 
       + D_{q_d}^{t - \Delta t } \Bigg) \nonumber \\
  & &
    + \sum \frac{\langle M \rangle}{M_v}
     \Bigg( - \frac{\bar{q}_v}{\bar{\rho}} \sum M_{fall}^{\tau} 
       + M_{src}^{\tau}(q_v) D_{q_v}^{t - \Delta t }
     \Bigg)    
   \Bigg\}
   \Bigg]_{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 Turb}.{u} \right]_{i(u),k}^{t-\Delta t}
   + \left[{\rm Diff}.u \right]_{i(u),k}^{t-\Delta t}, 
  \Deqlab{uwpi:Fu2}
  \\
&&  F_{w,i,k(w)}^{t} =
   - \left[ {\rm Adv}.w \right]_{i(u),k}^{t}
   + \left[ {\rm Buoy} \right]_{i, k(w)}^{t} 
   + \left[ {\rm Turb}.{w} \right]_{i,k(w)}^{t - \Delta t}
   + \left[ {\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t}. \nonumber \\
  \Deqlab{uwpi:Fw2}
 \end{eqnarray}
%
ここで, Adv は移流項, Turb は粘性拡散項, 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} \nonumber \\
   &&
   + \overline{u}_{i(u),k} \left[\DP{u}{x}\right]_{i(u),k}^{t} 
   + w_{i(u),k}^{t} \left[\DP{\overline{u}}{z}\right]_{i(u),k},  \\
  \left[ {\rm Adv}.{w} \right]_{i,k(w)}^{t} 
   &=&
    u_{i,k(w)}^{t} \left[\DP{w}{x}\right]_{i,k(w)}^{t}
   + w_{i,k(w)}^{t} \left[\DP{w}{z}\right]_{i,k(w)}^{t} \nonumber \\
   &&
   + \overline{u}_{i,k(w)} \left[\DP{w}{x}\right]_{i,k(w)}^{t}
 \end{eqnarray}
%
であり, 浮力項は, 
%
 \begin{eqnarray}
  [{\rm Buoy}]^{t}_{i,k(w)} 
&=& g \frac{\theta_{i,k(w)}^{t}}{\overline{\theta}_{i,k(w)}}
 \nonumber \\
&&  
 + g 
  \frac{[q_d]_{i,k(w)}^{t} \langle M \rangle / M_d 
        + \sum [q_v]_{i,k(w)}^{t} \langle M \rangle / M_v }
       {[\overline{q_d}]_{i,k(w)} \langle M \rangle / M_d 
        + \sum [\overline{q_v}]_{i,k(w)}^{t} \langle M \rangle / M_v }
%   = g \left[ \frac{\theta}{\overline{\theta}} \right]_{i,k(w)}^{t}
 \end{eqnarray}
%
であり, 粘性拡散項は, 
%
 \begin{eqnarray}
&&  \left[ {\rm Turb}.{u} \right]_{i(u),k}^{t - \Delta t} \nonumber \\
&=& 
  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 Turb}.{w} \right]_{i,k(w)}^{t - \Delta t} \nonumber \\
&=& 
  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} 
%
とする
\footnote{
\Deqref{nu} は CreSS のドキュメントにも現れるが, その導出方法は良く分か
らない. 

}. 
