主成分凝結大気における 2 次元非静力学モデルの離散化

山下達也

YAMASHITA Tatsuya 北海道大学理学院宇宙理学専攻

2011 年 12 月 14 日

\underline{2 大気主成分凝結対流離散化}

目次

第 1 章離散化 3 1.1 離散化の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 格子の配置 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 空間方向の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.3.1 平均操作 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.2 空間微分の離散化 . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.3 準圧縮方程式系の空間方向の離散化 . . . . . . . . . . . . . 7

1.4 時間方向の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.1 モード別時間分割法 . . . . . . . . . . . . . . . . . . . . . . 9 1.4.2 音波減衰項 . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4.3 数値粘性項 . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.4.4 準圧縮方程式系の時間方向の離散化 . . . . . . . . . . . . . 11 1.4.5 Robert, Asselin の時間フィルター . . . . . . . . . . . . . . . 15

参考文献 17

付録 A: 負の雲密度の発生を防ぐ処置に関して 19

\underline{大気主成分凝結対流離散化 3}

第1章離散化

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

1.1 離散化の概要

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

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

1.2 格子の配置

本モデルでは水平方向, 鉛直方向の格子点配置方法として Arakawa C グリッド, Lorenz グリッドをそれぞれ採用している. Arakawa C グリッドとは水平方向のベクトル量とスカラー量を半格子ずらして配置する格子点配置方法のことを言う. Arakawa C グリッドは重力波を表現するのに適しているとされている (Arakawa and Lamb, 1977). Lorenz グリッドとは鉛直方向のベクトル量とスカラー量を半格子ずらして配置する格子点配置方法のことを言う. スカラー量の格子点を (i, k), ベクトル量の水平成分に対する格子点を (i(u), k), ベクトル量の鉛直成分に対する格子点を (i, k(w)), 格子の角に相当する点を (i(u), k(w)) のように表すことにすると, 格子点の配置は図 1.1 のように表される.

図 1.1: 格子点の配置. 杉山他 (2006) より引用した.

1.3 空間方向の離散化

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

1.3.1 平均操作

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

以下, 計算に必要な平均操作を示す. 但し ϕ, H, V はそれぞれスカラー量, ベクトル量の水平成分, ベクトル量の鉛直成分を表す. また下付き添字は格子点位置を表している.

ϕ_{i(u),k} = \frac{ϕ_{i+1,k} + ϕ_{i,k}}{2} ,  (1.1)
ϕ_{i,k(w)} = \frac{ϕ_{i,k+1} + ϕ_{i,k}}{2} ,  (1.2)
ϕ_{i(u),k(w)} = \frac{ϕ_{i,k} + ϕ_{i+1,k} + ϕ_{i,k+1} + ϕ_{i+1,k+1}}{4} ,  (1.3)
H_{i,k} = \frac{H_{i(u),k} + H_{i−1(u),k}}{2} ,  (1.4)
H_{i,k(w)} = \frac{H_{i(u),k+1} + H_{i−1(u),k+1} + H_{i(u),k} + u_{i−1(u),k}}{4} ,  (1.5)
H_{i(u),k(w)} = \frac{H_{i(u),k+1} + H_{i(u),k}}{2} ,  (1.6)
V_{i,k} = \frac{V_{i,k(w)} + V_{i,k−1(w)}}{2} ,  (1.7)
V_{i(u),k} = \frac{V_{i+1,k(w)} + V_{i,k(w)} + V_{i+1,k−1(w)}, V_{i,k−1(w)}}{4} ,  (1.8)
V_{i(u),k(w)} = \frac{V_{i+1,k(w)} + V_{i,k(w)}}{2} .  (1.9)

1.3.2 空間微分の離散化

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

\underline{6 大気主成分凝結対流離散化} 以下に 2 次精度の中心差分を用いた微分操作を示す. 但し ψ は格子の角に相当する点で評価している変数を表す.

(\frac{∂ϕ}{∂x})_{i(u),k} = \frac{ϕ_{i+1,k} − ϕ_{i,k}}{Δx} ,  (1.10)
(\frac{∂ϕ}{∂z})_{i,k(w)} = \frac{ϕ_{i,k+1} − ϕ_{i,k}}{Δz} ,  (1.11)
(\frac{∂H}{∂x} )_{i,k} = \frac{H_{i(u),k} − H_{i−1(u),k}}{Δx} ,  (1.12)
(\frac{∂H}{∂z} )_{i(u),k(w)} = \frac{H_{i(u),k+1} − H_{i(u),k}}{Δz} ,  (1.13)
(\frac{∂V}{∂x} )_{i(u),k(w)} = \frac{V_{i+1,k(w)} − V_{i,k(w)}}{Δx} ,  (1.14)
(\frac{∂V}{∂z} )_{i,k} = \frac{V_{i,k(w)} − V_{i,k−1(w)}}{Δx} ,  (1.15)
(\frac{∂ψ}{∂x} )_{i,k(w)} = \frac{ψ_{i(u),k(w)} − ψ_{i−1(u),k(w)}}{Δx} ,  (1.16)
(\frac{∂ψ}{∂z} )_{i(u),k} = \frac{ψ_{i(u),k(w)} − ψ_{i(u),k−1(w)}}{Δz} ,  (1.17)
(\frac{∂ϕ}{∂z})_{i,k} = \frac{ϕ_{i,k(w)} − ϕ_{i,k−1(w)}}{Δz} ,  (1.18)
(\frac{∂(ϕH)}{∂x})_{i,k} = \frac{ϕ_{i(u),k}H_{i(u),k} − ϕ_{i−1(u),k}H_{i−1(u),k}}{Δx} ,  (1.19)
(\frac{∂(ϕV)}{∂z})_{i,k} = \frac{ϕ_{i,k(w)}V_{i,k(w)} − ϕ_{i,k−1(w)}V_{i,k−1(w)}}{Δz} .  (1.20)

以下, 4 次精度の中心差分を用いた微分操作を示す.

(\frac{∂ϕ}{∂x})_{i(u),k} = \frac{9}{8} \frac{ϕ_{i+1,k} − ϕ_{i,k}}{Δx} − \frac{1}{24} \frac{ϕ_{i+2,k} − ϕ_{i−1,k}}{Δz} ,  (1.21)
(\frac{∂ϕ}{∂z})_{i,k(w)} = \frac{9}{8} \frac{ϕ_{i,k+1} − ϕ_{i,k}}{Δx} − \frac{1}{24} \frac{ϕ_{i,k+2} − ϕ_{i,k−1}}{Δz} ,  (1.22)
(\frac{∂H}{∂x} )_{i,k} = \frac{9}{8} \frac{H_{i(u),k} − H_{i−1(u),k}}{Δx} − \frac{1}{24} \frac{H_{i+1(u),k} − H_{i−2(u),k}}{Δz} , (1.23)
(\frac{∂H}{∂z} )_{i(u),k(w)} = \frac{9}{8} \frac{H_{i(u),k+1} − H_{i(u),k}}{Δx} − \frac{1}{24} \frac{H_{i(u),k+2} − H_{i(u),k−1}}{Δz} , (1.24)
(\frac{∂V}{∂x} )_{i(u),k(w)} = \frac{9}{8} \frac{V_{i+1,k(w)} − V_{i,k(w)}}{Δx} − \frac{1}{24} \frac{V_{i+2,k(w)} − V_{i−1,k(w)}}{Δz} , (1.25)

\overline{2011/12/14(山下達也)}

\underline{大気主成分凝結対流離散化 7}

(\frac{∂V}{∂z} )_{i,k} = \frac{9}{8} \frac{V_{i,k(w)} − V_{i,k−1(w)}}{Δx} − \frac{1}{24} \frac{V_{i,k+1(w)} − V_{i,k−2(w)}}{Δz} ,  (1.26)
(\frac{∂ψ}{∂x} )_{i,k(w)} = \frac{9}{8} \frac{ψ_{i(u),k(w)} − ψ_{i−1(u),k(w)}}{Δx} − \frac{1}{24} \frac{ψ_{i+1(u),k(w)} − ψ_{i−2(u),k(w)}}{Δz} ,

(1.27)

(\frac{∂ψ}{∂z} )_{i(u),k} = \frac{9}{8} \frac{ψ_{i(u),k(w)} − ψ_{i(u),k−1(w)}}{Δx} − \frac{1}{24} \frac{ψ_{i(u),k+1(w)} − ψ_{i(u),k−2(w)}}{Δz} ,

(1.28)

(\frac{∂ϕ}{∂z})_{i,k} = \frac{9}{8} \frac{ϕ_{i,k(w)} − ϕ_{i,k−1(w)}}{Δx} − \frac{1}{24} \frac{ϕ_{i,k+1(w)} − ϕ_{i,k−2(w)}}{Δz} ,  (1.29)
(\frac{∂(ϕH)}{∂x})_{i,k} = \frac{9}{8} \frac{ϕ_{i(u),k}H_{i(u),k} − ϕ_{i−1(u),k}H_{i−1(u),k}}{Δx}
−\frac{1}{24} \frac{ϕ_{i+1(u),k}H_{i+1(u),k} − ϕ_{i−2(u),k}H_{i−2(u),k}}{Δz} ,  (1.30)
(\frac{∂(ϕV)}{∂z})_{i,k} = \frac{9}{8} \frac{ϕ_{i,k(w)}V_{i,k(w)} − ϕ_{i,k−1(w)}V_{i,k−1(w)}}{Δx}
−\frac{1}{24} \frac{ϕ_{i,k+1(w)}V_{i,k+1(w)} − ϕ_{i,k−2(w)}V_{i,k−2(w)}}{Δz} .  (1.31)

1.3.3 準圧縮方程式系の空間方向の離散化

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

\frac{∂u_{i(u),k}}{∂t} = −u_{i(u),k} (\frac{∂u}{∂x})_{i(u),k} − w_{i(u),k} (\frac{∂u}{∂z})_{i(u),k}
−c_{pv}\overline{θ}_{i(u),k} (\frac{∂Π}{∂x})_{i(u),k} + [D_{u}]_{i(u),k} .  (1.32)

鉛直方向の運動方程式

\frac{∂w_{i,k(w)}}{∂t} = −u_{i,k(w)} (\frac{∂w}{∂x})_{i,k(w)} − w_{i,k(w)} (\frac{∂w}{∂z})_{i,k(w)}

−c_{pv}\overline{θ}_{i,k(w)} (\frac{∂Π}{∂z})_{i,k(w)} + [D_{w}]_{i,k(w)} + g (\frac{θ′}{\overline{θ}} − \frac{R}{p_{0}} \frac{\overline{θ}}{\overline{Π}^{c_{vv}/R}} ρ_{s})_{i,k(w)} . (1.33)

\underline{8 大気主成分凝結対流離散化} 圧力方程式

\frac{∂Π_{i,k}}{∂t} = −\frac{\overline{c}_{i,k}^{2}}{c_{pv}\overline{ρ}_{i,k}\overline{θ}_{i,k}^{2}} [(\frac{∂}{∂x} (\overline{ρ}\overline{θ}u))_{i,k} + (\frac{∂}{∂z} (\overline{ρ}\overline{θ}w))_{i,k}]

+\frac{\overline{c}_{i,k}^{2}}{c_{pv}\overline{ρ}_{i,k}\overline{θ}_{i,k}} (\frac{L}{c_{pv}\overline{T}_{i,k}} − 1) [M_{cond}]_{i,k}

+\frac{\overline{c}_{i,k}^{2}}{c_{pv}\overline{θ}_{i,k}^{2}\overline{Π}_{i,k}} ([Q_{rad}]_{i,k} + [Q_{dis}]_{i,k}) .  (1.34)

雲密度の式

\frac{∂ρ_{s}_{i,k}}{∂t} = − [\frac{∂}{∂x} (ρ_{s}u)]_{i,k} − [\frac{∂}{∂z} (ρ_{s}w)]_{i,k} + [M_{cond}]_{i,k} + [M_{fall}]_{i,k}
+ [D_{ρ_{s}}]_{i,k} .  (1.35)

ここで

[M_{cond}]_{i,k} = \frac{4πr_{d}_{i,k}\overline{ρ}_{i,k}N^{∗}}{R_{h}_{i,k}} (S_{i,k} − 1),  (1.36)
R_{h}_{i,k} = \frac{L^{2}}{kR[θ^{2}]_{i,k}[Π^{2}]_{i,k}} ,  (1.37)
S_{i,k} = \frac{p_{0}^{\frac{cpv}{R}} Π_{i,k}}{[p_{∗}]_{i,k}} ,  (1.38)
[M_{fall}]_{i,k} = [\frac{∂}{∂z} (ρ_{s}V_{term})]_{i,k} ,  (1.39)
V_{term}_{i,k} = [C_{sc}]_{i,k} \frac{2 [r_{d}^{2}]_{i,k} gρ_{I}}{9η_{i,k}} ,  (1.40)
r_{d}_{i,k} = (r_{aero}^{3} + \frac{3ρ_{s}_{i,k}}{4ρ_{I}π\overline{ρ}_{i,k}N^{∗}} )^{1/3} ,  (1.41)
[C_{sc}]_{i,k} = 1 + \frac{4}{3} [K_{n}]_{i,k} ,  (1.42)
[K_{n}]_{i,k} = \frac{λ_{i,k}}{r_{d}_{i,k}} ,  (1.43)
λ_{i,k} = \frac{k_{B}\overline{θ}_{i,k}}{\sqrt{2}πσ^{2}p_{0}\overline{Π}_{i,k}^{c_{vv}/R}} ,  (1.44)
η_{i,k} = η_{ref} (\frac{T_{ref} + C_{CO_{2}}}{\overline{θ}_{i,k}\overline{Π}_{i,k} + C_{CO_{2}}} ) (\frac{\overline{θ}_{i,k}\overline{Π}_{i,k}}{T_{ref}} )^{3/2}  (1.45)

である.

\overline{2011/12/14(山下達也)}

熱力学の式

\frac{∂θ_{i,k}}{∂t} = −u_{i,k} (\frac{∂θ}{∂x} )_{i,k} − w_{i,k} (\frac{∂θ}{∂z} )_{i,k}
+\frac{1}{\overline{Π}_{i,k}} \frac{L[M_{cond}]_{i,k}}{\overline{ρ}_{i,k}c_{pv}} + \frac{1}{\overline{Π}_{i,k}} ([Q_{rad}]_{i,k} + [Q_{dis}]_{i,k}) + [D_{θ}]_{i,k} . (1.46)

乱流拡散係数の式

\frac{∂K_{m}_{i,k}}{∂t} = −u_{i,k} (\frac{∂K_{m}}{∂x} )_{i,k} − w_{i,k} (\frac{∂K_{m}}{∂z} )_{i,k} − \frac{3gC_{m}^{2}l^{2}}{2\overline{θ}_{i,k}} (\frac{∂θ′}{∂z} )_{i,k}

+C_{m}^{2}l^{2} [(\frac{∂u}{∂x})_{i,k}^{2} + (\frac{∂w}{∂z})_{i,k}^{2}] + \frac{C_{m}^{2}l^{2}}{2} (\frac{∂u}{∂z} + \frac{∂w}{∂x})_{i,k}^{2}

−\frac{K_{m}_{i,k}}{3} (\frac{∂u}{∂x} + \frac{∂w}{∂z})_{i,k} + \frac{1}{2} [(\frac{∂^{2}K_{m}^{2}}{∂x^{2}} )_{i,k} + (\frac{∂^{2}K_{m}^{2}}{∂z^{2}} )_{i,k}]

+ (\frac{∂K_{m}}{∂x} )_{i,k}^{2} + (\frac{∂K_{m}}{∂z} )_{i,k}^{2} − \frac{C_{ε}}{2C_{m}l^{2}} K_{m}_{i,k}^{2}.  (1.47)

1.4 時間方向の離散化

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

1.4.1 モード別時間分割法

一般に安定に計算を進める為には少なくとも CFL 条件を満たしている必要がある. 例えば簡単な例として 1 次元移流方程式

\frac{∂ϕ}{∂t} + U \frac{∂ϕ}{∂x} = 0  (1.48)

を考えると, CFL 条件は

\frac{UΔt}{Δx} ≤ 1  (1.49)

と表される. 但し U, Δt はそれぞれ移流の速さ, 時間ステップである. 系において様々な速度スケールの現象が生じている場合, 最も速い速度スケールを持つ現象が CFL 条件に制約を加えることになる.

準圧縮方程式は音波を解に含んでいる. 本研究では対流に着目しているので, 音波自体はあまり重要ではない. しかし音波の位相速度は対流の速度スケールに比べて 10 倍程度大きい. 従ってたとえ対流のみに着目しようとしても, 計算を安定に進めるために時間ステップを小さくとらなければならなくなり, 計算のコストが高くなってしまう. そこで計算の効率化を図るためにモード別時間分割法を採用する. モード別時間分割法とは時間ステップを 2 種類用意し, 短い方の時間ステップで音波に関連する項を解き, 長い時間ステップで音波に関連しない移流項や拡散項を解くという方法である. 短い時間ステップで時間積分を行なっている間は長い時間ステップで評価する項の値は一定とみなして計算を行なう. モード別時間分割法の概念図を図 1.2 に示す.

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

1.4.2 音波減衰項

モード別時間分割法を用いると音波についての CFL 条件を満たしているにもかかわらず計算不安定を起こすことがある (Skamarock and Klemp, 1992). この計算不安定を抑制する為に, 運動方程式 (1.32), (1.33) の Π′ を

Π′ − αDiv  (1.50)

で置き換える. 但し Div ≡ (∂u/∂x) + (∂w/∂z) である. 音波を選択的に減衰させる為には, 係数 α を適切な値に設定する必要がある. 本モデルでは北守 (2006) に従い

α = 5.0 × 10^{−7} × \frac{min[(Δx)^{2}, (Δz)^{2}]}{Δτ}  (1.51)

図 1.2: モード別時間分割法の概念図. 北守 (2006) より引用した.

とする.

1.4.3 数値粘性項

移流項の空間微分を中心差分を用いて離散化すると計算不安定が生じることがある. この計算不安定を抑制する為に運動方程式, 熱力学の式, 雲密度の式の移流項に人工的な数値粘性項を加える. 即ち任意の予報変数 ϕ に関する方程式に数値粘性項

ν_{H} \frac{∂^{2}ϕ}{∂x^{2}} + ν_{V} \frac{∂^{2}ϕ}{∂z^{2}}  (1.52)

を付加する. 係数 ν_{H}, ν_{V} については

ν_{H} = ν_{V} = 1.0 × 10^{−4} × \frac{min[(Δx)^{2}, (Δz)^{2}]}{Δt}  (1.53)

とする.

1.4.4 準圧縮方程式系の時間方向の離散化

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

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

水平方向の運動方程式の離散化

(1.32) を離散化すると

u_{i(u),k}^{τ+Δτ} = u_{i(u),k}^{τ} − (c_{pv}\overline{θ}\frac{∂[Π′]^{τ}}{∂x} )_{i,u(k)} Δτ + F_{u}_{i(u),k}^{t}Δτ  (1.54)

となる. F_{u} は音波, 凝結に関連しない項で

F_{u}_{i(u),k}^{t} = − (u^{t} \frac{∂u^{t}}{∂x} )_{i(u),k} − (w^{t} \frac{∂u^{t}}{∂z} )_{i(u),k} + [D_{u}]_{i(u),k}^{t−Δt}  (1.55)

である. 乱流拡散項は坪木・榊原 (2001) と同様に全て時刻 t−Δt での値で評価する.

鉛直方向の運動方程式と圧力方程式の離散化

本モデルでは HE-VI 法を用いるので, 鉛直方向の運動方程式と圧力方程式を連立させて解く. (1.33) において音波減衰項, 圧力勾配項はそれぞれ前進差分, 後退差分で離散化する. (1.34) において水平方向のフラックス項, 鉛直方向のフラックス項はそれぞれ前進差分, 後退差分で離散化する. (1.33) 及び (1.34) を離散化すると

w_{i,k(w)}^{τ+Δτ} = w_{i,k(w)}^{τ} − c_{pv}\overline{θ}_{i,k(w)} [β \frac{∂[Π′]^{τ+Δτ}}{∂z} + (1 − β)\frac{∂[Π′]^{τ}}{∂z} ]_{i,k(w)} Δτ
+F_{w}_{i,k(w)}^{t}Δτ,  (1.56)
[Π′]_{i,k}^{τ+Δτ} = [Π′]_{i,k}^{τ} − [\frac{\overline{c}^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} \frac{∂}{∂x} (\overline{ρ}\overline{θ}u^{τ+Δτ})]_{i,k} Δτ

−β [\frac{\overline{c}^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} \frac{∂}{∂z} (\overline{ρ}\overline{θ}w^{τ+Δτ})]_{i,k} Δτ − (1 − β) [\frac{\overline{c}^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} \frac{∂}{∂z} (\overline{ρ}\overline{θ}w^{τ})]_{i,k} Δτ

+ [\frac{\overline{c}^{2}}{c_{pv}\overline{θ}} (\frac{L}{c_{pv}\overline{T}} − 1) \frac{M_{cond}^{τ}}{\overline{ρ}} ]_{i,k} Δτ

+ [\frac{\overline{c}^{2}}{c_{pv}\overline{θ}^{2}\overline{Π}}(Q_{rad} + Q_{dis})]_{i,k} Δτ  (1.57)

となる. 但し β は後退差分における重み係数を表し, クランク・ニコルソン法の場合 β = 0.5 とする. また F_{w} は音波, 凝結に関連しない項で

F_{w}_{i,k(w)}^{t} = − (u^{t} \frac{∂w^{t}}{∂x} )_{i,k(w)} − (w^{t} \frac{∂w^{t}}{∂z} )_{i,k(w)} + g \frac{[θ′]_{i,k(w)}^{t}}{\overline{θ}_{i,k(w)}} − g \frac{R}{p_{0}} \begin{matrix}\overline{θ} & i,k(w) \\ \overline{Π} & {c_{vv}/R \atop i,k(w)} \end{matrix} ρ_{s}_{i,k(w)}^{t}
+ [D_{w}]_{i,k(w)}^{t−Δt}  (1.58)

である. (1.57) に (1.54), (1.56) を代入して u^{τ+Δτ}, w^{τ+Δτ} を消去し, 式を整理すると以下のような行列で表記出来る.

(\begin{matrix}A_{1} & B_{2} & ⋯ & 0 & 0 \\ C_{1} & A_{2} & ⋯ & 0 & 0 \\ : & : & & : & : \\ 0 & 0 & ⋯ & A_{km−1} & B_{km} \\ 0 & 0 & ⋯ & C_{km−1} & A_{km} \end{matrix}) (\begin{matrix}Π′_{1,1} & Π′_{2,1} & ⋯ & Π′_{km−1,1} & Π′_{km,1} \\ Π′_{1,2} & Π′_{2,2} & ⋯ & Π′_{km−1,2} & Π′_{km,2} \\ : & : & & : & : \\ Π′_{1,km−1} & Π′_{2,km−1} & ⋯ & Π′_{km−1,km−1} & Π′_{km,km−1} \\ Π′_{1,km} & Π′_{2,km} & ⋯ & Π′_{km−1,km} & Π′_{km,km} \end{matrix})^{τ+Δτ}

= (\begin{matrix}D_{1,1} & D_{2,1} & ⋯ & D_{km−1,1} & D_{km,1} \\ D_{1,2} & D_{2,2} & ⋯ & D_{km−1,2} & D_{km,2} \\ : & : & & : & : \\ D_{1,km−1} & D_{2,km−1} & ⋯ & D_{km−1,km−1} & D_{km,km−1} \\ D_{1,km} & D_{2,km} & ⋯ & D_{km−1,km} & D_{km,km} \end{matrix})^{τ}  (1.59)

\underline{大気主成分凝結対流離散化 13} ここで km は鉛直方向の格子点の総数である. (1.59) の左辺の係数行列の各成分は以下のように表される.

A_{1} = 1 + β^{2} [\frac{\overline{c}^{2}(Δτ)^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{1} \frac{1}{(Δz)^{2}} (c_{pv}\overline{ρ}\overline{θ}^{2})_{i,1(w)},  (1.60)
A_{k} = 1 + β^{2} [\frac{\overline{c}^{2}(Δτ)^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{k} \frac{1}{(Δz)^{2}} [(c_{pv}\overline{ρ}\overline{θ}^{2})_{i,k(w)} + (c_{pv}\overline{ρ}\overline{θ}^{2})_{i,k−1(w)}],
(k = 2, 3, …, km − 1)  (1.61)
A_{km} = 1 + β^{2} [\frac{\overline{c}^{2}(Δτ)^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{km} \frac{1}{(Δz)^{2}} (c_{pv}\overline{ρ}\overline{θ}^{2})_{i,km−1(w)},  (1.62)
B_{k} = −β^{2} [\frac{\overline{c}^{2}(Δτ)^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{k−1} \frac{1}{(Δz)^{2}} (c_{pv}\overline{ρ}\overline{θ}^{2})_{i,k−1(w)},  (1.63)
C_{k} = −β^{2} [\frac{\overline{c}^{2}(Δτ)^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{k+1} \frac{1}{(Δz)^{2}} (c_{pv}\overline{ρ}\overline{θ}^{2})_{i,k(w)},  (1.64)
D_{i,1} = [Π′]_{i,1}^{τ} − (1 − β) [\frac{\overline{c}^{2}Δτ}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{1} [\frac{∂(\overline{ρ}\overline{θ}w^{τ})}{∂z}]_{i,1} − [\frac{\overline{c}^{2}Δτ}{c_{pv}\overline{θ}} ]_{1} [\frac{∂u^{τ+Δτ}}{∂x} ]_{i,1} + F_{i,1}
−β [\frac{\overline{c}^{2}(Δτ)^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{1} \frac{1}{Δz} (c_{pv}\overline{ρ}\overline{θ}^{2})_{i,0(w)}E_{i,0(w)},  (1.65)
D_{i,k} = [Π′]_{i,k}^{τ} − (1 − β) [\frac{\overline{c}^{2}Δτ}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{k} [\frac{∂(\overline{ρ}\overline{θ}w^{τ})}{∂z}]_{i,k} − [\frac{\overline{c}^{2}Δτ}{c_{pv}\overline{θ}} ]_{k} [\frac{∂u^{τ+Δτ}}{∂x} ]_{i,k} + F_{i,k},
(k = 2, 3, …, km − 1)  (1.66)
D_{i,km} = [Π′]_{i,km}^{τ} − (1 − β) [\frac{\overline{c}^{2}Δτ}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{km} [\frac{∂(\overline{ρ}\overline{θ}w^{τ})}{∂z}]_{i,km} − [\frac{\overline{c}^{2}Δτ}{c_{pv}\overline{θ}} ]_{km} [\frac{∂u^{τ+Δτ}}{∂x} ]_{i,km}
+F_{i,km} + β [\frac{\overline{c}^{2}(Δτ)^{2}}{c_{pv}\overline{ρ}\overline{θ}^{2}} ]_{km} \frac{1}{Δz} (c_{pv}\overline{ρ}\overline{θ}^{2})_{i,km(w)}E_{i,km(w)}.  (1.67)

但し

E_{i,k(w)} = [\frac{∂ (αDiv)^{τ}}{∂z} ]_{i,k(w)} − (1 − β) [\frac{∂[Π′]^{τ}}{∂z} ]_{i,k(w)} + [\frac{F_{w}^{t}}{c_{pv}\overline{θ}}]_{i,k(w)} ,  (1.68)
F_{i,k} = −β (\frac{\overline{c}^{2}Δτ}{c_{pv}\overline{ρ}\overline{θ}})_{i,k}

[\frac{∂}{∂z} \{\overline{ρ}\overline{θ}(w^{τ} − c_{pv}\overline{θ}Δτ \{(1 − β)\frac{∂[Π′]^{τ}}{∂z} − \frac{∂ (αDiv)^{τ}}{∂z} \} + F_{w}^{t}Δτ)\}_{i,k(w)}

+\frac{\overline{c}^{2}}{c_{pv}\overline{θ}} (\frac{L}{c_{pv}\overline{T}} − 1) \frac{M_{cond}^{τ}}{\overline{ρ}} Δτ + \frac{\overline{c}^{2}}{c_{pv}\overline{θ}^{2}\overline{Π}}(Q_{rad} + Q_{dis})Δτ]_{i,k}  (1.69)

\overline{2011/12/14(山下達也)}

\underline{14 大気主成分凝結対流離散化} である. (1.59) の左辺の係数行列は 3 重対角行列となっているので, Thomas 法を用いて時刻 τ + Δτ での Π′ の値を求めることができる (Thomas, 1949). 本モデルでは計算ライブラリ LAPACK を用いて (1.59) を解いている. 雲密度の式の離散化雲密度の時間積分を行なうに当たっては, 先ず音波・凝結に関連しない項と雲粒落下項を計算し, 次に負の雲密度の発生を防ぐ処置を行い, 最後に音波・凝結に関連する項を計算する^{1 )}. 即ち, (1.35) を離散化すると以下のように書ける.

ρ_{s}_{i,k}^{∗} = ρ_{s}_{i,k}^{τ} + [M_{fall}]_{i,k}^{τ} Δτ + [F_{ρ_{s}}]_{i,k}^{t} Δτ,  (1.70)
ρ_{s}_{i,k}^{∗∗} = FillNeg.[ρ_{s}_{i,k}^{∗}],  (1.71)
ρ_{s}_{i,k}^{τ+Δτ} = ρ_{s}_{i,k}^{∗∗} + [M_{cond}]_{i,k}^{∗∗} Δτ.  (1.72)

ここで FillNeg.[ρ_{s}_{i,k}^{∗}] は負の雲密度の発生を防ぐ為の処置を表し, ρ_{s}_{i,k}^{∗∗} は処理を行なった後の雲密度の暫定値, F_{ρ_{s}} は音波・凝結に関連しない項であり,

[F_{ρ_{s}}]_{i,k}^{t} = − [\frac{∂}{∂x} (ρ_{s}^{t}u^{t})]_{i,k} − [\frac{∂}{∂z} (ρ_{s}^{t}w^{t})]_{i,k} + D_{ρ_{s}} _{i,k}^{t}.  (1.73)

である. 熱力学の式の離散化 (1.46) を離散化すると

θ′_{i,k}^{τ+Δτ} = θ′_{i,k}^{τ} + F_{θ}_{i,k}^{t}Δτ + \frac{1}{\overline{Π}_{i,k}} \frac{L [M_{cond}]_{i,k}^{∗∗}}{c_{pv}\overline{ρ}_{i,k}} Δτ  (1.74)

となる. F_{θ} は音波・凝結に関連しない項であり

F_{θ}_{i,k}^{t} = −u_{i,k}^{t} (\frac{∂θ′^{t}}{∂x} )_{i,k} − w_{i,k}^{t} (\frac{∂θ′^{t}}{∂z} )_{i,k} − w_{i,k}^{t} (\frac{∂\overline{θ}^{t}}{∂z} )_{i,k}
+\frac{1}{\overline{Π}_{i,k}} (Q_{rad} + Q_{dis})_{i,k}^{t} + D_{θ}_{i,k}^{t−Δt}  (1.75)

である. 乱流拡散係数の式の離散化 (1.47) をリープフロッグ法を用いて離散化すると

K_{m}_{i,k}^{t+Δt} = K_{m}_{i,k}^{t−Δt} + 2Δt[F_{K_{m}}]_{i,k}^{t}  (1.76)

1 )負の雲密度の発生を防ぐ処置の具体的な手順については, 付録 A を参照されたい.

\overline{2011/12/14(山下達也)}

となる. ここで

[F_{K_{m}}]_{i,k}^{t} = − [u_{i(u),k} (\frac{∂K_{m}}{∂x} )_{i(u),k}]_{i,k}^{t} − [w_{i,k(w)} (\frac{∂K_{m}}{∂z} )_{i,k(w)}]_{i,k}^{t}

− [\frac{3gC_{m}^{2}l^{2}}{2\overline{θ}} (\frac{∂θ′}{∂z} )_{i,k(w)}]_{i,k}^{t−Δt}

+C_{m}^{2}l^{2} [(\frac{∂u}{∂x})^{2}]_{i,k}^{t−Δt} + C_{m}^{2}l^{2} [(\frac{∂w}{∂z})^{2}]_{i,k}^{t−Δt}

+\frac{C_{m}^{2}l^{2}}{2} [\{(\frac{∂u}{∂z})_{i(u),k(w)}\}_{i,k}^{t−Δt} + \{(\frac{∂w}{∂x})_{i(u),k(w)}\}_{i,k}^{t−Δt}]^{2}

−\frac{1}{3}[K_{m}]_{i,k}^{t−Δt} [(\frac{∂u}{∂x})_{i,k}^{t−Δt} + (\frac{∂w}{∂z})_{i,k}^{t−Δt}]

+\frac{1}{2} [\{\frac{∂}{∂x} (\frac{∂K_{m}^{2}}{∂x} )_{i(u),k}\}_{i,k}^{t−Δt} + \{\frac{∂}{∂z} (\frac{∂K_{m}^{2}}{∂z} )_{i,k(w)}\}_{i,k}^{t−Δt}]

+ [\{(\frac{∂K_{m}}{∂x} )^{2}\}_{i(u),k}]_{i,k}^{t−Δt} + [\{(\frac{∂K_{m}}{∂z} )^{2}\}_{i,k(w)}]_{i,k}^{t−Δt}

−\frac{C_{ε}}{2C_{m}l^{2}} [K_{m}^{2}]_{i,k}^{t−Δt}  (1.77)

である. 坪木・榊原 (2001) と同様に移流項を時刻 t で, それ以外の項を時刻 t − Δt で評価した.

1.4.5 Robert, Asselin の時間フィルター

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

時間フィルター適用前の変数を ϕ, 時間フィルター適用後の変数を\overline{ϕ} とすると

\overline{ϕ^{t}} = ϕ^{t} + μ_{a}(\overline{ϕ^{t−Δt}} − 2ϕ^{t} + ϕ^{t+Δt})  (1.78)

\underline{16 大気主成分凝結対流離散化} と表される. ここで μ_{a} はフィルターの強さを表す係数であり, 本モデルでは Klemp and Wilhelmson(1978), 坪木・榊原 (2001) と同様に μ_{a} = 0.1 とした.

\underline{大気主成分凝結対流参考文献 17}

参考文献

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

Asselin, R., 1972 : Frequency filter for time integrations. Monthly Weather Review, 100, 487 – 490

Deardorff, J. W., 1972: Numerical investigation of neutral and unstable planetary boundary layers. Journal of the Atmospheric Sciences, 29, 91 – 115

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

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

Lorenz, E. N., 1960 : Energy and numerical weather prediction. Tellus, 12, 364 – 373

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

Robert, A., 1966 : The integration of a low order spectral form of the primitive meteorological equations. Journal of the Meteorological Society of Japan, 44, 237 – 244

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

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

北守太一, 2006: 大気主成分の凝結を考慮した火星大気非静力学モデルの開発. 北 \underline{18 大気主成分凝結対流参考文献}

海道大学大学院地球惑星科学専攻修士論文, 99 pp http://www.gfd-dennou.org/arch/prepri/2006/hokudai/kitamo/paper/pub/main.pdf

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

坪木和久, 榊原篤志, 2001: CReSS ユーザーズガイド第 2 版. 253 pp http://cf.tokyo.rist.or.jp/CReSS.top.html \underline{大気主成分凝結対流負の雲密度の発生を防ぐ処置に関して 19}

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

本もデルでは中央差分を用いることで生じる負の雲密度を防ぐ為, 全質量が保存されるように負の雲密度の穴埋めを行なっている. 雲密度に対する穴埋めを行なう関数 FillNegativeMMC の中の式を具体的に書き下すと, 以下のようになる.

QSUMP N(i, k) = 1.0/

{0.75 ∗ [max(0.0, ρ_{s}(i − 1, k)) +max(0.0, ρ_{s}(i + 1, k)) +max(0.0, ρ_{s}(i, k − 1)) +max(0.0, ρ_{s}(i, k + 1))] +0.25 ∗ [max(0.0, ρ_{s}(i − 2, k)) +max(0.0, ρ_{s}(i + 2, k)) +max(0.0, ρ_{s}(i, k − 2))

+max(0.0, ρ_{s}(i, k + 2))]}  (A.1)
DQF ILL(i, k) = −min(0.0, ρ_{s}(i, k))

+max(0.0, ρ_{s}(i, k)) ∗{0.75 ∗ [min(0.0, ρ_{s}(i − 1, k)) ∗ QSUMP N(i − 1, k) +min(0.0, ρ_{s}(i + 1, k)) ∗ QSUMP N(i + 1, k) +min(0.0, ρ_{s}(i, k − 1)) ∗ QSUMP N(i, k − 1) +min(0.0, ρ_{s}(i, k + 1)) ∗ QSUMP N(i, k + 1)] +0.25 ∗ [min(0.0, ρ_{s}(i − 2, k)) ∗ QSUMP N(i − 2, k) +min(0.0, ρ_{s}(i + 2, k)) ∗ QSUMP N(i + 2, k) +min(0.0, ρ_{s}(i, k − 2)) ∗ QSUMP N(i, k − 2) +min(0.0, ρ_{s}(i, k + 2)) ∗ QSUMP N(i, k + 2)]}

(A.2)

ρ̂_{s}(i, k) = ρ_{s}(i, k) + DQF ILL(i, k).  (A.3)

\underline{20 大気主成分凝結対流負の雲密度の発生を防ぐ処置に関して}

但し, ρ_{s} は穴埋め前の雲密度, ρ̂_{s} は穴埋め後の雲密度, DQF ILL(i, k) は穴埋めする量である. 例えばある一点のみで雲密度が負となった場合, この関数は負の点がゼロとなるように隣接する点と 2 つ隣の点から雲を削って埋める. その際, 隣接する各点は埋める量全体の 3/16 を負担し, 2 つ隣の各点は埋める量全体の 1/16 を負担する.

以下では, 単純な場合として, ある一点 (i, k) = (M, N) のみ穴埋め前の雲密度が負である場合に, 穴埋め前後での全雲質量が保存されることを示す. ある一点 (i, k) = (M, N) のみ穴埋め前の雲密度が負であるとき, (A.1), (A.2) より (i, k) = (M − 2, N), (M − 1, N), (M, N), (M + 1, N), (M + 2, N), (M, N − 2), (M, N − 1), (M, N + 1), (M, N + 2) 以外の点では DQF ILL(i, k) = 0 となる. 以上の 9 点での DQF ILL(i, k) を書き下すと,

DQF ILL(M − 2, N) = 0.25 ∗ ρ_{s}(M − 2, N)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.4)

DQF ILL(M − 1, N) = 0.75 ∗ ρ_{s}(M − 1, N)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.5)

DQF ILL(M, N) = −ρ_{s}(M, N)

(A.6)

DQF ILL(M + 1, N) = 0.75 ∗ ρ_{s}(M + 1, N)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.7)

DQF ILL(M + 2, N) = 0.25 ∗ ρ_{s}(M + 2, N)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.8)

DQF ILL(M, N − 2) = 0.25 ∗ ρ_{s}(M, N − 2)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.9)

DQF ILL(M, N − 1) = 0.75 ∗ ρ_{s}(M, N − 1)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.10)

DQF ILL(M, N + 1) = 0.75 ∗ ρ_{s}(M, N + 1)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.11)

DQF ILL(M, N + 2) = 0.25 ∗ ρ_{s}(M, N + 2)ρ_{s}(M, N) ∗ QSUMP N(M, N),

(A.12)

となる. 但し

QSUMP N(M, N) = {0.75 ∗ [ρ_{s}(M − 1, N) + ρ_{s}(M + 1, N)
+ρ_{s}(M, N − 1) + ρ_{s}(M, N + 1)]
+0.25 ∗ [ρ_{s}(M − 2, N) + ρ_{s}(M + 2, N)
+ρ_{s}(M, N − 2) + ρ_{s}(M, N + 2)]}^{−1}  (A.13)

\underline{大気主成分凝結対流負の雲密度の発生を防ぐ処置に関して 21} である. 従って

∑_{i} ∑_{k} ρ̂_{s}(i, k) = ∑_{i} ∑_{k} [ρ_{s}(i, k) + DQF ILL(i, k)]
= ∑_{i} ∑_{k} ρ_{s}(i, k)

+DQF ILL(M − 2, N) + DQF ILL(M − 1, N) +DQF ILL(M, N) + DQF ILL(M + 1, N) +DQF ILL(M + 2, N) + DQF ILL(M, N − 2) +DQF ILL(M, N − 1) + DQF ILL(M, N + 1) +DQF ILL(M, N + 2)

= ∑_{i} ∑_{k} ρ_{s}(i, k) − ρ_{s}(M, N)

+ρ_{s}(M, N) ∗ QSUMP N(M, N) {0.75 ∗ [ρ_{s}(M − 1, N) + ρ_{s}(M + 1, N) +ρ_{s}(M, N + 1) + ρ_{s}(M, N + 1)] +0.25 ∗ [ρ_{s}(M − 2, N) + ρ_{s}(M + 2, N) +ρ_{s}(M, N + 2) + ρ_{s}(M, N + 2)]}

= ∑_{i} ∑_{k} ρ_{s}(i, k) − ρ_{s}(M, N)

+ρ_{s}(M, N) ∗ QSUMP N(M, N) ∗ QSUMP N(M, N)^{−1}

= ∑_{i} ∑_{k} ρ_{s}(i, k)  (A.14)

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