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)

鉛直方向の運動方程式

69#69 12#12 70#70

71#71 (1.33)

圧力方程式

72#72 12#12 73#73 74#74

75#75 (1.34)

雲密度の式

76#76 12#12 77#77

78#78 (1.35)

ここで

79#79 12#12 80#80 (1.36)

81#81 12#12 82#82 (1.37)

83#83 12#12 84#84 (1.38)

85#85 12#12 86#86 (1.39)

87#87 12#12 88#88 (1.40)

89#89 12#12 90#90 (1.41)

91#91 12#12 92#92 (1.42)

93#93 12#12 94#94 (1.43)

95#95 12#12 96#96 (1.44)

97#97 12#12 98#98 (1.45)

である.

熱力学の式

99#99 12#12 100#100

101#101 (1.46)

乱流拡散係数の式

102#102 12#12 103#103 104#104 105#105

106#106 (1.47)

1.4 時間方向の離散化

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

1.4.1 モード別時間分割法

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

107#107 (1.48)

を考えると, CFL 条件は

108#108 (1.49)

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

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

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

110#110

凝結に関連する項は音波にも移流にも直接関連しないので, 解くべき時間ステッ プは凝結の時間スケールによって決まると考えられる. 北守(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 の 111#111 を

112#112 (1.50)

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

115#115 (1.51)

とする.

1.4.3 数値粘性項

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

116#116 (1.52)

を付加する. 係数 117#117, 118#118 については

119#119 (1.53)

とする.

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

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

以下, 時間積分により求まる量を 122#122, 123#123, 時間積分によって得られている最新の物理量を 124#124, 125#125, 最新の物理量の 1 ステップ前の時刻での物理量を 126#126, 127#127 と表すこととする.

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

RisanAA を離散化すると

128#128 (1.54)

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

130#130 (1.55)

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

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

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

132#132 12#12 133#133

134#134 (1.56)

135#135 12#12 136#136 137#137 138#138

139#139 (1.57)

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

143#143 12#12 144#144

145#145 (1.58)

である. TrisanF に TrisanC, TrisanE を代入して 146#146, 147#147 を消去し, 式を整理す ると以下のような行列で表記出来る.

149#149 (1.59)

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

151#151 12#12 152#152 (1.60)

153#153 12#12 154#154

155#155 (1.61)

156#156 12#12 157#157 (1.62)

158#158 12#12 159#159 (1.63)

160#160 12#12 161#161 (1.64)

162#162 12#12 163#163

164#164 (1.65)

165#165 12#12 166#166

167#167 (1.66)

168#168 12#12 169#169

170#170 (1.67)

但し

171#171 12#12 172#172 (1.68)

173#173 12#12 174#174 175#175

176#176 (1.69)

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

雲密度の式の離散化

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

178#178 12#12 179#179 (1.70)

180#180 12#12 181#181 (1.71)

182#182 12#12 183#183 (1.72)

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

187#187 12#12 188#188 (1.73)

である.

熱力学の式の離散化

RisanACを離散化すると

189#189 (1.74)

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

191#191 12#12 192#192

193#193 (1.75)

である.

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

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

194#194 (1.76)

となる. ここで

195#195 12#12 196#196 197#197 198#198 199#199 200#200 201#201 202#202

203#203 (1.77)

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

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

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

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

210#210 (1.78)

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

Footnotes

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

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