修士論文

非静力学対流モデルの定式化に関する再検討: 地形を考慮した雲解像モデル

の開発へ向けて Reconsideration of the

formulation of a nonhydrostatic convection model: Toward a cloud resolving model with terrain-following coordinates

高橋こう子 TAKAHASHI Koko

北海道大学大学院理学研究科地球惑星科学専攻

地球惑星流体科学講座地球流体力学研究室

2004年1月30日

\underline{1}

目次

第 1 章はじめに 1

第 2 章乾燥大気における準圧縮モデルの定式化と離散化 3

2.1 準圧縮モデルの定式化 . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 変数と基本場の定義 . . . . . . . . . . . . . . . . . . . . . . 3 2.1.2 直交直線座標系における基礎方程式系 . . . . . . . . . . . . 4 2.1.3 サブグリッドスケールの拡散 . . . . . . . . . . . . . . . . . 5 2.1.4 地形に沿った座標系における基礎方程式系 . . . . . . . . . . 6

2.2 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 基礎方程式系の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3.1 離散化の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.2 音波減衰項 . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.3 数値粘性項 . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.4 Asselin の時間フィルタ . . . . . . . . . . . . . . . . . . . . . 12

第 3 章数値実験 14 3.1 数値実験の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 数値実験 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.2.1 音波の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2.2 音波減衰項の係数を決めるための計算 . . . . . . . . . . . . 21 3.2.3 数値粘性項の係数を決めるための計算 . . . . . . . . . . . . 22 3.2.4 サーマルの計算 . . . . . . . . . . . . . . . . . . . . . . . . . 28

第 4 章まとめと考察 37

謝辞 39

付録 A 乾燥大気における準圧縮モデルの定式化の補足 40 A.1 地形に沿った座標系 . . . . . . . . . . . . . . . . . . . . . . . . . . 40 A.2 サブグリッドスケールの拡散 . . . . . . . . . . . . . . . . . . . . . 42

付録 B 湿潤大気における準圧縮モデルの定式化 44 B.1 直交直線座標系における基礎方程式系 . . . . . . . . . . . . . . . . 44 B.1.1 変数の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 B.1.2 基本場の定義と分離 . . . . . . . . . . . . . . . . . . . . . . 45 B.1.3 基礎方程式系 . . . . . . . . . . . . . . . . . . . . . . . . . . 46

B.2 地形に沿った座標系における基礎方程式系 . . . . . . . . . . . . . . 47

B.2.1 基礎方程式系 . . . . . . . . . . . . . . . . . . . . . . . . . . 47 B.3 サブグリッドスケールの拡散 . . . . . . . . . . . . . . . . . . . . . 49 B.4 雲・降水の微物理過程 . . . . . . . . . . . . . . . . . . . . . . . . . 53

付録 C 乾燥大気における準圧縮モデルの離散化 56 C.1 空間微分の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 C.2 時間積分の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 C.3 乾燥大気における基礎方程式系の離散化 . . . . . . . . . . . . . . . 59

C.3.1 水平方向運動方程式 . . . . . . . . . . . . . . . . . . . . . . 60 C.3.2 鉛直方向運動方程式と圧力方程式 . . . . . . . . . . . . . . . 62 C.3.3 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 C.3.4 サブグリッドスケールの乱流運動エネルギー方程式 . . . . . 70

付録 D 湿潤大気における準圧縮モデルの離散化 72 D.1 空間微分の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 D.2 時間積分の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 D.3 湿潤大気における基礎方程式系の離散化 . . . . . . . . . . . . . . . 73

D.3.1 水平方向運動方程式 . . . . . . . . . . . . . . . . . . . . . . 73 D.3.2 鉛直方向運動方程式と圧力方程式 . . . . . . . . . . . . . . . 76 D.3.3 熱力学の式と混合比の方程式 . . . . . . . . . . . . . . . . . 80 D.3.4 サブグリッドスケールの乱流運動エネルギー方程式 . . . . . 83

D.4 湿潤飽和調節法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

参考文献 87

要旨

本研究では乾燥大気における 3 次元準圧縮モデルの定式化と離散化の再検討を行い, それを基に 2 次元準圧縮の数値モデルを作成し数値実験を行った. 直交直線座標系を用いた定式化は Klemp and Wilhelmson (1978) に従い, 地形に沿った座標系を用いたは定式化と斉藤 (1999) に従った. 離散化の方法は CReSS (坪木他, 2001) とほぼ同じである. 行った数値実験は, 音波の計算, 音波減衰項の係数を決めるための計算, 数値粘性項の係数を決めるための計算, サーマルの計算の 4 つである.

定式化と離散化の再検討を行ったところ, Klemp and Wilhelmson (1978) において示された離散化されたサブグリッドスケールの乱流運動エネルギー方程式には, 離散化する前の式には存在した 2 つの項が無視されていることがわかった. これら 2 つの項のうち, 1 つはサブグリッドスケールの運動は非圧縮であると考えることによって無視することができるものであった. もうひとつの項のオーダは無視されていない別の項と同じオーダであったため, 消去すべきではないことがわかった.

音波を計算する数値実験では, 初期擾乱として水平方向一様または鉛直方向一様なガウス型の圧力をあたえ, 音波の伝播する様子を調べた. その結果, 圧力水平・鉛直速度とエクスナー関数が音速で伝播することが確認された. 音波計算のための時間間隔を変えて鉛直方向の数値解法にともなう音波の減衰の様子を調べたところ, 時間間隔を短くすると音波の減衰が小さくなることが確かめられた.

音波減衰項を決める数値実験では, 初期擾乱としてガウス型の圧力を与え, 音波の伝播する様子を調べた. その結果, 音波減衰項の値は水平空間格子間隔が 1,000 m, 鉛直格子間隔が 500 m, 時間格子間隔が 1 s の場合, 0.025 が適当であることがわかった. 数値粘性項の係数を決めるための数値実験では, 数値粘性項の係数を変えた温位の移流を計算を行い, 適当な係数の値を調べた. その結果, 水平空間格子間隔が 1,000 m, 鉛直格子間隔が 500 m, 時間格子間隔が 1 s の場合にもっとも安定に計算できたのは, 水平方向の係数を 5.0, 鉛直方向の係数を 1.25 とした場合であった. 上記の音波減衰項の係数と数値粘性係数の値を用いてサーマルの数値実験を行った. 下部境界の中心にガウス型の温度擾乱を与えてサーマルの浮上する様子を調べた. 基本場の温度分布として, 高度 5,000 m まで乾燥断熱減率で減少し, それより上は等温となる分布を与えた. その結果, 高度 5,000 m までサーマルが上昇し, その高度で横に広がっていくのが確認できた.

\underline{第 1 章はじめに 1}

第1章はじめに

最近の計算機の発達により, 非静力学モデルの構築が活発に行われるようになった. 非静力学モデルとは, 水平スケールとして 2 〜 2,000 km (メソスケール), 鉛直スケールとして 10 km 前後の気象現象を対象としたものである. 非静力学モデルの例として, 気象研究所の非静力学モデル (NHM; Ikawa and Saito, 1991), オクラホマ大学の Adcanced Regional Prediction System (ARPS; Xue et. al., 1995), 名古屋大学と高度情報科学技術研究機構の CReSS (Cloud Resolving Storm Simulator; 坪木ら, 2001) などが挙げられる. メソスケールの現象に対しては, 鉛直方向において静水圧平衡が成り立つと仮定することができない.

非静力学モデルは音波を含まない非弾性系モデルと含む弾性系モデルに分類される. 非静力学モデルを音波を含まないモデルを考えたのは以下の理由による. 大気に存在する波のうち気象学的に重要なのはロスビー波や重力波であり, ふつう音波には興味がないからである. 近似なしの運動方程式を解くと, 解にロスビー波, 重力波および音波が出てくる. 数値計算する際, これら 3 つの波のうちもっとも位相速度が速い音波によって CFL 条件が決められる. 音波によって決められる時間間隔は非常に短いため, 効率的な数値計算ができない. 非弾性系モデルでは, 大気の圧縮性が無視できると仮定して音波を取り除いた.

非静力学モデルは音波によって 2 つに分類され, さらに連続の式の近似方法によって以下の 4 つに分けられる (斉藤, 1999). 非弾性系として非圧縮モデルと非弾性モデル, 弾性系として準圧縮モデルと完全圧縮系モデルである. 非圧縮モデルでは, 密度が空間的時間的に一定であると仮定し, 非圧縮性流体の連続の式を用いる. 非弾性モデルでは, 密度が空間時間変化しないと仮定する. 準圧縮モデルでは, 密度が基本場からの偏差が小さいと仮定し, 線形化を行った連続の式を用いる. 完全圧縮系モデルでは圧縮性流体の連続の式を用いる.

非圧縮モデルは主に基礎研究や浅い対流を数値計算する場合に用いられる. 非弾性モデルは研究用シミュレーションや再現実験に使われる. しかし, 密度の時間変化を考慮していないため, 密度変化がゆっくりな場合と誤差が大きくなる. 完全圧縮系モデルは最も精度が高いものの, 方程式が複雑であるため数値計算が困難である. 準圧縮モデルの連続の式の精度は線形化を行っているため, 非弾性モデルのものと同じ精度である. しかし, 連続の式に密度の時間変化項を含むため, 密度時間変化がゆっくりな場合にも用いることができる. 準圧縮モデルは非弾性モデルより研究に用いられている.

現在, 様々な準圧縮の数値モデルが開発されている. 例えば先に上げた NHM や ARPS である. 数値シミュレーションや再現実験を行う際, 多くの研究者が既存の数値モデルを用いている. 既存の数値モデルを使用する利点のひとつは, 非常に時間や手間ががかかるプログラミングをする必要がない点である. しかし, ある目的の研究のために開発された数値モデルの場合, その目的に沿った対象に対してしか使うことができないことが多い. 例えば地球の気象を対称にした数値モデルをそのまま地球以外の惑星にそのまま適用できるかは疑問である. 汎用性のある数値モデルであったとしても, 数値モデル内に導入されているパラメタリゼーションが明記されていないこともある. 既存の数値モデルが自分の研究目的にあうとは限らない. そこで, 自ら数値モデルを開発する必要がある. 自ら数値モデルを開発することは, その数値モデルの長所や欠点を把握できるという利点を持つ. 数値モデルの特徴を把握しているため, 新しい概念を導入するときも既存の数値モデルを改良する場合と比べてずっと容易である.

自ら数値モデルを開発することは有用であるものの, 準圧縮系モデルの場合, 基礎方程式だけでなく乱流モデルや微物理過程の定式化や離散化が複雑である. 基礎となる数理モデルの定式化と離散化がまとめられていれば, これを基に数値モデルの開発をさらに進めることができる. したがって, 本研究は準圧縮の数値モデルの開発に必要な基礎知識の集積を目的とした. 非静力学の数値モデル開発のための基礎知識となるために, 以下の 2 点を考慮した. 多くの準圧縮モデルは Klemp and Wilhelmson (1978) を基に開発されている. メソスケールの気象現象は地形の影響を受けやすいため, 地形の効果を考慮する必要がある.

本研究では Klemp and Wilhelmson (1978) を基に定式化と離散化について再考察する. 2.1 節では Klemp and Wilhelmson (1978) を基に直交直線座標系における定式化と, 斉藤 (1999) を基に地形に沿った座標系における定式化を示す. 2.2 節では境界条件を示し, 2.3 節では地形に沿った座標系における基礎方程式の離散化を示す. 第 3 章では, 具体的な値がわからない音波減衰項と数値粘性項の係数を決め, 音波とサーマルの計算結果を示す.

\underline{第 2 章乾燥大気における準圧縮モデルの定式化と離散化 3}

第2章{乾燥大気における準圧縮モデ \atop ルの定式化と離散化}

本章では, 乾燥大気における 3 次元準圧縮モデルの定式化と離散化について述べる. Klemp and Wilhelmson (1978) に従って, 直交直線座標系の定式化を示す. 次に, 斉藤 (1999) に従って地形に沿った座標系の定式化を示す.

2.1 準圧縮モデルの定式化

準圧縮モデルの定式化をまとめると以下のようになる:

• 力学過程の基礎方程式系は非静力学・準圧縮系である. 座標系は 3 次元の地形に沿った座標系である.

• 力学過程の予報変数は速度 3 成分, 圧力, 温位である.

• サブグリッドスケールの乱流モデルは乱流運動エネルギーを予報する 1.5 次のクロージャモデル.

2.1.1 変数と基本場の定義

力学過程の予報変数は x, y, z 方向の速度 u, v, w と無次元化した圧力 (エクスナー関数) Π(≡ (p/p_{0})^{R_{d}/c_{p}}) および温位 θ(≡ T (p_{0}/p)^{R_{d}/c_{p}}) の 5 つである. これら 5 つの予報変数は空間の座標 x, y, z と時間の座標 t である. ここで T は温度, p_{0} は地表面での圧力, c_{p} は定圧比熱, R_{d} は乾燥空気の気体定数である. 密度 ρ は乾燥空気の状態方程式

Π = (\frac{R_{d}}{p_{0}} ρθ)^{R_{d}/c_{v}}  (2.1)

から得られる.

変数と基礎方程式系を基本場とそこからの偏差に分離する. 分離する変数は u, v, w, Π, θ および ρ である. 基本場の変数を上付きバー (¯ ) で, 偏差の変数を上付きプライム (′) で表す. ある変数 ϕ は

ϕ = ϕ̄ + ϕ′  (2.2)

と分離される. 本論文では記述の簡便のため, 本数理モデルではエクスナー関数 Π のずれを π と表す.

基本場は水平一様 (z のみの関数) かつ静止している (ū = v̄ = w̄ = 0) 状態とする. 基本場において, 鉛直方向の運動方程式は静水圧平衡の式となる:

\frac{∂Π̄}{∂z} = −\frac{g}{c_{p}θ̄}.  (2.3)

基本場の密度は

ρ̄ = \frac{p_{0}}{R_{d}} \frac{Π̄^{c_{v}/R_{d}}}{θ̄}  (2.4)

である.

2.1.2 直交直線座標系における基礎方程式系

基礎方程式は x, y, z^{∗} 方向の運動方程式, 熱力学の式, 準圧縮系の圧力方程式である.

x, y, z^{∗} 方向の運動方程式は

\frac{du_{i}}{dt} + c_{p}θ̄\frac{∂π}{∂x_{i}} = δ_{i3}g [\frac{θ}{θ̄} − 1] − ε_{i3k}fu_{k} + D_{u_{i}}  (2.5)

である. i = 1, 2, 3 で, それぞれ x, y および z 成分を表す. f はコリオリパラメータである. D_{u_{i}} は 3 方向の運動方程式におけるサブグリッドスケールの乱流に伴う拡散項である. 以下では, サブグリッドスケールの乱流に伴う拡散項を単に拡散項という. 3 方向の運動方程式の拡散項の表現については 2.1.3 節で述べる.

圧力方程式は

\frac{∂π}{∂t} + \frac{c̄^{2}}{c_{p}ρ̄θ̄^{2}} \frac{∂}{∂x_{j}} (ρ̄θ̄u_{j}) = f_{π},  (2.6)
f_{π} = −u_{j} \frac{∂π}{∂x_{j}} + \frac{R_{d}π}{c_{v}} \frac{∂u_{j}}{∂x_{j}} + \frac{c^{2}}{c_{p}θ^{2}} \frac{dθ}{dt}  (2.7)

である. c は音波で以下のように表される:

c^{2} = \frac{c_{p}R_{d}}{c_{v}} Πθ.  (2.8)

熱力学の式は

\frac{dθ}{dt} = D_{θ}  (2.9)

である. D_{θ} は熱力学の式における拡散項であり, 熱力学の式の拡散項の表現については 2.1.3 節で述べる.

密度は乾燥空気の状態方程式 (2.1) より診断的に得られる:

ρ = \frac{p_{0}}{R_{d}} \frac{Π^{c_{v}/R_{d}}}{θ} .  (2.10)

2.1.3 サブグリッドスケールの拡散

レイノルズ平均を用いて方程式を, グリッドスケールの平均運動とサブグリッドスケールの乱流運動とに分離する. 直交直線座標系における 3 方向の運動方程式 (2.5) と熱力学の式 (2.9) の拡散項 D_{u_{j}}, D_{θ} は, それぞれ,

D_{u_{j}} = −\frac{∂}{∂x_{j}} \overline{u′′_{i}u′′_{j}},  (2.11)
D_{θ} = −\frac{∂}{∂x_{j}} \overline{u′′_{j}θ′′}.  (2.12)

\overline{u′′_{i}u′′_{j}}, \overline{u′′_{i}θ′′} は, それぞれ, 3 方向の運動方程式と熱力学の式のサブグリッドスケールの変動成分である. 各拡散項 D_{u_{j}}, D_{θ} は導出は付録 A.2 節に示す.

Klemp and Wilhelmson (1978) に従って, 変動成分を

\overline{u′′_{i}u′′_{j}} = −K_{m} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} ) + \frac{2}{3}δ_{ij}E,  (2.13)
\overline{u′′_{j}θ′′} = K_{h} \frac{∂θ}{∂x_{j}}  (2.14)

と仮定する. K_{m} は運動量渦拡散係数, K_{h} は熱渦拡散係数である. E は

E ≡ \frac{1}{2} (\overline{u′′} + \overline{v′′} + \overline{w′′})  (2.15)

で与えられるサブグリッドスケールの乱流運動エネルギーである.

K_{m} は以下の時間発展方程式から求められる:

\frac{dK_{m}}{dt} = −\frac{3g}{θ̄} \frac{∂θ}{∂z} + 2C_{m}^{2}l^{2} (\frac{∂u_{j}}{∂x_{j}} )^{2} + \frac{1}{2} |ε_{ijk}|C_{m}^{2}l^{2} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} )^{2}

−\frac{1}{3}δ_{ij}K_{m} (\frac{∂u_{j}}{∂x_{j}} ) + \frac{1}{2}δ_{ij} (\frac{∂^{2}K_{m}^{2}}{∂x^{2}} ) + δ_{ij} (\frac{∂K_{m}}{∂x_{j}} )^{2}

− (\frac{C_{ε}}{2C_{m}l^{2}} ) K^{2}.  (2.16)

ただし K_{m} に関する上式 (2.16) はサブグリッドスケールの乱流運動エネルギーの時間発展方程式から得られる. 式 (2.16) の導出は付録 A.2 に示す. E と K_{m} の関係は

E = (\frac{K_{m}}{C_{m}l} )^{2} .  (2.17)

K_{m} と K_{h} の関係は

K_{h} = 3K_{m}  (2.18)

とする (Deardorff, 1972). ここで l は混合距離で l = (ΔxΔyΔz)^{1/3} である. Deardorff (1975) によると C_{ε} = C_{m} = 0.2 である.

Klemp and Wilhelmson (1978) 中の離散化された運動量渦拡散係数に関する方程式には, (2.16) 中の右辺第 4 項と第 6 項

\frac{1}{3}δ_{ij}K_{m} (\frac{∂u_{j}}{∂x_{j}} ) ,  (2.19)
δ_{ij} (\frac{∂K_{m}}{∂x_{j}} )^{2}  (2.20)

が消去されている. 上の 2 つの項 (2.19) と (2.20) が消去された理由については第 4 章で考える.

2.1.4 地形に沿った座標系における基礎方程式系

図 2.1 のような地形に沿った座標 (x^{∗}, y^{∗}, z^{∗}) は以下のように定義される:

x^{∗} = x,
y^{∗} = y,
z^{∗} = z^{∗}(x, y, z)

かつ,

z^{∗} ≡ \frac{H \{z − z_{s}(x^{∗}, y^{∗})\}}{H − z_{s}(x^{∗}, y^{∗})}.  (2.21)

z_{s}(x^{∗}, y^{∗}) は地表面の高さ, H はモデル領域の高さとする. 地形に沿った座標系における速度 (u^{∗}, v^{∗}, ω) は

u^{∗} = u,  (2.22)
v^{∗} = v,  (2.23)
ω = G^{13}u + G^{23}v + \frac{1}{G^{1/2}} w  (2.24)

である (速度の導出は付録 A.1 節参照). ただし

G^{1/2} = \frac{H − z_{s}}{H} ,  (2.25)
G^{13} = \frac{z^{∗} − H}{H − z_{s}} \frac{∂z_{s}}{∂x} ,  (2.26)
G^{23} = \frac{z^{∗} − H}{H − z_{s}} \frac{∂z_{s}}{∂y}  (2.27)

である. 地形に沿った座標系における空間微分は付録 A.1 節の式 (A.10)〜(A.12) のように定義される.

図 2.1: 直線直交座標系における速度 (u, v, w) と地形に沿った座標系における速度 (u^{∗}, v^{∗}, ω).

地形に沿った座標系の静水圧平衡の式と基礎方程式は以下のようになる. 静水圧平衡の式は

\frac{∂Π̄}{∂z^{∗}} = −\frac{G^{1/2}g}{c_{p}θ̄} .  (2.28)

\underline{第 2 章乾燥大気における準圧縮モデルの定式化と離散化 8}

である. x, y, z^{∗} 方向の運動方程式はそれぞれ以下のように表される:

\frac{∂u}{∂t} + u\frac{∂u}{∂x} + v \frac{∂u}{∂y} + ω \frac{∂u}{∂z^{∗}} + c_{p}θ̄(\frac{∂π}{∂x} + G^{13} \frac{∂π}{∂z^{∗}} ) = fv + D_{u}, (2.29)
\frac{∂v}{∂t} + u\frac{∂v}{∂x} + v \frac{∂u}{∂y} + ω \frac{∂v}{∂z^{∗}} + c_{p}θ̄(\frac{∂π}{∂y} + G^{23} \frac{∂π}{∂z^{∗}} ) = −fu + D_{v}, (2.30)
\frac{∂w}{∂t} + u\frac{∂w}{∂x} + v \frac{∂w}{∂y} + ω \frac{∂w}{∂z^{∗}} + \frac{c_{p}θ̄}{G^{1/2}} \frac{∂π}{∂z^{∗}} = g [\frac{θ}{θ̄} − 1] + D_{w}.  (2.31)

上式 (2.29)〜(2.31) 中の拡散項 D_{u}, D_{v}, D_{w} はそれぞれ (2.37), (2.38) および (2.39) である. 圧力方程式は

\frac{∂π}{∂t} + \frac{c̄^{2}}{c_{p}ρ̄θ̄^{2}} \{\frac{∂ (ρ̄θ̄u)}{∂x} + \frac{∂ (ρ̄θ̄v)}{∂y}

+G^{13} \frac{∂ (ρ̄θ̄u)}{∂z^{∗}} + G^{23} \frac{∂ (ρ̄θ̄v)}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂ (ρ̄θ̄w)}{∂z^{∗}}\}

= f_{π},  (2.32)
f_{π} = −u\frac{∂π}{∂x} − v \frac{∂π}{∂y} − ω \frac{∂π}{∂z^{∗}}

+\frac{R_{d}π}{c_{v}} \{\frac{∂u}{∂x} + \frac{∂v}{∂y} + (G^{13} + G^{23} + \frac{1}{G^{1/2}} ) \frac{∂w}{∂z^{∗}} \}

+\frac{c^{2}}{c_{p}θ^{2}} \{\frac{∂θ}{∂t} + \frac{∂θ}{∂x} + \frac{∂θ}{∂y} + (G^{13} + G^{23} + \frac{1}{G^{1/2}} ) \frac{∂θ}{∂z^{∗}} \} . (2.33)

Klemp and Wilhelmson (1978) に従って f_{π} = 0 とする. f_{π} = 0 とすることは圧力方程式を線形化し, 非断熱項 ((2.33) の θ の項) を無視することを意味する. 熱力学の式は

\frac{∂θ}{∂t} + u\frac{∂θ}{∂x} + v \frac{∂θ}{∂y} + ω \frac{∂θ}{∂z^{∗}} = D_{θ}  (2.34)

である. 上式 (2.34) 中の拡散項 D_{θ} は (2.40) 密度は, 地形に沿った座標系における乾燥大気の状態方程式から診断的に得られる:

ρ = \frac{p_{0}}{R_{d}} \frac{Π^{c_{v}/R_{d}}}{θ} .  (2.35)

\underline{第 2 章乾燥大気における準圧縮モデルの定式化と離散化 9}

運動量渦拡散係数の時間発展方程式は

\frac{∂K_{m}}{∂t} = − (u\frac{∂K_{m}}{∂x} + v \frac{∂K_{m}}{∂y} + ω \frac{∂K_{m}}{∂z^{∗}} ) − \frac{3gC_{m}^{2}l^{2}}{2θ̄G^{1/2}} \frac{∂θ_{e}}{∂z^{∗}}

+C_{m}^{2}l^{2}\{(\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} )^{2} + (\frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} )^{2} + (\frac{1}{G^{1/2}} \frac{∂w}{∂z} )^{2}\}

+\frac{C_{m}^{2}l^{2}}{2} \{(\frac{∂v}{∂x} + G^{13} \frac{∂v}{∂z^{∗}} + \frac{∂u}{∂y} + G^{23} \frac{∂u}{∂z^{∗}} )^{2}

+ (\frac{∂w}{∂x} + G^{13} \frac{∂w}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂u}{∂z})^{2}

+ (\frac{∂w}{∂y} + G^{23} \frac{∂w}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂v}{∂z} )^{2}\}

−\frac{K_{m}}{3} (\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} + \frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂w}{∂z} )

+\frac{1}{2}\{(\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) (\frac{∂K_{m}^{2}}{∂x} + G^{13} \frac{∂K_{m}^{2}}{∂z^{∗}} )

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) (\frac{∂K_{m}^{2}}{∂y} + G^{23} \frac{∂K_{m}^{2}}{∂z^{∗}} ) + \frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} \}

+ (\frac{∂K_{m}}{∂x} + G^{13} \frac{∂K_{m}}{∂z^{∗}} )^{2} + (\frac{∂K_{m}}{∂y} + G^{23} \frac{∂K_{m}}{∂z^{∗}} )^{2} + (\frac{1}{G^{1/2}} \frac{∂K_{m}}{∂z} )^{2}

−\frac{C_{ε}}{2C_{m}l^{2}} K_{m}^{2}  (2.36)

である. 地形に沿った座標系における各方程式の拡散項は以下のようになる. x, y, z^{∗} 方向の運動方程式 (2.29)〜(2.31) 中の拡散項は, それぞれ,

D_{u} = 2 (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} )\}

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂u}{∂y} + \frac{∂v}{∂x} + G^{23} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂v}{∂z^{∗}} )\}

+\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} \{K_{m} (\frac{∂w}{∂x} + \frac{1}{G^{1/2}} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂w}{∂z^{∗}} )\}

−\frac{2}{3} (\frac{∂E}{∂x} + G^{13} \frac{∂E}{∂z^{∗}} + \frac{∂E}{∂y} + G^{23} \frac{∂E}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂E}{∂z^{∗}} )  (2.37)

\underline{第 2 章乾燥大気における準圧縮モデルの定式化と離散化 10}

D_{v} = (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂v}{∂x} + \frac{∂u}{∂y} + G^{23} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂v}{∂z^{∗}} )\}

+2 (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} )\}

+\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} \{K_{m} (\frac{∂w}{∂y} + \frac{1}{G^{1/2}} \frac{∂v}{∂z^{∗}} + G^{23} \frac{∂w}{∂z^{∗}} )\}

−\frac{2}{3} (\frac{∂E}{∂x} + G^{13} \frac{∂E}{∂z^{∗}} + \frac{∂E}{∂y} + G^{23} \frac{∂E}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂E}{∂z^{∗}} )  (2.38)
D_{w} = (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂w}{∂x} + \frac{1}{G^{1/2}} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂w}{∂z^{∗}} )\}

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂w}{∂y} + \frac{1}{G^{1/2}} \frac{∂v}{∂z^{∗}} + G^{23} \frac{∂w}{∂z^{∗}} )\}

+2\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} (K_{m} \frac{1}{G^{1/2}} \frac{∂w}{∂z^{∗}} )

−\frac{2}{3} (\frac{∂E}{∂x} + G^{13} \frac{∂E}{∂z^{∗}} + \frac{∂E}{∂y} + G^{23} \frac{∂E}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂E}{∂z^{∗}} ) .  (2.39)

熱力学の式 (2.34) 中の拡散項は

D_{ϕ} = (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) (K_{h} \frac{∂θ}{∂x} + K_{h}G^{13} \frac{∂θ}{∂z^{∗}} )

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) (K_{h} \frac{∂θ}{∂y} + K_{h}G^{23} \frac{∂θ}{∂z^{∗}} )

+\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} (\frac{K_{h}}{G^{1/2}} \frac{∂θ}{∂z^{∗}} ) .  (2.40)

2.2 境界条件

下部境界 (z^{∗} = 0) と上部境界 (z^{∗} = H) では境界を横切る流れがないとする. 上下境界における境界条件はそれぞれ

ω = 0  (z^{∗} = 0),  (2.41)
ω = u\frac{∂z_{s}}{∂x} + v \frac{∂z_{s}}{∂y}  (z^{∗} = H)  (2.42)

である. 側面境界では周期境界条件を用いる.

2.3 基礎方程式系の離散化

2.3.1 離散化の概要

基礎方程式の離散化は有限差分法を用いる. 離散化された方程式系の具体的な形は付録 C.3 節に詳述した. 本節では, 本数値モデルで用いる有限差分法について簡単にまとめる.

空間微分の離散化は中心差分を用いる. 格子構造はスタッガード格子とし, 水平方向は Arakawa C グリッド, 鉛直方向は Lorenz グリッドを用いる (付録 C.1 節参照).

時間微分の離散化は, 音波に関係する項を短い時間間隔 Δτ で解き, それ以外の項を長い時間間隔 Δt で解く時間分割法を用いる (付録 C.2 節参照). 短い時間間隔で解く方程式の離散化は前進差分を, 長い時間間隔で解く方程式の離散化は中心差分を用いる. 音波に関する項のうち, 水平方向を陽解法で, 鉛直方向を陰解法で解く. この方法を水平陽解-鉛直陰解法 (HE-VI 法) という.

計算不安定を抑えるため以下の 3 点を考慮した:

• 短い時間間隔で積分する運動方程式中に音波減衰項を加える (2.3.2 節).

• 数値粘性項を圧力方程式以外の移流項に加える (2.3.3 節).

• 中央差分で長い時間間隔で積分する方程式の解に Asselin の時間フィルタ (Asselin, 1972) を適用する (2.3.4 節).

これら 3 点について以下の 3 節で詳しく述べる.

2.3.2 音波減衰項

Skamarock and Klemp (1992) は, 時間分割法を短い時間間隔で解く運動方程式に適用した場合, 不安定な音波モードが非線形項の計算によって励起し成長することがあると指摘した. 不安定なモードを除去するために, 発散に重みをかけた音波減衰項 αDiv を人工的に加える. 具体的には運動方程式 (2.29)〜(2.29) 中の π を π − αDiv に置き換える. 音波減衰項は

αDiv = α (\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} + \frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂w}{∂z^{∗}} )  (2.43)

で与えられる. Skamarock and Klemp (1992) によると, 安定な計算のためには

α ≤ 5 × 10^{−6} × \frac{Δx^{2}}{Δτ}  (2.44)

でなければならない.

2.3.3 数値粘性項

本数値モデルでは移流項を 2 次精度の中心差分で離散化する. この場合, 人工的な粘性項を加えないと計算不安定を起こすことがある. そこで以下のような数値粘性項を, 運動方程式 (2.29)〜(2.31) および熱力学の式 (2.34) の移流項に加える:

Diff.ϕ = ν_{h} [\frac{∂^{2}ϕ}{∂x^{2}} + \frac{∂^{2}ϕ}{∂y^{2}} ] + ν_{v} \frac{∂^{2}ϕ}{∂z^{∗}^{2}} .  (2.45)

ただし ϕ = u, v, w, θ とし,

ν_{h} = α_{h} \frac{ΔxΔy}{Δt} ,  (2.46)
ν_{v} = α_{v} \frac{Δz^{∗2}}{Δt}  (2.47)

である. 坪木ら (2001) によると, 安定な計算のためには

ν_{h} ≤ \frac{1}{8} \frac{ΔxΔy}{Δt} ,  (2.48)
ν_{v} ≤ \frac{1}{8} \frac{Δz^{∗2}}{Δt}  (2.49)

でなければならない.

2.3.4 Asselin の時間フィルタ

中央差分を用いて長い時間間隔で積分する場合, 微小な物理モードの擾乱で強い計算モードが生じることがある. 計算モードを減衰させるため Asselin (1972) が示した時間フィルタを用いる. ∂ϕ/∂t = F という方程式を解く場合,

ϕ^{∗t+Δt} = ϕ^{t−Δt} + 2ΔtF^{∗t}  (2.50)

と離散化される. アスタリスク (*) 付変数はフィルタがかけられていない量とする. 予報変数 ϕ^{t+Δt} を求めた後,

ϕ^{t} = ϕ^{∗t} + μ_{a} (ϕ^{∗t+Δt} − 2ϕ^{∗t} + ϕ^{t−Δt})  (2.51)

のようにフィルタリングする. 次の時刻 t + Δt における値を求めるとき, (2.51) でフィルタリングした ϕ^{t} を ϕ^{t−Δt} として解く. Asselin (1972) によると, μ_{a} の条件は

0 < μ_{a} < \frac{1}{3}  (2.52)

である. ここでは μ_{a} = 0.1 とする.

\underline{第 3 章数値実験 14}

第3章数値実験

前章で示した数理モデルを基にして 2 次元非静力対流モデルを作成した. 本章では, 音波減衰項の係数 α と数値粘性項の係数 ν_{h}, ν_{v} の具体的な値を決めるために行った数値実験の結果を示す.

3.1 数値実験の概要

作成した数値モデルは以下の点で前章で示した数理モデルと異なる:

• x, z^{∗} 方向の 2 次元.

• 地形なし. ただし数値コードにはマップファクターを入れてある.

• 運動量拡散係数は定数.

基本場における地表面圧力を 10^{5} Pa, 基本場の地表面温度を 300 K とする. その他の設定は数値実験毎に異なるので各節で示す.

ソースコードは Fortran 90 で記述した. 圧力方程式を解く際, Fujitsu Fortran & C Package ver. 3 の SSL II (Scientific Subroutine Library II) を使用した. 使用した計算機の CPU は Intell Pentium III 600 Mhz, メモリは 256 MB である.

行った数値実験は以下の 4 つである:

1. 音波の計算.

2. 音波減衰項の係数を決めるのための計算.

3. 数値粘性項の係数を決めるのための計算.

4. サーマルの計算.

3.2 数値実験

3.2.1 音波の計算

準圧縮モデルでは解に音波が含まれる. HE-VI 法で音波が正しく計算できているかを確認するため,

• 音速で伝播する波が存在する (case 1-1, case 1-2).

• 短い時間間隔に依存して鉛直方向の音波が減衰される (case 1-1 と case 1-3 の比較).

かどうかを調べる. そのための数値実験として

case 1-1 : Δτ = 1 s, 初期擾乱として, x 方向一様で, 中心 x = 1.2 × 10^{4} m, 振幅 5.0 × 10^{2} Pa, 半値幅 1.0 × 10^{3} m のガウス型の圧力を与える.

case 1-2 : Δτ = 1 s, 初期擾乱として, z^{∗} 方向一様で, 中心 z^{∗} = 5.0 × 10^{4} m, 振幅 5.0 × 10^{2} Pa, 半値幅\sqrt{5.0} × 10^{3} m のガウス型の圧力を与える.

case 1-3 : Δτ = 0.1 s, 初期擾乱として, x 方向一様で, 中心 x = 1.2 × 10^{4} m, 振幅 5.0 × 10^{2} Pa, 半値幅 1.0 × 10^{3} m のガウス型の圧力を与える (case 1-1 と同じ).

を行う. case 1-1 〜 1-3 の数値実験を行う際の設定を表 3.1 に示す. ここでは音波のみを見るため, 重力加速度を 0 とし, 移流項, コリオリ項, 浮力項, 拡散項および音波減衰項は全て 0 とする.

表 3.1: case 1-1〜1-3 の数値実験における設定.

水平計算領域  : 2.4 ×10^{4} m
鉛直計算領域  : 1.0 ×10^{4} m
水平格子間隔 Δx : 1.0 ×10^{3} m
鉛直格子間隔 Δz^{∗} : 5.0×10^{2} m (z^{∗} 方向に等間隔)
短い時間間隔 Δτ : 1.0 s または 1.0 ×10^{−1} s
計算時間 t  : 1.0 ×10^{2} s
初期擾乱  : x 方向一様で, 中心 x = 1.2 × 10^{4} m, 振幅 5.0 × 10^{2} Pa,

半値幅 1.0 × 10^{3} m のガウス型の圧力

または

z^{∗} 方向一様で, 中心 z^{∗} = 5.0 × 10^{4} m, 振幅 5.0 × 10^{2} Pa, 半値幅\sqrt{5.0} × 10^{3} m のガウス型の圧力

基本場  : 等温 (300 K) , 等圧大気

図 3.1〜図 3.4 の実験結果から, 鉛直速度, 水平速度およびエクスナー関数の偏差がほぼ音速で伝播していることがわかる. 与えた音速は

c̄ = (\frac{c_{p}R_{d}}{c_{v}} Π̄θ̄)^{1/2}
∼ 339 m/s  (3.1)

である. ただし c_{p} = 1.952 × 10^{3} J K^{−1} kg^{−1}, c_{v} = 1.463 × 10^{3} J K^{−1} kg^{−1}, R_{d} = 2.87 × 10^{2} J K^{−1} kg^{−1} とした. 図 3.1 と図 3.2 では, それぞれ鉛直速度の偏差とエクスナー関数の偏差が t = 15 s 後に z^{∗} = 5.0 × 10^{3} m に達していることがわかる. これより見積もった音速の大きさは約 333 m/s である. 図 3.3 と図 3.4 では, それぞれ鉛直速度とエクスナー関数の偏差が共に t = 35 s 後に z^{∗} = 1.2 × 10^{4} m に達しているので, 音速の大きさは約 343 m/s である.

図 3.5〜図 3.8 の実験結果から, 鉛直方向の音波が減衰されることが確認できる. 短い時間間隔が 1 s である図 3.5 と図 3.6 の波の振幅は, 始めの 1 周期に比べて約 1/2 になっている. これに対し, 短い時間間隔を 1/10 にした図 3.7 と図 3.8 の波の振幅は, ほとんど減衰が見られない.

図 3.1: case 1-1 における音波の計算結果 . 鉛直速度の偏差の x = 1.2 × 10^{4} m における時間- z^{∗} 軸断面図.

図 3.2: case 1-1 における音波の計算結果 . エクスナー関数の偏差の x = 1.2 × 10^{4} m における時間- z^{∗} 軸断面図.

図 3.3: case 1-2 における音波の計算結果 . 鉛直速度の偏差の z^{∗} = 5.0 × 10^{3} m における時間- x 軸断面図.

図 3.4: case 1-2 における音波の計算結果 . エクスナー関数の偏差の z^{∗} = 5.0 × 10^{3} m における時間- x 軸断面図.

図 3.5: case 1-1 における音波の計算結果. (x, z^{∗}) = (2.0 × 10^{3}, 5.0 × 10^{3}) における鉛直速度の偏差の時間変化.

図 3.6: case 1-1 における音波の計算結果. (x, z^{∗}) = (2.0 × 10^{3}, 5.0 × 10^{3}) におけるエクスナー関数の偏差の時間変化.

図 3.7: case 1-3 における音波の計算結果. (x, z^{∗}) = (2.0 × 10^{3}, 5.0 × 10^{3}) における鉛直速度の偏差の時間変化.

図 3.8: case 1-3 における音波の計算結果. (x, z^{∗}) = (2.0 × 10^{3}, 5.0 × 10^{3}) におけるエクスナー関数の偏差の時間変化.

3.2.2 音波減衰項の係数を決めるための計算

音波減衰項の係数 α の条件は (2.44) で与えられる. 表 3.2 の示す設定においては Δx = 1.0 × 10^{3} m, Δτ = 1.0 s なので, 音波減衰項の係数の条件は

α ≤ 2.5  (3.2)

となる. 具体的な α を決めるために以下のように値を代入して数値計算を行う:

case 2-1 : α = 0 (音波減衰項がない場合).

case 2-2 : α = 2.5 × 10^{−2}.

case 2-1 と 2-2 の数値実験を行う際の設定を表 3.2 に示す. ここでは音波のみを見るため, 重力加速度を 0 とし, 移流項, コリオリ項, 浮力項, 拡散項および音波減衰項は全て 0 とする.

表 3.2: case 2-1〜2-2 の数値実験における設定.

水平計算領域  : 2.4 ×10^{4} m
鉛直計算領域  : 1.0 ×10^{4} m
水平格子間隔 Δx : 1.0 ×10^{3} m
鉛直格子間隔 Δz^{∗} : 5.0 ×10^{2} m (z^{∗} 方向に等間隔)
短い時間間隔 Δτ : 1 s
計算時間 t  : 1.0 ×10^{4}s
初期擾乱  : 中心 (x, z^{∗}) = (1.2 × 10^{4}, 5.0 × 10^{3}), 振幅 5.0 × 10^{2} Pa,

半振幅 x = \sqrt{5.0} × 10^{3}, z^{∗} = 1.0 × 10^{3} の 2 次元ガウス型の圧力 (図 3.9)

基本場  : 等温 (300 K) 大気, 静水圧平衡

図3.10〜図3.13 の実験結果から, 音波減衰項の係数として適当な値はα = 2.5×10^{−2} であることがわかる. 音波減衰項を加えない図 3.11 では t = s あたりでエクスナー関数が急増している. 一方, 図 3.13 から t = 9100 s あたりでも値が安定している.

図 3.9: 初期擾乱として与えた圧力をエクスナー関数に変換したもの. 与えた初期擾乱は, 中心 (x, z^{∗}) = (1.2 × 10^{4}, 5.0 × 10^{3}), 振幅 5.0 × 10^{2} Pa, 半振幅 x = \sqrt{5.0} × 10^{3}, z^{∗} = 1.0 × 10^{3} の 2 次元ガウス型の圧力である.

3.2.3 数値粘性項の係数を決めるための計算

数値粘性項の係数 ν_{h}, ν_{v} の条件は (2.48) と (2.49) である. 表 3.3 の示す設定においては

ν_{h} ≤ 1.25 × 10^{5},  (3.3)
ν_{v} ≤ 3.125 × 10^{4}  (3.4)

である. 具体的な ν_{h} と ν_{v} の値を決めるため, 以下のように値を代入して温位の移流を数値計算する:

case 3-1 : ν_{h} = 0, ν_{v} = 0 (数値粘性項がない場合).
case 3-2 : ν_{h} = 5.0 × 10^{3}, ν_{v} = 1.25 × 10^{3} .
case 3-3 : ν_{h} = 5.0 × 10^{2}, ν_{v} = 1.25 × 10^{2} .
case 3-4 : ν_{h} = 5.0, ν_{v} = 1.25 .

case 3-1 〜 3-4 の数値実験を行う際の設定を表 3.3 に示す. ここでの数値モデルはコリオリ項, 浮力項および拡散項を 0 とする.

図 3.14〜図 3.17 の結果から, 数値粘性項の係数として ν_{h} = 5.0, ν_{v} = 1.25 が適当であることがわかる. 数値粘性項を追加した図 3.15〜図 3.17 のうち, もっとも安定なものは図 3.17 である. 図 3.15 では t = 500 あたりで数値計算によって出てきた偏差が大きくなっている. 図 3.16 では t = 750 あたりで温位偏差の後ろに擾乱が出てきている. 図 3.16 と図 3.17 を比べると, 温位偏差の後ろに出ている擾乱が小さいのは図 3.17 である.

図 3.10: case 2-1 における音波の計算結果. エクスナー関数の偏差について, 時刻 t = 0 〜 250 s 間, z^{∗} = 5.0 × 10^{3} における時間- x 軸断面図.

図 3.11: case 2-1 における音波の計算結果. エクスナー関数の偏差について, 時刻 t = 9000 〜 9250 s 間, z^{∗} = 5.0 × 10^{3} における時間- x 軸断面図.

図 3.12: case 2-2 における音波の計算結果. エクスナー関数の偏差について, 時刻 t = 0 〜 250 s 間, z^{∗} = 5.0 × 10^{3} における時間- x 軸断面図.

図 3.13: case 2-2 における音波の計算結果. エクスナー関数の偏差について, 時刻 t = 9000 〜 9250 s 間, z^{∗} = 5.0 × 10^{3} における時間- x 軸断面図.

表 3.3: case 3-1〜3-2 の数値実験における設定.

水平計算領域  : 2.4 ×10^{4} m
鉛直計算領域  : 1.0 ×10^{4} m
水平格子間隔 Δx : 1.0 ×10^{3} m
鉛直格子間隔 Δz^{∗} : 5.0×10^{2} m (z^{∗} 方向に等間隔)
長い時間間隔 Δt : 1 s
計算時間 t  : 10,000 s
初期擾乱  : 中心 (x, z^{∗}) = (5.0 × 10^{3}, 5.0 × 10^{3}), 振幅 5.0 × 10^{2} Pa,

半振幅 x = \sqrt{5.0} × 10^{3}, z^{∗} = \sqrt{5.0} × 10^{3} の 2 次元ガウス型の温位

水平速度  10 m/s (一定値)
鉛直速度  0 m/s (一定値)

図 3.14: case 3-1 における温位の計算結果. 左から順に時刻 t = 0, t = 250, t = 500, t = 750, t = 1000 s における温位の偏差の時間-z^{∗} 軸断面図.

\underline{第 3 章数値実験 26}

図 3.15: case 3-2 における温位の計算結果. 左から順に時刻 t = 0, t = 250, t = 500 s における温位の偏差の時間-z^{∗} 軸断面図. 図 3.16: case 3-3 における温位の計算結果. 左から順に時刻 t = 0, t = 250, t = 500, t = 750, t = 1000 s における温位の偏差の時間-z^{∗} 軸断面図.

\underline{第 3 章数値実験 27}

{図 3.17: case 3-4 における温位の計算結果. 左から順に時刻 t = 0, t = 250, t = \atop 500, t = 750, t = 1000 s における温位の偏差の時間-z^{∗} 軸断面図.}

3.2.4 サーマルの計算

前節までに得られた安定計算に必要な音波減衰項と数値粘性項の係数を用いて, サーマルの数値実験を行う. サーマルの計算によって, 簡単な対流が再現できるかを確認する. サーマルの計算に用いた数値モデルには拡散項が加えられている. 運動量渦拡散係数として一定値 (10, 50, 100 m^{2}/s) を用いる:

case 4-1 : K_{m} = 10 m^{2}/s
case 4-2 : K_{m} = 50 m^{2}/s
case 4-3 : K_{m} = 100 m^{2}/s

case 4-1 〜 4-3 の数値実験を行う際の設定を表 3.4 に示す. ここでの数値モデルはコリオリ項のみ 0 とする.

表 3.4: case 4-1〜4-3 の数値実験における設定.

水平計算領域  : 2.4 ×10^{4} m
鉛直計算領域  : 1.0 ×10^{4} m
水平格子間隔 Δx  : 1.0 ×10^{3} m
鉛直格子間隔 Δz^{∗}  : 5.0×10^{2} m (z^{∗} 方向に等間隔)
短い時間間隔 Δτ  : 2 s
長い時間間隔 Δt  : 10 s
計算時間 t  : 2.0 ×10^{3} s
初期擾乱  : 中心 (x, z^{∗}) = (1.2 × 10^{4}, 0.0), 振幅 5.0 × 10^{2} m,

半振幅 x = \sqrt{5.0} × 10^{3}, z^{∗} = \sqrt{5.0} × 10^{3} の 2 次元ガウス型の温度

基本場  : z^{∗} = 5000 m まで乾燥断熱減率で温度減少,

z^{∗} = 5000 m より上は z^{∗} = 5000 m での温度で等温, 静水圧平衡

音波減衰項の係数 α  : 2.5 × 10^{−2}
水平数値粘性項の係数 ν_{h} : 5.0
鉛直数値粘性項の係数 ν_{v} : 1.25

図 3.18〜図 3.32 の結果から, 上昇するサーマルを表現できたことが確認できる. case 4-1 〜 4-3 のどれも t = 500 s あたりからサーマルが上昇し始めるのが確認できる. t = 800 s あたりでサーマルの先端が z^{∗} = 5000 m に達する. t = 800 s あたりでサーマルが潰れて横に広がっているのがわかる. これはサーマルが上昇するにつれてサーマル内の温度が下がり, 高度 z^{∗} = 5.0 × 10^{3} m あたりで周りの温度に等しくなるため, これ以上上昇しなくなる. 上昇しなくなったサーマルは横に広がっていく.

サーマルの数値実験において, t = 700 s あたり (図 3.20, 図 3.25 および図 3.30) でサーマルの下から擾乱が生じている. t = 1000 s である図 3.22, 図 3.27 および図 3.32 では, さらに大きくなっている.

図 3.18: case 4-1 における時刻 t = 0 s での温位の偏差の計算結果.

\underline{第 3 章数値実験 30}

図 3.19: case 4-1 における時刻 t = 500 s での温位の偏差の計算結果. 図 3.20: case 4-1 における時刻 t = 700 s での温位の偏差の計算結果.

\underline{第 3 章数値実験 31}

図 3.21: case 4-1 における時刻 t = 800 s での温位の偏差の計算結果. 図 3.22: case 4-1 における時刻 t = 1000 s での温位の偏差の計算結果.

\underline{第 3 章数値実験 32}

図 3.23: case 4-2 における時刻 t = 0 s での温位の偏差の計算結果. 図 3.24: case 4-2 における時刻 t = 500 s での温位の偏差の計算結果.

\underline{第 3 章数値実験 33}

図 3.25: case 4-2 における時刻 t = 700 s での温位の偏差の計算結果. 図 3.26: case 4-2 における時刻 t = 800 s での温位の偏差の計算結果.

\underline{第 3 章数値実験 34}

図 3.27: case 4-2 における時刻 t = 1000 s での温位の偏差の計算結果. 図 3.28: case 4-3 における時刻 t = 0 s での温位の偏差の計算結果.

\underline{第 3 章数値実験 35}

図 3.29: case 4-3 における時刻 t = 500 s での温位の偏差の計算結果. 図 3.30: case 4-3 における時刻 t = 700 s での温位の偏差の計算結果.

\underline{第 3 章数値実験 36}

図 3.31: case 4-3 における時刻 t = 800 s での温位の偏差の計算結果. 図 3.32: case 4-3 における時刻 t = 1000 s での温位の偏差の計算結果.

\underline{第 4 章まとめと考察 37}

第4章まとめと考察

本研究では, 地形に沿った座標系における非静力学モデルを作成するための基礎となるべき定式化と離散化を示した. 直交直線座標系における定式化は Klemp and Wilhelmson (1978) に従い, 地形に沿った座標系は斉藤 (1999) に従って定式化を示した. 離散化では計算不安定を抑えるため, 音波減衰項と数値粘性項を人工的に加え, Asselin (1972) の時間フィルターを用いた.

定式化の際, Klemp and Wilhelmson (1978) の数理モデルには消去されている項があることを見つけた. 本論文中のサブグリッドスケールの乱流運動エネルギー方程式中の項 (2.19) と (2.20) が, Klemp and Wilhelmson (1978) では離散化の段階で消去されている. また消去されたことに関する議論はない.

定式化と離散化を行った後, 乾燥大気について 2 次元の非静力学対流モデルを作成した. 作成した数値モデルは音波と簡単なサーマルの計算ができる. 作成した数値モデルを用いて数値実験を行った結果, 以下のことがわかった:

• Δx = 1, 000 m, Δτ = 1.0 s の場合, 音波減衰項の係数は α = 2.5 × 10^{−2} が適当である.

• Δx = 1, 000 m, Δτ = 1.0 s の場合, 数値粘性項の係数は ν_{h} = 5.0, ν_{v} = 1.25 が適当である.

ここで, Klemp and Wilhelmson (1978) が項 (2.19) と (2.20) を消去した理由について考察する. まず, 項 (2.19) について考える. 乱流運動エネルギー方程式はサブグリッドスケールでの運動を考えている. サブグリッドスケールでは密度の空間的時間的変化が小さいと考えると, 非圧縮性流体の連続の式を用いることができる. したがって, 項 (2.19) は連続の式によって消去される.

次に, 項 (2.20) に関してはオーダを見る. 代表的な長さスケールを L, 代表的な運動量渦拡散係数の大きさを K とする. 項 (2.20) のオーダは

δ_{ij} (\frac{∂K_{m}}{∂x_{j}} )^{2} ∼ O (\frac{K^{2}}{L^{2}} )  (4.1)

である. ここで, 式 (2.16) 中で Klemp and Wilhelmson (1978) には残されている第 5 項について, オーダを見積もると

\frac{1}{2} \frac{∂^{2}K_{m}^{2}}{∂x^{2}} ∼ O (\frac{K^{2}}{L^{2}} ) .  (4.2)

したがって, 項 (2.20) と残されている第 6 項 (4.2) のオーダは同じである. オーダを見ると項 (2.20) を消去すべきではない.

謝辞

本研究を進めるにあたって多くの指導と助言をして頂いた方々に深く感謝します.

指導教官である林祥介教授には本研究を始めるきっかけを頂き, 研究を進める上で多くの助言を頂きました. 九州大学理学部地球惑星科学科の中島健介助手には, 地理的距離があるにもかかわらず, 定式化から離散化, プログラミングに至るまで様々な指導と助言を頂きました. 北海道大学大学院地球環境科学研究科の石渡正樹助手にはプログラミングや修論の内容等に関して様々な助言を頂きました. 地球流体力学研究室の小高正嗣助手と同研究室博士課程の杉山耕一朗さんにはなかなか進まない研究に毎週のようにつき合って頂き, 多くの指導や助言をして頂きました. 同研究室博士課程の小西丈予さんには定式化やプログラミングで詰まったとき等に助言や励ましを頂きました. 九州大学大学院理学府地球惑星科学専攻地球流体力学研究分野博士課程の中野満寿男さんには修士論文を参考にさせて頂いただけでなく, 定式化等に関して助言も頂きました. 本研究の大きな支えとなった CReSS を開発された名古屋大学地球水循環研究センターの坪木和久さんと財団法人高度情報科学技術研究機構の榊原篤志さんには紙面をかりて御礼申し上げます.

北海道大学低温科学研究所寒冷海洋圏科学部門雲科学分野の方々や名古屋大学地球水循環研究センター局域水循環過程研究部門第 1 研究室の院生の方々には, 本研究を始めるにあたり多くの助言を頂きました. 地球流体力学研究室・惑星物理学研究室の方々には研究に関する議論だけでなく, 多くの励ましも頂きました. 本当にありがとうございました.

本研究の数値モデルを作成するにあたり, コンパイラには Fujitsu Fortran & C Package ver.3 を使わせて頂きました. また, 作図には地球流体電脳倶楽部の gtool4 Fortran90 ツール・ライブラリ (http://www.gfd-dennou.org/arch/gtool4/) を利用させて頂きました.

\underline{付録 A 乾燥大気における準圧縮モデルの定式化の補足 40}

付録A {乾燥大気における準圧縮モ \atop デルの定式化の補足}

本章では主に第 2 章の定式化で解説していないことについての詳細を述べる. 予報変数, 基本場の設定や記号などは本文中と同じである.

A.1 地形に沿った座標系

ここでは基礎方程式系を地形に沿った座標系に変換するために, 座標系, 速度および空間微分の定義をする.

斉藤 (1999) や CReSS (坪木ら, 2001) と同様に, 地形に沿った座標 (x^{∗}, y^{∗}, z^{∗}) を

x^{∗} = x,
y^{∗} = y,
z^{∗} = z^{∗}(x, y, z)

とする. z^{∗} は (2.21) で定義される (図 2.1).

地形に沿った座標系での速度を (u^{∗}, v^{∗}, ω) とすると, 直交直線座標系での速度は

u = u^{∗} \frac{∂x}{∂x^{∗}} + v^{∗} \frac{∂x}{∂y^{∗}} + ω \frac{∂x}{∂z^{∗}} ,  (A.1)
v = u^{∗} \frac{∂y}{∂x^{∗}} + v^{∗} \frac{∂y}{∂y^{∗}} + ω \frac{∂y}{∂z^{∗}} ,  (A.2)
w = u^{∗} \frac{∂z}{∂x^{∗}} + v^{∗} \frac{∂z}{∂y^{∗}} + ω \frac{∂z}{∂z^{∗}}  (A.3)

と書ける. 上式 (A.1)-(A.3) を u^{∗}, v^{∗}, ω について解くと,

G^{1/2}u^{∗} = uJ_{y^{∗}z^{∗}}^{yz} + vJ_{y^{∗}z^{∗}}^{zx} + wJ_{y^{∗}z^{∗}}^{xy},  (A.4)

G^{1/2}v^{∗} = uJ_{z^{∗}x^{∗}}^{yz} + vJ_{z^{∗}x^{∗}}^{zx} + wJ_{z^{∗}x^{∗}}^{xy},  (A.5)
G^{1/2}ω = uJ_{x^{∗}y^{∗}}^{yz} + vJ_{x^{∗}y^{∗}}^{zx} + wJ_{x^{∗}y^{∗}}^{xy}  (A.6)

を得る. ただし J は 2 次元のヤコビアンで, 例えば

J_{y^{∗}z^{∗}}^{yz} ≡ \frac{∂(y, z)}{∂(y^{∗}, z^{∗})}

= | \begin{matrix}\frac{∂y}{∂y^{∗}} & \frac{∂y}{∂z^{∗}} \\ \frac{∂z}{∂y^{∗}} & \frac{∂z}{∂z^{∗}} \end{matrix} | ,

G は 3 次元のヤコビアンで

G ≡ \frac{∂(x, y, z)}{∂(x^{∗}, y^{∗}, z^{∗})}

= | \begin{matrix}\frac{∂x}{∂x^{∗}} & \frac{∂x}{∂y^{∗}} & \frac{∂x}{∂z^{∗}} \\ \frac{∂y}{∂x^{∗}} & \frac{∂y}{∂y^{∗}} & \frac{∂y}{∂z^{∗}} \\ \frac{∂z}{∂x^{∗}} & \frac{∂z}{∂y^{∗}} & \frac{∂z}{∂z^{∗}} \end{matrix} |

と定義される. したがって (2.21) を代入して整理すると, 地形に沿った座標系での速度は

U = u,  (A.7)
V = v,  (A.8)
W = G^{13}u + G^{23}v + \frac{1}{G^{1/2}} w  (A.9)

と書ける.

直交直線座標系から地形に沿う座標系へ変換する場合, ある変数 ϕ の空間微分は

\frac{∂ϕ}{∂x} = \frac{∂ϕ}{∂x} + G^{13} \frac{∂ϕ}{∂z^{∗}} ,  (A.10)
\frac{∂ϕ}{∂y} = \frac{∂ϕ}{∂y} + G^{23} \frac{∂ϕ}{∂z^{∗}} ,  (A.11)
\frac{∂ϕ}{∂z} = \frac{1}{G^{1/2}} \frac{∂ϕ}{∂z^{∗}}  (A.12)

となる.

A.2 サブグリッドスケールの拡散

大気中の運動には様々なスケールがある. これらの運動には, 数値モデルで考える格子点で表現できるグリッドスケールの平均運動と, 表現できない乱流運動 (またはサブグリッドスケールの運動) がある. 乱流状態にある場合, 同時に存在する様々な時空間スケールの運動が相互作用しているため, サブグリッドスケールの運動も無視することはできない. これらの運動を分離するために, 予報変数 A(= u, v, w, π, θ, q_{v}, q_{c}, q_{r}) をグリッドスケールで表現できる平均成分 (\overline{A}) とそこからの変動成分 (A′′) に分ける:

A = \overline{A} + A′′.

平均の方法としてレイノルズ平均を用いる. このとき平均成分と変動成分は

\overline{A′′} = 0,
\overline{AB} = \overline{A} \overline{B} + \overline{A′′}\overline{B′′},
\overline{A′′\overline{B}} = 0

を満たす.

直交直線座標系における運動方程式, 連続の式に対してレイノルズ平均操作を行うと,

\frac{∂\overline{u}_{i}}{∂t} + \overline{u}_{j} \frac{∂\overline{u}_{i}}{∂x_{j}} + c_{p}Θ̄_{v} \frac{∂\overline{π}}{∂x_{i}} = δ_{i3}g [\frac{θ}{Θ̄} − 1]
−ε_{i3k}fu_{k} − \frac{∂}{∂x_{j}} \overline{u′′_{i}u′′_{j}},  (A.13)
\frac{∂\overline{ρ}}{∂t} + \frac{∂\overline{ρu}_{i}}{∂x_{i}} = 0  (A.14)

を得る. ただし上付バーの混同を避けるため本節のみ, 基本場での θ を Θ̄ と表す. 上式 (A.13) 中の右辺第 3 項 −\frac{∂}{∂x_{j}} \overline{u′′_{i}u′′_{j}} が式 (2.5) 中の D_{u_{i}} である. Klemp and

Wilhelmson (1978) では, 連続の式を変形した圧力方程式にサブグリッドスケールの拡散項 D_{π} がある.

熱力学の式に対して同様な操作を行うと,

\frac{∂\overline{θ}}{∂t} + \overline{u_{j}} \frac{∂\overline{θ}}{∂x_{j}} = −\frac{∂}{∂x_{j}} \overline{u′′_{j}θ′′}.  (A.15)

上式 (A.15) 中の −\frac{∂}{∂x_{j}} \overline{u′′_{j}θ′′} が, それぞれ式 (2.9) の拡散項 D_{θ} である.

(A.13) の右辺第 3 項に\overline{u′′_{i}u′′_{j}}, (A.15) の右辺に\overline{u′′_{j}ϕ′′} があるために, \overline{u}_{i}, \overline{π}, \overline{θ}, \overline{q}_{v}, \overline{q}_{c}, \overline{q}_{r} に関して方程式系が閉じていない. 方程式を閉じるために, Klemp and Wilhelmson (1978) では乱流運動エネルギーを用いた 1.5 次のクロージャモデルを使っている. 本数理モデルは Klemp and Wilhelmson (1978) に従う. このとき\overline{u′′_{i}u′′_{j}} と\overline{u′′_{j}θ′′} を以下のように仮定する:

\overline{u′′_{i}u′′_{j}} = −K_{m} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} ) + \frac{2}{3}δ_{ij}E,  (A.16)
\overline{u′′_{j}θ′′} = K_{h} \frac{∂θ}{∂x_{j}} .  (A.17)

E, K_{m} および K_{h} は本文中と同様にして求める.

\underline{付録 B 湿潤大気における準圧縮モデルの定式化 44}

付録B {湿潤大気における準圧縮モ \atop デルの定式化}

本章では湿潤大気における 3 次元準圧縮モデルの定式化を示す.

基礎方程式系は, 3 方向の運動方程式と熱力学の式に加え, 水蒸気混合比の式, 雲水混合比の式および雨水混合比の式である. x, y 方向の運動方程式は乾燥空気の場合と同じある. z または z^{∗} 方向の運動方程式, 熱力学の式およびサブグリッドスケールの乱流運動エネルギー方程式は蒸発・凝結の効果を考えるため, 乾燥空気の方程式と形が異なる.

B.1 直交直線座標系における基礎方程式系

本節では直交直線座標系における基礎方程式系を導出する.

B.1.1 変数の定義

湿潤大気における予報変数は

\begin{aligned}u & : & 速度の x 成分 \\ v & : & 速度の y 成分 \\ w & : & 速度の z 成分 \\ Π & : & 無次元化した圧力 (エクスナー関数) \\ θ & : & 温位 \\ q_{v} & : & 水蒸気の混合比 \\ q_{c} & : & 雲水の混合比 \\ q_{r} & : & 雨水の混合比 \end{aligned}

であり, 空間の変数 x, y, z と時間の変数 t である.

密度は湿潤空気の状態方程式

p = ρR_{d}T (1 + 0.61q_{v})  (B.1)

または

Π = (\frac{R_{d}}{p_{0}} ρθ_{v})^{R_{d}/c_{v}}  (B.2)

から診断的に与えられる. ただし本章での ρ は湿潤空気の密度とする. θ_{v} は仮温位で

θ_{v} = θ (1 + 0.61q_{v})  (B.3)

と表される.

B.1.2 基本場の定義と分離

基本場も乾燥大気と同様の状態を考える. 水平一様 (z のみの関数) な静止 (ū = v̄ = w̄ = 0) 状態とする. 基本場における無次元量の圧力と温位は, 乾燥大気と同様であり, 仮温位は

θ̄_{v} = θ̄(1 + 0.61q̄_{v})  (B.4)

である. 湿潤大気における基本場の静水圧平衡の式は

\frac{∂Π̄}{∂z} = −\frac{g}{c_{p}θ̄_{v}}  (B.5)

である.

乾燥大気の場合と同様に, 変数を基本場とそこからの偏差に分ける:

ϕ = ϕ̄ + ϕ′.

ただし ϕ = u, v, w, Π, θ, q_{v}, q_{v}, q_{r}, θ_{v}, ρ を表す. 以下では, 記述の簡便のため Π′ = π とし, 速度のみ上付きプライム (′) を省略する.

B.1.3 基礎方程式系

湿潤大気の場合, 直交直線座標系における基礎方程式は, 運動方程式, 熱力学の式, 準圧縮系の連続の式, 水蒸気混合比の式, 雲水混合比の式および雨水混合比の式である. ただし連続の式の代わりに圧力方程式を用いる.

運動方程式は

\frac{du_{i}}{dt} + c_{p}θ̄_{v} \frac{∂π}{∂x_{i}}
= δ_{i3}g [\frac{θ}{θ̄} − 1 + 0.61 (q_{v} − q̄_{v}) − q_{c} − q_{r}] − ε_{i3k}fu_{k} + D_{u_{i}} (B.6)

である. i = 1, 2, 3 で, それぞれ x, y および z 成分を表す. f はコリオリパラメータである. D_{u_{i}} はサブグリッドスケールの乱流に伴う拡散項であり B.3 節で述べる.

圧力方程式は

\frac{∂π}{∂t} + \frac{c̄^{2}}{c_{p}ρ̄θ̄_{v}^{2}} \frac{∂}{∂x_{j}} (ρ̄θ̄_{v}u_{j}) = f_{π},  (B.7)
f_{π} = −u_{j} \frac{∂π}{∂x_{j}} + \frac{R_{d}π}{c_{v}} \frac{∂u_{j}}{∂x_{j}} + \frac{c^{2}}{c_{p}θ_{v}^{2}} \frac{dθ_{v}}{dt}  (B.8)

である. D_{π} はサブグリッドスケールの乱流に伴う拡散項, c は音波で

c^{2} = \frac{c_{p}R_{d}}{c_{v}} Πθ_{v}  (B.9)

と表され, c̄ は基本場での音波である. 圧力方程式 (B.8) はエクスナー関数を微分したものと圧縮性流体の連続の式

\frac{∂ρ}{∂t} + \frac{∂}{∂x_{j}} (ρu_{j}) = 0  (B.10)

から得られる.

熱力学の式, 水蒸気混合比の式, 雲水混合比の式および雨水混合比の式は, それぞれ,

\frac{dθ}{dt} = D_{θ} + M_{θ},  (B.11)
\frac{dq_{v}}{dt} = D_{q_{v}} + M_{q_{v}},  (B.12)
\frac{dq_{c}}{dt} = D_{q_{c}} + M_{q_{c}},  (B.13)
\frac{dq_{r}}{dt} = D_{q_{r}} + M_{q_{r}}  (B.14)

である. それぞれの方程式中の M_{θ}, M_{q_{v}}, M_{q_{c}}, M_{q_{v}} は微物理過程に伴う生成・消滅項である. D_{θ}, D_{q_{v}}, D_{q_{c}}, D_{q_{v}} はサブグリッドスケールの乱流に伴う拡散項である. 生成・消滅項については B.4 節で, 拡散項については B.3 節で述べる.

密度は湿潤空気の状態方程式

ρ = \frac{p_{0}}{R_{d}} \frac{Π^{c_{v}/R_{d}}}{θ_{v}}  (B.15)

より与えられる.

B.2 地形に沿った座標系における基礎方程式系

本節では, 前節で導出した直交直線座標系での方程式系を地形に沿う座標系に変換する. 座標系の設定は乾燥大気の場合と同様である.

B.2.1 基礎方程式系

静水圧平衡の式および前節で導出した基礎方程式系 (B.6), (B.8), (B.11)〜(B.14) を, 直交直線座標系から地形に沿った座標系に変換する.

静水圧平衡の式は

\frac{∂Π̄}{∂z^{∗}} = −\frac{G^{1/2}g}{c_{p}θ̄_{v}}  (B.16)

となる.

運動方程式は

\frac{∂u}{∂t} + u\frac{∂u}{∂x} + v \frac{∂u}{∂y} + ω \frac{∂u}{∂z^{∗}} + c_{p}θ̄_{v} (\frac{∂π}{∂x} + G^{13} \frac{∂π}{∂z^{∗}} ) = fv + D_{u}, (B.17)
\frac{∂v}{∂t} + u\frac{∂v}{∂x} + v \frac{∂u}{∂y} + ω \frac{∂v}{∂z^{∗}} + c_{p}θ̄_{v} (\frac{∂π}{∂y} + G^{23} \frac{∂π}{∂z^{∗}} ) = −fu + D_{v}, (B.18)
\frac{∂w}{∂t} + u\frac{∂w}{∂x} + v \frac{∂w}{∂y} + ω \frac{∂w}{∂z^{∗}} + \frac{c_{p}θ̄_{v}}{G^{1/2}} \frac{∂π}{∂z^{∗}}
= g [\frac{θ}{θ̄} − 1 + 0.61 (q_{v} − q̄_{v}) − q_{c} − q_{r}] + D_{w} (B.19)

\underline{付録 B 湿潤大気における準圧縮モデルの定式化 48}

となる. 拡散項 D_{u}, D_{v}, D_{w} は B.3 節で述べる. 圧力方程式は

\frac{∂π}{∂t} + \frac{c̄^{2}}{c_{p}ρ̄θ̄_{v}^{2}} \{\frac{∂ (ρ̄θ̄_{v}u)}{∂x} + \frac{∂ (ρ̄θ̄_{v}v)}{∂y}

+G^{13} \frac{∂ (ρ̄θ̄_{v}u)}{∂z^{∗}} + G^{23} \frac{∂ (ρ̄θ̄_{v}v)}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂ (ρ̄θ̄_{v}w)}{∂z^{∗}}\}

= f_{π},  (B.20)
f_{π} = −u\frac{∂π}{∂x} − v \frac{∂π}{∂y} − ω \frac{∂π}{∂z^{∗}}

+\frac{R_{d}π}{c_{v}} \{\frac{∂u}{∂x} + \frac{∂v}{∂y} + (G^{13} + G^{23} + \frac{1}{G^{1/2}} ) \frac{∂w}{∂z^{∗}} \}

+\frac{c^{2}}{c_{p}θ_{v}^{2}} \{\frac{∂θ_{v}}{∂t} + \frac{∂θ_{v}}{∂x} + \frac{∂θ_{v}}{∂y} + (G^{13} + G^{23} + \frac{1}{G^{1/2}} ) \frac{∂θ_{v}}{∂z^{∗}} \} . (B.21)

Klemp and Wilhelmson (1978) に従い f_{π} = 0 とすると, 準圧縮モデルの圧力方程式になる. f_{π} = 0 とすることは線形化を行って非断熱項 ((B.21) の右辺第 3 項) を無視することを意味する. 熱力学の式, 水蒸気混合比の式, 雲水混合比の式および雨水混合比の式は, それぞれ,

\frac{∂θ}{∂t} + u\frac{∂θ}{∂x} + v \frac{∂θ}{∂y} + ω \frac{∂θ}{∂z^{∗}} = M_{θ} + D_{θ},  (B.22)
\frac{∂q_{v}}{∂t} + u\frac{∂q_{v}}{∂x} + v \frac{∂q_{v}}{∂y} + ω \frac{∂q_{v}}{∂z^{∗}} = M_{q_{v}} + D_{q_{v}},  (B.23)
\frac{∂q_{c}}{∂t} + u\frac{∂q_{c}}{∂x} + v \frac{∂q_{c}}{∂y} + ω \frac{∂q_{c}}{∂z^{∗}} = M_{q_{c}} + D_{q_{c}},  (B.24)
\frac{∂q_{r}}{∂t} + u\frac{∂q_{r}}{∂x} + v \frac{∂q_{r}}{∂y} + ω \frac{∂q_{r}}{∂z^{∗}} = M_{q_{r}} + D_{q_{r}}  (B.25)

と書ける. それぞれの方程式中の乱流項 D_{θ}, D_{q_{v}}, D_{q_{c}}, D_{q_{v}} は B.3 節で, 微物理項 M_{θ}, M_{q_{v}}, M_{q_{c}}, M_{q_{v}} は B.4 節で述べる. 密度は, 湿潤空気の状態方程式

ρ = \frac{p_{0}}{R_{d}} \frac{Π^{c_{v}/c_{p}}}{θ_{v}}  (B.26)

から診断的に得られる.

B.3 サブグリッドスケールの拡散

乾燥大気の場合と同様に, レイノルズ平均を用いて, グリッドスケールの運動とサブグリッドスケールの運動に分離する. 直交直線座標系における運動方程式, 連続の式に対してレイノルズ平均操作を行うと,

\frac{∂\overline{u}_{i}}{∂t} + \overline{u}_{j} \frac{∂\overline{u}_{i}}{∂x_{j}} + c_{p}Θ̄_{v} \frac{∂\overline{π}}{∂x_{i}} = δ_{i3}g [\frac{θ}{Θ̄} − 1 + 0.61 (q_{v} − Q̄_{v}) − q_{c} − q_{r}]
−ε_{i3k}fu_{k} − \frac{∂}{∂x_{j}} \overline{u′′_{i}u′′_{j}},  (B.27)
\frac{∂\overline{ρ}}{∂t} + \frac{∂\overline{ρu}_{i}}{∂x_{i}} = 0  (B.28)

を得る. ただし上付バーの混同を避けるため本節でのみ, 基本場での θ, θ_{v}, q_{v} をそれぞれ Θ̄, Θ̄_{v}, Q̄_{v} と表す.

また, 熱力学の式, 水蒸気・雲水・雨水の混合比の式に対しても同様な操作を行うと,

\frac{∂\overline{ϕ}}{∂t} + \overline{u_{j}} \frac{∂\overline{ϕ}}{∂x_{j}} = −\frac{∂}{∂x_{j}} \overline{u′′_{j}ϕ′′}.  (B.29)

ただし ϕ = θ, q_{v}, q_{c}, q_{r} である. 上式 (B.29) 中の −\frac{∂}{∂x_{j}} \overline{u′′_{j}ϕ′′} が, それぞれ式 (B.11)

〜(B.14) の拡散項 D_{ϕ} である.

(B.27) の右辺第 3 項に\overline{u′′_{i}u′′_{j}}, (B.29) の右辺に\overline{u′′_{j}ϕ′′} があるために, \overline{u}_{i}, \overline{π}, \overline{θ}, \overline{q}_{v}, \overline{q}_{c}, \overline{q}_{r} に関して方程式系が閉じていない. 方程式を閉じるために, 本数理モデルでは乾燥大気と同様に Klemp and Wilhelmson (1978) に従う. このとき\overline{u′′_{i}u′′_{j}} と\overline{u′′_{j}ϕ′′} を以下のように仮定する:

\overline{u′′_{i}u′′_{j}} = −K_{m} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} ) + \frac{2}{3}δ_{ij}E,  (B.30)
\overline{u′′_{j}ϕ′′} = K_{h} \frac{∂ϕ}{∂x_{j}} .  (B.31)

E, K_{m} および K_{h} を求める. E は以下の時間発展方程式から求められる:

\frac{dE}{dt} = g \{−AK_{h} \frac{∂θ_{e}}{∂z} + K_{h} (1 − γA) \frac{∂q_{v}}{∂z} + K_{h} \frac{∂q_{c}}{∂z} \}

+2K_{m} (\frac{∂u_{j}}{∂x_{j}} )^{2} + \frac{1}{2} |ε_{ijk}|K_{m} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} )^{2}

−\frac{3}{2}δ_{ij}E (\frac{∂u_{j}}{∂x_{j}} ) + \frac{∂}{∂x_{j}} (K_{m} \frac{∂E}{∂x_{j}} ) − (\frac{C_{ε}}{l} ) E^{\frac{3}{2}} ,  (B.32)

A = \frac{1}{θ̄} \{\frac{1 + \frac{1.61εLq_{v}}{R_{d}T}}{1 + \frac{εL^{2}q_{v}}{c_{p}R_{d}T^{2}}} \} .  (B.33)

ただし θ_{e} は相当温位, γ = L/ (c_{p}Π̄) , θ_{e} = θ (1 + Lq_{v}/(c_{p}T)) である.

E が得られると, K_{m} は (2.17) から求められる. Deardorff (1975) によると C_{ε} = C_{m} = 0.2 である. CReSS (坪木ら, 2001) では C_{ε} = 3.9 (最下層), C_{ε} = 0.93 (それ以外), C_{m} = 0.1 としている.

式 (B.32) を数値計算したとき E が負になると扱いが困る. 数値モデルでは, 関係式 (2.18) と (2.17) を用いて K_{m} のみにした式を用いる:

\frac{dK_{m}}{dt} = 3g \{−AK_{m} \frac{∂θ_{e}}{∂z} + K_{m} (1 − γA) \frac{∂q_{v}}{∂z} + K_{m} \frac{∂q_{c}}{∂z} \}

+2C_{m}^{2}l^{2} (\frac{∂u_{j}}{∂x_{j}} )^{2} + \frac{1}{2} |ε_{ijk}|C_{m}^{2}l^{2} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} )^{2} − \frac{1}{3}δ_{ij}K_{m} (\frac{∂u_{j}}{∂x_{j}} )

+\frac{1}{2}δ_{ij} (\frac{∂^{2}K_{m}^{2}}{∂x^{2}} ) + δ_{ij} (\frac{∂K_{m}}{∂x_{j}} )^{2} − \frac{C_{ε}}{2C_{m}l^{2}} K_{m}^{2}.  (B.34)

地形に沿った座標系における各方程式の拡散項は以下のようになる. x, y, z 方向の運動方程式の拡散項は, それぞれ,

D_{u} = 2 (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} )\}

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂u}{∂y} + \frac{∂v}{∂x} + G^{23} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂v}{∂z^{∗}} )\}

+\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} \{K_{m} (\frac{∂w}{∂x} + \frac{1}{G^{1/2}} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂w}{∂z^{∗}} )\}

−\frac{2}{3} (\frac{∂E}{∂x} + G^{13} \frac{∂E}{∂z^{∗}} + \frac{∂E}{∂y} + G^{23} \frac{∂E}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂E}{∂z^{∗}} ) ,  (B.35)

\underline{付録 B 湿潤大気における準圧縮モデルの定式化 51}

D_{v} = (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂v}{∂x} + \frac{∂u}{∂y} + G^{23} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂v}{∂z^{∗}} )\}

+2 (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} )\}

+\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} \{K_{m} (\frac{∂w}{∂y} + \frac{1}{G^{1/2}} \frac{∂v}{∂z^{∗}} + G^{23} \frac{∂w}{∂z^{∗}} )\}

−\frac{2}{3} (\frac{∂E}{∂x} + G^{13} \frac{∂E}{∂z^{∗}} + \frac{∂E}{∂y} + G^{23} \frac{∂E}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂E}{∂z^{∗}} ) ,  (B.36)
D_{w} = (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂w}{∂x} + \frac{1}{G^{1/2}} \frac{∂u}{∂z^{∗}} + G^{13} \frac{∂w}{∂z^{∗}} )\}

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) \{K_{m} (\frac{∂w}{∂y} + \frac{1}{G^{1/2}} \frac{∂v}{∂z^{∗}} + G^{23} \frac{∂w}{∂z^{∗}} )\}

+2\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} (K_{m} \frac{1}{G^{1/2}} \frac{∂w}{∂z^{∗}} )

−\frac{2}{3} (\frac{∂E}{∂x} + G^{13} \frac{∂E}{∂z^{∗}} + \frac{∂E}{∂y} + G^{23} \frac{∂E}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂E}{∂z^{∗}} )  (B.37)

となる. 熱力学の式, 水蒸気・雲水・雨水混合比の式中の拡散項は, それぞれ,

D_{ϕ} = (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) (K_{h} \frac{∂ϕ}{∂x} + K_{h}G^{13} \frac{∂ϕ}{∂z^{∗}} )

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) (K_{h} \frac{∂ϕ}{∂y} + K_{h}G^{23} \frac{∂ϕ}{∂z^{∗}} )

+\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} (\frac{K_{h}}{G^{1/2}} \frac{∂ϕ}{∂z^{∗}} )  (B.38)

となる.

乱流運動エネルギー方程式は

\frac{∂E}{∂t} + u\frac{∂E}{∂x} + v \frac{∂E}{∂y} + ω \frac{∂E}{∂z^{∗}}
= \frac{g}{G^{1/2}} \{−AK_{h} \frac{∂θ_{e}}{∂z^{∗}} + K_{h} (1 − γA) \frac{∂q_{v}}{∂z^{∗}} + K_{h} \frac{∂q_{c}}{∂z^{∗}} \}

+2K_{m} \{(\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} )^{2} + (\frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} )^{2} + (\frac{1}{G^{1/2}} \frac{∂w}{∂z^{∗}} )^{2}\}

+K_{m}\{(\frac{∂v}{∂x} + \frac{∂u}{∂y} + G^{13} \frac{∂v}{∂z^{∗}} + G^{23} \frac{∂u}{∂z^{∗}} )^{2}

+ (\frac{∂w}{∂y} + G^{23} \frac{∂w}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂v}{∂z^{∗}} )^{2}

+ (\frac{∂w}{∂x} + G^{13} \frac{∂w}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂u}{∂z^{∗}} )^{2}\}

−\frac{2}{3}E (\frac{∂u}{∂x} + \frac{∂v}{∂y} + G^{13} \frac{∂u}{∂z^{∗}} + G^{23} \frac{∂v}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂w}{∂z^{∗}} )

+ (\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) (K_{h} \frac{∂E}{∂x} + K_{h}G^{13} \frac{∂E}{∂z^{∗}} )

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) (K_{h} \frac{∂E}{∂y} + K_{h}G^{23} \frac{∂E}{∂z^{∗}} )

+\frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} (\frac{K_{m}}{G^{1/2}} \frac{∂E}{∂z^{∗}} ) − (\frac{C_{ε}}{l} ) E^{\frac{3}{2}}  (B.39)

である. 運動量渦拡散係数に関する方程式は

\frac{∂K_{m}}{∂t} = − (u\frac{∂K_{m}}{∂x} + v \frac{∂K_{m}}{∂y} + ω \frac{∂K_{m}}{∂z^{∗}} ) + \frac{3gC_{m}^{2}l^{2}}{2G^{1/2}} (−A\frac{∂θ_{e}}{∂z^{∗}} + \frac{∂q_{l}}{∂z^{∗}} )

+C_{m}^{2}l^{2}\{(\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} )^{2} + (\frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} )^{2} + (\frac{1}{G^{1/2}} \frac{∂w}{∂z} )^{2}\}

+\frac{C_{m}^{2}l^{2}}{2} \{(\frac{∂v}{∂x} + G^{13} \frac{∂v}{∂z^{∗}} + \frac{∂u}{∂y} + G^{23} \frac{∂u}{∂z^{∗}} )^{2}

+ (\frac{∂w}{∂x} + G^{13} \frac{∂w}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂u}{∂z})^{2}

+ (\frac{∂w}{∂y} + G^{23} \frac{∂w}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂v}{∂z} )^{2}\}

−\frac{K_{m}}{3} (\frac{∂u}{∂x} + G^{13} \frac{∂u}{∂z^{∗}} + \frac{∂v}{∂y} + G^{23} \frac{∂v}{∂z^{∗}} + \frac{1}{G^{1/2}} \frac{∂w}{∂z} )

+\frac{1}{2}\{(\frac{∂}{∂x} + G^{13} \frac{∂}{∂z^{∗}} ) (\frac{∂K_{m}^{2}}{∂x} + G^{13} \frac{∂K_{m}^{2}}{∂z^{∗}} )

+ (\frac{∂}{∂y} + G^{23} \frac{∂}{∂z^{∗}} ) (\frac{∂K_{m}^{2}}{∂y} + G^{23} \frac{∂K_{m}^{2}}{∂z^{∗}} ) + \frac{1}{G^{1/2}} \frac{∂}{∂z^{∗}} \}

+ (\frac{∂K_{m}}{∂x} + G^{13} \frac{∂K_{m}}{∂z^{∗}} )^{2} + (\frac{∂K_{m}}{∂y} + G^{23} \frac{∂K_{m}}{∂z^{∗}} )^{2} + (\frac{1}{G^{1/2}} \frac{∂K_{m}}{∂z} )^{2}

−\frac{C_{ε}}{2C_{m}l^{2}} K_{m}^{2}  (B.40)

となる.

B.4 雲・降水の微物理過程

水蒸気を含む空気塊が上昇すると, 水蒸気から雲水になり雲が形成され, さらに降水粒子へと成長すると降水が起こる. 降水過程には暖かい雨と冷たい雨の 2 種類がある. 暖かい雨は氷相過程を含まない降水過程である. このような雨をもらたす雲は全領域が 0^{∘}C 以上にあり, 暖かい雲といわれる. 冷たい雨は氷相過程を含む降水過程である. このような雨をもたらす雲は, 雲内の全領域または一部が 0^{∘}C 以下にあり, 冷たい雲といわれる. ここでは暖かい雨のみを考える.

暖かい雲のモデル化にはバルク法のパラメタリゼーションを用いる. 暖かい雲についてのバルク法では, 水物質に関する予報変数として水蒸気混合比 q_{v}, 雲水混合比 q_{c} および雨水混合比 q_{r} の 3 つを用いる. 雲・降水の微物理過程が関係するのはこれら 3 つの混合比の式 (B.23)〜 (B.25), および熱力学の式 (B.22) である. 直交直線座標系における方程式の雲・降水微物理過程に伴う生成・消滅項は, 以下のように記述される:

M_{θ} = −γ (\frac{dq_{vs}}{dt} + E_{r}) ,  (B.41)
M_{q_{v}} = \frac{dq_{vs}}{dt} + E_{r},  (B.42)
M_{q_{c}} = −\frac{dq_{vs}}{dt} − A_{r} − C_{r},  (B.43)
M_{q_{r}} = F_{r} − E_{r} + A_{r} + C_{r}.  (B.44)

q_{vs}, E_{r}, A_{r}, C_{r}, F_{r} は以下で述べる.

本モデルで考慮する雲・降水の微物理過程は以下の通りである (図 B.1). 水蒸気が直接雨水になる凝結する過程は無視する.

図 B.1: 考慮する微物理過程: 凝結または蒸発による水蒸気から雲水間の変換, 蒸発による雨水から水蒸気への変換, 雲水から雨水への変換, 雨水の鉛直フラックスである.

• 凝結または蒸発による水蒸気から雲水間の変換: \frac{dq_{vs}}{dt}

この変換を計算するときは, Soong and Ogura (1973) の湿潤飽和調節法 (D.4 節参照) を用いる. q_{vs} は飽和混合比であり, Tetens の公式を用いて,

q_{vs} = \frac{3.8}{p̄} exp (17.27\frac{Π̄θ − 273}{Π̄θ − 36})  (B.45)

と表される.

• 蒸発による雨水から水蒸気への変換: E_{r} この変換は Ogura and Takahashi (1971) と同様に

E_{r} = \frac{1}{ρ̄}\frac{(1 − \frac{q_{v}}{q_{vs}} ) C (ρ̄q_{r})^{0.525}}{5.4 × 10^{5} + \frac{2.55 × 10^{6}}{p̄q_{vs}}}  (B.46)

と表される. ただし C は ventilation factor で C = 1.6 + 124.9 (ρ̄q_{r})^{0.2046} である.

• 雲水から雨水への変換: A_{r}, C_{r}

併合成長による変換 A_{r} と衝突併合による変換 C_{r} は, Kessler (1969) のパラメタリゼーションを用いて

A_{r} = k_{1} (q_{c} − a) ,  (B.47)
C_{r} = k_{2}q_{c}q_{r}^{0.875}  (B.48)

と書ける. ただし k_{1} = 0.001 s^{−1}, a = 0.001 kg kg^{−1}, k_{2} = 2.2 s^{−1} である.

• 雨水の鉛直フラックス: F_{r}

F_{r} = \frac{1}{ρ̄}\frac{∂}{∂z} (ρ̄V q_{r}) .  (B.49)

V は雨水の終端速度である. Soong and Ogura (1973) の終端速度に密度の変化を加えて

V = 36.34 (ρ̄q_{r})^{0.1346} (\frac{ρ̄}{ρ_{0}} )^{− \frac{1}{2}}  (B.50)

のように与えられる. ただし ρ_{0} は基本場の地表面での密度である.

地形に沿った座標系における生成・消滅項は, それぞれ,

M_{θ} = −γ \{\frac{∂q_{vs}}{∂t} + u\frac{∂q_{vs}}{∂x} + v \frac{∂q_{vs}}{∂y} + ω \frac{∂q_{vs}}{∂z^{∗}} + E_{r}\} ,  (B.51)
M_{q_{v}} = \frac{∂q_{vs}}{∂t} + u\frac{∂q_{vs}}{∂x} + v \frac{∂q_{vs}}{∂y} + ω \frac{∂q_{vs}}{∂z^{∗}} + E_{r},  (B.52)
M_{q_{c}} = − (\frac{∂q_{vs}}{∂t} + u\frac{∂q_{vs}}{∂x} + v \frac{∂q_{vs}}{∂y} + ω \frac{∂q_{vs}}{∂z^{∗}} ) − A_{r} − C_{r},  (B.53)
M_{q_{r}} = F_{r} − E_{r} + A_{r} + C_{r}  (B.54)

である.

\underline{付録 C 乾燥大気における準圧縮モデルの離散化 56}

付録C {乾燥大気における準圧縮モ \atop デルの離散化}

本章では空間微分と時間積分および計算不安定を抑える方法と, 乾燥大気と湿潤大気における基礎方程式系の離散化について述べる. 本数値モデルの数値計算法は以下のとおりである:

• 空間微分は有限差分法で, スタッガード格子.

• 時間微分は, 音波に関係する項を短い時間間隔で解き, それ以外の項を長い時間間隔で解く時間分割法. 水平方向の運動方程式は陽解法, 鉛直方向の運動方程式と圧力方程式は陰解法.

• 移流項は 2 次精度の中心差分, 前進差分の数値粘性項を追加.

• 上下境界条件は w = 0, 側面境界は周期境界条件.

C.1 空間微分の離散化

空間微分は中心差分で離散化する. 格子点の構造と変数の配置を図 C.1 と図 C.2 に示す. 変数はスタッガードに配置され, 水平方向を Arakawa C グリッド, 鉛直方向を Lorenz グリッドとした. このとき, すべてのスカラー量 (気圧, 温位, 水蒸気混合比, 雲水混合比, 雨水混合比) は格子ボックスの中心に配置される. ベクトル量の水平成分 (速度 u, v) は, 鉛直方向をスカラー量と同じレベルにして, u を x 負方向に半格子ずらし, v を y 負方向に半格子ずらす. ベクトル量の鉛直成分 w は, 水平方向をスカラー量と同じレベルにして, z 負方向に半格子ずらす.

また, 座標変換に伴うヤコビアン G^{1/2} および G^{13}, G^{23} は以下のように配置される. G^{1/2} は格子ボックスの中心に配置される. G^{13} は, u の定義点から z 負方向に半格子ずらした点に配置される. G^{23} は, v の定義点からさらに z 負方向に半格子ず {図 C.1: (i, j, k) 番目の格子における変数のスタッガードな配置図. ベクトル量は \atop 半格子ずれ, スカラー量は格子の中心に配置される.} 図 C.2: 全計算領域における変数の配置図. 黒丸は半格子ずらしたベクトル量, × 印はスカラー量の位置を表す. 計算領域の外の格子は境界における微分で必要な 格子である.

らした点に配置される. 座標 x, y, z^{∗} についても半格子ずらし, ベクトル量と同じ位置に配置される.

このようなスタッガード格子を用いる場合, ある位置 (i, j, k) での値を求めるためには平均操作が必要である. この平均化のオペレータを以下のように定義する.

任意のスカラー量の変数 ϕ を空間の x 方向 (または y, z^{∗} 方向) に平均する. 2 次精度の中心差分に現れる 2 点の平均は, 例えば,

\overline{ϕ}_{i,j,k}^{x} ≡ \frac{1}{2} (ϕ_{i− \frac{1}{2} ,j,k} + ϕ_{i+ \frac{1}{2} ,j,k})  (C.1)

とする. ただし下付添字 i, j, k は側面からそれぞれ x, y, z 方向に i, j, k 番目の格子点であることを表す.

平均化の添字が 2 つある場合, 例えば x, z^{∗} 方向の平均化は

\overline{ϕ}_{i,j,k}^{xz^{∗}} ≡ \frac{1}{4} (ϕ_{i− \frac{1}{2} ,j,k− \frac{1}{2}} + ϕ_{i+ \frac{1}{2} ,j,k− \frac{1}{2}} + ϕ_{i− \frac{1}{2} ,j,k+ \frac{1}{2}} + ϕ_{i+ \frac{1}{2} ,j,k+ \frac{1}{2}} )  (C.2)

とする.

微分オペレータは

(∂_{x}ϕ)_{i,j,k} ≡ \frac{1}{Δx} (ϕ_{i+ \frac{Δx}{2} ,j,k} − ϕ_{i− \frac{Δx}{2} ,j,k})  (C.3)

と定義する. 平均化オペレータと同様に, 下付添字 i は側面から i 番目の格子点であることを表し, ある方向 x の格子間隔を Δx とする.

2 階の微分オペレータは以下のように定義する:

(∂_{x}^{2}ϕ)_{i,j,k} ≡ \frac{1}{Δx} (\frac{ϕ_{i+Δx,j,k} − ϕ_{i,j,k}}{Δx} − \frac{ϕ_{i,j,k} − ϕ_{i−Δx,j,k}}{Δx} )
= \frac{1}{Δx^{2}} (ϕ_{i+Δx,j,k} − 2ϕ_{i,j,k} + ϕ_{i−Δx,j,k}) .  (C.4)

直交直線座標での鉛直速度 w と地形に沿う座標系における鉛直速度 ω を同じ点で定義する:

ω = \overline{\overline{u}^{z^{∗}}G^{13}}^{x} + \overline{\overline{v}^{z^{∗}}G^{23}}^{y} + \frac{w}{\overline{G^{1/2}}^{z^{∗}}} .  (C.5)

C.2 時間積分の離散化

大気中には音波, 内部重力波やロスビー波が存在する. 普通, 音波は気象学的に興味がない. 数値計算を行う際音波が CFL 条件を決めるため, 時間間隔を非常に小さくしなければならないという問題が生じる. このような音波による問題を回避するために, 音波に関する項をとそれ以外の項を別の時間間隔で解く時間分割法 (またはモード分離時間積分法) で解く. さらに音波に関する項は鉛直方向のみ陰解法 (鉛直陰解法) で解く.

時間分割法 (図 C.3) では, 音波に関係する項 (u, v, w, π) を短い時間間隔 Δτ で解き, それ以外の項を長い時間間隔 Δt で解く. ただし Δt を n 分割したものを Δτ とする:

Δt = nΔτ.  (C.6)

また, 音波を含む方程式中の音波に関係しない項を長い時間間隔で解き, t から t + Δt までは t での値 (一定値) を用いる.

図 C.3: 時間分割法の概念図. 音波に関係しない量を長い時間間隔で, 音波に関係する量を短い時間間隔で解く. 短い時間間隔で解いているとき, 音波に関係しない量は一定値である.

C.3 乾燥大気における基礎方程式系の離散化

予報変数の τ + Δτ, t + Δt における値を数値的に求めるために, 第 2 章に示した基礎方程式系を離散化する. 2.3.2 節と 2.3.3 節で述べたように, 基礎方程式には計算不安定を抑えるために音波減衰項と数値粘性項を加えてある.

C.3.1 水平方向運動方程式

(i, j, k) 番目の格子について, 水平方向の運動方程式 (2.29) と (2.30) を陽解法で離散化する. これら方程式中の音波に関係する項を短い時間間隔 Δτ で, その他の項を長い時間間隔 Δt で離散化する.

式 (2.29) と (2.30) は, それぞれ,

\frac{u_{i− \frac{Δx}{2} ,j,k}^{τ+Δτ} − u_{i− \frac{Δx}{2} ,j,k}^{τ}}{Δτ}

= −c_{p}θ̄[\{(∂_{x}π)_{i− \frac{Δx}{2} ,j,k} − α (∂_{x}Div)_{i− \frac{Δx}{2} ,j,k}\}^{τ}

+\overline{G^{13}}^{z^{∗}} \{(∂_{z^{∗}} \overline{π}^{xz^{∗}})_{i− \frac{Δx}{2} ,j,k} − α (∂_{z^{∗}} \overline{Div}^{xz^{∗}})_{i− \frac{Δx}{2} ,j,k}\}^{τ} ]

−Adv.u_{i− \frac{Δx}{2} ,j,k}^{t} + Cori.u_{i− \frac{Δx}{2} ,j,k}^{t} + (D_{u})_{i− \frac{Δx}{2} ,j,k}^{t−Δt} ,  (C.7)

\frac{v_{i,j− \frac{Δy}{2} ,k}^{τ+Δτ} − v_{i,j− \frac{Δy}{2} ,k}^{τ}}{Δτ}

= −c_{p}θ̄[\{(∂_{y}π)_{i,j− \frac{Δy}{2} ,k} − α (∂_{x}Div)_{i,j− \frac{Δy}{2} ,k}\}^{τ}
+\overline{G^{23}}^{z^{∗}} \{(∂_{z^{∗}} \overline{π}^{yz^{∗}})_{i,j− \frac{Δy}{2} ,k} − α (∂_{z^{∗}} \overline{Div}^{yz^{∗}})_{i,j− \frac{Δy}{2} ,k}\}^{τ} ]
−Adv.v_{i,j− \frac{Δy}{2} ,k}^{t} + Cori.v_{i,j− \frac{Δy}{2} ,k}^{t} + (D_{v})_{i,j− \frac{Δy}{2} ,k}^{t−Δt}  (C.8)

と離散化される. ただし

Div_{i,j,k} = (∂_{x}u)_{i,j,k} + (∂_{y}v)_{i,j,k} + \overline{G^{13}}_{i,j,k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{i,j,k}
+\overline{G^{23}}_{i,j,k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{xz^{∗}})_{i,j,k} + \frac{1}{G_{i,j,k}^{1/2}} (∂_{z^{∗}}w)_{i,j,k}  (C.9)

であり, スカラー量と同じ位置にある.

水平方向の運動方程式 (C.7) と (C.8) をそれぞれ数値計算すると, u^{τ+Δτ} と v^{τ+Δτ} が得られる.

長い時間間隔 t で計算する項を離散化したものは以下のようになる. 前述したように, これらの項は t から t + Δt までの間, t における値 (一定値) を用いて計算される. 移流項は

Adv.u_{i− \frac{Δx}{2} ,j,k}^{t} = (u∂_{x}\overline{u}^{x})_{i− \frac{Δx}{2} ,j,k}^{t} + (\overline{v}^{xy}∂_{y}\overline{u}^{y})_{i− \frac{Δx}{2} ,j,k}^{t}

+ (\overline{ω}^{xz^{∗}}∂_{z^{∗}} \overline{u}^{z^{∗}})_{i− \frac{Δx}{2} ,j,k}^{t} + Diff.u_{i− \frac{Δx}{2} ,j,k}^{t},  (C.10)
Adv.v_{i,j− \frac{Δy}{2} ,k}^{t} = (\overline{u}^{xy}∂_{x}\overline{v}^{x})_{i,j− \frac{Δx}{2} ,k}^{t} + (v∂_{y}\overline{v}^{y})_{i,j− \frac{Δx}{2} ,k}^{t}
+ (\overline{ω}^{yz^{∗}}∂_{z^{∗}} \overline{v}^{z^{∗}})_{i,j− \frac{Δy}{2} ,k}^{t} + Diff.v_{i,j− \frac{Δx}{2} ,k}^{t}.  (C.11)

数値粘性項は

Diff.u_{i− \frac{Δx}{2} ,j,k}^{t} = ν_{2h} [∂_{x}^{2}u + ∂_{y}^{2}u]_{i− \frac{Δx}{2} ,j,k}^{t} + ν_{2v} [∂_{z^{∗}}^{2}u]_{i− \frac{Δx}{2} ,j,k}^{t} , (C.12)
Diff.v_{i,j− \frac{Δy}{2} ,k}^{t} = ν_{2h} [∂_{x}^{2}v + ∂_{y}^{2}v]_{i,j− \frac{Δy}{2} ,k}^{t} + ν_{2v} [∂_{z^{∗}}^{2}v]_{i,j− \frac{Δy}{2} ,k}^{t} . (C.13)

コリオリ項は

Cori.u_{i− \frac{Δx}{2} ,j,k}^{t} = [\overline{f\overline{v}^{y}}^{x}]_{i− \frac{Δx}{2} ,j,k}^{t} ,  (C.14)
Cori.v_{i,j− \frac{Δy}{2} ,k}^{t} = [\overline{f\overline{u}^{x}}^{y}]_{i,j− \frac{Δy}{2} ,k}^{t} .  (C.15)

拡散項 (右辺の下付添字と上付添字 t − Δt は省略した) は

D_{i− \frac{Δx}{2} ,j,k} = 2∂_{x} \{K_{m} (∂_{x}u) + K_{m}\overline{G^{13}}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{z^{∗}})\}

+∂_{y}\{\overline{K_{m}}^{xy} (∂_{y}u) + \overline{K_{m}}^{xy} (∂_{x}v) + \overline{K_{m}}^{xy}\overline{G^{23}}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{yz∗})

+\overline{K_{m}}^{xy}\overline{G^{13}}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{xz^{∗}})\}

+2\overline{G^{13}}^{z^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{xz^{∗}} (∂_{x}\overline{u}^{xz^{∗}}) + \overline{K_{m}}^{xz^{∗}}G^{13} (∂_{z^{∗}}u)\}

+\overline{G^{23}}^{xyz^{∗}}∂_{z^{∗}}\{\overline{K_{m}}^{xz^{∗}} (∂_{y}\overline{u}^{yz^{∗}}) + \overline{K_{m}}^{xz^{∗}}∂_{x}\overline{v}^{yz^{∗}}

+\overline{K_{m}}^{xz^{∗}} \overline{G^{23}}^{xy} (∂_{z^{∗}}u) + \overline{K_{m}}^{xz^{∗}} \overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{v}^{xy})\}

+\frac{1}{\overline{G^{1/2}}^{x}} ∂_{z^{∗}} \{\overline{K_{m}}^{xz^{∗}} (∂_{x}w) + \begin{matrix}\overline{K_{m}} & xz^{∗} \\ \overline{G^{1/2}} & xz^{∗} \end{matrix} (∂_{z^{∗}}u) + \overline{K_{m}}^{xz^{∗}}G^{13}∂_{z^{∗}} \overline{w}^{xz^{∗}} \}

−\frac{2}{3C_{m}^{2}l^{2}} \{(∂_{x}K_{m}^{2}) + \overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{xz^{∗}})

+ (∂_{y}\overline{K_{m}^{2}}^{xy}) + \overline{G^{23}}^{xz^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})
+\frac{1}{\overline{G^{1/2}}^{x}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{xz^{∗}})\},  (C.16)

D_{i,j− \frac{Δy}{2} ,k} = ∂_{x}\{\overline{K_{m}}^{xy} (∂_{x}v) + \overline{K_{m}}^{xy} (∂_{y}u)
+\overline{K_{m}}^{xy}\overline{G^{23}}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{yz^{∗}}) + \overline{K_{m}}^{xy}\overline{G^{13}}^{yz^{∗}}∂_{z^{∗}} (\overline{v}^{xz^{∗}})\}

+2∂_{y} \{K_{m} (∂_{y}v) + K_{m}\overline{G^{23}}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})\}

+\overline{G^{13}}^{xyz^{∗}}∂_{z^{∗}}\{\overline{K_{m}}^{yz^{∗}} (\overline{∂_{x}v}^{xz^{∗}}) + \overline{K_{m}}^{yz^{∗}} (∂_{y}\overline{u}^{xz^{∗}})

+\overline{K_{m}}^{yz^{∗}}G^{23} (∂_{z^{∗}}) \overline{u}^{xy} + \overline{K_{m}}^{yz^{∗}} \overline{G^{13}}^{xy} (∂_{z^{∗}}v)\}

+2\overline{G^{23}}^{z^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{yz^{∗}} (∂_{y}\overline{v}^{yz^{∗}}) + \overline{K_{m}}^{yz^{∗}}G^{23} (∂_{z^{∗}}v)\}

−\frac{2}{3C_{m}^{2}l^{2}} \{(∂_{x}\overline{K_{m}^{2}}^{xy}) + \overline{G^{13}}^{yz^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})

+ (∂_{y}K_{m}^{2}) + \overline{G^{23}}^{z^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})
+\frac{1}{\overline{G^{1/2}}^{y}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})\}.  (C.17)

C.3.2 鉛直方向運動方程式と圧力方程式

鉛直方向の運動方程式 (2.31) と圧力方程式 (2.33) を陰解法で離散化する. これら方程式中の音波に関係する項を短い時間間隔 Δτ で, その他の項を長い時間間隔 Δt で離散化する. 陰解法で離散化するために, w, π は重み β(≠ 0) をつけて τ と τ + Δτ で平均をしたものを用いる. 本数値モデルでは β = 1 とする.

鉛直方向の運動方程式は

\frac{w_{i,j,k− \frac{Δz^{∗}}{2}}^{τ+Δτ} − w_{i,j,k− \frac{Δz^{∗}}{2}}^{τ}}{Δτ}

= α (\frac{c_{p}θ̄}{\overline{G^{1/2}}^{z^{∗}}} )_{i,j,k− \frac{Δz^{∗}}{2}} (∂_{z^{∗}}Div)_{i,j,k− \frac{Δz^{∗}}{2}}^{τ}

− (\frac{c_{p}θ̄}{\overline{G^{1/2}}^{z^{∗}}} )_{i,j,k− \frac{Δz^{∗}}{2}} [β (∂_{z^{∗}}π)_{i,j,k− \frac{Δz^{∗}}{2}}^{τ+Δτ} + (1 − β) (∂_{z^{∗}}π)_{i,j,k− \frac{Δz^{∗}}{2}}^{τ} ]

−Adv.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + Buoy.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + (D_{w})_{i,j,k− \frac{Δz^{∗}}{2}}^{t−Δt} ,  (C.18)

圧力方程式は

\underline{付録 C 乾燥大気における準圧縮モデルの離散化 63}
\frac{π_{i,j,k}^{τ+Δτ} − π_{i,j,k}^{τ}}{Δτ} + \frac{1}{G^{1/2}} \frac{c̄^{2}}{c_{p}ρ̄θ̄^{2}} \{β (∂_{z^{∗}} ρ̄θ̄w)_{i,j,k}^{τ+Δτ} + (1 − β) (∂_{z^{∗}} ρ̄θ̄w)_{i,j,k}^{τ}\}
= −\frac{c̄^{2}}{c_{p}θ̄}[(∂_{x}u)_{i,j,k} + (∂_{y}v)_{i,j,k} + \overline{G^{13}}_{i,j,k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{i,j,k}
+\overline{G^{23}}_{i,j,k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{i,j,k}]^{τ+Δτ}(C.19)

となる. 長い時間間隔 t で計算する項について離散化したものを書き下す. 移流項は

Adv.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} = (\overline{u}^{xz^{∗}}∂_{x}\overline{w}^{x})_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + (\overline{v}^{yz^{∗}}∂_{y}\overline{w}^{y})_{i,j,k− \frac{Δz^{∗}}{2}}^{t}
+ (ω∂_{z^{∗}} \overline{w}^{z^{∗}})_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + Diff.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} .  (C.20)

数値粘性項は

Diff.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} = ν_{2h} [∂_{x}^{2}w + ∂_{y}^{2}w]_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + ν_{2v} [∂_{z^{∗}}^{2}w]_{i,j,k− \frac{Δz^{∗}}{2}}^{t} . (C.21)

浮力項は

Buoy.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} = g [\overline{\frac{θ}{θ̄}}^{z^{∗}} − 1]_{i,j,k− \frac{Δz^{∗}}{2}}^{t} .  (C.22)

拡散項 (右辺の下付添字と上付添字 t − Δt は省略) は

D_{i,j,k− \frac{Δz^{∗}}{2}} = ∂_{x} \{\overline{K_{m}}^{xz^{∗}} (∂_{x}w) + \begin{matrix}\overline{K_{m}} & xz^{∗} \\ \overline{G^{1/2}} & xz^{∗} \end{matrix} (∂_{z^{∗}} \overline{u}^{xz^{∗}}) + \overline{K_{m}}^{xz^{∗}}G^{13} (∂_{z^{∗}} \overline{w}^{xz^{∗}})\}

+∂_{y} \{\overline{K_{m}}^{yz^{∗}} (∂_{y}w) + \begin{matrix}\overline{K_{m}} & yz^{∗} \\ \overline{G^{1/2}} & yz^{∗} \end{matrix} (∂_{z^{∗}}v) + \overline{K_{m}}^{yz^{∗}}G^{23} (∂_{z^{∗}} \overline{w}^{yz^{∗}})\}

+\overline{G^{13}}^{x}∂_{z^{∗}} \{K_{m} (∂_{x}\overline{w}^{xz^{∗}}) + \frac{K_{m}}{G^{1/2}} (∂_{z^{∗}} \overline{u}^{xz^{∗}}) + K_{m}\overline{G^{13}}^{xz^{∗}} (∂_{z^{∗}}w)\}

+\overline{G^{23}}^{y}∂_{z^{∗}} \{K_{m} (∂_{y}\overline{w}^{yz^{∗}}) + \frac{K_{m}}{G^{1/2}} (∂_{z^{∗}} \overline{v}^{yz^{∗}}) + K_{m}\overline{G^{23}}^{yz^{∗}} (∂_{z^{∗}}w)\}

+\frac{2}{\overline{G^{1/2}}^{z^{∗}}} ∂_{z^{∗}} \{\frac{K_{m}}{G^{1/2}} (∂_{z^{∗}}w)\}

−\frac{2}{3C_{m}^{2}l^{2}} \{(∂_{x}\overline{K_{m}^{2}}^{xz^{∗}}) + \overline{G^{13}}^{x} (∂_{z^{∗}}K_{m}^{2})

+ (∂_{y}\overline{K_{m}^{2}}^{yz^{∗}}) + \overline{G^{23}}^{y} (∂_{z^{∗}}K_{m}^{2}) + \frac{1}{\overline{G^{1/2}}^{z^{∗}}} (∂_{z^{∗}}K_{m}^{2})\}. (C.23)

式 (C.18) と (C.19) から直接 w^{τ+Δτ}, π^{τ+Δτ} を求めることはできない. これら 2 式を連立して w^{τ+Δτ} を消去し, π^{τ+Δτ} を求め, 式 (C.18) から w^{τ+Δτ} を求める.

以下で具体的に連立する. ここからは予報変数毎に格子番号をつけ, i, j を省略する. また w と同じ位置にある (i, j, k) 番目の変数を w_{w(k)} などと表す. 離散化した鉛直方向の運動方程式式 (C.18) より,

w_{w(k)}^{τ+Δτ} = w_{w(k)}^{τ}

− (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} )_{w(k)} \{−α (∂_{z^{∗}}Div)_{w(k)}^{τ} + (1 − β) (∂_{z^{∗}}π)_{w(k)}^{τ}\}

+F.w_{w(k)}^{t}Δτ − β (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} )_{w(k)} (∂_{z^{∗}}π)_{w(k)}^{τ+Δτ} ,  (C.24)
F.w_{w(k)}^{t} = −Adv.w_{w(k)}^{t} + Buoy.w_{w(k)}^{t} + (D_{w})_{w(k)}^{t−Δt} .  (C.25)

離散化した圧力方程式 (D.13) より,

π_{k}^{τ+Δτ} + β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{k} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{k}^{τ+Δτ}
= π_{k}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{k} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{k}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄} [(∂_{x}u)_{k} + (∂_{y}v)_{k}

+\overline{G^{13}}_{k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{k} + \overline{G^{23}}_{k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{k}]^{τ+Δτ}(C.26)

である.

(C.24) を代入して w^{τ+Δτ} を消去すると,

π_{k}^{τ+Δτ} − β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{k} [∂_{z^{∗}} \{(\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} ) (∂_{z^{∗}}π)\}]_{k}^{τ+Δτ}
= π_{k}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{k} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{k}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄} [(∂_{x}u)_{k} + (∂_{y}v)_{k} + \overline{G^{13}}_{k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{k} + \overline{G^{23}}_{k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{k}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{k}

×[∂_{z^{∗}} ρ̄θ̄\{w^{τ} − (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{w(k)}.  (C.27)

左辺について z^{∗} 方向の微分をオペレータなしで書き下すと,

(左辺) = π_{k}^{τ+Δτ}

−β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{k}

×\frac{1}{Δz^{∗}} [(\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k+1)} \frac{1}{Δz^{∗}} (π_{k+1} − π_{k})

− (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k)} \frac{1}{Δz^{∗}} (π_{k} − π_{k−1})]^{τ+Δτ}

= π_{k}^{τ+Δτ}

−β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{k} \frac{1}{Δz^{∗2}}

×[(\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k+1)} π_{k+1} − \{(\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k+1)} + (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k)}\} π_{k}

+ (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k)} π_{k−1}]^{τ+Δτ}.  (C.28)

本数値モデルにおける上部と下部の境界条件は, それぞれ

ω = 0  (z^{∗} = H; k = km + 1),  (C.29)
ω = u\frac{∂z_{s}}{∂x} + v \frac{∂z_{s}}{∂y}  (z^{∗} = 0; k = 0)  (C.30)

とする. 上部境界における鉛直方向の運動方程式を離散化したものは

β (∂_{z^{∗}}π)_{w(km+1)}^{τ+Δτ} = α (∂_{z^{∗}}Div)_{w(km+1)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(km+1)}^{τ}
+ (\frac{G^{1/2}}{c_{p}θ̄} F.w^{t})_{w(km+1)} .  (C.31)

このときの圧力方程式は

(左辺) = π_{km}^{τ+Δτ}

+β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{km} \frac{1}{Δz^{∗2}} [(\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(km)} (π_{km} − π_{km−1})]^{τ+Δτ}

(右辺) = π_{km}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{km} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{km}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄} [(∂_{x}u)_{km} + (∂_{y}v)_{km}

+\overline{G^{13}}_{km}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{km} + \overline{G^{23}}_{km}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{km}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{km}

×[∂_{z^{∗}} ρ̄θ̄\{w^{τ} − (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{km}

+β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{km} \frac{1}{Δz^{∗}} (\frac{c_{p}θ̄^{2}}{G^{1/2}} )_{w(km+1)}

×{α (∂_{z^{∗}}Div)_{w(km+1)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(km+1)}^{τ}

+ (\frac{G^{1/2}}{c_{p}θ̄} F.w^{t})_{w(km+1)}\}  (C.32)

である. また, 下部境界における鉛直方向の運動方程式を離散化したものは

β (∂_{z^{∗}}π)_{w(0)}^{τ+Δτ} = α (∂_{z^{∗}}Div)_{w(0)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(0)}^{τ}
+ (\frac{G^{1/2}}{c_{p}θ̄} F.w^{t})_{w(0)} .  (C.33)

このときの圧力方程式は

(左辺) = π_{0}^{τ+Δτ}

−β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{0} \frac{1}{Δz^{∗2}} [(\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(1)} (π_{1} − π_{0})]^{τ+Δτ}

(右辺) = π_{0}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{0} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{0}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄} [(∂_{x}u)_{0} + (∂_{y}v)_{0}

+\overline{G^{13}}_{0}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{1} + \overline{G^{23}}_{0}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{0}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{0}

×[∂_{z^{∗}} ρ̄θ̄\{w^{τ} − (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{0}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{1} \frac{1}{Δz^{∗}} (\frac{c_{p}θ̄^{2}}{G^{1/2}} )_{w(0)}

×\{α (∂_{z^{∗}}Div)_{w(0)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(0)}^{τ} + (\frac{G^{1/2}}{c_{p}θ̄} F.w^{t})_{w(0)}\}(C.34)

である.

これらの式は以下のように行列を用いてまとめて書くことができる:

( \begin{matrix}A_{0} & B_{0} & & 0 \\ C_{0} & ... & ... \\ & ... & ... & B_{km−1} \\ 0 & & C_{km−1} & A_{km} \end{matrix} ) ( \begin{matrix}π_{0} & π_{0} & ⋯ & π_{0} \\ π_{1} & ... & ... & ⋮ \\ ⋮ & & ... & ⋮ \\ π_{km} & ⋯ & ⋯ & π_{km} \end{matrix} )^{τ+Δτ}

= ( \begin{matrix}D_{0} & D_{0} & ⋯ & D_{0} \\ D_{1} & ... & ... & ⋮ \\ ⋮ & & ... & ⋮ \\ D_{km} & ⋯ & ⋯ & D_{km} \end{matrix} )^{τ} .  (C.35)

この連立 1 次方程式の係数行列は 3 項行列なので, LU 分解を行って計算すると大幅に計算量を減らすことができる. こうして求められた π^{τ+Δτ} と鉛直方向の運動方程式 (C.18) から w^{τ+Δτ} を求める. ただし

A_{k} = 1 + β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{k} \frac{1}{Δz^{∗2}} \{(\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k+1)} + (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k)}\} ,
(k = 1, 3, ⋯, km − 1),
A_{0} = 1 + β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{0} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(1)} ,

\underline{付録 C 乾燥大気における準圧縮モデルの離散化 68}

A_{km} = 1 + β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{km} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(km)} ,
B_{k} = −β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{k} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k+1)} ,
C_{k} = −β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{k+1} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}^{2}}{G^{1/2}} )_{w(k+1)} ,  (C.36)
D_{k} = π_{k}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{k} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{k}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄} [(∂_{x}u)_{k} + (∂_{y}v)_{k} + \overline{G^{13}}_{k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{k}

+\overline{G^{23}}_{k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{k}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{k}

×[∂_{z^{∗}} ρ̄θ̄\{w^{τ} − (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{k},

(k = 1, 3, ⋯, km − 3),  (C.37)
D_{0} = π_{0}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{0} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{0}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄} [(∂_{x}u)_{0} + (∂_{y}v)_{0}

+\overline{G^{13}}_{0}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{0} + \overline{G^{23}}_{1}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{0}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{0}

×[∂_{z^{∗}} ρ̄θ̄\{w^{τ} − (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{0}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{0} \frac{1}{Δz^{∗}} (\frac{c_{p}θ̄^{2}}{G^{1/2}} )_{w(0)}

×\{α (∂_{z^{∗}}Div)_{w(0)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(0)}^{τ} + (\frac{G^{1/2}}{c_{p}θ̄} F.w^{t})_{w(0)}\}, (C.38)

D_{km−2} = π_{km−2}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{km−2} \{∂_{z^{∗}} (ρ̄θ̄w)\}_{km−2}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄} [(∂_{x}u)_{km−2} + (∂_{y}v)_{km−2}

+\overline{G^{13}}_{km}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{km−2} + \overline{G^{23}}_{km}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{km−2}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄^{2}} )_{km−2}

×[∂_{z^{∗}} ρ̄θ̄\{w^{τ} − (\frac{c_{p}θ̄Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{km−2}

+β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄^{2}} )_{km−2} \frac{1}{Δz^{∗}} (\frac{c_{p}θ̄^{2}}{G^{1/2}} )_{w(km−1)}

×{α (∂_{z^{∗}}Div)_{w(km−1)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(km−1)}^{τ}

+ (\frac{G^{1/2}}{c_{p}θ̄} F.w^{t})_{w(km−1)}\}  (C.39)

C.3.3 熱力学の式

熱力学の式 (2.34) を中心差分で離散化する. 時間間隔は長い時間間隔 Δt である. 本節では変数の下付添字 i, j, k を省略する.

熱力学の式は

\frac{θ^{t+Δt} − θ^{t−Δt}}{2Δt} = −Adv.θ^{t} + D_{θ}^{t−Δt}.  (C.40)

となる.

前節までに求めた u^{τ+Δτ}, v^{τ+Δτ}, w^{τ+Δτ} および π^{τ+Δτ} を用いると θ^{τ+Δτ} を求めることができる. 2.3.4 節で述べたような計算不安定モードを抑えるため, Asselin (1972) の時間フィルタを用いる.

上式 (C.40) 中の各項は以下のように離散化される; 移流項は

Adv.θ^{t} = [\overline{u∂_{x}θ}^{x} + \overline{v∂_{y}θ}^{y} + \overline{ω∂_{z^{∗}}θ}^{z^{∗}}]^{t} + Diff.θ^{t}.  (C.41)

数値粘性項は

Diff.θ^{t} = ν_{2h} [∂_{x}^{2}θ + ∂_{y}^{2}θ]^{t} + ν_{2v} [∂_{z^{∗}}^{2}θ]^{t} .  (C.42)

拡散項は

D_{θ}^{t−Δt} = 3∂_{x} \{\overline{K_{m}}^{x} (∂_{x}θ) + \overline{K_{m}}^{x}\overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{θ}^{xz^{∗}})\}^{t−Δt}
+3∂_{y} \{\overline{K_{m}}^{y} (∂_{y}θ) + \overline{K_{m}}^{y}\overline{G^{23}}^{z^{∗}} (∂_{z^{∗}} \overline{θ}^{yz^{∗}})\}^{t−Δt}
+3\overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{x}\overline{θ}^{xz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{13}}^{x} (∂_{z^{∗}}θ)\}^{t−Δt}
+3\overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{y}\overline{θ}^{yz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{23}}^{y} (∂_{z^{∗}}θ)\}^{t−Δt}
+\frac{3}{G^{1/2}} ∂_{z^{∗}} \{\frac{K_{m}}{\overline{G^{1/2}}^{z^{∗}}} (∂_{z^{∗}}θ)\}^{t−Δt}  (C.43)

となる.

C.3.4 サブグリッドスケールの乱流運動エネルギー方程式

サブグリッドスケールの乱流運動エネルギー方程式を運動量渦拡散係数で書き換えた式 (B.40) は長い時間間隔 t で計算される. ここでは変数の下付添字 i, j, k は省略する.

\frac{K_{m}^{t+Δt} − K_{m}^{t−Δt}}{2Δt} = −Adv.K_{m}^{t} + \frac{3gC_{m}^{2}l^{2}}{2G^{1/2}θ̄} (−∂_{z^{∗}} \overline{θ}^{z^{∗}})^{t}

+C_{m}^{2}l^{2}\{ (∂_{x}u + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \overline{u}^{xz^{∗}})^{2}

+ (∂_{y}v + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \overline{v}^{yz^{∗}})^{2} + (\frac{1}{G^{1/2}} ∂_{z^{∗}}w)^{2}\}^{t}

+\frac{C_{m}^{2}l^{2}}{2} \{(∂_{x}\overline{v}^{xy} + ∂_{y}\overline{u}^{xy} + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \overline{v}^{yz^{∗}} + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \overline{u}^{xz^{∗}})^{2}

+ (∂_{x}w + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}}w + \frac{1}{G^{1/2}} ∂_{z^{∗}} \overline{u}^{xz^{∗}})^{2}

+ (∂_{y}\overline{w}^{yz^{∗}} + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}}w + \frac{1}{G^{1/2}} ∂_{z^{∗}} \overline{v}^{yz^{∗}})^{2}\}^{t}

+\frac{1}{2}\{∂_{x}^{2}K_{m}^{2} + ∂_{x} (\overline{G^{13}}^{z^{∗}}∂_{z^{∗}} \overline{K_{m}^{2}}^{xz^{∗}})

+∂_{z^{∗}} (\overline{G^{13}}^{x}∂_{x}\overline{K_{m}^{2}}^{xz^{∗}}) + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} (\overline{G^{13}}^{x}∂_{z^{∗}}K_{m}^{2})

+∂_{y}^{2}K_{m}^{2} + ∂_{y} (\overline{G^{23}}^{z^{∗}}∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})

+∂_{z^{∗}} (\overline{G^{23}}^{y}∂_{y}\overline{K_{m}^{2}}^{yz^{∗}}) + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} (\overline{G^{23}}^{y}∂_{z^{∗}}K_{m}^{2})

+\frac{1}{G^{1/2}} ∂_{z^{∗}} (\overline{\frac{1}{G^{1/2}}}^{∗}∂_{z^{∗}}K_{m}^{2})\}^{t}

+\{(∂_{x}\overline{K_{m}}^{x} + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \overline{K_{m}}^{z^{∗}})^{2}

+ (∂_{y}\overline{K_{m}}^{y} + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \overline{K_{m}}^{z^{∗}})^{2}

+ (\frac{1}{G^{1/2}} ∂_{z^{∗}} \overline{K_{m}}^{z^{∗}})^{2}\}^{t}

−\frac{C_{ε}}{2C_{m}l^{2}} K_{m}^{2}  (C.44)

となる. 移流項と数値粘性項はそれぞれ以下のようになる:

Adv.K_{m}^{t} = [\overline{u∂_{x}K_{m}}^{x} + \overline{v∂_{y}K_{m}}^{y} + \overline{ω∂_{z^{∗}}K_{m}}^{z^{∗}}]^{t} + Diff.K_{m}^{t},  (C.45)
Diff.K_{m}^{t} = ν_{2h} [∂_{x}^{2}K_{m} + ∂_{y}^{2}K_{m}]^{t} + ν_{2v} [∂_{z^{∗}}^{2}K_{m}]^{t} .  (C.46)

\underline{付録 D 湿潤大気における準圧縮モデルの離散化 72}

付録D {湿潤大気における準圧縮モ \atop デルの離散化}

本章では空間微分と時間積分および計算不安定を抑える方法と, 乾燥大気と湿潤大気における基礎方程式系の離散化について述べる. 本数値モデルの離散化は以下のとおりである:

• 空間微分の離散化は中心差分法を用いる. 変数の配置はスタッガード格子を用いる.

• 時間微分の離散化は, 音波に関係する項を短い時間間隔で解き, それ以外の項を長い時間間隔で解く時間分割法. 水平方向の音波に関する項は陽解法, 鉛直方向の音波に関する項は陰解法.

• 上下境界条件は w = 0, 側面境界は周期境界条件.

D.1 空間微分の離散化

湿潤大気についての空間微分の離散化は, 乾燥大気の場合と同様なので, C.1 節を参照のこと.

D.2 時間積分の離散化

湿潤大気についての時間積分の離散化は, 乾燥大気の場合と同様なので, C.2 節を参照のこと.

D.3 湿潤大気における基礎方程式系の離散化

予報変数の τ + Δτ, t + Δt での値を求めるために, 付録 B 章で求めた湿潤大気における基礎方程式系を離散化する. 2.3.2 節と 2.3.3 節で述べたように, 基礎方程式には計算不安定を抑えるために音波減衰項と数値粘性項を加えてある.

D.3.1 水平方向運動方程式

(i, j, k) 番目の格子について, 水平方向の運動方程式 (B.17) と (B.18) を陽解法で離散化する. これら方程式中の音波に関係する項は, 短い時間間隔 Δτ で離散化される.

式 (B.17) と (B.18) は, それぞれ,

\frac{u_{i− \frac{Δx}{2} ,j,k}^{τ+Δτ} − u_{i− \frac{Δx}{2} ,j,k}^{τ}}{Δτ}

= −c_{p}θ̄_{v}[\{(∂_{x}π)_{i− \frac{Δx}{2} ,j,k} − α (∂_{x}Div)_{i− \frac{Δx}{2} ,j,k}\}^{τ}

+\overline{G^{13}}^{z^{∗}} \{(∂_{z^{∗}} \overline{π}^{xz^{∗}})_{i− \frac{Δx}{2} ,j,k} − α (∂_{z^{∗}} \overline{Div}^{xz^{∗}})_{i− \frac{Δx}{2} ,j,k}\}^{τ} ]

−Adv.u_{i− \frac{Δx}{2} ,j,k}^{t} + Cori.u_{i− \frac{Δx}{2} ,j,k}^{t} + (D_{u})_{i− \frac{Δx}{2} ,j,k}^{t−Δt} ,  (D.1)

\frac{v_{i,j− \frac{Δy}{2} ,k}^{τ+Δτ} − v_{i,j− \frac{Δy}{2} ,k}^{τ}}{Δτ}

= −c_{p}θ̄_{v}[\{(∂_{y}π)_{i,j− \frac{Δy}{2} ,k} − α (∂_{x}Div)_{i,j− \frac{Δy}{2} ,k}\}^{τ}
+\overline{G^{23}}^{z^{∗}} \{(∂_{z^{∗}} \overline{π}^{yz^{∗}})_{i,j− \frac{Δy}{2} ,k} − α (∂_{z^{∗}} \overline{Div}^{yz^{∗}})_{i,j− \frac{Δy}{2} ,k}\}^{τ} ]
−Adv.v_{i,j− \frac{Δy}{2} ,k}^{t} + Cori.v_{i,j− \frac{Δy}{2} ,k}^{t} + (D_{v})_{i,j− \frac{Δy}{2} ,k}^{t−Δt}  (D.2)

と離散化される. ただし

Div_{i,j,k} = (∂_{x}u)_{i,j,k} + (∂_{y}v)_{i,j,k} + \overline{G^{13}}_{i,j,k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{i,j,k}
+\overline{G^{23}}_{i,j,k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{xz^{∗}})_{i,j,k} + \frac{1}{G_{i,j,k}^{1/2}} (∂_{z^{∗}}w)_{i,j,k}  (D.3)

であり, スカラー量と同じ位置にある.

水平方向の運動方程式 (D.1) と (D.2) を計算すると, u^{τ+Δτ} と v^{τ+Δτ} が得られる.

長い時間間隔 t で計算する項を離散化したものは以下のようになる. 前述したように, これらの項は t から t + Δt までの間, t における値 (一定値) を用いて計算される. 移流項は

Adv.u_{i− \frac{Δx}{2} ,j,k}^{t} = (u∂_{x}\overline{u}^{x})_{i− \frac{Δx}{2} ,j,k}^{t} + (\overline{v}^{xy}∂_{y}\overline{u}^{y})_{i− \frac{Δx}{2} ,j,k}^{t}
+ (\overline{ω}^{xz^{∗}}∂_{z^{∗}} \overline{u}^{z^{∗}})_{i− \frac{Δx}{2} ,j,k}^{t} + Diff.u_{i− \frac{Δx}{2} ,j,k}^{t},  (D.4)
Adv.v_{i,j− \frac{Δy}{2} ,k}^{t} = (\overline{u}^{xy}∂_{x}\overline{v}^{x})_{i,j− \frac{Δx}{2} ,k}^{t} + (v∂_{y}\overline{v}^{y})_{i,j− \frac{Δx}{2} ,k}^{t}
+ (\overline{ω}^{yz^{∗}}∂_{z^{∗}} \overline{v}^{z^{∗}})_{i,j− \frac{Δy}{2} ,k}^{t} + Diff.v_{i,j− \frac{Δx}{2} ,k}^{t}.  (D.5)

数値粘性項は

Diff.u_{i− \frac{Δx}{2} ,j,k}^{t} = ν_{2h} [∂_{x}^{2}u + ∂_{y}^{2}u]_{i− \frac{Δx}{2} ,j,k}^{t} + ν_{2v} [∂_{z^{∗}}^{2}u]_{i− \frac{Δx}{2} ,j,k}^{t} ,  (D.6)
Diff.v_{i,j− \frac{Δy}{2} ,k}^{t} = ν_{2h} [∂_{x}^{2}v + ∂_{y}^{2}v]_{i,j− \frac{Δy}{2} ,k}^{t} + ν_{2v} [∂_{z^{∗}}^{2}v]_{i,j− \frac{Δy}{2} ,k}^{t} .  (D.7)

コリオリ項は

Cori.u_{i− \frac{Δx}{2} ,j,k}^{t} = [\overline{f\overline{v}^{y}}^{x}]_{i− \frac{Δx}{2} ,j,k}^{t} ,  (D.8)
Cori.v_{i,j− \frac{Δy}{2} ,k}^{t} = [\overline{f\overline{u}^{x}}^{y}]_{i,j− \frac{Δy}{2} ,k}^{t} .  (D.9)

拡散項 (右辺の下付添字と上付添字 t − Δt は省略した) は

\underline{付録 D 湿潤大気における準圧縮モデルの離散化 75}

D_{i− \frac{Δx}{2} ,j,k} = 2∂_{x} \{K_{m} (∂_{x}u) + K_{m}\overline{G^{13}}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{z^{∗}})\}

+∂_{y}\{\overline{K_{m}}^{xy} (∂_{y}u) + \overline{K_{m}}^{xy} (∂_{x}v) + \overline{K_{m}}^{xy}\overline{G^{23}}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{yz∗})

+\overline{K_{m}}^{xy}\overline{G^{13}}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{xz^{∗}})\} +2\overline{G^{13}}^{z^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{xz^{∗}} (∂_{x}\overline{u}^{xz^{∗}}) + \overline{K_{m}}^{xz^{∗}}G^{13} (∂_{z^{∗}}u)\} +\overline{G^{23}}^{xyz^{∗}}∂_{z^{∗}}\{\overline{K_{m}}^{xz^{∗}} (∂_{y}\overline{u}^{yz^{∗}}) + \overline{K_{m}}^{xz^{∗}}∂_{x}\overline{v}^{yz^{∗}}

+\overline{K_{m}}^{xz^{∗}} \overline{G^{23}}^{xy} (∂_{z^{∗}}u) + \overline{K_{m}}^{xz^{∗}} \overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{v}^{xy})\}

+\frac{1}{\overline{G^{1/2}}^{x}} ∂_{z^{∗}} \{\overline{K_{m}}^{xz^{∗}} (∂_{x}w) + \begin{matrix}\overline{K_{m}} & xz^{∗} \\ \overline{G^{1/2}} & xz^{∗} \end{matrix} (∂_{z^{∗}}u) + \overline{K_{m}}^{xz^{∗}}G^{13}∂_{z^{∗}} \overline{w}^{xz^{∗}} \}

−\frac{2}{3C_{m}^{2}l^{2}} \{(∂_{x}K_{m}^{2}) + \overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{xz^{∗}})

+ (∂_{y}\overline{K_{m}^{2}}^{xy}) + \overline{G^{23}}^{xz^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})

+\frac{1}{\overline{G^{1/2}}^{x}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{xz^{∗}})\},  (D.10)
D_{i,j− \frac{Δy}{2} ,k} = ∂_{x}\{\overline{K_{m}}^{xy} (∂_{x}v) + \overline{K_{m}}^{xy} (∂_{y}u)

+\overline{K_{m}}^{xy}\overline{G^{23}}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{yz^{∗}}) + \overline{K_{m}}^{xy}\overline{G^{13}}^{yz^{∗}}∂_{z^{∗}} (\overline{v}^{xz^{∗}})\} +2∂_{y} \{K_{m} (∂_{y}v) + K_{m}\overline{G^{23}}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})\} +\overline{G^{13}}^{xyz^{∗}}∂_{z^{∗}}\{\overline{K_{m}}^{yz^{∗}} (\overline{∂_{x}v}^{xz^{∗}}) + \overline{K_{m}}^{yz^{∗}} (∂_{y}\overline{u}^{xz^{∗}})

+\overline{K_{m}}^{yz^{∗}}G^{23} (∂_{z^{∗}}) \overline{u}^{xy} + \overline{K_{m}}^{yz^{∗}} \overline{G^{13}}^{xy} (∂_{z^{∗}}v)\}

+2\overline{G^{23}}^{z^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{yz^{∗}} (∂_{y}\overline{v}^{yz^{∗}}) + \overline{K_{m}}^{yz^{∗}}G^{23} (∂_{z^{∗}}v)\}

−\frac{2}{3C_{m}^{2}l^{2}} \{(∂_{x}\overline{K_{m}^{2}}^{xy}) + \overline{G^{13}}^{yz^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})

+ (∂_{y}K_{m}^{2}) + \overline{G^{23}}^{z^{∗}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})

+\frac{1}{\overline{G^{1/2}}^{y}} (∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})\}.  (D.11)

D.3.2 鉛直方向運動方程式と圧力方程式

鉛直方向の運動方程式 (B.19) と圧力方程式 (B.21) を陰解法で離散化する. これら方程式は短い時間間隔 Δτ で離散化される. 鉛直陰解法では w, π は重み β(≠ 0) をつけて τ と τ + Δτ で平均をしたものを用いる.

鉛直方向の運動方程式は

\frac{w_{i,j,k− \frac{Δz^{∗}}{2}}^{τ+Δτ} − w_{i,j,k− \frac{Δz^{∗}}{2}}^{τ}}{Δτ}

= α (\frac{c_{p}θ̄_{v}}{\overline{G^{1/2}}^{z^{∗}}} )_{i,j,k− \frac{Δz^{∗}}{2}} (∂_{z^{∗}}Div)_{i,j,k− \frac{Δz^{∗}}{2}}^{τ}

− (\frac{c_{p}θ̄_{v}}{\overline{G^{1/2}}^{z^{∗}}} )_{i,j,k− \frac{Δz^{∗}}{2}} [β (∂_{z^{∗}}π)_{i,j,k− \frac{Δz^{∗}}{2}}^{τ+Δτ} + (1 − β) (∂_{z^{∗}}π)_{i,j,k− \frac{Δz^{∗}}{2}}^{τ} ]

−Adv.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + Buoy.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + (D_{w})_{i,j,k− \frac{Δz^{∗}}{2}}^{t−Δt} ,  (D.12)

圧力方程式は

\frac{π_{i,j,k}^{τ+Δτ} − π_{i,j,k}^{τ}}{Δτ} + \frac{1}{G^{1/2}} \frac{c̄^{2}}{c_{p}ρ̄θ̄_{v}^{2}} \{β (∂_{z^{∗}} ρ̄θ̄_{v}w)_{i,j,k}^{τ+Δτ} + (1 − β) (∂_{z^{∗}} ρ̄θ̄_{v}w)_{i,j,k}^{τ}\}
= −\frac{c̄^{2}}{c_{p}θ̄_{v}} [(∂_{x}u)_{i,j,k} + (∂_{y}v)_{i,j,k} + \overline{G^{13}}_{i,j,k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{i,j,k}
+\overline{G^{23}}_{i,j,k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{i,j,k}]^{τ+Δτ}(D.13)

となる.

長い時間間隔 t で計算する項について離散化したものを書き下す. 移流項は

Adv.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} = (\overline{u}^{xz^{∗}}∂_{x}\overline{w}^{x})_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + (\overline{v}^{yz^{∗}}∂_{y}\overline{w}^{y})_{i,j,k− \frac{Δz^{∗}}{2}}^{t}
+ (ω∂_{z^{∗}} \overline{w}^{z^{∗}})_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + Diff.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} .  (D.14)

数値粘性項は

Diff.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} = ν_{2h} [∂_{x}^{2}w + ∂_{y}^{2}w]_{i,j,k− \frac{Δz^{∗}}{2}}^{t} + ν_{2v} [∂_{z^{∗}}^{2}w]_{i,j,k− \frac{Δz^{∗}}{2}}^{t} . (D.15)

浮力項は

Buoy.w_{i,j,k− \frac{Δz^{∗}}{2}}^{t} = g [\overline{\frac{θ}{θ̄}}^{z^{∗}} − 1 + 0.61 (\overline{q_{v}}^{z^{∗}} − \overline{q̄_{v}}^{z^{∗}}) − \overline{q_{c}}^{z^{∗}} − \overline{q_{r}}^{z^{∗}} ]_{i,j,k− \frac{Δz^{∗}}{2}}^{t} . (D.16)

拡散項 (右辺の下付添字と上付添字 t − Δt は省略) は

D_{i,j,k− \frac{Δz^{∗}}{2}} = ∂_{x} \{\overline{K_{m}}^{xz^{∗}} (∂_{x}w) + \begin{matrix}\overline{K_{m}} & xz^{∗} \\ \overline{G^{1/2}} & xz^{∗} \end{matrix} (∂_{z^{∗}} \overline{u}^{xz^{∗}}) + \overline{K_{m}}^{xz^{∗}}G^{13} (∂_{z^{∗}} \overline{w}^{xz^{∗}})\}

+∂_{y} \{\overline{K_{m}}^{yz^{∗}} (∂_{y}w) + \begin{matrix}\overline{K_{m}} & yz^{∗} \\ \overline{G^{1/2}} & yz^{∗} \end{matrix} (∂_{z^{∗}}v) + \overline{K_{m}}^{yz^{∗}}G^{23} (∂_{z^{∗}} \overline{w}^{yz^{∗}})\}

+\overline{G^{13}}^{x}∂_{z^{∗}} \{K_{m} (∂_{x}\overline{w}^{xz^{∗}}) + \frac{K_{m}}{G^{1/2}} (∂_{z^{∗}} \overline{u}^{xz^{∗}}) + K_{m}\overline{G^{13}}^{xz^{∗}} (∂_{z^{∗}}w)\}

+\overline{G^{23}}^{y}∂_{z^{∗}} \{K_{m} (∂_{y}\overline{w}^{yz^{∗}}) + \frac{K_{m}}{G^{1/2}} (∂_{z^{∗}} \overline{v}^{yz^{∗}}) + K_{m}\overline{G^{23}}^{yz^{∗}} (∂_{z^{∗}}w)\}

+\frac{2}{\overline{G^{1/2}}^{z^{∗}}} ∂_{z^{∗}} \{\frac{K_{m}}{G^{1/2}} (∂_{z^{∗}}w)\}

−\frac{2}{3C_{m}^{2}l^{2}} \{(∂_{x}\overline{K_{m}^{2}}^{xz^{∗}}) + \overline{G^{13}}^{x} (∂_{z^{∗}}K_{m}^{2})

+ (∂_{y}\overline{K_{m}^{2}}^{yz^{∗}}) + \overline{G^{23}}^{y} (∂_{z^{∗}}K_{m}^{2}) + \frac{1}{\overline{G^{1/2}}^{z^{∗}}} (∂_{z^{∗}}K_{m}^{2})\}. (D.17)

式 (D.12) と (D.13) から直接 w^{τ+Δτ}, π^{τ+Δτ} を求めることはできない. これらの式を連立して w^{τ+Δτ} を消去し, π^{τ+Δτ} を求める. 次に式 (D.12) から w^{τ+Δτ} を求める. 連立の方法は乾燥大気の場合と同様なので省略する. 圧力に関する行列は以下のようになる:

( \begin{matrix}A_{0} & B_{0} & & 0 \\ C_{0} & ... & ... \\ & ... & ... & B_{km−1} \\ 0 & & C_{km−1} & A_{km} \end{matrix} ) ( \begin{matrix}π_{0} & π_{0} & ⋯ & π_{0} \\ π_{1} & ... & ... & ⋮ \\ ⋮ & & ... & ⋮ \\ π_{km} & ⋯ & ⋯ & π_{km} \end{matrix} )^{τ+Δτ}

= ( \begin{matrix}D_{0} & D_{0} & ⋯ & D_{0} \\ D_{1} & ... & ... & ⋮ \\ ⋮ & & ... & ⋮ \\ D_{km} & ⋯ & ⋯ & D_{km} \end{matrix} )^{τ} .  (D.18)

この連立 1 次方程式の係数行列は 3 項行列なので, LU 分解を行って計算すると大幅に計算量を減らすことができる. こうして求められた π^{τ+Δτ} と鉛直方向の運動方程式 (D.12) から w^{τ+Δτ} を求める. A_{n}, B_{n}, C_{n}(n = 0, 1, ⋯, n) は以下のようである:

A_{k} = 1 + β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄_{v}^{2}} )_{k} \frac{1}{Δz^{∗2}} \{(\frac{ρ̄\overline{θ}_{v}^{2}}{G^{1/2}} )_{w(k+1)} + (\frac{ρ̄\overline{θ}_{v}^{2}}{G^{1/2}} )_{w(k)}\} ,(D.19)
(k = 1, 3, ⋯, km − 1)  (D.20)

\underline{付録 D 湿潤大気における準圧縮モデルの離散化 78}

A_{0} = 1 + β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄_{v}^{2}} )_{0} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}_{v}^{2}}{G^{1/2}} )_{w(1)} ,  (D.21)
A_{km} = 1 + β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄_{v}^{2}} )_{km} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}_{v}^{2}}{G^{1/2}} )_{w(km)} ,  (D.22)
B_{k} = −β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄_{v}^{2}} )_{k} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}_{v}^{2}}{G^{1/2}} )_{w(k+1)} ,  (D.23)
C_{k} = −β^{2} (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄_{v}^{2}} )_{k+1} \frac{1}{Δz^{∗2}} (\frac{ρ̄\overline{θ}_{v}^{2}}{G^{1/2}} )_{w(k+1)} ,  (D.24)
D_{k} = π_{k}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄_{v}^{2}} )_{k} \{∂_{z^{∗}} (ρ̄θ̄_{v}w)\}_{k}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄_{v}} [(∂_{x}u)_{k} + (∂_{y}v)_{k} + \overline{G^{13}}_{k}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{k}

+\overline{G^{23}}_{k}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{k}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄_{v}^{2}} )_{k}

×[∂_{z^{∗}} ρ̄θ̄_{v}\{w^{τ} − (\frac{c_{p}θ̄_{v}Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{k},

(k = 1, 3, ⋯, km − 3),  (D.25)

\underline{付録 D 湿潤大気における準圧縮モデルの離散化 79}

D_{0} = π_{0}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄_{v}^{2}} )_{0} \{∂_{z^{∗}} (ρ̄θ̄_{v}w)\}_{0}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄_{v}} [(∂_{x}u)_{0} + (∂_{y}v)_{0}

+\overline{G^{13}}_{0}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{0} + \overline{G^{23}}_{1}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{0}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄_{v}^{2}} )_{0}

×[∂_{z^{∗}} ρ̄θ̄_{v}\{w^{τ} − (\frac{c_{p}θ̄_{v}Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{0}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄_{v}^{2}} )_{0} \frac{1}{Δz^{∗}} (\frac{c_{p}θ̄_{v}^{2}}{G^{1/2}} )_{w(0)}

×\{α (∂_{z^{∗}}Div)_{w(0)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(0)}^{τ} + (\frac{G^{1/2}}{c_{p}θ̄_{v}} F.w^{t})_{w(0)}\}, (D.26)
D_{km−2} = π_{km−2}^{τ} − (1 − β) (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄_{v}^{2}} )_{km−2} \{∂_{z^{∗}} (ρ̄θ̄_{v}w)\}_{km−2}^{τ}

−\frac{c̄^{2}Δτ}{c_{p}θ̄_{v}} [(∂_{x}u)_{km−2} + (∂_{y}v)_{km−2}

+\overline{G^{13}}_{km}^{xz^{∗}} (∂_{z^{∗}} \overline{u}^{xz^{∗}})_{km−2} + \overline{G^{23}}_{km}^{yz^{∗}} (∂_{z^{∗}} \overline{v}^{yz^{∗}})_{km−2}]^{τ+Δτ}

−β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ}{c_{p}ρ̄θ̄_{v}^{2}} )_{km−2}

×[∂_{z^{∗}} ρ̄θ̄_{v}\{w^{τ} − (\frac{c_{p}θ̄_{v}Δτ}{\overline{G^{1/2}}^{z^{∗}}} ) \{−α (∂_{z^{∗}}Div)^{τ} + (1 − β) (∂_{z^{∗}}π)^{τ}\}

−F.w^{t}Δτ }]_{km−2}

+β (\frac{1}{G^{1/2}} \frac{c̄^{2}Δτ^{2}}{ρ̄θ̄_{v}^{2}} )_{km−2} \frac{1}{Δz^{∗}} (\frac{c_{p}θ̄_{v}^{2}}{G^{1/2}} )_{w(km−1)}

×{α (∂_{z^{∗}}Div)_{w(km−1)}^{τ} − (1 − β) (∂_{z^{∗}}π)_{w(km−1)}^{τ}

+ (\frac{G^{1/2}}{c_{p}θ̄_{v}} F.w^{t})_{w(km−1)}\}  (D.27)

D.3.3 熱力学の式と混合比の方程式

熱力学の式 (B.22) と水蒸気・雲水および雨水の混合比の方程式 (B.23)〜(B.25) を中央差分で離散化する. 時間間隔は長い時間間隔 Δt である. 本節では変数の下付添字 i, j, k を省略する.

熱力学の式は

\frac{θ^{t+Δt} − θ^{t−Δt}}{2Δt} = −Adv.θ^{t} + M_{θ}^{t} + D_{θ}^{t−Δt}.  (D.28)

水蒸気混合比の方程式は

\frac{q_{v}^{t+Δt} − q_{v}^{t−Δt}}{2Δt} = −Adv.q_{v}^{t} + M_{q_{v}}^{t} + D_{q_{v}}^{t−Δt}.  (D.29)

本モデルでは過飽和がない場合を考える. 求めた水蒸気の混合比と飽和混合比の値が q_{v}^{t+Δt} > q_{vs}^{t+Δt} となった場合, 湿潤飽和調節法 (D.4 参照) を用いる.

雲水混合比の方程式は

\frac{q_{c}^{t+Δt} − q_{c}^{t−Δt}}{2Δt} = −Adv.q_{c}^{t} + M_{q_{c}}^{t} + D_{q_{c}}^{t−Δt},  (D.30)

雨水混合比の方程式は

\frac{q_{r}^{t+Δt} − q_{r}^{t−Δt}}{2Δt} = −Adv.q_{r}^{t} + M_{q_{r}}^{t} + D_{q_{r}}^{t−Δt}  (D.31)

となる.

これまでに求めた u^{τ+Δτ}, v^{τ+Δτ}, w^{τ+Δτ} および π^{τ+Δτ} を用いてθ^{τ+Δτ}, q_{v}^{τ+Δτ}, q_{c}^{τ+Δτ}, q_{r}^{τ+Δτ} を順に求める. ただし, 計算の際, 長い時間間隔における中央差分での計算モードを抑えるため, Asselin (1972) の時間フィルタ (2.3.4 節参照) を用いる.

上式 (D.28)〜 (D.31) 中の各項はそれぞれ以下のようになる; 移流項は

Adv.θ^{t} = [\overline{u∂_{x}θ}^{x} + \overline{v∂_{y}θ}^{y} + \overline{ω∂_{z^{∗}}θ}^{z^{∗}}]^{t} + Diff.θ^{t},  (D.32)
Adv.q_{v}^{t} = [\overline{u∂_{x}q_{v}}^{x} + \overline{v∂_{y}q_{v}}^{y} + \overline{ω∂_{z^{∗}}q_{v}}^{z^{∗}}]^{t} + Diff.q_{v}^{t},  (D.33)
Adv.q_{c}^{t} = [\overline{u∂_{x}q_{c}}^{x} + \overline{v∂_{y}q_{c}}^{y} + \overline{ω∂_{z^{∗}}q_{c}}^{z^{∗}}]^{t} + Diff.q_{c}^{t},  (D.34)
Adv.q_{r}^{t} = [\overline{u∂_{x}q_{r}}^{x} + \overline{v∂_{y}q_{r}}^{y} + \overline{ω∂_{z^{∗}}q_{r}}^{z^{∗}}]^{t} + Diff.q_{r}^{t}.  (D.35)

数値粘性項は

Diff.θ^{t} = ν_{2h} [∂_{x}^{2}θ + ∂_{y}^{2}θ]^{t} + ν_{2v} [∂_{z^{∗}}^{2}θ]^{t} ,  (D.36)
Diff.q_{v}^{t} = ν_{2h} [∂_{x}^{2}q_{v} + ∂_{y}^{2}q_{v}]^{t} + ν_{2v} [∂_{z^{∗}}^{2}q_{v}]^{t} ,  (D.37)
Diff.q_{c}^{t} = ν_{2h} [∂_{x}^{2}q_{c} + ∂_{y}^{2}q_{c}]^{t} + ν_{2v} [∂_{z^{∗}}^{2}q_{c}]^{t} ,  (D.38)
Diff.q_{r}^{t} = ν_{2h} [∂_{x}^{2}q_{r} + ∂_{y}^{2}q_{r}]^{t} + ν_{2v} [∂_{z^{∗}}^{2}q_{r}]^{t} .  (D.39)

生成・消滅項は

M_{θ}^{t} = −γ (\frac{q_{vs}^{t+Δt} − q_{vs}^{t−Δt}}{2Δt} + Adv.q_{vs}^{t} + E_{r}^{t}) ,  (D.40)
M_{q_{v}}^{t} = \frac{q_{vs}^{t+Δt} − q_{vs}^{t−Δt}}{2Δt} + Adv.q_{vs}^{t−Δt} + E_{r}^{t},  (D.41)
M_{q_{c}}^{t} = −\frac{q_{vs}^{t+Δt} − q_{vs}^{t−Δt}}{2Δt} − Adv.q_{vs}^{t} − A_{r}^{t} − C_{r}^{t},  (D.42)
M_{q_{r}}^{t} = F_{r}^{t} − E_{r}^{t} + A_{r}^{t} + C_{r}^{t}.  (D.43)

また,

Adv.q_{vs}^{t} = [\overline{u∂_{x}q_{vs}}^{x} + \overline{v∂_{y}q_{vs}}^{y} + \overline{ω∂_{z^{∗}}q_{vs}}^{z^{∗}}]^{t} + Diff.q_{vs}^{t},  (D.44)
q_{vs} = \frac{3.8}{p̄} exp (17.27\frac{Π̄θ − 273}{Π̄θ − 36})
≃ q_{vs}^{∗} [1 + \frac{4093aΠ̄}{(Π̄θ^{∗} − 36)^{2}} (θ^{t+Δt} − θ^{∗})] ,  (D.45)

F_{r}^{t} = [\frac{1}{G^{1/2}ρ̄}∂_{z^{∗}} (ρ̄V q_{r})]^{t} ,  (D.46)
V = 36.34 (ρ̄q_{r})^{0.1346} (\frac{ρ̄}{ρ_{0}} )^{− \frac{1}{2}}  (D.47)

である. ただし q_{vs}^{∗} は θ = θ^{∗} のときの q_{vs} の値である. 式 (D.45) の近似は Tetens の式より θ = θ^{∗} について線形化を行った.

拡散項は以下のようになる:

D_{θ}^{t−Δt} = 3∂_{x} \{\overline{K_{m}}^{x} (∂_{x}θ) + \overline{K_{m}}^{x}\overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{θ}^{xz^{∗}})\}^{t−Δt}

+3∂_{y} \{\overline{K_{m}}^{y} (∂_{y}θ) + \overline{K_{m}}^{y}\overline{G^{23}}^{z^{∗}} (∂_{z^{∗}} \overline{θ}^{yz^{∗}})\}^{t−Δt}

+3\overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{x}\overline{θ}^{xz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{13}}^{x} (∂_{z^{∗}}θ)\}^{t−Δt}

+3\overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{y}\overline{θ}^{yz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{23}}^{y} (∂_{z^{∗}}θ)\}^{t−Δt}

+\frac{3}{G^{1/2}} ∂_{z^{∗}} \{\frac{K_{m}}{\overline{G^{1/2}}^{z^{∗}}} (∂_{z^{∗}}θ)\}^{t−Δt} ,  (D.48)
D_{q_{v}}^{t−Δt} = 3∂_{x} \{\overline{K_{m}}^{x} (∂_{x}q_{v}) + \overline{K_{m}}^{x}\overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{q_{v}}^{xz^{∗}})\}^{t−Δt}

+3∂_{y} \{\overline{K_{m}}^{y} (∂_{y}q_{v}) + \overline{K_{m}}^{y}\overline{G^{23}}^{z^{∗}} (∂_{z^{∗}} \overline{q_{v}}^{yz^{∗}})\}^{t−Δt}

+3\overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{x}\overline{q_{v}}^{xz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{13}}^{x} (∂_{z^{∗}}q_{v})\}^{t−Δt}

+3\overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{y}\overline{q_{v}}^{yz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{23}}^{y} (∂_{z^{∗}}q_{v})\}^{t−Δt}

+\frac{3}{G^{1/2}} ∂_{z^{∗}} \{\frac{K_{m}}{\overline{G^{1/2}}^{z^{∗}}} (∂_{z^{∗}}q_{v})\}^{t−Δt} ,  (D.49)
D_{q_{c}}^{t−Δt} = 3∂_{x} \{\overline{K_{m}}^{x} (∂_{x}q_{c}) + \overline{K_{m}}^{x}\overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{q_{c}}^{xz^{∗}})\}^{t−Δt}

+3∂_{y} \{\overline{K_{m}}^{y} (∂_{y}q_{c}) + \overline{K_{m}}^{y}\overline{G^{23}}^{z^{∗}} (∂_{z^{∗}} \overline{q_{v}}^{yz^{∗}})\}^{t−Δt}

+3\overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{x}\overline{q_{c}}^{xz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{13}}^{x} (∂_{z^{∗}}q_{c})\}^{t−Δt}

+3\overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{y}\overline{q_{c}}^{yz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{23}}^{y} (∂_{z^{∗}}q_{c})\}^{t−Δt}

+\frac{3}{G^{1/2}} ∂_{z^{∗}} \{\frac{K_{m}}{\overline{G^{1/2}}^{z^{∗}}} (∂_{z^{∗}}q_{c})\}^{t−Δt} ,  (D.50)
D_{q_{r}}^{t−Δt} = 3∂_{x} \{\overline{K_{m}}^{x} (∂_{x}q_{r}) + \overline{K_{m}}^{x}\overline{G^{13}}^{z^{∗}} (∂_{z^{∗}} \overline{q_{r}}^{xz^{∗}})\}^{t−Δt}

+3∂_{y} \{\overline{K_{m}}^{y} (∂_{y}q_{r}) + \overline{K_{m}}^{y}\overline{G^{23}}^{z^{∗}} (∂_{z^{∗}} \overline{q_{r}}^{yz^{∗}})\}^{t−Δt}

+3\overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{x}\overline{q_{r}}^{xz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{13}}^{x} (∂_{z^{∗}}q_{r})\}^{t−Δt}

+3\overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \{\overline{K_{m}}^{z^{∗}} (∂_{y}\overline{q_{r}}^{yz^{∗}}) + \overline{K_{m}}^{z^{∗}} \overline{G^{23}}^{y} (∂_{z^{∗}}q_{r})\}^{t−Δt}

+\frac{3}{G^{1/2}} ∂_{z^{∗}} \{\frac{K_{m}}{\overline{G^{1/2}}^{z^{∗}}} (∂_{z^{∗}}q_{r})\}^{t−Δt} .  (D.51)

D.3.4 サブグリッドスケールの乱流運動エネルギー方程式

サブグリッドスケールの乱流運動エネルギー方程式を運動量渦拡散係数で書き換えた式 (B.40) は長い時間間隔 t で計算される. ここでは変数の下付添字 i, j, k は省略する.

\frac{K_{m}^{t+Δt} − K_{m}^{t−Δt}}{2Δt} = −Adv.K_{m}^{t} + \frac{3gC_{m}^{2}l^{2}}{2G^{1/2}} (−A∂_{z^{∗}} \overline{θ_{e}}^{z^{∗}} + ∂_{z^{∗}} \overline{q_{l}}^{z^{∗}})^{t}
+C_{m}^{2}l^{2}\{ (∂_{x}u + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \overline{u}^{xz^{∗}})^{2}

+ (∂_{y}v + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \overline{v}^{yz^{∗}})^{2} + (\frac{1}{G^{1/2}} ∂_{z^{∗}}w)^{2}\}^{t}

+\frac{C_{m}^{2}l^{2}}{2} \{(∂_{x}\overline{v}^{xy} + ∂_{y}\overline{u}^{xy} + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \overline{v}^{yz^{∗}} + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \overline{u}^{xz^{∗}})^{2}

+ (∂_{x}w + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}}w + \frac{1}{G^{1/2}} ∂_{z^{∗}} \overline{u}^{xz^{∗}})^{2}

+ (∂_{y}\overline{w}^{yz^{∗}} + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}}w + \frac{1}{G^{1/2}} ∂_{z^{∗}} \overline{v}^{yz^{∗}})^{2}\}^{t}

+\frac{1}{2}\{∂_{x}^{2}K_{m}^{2} + ∂_{x} (\overline{G^{13}}^{z^{∗}}∂_{z^{∗}} \overline{K_{m}^{2}}^{xz^{∗}})

+∂_{z^{∗}} (\overline{G^{13}}^{x}∂_{x}\overline{K_{m}^{2}}^{xz^{∗}}) + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} (\overline{G^{13}}^{x}∂_{z^{∗}}K_{m}^{2})

+∂_{y}^{2}K_{m}^{2} + ∂_{y} (\overline{G^{23}}^{z^{∗}}∂_{z^{∗}} \overline{K_{m}^{2}}^{yz^{∗}})

+∂_{z^{∗}} (\overline{G^{23}}^{y}∂_{y}\overline{K_{m}^{2}}^{yz^{∗}}) + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} (\overline{G^{23}}^{y}∂_{z^{∗}}K_{m}^{2})

+\frac{1}{G^{1/2}} ∂_{z^{∗}} (\overline{\frac{1}{G^{1/2}}}^{∗}∂_{z^{∗}}K_{m}^{2})\}^{t}

+\{(∂_{x}\overline{K_{m}}^{x} + \overline{G^{13}}^{xz^{∗}}∂_{z^{∗}} \overline{K_{m}}^{z^{∗}})^{2}
+ (∂_{y}\overline{K_{m}}^{y} + \overline{G^{23}}^{yz^{∗}}∂_{z^{∗}} \overline{K_{m}}^{z^{∗}})^{2}

+ (\frac{1}{G^{1/2}} ∂_{z^{∗}} \overline{K_{m}}^{z^{∗}})^{2}\}^{t}

−\frac{C_{ε}}{2C_{m}l^{2}} K_{m}^{2}  (D.52)

となる. 移流項と数値粘性項はそれぞれ以下のようになる:

Adv.K_{m}^{t} = [\overline{u∂_{x}K_{m}}^{x} + \overline{v∂_{y}K_{m}}^{y} + \overline{ω∂_{z^{∗}}K_{m}}^{z^{∗}}]^{t} + Diff.K_{m}^{t},  (D.53)
Diff.K_{m}^{t} = ν_{2h} [∂_{x}^{2}K_{m} + ∂_{y}^{2}K_{m}]^{t} + ν_{2v} [∂_{z^{∗}}^{2}K_{m}]^{t} .  (D.54)

D.4 湿潤飽和調節法

本モデルでは過飽和なしの場合を考える. 過飽和しないように温位と混合比を計算するために, 湿潤飽和調節法 (Soong and Ogura, 1973) を用いる. 以下でその方法について述べる.

先に示したように, 温位, 水蒸気と雲水の混合比の方程式は

\frac{dϕ}{dt} = M_{ϕ} + D_{ϕ}  (D.55)

と書ける. ただし ϕ = θ, q_{v}, q_{c} であり, D_{ϕ} は拡散項, 微物理過程の生成・消滅項は

M_{θ} = −γ \{\frac{dq_{vs}}{dt} + E_{r}\} ,  (D.56)
M_{q_{v}} = \frac{dq_{vs}}{dt} + E_{r},  (D.57)
M_{q_{c}} = −\frac{dq_{vs}}{dt} − A_{r} − C_{r},  (D.58)
q_{vs} = \frac{3.8}{p̄} exp (17.27\frac{Π̄θ − 273}{Π̄θ − 36})  (D.59)

である.

ここで湿潤調節法に必要な方程式を求める. (温位式)+γ×(水蒸気混合比の式) より

\frac{dθ}{dt} + γ \frac{dq_{v}}{dt} = D_{θ} + D_{q_{v}},  (D.60)

(水蒸気混合比の式)+(雲水混合比の式) より

\frac{dq_{v}}{dt} + \frac{dq_{c}}{dt} = D_{q_{v}} + D_{q_{c}} + E_{r} − A_{r} − C_{r}  (D.61)

である.

予報方程式 (D.55) から計算された t + Δt での各値を仮に θ^{∗}, q_{v}^{∗}, q_{c}^{∗} と表す. ここでは過飽和しないという条件で考えているので, q_{v}^{∗} > q_{vs}^{∗} ならば湿潤飽和調節を行って t + Δt での各値を求める. q_{v}^{∗} ≤ q_{vs}^{∗} ならば予報方程式から得られた値を用いる. 以下で湿潤飽和調節法を解説する.

式 (D.60) から

\frac{θ^{t+Δt} − θ^{t−Δt}}{2Δt} + γ \frac{q_{v}^{t+Δt} − q_{v}^{t−Δt}}{2Δt} = D_{θ}^{t−Δt} + D_{q_{v}}^{t−Δt},  (D.62)
\frac{θ^{∗} − θ^{t−Δt}}{2Δt} + γ \frac{q_{v}^{∗} − q_{v}^{t−Δt}}{2Δt} = D_{θ}^{t−Δt} + D_{q_{v}}^{t−Δt}.  (D.63)

\underline{付録 D 湿潤大気における準圧縮モデルの離散化 85}

辺々引くと

θ^{t+Δt} + γq_{v}^{t+Δt} = θ^{∗} + γq_{v}^{∗}  (D.64)

が成り立つ. 過飽和しないので, 飽和に達した後の t + Δt における水蒸気の混合比と飽和混合比は等しくなければならない:

q_{v}^{t+Δt} = q_{vs}^{t+Δt}
= b exp [\frac{a (Π̄θ^{t+Δt} − 273)}{Π̄θ^{t+Δt} − 36}]
≈ q_{vs}^{∗} \{1 + \frac{273aΠ̄}{(Π̄θ^{∗} − 36)^{2}} (θ^{t+Δt} − θ^{∗})\} .  (D.65)

ただし a = 7.5 ln 10, b = 3.8/p̄ である. したがって, 式 (D.64) は

θ^{t+Δt} = θ^{∗} + γ (q_{v}^{∗} − q_{vs}^{t+Δt})
= θ^{∗} + γr_{1} (q_{v}^{∗} − q_{vs}^{∗})  (D.66)

となる. ただし

r_{1} = [1 + γq_{vs}^{∗} \frac{273aΠ̄}{(Π̄θ^{∗} − 36)^{2}} ]^{−1} .  (D.67)

式 (D.64) と (D.66) から

q_{v}^{t+Δt} = q_{v}^{∗} + \frac{1}{γ} (θ^{∗} − θ^{t+Δt})  (D.68)

また, (D.64) を求めたのと同様に, (D.61) より

q_{v}^{t+Δt} + q_{c}^{t+Δt} = q_{v}^{∗} + q_{c}^{∗}  (D.69)

が得られる. 式 (D.66), (D.68), (D.69) より暫定的に θ, q_{v}, q_{c} の値を求める. q_{c}^{t+Δt} > 0, つまり雲水ができる場合, 式 (D.66), (D.68), (D.69) を繰り返し用いて値が収束するまで計算する^{1} q_{c}^{t+Δt} ≤ 0, つまり雲水ができない場合,

1凝結して周りの温度が上がるため, 蒸発が起こるのを繰り返す. この状態が落ち着くまで計算する.

q_{c}^{t+Δt} = 0,  (D.70)
q_{v}^{t+Δt} = q_{v}^{∗} + q_{c}^{∗},  (D.71)
θ^{t+Δt} = θ^{∗} − γq_{c}^{∗}  (D.72)

を用いる. ここで式 (D.71) は式 (D.69) から, 式 (D.72) は (D.68) と (D.71) から求められる.

参考文献

Asselin, R., 1972: Frequency filter for time integrations. Mon. Weather Rev., 100, 487-490.

Ćurić, M., D. Janc, D. Vujović, V. Vučković, 2003: The effects fo a river valley on an isolated cumulonimbus cloud development. Atmos. Res., 66, 123-139.

Ćurić, M., Agee, E.M., Asai, T. (Eds.), 1982: The development of the cumulonimbus clouds which moves along a valley. Clouds Dynamics. Reidel, Dordrecht, 259-272.

Deardorff, J. W., 1972: Numerical investigation of neutral and unstable planetary boundary layers. J. Atmos. Sci., 29, 91-115.

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

Houze, R. A., 1993: Cloud dynamics. Academic Press.

Ikawa, M. and K. Saito, 1991: Description of a nonhydrostatic model developed at the Forecast Research Department of the MRI. Technical Reports of the MRI, 28, 238pp.

Kessler, E., 1969: On the Distribution and Continuity of Water Substance in Atmospheric Circuration. Meteor. Monogr., Amer. Meteor.Soc., 32, 84 pp.

Klemp J. B. and Robert B. Wilhelmson, 1978: The simulation of three-dimensional convective storm dynamics. J. Atmos. Sci., 35, 1070-1096.

Lee, B.D. and R.B. Wilhelmson, 1997: The numerical simulation of nonsupercell tornadogenesis.

Mellor, G., and T. Yamada, 1974: A hierarchy of turbulence closure models for atmospheric boundary layers. J. Atmos. Sci., 31, 1791-1806.

中野満寿男, 2003: 鉛直シアーが台風初期渦形成に及ぼす影響. 九州大学大学院理学府地球惑星科学専攻修士論文.

小倉義光, 1997: メソ気象の基礎理論, 東京大学出版会, 215pp.

Ogura, Y. and N. A. Phillips, 1962; Scale analysis of deep and shallow water convection in the atmosphere. J. Atmos. Sci., 19, 173-179.

Ogura, Y., and T. Takahashi, 1971: Numerical simulation of the life cycle of a thunderstorm cell. Mon. Wea. Rev., 99, 895-911.

Ogura, Y. and M. Yoshizaki, 1988: Numerical study of orographic-convective precipitation over the eastern Arabic Sea and the Ghats Mountains during the summer monsoon. J. Atoms. Sci., 45, 2097-2122.

大宮司久明, 三宅裕, 吉澤徴, 1998: 乱流の数値流体力学—モデルと計算法. 東京大学出版会, 652pp.

斉藤和雄編, 1999: 非静力学モデル. 気象研究ノート第 196 号, 日本気象学会,195pp.

Skamarock, Wiliam C. and Joseph B. Klemp, 1992: The stability of time-split methods for the hydrostatic and the nonhydrostatic elastic equation. Mon. Wea. Rev., 120, 2109-2127.

Soong, S-T., and Y. Ogura, 1973: A comparison between axi-symmetric and slab-symmetric cumulus cloud models. J. Atmos. Sci., 30, 879-893

Takahashi, T. and T. Kawano, 1998: Numerical sensitivity study of rainband precipitation and evolution. J. Atmos. Soc., 55, 57-87.

豊田英司, 石渡正樹, 竹広真一, 林祥介, gtool4 開発グループ, 2002: gtool4 Fortran90 ツール・ライブラリ, http://www.gfd-dennou.org/arch/gtool4/, 地球流体電脳倶楽部.

坪木和久, 榊原篤志, 2001: CReSS ユーザーズガイド第 2 版.

http://www.tokyo.rist.or.jp/CReSS Fujin/CReSS.top.html

Wilhelmson, R. B., 1977: On the thermodyanmics equation for deep convetion. Mon. Wea. Rev., 105, 545-547.

Wilhelmson, R., and Y. Ogura, 1972: The pressure perturbation and the numerical modeling of a cloud. J. Atmos. Sci., 29, 1295-1307.

Wilhelmson, R. B., and Y. Ogura, 1977: The pressure perturbation and the numerical modleing of a cloud. J. Atmos. Sci., 29, 1295-1307.

Xue, M., K. K. Droegemeier, V. Wong, A. Shapiro and K. Brewster, 1995: Advanced Regional Prediction System, Version 4.0. Center for Analysis and Prediction of Storms, University of Oklahoma, 380pp.