<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 山下 達也 改訂


%\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

\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[2011/12/14]{2011 年 12 月 14 日}
\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}. 
     \Deqlab{RisanAA}
\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}
%
{\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}
\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}.
 \Deqlab{RisanAC}
\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}
    - \frac{3 g {C_m}^2 l^2}{2 \overline{\theta}_{i,k}}
    \left( \DP{\theta^{\prime}}{z} \right)_{i,k} \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}
とする. 


\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
  + {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} \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; + {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} の導出については北守(2006) の付録 B を参照され
% たい. }. 
%
\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_{km - 1} &amp; B_{km} \\
  0 &amp; 0 &amp; \cdots &amp; C_{km - 1} &amp; A_{km} 
 \end{pmatrix}
 \begin{pmatrix}
  \Pi^{\prime}_{1,1} &amp; \Pi^{\prime}_{2,1} &amp; \cdots &amp; \Pi^{\prime}_{km-1,1} &amp;
  \Pi^{\prime}_{km,1} \\
  \Pi^{\prime}_{1,2} &amp; \Pi^{\prime}_{2,2} &amp; \cdots &amp; \Pi^{\prime}_{km-1,2} &amp;
  \Pi^{\prime}_{km,2} \\
  \colon &amp; \colon &amp; &amp; \colon &amp; \colon \\
  \Pi^{\prime}_{1,km-1} &amp; \Pi^{\prime}_{2,km-1} &amp; \cdots &amp;
  \Pi^{\prime}_{km-1,km-1} &amp; \Pi^{\prime}_{km,km-1} \\
  \Pi^{\prime}_{1,km} &amp; \Pi^{\prime}_{2,km} &amp; \cdots &amp;
  \Pi^{\prime}_{km-1,km} &amp; \Pi^{\prime}_{km,km} 
 \end{pmatrix} ^{\tau + \Delta \tau}
 \nonumber
\end{eqnarray}
\begin{eqnarray}
=
 \begin{pmatrix}
  D_{1,1} &amp; D_{2,1} &amp; \cdots &amp; D_{km-1,1} &amp; D_{km,1} \\
  D_{1,2} &amp; D_{2,2} &amp; \cdots &amp; D_{km-1,2} &amp; D_{km,2} \\
  \colon &amp; \colon &amp; &amp; \colon &amp; \colon \\
  D_{1,km-1} &amp; D_{2,km-1} &amp; \cdots &amp; D_{km-1,km-1} &amp; D_{km,km-1} \\
  D_{1,km} &amp; D_{2,km} &amp; \cdots &amp; D_{km-1,km} &amp; D_{km,km}
 \end{pmatrix} ^{\tau}
  \Deqlab{TrisanH}
\end{eqnarray}
ここで $km$ は鉛直方向の格子点の総数である. 
\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 , km -1 ) \Deqlab{TrisanJ} \\
A_{km} &amp;=&amp; 1 + \beta^2 \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{km}
        \Dinv{(\Delta z)^2} 
        ( c_{pv} \overline{\rho} \overline{\theta}^2)_{i,km-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 , km -1 )
  \Deqlab{TrisanO} \\ 
D_{i,km}
 &amp;=&amp; [ \Pi^{\prime} ]_{i,km}^{\tau}
    - (1 - \beta) \left[ \frac{\overline{c}^2 \Delta \tau}{c_{pv}
          \overline{\rho} \overline{\theta}^2} \right]_{km}
    \left[ \DP{(\overline{\rho} \overline{\theta} w^{\tau})}{z}
    \right]_{i,km} 
    - \left[ \frac{\overline{c}^2 \Delta \tau}{c_{pv}
          \overline{\theta} } \right]_{km}
      \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,km} \nonumber \\
 &amp; &amp;   + F_{i,km}
     + \beta \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_{pv}
        \overline{\rho} \overline{\theta}^2} \right]_{km}
      \Dinv{\Delta z} (c_{pv} \overline{\rho} \overline{\theta}^2)_{i,km(w)}
      E_{i,km(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}} 
  \right)_{i,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\}_{i,k(w)}
   \nonumber \\
 &amp; &amp; 
   + \frac{\overline{c}^2}{c_{pv} \overline{\theta}} 
   \left(\frac{L}{c_{pv} \overline{T}} -1 \right) 
   \frac{M_{cond}^{\tau}}{\overline{\rho}}
 \Delta \tau
 + \frac{ \overline{c}^2 }{c_{pv} \overline{\theta}^2
       \overline{\Pi} } ( Q_{rad} + Q_{dis} ) \Delta \tau
 \Big]_{i,k}
  \Deqlab{TrisanR}
\end{eqnarray}
である. 
\Deqref{TrisanH} の左辺の係数行列は 3 重対角行列となっているので, Thomas
法を用いて時刻 $ \tau + \Delta \tau$ での $\Pi^{\prime}$ の値を求めるこ
とができる(Thomas, 1949). 
本モデルでは計算ライブラリ LAPACK を用いて \Deqref{TrisanH} を解いている. 

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

雲密度の時間積分を行なうに当たっては, 
先ず音波・凝結に関連しない項と雲粒落下項を計算し, 
次に負の雲密度の発生を防ぐ処置を行い, 
最後に音波・凝結に関連する項を計算する
\footnote{
負の雲密度の発生を防ぐ処置の具体的な手順については, 付録 A を参照された
い. 
}. 
即ち, \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}
 \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;
- \left[
\frac{3 g {C_m}^2 l^2}{2 \overline{\theta}} 
      \left( \DP{\theta^{\prime}}{z} \right)_{i,k(w)} 
 \right]_{i,k}^{t - \Delta t} \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: 負の雲密度の発生を防ぐ処置に関して}

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

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

\setcounter{equation}{0}
\setcounter{chapter}{1}
\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>