\underline{球面浅水方程式モデル : 渦度発散型目次 1}

SPMODEL サンプルプログラム

球面浅水方程式モデル : 渦度発散型

spshallow-zd *.f90

竹広真一

2004 年 10 月 7 日

目次

1 概要 2

2 支配方程式系 3 2.1 支配方程式系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 渦度・発散を用いた支配方程式系 . . . . . . . . . . . . . . . . . . . 3 2.3 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3 水平離散化 6

4 時間積分 7 4.1 陽解法 (オイラースキーム) . . . . . . . . . . . . . . . . . . . . . . . 7 4.2 セミインプリシット & リープフロッグスキーム . . . . . . . . . . . 7

4.2.1 時間フィルター . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.3 リープフロッグ & クランク・ニコルソン法 . . . . . . . . . . . . . 10 4.4 ルンゲ・クッタ & クランク・ニコルソン法 . . . . . . . . . . . . . 11

5 使用モジュールとその他の設定 14

6 Williamson et al. (1992) の標準実験 15 6.1 テストケース 1: 極を越える cos 型の山の移流 . . . . . . . . . . . . 15 6.2 テストケース 2: 定常な非線形帯状地衡流 . . . . . . . . . . . . . . . 16 6.3 テストケース 3: 局在した非線形帯状地衡流 . . . . . . . . . . . . . 17 6.4 テストケース 4: 移動する低気圧を持つ強制された流れ . . . . . . . 19 6.5 テストケース 5: 孤立した山を越える帯状流 . . . . . . . . . . . . . 21 \underline{球面浅水方程式モデル : 渦度発散型目次 2} 6.6 テストケース 6: ロスビー・ハウルヴィッツ波 . . . . . . . . . . . . 22 7 参考文献 23

spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

\underline{球面浅水方程式モデル : 渦度発散型  1 概要  3}

1 概要

SPMODEL サンプルプログラム『shallow-zd *.f90』に用いられている基礎方程式と境界条件, および, このプログラムを用いた数値実験の方法について解説する. 基礎方程式は, 渦度と発散を用いて表現した球面上の浅水方程式である. 計算はスペクトル法を用いて行い, 展開関数は球面調和関数である. 波数切断は三角切断である. スペクトル変換と逆変換および微分演算には SPMODEL ライブラリ (spml) を用いている. 時間積分はオイラー法, リープフロッグ法と台形公式を組み合わせたセミインプリシット法, およびリープフロッグ法とクランク・ニコルソン法を組み合わせた方法, および 4 次のルンゲ・クッタ法とクランク・ニコルソン法を組み合わせた方法について定式化を記述する.

数値実験として Williamson et al. (1992) によって提案された, 球面浅水方程式モデルの標準試験の実験を行う.

プログラム名

shallow-zd euler *.f90 : オイラー法
shallow-zd lfsi *.f90  : セミインプリシット法
shallow-zd lfcn *.f90  : リープフロッグ & クランク・ニコルソン法
shallow-zd rn4cn *.f90 : ルンゲ・クッタ & クランク・ニコルソン法

プログラム取得元

http://www.gfd-dennou.org/arch/spmodel/2d-sphere-w/shallow-zd/SIGEN.htm

SPMODEL サンプルプログラム目次

http://www.gfd-dennou.org/arch/spmodel/html/sample.htm

SPMODEL の使い方

http://www.gfd-dennou.org/arch/spmodel/

\underline{球面浅水方程式モデル : 渦度発散型  2 支配方程式系  4}

2 支配方程式系

ここでは支配方程式系と境界条件を記す.

2.1 支配方程式系

支配方程式系は球面浅水方程式系である.

2.2 渦度・発散を用いた支配方程式系

球面浅水方程式系を渦度と発散を用いて表した式はつぎのとおりである.

\frac{∂ζ}{∂t} = −\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)U}{∂λ} − \frac{1}{a} \frac{∂(f + ζ)V}{∂μ}
−K_{m} [(−1)^{N_{m}}∇^{2N_{m}} − (\frac{2}{a^{2}} )^{2N_{m}}] ζ,  (1)
\frac{∂D}{∂t} = \frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)V}{∂λ} − \frac{1}{a} \frac{∂(f + ζ)U}{∂μ} − ∇^{2} [g(h + h_{s}) + E]
−K_{m} [(−1)^{N_{m}}∇^{2N_{m}} − (\frac{2}{a^{2}} )^{2N_{m}}] D,  (2)
\frac{∂h}{∂t} = −\frac{1}{a(1 − μ^{2})} \frac{∂(hU)}{∂λ} − \frac{1}{a} \frac{∂(hV)}{∂μ} − (−1)^{N_{h}}K_{h}∇^{2N_{h}}h.  (3)

各記号の定義は表 1 に示した. 渦度 ζ と発散 D , 単位質量あたりの運動エネルギー E, 演算子 ∇^{2} は以下のように定義される.

ζ = \frac{1}{a(1 − μ^{2})} \frac{∂V}{∂λ} − \frac{1}{a} \frac{∂U}{∂μ} ,  (4)
D = \frac{1}{a(1 − μ^{2})} \frac{∂U}{∂λ} + \frac{1}{a} \frac{∂V}{∂μ} ,  (5)
E = \frac{U^{2} + V^{2}}{2(1 − μ^{2})},  (6)
∇^{2}A = \frac{1}{a^{2}(1 − μ^{2})} \frac{∂^{2}A}{∂λ^{2}} + \frac{1}{a^{2}} \frac{∂}{∂μ} [(1 − μ^{2})\frac{∂A}{∂μ}] .  (7)

ここで A は任意のスカラーである. 流線関数 ψ と速度ポテンシャル χ を導入すると,

U = −\frac{1 − μ^{2}}{a} \frac{∂ψ}{∂μ} + \frac{1}{a} \frac{∂χ}{∂λ},  (8)
spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

\underline{球面浅水方程式モデル : 渦度発散型  2 支配方程式系  5}
V = \frac{1}{a} \frac{∂ψ}{∂λ} + \frac{1 − μ^{2}}{a} \frac{∂χ}{∂μ},  (9)
ζ = ∇^{2}ψ,  (10)
D = ∇^{2}χ  (11)

と表される.

記号変数/物理定数
λ  経度
ϕ  緯度
μ  sin 緯度 (≡ sin ϕ)
t  時間
u  経度方向流速
v  緯度方向流速
U  u cos ϕ
V  v cos ϕ
h  流体の全厚さ
ζ  渦度
D  発散
E  運動エネルギー
ψ  流線関数
χ  速度ポテンシャル
f  2Ω sin ϕ
Ω  自転角速度
g  重力加速度
a  惑星半径
K_{m}  水平粘性係数
K_{h}  水平拡散係数
N_{m}  超粘性の次数 (水平ラプラシアンの階数)
N_{h}  超拡散の次数 (水平ラプラシアンの階数)

表 1: 変数, 物理定数の定義

2.3 境界条件

境界条件は緯度, 経度方向ともに周期境界条件である.

ζ(λ + 2π, ϕ) = ζ(λ, ϕ),  (12)
spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

\underline{球面浅水方程式モデル : 渦度発散型  2 支配方程式系  6}
ζ(λ, ϕ + 2π) = ζ(λ, ϕ).  (13)

などと表される. 後の時間積分スキームの構成のために, h を平均成分\overline{h} とそれからの擾乱 h′ を用いて h = \overline{h} + h′(λ, μ, t) と表しておくと,

\frac{∂ζ}{∂t} = −\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)U}{∂λ} − \frac{1}{a} \frac{∂(f + ζ)V}{∂μ}
−K_{m} [(−1)^{N_{m}}∇^{2N_{m}} − (\frac{2}{a^{2}} )^{2N_{m}}] ζ,  (14)
\frac{∂D}{∂t} = \frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)V}{∂λ} − \frac{1}{a} \frac{∂(f + ζ)U}{∂μ} − ∇^{2} [g(h′ + h_{s}) + E]
−K_{m} [(−1)^{N_{m}}∇^{2N_{m}} − (\frac{2}{a^{2}} )^{2N_{m}}] D,  (15)
\frac{∂h}{∂t} = −\frac{1}{a(1 − μ^{2})} \frac{∂(h′U)}{∂λ} − \frac{1}{a} \frac{∂(h′V)}{∂μ} − h̄D − (−1)^{N_{h}}K_{h}∇^{2N_{h}}h. (16)

\underline{球面浅水方程式モデル : 渦度発散型  3 水平離散化  7}

3 水平離散化

この節では方程式系の空間離散化について説明, プログラム内で実際に用いられている方程式を記述する.

格子点を経度方向には等間隔, 緯度方向にはガウス緯度上に置き, 物理量 A を格子点上の値で定義したものを A_{i,j} と表すことにする. 添字 i, j は格子点 (λ_{i}, μ_{j}) 上の値であることを示す.

A_{i,j} に対する離散スペクトル変換を Â_{m,n} であらわす. すなわち,

A_{i,j} = ∑_{n=0}^{M} ∑_{m=−n}^{n} P_{n}^{m}(μ)e^{imλ}Â_{m,n},  (17)
Â_{m,n} = \frac{1}{2I} ∑_{j=1}^{J} ∑_{i=1}^{I} A_{i,j}e^{−imλ_{i}}P_{n}^{m}(μ_{j})w_{j},  (18)

支配方程式 (14),(15), (16), をスペクトル変換すると

\frac{∂ζ̂_{m,n}}{∂t} = −  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ_{)}U}{∂λ} ]_{m,n} − ˆ [\frac{1}{a} \frac{∂(f + ζ)V}{∂μ} ]_{m,n} − K_{m_{n,Nm}} ζ̂_{m,n},(19)
\frac{∂D̂_{m,n}}{∂t} = ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)V}{∂λ} ]_{m,n} − ˆ [\frac{1}{a} \frac{∂(f + ζ)U}{∂μ} ]_{m,n}
−g∇^{2}ĥ′_{m,n} − ˆ [∇^{2} {gh_{s} + E}]_{m,n} − K_{m_{n,Nm}} D̂_{m,n},  (20)
\frac{∂ĥ_{m,n}}{∂t} = −  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(h′U)}{∂λ}]_{m,n} − ˆ [\frac{1}{a} \frac{∂(h′V)}{∂μ}]_{m,n} − h̄D̂_{m,n} − K_{h_{n,N_{h}}} ĥ_{m,n}. (21)

ここで各散逸項に対して, 球面調和函数がラプラシアンの固有関数であることを用いて ˆ (∇^{2}A)_{m,n} = −n(n + 1)Ã_{m,n}/a^{2} の関係を散逸項に対して適用し,

K_{m_{n,Nm}} = K_{m} [(−1)^{N_{m}} (\frac{−n(n + 1)}{a^{2}})^{N_{m}} − (\frac{2}{a^{2}} )^{N_{m}} ]
= K_{m} [(\frac{−[−n(n + 1)]}{a^{2}})^{N_{m}} − (\frac{2}{a^{2}} )^{N_{m}} ] ,
K_{h_{n,N_{h}}} = (−1)^{N_{h}}K_{h} (\frac{−n(n + 1)}{a^{2}})^{N_{h}} = K_{h} (\frac{−[−n(n + 1)]}{a^{2}})^{N_{h}} .

として表している.

非線形項のスペクトル変換は, 格子点上での値を用いて非線形項を評価し, それをスペクトル変換する (変換法) を適用することにする.

\underline{球面浅水方程式モデル : 渦度発散型  4 時間積分  8}

4 時間積分

ここでは時間積分法について記述し, プログラム内で実際に用いられている方程式を記述する. 以下では Δt を時間格子間隔, 時刻 t, t + Δt における ζ̂_{m,n} の値をそれぞれ ζ̂_{m,n}^{τ}, ζ̂_{m,n}^{τ+1} 等と表す.

4.1 陽解法 (オイラースキーム)

時間積分を 1 次の陽解法 (オイラースキーム) で解く場合の定式化を記す. 時間変化項を全て時間ステップ τ の値で評価し, 時間変化を τ + 1 と τ の差分で見積もれば, (19), (20), (21) はつぎのように時間方向に差分化される.

ζ̂_{m,n}^{τ+1} = ζ̂_{m,n}^{τ} + (Δt) × \{−  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ_{)}U}{∂λ} ]_{m,n} − ˆ [\frac{1}{a} \frac{∂(f + ζ)V}{∂μ} ]_{m,n} − K_{m_{n,Nm}} ζ̂_{m,n}\} ,

(22)

D̂_{m,n}^{τ+1} = D̂_{m,n}^{τ} + (Δt) × \{  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)V}{∂λ} ]_{m,n} − ˆ [\frac{1}{a} \frac{∂(f + ζ)U}{∂μ} ]_{m,n}
−g∇^{2}ĥ′_{m,n} − ˆ [∇^{2} {gh_{s} + E}]_{m,n} − K_{m_{n,Nm}} D̂_{m,n}} ,

(23)

ĥ_{m,n}^{τ+1} = ĥ_{m,n}^{τ} + (Δt) × \{−  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(h′U)}{∂λ}]_{m,n} − ˆ [\frac{1}{a} \frac{∂(h′V)}{∂μ}]_{m,n} − h̄D̂_{m,n} − K_{h_{n,N_{h}}} ĥ_{m,n}\} .

(24)

4.2 セミインプリシット & リープフロッグスキーム

時間積分にセミインプリシット法を用いる場合の定式化を行なう. 重力波に関係する線形項は陰的に解き, その他の項は陽的に解く (Krishnamurti, et al., 1998).

(19) の時間積分は陽的解法を用いて行う. 摩擦項に対しては Euler 法, その他の項に対しては Leapfrog 法を適用する.

ζ̂_{m,n}^{τ+1} = ζ̂_{m,n}^{τ−1} + 2Δt × \{−  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ_{)}U}{∂λ} ]_{m,n}^{τ}

−  ˆ [\frac{1}{a} \frac{∂(f + ζ)V}{∂μ} ]_{m,n}^{τ} − K_{m_{n,Nm}} ζ̂_{m,n}^{τ−1}\} .  (25)

(20), (21) は時間方向に離散化する前に重力波成分の分離を行う. (20), (21) を書き直すと

\frac{∂D̂_{m,n}}{∂t} = (\frac{∂D̂_{m,n}}{∂t} )_{NG} − g∇^{2}ĥ′_{m,n} − K_{m_{n,Nm}} D̂_{m,n}  (26)
\frac{∂ĥ′_{m,n}}{∂t} = (\frac{∂ĥ′_{m,n}}{∂t} )_{NG} − \overline{h}D̂_{m,n} − K_{h_{n,N_{h}}} ĥ′_{i,j}  (27)

(26), (27) のそれぞれ右辺第 2 項が重力波に関係する線形項である. ここで,

(\frac{∂D̂_{m,n}}{∂t} )_{NG} (ζ, D, U, V)

= ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)V}{∂λ} ]_{m,n} − ˆ [\frac{1}{a} \frac{∂(f + ζ)U}{∂μ} ]_{m,n} − ˆ [∇^{2}(gh_{s} + E)]_{m,n}, (28)
(\frac{∂ĥ′_{m,n}}{∂t} )_{NG} (ζ, D, U, V) = −  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(h′U)}{∂λ}]_{m,n} − ˆ [\frac{1}{a} \frac{∂(h′V)}{∂μ}]_{m,n} (29)

は重力波に関係しない (非線形) 項である.

セミインプリシット法を適用するために, つぎの時間平均演算子を導入する.

D^{τ∗} = \frac{D^{τ+1} + D^{τ−1}}{2} ,  (30)

(26), (27) の時間変化を中央差分で評価する. 一方右辺は重力波に関係する項を τ + 1 と τ − 1 の時間平均の値, 摩擦項はオイラースキームを適用すべく τ − 1 での値, それ以外の非線形項を τ での値で評価する. さらに左辺の時間変化項を上の演算子を用いて表せば,

\frac{D̂_{m,n}^{τ+1} − D̂_{m,n}^{τ−1}}{2Δt} = \frac{D̂_{m,n}^{τ∗} − D̂_{m,n}^{τ−1}}{Δt} = (\frac{∂D_{m,n}}{∂τ} )_{NG}^{τ} − g∇^{2}ĥ_{m,n}^{′τ∗} − K_{m_{n,Nm}} D̂_{m,n}^{τ−1},
\frac{ĥ_{m,n}^{′τ+1} − ĥ_{m,n}^{′τ−1}}{2Δt} = \frac{ĥ_{m,n}^{′τ∗} − ĥ_{m,n}^{′τ−1}}{Δt} = (\frac{∂h′_{m,n}}{∂τ} )_{NG}^{τ} − \overline{h}D̂_{m,n}^{τ∗} − K_{h_{n,N_{h}}} ĥ_{n,m}^{′τ−1},

すなわち

D̂_{m,n}^{τ∗} = D̂_{m,n}^{τ−1} + (Δt) [(\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} − g∇^{2}ĥ_{m,n}^{′τ∗} − K_{m_{n,Nm}} D̂_{m,n}^{τ−1}] , . (31)
ĥ_{m,n}^{τ∗} = ĥ_{m,n}^{′τ−1} + (Δt) [(\frac{∂ĥ′_{m,n}}{∂t} )_{NG}^{τ} − \overline{h}D̂_{m,n}^{τ∗} − K_{h_{n,N_{h}}} ĥ_{n,m}^{′τ−1}] ,  (32)

spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

(31), (32) から ĥ_{m,n}^{τ∗} を消去する. そのために

ĥ_{m,n}^{τ+} = ĥ_{m,n}^{′τ−1} + (Δt) [(\frac{∂ĥ′_{m,n}}{∂t} )_{NG}^{τ} − K_{h_{n,N_{h}}} ĥ_{n,m}^{′τ−1}]  (33)
= [1 − (Δt)K_{h_{n,N_{h}}}]ĥ_{m,n}^{′τ−1} + (Δt) (\frac{∂ĥ′_{m,n}}{∂t} )_{NG}^{τ}  (34)

と定義すると,

ĥ_{m,n}^{τ∗} = ĥ_{m,n}^{τ+} − (Δt)\overline{h}D̂_{m,n}^{τ∗},  (35)

であり,

D̂_{m,n}^{τ∗} = D̂_{m,n}^{τ−1} + (Δt) [(\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} − K_{m_{n,Nm}} D̂_{m,n}^{τ−1}] − (Δt)g∇^{2}ĥ_{m,n}^{′τ∗}
= D̂_{m,n}^{τ−1} + (Δt) [(\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} − K_{m_{n,Nm}} D̂_{m,n}^{τ−1}] − (Δt)g∇^{2}[ĥ_{m,n}^{τ+} − (Δt)\overline{h}D̂_{m,n}^{τ∗}],
D̂_{m,n}^{τ∗} − (Δt)^{2}g\overline{h}∇^{2}D̂_{m,n}^{τ∗}
= (1 − K_{m_{n,Nm}}Δt)D̂_{m,n}^{τ−1} + (Δt) (\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} − (Δt)g∇^{2}ĥ_{m,n}^{τ+}.

左辺に ˆ (∇^{2}A)_{m,n} = −n(n + 1)Ã_{m,n}/a^{2} の関係を再び用いることで D̂_{m,n}^{τ∗} について解くことができて,

D̂_{m,n}^{τ∗} = \frac{1}{1 + (Δt)^{2}g\overline{h}\frac{n(n+1)}{a^{2}}} \{(1 − K_{m_{n,Nm}}Δt)D̂_{m,n}^{τ−1} + (Δt) (\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} − (Δt)g∇^{2}ĥ_{m,n}^{τ+}\}

(36)

これを (35) に代入して ĥ_{m,n}^{τ∗} が求められる.

それぞれの時間平均が求められたら, 時間平均の定義式 (30) からつぎの時間ステップの値が求められる.

以上の手順をまとめると,

1. 重力波以外の時間変化項から h′ を見積もる.

ĥ_{m,n}^{τ+} = [1 − (Δt)K_{h_{n,N_{h}}}]ĥ_{m,n}^{′τ−1} + (Δt) (\frac{∂ĥ′_{m,n}}{∂t} )_{NG}^{τ}  (37)

2. 発散 D に関してセミインプリシット時間積分を実行する.

D̂_{m,n}^{τ∗} = \frac{1}{1 + (Δt)^{2}g\overline{h} · n(n + 1)/a^{2}}
× \{(1 − K_{m_{n,Nm}}Δt)D̂_{m,n}^{τ−1} + (Δt) (\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} − (Δt)g∇^{2}ĥ_{m,n}^{τ+}\} . (38)

3. 変位擾乱 h′ に対してセミインプリシット時間積分を実行する.

ĥ_{m,n}^{τ∗} = ĥ_{m,n}^{τ+} − (Δt)\overline{h}D̂_{m,n}^{τ∗},  (39)

4. つぎの時間ステップの値を求める

D_{m,n}^{τ+1} = 2D_{m,n}^{τ∗} − D_{m,n}^{τ−1},  (40)
h_{m,n}^{′τ+1} = 2h_{m,n}^{′τ∗} − h_{m,n}^{′τ−1}.  (41)

4.2.1 時間フィルター

Leap frog スキームを用いる場合には, 計算モードの増幅を抑制するため Asselin (1972) の時間フィルターを適用する. フィルターの定義式は

A^{τ∗} = A^{τ} + \frac{1}{2}γ(A^{τ+1} − 2A^{τ} + A^{τ−1})  (42)

ここで A^{τ∗} が修正された τ での A の値, γ はフィルターの係数である. この時間フィルターを, とあるタイムステップ間隔で適用する.

4.3 リープフロッグ & クランク・ニコルソン法

時間積分にリープフロッグ法とクランク・ニコルソン法を組み合わせて用いる場合の定式化を行なう^{1}.

(19) の時間積分は陽的解法を用いて行う. 摩擦項に対しては Euler 法, その他の項に対しては Leapfrog 法を適用する.

ζ̂_{m,n}^{τ+1} = ζ̂_{m,n}^{τ−1} + 2Δt × \{−  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ_{)}U}{∂λ} ]_{m,n}^{τ}
−  ˆ [\frac{1}{a} \frac{∂(f + ζ)V}{∂μ} ]_{m,n}^{τ} − K_{m_{n,Nm}} ζ̂_{m,n}^{τ−1}\} .  (43)

1このスキームは第 4.2 節でのものと計算手順が異なるだけで, 本質的にはおなじであると考えられる (要確認).

spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

(26), (27) には, 重力波に関係するそれぞれの式の右辺第 2 項に対して台形型の陰解法 (Clank-Nicolson 法), 摩擦項に対してはオイラー法, その他の項に対してはリープフロッグ法を適用する.

\frac{D̂_{m,n}^{τ+1} − D̂_{m,n}^{τ−1}}{2Δt} = (\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} + \frac{n(n + 1)}{a^{2}}g \frac{ĥ_{m,n}^{′τ+1} + ĥ_{m,n}^{′τ−1}}{2} − K_{m_{n,Nm}} D̂_{m,n}^{τ−1},(44)
\frac{ĥ_{m,n}^{′τ+1} − ĥ_{m,n}^{′τ−1}}{2Δt} = (\frac{∂ĥ′_{m,n}}{∂t} )_{NG}^{τ} − \overline{h}\frac{D̂_{m,n}^{τ+1} + D̂_{m,n}^{τ−1}}{2} − K_{h_{n,N_{h}}} ĥ_{n,m}^{′τ−1}.  (45)

これらを D̂_{m,n}^{τ+1}, ĥ_{m,n}^{′τ+1} について解くと以下の式が得られる (Krishnamurti et al., 1998; Satoh, 2004).

(1 + Δt^{2}g\overline{h}\frac{n(n + 1)}{a^{2}}) D̂_{m,n}^{τ+1} = (1 − Δt^{2}g\overline{h}\frac{n(n + 1)}{a^{2}} − 2ΔtK_{m_{n,Nm}} ) D̂_{m,n}^{τ−1}
+2ΔtG_{1}(ζ, D, U, V),  (46)
(1 + Δt^{2}g\overline{h}\frac{n(n + 1)}{a^{2}}) ĥ_{m,n}^{′τ+1} = (1 − Δt^{2}g\overline{h}\frac{n(n + 1)}{a^{2}} − 2ΔtK_{h_{n,N_{h}}} ) ĥ_{m,n}^{′τ−1}
+2ΔtG_{2}(ζ, D, U, V)  (47)

ここで,

G_{1}(ζ, D, U, V) = (\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} + \frac{n(n + 1)}{a^{2}}g [(1 − ΔtK_{h_{n,p}})ĥ_{m,n}^{′τ−1} + (\frac{∂ĥ′_{m,n}}{∂t} )_{NG}^{τ} (Δt)] ,

(48)

G_{2}(ζ, D, U, V) = (\frac{∂ĥ′_{m,n}}{∂t} )_{NG}^{τ} − \overline{h}(1 − ΔtK_{m_{n,p}})D̂_{m,n}^{τ−1} − \overline{h} (\frac{∂D̂_{m,n}}{∂t} )_{NG}^{τ} (Δt).  (49)

第 4.2 節と同様, 計算モードを抑えるための時間フィルター (第 4.2.1 節) を, とあるタイムステップ間隔で適用する.

4.4 ルンゲ・クッタ & クランク・ニコルソン法

時間積分の数値スキームとして, 散逸項にはクランク・ニコルソン法, それ以外の項を 4 次のルンゲ・クッタ法を組み合わせて用いる場合の定式化を行なう.

まず散逸項以外の項に関してルンゲ・クッタ法で積分する. (19), (20), (21) から散逸項以外の時間変化を

(\frac{∂ζ̂_{m,n}}{∂t} )_{ND} [ζ̂_{m,n}, D̂_{m,n}, ĥ′_{m,n}] = −  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)U}{∂λ} ]_{m,n} − ˆ [\frac{1}{a} \frac{∂(f + ζ)V}{∂μ} ]_{m,n}

(50)

(\frac{∂D̂_{m,n}}{∂t} )_{ND} [ζ̂_{m,n}, D̂_{m,n}, ĥ′_{m,n}] = ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(f + ζ)V}{∂λ} ]_{m,n} − ˆ [\frac{1}{a} \frac{∂(f + ζ)U}{∂μ} ]_{m,n}
−g∇^{2}ĥ_{m,n}^{′τ} − ˆ [∇^{2} {gh_{s} + E}]_{m,n}  (51)
(\frac{∂ĥ_{m,n}}{∂t} )_{ND} [ζ̂_{m,n}, D̂_{m,n}, ĥ′_{m,n}] = −  ˆ [\frac{1}{a(1 − μ^{2})} \frac{∂(h′U)}{∂λ}]_{m,n} − ˆ [\frac{1}{a} \frac{∂(h′V)}{∂μ}]_{m,n} − h̄D̂_{m,n}

(52)

と定義する. ここで非線形項で用いられている格子点で定義された変数は各スペクトル変数 ζ̂_{m,n}, D̂_{m,n}, ĥ′_{m,n} からスペクトル逆変換して求められる.

ルンゲ・クッタ法による積分の手順はつぎのようになる. 渦度方程式だけを例に手順を記す. 第 1 段は

(\frac{∂ζ̂_{m,n}}{∂t} )^{k1} = (\frac{∂ζ̂_{m,n}}{∂t} )_{ND} [ζ̂_{m,n}^{τ}, D̂_{m,n}^{τ}, ĥ′_{m,n}^{τ}],  (53)
ζ̂_{m,n}^{τ+1/2} = ζ̂_{m,n}^{τ} + (\frac{∂ζ̂_{m,n}}{∂t} )^{k1} · \frac{(Δt)}{2},  (54)

第 2 段は

(\frac{∂ζ̂_{m,n}}{∂t} )^{k2} = (\frac{∂ζ̂_{m,n}}{∂t} )_{ND} [ζ̂_{m,n}^{τ+1/2}, D̂_{m,n}^{τ+1/2}, ĥ′_{m,n}^{τ+1/2}],  (55)
ζ̂_{m,n}^{∗τ+1/2} = ζ̂_{m,n}^{τ} + (\frac{∂ζ̂_{m,n}}{∂t} )^{k2} · \frac{(Δt)}{2},  (56)

第 3 段は

(\frac{∂ζ̂_{m,n}}{∂t} )^{k3} = (\frac{∂ζ̂_{m,n}}{∂t} )_{ND} [ζ̂_{m,n}^{∗τ+1/2}, D̂_{m,n}^{∗τ+1/2}, ĥ′_{m,n}^{∗τ+1/2}],  (57)
ζ̂_{m,n}^{∗τ+1} = ζ̂_{m,n}^{τ} + (\frac{∂ζ̂_{m,n}}{∂t} )^{k3} · (Δt),  (58)

第 4 段は

(\frac{∂ζ̂_{m,n}}{∂t} )^{k4} = (\frac{∂ζ̂_{m,n}}{∂t} )_{ND} [ζ̂_{m,n}^{∗τ+1}, D̂_{m,n}^{∗τ+1}, ĥ′_{m,n}^{∗τ+1}],  (59)

最終的に

ζ̂_{m,n}^{τ+1,ND} = ζ̂_{m,n}^{τ} + (Δt) [\frac{1}{6} (\frac{∂ζ̂_{m,n}}{∂t} )^{k1} + \frac{1}{3} (\frac{∂ζ̂_{m,n}}{∂t} )^{k2} + \frac{1}{3} (\frac{∂ζ̂_{m,n}}{∂t} )^{k3} + \frac{1}{6} (\frac{∂ζ̂_{m,n}}{∂t} )^{k4}] ,

(60)

以上の手順を発散 D̂_{m,n}, 変位擾乱 ĥ′_{m,n} についても平行して計算していくことで, 非散逸項の積分が完了する.

ルンゲ・クッタ法による積分ののち, 散逸項からの修正を行なう. 渦度を例に取ると,

\frac{ζ_{m,n}^{τ+1} − ζ_{m,n}^{τ+1,ND}}{Δt} = −K_{m_{n,Nm}} \frac{ζ_{m,n}^{τ+1} + ζ_{m,n}^{τ+1,ND}}{2}

したがって

ζ_{m,n}^{τ+1} = \frac{1 − K_{m_{n,Nm}}(Δt)/2}{1 + K_{m_{n,Nm}}(Δt)/2} · ζ_{m,n}^{τ+1,ND},  (61)

と修正される. 発散 D̂_{m,n}, 変位擾乱 ĥ′_{m,n} についても同様に計算する.

\underline{{球面浅水方程式モデル : 渦度発散型 \atop 5 使用モジュールとその他の設定 15}}

5 使用モジュールとその他の設定

スペクトル変換と逆変換, 微分演算は SPMODEL ライブラリ (spml) の w module に含まれる関数を用いて行う. スペクトル変換と逆変換の際の数値積分は台形公式を用いて行う. spml が下位で使用する ISPACK の仕様から, 格子点数 I は偶数で, かつ I/2 = 2^{a}3^{b}5^{c} (a, b, c は 0 または整数) でなければならない. 非線形項の計算によって生じるエリアジングを防ぐため, 格子点数 I, J は切断波数 M に対し I > 3M + 1, J > 3M/2 を満たすように与える.

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 16

6 Williamson et al. (1992) の標準実験

作成したプログラムのチェックのために, Williamson et al. (1992) によって提案された球面浅水モデルの標準実験を行なった. 以下ではその実験設定について説明する.

6.1 テストケース 1: 極を越える cos 型の山の移流

テストケース 1 は剛体回転の流れ場に cosine 型の山を置き移流させる実験である. これは球面浅水方程式系において, u, v を定常とし h の式を移流方程式として扱う. u, v の初期値は

u = u_{0}(cos ϕ cos α + sin ϕ cos λ sin α),  (62)
v = −u_{0} sin λ sin α  (63)

と与える. 流線関数と速度ポテンシャルの初期値は

ψ = −au_{0}(sin ϕ cos α − cos λ cos ϕ sin α),  (64)
χ = 0  (65)

である. ここで α は球面上での流れの方向を指定するパラメータである. 初期の山の分布は

h = \{ \begin{matrix}(h_{0}/2)(1 + cos(πr/R)) & (r < R) \\ 0 & (r ≥ R) \end{matrix}  (66)

とする. ここで

r = a arccos(sin ϕ_{c} sin ϕ + cos ϕ_{c} cos ϕ cos(λ − λ_{c}))  (67)

である. 与えたパラメータは表 6.1 にまとめる. 計算時間は 12 日間で 1 ステップ毎のデータを出力する.

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 17

表 2: テストケース 1 で用いたパラメータ\begin{matrix}パラメータ & 値 \\ α & 0, 0.05, π/2 − 0.05, π/2 \\ u_{0} & 2πa/(12 日) msec^{−1} \\ h_{0} & 1000 m \\ λ_{c} & 3π/2 \\ ϕ_{c} & 0 \\ R & a/3 \end{matrix}

6.2 テストケース 2: 定常な非線形帯状地衡流

テストケース 2 では非線形球面浅水方程式系の定常解として得られる剛体回転の流れ場と, それに地衡流平衡する厚さ分布を初期に与える.

u = u_{0}(cos ϕ cos α + sin ϕ cos λ sin α),  (68)
v = −u_{0} sin λ sin α,  (69)
h = h_{0} − (aΩu_{0} + \frac{u_{0}^{2}}{2} ) (− cos λ cos ϕ sin α + sin ϕ cos α)^{2}/g.  (70)

流線関数と速度ポテンシャルの初期値はケース 1 と同様に

ψ = −au_{0}(sin ϕ cos α − cos λ cos ϕ sin α),  (71)
χ = 0  (72)

である.

この計算ではコリオリパラメータ f は球面上での流れの方向を指定するパラメータ α に依存する.

f = 2Ω(− cos λ cos ϕ sin α + sin ϕ cos α).  (73)

与えたパラメータは表 6.2 にまとめる. 計算時間は 5 日間で, テストケース 1 と同様に 1 ステップ毎のデータを出力する.

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 18

表 3: テストケース 2 で用いたパラメータ\begin{matrix}パラメータ & 値 \\ α & 0, 0.05, π/2 − 0.05, π/2 \\ u_{0} & 2πa/(12 日) msec^{−1} \\ h_{0} & 2.94 × 10^{4}/g m \end{matrix}

6.3 テストケース 3: 局在した非線形帯状地衡流

この実験はケース 2 と同様に定常解として得られる帯状な流れを与える. この解では限られた場所でのみ流れが存在する. 流れの軸を変えた場合の初期値の表現が難しいので, まず流れが経度座標と平行な系 (λ′, ϕ′) での初期値を u′, v′, h′ と表し, 座標軸を α だけ回転した系 (λ, ϕ) における u, v, h へ変換する.

変換前の u′, v′, h′ は以下のように与えられる.

u′ = u_{0}b(x)b(x_{e} − x)e^{4/x_{e}},  (74)
v′ = 0,  (75)
h′ = h_{0} − \frac{a}{g} ∫_{−π/2}^{ϕ′} (2Ω sin τ + \frac{u′(τ) tan τ}{a} ) u′(τ) dτ.  (76)

ここで

b(x) = \{ \begin{matrix}0 & x ≤ 0 \\ e^{−x^{−1}} & 0 < x \end{matrix}

x = x_{e} \frac{ϕ′ − ϕ_{b}}{ϕ_{e} − ϕ_{b}}

である.

u′, v′, h′ から u, v, h への変換は以下のように行う. まず座標は

sin ϕ′ = sin ϕ cos α − cos α cos λ sin α,  (77)
sin λ′ cos θ′ = sin λ cos ϕ  (78)

にしたがって変換する. λ′ を決める際には

λ′ = \{ \begin{matrix}λ′_{p} & cos α cos λ cos ϕ + sin α sin ϕ ≥ 0 \\ π − λ′_{p} & cos α cos λ cos ϕ + sin α sin ϕ < 0 \end{matrix}

とする. ここで λ′_{p} は λ の主値である. u, v は

v cos ϕ = −u′ sin α sin λ′,  (79)
u cos λ = v sin ϕ sin λ + u′ cos λ′.  (80)

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 19

より求める. h は h′ の式を回転後の座標系で数値積分することで得られる. コリオリパラメータ f はテストケース 3 の場合と同様に

f = 2Ω(− cos λ cos ϕ sin α + sin ϕ cos α)  (81)

とする. 与えたパラメータは表 6.3 にまとめる. 計算時間は 5 日間で, テストケース 1 と同様に 1 ステップ毎のデータを出力する.

表 4: テストケース 3 で用いたパラメータ\begin{matrix}パラメータ & 値 \\ α & 0, π/3 \\ u_{0} & 2πa/(12 日) msec^{−1} \\ h_{0} & 2.94 × 10^{4}/g m \\ ϕ_{b} & −pi/6 \\ ϕ_{e} & pi/2 \\ x_{e} & 0.3 \end{matrix}

spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 20

6.4 テストケース 4: 移動する低気圧を持つ強制された流れ

このケースでは u, v, h を予報変数とする形で表された浅水方程式系に強制項を加える. 強制項は

F_{u} = \frac{dũ}{dt} − \frac{ũṽ tan ϕ}{a} − fṽ + \frac{g}{a cos ϕ} \frac{∂h̃}{∂λ},  (82)
F_{v} = \frac{dṽ}{dt} + \frac{ũṽ tan ϕ}{a} + fṽ + \frac{g}{a} \frac{∂h̃}{∂λ},  (83)
F_{h} = \frac{dh̃}{dt} + \frac{h̃}{a cos ϕ} [\frac{∂ũ}{∂λ} + \frac{∂ṽ cos ϕ}{∂ϕ}]  (84)

として与えられる. ここで

\frac{d}{dt} = \frac{∂}{∂t} + \frac{u}{a cos ϕ} \frac{∂}{∂λ} + \frac{v}{a} \frac{∂}{∂ϕ}

である. 渦度方程式 (1) と発散方程式 (2) における強制 F_{ζ}, F_{D} は

F_{ζ} = \frac{1}{a(1 − μ^{2})} \frac{∂(F_{v} cos ϕ)}{∂λ} − \frac{1}{a} \frac{∂(F_{u} cos ϕ)}{∂μ},  (85)
F_{D} = \frac{1}{a(1 − μ^{2})} \frac{∂(F_{u} cos ϕ)}{∂λ} + \frac{1}{a} \frac{∂(F_{v} cos ϕ)}{∂μ}  (86)

より計算される. ũ, ṽ, h̃ は

ũ = \overline{u} − \frac{1}{a} \frac{∂\overline{ψ}}{∂ϕ} ,  (87)
ṽ = \frac{1}{a cos ϕ} \frac{∂\overline{ψ}}{∂λ} ,  (88)
h̃ = \overline{h} + \frac{f}{g} \overline{ψ}  (89)

で与えられ, ここで

\overline{u} = u_{0} sin^{14}(2ϕ),  (90)
\overline{h} = h_{0} − \frac{a}{g} ∫_{−π/2}^{ϕ} (2Ω sin τ + \frac{\overline{u}(τ) tan τ}{a} ) \overline{u}(τ) dτ,  (91)
\overline{ψ} = ψ_{0} exp (−σ \frac{1 − C}{1 + C} ) ,  (92)
ψ_{0} = −0.03\frac{gh_{0}}{2Ω sin(π/4)},  (93)
C = sin ϕ_{0} sin ϕ + cos ϕ_{0} cos ϕ (λ − \frac{u_{0}}{a} t − λ_{0})  (94)

spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 21

である. 与えたパラメータは表 6.4 にまとめる. 計算時間は 5 日間で, テストケース 1 と同様に 1 ステップ毎のデータを出力する.

表 5: テストケース 4 で用いたパラメータ\begin{matrix}パラメータ & 値 \\ u_{0} & 20, 40 (msec^{−1}) \\ h_{0} & 10^{5}/g m \\ σ & (12.74244)^{2} \\ λ_{0} & 0 \\ ϕ_{0} & π/4 \end{matrix}

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 22

6.5 テストケース 5: 孤立した山を越える帯状流

これは中緯度に孤立した山を置き, それによる流体の応答を調べる実験である. 底面の地形分布 h_{s} は以下のように与える.

h_{s} = h_{s0} (1 − \frac{1}{R}\sqrt{min[R^{2}, (λ − λ_{c})^{2} + (ϕ − ϕ_{c})^{2}]})  (95)

初期の場はケース 2 の α = 0 の場合と同じである.

u = u_{0} cos ϕ,  (96)
v = 0,  (97)
h = h_{0} − (aΩu_{0} + \frac{u_{0}^{2}}{2} ) (sin ϕ)^{2}/g.  (98)

流線関数と速度ポテンシャルの初期値は

ψ = −u_{0}a sin θ,  (99)
χ = 0  (100)

である. これより渦度と発散の初期値は

ζ = \frac{2u_{0}}{a} sin θ,  (101)
D = 0  (102)

である.

与えたパラメータは表 6.6 にまとめる. 計算時間は 15 日間で, テストケース 1 と同様に 1 ステップ毎のデータを出力する.

表 6: テストケース 5 で用いたパラメータ\begin{matrix}パラメータ & 値 \\ u_{0} & 20 msec^{−1} \\ h_{0} & 5960 m \\ h_{s0} & 2000 m \\ R & π/9 \\ λ_{c} & π/2 \\ ϕ_{c} & π/6 \end{matrix}

spshallow-zd.tex  2004 年 10 月 7 日 (竹広真一)

球面浅水方程式モデル : 渦度発散型

6 Williamson et al. (1992) の標準実験 23

6.6 テストケース 6: ロスビー・ハウルヴィッツ波

初期に非発散順圧渦度方程式の解であるロスビー・ハウルヴィッツ波を表現する流線関数を与える.

ψ = −a^{2}ω sin ϕ + a^{2}K cos^{R} ϕ cos Rλ,  (103)
χ = 0.  (104)

初期の速度場は

u = aω cos ϕ + aK cos^{R−1} ϕ(R sin^{2} ϕ − cos^{2} ϕ) cos Rλ,  (105)
v = −aKR cos^{R−1} ϕ sin ϕ sin Rλ  (106)

であり, 渦度と発散の初期値は

ζ = 2ω sin ϕ − K sin ϕ cos^{R} ϕ(R^{2} + 3R + 2) cos Rλ,  (107)
D = 0  (108)

である.

初期の h は

h = h_{0} + [a^{2}A(ϕ) + a^{2}B(ϕ) cos Rλ + a^{2}C(ϕ) cos 2Rλ]/g  (109)

で与えられる. ここで

A(ϕ) = \frac{ω}{2} (2Ω + ω) cos^{2} ϕ
+\frac{1}{4}K^{2} cos^{2R}[(R + 1) cos^{2} ϕ + (2R^{2} − R − 2) − 2R^{2}cos^{−2}ϕ], (110)
B(ϕ) = \frac{2(Ω + ω)K}{(R + 1)(R + 2)} cos^{R} ϕ[(R^{2} + 2R + 2) − (R + 1)^{2} cos^{2} ϕ], (111)
C(ϕ) = \frac{1}{4}K^{2} cos^{2R} ϕ[(R + 1) cos^{2} ϕ − (R + 2)]  (112)

である.

与えたパラメータは表 6.6 にまとめる. 計算時間は 15 日間で, テストケース 1 と同様に 1 ステップ毎のデータを出力する.

\underline{球面浅水方程式モデル : 渦度発散型  7 参考文献  24}

表 7: テストケース 5 で用いたパラメータ\begin{matrix}パラメータ & 値 \\ h_{0} & 8000 m \\ R & 4 \\ ω & 7.848 × 10^{−6} sec^{−1} \\ K & 7.848 × 10^{−6} sec^{−1} \end{matrix}

7 参考文献

Asselin, R. A., 1972: Frequency filter for time integrations. Mon. Wea. Rev., 100, 487–490.

石岡圭一, 2004: スペクトル法による数値計算入門, 東京大学出版会, 232pp.

Krishnamurti, T. N., H. S. Bedi, and V. M. Hardiker, 1998: An introduction to global spectral modeling. Ocford Univ. Press, 253 pp.

竹広真一, 石岡圭一, 柿並義宏, 西澤誠也, 森川靖大, 小高正嗣, 石渡正樹, 林祥介, SPMODEL 開発グループ, 2005: 階層的地球流体力学スペクトルモデル集 (SPMODEL), http://www.gfd-dennou.org/library/spmodel/, 地球流体電脳倶楽部.

Williamson, D. L., J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber, 1992: A Standard Test Set for Numerical Approximations to the Shallow Water Equation in Spherical Geometry. J. Comput. Phys., 102, 211–224.

\underline{球面浅水方程式モデル : 渦度発散型謝辞 25}

謝辞

本資源は, 地球流体電脳倶楽部のインターネット上での学術知識の集積と活用の実験の一環として

http://www.gfd-dennou.org/arch/spmodel/

において公開されているものである (\begin{matrix}c \\ ◯ \end{matrix}地球流体電脳倶楽部スペクトルモデルプロジェクト spmodel@gfd-dennou.org 2002. ). 本資源は, 著作者の諸権利に抵触しない (迷惑をかけない) 限りにおいて自由に利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめることをお願いする (無保証無責任原則).

本資源に含まれる元資源提供者(図等の版元等を含む)からは, 直接的な形での WEB 上での著作権または使用許諾を得ていない場合があるが, 勝手ながら, 「未来の教育」のための実験という学術目的であることをご理解いただけるものと信じ, 学術標準の引用手順を守ることで諸手続きを略させていただいている. 本資源の利用者には, この点を理解の上, 注意して扱っていただけるようお願いする. 万一, 不都合のある場合には

spmodel@gfd-dennou.org

まで連絡していただければ幸いである.