A 離散化

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

1.1 離散化の概要

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

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

1.2 格子の配置

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

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

7#7

1.3 空間方向の離散化

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

1.3.1 平均操作

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

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

11#11 12#12 13#13 (1.1)

14#14 12#12 15#15 (1.2)

16#16 12#12 17#17 (1.3)

18#18 12#12 19#19 (1.4)

20#20 12#12 21#21 (1.5)

22#22 12#12 23#23 (1.6)

24#24 12#12 25#25 (1.7)

26#26 12#12 27#27 (1.8)

28#28 12#12 29#29 (1.9)

1.3.2 空間微分の離散化

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

以下に 2 次精度の中心差分を用いた微分操作を示す. 但し 30#30 は格子の角に相当する点で評価している変数を表す.

31#31 12#12 32#32 (1.10)

33#33 12#12 34#34 (1.11)

35#35 12#12 36#36 (1.12)

37#37 12#12 38#38 (1.13)

39#39 12#12 40#40 (1.14)

41#41 12#12 42#42 (1.15)

43#43 12#12 44#44 (1.16)

45#45 12#12 46#46 (1.17)

47#47 12#12 48#48 (1.18)

49#49 12#12 50#50 (1.19)

51#51 12#12 52#52 (1.20)

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

31#31 12#12 53#53 (1.21)

33#33 12#12 54#54 (1.22)

35#35 12#12 55#55 (1.23)

37#37 12#12 56#56 (1.24)

39#39 12#12 57#57 (1.25)

41#41 12#12 58#58 (1.26)

43#43 12#12 59#59 (1.27)

45#45 12#12 60#60 (1.28)

47#47 12#12 61#61 (1.29)

49#49 12#12 62#62

63#63 (1.30)

51#51 12#12 64#64

65#65 (1.31)

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

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

水平方向の運動方程式

66#66 12#12 67#67

68#68 (1.32)

ここで

12#12 70#70 71#71

72#72 (1.33)

73#73 (1.34)

である.

鉛直方向の運動方程式

74#74 12#12 75#75

76#76 (1.35)

ここで

12#12 78#78 79#79

80#80 (1.36)

である.

圧力方程式

81#81 12#12 82#82 83#83

84#84 (1.37)

雲密度の式

85#85 12#12 86#86

87#87 (1.38)

ここで

88#88 12#12 89#89 (1.39)

90#90 12#12 91#91 (1.40)

92#92 12#12 93#93 (1.41)

94#94 12#12 95#95 (1.42)

96#96 12#12 97#97 (1.43)

98#98 12#12 99#99 (1.44)

100#100 12#12 101#101 (1.45)

102#102 12#12 103#103 (1.46)

104#104 12#12 105#105 (1.47)

106#106 12#12 107#107 (1.48)

108#108 12#12 109#109 (1.49)

である.

熱力学の式

110#110 12#12 111#111

112#112 (1.50)

ここで

113#113 12#12 114#114 (1.51)

115#115 12#12 116#116 (1.52)

である.

乱流拡散係数の式

117#117 12#12 118#118 119#119 120#120 121#121

122#122 (1.53)

1.4 時間方向の離散化

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

1.4.1 モード別時間分割法

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

123#123 (1.54)

を考えると, CFL 条件は

124#124 (1.55)

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

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

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

126#126

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

1.4.2 音波減衰項

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

128#128 (1.56)

で置き換える. 但し 129#129 である. 音波を選択的に減衰させる為には, 係数 130#130 を適切な値に設定する必要が ある. 本モデルでは北守(2006)に従い

131#131 (1.57)

とする ^{1 )}.

1.4.3 数値粘性項

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

132#132 (1.58)

を付加する. 係数 133#133, 134#134 については

135#135 (1.59)

とする.

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

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

以下, 時間積分により求まる量を 138#138, 139#139, 時間積分によって得られている最新の物理量を 140#140, 141#141, 最新の物理量の 1 ステップ前の時刻での物理量を 142#142, 143#143 と表すこととする.

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

式 RisanAA を離散化すると

144#144 (1.60)

となる. 145#145 は音波, 凝結に関連しない項で

146#146 (1.61)

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

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

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

式 RisanAB 及び 式 RisanAD を離散化すると

148#148 12#12 149#149 150#150

151#151 (1.62)

152#152 12#12 153#153 154#154 155#155

156#156 (1.63)

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

160#160 12#12 161#161

162#162 (1.64)

である. 式 TrisanF に 式 TrisanC, TrisanE を代入して 163#163, 164#164 を消去し, 式を整理す ると以下のような行列で表記出来る ^{2 )}.

166#166 (1.65)

ここで 167#167 は鉛直方向の格子点の総数である. 式 TrisanH の左辺の係数行列の各成分は以下のように表される.

168#168 12#12 169#169 (1.66)

170#170 12#12 171#171

172#172 (1.67)

173#173 12#12 174#174 (1.68)

175#175 12#12 176#176 (1.69)

177#177 12#12 178#178 (1.70)

179#179 12#12 180#180

181#181 (1.71)

182#182 12#12 183#183

184#184 (1.72)

185#185 12#12 186#186

187#187 (1.73)

但し

188#188 12#12 189#189 (1.74)

190#190 12#12 191#191 192#192 193#193

194#194 (1.75)

である. 式 TrisanH の左辺の係数行列は 3 重対角行列となっているので, Thomas 法を用いて時刻 195#195 での 127#127 の値を求 めることができる(Thomas, 1949). 本モデルでは計算ライブラリ LAPACK を用いて 式 TrisanH を解いて いる.

雲密度の式の離散化

雲密度の時間積分を行なうに当たっては, 先ず音波・凝結に関連しない項と雲粒落下項を計算し, 次に負の雲密度の発生を防ぐ処置を行い, 最後に音波・凝結に関連する項を計算する ^{3 )}. 即ち, 式 RisanAE を離散化すると以下のように書ける.

196#196 12#12 197#197 (1.76)

198#198 12#12 199#199 (1.77)

200#200 12#12 201#201 (1.78)

ここで 202#202 は負の雲密度の発生を防ぐ為 の処置を表し, 203#203 は処理を行なった後の雲密度の暫定値, 204#204 は音波・凝結に関連しない項であり,

205#205 12#12 206#206 (1.79)

である.

熱力学の式の離散化

式 RisanACを離散化すると

207#207 (1.80)

となる. 208#208 は音波・凝結に関連しない項であり

209#209 12#12 210#210

211#211 (1.81)

である.

乱流拡散係数の式の離散化

式 RisanAFをリープフロッグ法を用いて離散化すると

212#212 (1.82)

となる. ここで

213#213 12#12 214#214 215#215 216#216 217#217 218#218 219#219 220#220

221#221 (1.83)

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

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

1.4 節で述べた通り, 長い時間ステップの計算ではリープフロッグ法を用いてい る. リープフロッグ法では時刻 147#147 の値を用いて時刻 223#223 の 値を求める為, 隣接する時間ステップ間で物理量の値に大きな食い違いや振動が 生じる恐れがある. この問題を回避する為, 長い時間ステップの計算を 1 回行 なう度に Robert(1966), Asselin(1972) が考案した時間フィルターを適用する. 本モデルでは 136#136, 137#137, 224#224, 127#127, 225#225, 226#226 の時間積分に対してフィルターをかける.

時間フィルター適用前の変数を 8#8, 時間フィルター適用後の変数を 227#227 とすると

228#228 (1.84)

と表される. ここで 229#229 はフィルターの強さを表す係数であり, 本モデルでは Klemp and Wilhelmson(1978), 坪木・榊原(2001) と同様に 230#230 とした.

Footnotes

... とする^{1 )}

基本場や計算領域の鉛直方向の長さを変更すると, 計算が即破綻してしまうこと がある. その際, 経験的には 130#130 の値を適宜変えてみると計算がうまく進む場合が 多い.

... ると以下のような行列で表記出来る^{2 )}

式 TrisanH - 式 TrisanR の導出については付 録 A を参照されたい.

... 最後に音波・凝結に関連する項を計算する^{3 )}

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