%修論 tex ファイル %作成開始日 2008/11/16 %作成者 山下 達也 %\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 --- 火星大気中の CO$_2$ 氷雲を伴う対流に関する \\ \Large 基礎的数値実験 --- \\ \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[2009/02/02]{2009 年 2 月 2 日} \Dfile{} \maketitle %\begin{titlepage} % %\begin{center} %{\huge 火星湿潤対流(仮)} %\end{center} % %\begin{center} %Martian Moist Convection \\ %\end{center} % %\begin{center} %\Large{ 北海道大学理学院宇宙理学専攻 \\ 惑星物理学研究室修士課程 2 年} %\end{center} % %\begin{center} %Planetary Physics Laboratory, %Department of Cosmosciences, %Faculty of Science, Hokkaido University %\end{center} % %\begin{center} %\huge{ 山下 達也} %\end{center} % %\begin{center} %YAMASHITA Tatsuya %\end{center} % %\begin{center} %\Large{ 2009年2月2日} %\end{center} % %\end{titlepage} \markboth{離散化} {離散化} \chapter{離散化} 第 3 章では主成分凝結を考慮した 2 次元準圧縮方程式系の離散化について述べ る. \section{離散化の概要} 第 3.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{空間方向の離散化} 第 3.3 節では空間微分の離散化方法とその為に必要な平均操作について説明し た上で, 準圧縮方程式系の空間方向の離散化について述べる. \subsection{平均操作} 第 3.2 節で述べたようにスカラー量の格子点とベクトル量の格子点は互いに半 格子ずつずれている. 数値計算を行なう上でベクトル量をスカラー量の格子点で評価したり, 或いはス カラー量をベクトル量の格子点で評価する必要がある. その際, 平均操作を行なうことによって半格子ずれた点での値を評価することと する. 以下, 計算に必要と平均操作を示す. 但し $\phi$, $H$, $V$ はそれぞれスカラー量, ベクトル量の水平成分, ベクト ル量の鉛直成分を表す. また下付き添字は格子点位置を表している. \begin{eqnarray} \phi_{i(u), k} &=& \frac{\phi_{i+1, k} + \phi_{i, k}}{2}, \Deqlab{RisanA} \\ \phi_{i, k(w)} &=& \frac{\phi_{i, k+1} + \phi_{i, k}}{2}, \Deqlab{RisanB} \\ \phi_{i(u), k(w)} &=& \frac{\phi_{i, k} + \phi_{i+1, k} + \phi_{i, k+1} + \phi_{i+1, k+1}}{4}, \Deqlab{RisanC} \\ H_{i, k} &=& \frac{H_{i(u), k} + H_{i-1(u), k}}{2}, \Deqlab{RisanD} \\ H_{i, k(w)} &=& \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)} &=& \frac{H_{i(u), k+1} + H_{i(u), k}}{2}, \Deqlab{RisanF} \\ V_{i, k} &=& \frac{V_{i, k(w)} + V_{i, k-1(w)}}{2}, \Deqlab{RisanG} \\ V_{i(u), k} &=& \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)} &=& \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} &=& \frac{\phi_{i+1, k} - \phi_{i, k}}{\Delta x}, \Deqlab{RisanJ} \\ \left( \DP{\phi}{z} \right)_{i, k(w)} &=& \frac{\phi_{i, k+1} - \phi_{i, k}}{\Delta z}, \Deqlab{RisanK} \\ \left( \DP{H}{x} \right)_{i,k} &=& \frac{H_{i(u), k} - H_{i-1(u), k}}{\Delta x}, \Deqlab{RisanL} \\ \left( \DP{H}{z} \right)_{i(u), k(w)} &=& \frac{H_{i(u), k+1} - H_{i(u), k}}{\Delta z}, \Deqlab{RisanM} \\ \left( \DP{V}{x} \right)_{i(u), k(w)} &=& \frac{V_{i+1, k(w)} - V_{i, k(w)}}{\Delta x}, \Deqlab{RisanN} \\ \left( \DP{V}{z} \right)_{i, k} &=& \frac{V_{i, k(w)} - V_{i, k-1(w)}}{\Delta x}, \Deqlab{RisanO} \\ \left( \DP{\psi}{x} \right)_{i, k(w)} &=& \frac{\psi_{i(u), k(w)} - \psi_{i-1(u), k(w)}}{\Delta x}, \Deqlab{RisanP} \\ \left( \DP{\psi}{z} \right)_{i(u), k} &=& \frac{\psi_{i(u), k(w)} - \psi_{i(u), k-1(w)}}{\Delta z}. \Deqlab{RisanQ} \end{eqnarray} 以下, 4 次精度の中心差分を用いた微分操作を示す. \begin{eqnarray} \left( \DP{\phi}{x} \right)_{i(u), k} &=& \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)} &=& \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} &=& \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)} &=& \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)} &=& \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} \\ \left( \DP{V}{z} \right)_{i, k} &=& \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)} &=& \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} \notag \\ \end{eqnarray} \begin{eqnarray} \left( \DP{\psi}{z} \right)_{i(u), k} &=& \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} \notag \\ \end{eqnarray} \subsection{準圧縮方程式系の空間方向の離散化} 第 3.1 節, 第 3.2 節の結果を用いて準圧縮方程式系を空間方向に離散化すると, 以下のように書ける. \begin{eqnarray} \DP{u_{i(u), k}}{t} &=& - u_{i(u), k} \left( \DP{u}{x} \right)_{i(u), k} - w_{i(u), k} \left( \DP{u}{z} \right)_{i(u), k} \notag \\ & & - c_p \overline{\theta}_{i(u), k} \left( \DP{\Pi}{x} \right)_{i(u), k} + \left[ D_u \right]_{i(u), k}, \Deqlab{RisanAA} \\ \DP{w_{i, k(w)}}{t} &=& - u_{i, k(w)} \left( \DP{w}{x} \right)_{i, k(w)} - w_{i, k(w)} \left( \DP{u}{z} \right)_{i, k(w)} \notag \\ & & - c_p \overline{\theta}_{i, k(w)} \left( \DP{\Pi}{z} \right)_{i, k(w)} + \left[ D_w \right]_{i, k(w)}, \Deqlab{RisanAB} \\ \DP{\theta_{i,k}}{t} &=& - u_{i,k} \left( \DP{\theta}{x} \right)_{i,k} - w_{i,k} \left( \DP{\theta}{z} \right)_{i,k} \notag \\ & & + \frac{\overline{\theta}_{i,k}}{\overline{T}_{i,k}} \frac{L [M_{cond}]_{i,k}}{\overline{\rho}_{i,k} c_p} + \frac{\overline{\theta}_{i,k}}{\overline{T}_{i,k}} \left( [Q_{rad}]_{i,k} + [Q_{dis}]_{i,k} \right) + \left[ D_{\theta} \right]_{i,k}, \Deqlab{RisanAC} \\ \DP{\Pi_{i,k}}{t} &=& - \frac{\overline{c}^2_{i,k}}{c_p \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] \notag \\ & & + \frac{\overline{c}^2_{i,k}}{c_p \overline{\rho}_{i,k} \overline{\theta}_{i,k}} \left( \frac{L}{c_p \overline{T}_{i,k} } -1 \right) \left[ M_{cond} \right]_{i,k}, \Deqlab{RisanAD} \\ \DP{{\rho_s}_{i,k}}{t} &=& - \left( \rho_s u \right)_{i,k} - \left( \rho_s w \right)_{i,k} + \left[ M_{cond} \right]_{i,k} + \left[ M_{fall} \right]_{i,k} + \left[ D_{\rho_s} \right]_{i,k}. \Deqlab{RisanAE} \end{eqnarray} \section{時間方向の離散化} 第 3.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{AA}, \Deqref{AB} の $\Pi^{\prime}$ を \begin{equation} \Pi^{\prime} - \alpha \left( \DP{u}{x} + \DP{w}{z} \right) \Deqlab{SoundDecay1} \end{equation} で置き換える. 音波を選択的に減衰させる為には, 係数 $\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 = 5.0 \times 10^{-4} \times \frac{{\rm min} [(\Delta x)^2, (\Delta z)^2]}{\Delta t} \Deqlab{NumDisp2} \end{equation} とする \footnote{ 北守(2006)には $\nu_H = \nu_V = 1.0 \times 10^{-4} \times {\rm min} [(\Delta x)^2, (\Delta z)^2] / \Delta t$ とすると記述されているが, 北守(2006)の計算設定並びに結果をまとめたサイト \\ http://www.gfd-dennou.org/arch/deepconv/arare/sample/2006-02-09\_kitamo/ \\ によると \Deqref{NumDisp2} の値を用いたようである. }. \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{AA} を離散化すると \begin{equation} u_{i(u),k}^{\tau + \Delta \tau} = u_{i(u),k}^{\tau} - \left( c_p \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{AB} において音波減衰項, 圧力勾配項はそれぞれ前進差分, 後退差分で 離散化する. \Deqref{AH} において水平方向のフラックス項, 鉛直方向のフラックス項はそれ ぞれ前進差分, 後退差分で離散化する. \Deqref{AB} 及び \Deqref{AH} を離散化すると \begin{eqnarray} w_{i,k(w)}^{\tau + \Delta \tau} &=& w_{i,k(w)}^{\tau} - c_p \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 \notag \\ & & + {F_w}_{i,k(w)}^{t} \Delta \tau, \Deqlab{TrisanE} \\ \left[ \Pi^{\prime} \right]_{i,k}^{\tau + \Delta \tau} &=& \left[ \Pi^{\prime} \right]_{i,k}^{\tau} - \left[ \frac{ \overline{c}^2 }{c_p \overline{\rho} \overline{\theta}^2 } \DP{}{x} (\overline{\rho} \overline{\theta} u^{\tau + \Delta \tau})\right]_{i,k} \Delta \tau \notag \\ & & - \beta \left[ \frac{ \overline{c}^2 }{c_p \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_p \overline{\rho} \overline{\theta}^2 } \DP{}{z} (\overline{\rho} \overline{\theta} w^{\tau})\right]_{i,k} \Delta \tau \notag \\ & & + \left[ \frac{\overline{c}^2}{c_p \overline{\theta}} \left( \frac{L}{c_p \overline{T}} - 1 \right) \frac{M_{cond}^{\tau}}{\overline{\rho}} \right]_{i,k} \Delta \tau \Deqlab{TrisanF} \end{eqnarray} となる. 但し $\beta$ は後退差分における重み係数を表し, クランク・ニコルソン法の 場合 $\beta = 0.5$ とする. また $F_w$ は音波, 凝結に関連しない項で \begin{equation} {F_w}_{i,k(w)}^{t} = - \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)}} + \left[ D_w \right]_{i,k(w)}^{t - \Delta t} \Deqlab{TrisanG} \end{equation} である. \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 & B_2 & \cdots & 0 & 0 \\ C_1 & A_2 & \cdots & 0 & 0 \\ \colon & \colon & & \colon & \colon \\ 0 & 0 & \cdots & A_{km - 1} & B_{km} \\ 0 & 0 & \cdots & C_{km - 1} & A_{km} \end{pmatrix} \begin{pmatrix} \Pi^{\prime}_{1,1} & \Pi^{\prime}_{2,1} & \cdots & \Pi^{\prime}_{km-1,1} & \Pi^{\prime}_{km,1} \\ \Pi^{\prime}_{1,2} & \Pi^{\prime}_{2,2} & \cdots & \Pi^{\prime}_{km-1,2} & \Pi^{\prime}_{km,2} \\ \colon & \colon & & \colon & \colon \\ \Pi^{\prime}_{1,km-1} & \Pi^{\prime}_{2,km-1} & \cdots & \Pi^{\prime}_{km-1,km-1} & \Pi^{\prime}_{km,km-1} \\ \Pi^{\prime}_{1,km} & \Pi^{\prime}_{2,km} & \cdots & \Pi^{\prime}_{km-1,km} & \Pi^{\prime}_{km,km} \end{pmatrix} ^{\tau + \Delta \tau} \notag \end{eqnarray} \begin{eqnarray} = \begin{pmatrix} D_{1,1} & D_{2,1} & \cdots & D_{km-1,1} & D_{km,1} \\ D_{1,2} & D_{2,2} & \cdots & D_{km-1,2} & D_{km,2} \\ \colon & \colon & & \colon & \colon \\ D_{1,km-1} & D_{2,km-1} & \cdots & D_{km-1,km-1} & D_{km,km-1} \\ D_{1,km} & D_{2,km} & \cdots & D_{km-1,km} & D_{km,km} \end{pmatrix} \Deqlab{TrisanH} \end{eqnarray} ここで $km$ は鉛直方向の格子点の総数である. \Deqref{TrisanH} の左辺の係数行列の各成分は以下のように表される. \begin{eqnarray} A_1 &=& 1 + \beta^2 \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_p \overline{\rho} \overline{\theta}^2} \right]_1 \Dinv{(\Delta z)^2} ( c_p \overline{\rho} \overline{\theta}^2)_{i,1(w)}, \Deqlab{TrisanI} \\ A_k &=& 1 + \beta^2 \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_p \overline{\rho} \overline{\theta}^2} \right]_k \Dinv{(\Delta z)^2} [ ( c_p \overline{\rho} \overline{\theta}^2)_{i,k(w)} + ( c_p \overline{\rho} \overline{\theta}^2)_{i,k-1(w)} ], \notag \\ & & (k = 2, 3, \dots , km -1 ) \Deqlab{TrisanJ} \\ A_{km} &=& 1 + \beta^2 \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_p \overline{\rho} \overline{\theta}^2} \right]_{km} \Dinv{(\Delta z)^2} ( c_p \overline{\rho} \overline{\theta}^2)_{i,km-1(w)}, \Deqlab{TrisanK} \\ B_k &=& - \beta^2 \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_p \overline{\rho} \overline{\theta}^2} \right]_{k-1} \Dinv{(\Delta z)^2} ( c_p \overline{\rho} \overline{\theta}^2)_{i,k-1(w)}, \Deqlab{TrisanL} \\ C_k &=& - \beta^2 \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_p \overline{\rho} \overline{\theta}^2} \right]_{k+1} \Dinv{(\Delta z)^2} ( c_p \overline{\rho} \overline{\theta}^2)_{i,k(w)}, \Deqlab{TrisanM} \\ D_{i,1} &=& [ \Pi^{\prime} ]_{i,1}^{\tau} - (1 - \beta) \left[ \frac{\overline{c}^2 \Delta \tau}{c_p \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_p \overline{\rho} } \right]_{1} \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,1} + F_{i,1} \notag \\ & & - \beta \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_p \overline{\rho} \overline{\theta}^2} \right]_{1} \Dinv{\Delta z} (c_p \overline{\rho} \overline{\theta}^2)_{i,0(w)} E_{i,0(w)}, \Deqlab{TrisanN} \end{eqnarray} \begin{eqnarray} D_{i,k} &=& [ \Pi^{\prime} ]_{i,k}^{\tau} - (1 - \beta) \left[ \frac{\overline{c}^2 \Delta \tau}{c_p \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_p \overline{\rho} } \right]_{k} \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,k} + F_{i,k}, \notag \\ & & (k = 2, 3, \dots , km -1 ) \Deqlab{TrisanO} \\ D_{i,km} &=& [ \Pi^{\prime} ]_{i,km}^{\tau} - (1 - \beta) \left[ \frac{\overline{c}^2 \Delta \tau}{c_p \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_p \overline{\rho} } \right]_{km} \left[ \DP{u^{\tau + \Delta \tau}}{x} \right]_{i,km} \notag \\ & & + F_{i,km} + \beta \left[ \frac{\overline{c}^2 (\Delta \tau)^2}{c_p \overline{\rho} \overline{\theta}^2} \right]_{km} \Dinv{\Delta z} (c_p \overline{\rho} \overline{\theta}^2)_{i,km(w)} E_{i,km(w)}. \Deqlab{TrisanP} \end{eqnarray} 但し \begin{eqnarray} E_{i,k(w)} &=& - (1 - \beta )\left[ \DP{[\Pi^{\prime}]^{\tau}}{z} \right]_{i,k(w)} + \left[ \frac{F^{t}_{w}}{c_p \overline{\theta}} \right]_{i,k(w)}, \Deqlab{TrisanQ} \\ F_{i,k} &=& \Big[ - \Delta \tau \beta \frac{\overline{c}^2}{c_p \overline{\rho} \overline{\theta}} \DP{}{z} \Big\{ \overline{\rho} \overline{\theta} \Big( w^{\tau} - \Delta \tau (1 - \beta) c_p \overline{\theta} \DP{[\Pi^{\prime}]^{\tau}}{z} + F_{w}^{t} \Delta \tau \Big) \Big\} \notag \\ & & + \frac{\overline{c}^2}{c_p \overline{\theta}} \left(\frac{L}{c_p \overline{T}} -1 \right) \frac{M_{cond}^{\tau}}{\overline{\rho}} \Delta \tau \Big] \Deqlab{TrisanR} \end{eqnarray} である. \Deqref{TrisanH} の左辺の係数行列は 3 重対角行列となっているので, Thomas 法を用いて時刻 $ \tau + \Delta \tau$ での $\Pi^{\prime}$ の値を求めるこ とができる(Thomas, 1949). 本モデルでは計算ライブラリ LAPACK を用いて \Deqref{TrisanH} を求めている. {\bf 熱力学の式の離散化} \Deqref{AC}を離散化すると \begin{equation} {\theta^{\prime}}_{i,k}^{\tau + \Delta \tau} = {\theta^{\prime}}_{i,k}^{\tau} + {F_{\theta}}_{i,k}^{t} \Delta \tau + \frac{\overline{{\theta^{\prime}}}_{i,k}}{\overline{T}_{i,k}} \frac{L {M_{cond}^{\tau}}_{i,k}}{c_p \overline{\rho}_{i,k}} \Delta \tau \Deqlab{TrisanS} \end{equation} となる. ここで $F_{\theta}$ は音波・凝結に関連しない項であり \begin{eqnarray} {F_{\theta}}_{i,k}^{t} &=& - 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} \notag \\ & & + \frac{\overline{{\theta^{\prime}}}_{i,k}}{\overline{T}_{i,k}} ( Q_{rad} + Q_{dis})_{i,k}^{t} + {D_{\theta}}_{i,k}^{t - \Delta t} \Deqlab{TrisanT} \end{eqnarray} と表される. {\bf 雲密度の式の離散化} \Deqref{AE}を離散化すると \begin{equation} {\rho_s}_{i,k}^{\tau + \Delta \tau} = {\rho_s}_{i,k}^{\tau} + {M_{cond}}_{i,k}^{\tau} \Delta \tau + {F_{\rho_s}}_{i,k}^{t} \Delta \tau \Deqlab{TrisanU} \end{equation} となる. ここで $F_{\rho_s}$ は音波・凝結に関連しない項であり, \begin{equation} {F_{\rho_s}}_{i,k}^{t} = - \left[ \DP{}{x} ({\rho_s}^t u^t) \right]_{i,k} - \left[ \DP{}{z} ({\rho_s}^t w^t) \right]_{i,k} + {M_{fall}}_{i,k}^{t} + {D_{\rho_s}}_{i,k}^{t} \Deqlab{TrisanV} \end{equation} と表される. {\bf 乱流拡散係数の式の離散化} \Deqref{TurbG}をリープフロッグ法を用いて離散化すると \begin{equation} {K_m}_{i,k}^{t+ \Delta t} = {K_m}_{i,k}^{t+ \Delta t} + 2 \Delta [F_{K_m}]_{i,k}^{t} \Deqlab{TrisanW} \end{equation} となる. ここで \begin{eqnarray} [F_{K_m}]_{i,k}^{t} &=& - \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 \notag \\ & & - \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} \notag \\ & & + {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} \notag \\ & & + \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 \notag \\ & & - \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] \notag \\ & & + \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] \notag \\ & & + \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} \notag \\ & & - \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 の時間フィルター} 3.4 節で述べた通り, 長い時間ステップの計算ではリープフロッグ法を用いてい る. リープフロッグ法では時刻 $t - \Delta t$ の値を用いて時刻 $t + \Delta$ の 値を求める為, 隣接する時間ステップ間で物理量の値に大きな食い違いや振動が 生じる恐れがある. この問題を回避する為, 長い時間ステップの計算を 1 回行 なう度に Robert(1966), Asselin(1972) が考案した時間フィルターを適用する. 本モデルでは $u$, $w$, $\theta^{\prime}$, $\Pi^{\prime}$, $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$ とした. \end{document}