<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">% arare4 火星主成分凝結大気用離散化ドキュメント
%
% 履歴: 2008/11/16 山下 達也 作成開始
%       2011/09/05 山下 達也 改訂
%       2011/09/09 山下 達也 改訂
%       2011/10/26 山下 達也 改訂
%       2011/12/14 山下 達也 改訂
%       2012/03/27 山下 達也 改訂
%       2012/09/09 山下 達也 改訂


%\documentclass[a4j,12pt]{jreport}
\documentclass[a4j,12pt,twoside]{jreport}

\usepackage{Dennou6}         %電脳スタイルを書式として使う
\usepackage{amsmath}
\usepackage{times}
\usepackage{url}
\usepackage{bm}
%\usepackage[multiple]{footmisc}  %同じ場所に複数脚注が入る場合，コンマを
				 %入れる

\pagestyle{Dmyheadings}
%\Dauthor{山下 達也}
%\Ddate{2009/02/02}
%\Dtitle{火星湿潤対流(仮)}
%\Dfile{}

\setlength{\parskip}{2ex}

\setcounter{chapter}{0}

\def\thesection{\arabic{chapter}.\arabic{section}}
\def\thesubsection{\arabic{chapter}.\arabic{section}.\arabic{subsection}}
%\def\thesubsubsection{\arabic{chapter}.\arabic{section}.\arabic{subsection}.\ararbic{subsubsection}}
%\def\theequation{\arabic{chapter}.\arabic{section}.\arabic{equation}}
\def\theequation{\arabic{chapter}.\arabic{equation}}
%\def\thepage{}
\def\thepage{\arabic{page}}
%\def\thetable{\arabic{chapter}.\arabic{section}.\arabic{table}}
\def\thetable{\arabic{chapter}.\arabic{table}}
\def\thefootnote{\arabic{footnote} )}

\makeatletter

\def\sech{\mathop {\operator@font sech}\nolimits}
\def\cosech{\mathop {\operator@font cosech}\nolimits}
\def\Tr{\mathop {\operator@font Tr}\nolimits}
\def\Im{\mathop {\operator@font Im}\nolimits}
\def\Re{\mathop {\operator@font Re}\nolimits}
\def\Res{\mathop {\operator@font Res}\nolimits}

\makeatother

% eqnarray で改ページを許容する
\allowdisplaybreaks


\begin{document}

\Dtitle[大気主成分凝結対流]{ 
\Huge 主成分凝結大気における 2 次元非静力学モデルの離散化 \\
%\Large Two dimensional moist convection including condensation of major
%atmospheric component \\
%\large --- Basic numerical experiments on CO$_2$ ice cloud convection in
%Martian atmosphere ---
}
\Dauthor[山下 達也]{ 
\Huge
山下 達也 \\
YAMASHITA Tatsuya \\
北海道大学理学院宇宙理学専攻 \\ 
%惑星物理学研究室修士課程 2 年 \\
%\large
%Planetary Physics Laboratory,
%Department of Cosmosciences, \\
%\large
%Faculty of Science, Hokkaido University
}
\Ddate[2012/09/09]{2012 年 9 月 9 日}
\Dfile{}

\maketitle

\tableofcontents

\markboth{離散化}
         {離散化}

\chapter{離散化}

第 1 章では主成分凝結を考慮した 2 次元準圧縮方程式系の離散化について述べ
る. 

\section{離散化の概要}

第 1.1 節では離散化の概要について述べ, 詳細の説明は次節以降で行なうもの
とする. 

本モデルにおける格子点の配置方法として, 水平方向には Arakawa C グリッド
(Arakawa and Lamb, 1977), 鉛直方向には Lorenz グリッド(Lorenz, 1960)を
採用する. 
空間方向の離散化は 2 次精度又は 4 次精度の中心差分を用いて行ない, 時間方
向の離散化はモード別時間分割法を用いて行なう. 
音波・凝結に関する項は短い時間ステップ $\Delta \tau$ で離散化し, それ以
外の項は長い時間ステップ $\Delta t$ で離散化する. 
音波に関連する項の離散化には HE-VI 法を用い, 水平方向の運動方程式は前進
差分, 鉛直方向の運動方程式及び圧力方程式は後退差分(クランク・ニコルソン
法)で評価する. 
その他の項の離散化にはリープフロッグ法を用いる. 

\section{格子の配置}

本モデルでは水平方向, 鉛直方向の格子点配置方法として Arakawa C グリッド,
Lorenz グリッドをそれぞれ採用している. 
Arakawa C グリッドとは水平方向のベクトル量とスカラー量を半格子ずらして配
置する格子点配置方法のことを言う. 
Arakawa C グリッドは重力波を表現するのに適しているとされている(Arakawa
and Lamb, 1977). 
Lorenz グリッドとは鉛直方向のベクトル量とスカラー量を半格子ずらして配置
する格子点配置方法のことを言う. 
スカラー量の格子点を $(i,k)$, ベクトル量の水平成分に対する格子点を $(i(u),
k)$, ベクトル量の鉛直成分に対する格子点を $(i, k(w))$, 格子の角に相当す
る点を $(i(u), k(w))$ のように表すことにすると, 格子点の配置は
\Dfigref{grid}のように表される. 
%
%
%
\begin{figure}
 \begin{center}
%  \includegraphics[width=9.0cm]{ps/grid.ps}
  \Depsf[90mm]{ps/grid.ps}
  \caption{格子点の配置. 杉山他(2006) より引用した. \Dfiglab{grid}
}
 \end{center}
  \end{figure}
%
%
%
\section{空間方向の離散化}

第 1.3 節では空間微分の離散化方法とその為に必要な平均操作について説明し
た上で, 準圧縮方程式系の空間方向の離散化について述べる. 

\subsection{平均操作}

第 1.2 節で述べたようにスカラー量の格子点とベクトル量の格子点は互いに半
格子ずつずれている. 
数値計算を行なう上でベクトル量をスカラー量の格子点で評価したり, 或いはス
カラー量をベクトル量の格子点で評価する必要がある. 
その際, 平均操作を行なうことによって半格子ずれた点での値を評価することと
する. 

以下, 計算に必要な平均操作を示す. 
但し $\phi$, $H$, $V$ はそれぞれスカラー量, ベクトル量の水平成分, ベクト
ル量の鉛直成分を表す. 
また下付き添字は格子点位置を表している. 
\begin{eqnarray}
\phi_{i(u), k} &amp;=&amp; \frac{\phi_{i+1, k} + \phi_{i, k}}{2},
 \Deqlab{RisanA} \\
\phi_{i, k(w)} &amp;=&amp; \frac{\phi_{i, k+1} + \phi_{i, k}}{2},
 \Deqlab{RisanB} \\
\phi_{i(u), k(w)} &amp;=&amp; \frac{\phi_{i, k} + \phi_{i+1, k} + \phi_{i, k+1}
 + \phi_{i+1, k+1}}{4}, \Deqlab{RisanC} \\
H_{i, k} &amp;=&amp; \frac{H_{i(u), k} + H_{i-1(u), k}}{2}, \Deqlab{RisanD} \\
H_{i, k(w)} &amp;=&amp; \frac{H_{i(u), k+1} + H_{i-1(u), k+1} + H_{i(u), k} +
 u_{i-1(u), k}}{4}, \Deqlab{RisanE} \\
H_{i(u), k(w)} &amp;=&amp; \frac{H_{i(u), k+1} + H_{i(u), k}}{2},
 \Deqlab{RisanF} \\
V_{i, k} &amp;=&amp; \frac{V_{i, k(w)} + V_{i, k-1(w)}}{2}, \Deqlab{RisanG} \\
V_{i(u), k} &amp;=&amp; \frac{V_{i+1, k(w)} + V_{i, k(w)} + V_{i+1, k-1(w)},
 V_{i, k-1(w)}}{4}, \Deqlab{RisanH} \\
V_{i(u), k(w)} &amp;=&amp; \frac{V_{i+1, k(w)} + V_{i, k(w)}}{2}. \Deqlab{RisanI} 
\end{eqnarray}

\subsection{空間微分の離散化}

空間微分の離散化について述べる. 
音波に関連する項の空間微分については 2 次精度の中心差分を用い, その他の
項の空間微分については 4 次精度の中心差分を用いる. 

以下に 2 次精度の中心差分を用いた微分操作を示す. 
但し $\psi$ は格子の角に相当する点で評価している変数を表す. 
\begin{eqnarray}
\left( \DP{\phi}{x} \right)_{i(u), k} &amp;=&amp; \frac{\phi_{i+1, k} - \phi_{i,
 k}}{\Delta x}, \Deqlab{RisanJ} \\
\left( \DP{\phi}{z} \right)_{i, k(w)} &amp;=&amp; \frac{\phi_{i, k+1} - \phi_{i,
 k}}{\Delta z}, \Deqlab{RisanK} \\
\left( \DP{H}{x} \right)_{i,k} &amp;=&amp; \frac{H_{i(u), k} - H_{i-1(u),
 k}}{\Delta x}, \Deqlab{RisanL} \\
\left( \DP{H}{z} \right)_{i(u), k(w)} &amp;=&amp; \frac{H_{i(u), k+1} - H_{i(u),
 k}}{\Delta z}, \Deqlab{RisanM} \\
\left( \DP{V}{x} \right)_{i(u), k(w)} &amp;=&amp; \frac{V_{i+1, k(w)} - V_{i,
 k(w)}}{\Delta x}, \Deqlab{RisanN} \\
\left( \DP{V}{z} \right)_{i, k} &amp;=&amp; \frac{V_{i, k(w)} - V_{i,
 k-1(w)}}{\Delta x}, \Deqlab{RisanO} \\
\left( \DP{\psi}{x} \right)_{i, k(w)} &amp;=&amp; \frac{\psi_{i(u), k(w)} -
 \psi_{i-1(u), k(w)}}{\Delta x}, \Deqlab{RisanP} \\
\left( \DP{\psi}{z} \right)_{i(u), k} &amp;=&amp; \frac{\psi_{i(u), k(w)} -
 \psi_{i(u), k-1(w)}}{\Delta z}, \Deqlab{RisanQ} \\
\left( \DP{\phi}{z} \right)_{i, k} &amp;=&amp; \frac{\phi_{i, k(w)} - \phi_{i,
 k-1(w)}}{\Delta z}, \Deqlab{RisanQ2} \\
\left( \DP{ (\phi H)}{x} \right)_{i, k} 
 &amp;=&amp; 
  \frac{ \phi_{i(u), k} H_{i(u), k} - \phi_{i-1(u), k} H_{i-1(u), k} }
  {\Delta x}, \Deqlab{RisanQ3} \\
\left( \DP{ (\phi V)}{z} \right)_{i, k} 
 &amp;=&amp; 
  \frac{ \phi_{i, k(w)} V_{i, k(w)} - \phi_{i, k-1(w)} V_{i, k-1(w)} }
  {\Delta z}. \Deqlab{RisanQ4}
\end{eqnarray}

以下, 4 次精度の中心差分を用いた微分操作を示す. 
\begin{eqnarray}
\left( \DP{\phi}{x} \right)_{i(u), k} 
&amp;=&amp; \frac{9}{8} \frac{\phi_{i+1,k} - \phi_{i,k}}{\Delta x} 
   - \Dinv{24} \frac{\phi_{i+2,k} - \phi_{i-1,k}}{\Delta z}, 
   \Deqlab{RisanR} \\
\left( \DP{\phi}{z} \right)_{i, k(w)} 
&amp;=&amp; \frac{9}{8} \frac{\phi_{i,k+1} - \phi_{i,k}}{\Delta x} - \Dinv{24}
\frac{\phi_{i,k+2} - \phi_{i,k-1}}{\Delta z}, 
   \Deqlab{RisanS} \\
\left( \DP{H}{x} \right)_{i,k} 
&amp;=&amp; \frac{9}{8} \frac{H_{i(u), k} - H_{i-1(u), k}}{\Delta x} - \Dinv{24}
\frac{H_{i+1(u), k} - H_{i-2(u), k}}{\Delta z}, 
   \Deqlab{RisanT} \\
\left( \DP{H}{z} \right)_{i(u), k(w)} 
&amp;=&amp; \frac{9}{8} \frac{H_{i(u), k+1} - H_{i(u), k}}{\Delta x} - \Dinv{24}
\frac{H_{i(u), k+2} - H_{i(u), k-1}}{\Delta z}, 
   \Deqlab{RisanU} \\
\left( \DP{V}{x} \right)_{i(u), k(w)} 
&amp;=&amp; \frac{9}{8} \frac{V_{i+1, k(w)} - V_{i, k(w)}}{\Delta x} - \Dinv{24}
\frac{V_{i+2, k(w)} - V_{i-1, k(w)}}{\Delta z}, 
   \Deqlab{RisanV}
\end{eqnarray}
\begin{eqnarray}
\left( \DP{V}{z} \right)_{i, k} 
&amp;=&amp; \frac{9}{8} \frac{V_{i, k(w)} - V_{i, k-1(w)}}{\Delta x} - \Dinv{24}
\frac{V_{i, k+1(w)} - V_{i, k-2(w)}}{\Delta z}, 
   \Deqlab{RisanW} \\
\left( \DP{\psi}{x} \right)_{i, k(w)} 
&amp;=&amp; \frac{9}{8} \frac{\psi_{i(u), k(w)} - \psi_{i-1(u), k(w)} }{\Delta
x} - \Dinv{24} \frac{\psi_{i+1(u), k(w)} - \psi_{i-2(u), k(w)}}{\Delta
z}, \Deqlab{RisanX} \nonumber \\
\\
\left( \DP{\psi}{z} \right)_{i(u), k} 
&amp;=&amp; \frac{9}{8} \frac{\psi_{i(u), k(w)} - \psi_{i(u), k-1(w)} }{\Delta
x} - \Dinv{24} \frac{\psi_{i(u), k+1(w)} - \psi_{i(u), k-2(w)}}{\Delta
z}, \Deqlab{RisanY} \nonumber \\
\\
\left( \DP{\phi}{z} \right)_{i, k} 
&amp;=&amp;
\frac{9}{8} \frac{\phi_{i, k(w)} - \phi_{i, k-1(w)}}{\Delta x} - \Dinv{24}
\frac{\phi_{i, k+1(w)} - \phi_{i, k-2(w)}}{\Delta z},
   \Deqlab{RisanZ} \\
\left( \DP{ (\phi H)}{x} \right)_{i, k} 
&amp;=&amp;
\frac{9}{8} 
 \frac{ \phi_{i(u), k} H_{i(u), k} 
  - \phi_{i-1(u), k} H_{i-1(u), k}}{\Delta x} 
  \nonumber \\
&amp; &amp;
- \Dinv{24}
 \frac{\phi_{i+1(u), k} H_{i+1(u), k} 
  - \phi_{i-2(u), k} H_{i-2(u), k} }{\Delta z},
   \Deqlab{RisanZ-2} \\
\left( \DP{ (\phi V)}{z} \right)_{i, k} 
&amp;=&amp;
\frac{9}{8} 
 \frac{\phi_{i, k(w)} V_{i, k(w)} 
  - \phi_{i, k-1(w)} V_{i, k-1(w)} }{\Delta x} 
  \nonumber \\
&amp; &amp;
- \Dinv{24}
 \frac{\phi_{i, k+1(w)} V_{i, k+1(w)} 
  - \phi_{i, k-2(w)} V_{i, k-2(w)}}{\Delta z}.
   \Deqlab{RisanZ-3}
\end{eqnarray}

\subsection{準圧縮方程式系の空間方向の離散化}

1.1 節, 1.2 節の結果を用いて主成分凝結を考慮した準圧縮方程式系を空間方向
に離散化すると, 以下のように書ける. 

{\bf 水平方向の運動方程式}
%
\begin{eqnarray}
\DP{u_{i(u), k}}{t} 
 &amp;=&amp; - u_{i(u), k} \left( \DP{u}{x} \right)_{i(u), k}
     - w_{i(u), k} \left( \DP{u}{z} \right)_{i(u), k} \nonumber \\
 &amp; &amp; - c_{pv} \overline{\theta}_{i(u), k} 
     \left( \DP{\Pi}{x} \right)_{i(u), k} 
     + \left[ D_u \right]_{i(u), k}
     + \left[ D_u^{sfc} \right]_{i(u), k}. 
     \Deqlab{RisanAA}
\end{eqnarray}
ここで
\begin{eqnarray}
&amp;&amp;
\left[ D_u \right]_{i(u), k}
\nonumber \\
&amp;=&amp;
  2 \left[
      \DP{}{x}\left\{
          \left( K_{m} \right)_{i,k} \left( \DP{u}{x} \right)_{i,k}
       \right\}
    \right]_{i(u),k}
\nonumber \\
 &amp;&amp;
   + \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}
\nonumber \\
 &amp;&amp;
   - \frac{2}{3 C_{m}^{2} l^{2}}
     \left( \DP{ K_{m}^{2} }{x} \right)_{i(u),k},
   \\
&amp;&amp;
\left[ D_u^{sfc} \right]_{i(u), k}
 =
  \delta_{1k}
  \left[ \DP{}{z} (C_D V_{sfc} u_{z_1}) \right]_{i(u),k}
\end{eqnarray}
である. 


%
{\bf 鉛直方向の運動方程式}
%
\begin{eqnarray}
\DP{w_{i, k(w)}}{t} 
 &amp;=&amp; - u_{i, k(w)} \left( \DP{w}{x} \right)_{i, k(w)}
     - w_{i, k(w)} \left( \DP{w}{z} \right)_{i, k(w)} \nonumber \\
 &amp; &amp; - c_{pv} \overline{\theta}_{i, k(w)} 
     \left( \DP{\Pi}{z} \right)_{i, k(w)} 
     + \left[ D_w \right]_{i, k(w)}
     + g 
     \left( 
      \frac{\theta^{\prime}}{\overline{\theta}}
      - \frac{R}{p_0} 
      \frac{\overline{\theta}}{\overline{\Pi}^{c_{vv} / R}}
      \rho_s
     \right)_{i, k(w)}. 
     \Deqlab{RisanAB}
\end{eqnarray}
ここで
\begin{eqnarray}
&amp;&amp;
\left[ D_w \right]_{i, k(w)}
\nonumber \\
&amp;=&amp;
  2 \left[
      \DP{}{z}\left\{
          \left( K_{m} \right)_{i,k} \left( \DP{w}{z} \right)_{i,k}
       \right\}
    \right]_{i,k(w)}
\nonumber \\
 &amp;&amp;
   +\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)}
\nonumber \\
 &amp;&amp;
   - \frac{2}{3 C_{m}^{2} l^{2}}
     \left( \DP{ K_{m}^{2} }{z} \right)_{i,k(w)}
\end{eqnarray}
である. 

%
{\bf 圧力方程式}
%
\begin{eqnarray}
\DP{\Pi_{i,k}}{t}
 &amp;=&amp; - \frac{\overline{c^2}_{i,k}}{c_{pv} \overline{\rho}_{i,k}
 \overline{\theta}^2_{i,k}} 
 \left[ \left( \DP{}{x} (\overline{\rho} \overline{\theta} u)
        \right)_{i,k} +
        \left( \DP{}{z} (\overline{\rho} \overline{\theta} w)
        \right)_{i,k} \right] \nonumber \\
 &amp; &amp;
 + \frac{\overline{c^2}_{i,k}}{c_{pv} \overline{\rho}_{i,k}
 \overline{\theta}_{i,k}} 
 \left( \frac{L}{c_{pv} \overline{T}_{i,k} } -1  \right) 
  \left[ M_{cond} \right]_{i,k} \nonumber \\
 &amp; &amp;
 + \frac{\overline{c^2}_{i,k}}{c_{pv} \overline{\theta}_{i,k}^2
  \overline{\Pi}_{i,k}}
  \left( [Q_{rad}]_{i,k} + [Q_{dis}]_{i,k} \right). 
  \Deqlab{RisanAD}
\end{eqnarray}
%
{\bf 雲密度の式}
%
\begin{eqnarray}
\DP{{\rho_s}_{i,k}}{t} 
 &amp;=&amp; - \left[ \DP{}{x} \left( \rho_s u \right) \right]_{i,k} 
     - \left[ \DP{}{z} \left( \rho_s w \right) \right]_{i,k} 
     + \left[ M_{cond} \right]_{i,k}
     + \left[ M_{fall} \right]_{i,k}
     \nonumber \\
 &amp; &amp; + \left[ D_{\rho_s} \right]_{i,k}. 
  \Deqlab{RisanAE}
\end{eqnarray}
ここで
\begin{eqnarray}
\left[ M_{cond} \right]_{i,k}
 &amp;=&amp; \frac{ 4 \pi {r_d}_{i,k} \overline{\rho}_{i,k} N^{*} }
   {{R_h}_{i,k}} (S_{i,k} - 1), 
 \Deqlab{TrisanS2} \\
{R_h}_{i,k}
 &amp;=&amp; \frac{L^2}{k R [\theta^2]_{i,k} [\Pi^2]_{i,k}}, 
 \Deqlab{TrisanS3} \\
S_{i,k}
 &amp;=&amp; 
  \frac{{p_0}^{\frac{c_{pv}}{R}} \Pi_{i,k}}
  {\left[ p_{\ast} \right]_{i,k}},
 \Deqlab{TrisanS4} \\
\left[ M_{fall} \right]_{i,k}
 &amp;=&amp; 
  \left[ 
   \DP{}{z} \left( \rho_s V_{term} \right) 
  \right]_{i,k}
% =
%   \DP{}{z} \left( \rho_s V_{term} \right)_{i,k(w)}
  \Deqlab{TrisanV2}, \\
{V_{term}}_{i,k} 
 &amp;=&amp; 
  \left[ C_{sc} \right]_{i,k}
  \frac{2 \left[ r_d^2 \right]_{i,k} g \rho_I}{9 \eta_{i,k} }
  \Deqlab{TrisanV3}, \\
{r_d}_{i,k}
 &amp;=&amp; 
 \left(
  r_{aero}^3 + \frac{3 {\rho_s}_{i,k}}
  {4 \rho_I \pi \overline{\rho}_{i,k} N^{*}}
 \right)^{1/3}
  \Deqlab{TrisanV3-2}, \\
\left[ C_{sc} \right]_{i,k}
 &amp;=&amp; 
  1 + \frac{4}{3} \left[ K_n \right]_{i,k}, 
  \Deqlab{TrisanV4} \\
\left[ K_n \right]_{i,k}
 &amp;=&amp; 
  \frac{\lambda_{i,k}}{{r_d}_{i,k}}, 
  \Deqlab{TrisanV4-2} \\
\lambda_{i,k}
 &amp;=&amp;
  \frac{k_B \overline{\theta}_{i,k}}
  {\sqrt{2} \pi \sigma^2 p_0 \overline{\Pi}^{c_{vv}/R}_{i,k}}, 
  \Deqlab{TrisanV5} \\
\eta_{i,k}
 &amp;=&amp; 
 \eta_{ref}
 \left( 
  \frac{T_{ref} + C_{CO_2}}
       {\overline{\theta}_{i,k} \overline{\Pi}_{i,k} + C_{CO_2}} 
 \right)
 \left( 
  \frac{\overline{\theta}_{i,k} \overline{\Pi}_{i,k} }{T_{ref}} 
 \right)^{3/2}, 
  \Deqlab{TrisanV6} \\
\left[ D_{\rho_s} \right]_{i,k}
 &amp;=&amp;
  \left[
    u_{i(u),k}  \left( \DP{\rho_s}{x} \right)_{i(u),k}
  \right]_{i,k}
   +
  \left[
    w_{i,k(w)} \left( \DP{\rho_s}{z} \right)_{i,k(w)}
  \right]_{i,k}
\end{eqnarray}
である. 

%
{\bf 熱力学の式}
%
\begin{eqnarray}
\DP{\theta_{i,k}}{t}
 &amp;=&amp; - u_{i,k} \left( \DP{\theta}{x} \right)_{i,k}
     - w_{i,k} \left( \DP{\theta}{z} \right)_{i,k} \nonumber \\
 &amp; &amp; + \Dinv{\overline{\Pi}_{i,k}} \frac{L
  [M_{cond}]_{i,k}}{\overline{\rho}_{i,k} c_{pv}}
  + \Dinv{\overline{\Pi}_{i,k}}
   \left( [Q_{rad}]_{i,k} + [Q_{dis}]_{i,k} \right) 
  + \left[ D_{\theta} \right]_{i,k}
  + \left[ D_{\theta}^{sfc} \right]_{i,k}.
 \Deqlab{RisanAC}
\end{eqnarray}
ここで
\begin{eqnarray}
\left[ D_{\theta} \right]_{i,k}
 &amp;=&amp;
  \left[
    u_{i(u),k}  \left( \DP{\theta}{x} \right)_{i(u),k}
  \right]_{i,k}
   +
  \left[
    w_{i,k(w)} \left( \DP{\theta}{z} \right)_{i,k(w)}
  \right]_{i,k}, 
  \\
\left[ D_{\theta}^{sfc} \right]_{i,k}
 &amp;=&amp;
  \delta_{1k} \left\{ 
   \DP{}{z} [C_D V_{sfc} \rho (T_{z_1} - T_{sfc})]
  \right\}_{i,k}
\end{eqnarray}
である. 

%
{\bf 乱流拡散係数の式}
%
\begin{eqnarray}
\DP{{K_m}_{i,k}}{t} 
&amp;=&amp; - u_{i,k} \left( \DP{{K_m}}{x} \right)_{i,k} 
    - w_{i,k} \left( \DP{{K_m}}{z} \right)_{i,k} \nonumber \\
&amp; &amp; - \frac{3 g {C_m}^2 l^2}{2}
    \left[ 
     \Dinv{\overline{\theta}_{i,k}}
     \left( \DP{\theta}{z} \right)_{i,k}
     - \frac{R}{p_0} 
       \left( 
        \frac{\overline{\theta}}{\overline{\Pi}^{c_v/R}} 
       \right)_{i,k}
       \left( \DP{\rho_s}{z} \right)_{i,k}
    \right]
     \nonumber \\
&amp; &amp; + {C_m}^2 l^2 \left[ \left( \DP{u}{x} \right)^2_{i,k} 
       + \left( \DP{w}{z} \right)^2_{i,k} \right] 
    + \frac{{C_m}^2 l^2}{2} \left( \DP{u}{z} + \DP{w}{x} \right)^2_{i,k}
    \nonumber \\
&amp; &amp; - \frac{{K_m}_{i,k}}{3} \left( \DP{u}{x} + \DP{w}{z} \right)_{i,k}
    + \Dinv{2} \left[ 
    \left( \DP[2]{{K_m}^2}{x} \right)_{i,k}
    + \left( \DP[2]{{K_m}^2}{z} \right)_{i,k} \right]
    \nonumber \\
&amp; &amp; + \left( \DP{{K_m}}{x} \right)^2_{i,k} 
    + \left( \DP{{K_m}}{z} \right)^2_{i,k}
    - \frac{C_{\varepsilon}}{2 C_m l^2} {K_m}^2_{i,k}. 
 \Deqlab{RisanAF}
\end{eqnarray}

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

第 1.4 節では準圧縮方程式の時間方向の離散化について述べる. 

\subsection{モード別時間分割法}

一般に安定に計算を進める為には少なくとも CFL 条件を満たしている必要がある. 
例えば簡単な例として 1 次元移流方程式
\begin{equation}
\DP{\phi}{t} + U \DP{\phi}{x} = 0 \Deqlab{TrisanA}
\end{equation}
を考えると, CFL 条件は
\begin{equation}
\frac{U \Delta t}{\Delta x} \leq 1 \Deqlab{TrisanB}
\end{equation}
と表される. 
但し $U$, $\Delta t$ はそれぞれ移流の速さ, 時間ステップである. 
系において様々な速度スケールの現象が生じている場合, 最も速い速度スケール
を持つ現象が CFL 条件に制約を加えることになる. 

準圧縮方程式は音波を解に含んでいる. 
本研究では対流に着目しているので, 音波自体はあまり重要ではない. 
しかし音波の位相速度は対流の速度スケールに比べて 10 倍程度大きい. 
従ってたとえ対流のみに着目しようとしても, 計算を安定に進めるために時間ス
テップを小さくとらなければならなくなり, 計算のコストが高くなってしまう. 
そこで計算の効率化を図るためにモード別時間分割法を採用する. 
モード別時間分割法とは時間ステップを 2 種類用意し, 短い方の時間ステップ
で音波に関連する項を解き, 長い時間ステップで音波に関連しない移流項や拡散
項を解くという方法である. 
短い時間ステップで時間積分を行なっている間は長い時間ステップで評価する項
の値は一定とみなして計算を行なう. 
モード別時間分割法の概念図を \Dfigref{split} に示す. 
%
%
%
\begin{figure}[b]
 \begin{center}
%  \includegraphics[width=7.8cm]{ps/timesplit.ps}
  \Depsf[78mm]{ps/timesplit.ps}
  \caption{モード別時間分割法の概念図. 
  北守(2006) より引用した. \Dfiglab{split}
}
 \end{center}
  \end{figure}
%
%
%


凝結に関連する項は音波にも移流にも直接関連しないので, 解くべき時間ステッ
プは凝結の時間スケールによって決まると考えられる. 
北守(2006)は Odaka et al.(1998) の火星乾燥対流の実験結果をもとに CFL 条
件を満たす長い時間ステップと短い時間ステップの最大値をそれぞれ 5.0 [s],
0.5 [s] と見積もった. 
また北守(2006)は流れの存在しない火星大気での拡散成長についての数値計算を
行ない, 凝結の時間スケールが 1 -- 20 [s] 程度となることを見出した. 
そこで北守(2006)同様, 凝結に関連する項を短い時間ステップで解くこととする. 
本文書では長い時間ステップを $\Delta t$, 短い時間ステップを $\Delta
\tau$ と書くことにする. 

\subsection{音波減衰項}

モード別時間分割法を用いると音波についての CFL 条件を満たしているにもか
かわらず計算不安定を起こすことがある(Skamarock and Klemp, 1992). 
この計算不安定を抑制する為に, 運動方程式 \Deqref{RisanAA}, 
\Deqref{RisanAB} の
$\Pi^{\prime}$ を
\begin{equation}
\Pi^{\prime} - \alpha {\rm Div}
 \Deqlab{SoundDecay1}
\end{equation}
で置き換える. 
但し 
${\rm Div} \equiv (\partial u / \partial x) + (\partial w / \partial z)$
である. 
音波を選択的に減衰させる為には, 係数 $\alpha$ を適切な値に設定する必要が
ある. 
本モデルでは北守(2006)に従い
\begin{equation}
\alpha = 5.0 \times 10^{-7} \times \frac{{\rm min} [(\Delta x)^2, (\Delta
 z)^2]}{\Delta \tau} \Deqlab{SoundDecay2}
\end{equation}
とする
\footnote{
基本場や計算領域の鉛直方向の長さを変更すると, 計算が即破綻してしまうこと
がある. 
その際, 経験的には $\alpha$ の値を適宜変えてみると計算がうまく進む場合が
多い. 
}. 


\subsection{数値粘性項}

移流項の空間微分を中心差分を用いて離散化すると計算不安定が生じることがあ
る. 
この計算不安定を抑制する為に運動方程式, 熱力学の式, 雲密度の式の移流項に
人工的な数値粘性項を加える. 
即ち任意の予報変数 $\phi$ に関する方程式に数値粘性項
\begin{equation}
\nu_H \DP[2]{\phi}{x} + \nu_V \DP[2]{\phi}{z} \Deqlab{NumDisp1}
\end{equation}
を付加する. 
係数 $\nu_H$, $\nu_V$ については
\begin{equation}
\nu_H = \nu_V = 1.0 \times 10^{-4} \times \frac{{\rm min} [(\Delta x)^2,
 (\Delta z)^2]}{\Delta t} \Deqlab{NumDisp2}
\end{equation}
とする. 


\subsection{準圧縮方程式系の時間方向の離散化}

音波及び凝結に関する項は短い時間ステップ $\Delta \tau$ で, それ以外の項
は長い時間ステップ $\Delta t$ で離散化する. 
音波に関連する項の離散化には HE-VI (horizontally explicit - vertically
implicit)法を用いる. 
即ち $u$ の式は前進差分, $w$, $\Pi^{\prime}$ の式は後退差分(クランク・ニ
コルソン法)で離散化する. 
音波に関連しない項の項はリープフロッグ法で離散化する. 


以下, 時間積分により求まる量を $\phi^{t + \Delta t}$, $\phi^{\tau +
\Delta \tau}$, 時間積分によって得られている最新の物理量を $\phi^{t}$,
$\phi^{\tau}$, 最新の物理量の 1 ステップ前の時刻での物理量を $\phi^{t -
\Delta t}$, $\phi^{\tau - \Delta \tau}$ と表すこととする. 


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


式 \Deqref{RisanAA} を離散化すると
\begin{equation}
u_{i(u),k}^{\tau + \Delta \tau}
= u_{i(u),k}^{\tau} 
  - \left( c_{pv} \overline{\theta} \DP{[\Pi^{\prime}]^{\tau}}{x}
    \right)_{i,u(k)} \Delta \tau
  + \left( c_{pv} \overline{\theta} \DP{[{\rm Div}]^{\tau}}{x}
    \right)_{i,u(k)} \Delta \tau
  + {F_u}_{i(u),k}^{t} \Delta \tau 
  \Deqlab{TrisanC}
\end{equation}
となる. 
$F_u$ は音波, 凝結に関連しない項で
\begin{equation}
{F_u}_{i(u),k}^{t}
= - \left( u^t \DP{u^t}{x} \right)_{i(u), k}
  - \left( w^t \DP{u^t}{z} \right)_{i(u), k}
  + \left[ D_u \right]_{i(u),k}^{t - \Delta t} 
  + \left[ D_u^{sfc} \right]_{i(u),k}^{t}
  \Deqlab{TrisanD}
\end{equation}
である. 
乱流拡散項は坪木・榊原(2001)と同様に全て時刻 $t - \Delta t$ での値で評価
する. 


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


本モデルでは HE-VI 法を用いるので, 鉛直方向の運動方程式と圧力方程式を連
立させて解く. 
式 \Deqref{RisanAB} において音波減衰項, 圧力勾配項はそれぞれ前進差分, 後
退差分で離散化する. 
式 \Deqref{RisanAD} においてフラックス項は後退差分で離散化する. 
上下境界には剛体壁が存在するものとし, 運動学的境界条件を適用する. 

式 \Deqref{RisanAB} 及び 式 \Deqref{RisanAD} を離散化すると
\begin{eqnarray}
w_{i,k(w)}^{\tau + \Delta \tau}
&amp;=&amp; w_{i,k(w)}^{\tau} 
    - c_{pv} \overline{\theta}_{i,k(w)}
    \left[ \beta \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} + (1 -
     \beta) \DP{ [\Pi^{\prime}]^{\tau} }{z} \right]_{i,k(w)} \Delta \tau
    \nonumber \\
&amp; &amp; + c_{pv} \overline{\theta}_{i,k(w)}
    \left[ \DP{(\alpha {\rm Div})}{z} \right]_{i,k(w)} \Delta \tau
    \nonumber \\
&amp; &amp; + {F_w}_{i,k(w)}^{t} \Delta \tau, \Deqlab{TrisanE}
\end{eqnarray}
\begin{eqnarray}
\left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} 
&amp;=&amp; \left[ \Pi^{\prime} \right]_{i,k}^{\tau}
    - \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\rho}
       \overline{\theta}^2 } \DP{}{x} (\overline{\rho} \overline{\theta}
       u^{\tau + \Delta \tau})\right]_{i,k} \Delta \tau \nonumber \\
&amp; &amp;  - \beta \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\rho}
       \overline{\theta}^2 } \DP{}{z} (\overline{\rho} \overline{\theta}
       w^{\tau + \Delta \tau})\right]_{i,k} \Delta \tau 
     - (1 - \beta) \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\rho}
       \overline{\theta}^2 } \DP{}{z} (\overline{\rho} \overline{\theta}
       w^{\tau})\right]_{i,k} \Delta \tau  \nonumber \\
&amp; &amp; + \left[ \frac{\overline{c^2}}{c_{pv} \overline{\theta}} 
       \left( \frac{L}{c_{pv} \overline{T}} - 1 \right)
       \frac{M_{cond}^{\tau}}{\overline{\rho}} \right]_{i,k} \Delta \tau
      \nonumber \\
&amp; &amp; + \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\theta}^2
       \overline{\Pi} } ( Q_{rad} + Q_{dis} ) 
      \right]_{i,k} \Delta \tau
      \Deqlab{TrisanF}
\end{eqnarray}
となる. 
但し $\beta$ は後退差分における重み係数を表し, クランク・ニコルソン法の
場合 $\beta = 0.5$ とする. 
また $F_w$ は音波, 凝結に関連しない項で
\begin{eqnarray}
{F_w}_{i,k(w)}^{t}
 &amp;=&amp; 
  - \left( u^t \DP{w^t}{x} \right)_{i, k(w)}
  - \left( w^t \DP{w^t}{z} \right)_{i, k(w)}
  + g 
    \frac{ [\theta^{\prime}]_{i,k(w)}^{t}}{\overline{\theta}_{i,k(w)}} 
  - g 
    \frac{R}{p_0} 
     \frac{\overline{\theta}_{i, k(w)}}
     { \overline{\Pi}^{c_{vv} / R}_{i, k(w)} }
     [\rho_s]_{i, k(w)}^{t}
  \nonumber \\
 &amp; &amp;
  + \left[ D_w \right]_{i,k(w)}^{t - \Delta t}
 \Deqlab{TrisanG}
\end{eqnarray}
である. 
式 \Deqref{TrisanF} に 式 \Deqref{TrisanC}, \Deqref{TrisanE} を代入して
$u^{\tau + \Delta \tau}$, $w^{\tau + \Delta \tau}$ を消去し, 式を整理す
ると以下のような行列で表記出来る
\footnote{式 \Deqref{TrisanH} -- 式 \Deqref{TrisanR} の導出については付
録 A を参照されたい. }. 

\begin{eqnarray}
 \begin{pmatrix}
  A_1 &amp; B_2 &amp; \cdots &amp; 0 &amp; 0 \\
  C_1 &amp; A_2 &amp; \cdots &amp; 0 &amp; 0 \\
  \colon &amp; \colon &amp;  &amp; \colon &amp; \colon \\
  0 &amp; 0 &amp; \cdots &amp; A_{k_m - 1} &amp; B_{k_m} \\
  0 &amp; 0 &amp; \cdots &amp; C_{k_m - 1} &amp; A_{k_m} 
 \end{pmatrix}
 \begin{pmatrix}
  \Pi^{\prime}_{1,1} &amp; \Pi^{\prime}_{2,1} &amp; \cdots &amp; \Pi^{\prime}_{k_m-1,1} &amp;
  \Pi^{\prime}_{k_m,1} \\
  \Pi^{\prime}_{1,2} &amp; \Pi^{\prime}_{2,2} &amp; \cdots &amp; \Pi^{\prime}_{k_m-1,2} &amp;
  \Pi^{\prime}_{k_m,2} \\
  \colon &amp; \colon &amp; &amp; \colon &amp; \colon \\
  \Pi^{\prime}_{1,k_m-1} &amp; \Pi^{\prime}_{2,k_m-1} &amp; \cdots &amp;
  \Pi^{\prime}_{k_m-1,k_m-1} &amp; \Pi^{\prime}_{k_m,k_m-1} \\
  \Pi^{\prime}_{1,k_m} &amp; \Pi^{\prime}_{2,k_m} &amp; \cdots &amp;
  \Pi^{\prime}_{k_m-1,k_m} &amp; \Pi^{\prime}_{k_m,k_m} 
 \end{pmatrix} ^{\tau + \Delta \tau}
 \nonumber
\end{eqnarray}
\begin{eqnarray}
=
 \begin{pmatrix}
  D_{1,1} &amp; D_{2,1} &amp; \cdots &amp; D_{k_m-1,1} &amp; D_{k_m,1} \\
  D_{1,2} &amp; D_{2,2} &amp; \cdots &amp; D_{k_m-1,2} &amp; D_{k_m,2} \\
  \colon &amp; \colon &amp; &amp; \colon &amp; \colon \\
  D_{1,k_m-1} &amp; D_{2,k_m-1} &amp; \cdots &amp; D_{k_m-1,k_m-1} &amp; D_{k_m,k_m-1} \\
  D_{1,k_m} &amp; D_{2,k_m} &amp; \cdots &amp; D_{k_m-1,k_m} &amp; D_{k_m,k_m}
 \end{pmatrix} ^{\tau}
  \Deqlab{TrisanH}
\end{eqnarray}
ここで $k_m$ は鉛直方向の格子点の総数である. 
式 \Deqref{TrisanH} の左辺の係数行列の各成分は以下のように表される. 
\begin{eqnarray}
A_1 &amp;=&amp; 1 + \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_1
        \Dinv{(\Delta z)^2} 
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,1(w)},
	\Deqlab{TrisanI} \\
A_k &amp;=&amp; 1 + \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
           \overline{\rho} \overline{\theta}^2} \right]_k
        \Dinv{(\Delta z)^2} 
        [ ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k(w)} 
         + ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k-1(w)} 
        ], \nonumber \\
    &amp; &amp;  (k = 2, 3, \dots , k_m -1 ) \Deqlab{TrisanJ} \\
A_{k_m} &amp;=&amp; 1 + \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k_m}
        \Dinv{(\Delta z)^2} 
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k_m-1(w)},
	\Deqlab{TrisanK} \\
B_k &amp;=&amp; - \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k-1}
        \Dinv{(\Delta z)^2}
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k-1(w)},
        \Deqlab{TrisanL} \\
C_k &amp;=&amp; - \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k+1}
        \Dinv{(\Delta z)^2}
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k(w)},
        \Deqlab{TrisanM} \\
D_{i,1}
 &amp;=&amp; [ \Pi^{\prime} ]_{i,1}^{\tau}
    - (1 - \beta) \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{1}
    \left[ \DP{(\overline{\rho} \overline{\theta} w^{\tau})}{z}
    \right]_{i,1} 
    - \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\theta} } \right]_{1}
      \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,1} 
        + F_{i,1} \nonumber \\
 &amp; &amp; 
     - \beta \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
        \overline{\rho} \overline{\theta}^2} \right]_{1}
      \Dinv{\Delta z} (c_{pv} \overline{\rho} \overline{\theta}^2)_{i,0(w)}
      E_{i,0(w)}, \Deqlab{TrisanN} 
\end{eqnarray} 
\begin{eqnarray}
D_{i,k}
 &amp;=&amp; [ \Pi^{\prime} ]_{i,k}^{\tau}
    - (1 - \beta) \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k}
    \left[ \DP{(\overline{\rho} \overline{\theta} w^{\tau})}{z}
    \right]_{i,k} 
    - \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\theta} } \right]_{k}
      \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,k} 
    + F_{i,k}, \nonumber \\
 &amp; &amp; (k = 2, 3, \dots , k_m -1 )
  \Deqlab{TrisanO} \\ 
D_{i,k_m}
 &amp;=&amp; [ \Pi^{\prime} ]_{i,k_m}^{\tau}
    - (1 - \beta) \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k_m}
    \left[ \DP{(\overline{\rho} \overline{\theta} w^{\tau})}{z}
    \right]_{i,k_m} 
    - \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\theta} } \right]_{k_m}
      \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,k_m} \nonumber \\
 &amp; &amp;   + F_{i,k_m}
     + \beta \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
        \overline{\rho} \overline{\theta}^2} \right]_{k_m}
      \Dinv{\Delta z} (c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k_m(w)}
      E_{i,k_m(w)}. \Deqlab{TrisanP}
\end{eqnarray}
但し
\begin{eqnarray}
E_{i,k(w)} 
 &amp;=&amp;  
  \left[ \DP{\left( \alpha {\rm Div} \right)^{\tau}}{z} 
  \right]_{i,k(w)}
  - (1 - \beta )
    \left[ \DP{[\Pi^{\prime}]^{\tau}}{z} \right]_{i,k(w)} 
  + \left[ \frac{F^{t}_{w}}{c_{pv} \overline{\theta}} 
    \right]_{i,k(w)}, 
  \Deqlab{TrisanQ} \\
F_{i,k} 
 &amp;=&amp;  - \beta 
  \left( 
   \frac{\overline{c^2} \Delta \tau}
   {c_{pv} \overline{\rho} \overline{\theta}^2} 
  \right)_{k}
 \nonumber \\
 &amp; &amp;
  \Big[
  \DP{}{z} 
   \Big\{ \overline{\rho} \overline{\theta}
    \Big( w^{\tau} - c_{pv} \overline{\theta} \Delta \tau 
    \left\{ 
     (1 - \beta) \DP{[\Pi^{\prime}]^{\tau}}{z} 
     - \DP{\left( \alpha {\rm Div} \right)^{\tau}}{z} 
    \right\}
    + F_{w}^{t} \Delta \tau \Big) \Big\} \Big]_{i,k}
   \nonumber \\
 &amp; &amp; 
   + \left[
      \frac{\overline{c^2} \Delta \tau }
           {c_{pv} \overline{\rho} \overline{\theta}} 
      \left( \frac{L}{c_{pv} \overline{T}} -1 \right) 
     \right]_{k}
     [M_{cond}^{\tau}]_{i,k}
   \nonumber \\
 &amp; &amp; 
   + \left[
     \frac{ \overline{c^2} \Delta \tau }{c_{pv} \overline{\theta}^2
       \overline{\Pi} } \right]_{k}
     \left[ ( Q_{rad}^{\tau} + Q_{dis}^{\tau} ) \right]_{i,k}
  \Deqlab{TrisanR}
\end{eqnarray}
である. 
式 \Deqref{TrisanH} の左辺の係数行列は 3 重対角行列となっているので,
Thomas 法を用いて時刻 $ \tau + \Delta \tau$ での $\Pi^{\prime}$ の値を求
めることができる(Thomas, 1949). 
本モデルでは計算ライブラリ LAPACK を用いて 式 \Deqref{TrisanH} を解いて
いる. 

{\bf 雲密度の式の離散化}

雲密度の時間積分を行なうに当たっては, 
先ず音波・凝結に関連しない項と雲粒落下項を計算し, 
次に負の雲密度の発生を防ぐ処置を行い, 
最後に音波・凝結に関連する項を計算する
\footnote{
負の雲密度の発生を防ぐ処置の具体的な手順については, 付録 B を参照された
い. 
}. 
即ち, 式 \Deqref{RisanAE} を離散化すると以下のように書ける. 
\begin{eqnarray}
{\rho_s}_{i,k}^{*}
 &amp;=&amp; 
  {\rho_s}_{i,k}^{\tau} 
  + \left[ M_{fall} \right]_{i,k}^{\tau} \Delta \tau
  + \left[ F_{\rho_s} \right]_{i,k}^{t} \Delta \tau, 
  \Deqlab{TrisanU-1} \\
{\rho_s}_{i,k}^{**}
 &amp;=&amp; 
  {\rm FillNeg.} [{\rho_s}_{i,k}^{*}], 
  \Deqlab{TrisanU-2} \\
{\rho_s}_{i,k}^{\tau + \Delta \tau}
 &amp;=&amp; 
  {\rho_s}_{i,k}^{**} 
  + \left[ M_{cond} \right]_{i,k}^{**} \Delta \tau.
  \Deqlab{TrisanU-3}
\end{eqnarray}
ここで ${\rm FillNeg.} [{\rho_s}_{i,k}^{*}]$ は負の雲密度の発生を防ぐ為
の処置を表し, ${\rho_s}_{i,k}^{**}$ は処理を行なった後の雲密度の暫定値,
$F_{\rho_s}$ は音波・凝結に関連しない項であり, 
\begin{eqnarray}
\left[ F_{\rho_s} \right]_{i,k}^{t}
 &amp;=&amp; 
  - \left[ \DP{}{x} ({\rho_s}^t u^t) \right]_{i,k}
  - \left[ \DP{}{z} ({\rho_s}^t w^t) \right]_{i,k}
  + {D_{\rho_s}}_{i,k}^{t}. 
  \Deqlab{TrisanV}
\end{eqnarray}
である. 


{\bf 熱力学の式の離散化}


式 \Deqref{RisanAC}を離散化すると
\begin{equation}
{\theta^{\prime}}_{i,k}^{\tau + \Delta \tau} 
= {\theta^{\prime}}_{i,k}^{\tau} + {F_{\theta}}_{i,k}^{t}
\Delta \tau +  
  \Dinv{\overline{\Pi}_{i,k}} 
   \frac{L \left[ M_{cond} \right]^{**}_{i,k}}
   {c_{pv} \overline{\rho}_{i,k}} \Delta \tau 
 \Deqlab{TrisanS}
\end{equation}
となる. 
$F_{\theta}$ は音波・凝結に関連しない項であり
\begin{eqnarray}
{F_{\theta}}_{i,k}^{t}
&amp;=&amp; - u_{i,k}^{t} \left( \DP{{\theta^{\prime}}^t}{x} \right)_{i,k}
    - w_{i,k}^{t} \left( \DP{{\theta^{\prime}}^t}{z} \right)_{i,k}
    - w_{i,k}^{t} \left( \DP{\overline{\theta}^t}{z} \right)_{i,k}
    \nonumber \\
&amp; &amp; + \Dinv{\overline{\Pi}_{i,k}} 
    ( Q_{rad} + Q_{dis})_{i,k}^{t} 
    + [{D_{\theta}}]_{i,k}^{t - \Delta t}
    + [{D_{\theta}^{sfc}}]_{i,k}^{t}
 \Deqlab{TrisanT}
\end{eqnarray}
である. 


{\bf 乱流拡散係数の式の離散化}


式 \Deqref{RisanAF}をリープフロッグ法を用いて離散化すると
\begin{equation}
{K_m}_{i,k}^{t+ \Delta t}
= {K_m}_{i,k}^{t- \Delta t} + 2 \Delta t [F_{K_m}]_{i,k}^{t}
\Deqlab{TrisanW} 
\end{equation}
となる. 
ここで
\begin{eqnarray}
[F_{K_m}]_{i,k}^{t}
&amp;=&amp; 
- \left[ u_{i(u),k} \left( \DP{K_m}{x} \right)_{i(u),k} \right]_{i,k}^t
- \left[ w_{i,k(w)} \left( \DP{K_m}{z} \right)_{i,k(w)} \right]_{i,k}^t
 \nonumber \\
&amp; &amp;
- \frac{3 g {C_m}^2 l^2}{2}
    \left[ 
     \Dinv{\overline{\theta}}
     \left( \DP{\theta}{z} \right)_{i,k(w)}
     - \frac{R}{p_0} 
       \left( 
        \frac{\overline{\theta}}{\overline{\Pi}^{c_v/R}} 
       \right)
       \left( \DP{\rho_s}{z} \right)_{i,k(w)}
    \right]^{t - \Delta t}_{i,k} 
 \nonumber \\
&amp; &amp;
+ {C_m}^2 l^2 \left[ \left( \DP{u}{x} \right)^2 \right]_{i,k}^{t -
\Delta t}
+ {C_m}^2 l^2 \left[ \left( \DP{w}{z} \right)^2 \right]_{i,k}^{t -
\Delta t} \nonumber \\
&amp; &amp;
+ \frac{{C_m}^2 l^2}{2} \left[ 
    \left\{ \left( \DP{u}{z} \right)_{i(u),k(w)}
       \right\}_{i,k}^{t - \Delta t} 
    + \left\{ \left( \DP{w}{x} \right)_{i(u),k(w)}
       \right\}_{i,k}^{t - \Delta t} 
\right]^2 \nonumber \\
&amp; &amp;
- \Dinv{3} [K_m]_{i,k}^{t - \Delta t}
  \left[ \left( \DP{u}{x} \right)_{i,k}^{t - \Delta t} 
        + \left( \DP{w}{z} \right)_{i,k}^{t - \Delta t} 
\right] \nonumber \\
&amp; &amp;
+ \Dinv{2} \left[
 \left\{ \DP{}{x} \left( \DP{K_m^2}{x} \right)_{i(u),k} 
   \right\}_{i,k}^{t - \Delta t}
 + \left\{ \DP{}{z} \left( \DP{K_m^2}{z} \right)_{i,k(w)} 
   \right\}_{i,k}^{t - \Delta t} \right] \nonumber \\
&amp; &amp;
 + \left[ \left\{ \left( \DP{K_m}{x} \right)^2 \right\}_{i(u),k}
\right]_{i,k}^{t - \Delta t}
 + \left[ \left\{ \left( \DP{K_m}{z} \right)^2 \right\}_{i,k(w)}
\right]_{i,k}^{t - \Delta t} \nonumber \\
&amp; &amp;
 - \frac{C_{\varepsilon}}{2 C_m l^2} [K_m^2]_{i,k}^{t - \Delta t}
\Deqlab{TrisanX}
\end{eqnarray}
である. 
坪木・榊原(2001)と同様に移流項を時刻 $t$ で, それ以外の項を時刻 $t -
\Delta t$ で評価した. 


\subsection{Robert, Asselin の時間フィルター}

1.4 節で述べた通り, 長い時間ステップの計算ではリープフロッグ法を用いてい
る. 
リープフロッグ法では時刻 $t - \Delta t$ の値を用いて時刻 $t + \Delta$ の
値を求める為, 隣接する時間ステップ間で物理量の値に大きな食い違いや振動が
生じる恐れがある. 
この問題を回避する為, 長い時間ステップの計算を 1 回行
なう度に Robert(1966), Asselin(1972) が考案した時間フィルターを適用する. 
本モデルでは $u$, $w$, $\theta^{\prime}$, $\Pi^{\prime}$, $\rho_s$, 
$K_m$ の時間積分に対してフィルターをかける. 

時間フィルター適用前の変数を $\phi$, 時間フィルター適用後の変数を
$\overline{\phi}$ とすると
\begin{equation}
\overline{\phi^t} = \phi^t + \mu_a (\overline{\phi^{t - \Delta t}} - 2
 \phi^t + \phi^{t + \Delta t})
\end{equation}
と表される. 
ここで $\mu_a$ はフィルターの強さを表す係数であり, 本モデルでは Klemp
and Wilhelmson(1978), 坪木・榊原(2001) と同様に $\mu_a = 0.1$ とした. 

\clearpage

% \cleardoublepage

\chapter*{参考文献}

\addcontentsline{toc}{chapter}{参考文献}

\markboth{参考文献}
         {参考文献}

Arakawa, A., Lamb, V. R., 1977 :
Computational design of the basic dynamical processes of the UCLA
general circulation model. 
Methods in Computational Physics, $\bf{17}$, Academic Press, 173 -- 265

Asselin, R., 1972 :
Frequency filter for time integrations.
Monthly Weather Review, $\bf{100}$, 487 -- 490

Deardorff, J. W., 1972: 
Numerical investigation of neutral and unstable planetary boundary
layers. 
Journal of the Atmospheric Sciences, $\bf{29}$, 91 -- 115

Deardorff, J. W., 1975:
The development of boundary-layer turbulence models for use in studying
the severe storm environment. 
Proceeding of the SESAME meeting, Boulder, NOAA-ERL, 251 -- 261

Klemp, J. B., Wilhelmson, R. B., 1978: 
The simulation of three-dimensional convective storm dynamics. 
Journal of the Atmospheric Sciences, $\bf{35}$, 1070 -- 1096

Lorenz, E. N., 1960 :
Energy and numerical weather prediction. 
Tellus, $\bf{12}$, 364 -- 373

Odaka, M., Nakajima, K., Takehiro, S., Ishiwatari, M., Hayashi, Y.-Y.,
1998: 
A numerical study of the Martian atmospheric convection with a
two-dimensional anelastic model. 
Earth Planet Space, $\bf{50}$, 431 -- 437

Robert, A., 1966 :
The integration of a low order spectral form of the primitive
meteorological equations. 
Journal of the Meteorological Society of Japan, $\bf{44}$, 237 -- 244

Skamarock, W. C., Klemp, J. B., 1992: 
The stability of time-split numerical methods for the hydrostatic and
the nonhydrostatic elastic equations. 
Monthly Weather Review, $\bf{120}$, 2109 -- 2127

Thomas, L. H., 1949 :
Elliptic problems in linear difference equation over a network. 
Watson Scientific Computing Laboratory Report, Columbia University, New
York, 

北守 太一, 2006: 
大気主成分の凝結を考慮した火星大気非静力学モデルの開発.  
北海道大学大学院地球惑星科学専攻修士論文, 99 pp \\
http://www.gfd-dennou.org/arch/prepri/2006/hokudai/kitamo/paper/pub/main.pdf

杉山 耕一朗, 北守 太一, 小高 正嗣, 2006:
2 次元非静力学モデルの離散化. \\
arare/deepconv 離散化マニュアル, 42 pp \\
http://www.gfd-dennou.org/library/deepconv/arare/arare4\_current/doc/dai2bu/
\\ pub/main.pdf

坪木 和久, 榊原 篤志, 2001: 
CReSS ユーザーズガイド第 2 版. 253 pp \\
http://cf.tokyo.rist.or.jp/CReSS.top.html

\clearpage

%%%%%%%% 付録 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter*{付録 A: 式 \Deqref{TrisanH} -- 式 \Deqref{TrisanR} の導出}

\addcontentsline{toc}{chapter}
{付録 A: 式 \Deqref{TrisanH} -- 式 \Deqref{TrisanR} の導出 }

\markboth{式 \Deqref{TrisanH} -- 式 \Deqref{TrisanR} の導出}
         {式 \Deqref{TrisanH} -- 式 \Deqref{TrisanR} の導出}

\setcounter{equation}{0}
\setcounter{chapter}{1}
\renewcommand{\thechapter}{\Alph{chapter}}
\renewcommand{\theequation}{\thechapter.\arabic{equation}}

式 \Deqref{TrisanE}, \Deqref{TrisanF} から
式 \Deqref{TrisanH} -- 式 \Deqref{TrisanR} を導出する. 
式 \Deqref{TrisanF} に式 \Deqref{TrisanE} を代入して 
$w^{\tau + \Delta \tau}$ を消去すると, 
\begin{eqnarray}
\left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} 
&amp;=&amp; \left[ \Pi^{\prime} \right]_{i,k}^{\tau}
    - \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\theta} } 
      \DP{}{x} ( u^{\tau + \Delta \tau})\right]_{i,k} \Delta \tau 
      \nonumber \\
&amp; &amp;  - \beta \Bigg[ 
       \frac{ \overline{c^2} }{c_{pv} \overline{\rho}
       \overline{\theta}^2 }
       \DP{}{z}
       \Bigg( 
        \overline{\rho} \overline{\theta}
        \Bigg\{
         w^{\tau} 
         - c_{pv} \overline{\theta}
         \left[ 
          \beta \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} 
          + (1 - \beta) \DP{ [\Pi^{\prime}]^{\tau} }{z} \right] 
         \Delta \tau
         \nonumber \\
&amp; &amp;      + c_{pv} \overline{\theta}
         \left[ \DP{(\alpha {\rm Div})}{z} \right] \Delta \tau
         + {F_w}^{t} \Delta \tau \Bigg\}
       \Bigg) 
       \Bigg]_{i,k} \Delta \tau 
      \nonumber \\
&amp; &amp;  - (1 - \beta) \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\rho}
       \overline{\theta}^2 } \DP{}{z} (\overline{\rho} \overline{\theta}
       w^{\tau})\right]_{i,k} \Delta \tau  \nonumber \\
&amp; &amp; + \left[ \frac{\overline{c^2}}{c_{pv} \overline{\theta}} 
       \left( \frac{L}{c_{pv} \overline{T}} - 1 \right)
       \frac{M_{cond}^{\tau}}{\overline{\rho}} \right]_{i,k} \Delta \tau
      \nonumber \\
&amp; &amp; + \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\theta}^2
       \overline{\Pi} } ( Q_{rad} + Q_{dis} ) 
      \right]_{i,k} \Delta \tau
      \Deqlab{TrisanG-Deriv1}
\end{eqnarray}
となる. 
式 \Deqref{TrisanG-Deriv1} において, 未知の量である $\left[ \Pi^{\prime}
\right]^{\tau + \Delta \tau}$ が含まれる項を左辺にまとめると, 
\begin{eqnarray}
&amp; &amp;
\left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \left[ 
    \DP{}{z} 
     \left\{ 
      (c_{pv} \overline{\rho} \overline{\theta}^2)
      \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z}
     \right\} 
   \right]_{i,k}
  \nonumber \\
&amp;=&amp; \left[ \Pi^{\prime} \right]_{i,k}^{\tau}
    - \left[ 
       \frac{ \overline{c^2} }{c_{pv} \overline{\theta} } 
      \right]_{k}
      \left[  
       \DP{u^{\tau + \Delta \tau}}{x} 
      \right]_{i,k} 
      \Delta \tau 
      \nonumber \\
&amp; &amp; - \beta 
     \left[
      \frac{\overline{c^2} \Delta \tau}
      {c_{pv} \overline{\rho} \overline{\theta}^2} 
     \right]_{k}
 \nonumber \\
 &amp; &amp;
  \Big[
  \DP{}{z} 
   \Big\{ \overline{\rho} \overline{\theta}
    \Big( w^{\tau} - c_{pv} \overline{\theta} \Delta \tau 
    \left\{ 
     (1 - \beta) \DP{[\Pi^{\prime}]^{\tau}}{z} 
     - \DP{\left( \alpha {\rm Div} \right)^{\tau}}{z} 
    \right\}
    + F_{w}^{t} \Delta \tau \Big) \Big\} \Big]_{i,k}
   \nonumber \\
 &amp; &amp;  - (1 - \beta) \left[ \frac{ \overline{c^2} }{c_{pv} \overline{\rho}
       \overline{\theta}^2 } \DP{}{z} (\overline{\rho} \overline{\theta}
       w^{\tau})\right]_{i,k} \Delta \tau  
   \nonumber \\
 &amp; &amp; 
   + \left[
      \frac{\overline{c^2} \Delta \tau }
           {c_{pv} \overline{\rho} \overline{\theta}} 
      \left( \frac{L}{c_{pv} \overline{T}} -1 \right) 
     \right]_{k}
     [M_{cond}^{\tau}]_{i,k}
   + \left[
     \frac{ \overline{c^2} \Delta \tau }{c_{pv} \overline{\theta}^2
       \overline{\Pi} } \right]_{k}
     \left[ ( Q_{rad}^{\tau} + Q_{dis}^{\tau} ) \right]_{i,k}
      \Deqlab{TrisanG-Deriv2}
\end{eqnarray}
となる. 
ここで
\begin{eqnarray}
F_{i,k}
 &amp;=&amp;
    - \beta 
     \left[
      \frac{\overline{c^2} \Delta \tau}
      {c_{pv} \overline{\rho} \overline{\theta}^2} 
     \right]_{k}
 \nonumber \\
 &amp; &amp;
  \Big[
  \DP{}{z} 
   \Big\{ \overline{\rho} \overline{\theta}
    \Big( w^{\tau} - c_{pv} \overline{\theta} \Delta \tau 
    \left\{ 
     (1 - \beta) \DP{[\Pi^{\prime}]^{\tau}}{z} 
     - \DP{\left( \alpha {\rm Div} \right)^{\tau}}{z} 
    \right\}
    + F_{w}^{t} \Delta \tau \Big) \Big\} \Big]_{i,k}
   \nonumber \\
 &amp; &amp; 
   + \left[
      \frac{\overline{c^2} \Delta \tau }
           {c_{pv} \overline{\rho} \overline{\theta}} 
      \left( \frac{L}{c_{pv} \overline{T}} -1 \right) 
     \right]_{k}
     [M_{cond}^{\tau}]_{i,k}
   + \left[
     \frac{ \overline{c^2} \Delta \tau }{c_{pv} \overline{\theta}^2
       \overline{\Pi} } \right]_{k}
     \left[ ( Q_{rad}^{\tau} + Q_{dis}^{\tau} ) \right]_{i,k}
      \Deqlab{TrisanG-Deriv3}
\end{eqnarray}
と置くと, 
\begin{eqnarray}
&amp; &amp;
\left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \left[ 
    \DP{}{z} 
     \left\{ 
      (c_{pv} \overline{\rho} \overline{\theta}^2)
      \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z}
     \right\} 
   \right]_{i,k}
  \nonumber \\
&amp;=&amp; \left[ \Pi^{\prime} \right]_{i,k}^{\tau}
    - (1 - \beta) 
      \left[ \frac{ \overline{c^2} \Delta \tau }
                  {c_{pv} \overline{\rho}
             \overline{\theta}^2 } \right]_k
      \left[ \DP{}{z} (\overline{\rho} \overline{\theta}
             w^{\tau}) \right]_{i,k} 
    - \left[ 
       \frac{ \overline{c^2} \Delta \tau }{c_{pv} \overline{\theta} } 
      \right]_{k}
      \left[  
       \DP{u^{\tau + \Delta \tau}}{x} 
      \right]_{i,k} 
    + F_{i,k}
  \nonumber \\
      \Deqlab{TrisanG-Deriv4}
\end{eqnarray}
となる. 

以下, $k=1$, $k=k_m$, それ以外の 3 つの場合に分けて, 式
\Deqref{TrisanG-Deriv4} を鉛直方向に差分化した式を書き下す. 
鉛直微分は 2 次の中心差分で評価する. 
先ず $k=1$ の場合について考える. 
下部境界条件 $w_{i,1(w)} = 0$ を適用すると, 式 \Deqref{TrisanE} より
\begin{eqnarray}
0
&amp;=&amp; 
    - c_{pv} \overline{\theta}_{i,0(w)}
    \left[ 
     \beta \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} 
     + (1 - \beta) \DP{ [\Pi^{\prime}]^{\tau} }{z} \right]_{i,0(w)} 
    \nonumber \\
&amp; &amp; + c_{pv} \overline{\theta}_{i,0(w)}
    \left[ \DP{(\alpha {\rm Div})}{z} \right]_{i,0(w)}
    + {F_w}_{i,0(w)}^{t}, 
    \nonumber \\
\beta 
 \left[ \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} \right]_{i,0(w)} 
&amp;=&amp; 
  \left[ \DP{\left( \alpha {\rm Div} \right)^{\tau}}{z} 
  \right]_{i,0(w)}
  - (1 - \beta )
    \left[ \DP{[\Pi^{\prime}]^{\tau}}{z} \right]_{i,0(w)} 
  + \left[ \frac{F^{t}_{w}}{c_{pv} \overline{\theta}} 
    \right]_{i,0(w)}
    \Deqlab{TrisanG-Deriv5}
\end{eqnarray}
が成り立つ. 
ここで
\begin{eqnarray}
E_{i,k(w)} 
 &amp;=&amp;  
  \left[ \DP{\left( \alpha {\rm Div} \right)^{\tau}}{z} 
  \right]_{i,k(w)}
  - (1 - \beta )
    \left[ \DP{[\Pi^{\prime}]^{\tau}}{z} \right]_{i,k(w)} 
  + \left[ \frac{F^{t}_{w}}{c_{pv} \overline{\theta}} 
    \right]_{i,k(w)} 
    \Deqlab{TrisanG-Deriv6}
\end{eqnarray}
と置くと, 式 \Deqref{TrisanG-Deriv5} は以下のように書き換えられる. 
\begin{eqnarray}
\beta 
 \left[ \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} \right]_{i,0(w)} 
 = E_{i,0(w)}.
    \Deqlab{TrisanG-Deriv7}
\end{eqnarray}
$k=1$ の場合, 式 \Deqref{TrisanG-Deriv4} の左辺は
式 \Deqref{TrisanG-Deriv7} を用いると
\begin{eqnarray}
&amp; &amp;
\left[ \Pi^{\prime} \right]_{i,1}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{1}
   \left[ 
    \DP{}{z} 
     \left\{ 
      (c_{pv} \overline{\rho} \overline{\theta}^2)
      \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z}
     \right\} 
   \right]_{i,1}
  \nonumber \\
&amp;=&amp;
\left[ \Pi^{\prime} \right]_{i,1}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{1}
  \nonumber \\
&amp; &amp;
   \Dinv{\Delta z}
   \left\{ 
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{1(w)}
     \left[ \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} \right]_{i,1(w)}
     - (c_{pv} \overline{\rho} \overline{\theta}^2)_{0(w)}
       \left[ \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} \right]_{i,0(w)}
   \right\}
  \nonumber \\
&amp;=&amp;
\left[ \Pi^{\prime} \right]_{i,1}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{1}
  \nonumber \\
&amp; &amp;
   \Dinv{\Delta z}
   \left\{ 
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{1(w)}
     \Dinv{\Delta z}
     \left( 
      \left[ \Pi^{\prime} \right]_{i,2}^{\tau + \Delta \tau}
      - \left[ \Pi^{\prime} \right]_{i,1}^{\tau + \Delta \tau}
     \right)
     - (c_{pv} \overline{\rho} \overline{\theta}^2)_{0(w)}
       \Dinv{\beta} E_{i,0(w)}
   \right\}
  \nonumber \\
&amp;=&amp;
- \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{1}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{1(w)}
     \Dinv{(\Delta z)^2}
  \left[ \Pi^{\prime} \right]_{i,2}^{\tau + \Delta \tau} 
  \nonumber \\
&amp; &amp;
 + \left\{ 
    1 +
    \beta^2 
     \left[ 
      \frac{\overline{c^2} \Delta \tau}
      {c_{pv} \overline{\rho} \overline{\theta}^2} 
     \right]_{1}
     \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{1(w)}
   \right\}
     \left[ \Pi^{\prime} \right]_{i,1}^{\tau + \Delta \tau}
  \nonumber \\
&amp; &amp;
 + \beta
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{1}
   \Dinv{\Delta z}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{0(w)}
       E_{i,0(w)}
      \Deqlab{TrisanG-Deriv8}
\end{eqnarray}
となるので, 
\begin{eqnarray}
&amp; &amp;
- \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{1}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{1(w)}
     \Dinv{(\Delta z)^2}
  \left[ \Pi^{\prime} \right]_{i,2}^{\tau + \Delta \tau} 
  \nonumber \\
&amp; &amp;
 + \left\{ 
    1 +
    \beta^2 
     \left[ 
      \frac{\overline{c^2} \Delta \tau}
      {c_{pv} \overline{\rho} \overline{\theta}^2} 
     \right]_{1}
     \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{1(w)}
   \right\}
     \left[ \Pi^{\prime} \right]_{i,1}^{\tau + \Delta \tau}
  \nonumber \\
&amp;=&amp; \left[ \Pi^{\prime} \right]_{i,1}^{\tau}
    - (1 - \beta) 
      \left[ \frac{ \overline{c^2} \Delta \tau }
                  {c_{pv} \overline{\rho}
             \overline{\theta}^2 } \right]_1
      \left[ \DP{}{z} (\overline{\rho} \overline{\theta}
             w^{\tau}) \right]_{i,1} 
    - \left[ 
       \frac{ \overline{c^2} \Delta \tau }{c_{pv} \overline{\theta} } 
      \right]_{1}
      \left[  
       \DP{u^{\tau + \Delta \tau}}{x} 
      \right]_{i,1} 
    + F_{i,1}
  \nonumber \\
&amp; &amp;
 + \beta
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{1}
   \Dinv{\Delta z}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{0(w)}
       E_{i,0(w)}
      \Deqlab{TrisanG-Deriv9}
\end{eqnarray}
となる. 

次に $k = k_m$ の場合について考える. 
下部境界条件 $w_{i,k_m(w)} = 0$ を適用すると, 式 \Deqref{TrisanE} より
\begin{eqnarray}
0
&amp;=&amp; 
    - c_{pv} \overline{\theta}_{i,k_m(w)}
    \left[ \beta \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} + (1 -
     \beta) \DP{ [\Pi^{\prime}]^{\tau} }{z} \right]_{i,k_m(w)}
    \nonumber \\
&amp; &amp; + c_{pv} \overline{\theta}_{i,k_m(w)}
    \left[ \DP{(\alpha {\rm Div})}{z} \right]_{i,k_m(w)}
    \nonumber \\
&amp; &amp; + {F_w}_{i,k_m(w)}^{t}, 
    \nonumber \\
E_{i,k_m(w)}
&amp;=&amp;
\beta 
 \left[ \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} \right]_{i,k_m(w)} 
    \nonumber \\
&amp;=&amp; 
  \left[ \DP{\left( \alpha {\rm Div} \right)^{\tau}}{z} 
  \right]_{i,k_m(w)}
  - (1 - \beta )
    \left[ \DP{[\Pi^{\prime}]^{\tau}}{z} \right]_{i,k_m(w)} 
  + \left[ \frac{F^{t}_{w}}{c_{pv} \overline{\theta}} 
    \right]_{i,k_m(w)}
      \Deqlab{TrisanG-Deriv10}
\end{eqnarray}
となる. 
$k = k_m$ の場合, 
式 \Deqref{TrisanG-Deriv4} の左辺は
式 \Deqref{TrisanG-Deriv10} を用いると
\begin{eqnarray}
&amp; &amp;
\left[ \Pi^{\prime} \right]_{i,k_m}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
   \left[ 
    \DP{}{z} 
     \left\{ 
      (c_{pv} \overline{\rho} \overline{\theta}^2)
      \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z}
     \right\} 
   \right]_{i,k_m}
  \nonumber \\
&amp;=&amp;
\left[ \Pi^{\prime} \right]_{i,k_m}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
  \nonumber \\
&amp; &amp;
   \Dinv{\Delta z}
   \left\{ 
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m(w)}
     \left[ 
      \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} 
     \right]_{i,k_m(w)}
     - (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m - 1(w)}
       \left[ 
        \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} 
       \right]_{i,k_m - 1(w)}
   \right\}
  \nonumber \\
&amp;=&amp;
\left[ \Pi^{\prime} \right]_{i,k_m}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
  \nonumber \\
&amp; &amp;
   \Dinv{\Delta z}
   \left\{ 
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m(w)}
     \Dinv{\beta} F_{i,k_m(w)}
     - (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m - 1(w)}
     \Dinv{\Delta z}
     \left[ 
     [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k_m}
       - [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k_m - 1}
     \right]
   \right\}
  \nonumber \\
&amp;=&amp;
  \left\{ 
   1 + 
   \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m - 1(w)}
  \right\}
     [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k_m}
  \nonumber \\
&amp; &amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m - 1(w)}
       [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k_m - 1}
  \nonumber \\
&amp; &amp;
 - \beta
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m(w)}
     F_{i,k_m(w)}
      \Deqlab{TrisanG-Deriv11}
\end{eqnarray}
となるので, 
\begin{eqnarray}
&amp; &amp;
  \left\{ 
   1 + 
   \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m - 1(w)}
  \right\}
     [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k_m}
  \nonumber \\
&amp; &amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k_m}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m - 1(w)}
       [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k_m - 1}
  \nonumber \\
&amp;=&amp; \left[ \Pi^{\prime} \right]_{i,k_m}^{\tau}
    - (1 - \beta) 
      \left[ \frac{ \overline{c^2} \Delta \tau }
                  {c_{pv} \overline{\rho}
             \overline{\theta}^2 } \right]_{k_m}
      \left[ \DP{}{z} (\overline{\rho} \overline{\theta}
             w^{\tau}) \right]_{i,k_m} 
    - \left[ 
       \frac{ \overline{c^2} \Delta \tau }{c_{pv} \overline{\theta} } 
      \right]_{k_m}
      \left[  
       \DP{u^{\tau + \Delta \tau}}{x} 
      \right]_{i,k_m} 
  \nonumber \\
&amp; &amp;
    + F_{i,k}
    + \beta
      \left[ 
        \frac{\overline{c^2} \Delta \tau}
        {c_{pv} \overline{\rho} \overline{\theta}^2} 
      \right]_{k_m}
        (c_{pv} \overline{\rho} \overline{\theta}^2)_{k_m(w)}
        F_{i,k_m(w)}
      \Deqlab{TrisanG-Deriv12}
\end{eqnarray}
となる. 

更に $k = 1, k_m$ 以外の場合について考える. 
式 \Deqref{TrisanG-Deriv4} の左辺は
\begin{eqnarray}
&amp; &amp;
\left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} 
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \left[ 
    \DP{}{z} 
     \left\{ 
      (c_{pv} \overline{\rho} \overline{\theta}^2)
      \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z}
     \right\} 
   \right]_{i,k}
  \nonumber \\
&amp;=&amp;
\left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} 
  \nonumber \\
&amp; &amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{\Delta z}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k(w)}
     \left[ 
      \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} 
     \right]_{i,k(w)}
  \nonumber \\
&amp; &amp;
 + \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{\Delta z}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k - 1(w)}
       \left[ 
        \DP{[\Pi^{\prime}]^{\tau + \Delta \tau}}{z} 
       \right]_{i,k - 1(w)}
  \nonumber \\
&amp;=&amp;
\left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} 
  \nonumber \\
&amp; &amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k(w)}
     \left[ 
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k+1}
      - [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k}
     \right]
  \nonumber \\
&amp; &amp;
 + \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k - 1(w)}
     \left[ 
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k}
      - [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k-1}
     \right]
  \nonumber \\
&amp;=&amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k(w)}
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k+1}
  \nonumber \\
&amp; &amp;
 + \left\{
   1
   + \beta^2 
     \left[ 
       \frac{\overline{c^2} \Delta \tau}
       {c_{pv} \overline{\rho} \overline{\theta}^2} 
     \right]_{k}
     \Dinv{(\Delta z)^2}
      \left[ 
       (c_{pv} \overline{\rho} \overline{\theta}^2)_{k(w)}
       + (c_{pv} \overline{\rho} \overline{\theta}^2)_{k - 1(w)}
      \right]
   \right\}
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k}
  \nonumber \\
&amp; &amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k - 1(w)}
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k-1}
      \Deqlab{TrisanG-Deriv13}
\end{eqnarray}
となるので, 
\begin{eqnarray}
&amp; &amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k(w)}
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k+1}
  \nonumber \\
&amp; &amp;
 + \left\{
   1
   + \beta^2 
     \left[ 
       \frac{\overline{c^2} \Delta \tau}
       {c_{pv} \overline{\rho} \overline{\theta}^2} 
     \right]_{k}
     \Dinv{(\Delta z)^2}
      \left[ 
       (c_{pv} \overline{\rho} \overline{\theta}^2)_{k(w)}
       + (c_{pv} \overline{\rho} \overline{\theta}^2)_{k - 1(w)}
      \right]
   \right\}
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k}
  \nonumber \\
&amp; &amp;
 - \beta^2 
   \left[ 
     \frac{\overline{c^2} \Delta \tau}
     {c_{pv} \overline{\rho} \overline{\theta}^2} 
   \right]_{k}
   \Dinv{(\Delta z)^2}
     (c_{pv} \overline{\rho} \overline{\theta}^2)_{k - 1(w)}
      [\Pi^{\prime}]^{\tau + \Delta \tau}_{i,k-1}
  \nonumber \\
&amp;=&amp;
\left[ \Pi^{\prime} \right]_{i,k}^{\tau}
    - (1 - \beta) 
      \left[ \frac{ \overline{c^2} \Delta \tau }
                  {c_{pv} \overline{\rho}
             \overline{\theta}^2 } \right]_k
      \left[ \DP{}{z} (\overline{\rho} \overline{\theta}
             w^{\tau}) \right]_{i,k} 
    - \left[ 
       \frac{ \overline{c^2} \Delta \tau }{c_{pv} \overline{\theta} } 
      \right]_{k}
      \left[  
       \DP{u^{\tau + \Delta \tau}}{x} 
      \right]_{i,k} 
    + F_{i,k}
  \nonumber \\
      \Deqlab{TrisanG-Deriv14}
\end{eqnarray}
となる. 
式 \Deqref{TrisanG-Deriv9}, \Deqref{TrisanG-Deriv12},
\Deqref{TrisanG-Deriv14} を行列形式で表記すると, 
\begin{eqnarray}
 \begin{pmatrix}
  A_1 &amp; B_2 &amp; \cdots &amp; 0 &amp; 0 \\
  C_1 &amp; A_2 &amp; \cdots &amp; 0 &amp; 0 \\
  \colon &amp; \colon &amp;  &amp; \colon &amp; \colon \\
  0 &amp; 0 &amp; \cdots &amp; A_{k_m - 1} &amp; B_{k_m} \\
  0 &amp; 0 &amp; \cdots &amp; C_{k_m - 1} &amp; A_{k_m} 
 \end{pmatrix}
 \begin{pmatrix}
  \Pi^{\prime}_{1,1} &amp; \Pi^{\prime}_{2,1} &amp; \cdots &amp; \Pi^{\prime}_{k_m-1,1} &amp;
  \Pi^{\prime}_{k_m,1} \\
  \Pi^{\prime}_{1,2} &amp; \Pi^{\prime}_{2,2} &amp; \cdots &amp; \Pi^{\prime}_{k_m-1,2} &amp;
  \Pi^{\prime}_{k_m,2} \\
  \colon &amp; \colon &amp; &amp; \colon &amp; \colon \\
  \Pi^{\prime}_{1,k_m-1} &amp; \Pi^{\prime}_{2,k_m-1} &amp; \cdots &amp;
  \Pi^{\prime}_{k_m-1,k_m-1} &amp; \Pi^{\prime}_{k_m,k_m-1} \\
  \Pi^{\prime}_{1,k_m} &amp; \Pi^{\prime}_{2,k_m} &amp; \cdots &amp;
  \Pi^{\prime}_{k_m-1,k_m} &amp; \Pi^{\prime}_{k_m,k_m} 
 \end{pmatrix} ^{\tau + \Delta \tau}
 \nonumber
\end{eqnarray}
\begin{eqnarray}
=
 \begin{pmatrix}
  D_{1,1} &amp; D_{2,1} &amp; \cdots &amp; D_{k_m-1,1} &amp; D_{k_m,1} \\
  D_{1,2} &amp; D_{2,2} &amp; \cdots &amp; D_{k_m-1,2} &amp; D_{k_m,2} \\
  \colon &amp; \colon &amp; &amp; \colon &amp; \colon \\
  D_{1,k_m-1} &amp; D_{2,k_m-1} &amp; \cdots &amp; D_{k_m-1,k_m-1} &amp; D_{k_m,k_m-1} \\
  D_{1,k_m} &amp; D_{2,k_m} &amp; \cdots &amp; D_{k_m-1,k_m} &amp; D_{k_m,k_m}
 \end{pmatrix} ^{\tau},
  \Deqlab{TrisanG-Deriv15}
\end{eqnarray}
\begin{eqnarray}
A_1 &amp;=&amp; 1 + \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_1
        \Dinv{(\Delta z)^2} 
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,1(w)},
  \Deqlab{TrisanG-Deriv16} \\
A_k &amp;=&amp; 1 + \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
           \overline{\rho} \overline{\theta}^2} \right]_k
        \Dinv{(\Delta z)^2} 
        [ ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k(w)} 
         + ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k-1(w)} 
        ], \nonumber \\
    &amp; &amp;  (k = 2, 3, \dots , k_m -1 ) 
  \Deqlab{TrisanG-Deriv17} \\
A_{k_m} &amp;=&amp; 1 + \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k_m}
        \Dinv{(\Delta z)^2} 
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k_m-1(w)},
  \Deqlab{TrisanG-Deriv18} \\
B_k &amp;=&amp; - \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k-1}
        \Dinv{(\Delta z)^2}
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k-1(w)},
  \Deqlab{TrisanG-Deriv19} \\
C_k &amp;=&amp; - \beta^2 \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k+1}
        \Dinv{(\Delta z)^2}
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k(w)},
  \Deqlab{TrisanG-Deriv20} \\
D_{i,1}
 &amp;=&amp; [ \Pi^{\prime} ]_{i,1}^{\tau}
    - (1 - \beta) \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{1}
    \left[ \DP{(\overline{\rho} \overline{\theta} w^{\tau})}{z}
    \right]_{i,1} 
    - \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\theta} } \right]_{1}
      \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,1} 
        + F_{i,1} \nonumber \\
 &amp; &amp; 
     - \beta \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
        \overline{\rho} \overline{\theta}^2} \right]_{1}
      \Dinv{\Delta z} (c_{pv} \overline{\rho} \overline{\theta}^2)_{i,0(w)}
      E_{i,0(w)}, 
  \Deqlab{TrisanG-Deriv21}
\end{eqnarray} 
\begin{eqnarray}
D_{i,k}
 &amp;=&amp; [ \Pi^{\prime} ]_{i,k}^{\tau}
    - (1 - \beta) \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k}
    \left[ \DP{(\overline{\rho} \overline{\theta} w^{\tau})}{z}
    \right]_{i,k} 
    - \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\theta} } \right]_{k}
      \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,k} 
    + F_{i,k}, \nonumber \\
 &amp; &amp; (k = 2, 3, \dots , k_m -1 )
  \Deqlab{TrisanG-Deriv22} \\
D_{i,k_m}
 &amp;=&amp; [ \Pi^{\prime} ]_{i,k_m}^{\tau}
    - (1 - \beta) \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{k_m}
    \left[ \DP{(\overline{\rho} \overline{\theta} w^{\tau})}{z}
    \right]_{i,k_m} 
    - \left[ \frac{\overline{c^2} \Delta \tau}{c_{pv}
          \overline{\theta} } \right]_{k_m}
      \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,k_m} \nonumber \\
 &amp; &amp;   + F_{i,k_m}
     + \beta \left[ \frac{\overline{c^2} (\Delta \tau)^2}{c_{pv}
        \overline{\rho} \overline{\theta}^2} \right]_{k_m}
      \Dinv{\Delta z} (c_{pv} \overline{\rho} \overline{\theta}^2)_{i,k_m(w)}
      E_{i,k_m(w)}. 
  \Deqlab{TrisanG-Deriv23}
\end{eqnarray}
が得られる. 

\clearpage

\chapter*{付録 B: 負の雲密度の発生を防ぐ処置に関して}

\addcontentsline{toc}{chapter}{付録 B: 負の雲密度の発生を防ぐ処置に関して}

\markboth{負の雲密度の発生を防ぐ処置に関して}
         {負の雲密度の発生を防ぐ処置に関して}

\setcounter{equation}{0}
\setcounter{chapter}{2}
\renewcommand{\thechapter}{\Alph{chapter}}
\renewcommand{\theequation}{\thechapter.\arabic{equation}}

%\subsection*{fillnegative の保存性の検証}

本もデルでは中央差分を用いることで生じる負の雲密度を防ぐ為, 全質量が保存
されるように負の雲密度の穴埋めを行なっている. 
雲密度に対する穴埋めを行なう関数 FillNegativeMMC の中の式を具体的に書き
下すと, 以下のようになる. 
\begin{eqnarray}
QSUMPN(i,k) 
 &amp;=&amp; 
  1.0 / 
   \nonumber \\
 &amp; &amp;
    \{
     0.75 \ast [  
      {\rm max}( 0.0, \rho_s(i-1,k) )
       \nonumber \\ 
 &amp; &amp;  + {\rm max}( 0.0, \rho_s(i+1,k) )  
       \nonumber \\
 &amp; &amp;  + {\rm max}( 0.0, \rho_s(i,k-1) ) 
       \nonumber \\
 &amp; &amp;  + {\rm max}( 0.0, \rho_s(i,k+1) ) 
      ] 
       \nonumber \\                      
 &amp; &amp;  + 0.25 \ast [   
      {\rm max}( 0.0, \rho_s(i-2,k) )  
       \nonumber \\
 &amp; &amp;  + {\rm max}( 0.0, \rho_s(i+2,k) )  
       \nonumber \\
 &amp; &amp;  + {\rm max}( 0.0, \rho_s(i,k-2) )  
       \nonumber \\
 &amp; &amp;  + {\rm max}( 0.0, \rho_s(i,k+2) )  
      ]
    \} 
     \Deqlab{A} \\
DQFILL(i,k) 
 &amp;=&amp;                            
   - {\rm min}( 0.0, \rho_s(i,k) )
       \nonumber \\ 
 &amp; &amp; + {\rm max}( 0.0, \rho_s(i,k) )  
       \nonumber \\           
 &amp; &amp;   * \{
      0.75 * [
        {\rm min}( 0.0, \rho_s(i-1,k) ) * QSUMPN(i-1,k)            
       \nonumber \\ 
 &amp; &amp;    + {\rm min}( 0.0, \rho_s(i+1,k) ) * QSUMPN(i+1,k)            
       \nonumber \\ 
 &amp; &amp;    + {\rm min}( 0.0, \rho_s(i,k-1) ) * QSUMPN(i,k-1)            
       \nonumber \\ 
 &amp; &amp;    + {\rm min}( 0.0, \rho_s(i,k+1) ) * QSUMPN(i,k+1)            
      ]                        
       \nonumber \\ 
 &amp; &amp; + 0.25 * [
        {\rm min}( 0.0, \rho_s(i-2,k) ) * QSUMPN(i-2,k)            
       \nonumber \\ 
 &amp; &amp;    + {\rm min}( 0.0, \rho_s(i+2,k) ) * QSUMPN(i+2,k)            
       \nonumber \\ 
 &amp; &amp;    + {\rm min}( 0.0, \rho_s(i,k-2) ) * QSUMPN(i,k-2)            
       \nonumber \\ 
 &amp; &amp;    + {\rm min}( 0.0, \rho_s(i,k+2) ) * QSUMPN(i,k+2)            
      ]                        
                 \} 
      \nonumber \\
     \Deqlab{B} \\
\hat{\rho_s} (i,k) &amp;=&amp; \rho_s (i,k) + DQFILL (i,k)
     \Deqlab{C}. 
\end{eqnarray}
但し, $\rho_s$ は穴埋め前の雲密度, $\hat{\rho_s}$ は穴埋め後の雲密度,
$DQFILL(i,k)$ は穴埋めする量である. 
例えばある一点のみで雲密度が負となった場合, この関数は負の点がゼロとなる
ように隣接する点と 2 つ隣の点から雲を削って埋める. 
その際, 隣接する各点は埋める量全体の $3/16$ を負担し, 2 つ隣の各点は埋め
る量全体の $1/16$ を負担する. 

以下では, 単純な場合として, ある一点 $(i,k) = (M,N)$ のみ穴埋め前の雲密
度が負である場合に, 穴埋め前後での全雲質量が保存されることを示す. 
ある一点 $(i,k) = (M,N)$ のみ穴埋め前の雲密度が負であるとき, 式 \Deqref{A},
\Deqref{B} より
$(i,k) = (M-2,N), (M-1,N), (M,N), (M+1,N), (M+2,N), (M,N-2), (M,N-1),
(M,N+1), (M,N+2)$ 以外の点では $DQFILL(i,k) = 0$ となる. 
以上の 9 点での $DQFILL(i,k)$ を書き下すと, 
\begin{eqnarray}
DQFILL (M-2,N)
 &amp;=&amp; 0.25 * \rho_s(M-2,N) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{D} \\
DQFILL (M-1,N)
 &amp;=&amp; 0.75 * \rho_s(M-1,N) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{E} \\
DQFILL (M,N)
 &amp;=&amp; - \rho_s (M,N)
  \nonumber \\ \Deqlab{F} \\
DQFILL (M+1,N)
 &amp;=&amp; 0.75 * \rho_s(M+1,N) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{G} \\
DQFILL (M+2,N)
 &amp;=&amp; 0.25 * \rho_s(M+2,N) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{H} \\
DQFILL (M,N-2)
 &amp;=&amp; 0.25 * \rho_s(M,N-2) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{I} \\
DQFILL (M,N-1)
 &amp;=&amp; 0.75 * \rho_s(M,N-1) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{J} \\
DQFILL (M,N+1)
 &amp;=&amp; 0.75 * \rho_s(M,N+1) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{K} \\
DQFILL (M,N+2)
 &amp;=&amp; 0.25 * \rho_s(M,N+2) \rho_s(M,N) * QSUMPN(M,N), 
  \nonumber \\ \Deqlab{L}
\end{eqnarray}
となる. 
但し
\begin{eqnarray}
QSUMPN (M,N)
 &amp;=&amp;
    \{
     0.75 \ast [  
      \rho_s(M-1,N)
      + \rho_s(M+1,N)  
       \nonumber \\
 &amp; &amp;  + \rho_s(M,N-1) 
      + \rho_s(M,N+1)
      ] 
       \nonumber \\                      
 &amp; &amp;  + 0.25 \ast [   
      \rho_s(M-2,N)  
      + \rho_s(M+2,N)
       \nonumber \\
 &amp; &amp;  + \rho_s(M,N-2)  
      + \rho_s(M,N+2)  
      ]
    \}^{-1}
 \Deqlab{M}
\end{eqnarray}
である. 
従って
\begin{eqnarray}
\sum_{i} \sum_{k} \hat{\rho_s} (i, k)
 &amp;=&amp;
  \sum_{i} \sum_{k} 
   \left[ \rho_s (i,k) + DQFILL (i,k) \right]
   \nonumber \\
 &amp;=&amp;
  \sum_{i} \sum_{k} \rho_s (i,k)
   \nonumber \\
 &amp; &amp;
   + DQFILL (M-2,N)
   + DQFILL (M-1,N)
   \nonumber \\
 &amp; &amp;
   + DQFILL (M,N)
   + DQFILL (M+1,N)
   \nonumber \\
 &amp; &amp;
   + DQFILL (M+2,N)
   + DQFILL (M,N-2)
   \nonumber \\
 &amp; &amp;
   + DQFILL (M,N-1)
   + DQFILL (M,N+1)
   \nonumber \\
 &amp; &amp;
   + DQFILL (M,N+2)
   \nonumber \\
 &amp;=&amp;
   \sum_{i} \sum_{k} \rho_s (i,k)
   - \rho_s (M,N)
    \nonumber \\
 &amp; &amp;
   + \rho_s (M,N) * QSUMPN(M,N) 
    \nonumber \\
 &amp; &amp;
   \{
     0.75 
     * [ \rho_s(M-1,N) + \rho_s(M+1,N) 
    \nonumber \\
 &amp; &amp;     + \rho_s(M,N+1) + \rho_s(M,N+1) ]
    \nonumber \\
 &amp; &amp; + 0.25 
     * [ \rho_s(M-2,N) + \rho_s(M+2,N) 
    \nonumber \\
 &amp; &amp;     + \rho_s(M,N+2) + \rho_s(M,N+2) ] 
   \}  
    \nonumber \\
 &amp;=&amp; 
     \sum_{i} \sum_{k} \rho_s (i,k)
     - \rho_s (M,N)
    \nonumber \\
 &amp; &amp;
   + \rho_s (M,N) * QSUMPN(M,N) * QSUMPN(M,N)^{-1} 
    \nonumber \\
 &amp;=&amp;
     \sum_{i} \sum_{k} \rho_s (i,k)
  \Deqlab{N} 
\end{eqnarray}
となるので, 穴埋め前後での質量は相等しい. 

同様に雲密度が負である点が複数ある場合についても,穴埋め前後で質量は保存
されると考えられる. 


\end{document}
</pre></body></html>