DCPAM4 第2部離散化

地球流体電脳倶楽部

2008 年 6 月 10 日 (DCPAM4-20080609-1)

\underline{ii DCPAM4 第 2 部離散化}

目次

1  この文書について  1

1.1 この文書について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2  座標系・変換公式  2

2.1 座標系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1.1 水平格子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1.2 鉛直レベル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.2 水平スペクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2.1 水平スペクトルの基底の導入 . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2.2 波数切断 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2.3 離散化したスペクトルの基底の直交性 . . . . . . . . . . . . . . . . . . . . 5 2.2.4 格子点値とスペクトルの係数との変換法 . . . . . . . . . . . . . . . . . . . 5 2.2.5 内挿公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2.6 空間微分の評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3  支配方程式・力学過程  7

3.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.2 鉛直離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.2.1 連続の式, 鉛直速度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2.2 静水圧の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2.3 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2.4 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2.5 水蒸気の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.3 水平離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.1 連続の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.2 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.3 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3.4 水蒸気の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.4 時間積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.4.1 leap frog による時間積分と時間フィルター . . . . . . . . . . . . . . . . . . 13 3.4.2 semi-implicit 時間積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4  参考文献  16

4.1 座標系・変換公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.2 力学過程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.3 雲 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.4 放射 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2008/06/10(地球流体電脳倶楽部)  discretization.tex

4.5 鉛直輸送 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.6 地表面過程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

A  座標系・変換公式  1

A.1 球面調和函数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 A.1.1 定義と性質 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 A.1.2 球面調和函数の空間微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 A.1.3 コメント — 全波数について . . . . . . . . . . . . . . . . . . . . . . . . . . 5 A.1.4 グラフ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

A.2 微分公式, ＧＣＭの変数の微分関係式 . . . . . . . . . . . . . . . . . . . . . . . . . 7 A.2.1 スカラー量の微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 A.2.2 ベクトル量の微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 A.2.3 発散 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 A.2.4 渦度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 A.2.5 速度ポテンシャル, 流線関数と (U, V) . . . . . . . . . . . . . . . . . . . . . 8

A.3 Legendre 函数 P_{n} の性質 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 A.3.1 多項式と Legendre 函数の積の積分 . . . . . . . . . . . . . . . . . . . . . . 9 A.3.2 Legendre 函数の零点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

A.4 積分評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 A.4.1 Gauss の台形公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 A.4.2 Gauss-Legendre の公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

A.5 球面調和函数の離散的直交関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 A.6 スペクトルの係数と格子点値とのやり取り . . . . . . . . . . . . . . . . . . . . . . 17

A.6.1 スペクトルの係数と格子点値との値のやり取り . . . . . . . . . . . . . . . . 17 A.6.2 スペクトルの係数と格子点値との値のやり取り〜東西微分編 . . . . . . . . 17 A.6.3 スペクトルの係数と格子点値との値のやり取り〜南北微分編 . . . . . . . . 18 A.6.4 速度の格子点値から発散・渦度のスペクトルの係数への変換 . . . . . . . . 21 A.6.5 χ, ψ のスペクトルの係数から速度の格子点値への変換 . . . . . . . . . . . . 21

A.7 スペクトルの係数同士の関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 A.8 波数切断 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

A.8.1 波数切断の仕方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 A.8.2 切断波数の決め方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

A.9 スペクトルモデルと差分モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

B  謝辞  31
B  謝辞  31

B.1 ライセンス規定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 B.2 使用上の注意 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 B.3 開発グループメンバー . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

B.3.1 2007 年度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

discretization.tex  2008/06/10(地球流体電脳倶楽部)

第1章この文書について

1.1 この文書について

この文書は, 地球流体電脳倶楽部惑星大気モデル (DCPAM) の離散化方法を解説したものである. DCPAM は開発中であり, 本文書の内容とソースコードとは必ずしも一致しない.

2005/04/04 石渡正樹

第2章座標系・変換公式

2.1 座標系

ここでは水平格子点, 鉛直レベルのとり方を記す. さらに, 力学過程の時間積分において使用する水平スペクトルを定義し, 格子点値とスペクトルの係数との変換則を記す.

2.1.1 水平格子

水平方向の格子点の位置は, Gauss 緯度（格子点数 J 個^{1}） , 等間隔の経度（同 I 個）である.

• Gauss 緯度 Gauss 緯度を J 次の Legendre 函数 P_{J}(sin ϕ) の零点 ϕ_{j}（j = 1, 2, 3, ⋯, J）として定義する. 順番としては, \frac{π}{2} > ϕ_{1} > ϕ_{2} > ⋯ > ϕ_{J} > −\frac{π}{2} とする^{2}. なお以後, sin ϕ = μ と書くことがある.

• 経度方向の格子点経度方向の格子点の位置を

λ_{i} = \frac{2π（i − 1)}{I} （i = 1, 2, ⋯, I）  (2.2)

ととる.

2.1.2 鉛直レベル

Arakawa and Suarez(1983) のスキームを用いる. とり方は以下のとおりである^{3}.

下の層から上へと層の番号をつける. 整数レベルと半整数レベルを定義する^{4}. 半整数レベルでの σ の値 σ_{k−1/2}(k = 1, 2, ⋯, K) を定義する. ただし, レベル\frac{1}{2} は下端 (σ = 1), レベル K + \frac{1}{2} は上端 (σ = 0) とする. 整数レベルの σ の値 σ_{k}(k = 1, 2, … K) は次の式から求める.

[\frac{d}{dμ} \{(1 − μ^{2}) \frac{d}{dμ} \} + J(J + 1)] P_{J} (μ) = 0  (2.1)

91/10/31 保坂征宏

1以下, J は偶数とする. 現在の東大版大気 GCM では, （Gauss 緯度としてとる場合には） J は偶数でなければならない.

2J 次の Legendre 函数 P_{J} (μ) はを満たす J 次多項式であり, P_{J} (μ) の零点は全て −1 < μ < 1 にある. その理由について詳しくは「Legendre 函数 P_{n} の性質」参照.

なお, Gauss 緯度は近似的には sin^{−1} (cos \frac{j − 1/2}{J} π) で与えられる.

3なぜこうするとよいのかについては未調査. 保存量と関連するらしい.

4物理量により, 整数レベルで定義されるものと, 半整数レベルで定義されるものがある.

\underline{DCPAM4 第 2 部離散化  2.1 座標系  3}

σ_{k} = \{\frac{1}{1 + κ} (\frac{σ_{k−1/2}^{κ+1} − σ_{k+1/2}^{κ+1}}{σ_{k−1/2} − σ_{k+1/2}} )\}^{1/κ}  (2.3)

ただし, κ = \frac{R}{C_{p}} である. ここで, R は乾燥空気の気体定数, C_{p} は乾燥空気の等圧比熱である^{5} また, レベル加重Δσは以下のように定義される.

Δσ_{k} ≡ σ_{k−1/2} − σ_{k+1/2} (1 ≤ k ≤ K)  (2.4)
Δσ_{1/2} ≡ σ_{1/2} − σ_{1} = 1 − σ_{1}  (2.5)
Δσ_{K+1/2} ≡ σ_{K} − σ_{K+1/2} = σ_{K}  (2.6)

[figure: ・; ・; ・; \frac{1}{2}; 1; \frac{3}{2}; 2; \frac{5}{2}; K − \frac{1}{2}; K; K + \frac{1}{2}; σ = 1; σ = 0]

5いずれも定数としている.

zahyou/zahyou.tex  2008/06/10(地球流体電脳倶楽部)

\underline{4 DCPAM4 第 2 部離散化第 2 章座標系・変換公式}

2.2 水平スペクトル

ここでは, 力学過程の時間積分での計算において用いるスペクトルを導入し, 格子点での値とスペクトルの係数とのやり取りの公式を示す.

2.2.1 水平スペクトルの基底の導入

格子点上の点で定義された物理量は, 格子点上でのみ値を持つ（以下このことを, 「離散化した」と呼ぶ）球面調和函数の和の形で表現される. また, 各格子点における物理量の水平微分を評価するために, (λ, ϕ) 面で定義された（以下, 「連続系の」と呼ぶ）球面調和函数系で内挿して得られる関数を用いる. ここではその球面調和函数を導入する. なお, 簡単のために, 連続系の球面調和函数のみを陽に記す. 離散系の球面調和函数は連続系の球面調和函数に格子点の座標を代入したものから構成される.

(λ, ϕ) 面において, 球面調和函数 Y_{n}^{m}(λ, ϕ) は次のように定義される.

Y_{n}^{m}(λ, ϕ) ≡ P_{n}^{m}(sin ϕ) exp(imλ)  (2.7)

ただし, m, n は 0 ≤ |m| ≤ n を満たす整数であり, P_{n}^{m}(sin ϕ) は 2 で規格化された Legendre 函数・陪函数

P_{n}^{m}(μ) ≡ \sqrt{\frac{(2n + 1)(n − |m|)!}{(n + |m|)!}} \frac{(1 − μ^{2})^{\frac{|m|}{2}}}{2^{n}n!} \frac{d^{n+|m|}}{dμ^{n+|m|}} (μ^{2} − 1)^{n}  (2.8)
∫_{−1}^{1} P_{n}^{m}(μ)P_{n′}^{m}(μ)dμ = 2δ_{nn′}  (2.9)

である. なお, P_{n}^{0} を P_{n} とも書く.

2.2.2 波数切断

波数切断は三角形切断（T）または平行四辺形切断（R）とする. M , N は三角形切断, 平行四辺形切断のときについてそれぞれ以下のとおりである. ただし, 切断波数を N_{tr} とする.

• 三角形切断の場合

M = N_{tr}, N = N_{tr} , I ≥ 3N_{tr} + 1, かつ J ≥ \frac{3N_{tr} + 1}{2}.

自由度は, (N_{tr} + 1)^{2} である.

• 平行四辺形切断の場合 M = N_{tr}, N(m) = N_{tr} + |m| , I ≥ 3N_{tr} + 1, かつ J ≥ 3N_{tr} + 1. 自由度は, (2N_{tr} + 1)(N_{tr} + 1) である.

よく用いられる値の例としては, T42 の場合 I = 128, J = 64 , R21 の場合 I = 64, J = 64 がある.

なお, 球面調和函数についてはリファレンス「球面調和函数」を, 波数切断についてはリファレンス「波数切断」を参照せよ.

91/10/31 保坂征宏 2008/06/10(地球流体電脳倶楽部)  zahyou/spctl.tex

\underline{DCPAM4 第 2 部離散化  2.2 水平スペクトル  5}

2.2.3 離散化したスペクトルの基底の直交性

離散化した Legendre 函数と三角関数は次の直交条件を満たす^{6}.

∑_{j=1}^{J} P_{n}^{m}(μ_{j})P_{n′}^{m}(μ_{j})w_{j} = δ_{nn′}  (2.10)

∑_{i=1}^{I} exp(imλ_{i}) exp(−im′λ_{i}) = Iδ_{mm′}  (2.11)

(2.12)

ここで w_{j} は Gauss 荷重で, w_{j} ≡ \frac{(2J − 1)(1 − sin^{2} ϕ_{j})}{(JP_{J−1}(sin ϕ_{j}))^{2}} である.

2.2.4 格子点値とスペクトルの係数との変換法

物理量 A の格子点 (λ_{i}, ϕ_{j}) （ただし i = 1, 2, ⋯, I, j = 1, 2, ⋯, J）での値 A_{ij} = A(λ_{i}, ϕ_{j}) とスペクトル空間での Y_{n}^{m} （ただし m = −M, ⋯, M, n = |m|, ⋯, N(m) ）の係数 Ã_{n}^{m} とは次の変換則に従う^{7}.

A_{ij} ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} Ã_{n}^{m}Y_{n}^{m}(λ_{i}, ϕ_{j})  (2.13)

Ã_{n}^{m} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} A_{ij}Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j}  (2.14)

A が実数であることを用いると, (Ã_{n}^{m} exp(imλ))^{∗} = Ã_{n}^{−m} exp(−imλ) なので, m については負でない整数の範囲で和をとることができる^{8}. ただし, A_{n}^{m} の定義を修正していることに注意せよ.

A_{ij} = ∑_{m=0}^{M} ∑_{n=m}^{N} ℜÃ_{n}^{m}Y_{n}^{m}(λ_{i}, ϕ_{j})  (2.15)

Ã_{n}^{m} = \{ \begin{matrix}\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J}A_{ij}Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j} & m = 0, m ≤ n ≤ N \\ \frac{2}{I} ∑_{i=1}^{I} ∑_{j=1}^{J}A_{ij}Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j} & 1 ≤ m ≤ M, m ≤ n ≤ N \end{matrix}  (2.16)

6詳しくはリファレンス「球面調和函数の離散的直交関係」を参照せよ.

7正変換, 逆変換時の係数は consistent に与えてさえいれば問題がない. 現在の GCM では異なる動作をする 2 種類の FFT が用意されている. ここでは numerical recipy 版の方法に従っている. （もう一つは中村一版である.）仕様について詳しくは第３部の FFT に関する項目（「FFT99X」）を見られたい.

8さらに, 実際の計算手続きとしては, P_{n}^{m}(sin ϕ) が, n − m が偶数（even）の時 ϕ = 0 について対称, n − m が奇数（odd）の時 ϕ = 0 について反対称であることを考慮して演算回数を減らすことができる. すなわち, A_{ij} の計算では北半球のみについて南北対称成分 A_{ij}^{even} と反対称成分 A_{ij}^{odd} についてそれぞれ計算し, 南半球については Ai, J − j = A_{ij}^{even} −A_{ij}^{odd} とすればよい. また, A_{n}^{m} の計算においては, その対称性, 反対称性に基づいて A_{i,j} + A_{i,J−j} または A_{i,j} − A_{i,J−j} の一方を j について 1 から J/2 まで加えればよい.

zahyou/spctl.tex  2008/06/10(地球流体電脳倶楽部)

2.2.5 内挿公式

(λ, ϕ) 空間で定義される物理量 A(λ, ϕ) を格子点値 A_{ij} をもとに内挿する場合には, 変換公式を用いて A_{ij} から Ã_{n}^{m} を求めた上で,

A(λ, ϕ) ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} Ã_{n}^{m}Y_{n}^{m}(λ, ϕ)  (2.17)

として得る.

2.2.6 空間微分の評価

各格子点における空間微分値の評価は, 内挿公式を用いて得た連続関数の空間微分の格子点値で評価する.

• λ 微分

(\frac{∂f}{∂λ} )_{ij} ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} imf̃_{n}^{m}Y_{n}^{m}(λ_{i}, ϕ_{j})  (2.18)
˜ (\frac{∂f}{∂λ} )_{n}^{m} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} imf_{ij}Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j}  (2.19)

• μ 微分

(\frac{∂f}{∂μ} )_{ij} ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} f̃_{n}^{m} \frac{dP_{n}^{m}}{dμ} |_{j} exp(imλ_{i})  (2.20)
˜ (\frac{∂f}{∂μ} )_{n}^{m} = −\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} f_{ij} \frac{dP_{n}^{m}}{dμ} |_{j} exp(−imλ_{i})w_{j}  (2.21)

2008/06/10(地球流体電脳倶楽部)  zahyou/spctl.tex

第3章支配方程式・力学過程

3.1 はじめに

この章では支配方程式を離散化する. 空間に関する離散化である鉛直離散化と水平離散化の方法並びに時間に関する離散化について記す.

93/03/18 保坂征宏

\underline{8 DCPAM4 第 2 部離散化第 3 章支配方程式・力学過程}

3.2 鉛直離散化

ここでは支配方程式を鉛直方向に離散化する. Arakawa and Suarez(1983) に従って, 基礎方程式を鉛直方向に差分によって離散化する^{1}. 各方程式の離散化表現は次のようになる.

3.2.1 連続の式, 鉛直速度

\frac{∂π}{∂t} = − ∑_{k=1}^{K}(D_{k} + v_{k} · ∇π)Δσ_{k}  (3.1)

σ̇_{k−1/2} = −σ_{k−1/2} \frac{∂π}{∂t} − ∑_{l=k}^{K}(D_{l} + v_{l} · ∇π)Δσ_{l}  (3.2)
σ̇_{1/2} = σ̇_{K+1/2} = 0  (3.3)

3.2.2 静水圧の式

Φ_{1} = Φ_{s} + C_{p}(σ_{1}^{−κ} − 1)T_{v,1}  (3.4)
= Φ_{s} + C_{p}α_{1}T_{v,1}

Φ_{k} − Φ_{k−1} = C_{p} [(\frac{σ_{k−1/2}}{σ_{k}} )^{κ} − 1] T_{v,k} + C_{p} [1 − (\frac{σ_{k−1/2}}{σ_{k−1}} )^{κ}] T_{v,k−1}  (3.5)
= C_{p}α_{k}T_{v,k} + C_{p}β_{k−1}T_{v,k−1}

ここで,

α_{k} = (\frac{σ_{k−1/2}}{σ_{k}} )^{κ} − 1  (3.6)
β_{k} = 1 − (\frac{σ_{k+1/2}}{σ_{k}} )^{κ} .  (3.7)

3.2.3 運動方程式

\frac{∂ζ_{k}}{∂t} = \frac{1}{a(1 − μ^{2})} \frac{∂(VA)_{k}}{∂λ} − \frac{1}{a} \frac{∂(UA)_{k}}{∂μ} + D(ζ_{k})  (3.8)

93/03/18 沼口敦・保坂征宏

1このスキームは次のような特徴をもつ. ( 2005/4/4 石渡: 初の 4 つは不正確な表現かも )

• 全領域積分した質量を保存

• 全領域積分したエネルギーを保存

• 全領域積分の角運動量を保存

• 全質量積分した温位を保存

• 静水圧の式が local にきまる. (下層の高度が上層の温度に依存しない)

• 水平方向に一定の, ある特定の温度分布について, 静水圧の式が正確になり, 気圧傾度力が 0 になる.

• 等温位大気はいつまでも等温位に留まる 2008/06/10(地球流体電脳倶楽部)  rikigaku/entyoku.tex

\underline{DCPAM4 第 2 部離散化  3.2 鉛直離散化  9}
\frac{∂D}{∂t} = \frac{1}{a(1 − μ^{2})} \frac{∂(UA)_{k}}{∂λ} + \frac{1}{a} \frac{∂(VA)_{k}}{∂μ} − ∇_{σ}^{2}(Φ_{k} + C_{p}κ̂_{k}T̄_{k}π + (KE)_{k}) + D(D_{k})  (3.9)

ここで,

(UA)_{k} = (ζ_{k} + f)V_{k} − \frac{1}{2Δσ_{k}} [σ̇_{k−1/2}(U_{k−1} − U_{k}) + σ̇_{k+1/2}(U_{k} − U_{k+1})]
−\frac{C_{p}κ̂_{k}T′_{v,k}}{a} \frac{∂π}{∂λ} + F_{λk} cos φ  (3.10)

(3.11)

(VA)_{k} = −(ζ_{k} + f)U_{k} − \frac{1}{2Δσ_{k}} [σ̇_{k−1/2}(V_{k−1} − V_{k}) + σ̇_{k+1/2}(V_{k} − V_{k+1})]
−\frac{C_{p}κ̂_{k}T′_{v,k}}{a} (1 − μ^{2})\frac{∂π}{∂μ} + F_{φk} cos φ  (3.12)

ここで,

κ̂_{k} = \frac{σ_{k−1/2}(σ_{k−1/2}^{κ} − σ_{k}^{κ}) + σ_{k+1/2}(σ_{k}^{κ} − σ_{k+1/2}^{κ})}{σ_{k}^{κ}(σ_{k−1/2} − σ_{k+1/2})}
= \frac{σ_{k−1/2}α_{k} + σ_{k+1/2}β_{k}}{Δσ_{k}}  (3.13)
T_{v,k} = T_{k} − T̄_{k}  (3.14)
(KE)_{k} = \frac{U_{k}^{2} + V_{k}^{2}}{2(1 − μ^{2})}  (3.15)

3.2.4 熱力学の式

\frac{∂T_{k}}{∂t} = −\frac{1}{a(1 − μ^{2})} \frac{∂U_{k}T′_{k}}{∂λ} − \frac{1}{a} \frac{∂V_{k}T′_{k}}{∂μ} + Ĥ_{k}

+\frac{Q_{k}}{C_{p}} + D(T_{k}) + D′(v)

(3.16)

ここで,

Ĥ_{k} ≡ T′_{k}D_{k} − \frac{1}{Δσ_{k}} [σ̇_{k−1/2}(T̂_{k−1/2} − T_{k}) + σ̇_{k+1/2}(T_{k} − T̂_{k+1/2})]

+ {α_{k} [σ_{k−1/2}v_{k} · ∇π − ∑_{l=k}^{K}(D_{l} + v_{l} · ∇π)Δσ_{l}]

+ β_{k} [σ_{k+1/2}v_{k} · ∇π − ∑_{l=k+1}^{K}(D_{l} + v_{l} · ∇π)Δσ_{l}]\} \frac{1}{Δσ_{k}} T_{v,k}

= T′_{k}D_{k} − \frac{1}{Δσ_{k}} [σ̇_{k−1/2}(T̂_{k−1/2} − T_{k}) + σ̇_{k+1/2}(T_{k} − T̂_{k+1/2})]

+κ̂_{k}v_{k} · ∇πT_{v,k}

−α_{k} ∑_{l=k}^{K}(D_{l} + v_{l} · ∇π)Δσ_{l} \frac{T_{v,k}}{Δσ_{k}}

−β_{k} ∑_{l=k+1}^{K}(D_{l} + v_{l} · ∇π)Δσ_{l} \frac{T_{v,k}}{Δσ_{k}}  (3.17)
rikigaku/entyoku.tex  2008/06/10(地球流体電脳倶楽部)

T̂_{k−1/2} = \frac{[(\frac{σ_{k−1/2}}{σ_{k}} )^{κ} − 1] σ_{k−1}^{κ}T_{k} + [1 − (\frac{σ_{k−1/2}}{σ_{k−1}} )^{κ}] σ_{k}^{κ}T_{k−1}}{σ_{k−1}^{κ} − σ_{k}^{κ}}  (3.18)
= a_{k}T_{k} + b_{k−1}T_{k−1}  (3.19)

a_{k} = α_{k} [1 − (\frac{σ_{k}}{σ_{k−1}} )^{κ}]^{−1}  (3.20)

b_{k} = β_{k} [(\frac{σ_{k}}{σ_{k+1}} )^{κ} − 1]^{−1} .  (3.21)

3.2.5 水蒸気の式

\frac{∂q_{k}}{∂t} = −\frac{1}{a(1 − μ^{2})} \frac{∂U_{k}q_{k}}{∂λ} − \frac{1}{a} \frac{∂V_{k}q_{k}}{∂μ} + R_{k} + S_{q,k} + D(q_{k})  (3.22)
R_{k} = q_{k}D_{k} − \frac{1}{2Δσ_{k}} [σ̇_{k−1/2}(q_{k−1} − q_{k}) + σ̇_{k+1/2}(q_{k} − q_{k+1})]  (3.23)

2008/06/10(地球流体電脳倶楽部)  rikigaku/entyoku.tex

\underline{DCPAM4 第 2 部離散化  3.3 水平離散化  11}

3.3 水平離散化

ここでは支配方程式を水平離散化する. 水平方向の離散化はスペクトル変換法を用いる (Bourke, 1988). 経度, 緯度に関する微分の項は直交関数展開によって評価し, 一方, 非線型項は格子点上で計算する. 各方程式のスペクトル表現は以下のようになる.

3.3.1 連続の式

\frac{∂π̃_{n}^{m}}{∂t} = − ∑_{k=1}^{K}(D̃_{n}^{m})_{k}Δσ_{k}
+\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} Z_{ij}Y_{n}^{m∗}(λ_{i}, μ_{j})w_{j},  (3.24)

ここで,

Z ≡ − ∑_{k=1}^{K} v_{k} · ∇π.  (3.25)

3.3.2 運動方程式

\frac{∂ζ̃_{n}^{m}}{∂t} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} im(VA)_{ij}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}
+ \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J}(UA)_{ij}(1 − μ_{j}^{2})\frac{∂}{∂μ}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}
+ D̃_{M,n}^{m}ζ̃_{n}^{m},  (3.26)

\frac{∂D̃_{n}^{m}}{∂t} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} im(UA)_{ij}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}
− \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J}(VA)_{ij}(1 − μ_{j}^{2})\frac{∂}{∂μ}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}
− \frac{n(n + 1)}{a^{2}} \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J}(KE)_{ij}Y_{n}^{m∗}(λ_{i}, μ_{j})w_{j}
+ \frac{n(n + 1)}{a^{2}}(Φ_{n}^{m} + C_{p}κ̂_{k}T̄_{k}π_{n}^{m}) + D̃_{M,n}^{m}D̃_{n}^{m},  (3.27)

ただし,

D̃_{M,n}^{m} = −K_{HD} [(\frac{−n(n + 1)}{a^{2}})^{N_{D}/2} − (\frac{2}{a})^{N_{D}}] .  (3.28)

93/03/18 沼口敦・保坂征宏 rikigaku/suihei.tex  2008/06/10(地球流体電脳倶楽部)

\underline{12 DCPAM4 第 2 部離散化第 3 章支配方程式・力学過程}

3.3.3 熱力学の式

\frac{∂T̃_{n}^{m}}{∂t} = −\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} imU_{ij}T′_{ij}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}

+\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} V_{ij}T′_{ij}(1 − μ_{j}^{2})\frac{∂}{∂μ}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}

+\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} (Ĥ_{ij} + \frac{Q_{ij}}{C_{p}} ) Y_{n}^{m∗}(λ_{i}, μ_{j})w_{j}

+D̃_{H,n}^{m}T̃_{n}^{m}

+\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} D′_{ij}(v)Y_{n}^{m∗}(λ_{i}, μ_{j})w_{j},  (3.29)

ただし,

D̃_{H,n}^{m} = −K_{HD} (\frac{−n(n + 1)}{a^{2}})^{N_{D}/2} .  (3.30)

3.3.4 水蒸気の式

\frac{∂q̃_{n}^{m}}{∂t} = −\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} imU_{ij}q_{ij}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}

+\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} V_{ij}q_{ij}(1 − μ_{j}^{2})\frac{∂}{∂μ}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{w_{j}}{a(1 − μ_{j}^{2})}

+\frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} (R̂_{ij} + S_{q,ij}) Y_{n}^{m∗}(λ_{i}, μ_{j})w_{j}

+D̃_{H,n}^{m}q̃_{n}^{m}  (3.31)

2008/06/10(地球流体電脳倶楽部)  rikigaku/suihei.tex

\underline{DCPAM4 第 2 部離散化  3.4 時間積分  13}

3.4 時間積分

ここでは時間積分スキームについて記す.

時間差分スキームは基本的に leap frog である. ただし, 拡散項および物理過程の項は後方差分もしくは前方差分とする. 計算モードを抑えるために時間フィルター (Asselin, 1972) を用いる. さらに Δt を大きくとるために, 重力波の項に semi-implicit の手法を適用する (Bourke, 1988).

3.4.1 leap frog による時間積分と時間フィルター

移流項等の時間積分スキームとして leap frog を用いる. 水平拡散項には 2Δt の後方差分を使用する. 物理過程の項 (F_{λ}, F_{φ}, Q, S_{q}) には 2Δt の前方差分を使用する. ただし, 鉛直拡散の時間変化項の計算に関しては後方差分的な取扱いをする. また, 対流・凝結過程に関しては一度それらを考慮せずに時間積分して, その後で場を修正するという形式をとる (adjustment). 摩擦熱の項も補正として扱う. 各予報変数の総体を A と表すと,

Â^{t+Δt} = Ā^{t−Δt} + 2ΔtȦ_{adv} (A^{t}) + 2ΔtȦ_{dif} \binom{}{Â^{t+Δt}} + 2ΔtȦ_{phy} \binom{}{Ā^{t−Δt}}  (3.32)

Ȧ_{adv} は移流項等, Ȧ_{dif} は水平拡散項, Ȧ_{phy} は物理過程等による項である. Â^{t+Δt} には, 摩擦熱 (Ȧ_{dis}) および対流・凝結過程 (Ȧ_{cnd}) の補正が加えられ, A^{t+Δt} となる.

A^{t+Δt} = Â^{t+Δt} + 2ΔtȦ_{dis} \binom{}{Â^{t+Δt}} + 2ΔtȦ_{cnd} \binom{}{Â^{t+Δt}}  (3.33)

leap frog における計算モードの除去のために Asselin(1972) の時間フィルターを毎ステップ適用する. すなわち,

Ā^{t} = (1 − 2ϵ_{f})A^{t} + ϵ_{f} (Ā^{t−Δt} + A^{t+Δt})  (3.34)

と Ā を求める. ϵ_{f} としては標準的に 0.05 を使用する.

3.4.2 semi-implicit 時間積分

方程式系において, T = T̄_{k} であるような静止場を基本場とする線型重力波項とそれ以外の項 (添字 NG を付ける) に分離する. 鉛直方向のベクトル表現 D = {D_{k}}, T = {T_{k}} を用いて,

\frac{∂π}{∂t} = (\frac{∂π}{∂t} )_{NG} − C · D,  (3.35)
\frac{∂D}{∂t} = (\frac{∂D}{∂t} )_{NG} − ∇_{σ}^{2}(Φ_{S} + \underline{W} T + Gπ) + D_{M}D,  (3.36)

\frac{∂T}{∂t} = (\frac{∂T}{∂t} )_{NG} − \underline{h}D + D_{H}T ,  (3.37)

ここで, 非重力波項は,

(\frac{∂π}{∂t} )^{NG} = − ∑_{k=1}^{K} v_{k} · ∇πΔσ_{k} = Z_{k}  (3.38)

93/03/18 沼口敦・保坂征宏 rikigaku/jikan.tex  2008/06/10(地球流体電脳倶楽部)

σ̇_{k−1/2}^{NG} = −σ_{k−1/2} (\frac{∂π}{∂t} )^{NG} − ∑_{l=k}^{K} v_{l} · ∇πΔσ_{l}  (3.39)

(\frac{∂D}{∂t} )^{NG} = \frac{1}{a(1 − μ^{2})} \frac{∂(UA)_{k}}{∂λ} + \frac{1}{a} \frac{∂(VA)_{k}}{∂μ} −∇_{σ}^{2}(KE)_{k} −∇_{σ}^{2} ∑_{k=1}^{K} W_{kl}(T_{v,l} −T_{l})+D(D_{k}) (3.40)
(\frac{∂T_{k}}{∂t} )^{NG} ≡ −\frac{1}{a(1 − μ^{2})} \frac{∂U_{k}T′_{k}}{∂λ} − \frac{1}{a} \frac{∂V_{k}T′_{k}}{∂μ} + H_{k} + D(T_{k}) + D′(v)  (3.41)
H_{k} = T′_{k}D_{k}

−\frac{1}{Δσ_{k}} [σ̇_{k−1/2}(T̂′_{k−1/2} − T′_{k}) + σ̇_{k+1/2}(T′_{k} − T̂′_{k+1/2})]

−\frac{1}{Δσ_{k}} [σ̇_{k−1/2}^{NG}(ˆT̄_{k−1/2} − T̄_{k}) + σ̇_{k+1/2}^{NG}(T̄_{k} − ˆ T̄_{k+1/2})]

+κ̂_{k}T_{v,k}v_{k} · ∇π

−\frac{α_{k}}{Δσ_{k}} T_{v,k} ∑_{l=k}^{K} v_{l} · ∇πΔσ_{l} − \frac{β_{k}}{Δσ_{k}} T_{v,k} ∑_{l=k+1}^{K} v_{l} · ∇πΔσ_{l}

−\frac{α_{k}}{Δσ_{k}} T′_{v,k} ∑_{l=k}^{K} D_{l}Δσ_{l} − \frac{β_{k}}{Δσ_{k}} T′_{v,k} ∑_{l=k+1}^{K} D_{l}Δσ_{l}  (3.42)

また, 重力波項のベクトルおよび行列 (下線で表示) は,

C_{k} = Δσ_{k}  (3.43)
W_{kl} = C_{p}α_{l}δ_{k≥l} + C_{p}β_{l}δ_{k−1≥l}  (3.44)
G_{k} = κ̂_{k}C_{p}T̄_{k}  (3.45)
\underline{h} = \underline{Q}\underline{S} − \underline{R}  (3.46)
Q_{kl} = \frac{1}{Δσ_{k}} (ˆT̄_{k−1/2} − T̄_{k})δ_{k=l} + \frac{1}{Δσ_{k}} (T̄_{k} − ˆ T̄_{k+1/2})δ_{k+1=l}  (3.47)
S_{kl} = σ_{k−1/2}Δσ_{l} − Δσ_{l}δ_{k≤l}  (3.48)
R_{kl} = − (\frac{α_{k}}{Δσ_{k}} Δσ_{l}δ_{k≤l} + \frac{β_{k}}{Δσ_{k}} Δσ_{l}δ_{k+1≤l}) T̄_{k}.  (3.49)

ここで, 例えば δ_{k≤l} は, k ≤ l が成り立つとき 1, そうでないとき 0 となる関数である. 次のような表現を使用して,

δ_{t}A ≡ \frac{1}{2Δt} (A^{t+Δt} − A^{t−Δt})  (3.50)
\overline{A}^{t} ≡ \frac{1}{2} (A^{t+Δt} + A^{t−Δt})
= A^{t−Δt} + δ_{t}AΔt,  (3.51)

方程式系に semi-implicit 法を適用すると,

δ_{t}π = (\frac{∂π}{∂t} )_{NG} − C · \overline{D}^{t}  (3.52)
δ_{t}D = (\frac{∂D}{∂t} )_{NG} − ∇_{σ}^{2}(Φ_{S} + \underline{W} \overline{T}^{t} + G\overline{π}^{t}) + D_{M}(D^{t−Δt} + 2Δtδ_{t}D)  (3.53)
2008/06/10(地球流体電脳倶楽部)  rikigaku/jikan.tex

δ_{t}T = (\frac{∂T}{∂t} )_{NG} − \underline{h}\overline{D}^{t} + D_{H}(T^{t−Δt} + 2Δtδ_{t}T)  (3.54)

すると,

\{(1 − 2ΔtD_{H})(1 − 2ΔtD_{M})\underline{I} − (Δt)^{2}(\underline{W} \underline{h} + GC^{T})∇_{σ}^{2}\} \overline{D}^{t}

= (1 − 2ΔtD_{H})(1 − ΔtD_{M})D^{t−Δt} + (1 − 2ΔtD_{H})Δt (\frac{∂D}{∂t} )_{NG}

−Δt∇_{σ}^{2} \{(1 − 2ΔtD_{H})Φ_{S} + \underline{W} [(1 − 2ΔtD_{H})T^{t−Δt} + Δt (\frac{∂T}{∂t} )_{NG}]

+(1 − 2ΔtD_{H})G [π^{t−Δt} + Δt (\frac{∂π}{∂t} )_{NG}]\} .  (3.55)

球面調和関数展開を用いているので,

∇_{σ}^{2} = −\frac{n(n + 1)}{a^{2}}

であり上式を\overline{D_{n}^{m}}^{t} について解くことができる. その後,

D^{t+Δt} = 2\overline{D}^{t} − D^{t−Δt}  (3.56)

および, (3.52), (3.54) により t + Δt における値 Â^{t+Δt} が求められる.

rikigaku/jikan.tex  2008/06/10(地球流体電脳倶楽部)

第4章参考文献

4.1 座標系・変換公式

気象庁予報部, 1982 : スペクトル法による数値予報（その原理と実際） . 気象庁, 111pp. 寺沢寛一, 1983 : 自然科学者のための数学概論（増訂版） . 岩波書店, 711pp. Haltiner,G.J., Williams,R.T., 1980 : Numerical Prediction and Dynamic Meteorology. John

Wiley & Sons, 477pp. Arakawa A., Suarez M.J., 1983: Vertical differencing of the primitive equations in sigma cood-

inates. Mon. Weather Rev., 111, 34–45. 森口, 宇田川, 一松編 ,1956 : 岩波数学公式 I . 岩波書店, 318pp. 森口, 宇田川, 一松編 ,1960 : 岩波数学公式 III . 岩波書店, 310pp. 一松信, 1982 : 数値解析. 朝倉書店, 163pp. 森正武, 1984 : 数値解析法. 朝倉書店, 202pp.

4.2 力学過程

4.3 雲

4.4 放射

4.5 鉛直輸送

4.6 地表面過程

91/10/31 保坂征宏

第A章座標系・変換公式

A.1 球面調和函数

ここでは連続系での球面調和函数を定義し, スペクトル計算の理解に必要な性質を挙げ, 証明する.

まず球面調和函数を定義し, 次いで球面調和函数が完全直交系をなすことを主張する. このことにより, 球面上に分布するあらゆる連続関数が球面調和函数の重ね合わせで一意的に表されることになる.

球面調和函数は 2 次元ラプラシアンに関する固有関数であり, このために全波数という概念が生まれる. 参考までにこのことも記しておく.

さらに, 球面調和函数を空間微分した結果も書いておく.

1. 定義と性質（球面調和函数, Legendre 函数, Legendre 陪函数）

2. 空間微分

3. 全波数の概念

また, イメージをつかむために, ルジャンドル (陪) 関数のグラフを示す.

A.1.1 定義と性質

ここでは, 岩波公式集^{1}の Legendre 函数・陪函数 P̃_{n}^{m} , 2 で規格化した Legendre 函数・陪函数 P_{n}^{m}, 4π で規格化した球面調和函数 Y_{n}^{m} の順に定義する. さらにそれらの性質として, 従う微分方程式, 漸下式, 完全規格直交性について述べる.

岩波公式集の Legendre 函数・陪函数 P̃_{n}^{m}

• 定義岩波公式集によると Legendre 函数・陪函数 P̃_{n}^{m}(μ) は −1 ≤ μ ≤ 1 において次式で定義される（Rodrigues の公式） .

P̃_{n}^{m} ≡ \frac{(1 − μ^{2})^{\frac{|m|}{2}}}{2^{n}n!} \frac{d^{n+|m|}}{dμ^{n+|m|}} (μ^{2} − 1)^{n}  (A.1)

ただし, m, n は 0 ≤ |m| ≤ n を満たす整数である. Legendre 函数 P̃_{n}^{0} を P̃_{n} とも書く.

91/12/09 保坂征宏

2005/04/04 石渡正樹

1森口, 宇田川, 一松編「数学公式 III」 ,1960 を指す.

\underline{2 DCPAM4 第 2 部離散化付録 A 座標系・変換公式}

• Legendre 函数・陪函数の満たす方程式 P̃_{n}^{m}(μ) は次の方程式を満たす.

\frac{d}{dμ} \{(1 − μ^{2})\frac{d}{dμ} P̃_{n}^{m}\} + \{n(n + 1) − \frac{m^{2}}{1 − μ^{2}} \} P̃_{n}^{m} = 0  (A.2)

ただし, m, n は 0 ≤ |m| ≤ n を満たす整数である.

• Legendre 函数・陪函数の従う漸化式

P̃_{n}^{m}(μ) は次の漸化式に従う.

(n − |m| + 1)P̃_{n+1}^{m} − (2n + 1)μP̃_{n}^{m} + (n + |m|)P̃_{n−1}^{m} = 0  (A.3)

ただし, m, n は 0 ≤ |m| ≤ n − 1, または m = n = 0 を満たす整数である.

さらに, 次の関係式が成り立つ.

(1 − μ^{2})\frac{d}{dμ} P̃_{n}^{m} = (n + |m|)P̃_{n−1}^{m} − nμP̃_{n}^{m}  (A.4)

ただし, m, n は 0 ≤ |m| ≤ n − 1 を満たす整数である.

• 完全規格直交性 P̃_{n}^{m}(μ) (n = |m|, |m + 1, ⋯) は次の直交関係を満たす.

∫_{−1}^{1} P̃_{n}^{m}(μ)P̃_{n′}^{m}(μ)dμ = \frac{2}{2n + 1} \frac{(n + |m|)!}{(n − |m|)!}δ_{nn′}  (A.5)

ただし, m, n, n′ は 0 ≤ |m| ≤ n, n′ を満たす整数である.

−1 ≤ μ ≤ 1 で定義される連続関数 A(μ) は {P̃_{n}^{m}|n = |m|, |m + 1|, ⋯} を用いて

A(μ) = ∑_{n=|m|}^{∞} Ã_{n}^{m}P̃_{n}^{m}(μ),  (A.6)
Ã_{n}^{m} = \frac{2n + 1}{2} \frac{(n − |m|)!}{(n + |m|)!} ∫_{−1}^{1} A(μ)P̃_{n}^{m}(μ)dμ  (A.7)

と表される.

2 で規格化した Legendre 函数・陪函数 P_{n}^{m}

• 定義 2 で規格化した Legendre 函数・陪函数 P_{n}^{m}(μ) は −1 ≤ μ ≤ 1 において次式で定義される.

P_{n}^{m} ≡ \sqrt{\frac{(2n + 1)(n − |m|)!}{(n + |m|)!}}P̃_{n}^{m} = \sqrt{\frac{(2n + 1)(n − |m|)!}{(n + |m|)!}} \frac{(1 − μ^{2})^{\frac{|m|}{2}}}{2^{n}n!} \frac{d^{n+|m|}}{dμ^{n+|m|}} (μ^{2} − 1)^{n}, (A.8)

(A.9)

ただし, m, n は 0 ≤ |m| ≤ n を満たす整数である. Legendre 函数 P_{n}^{0} を P_{n} とも書く.

2008/06/10(地球流体電脳倶楽部)  zahyou/kyuumen.tex

\underline{DCPAM4 第 2 部離散化  A.1 球面調和函数  3}

• Legendre 函数・陪函数の満たす方程式 P_{n}^{m}(μ) は, 次の方程式を満たす.

\frac{d}{dμ} \{(1 − μ^{2})\frac{d}{dμ} P_{n}^{m}\} + \{n(n + 1) − \frac{m^{2}}{1 − μ^{2}} \} P_{n}^{m} = 0  (A.10)

ただし, m, n は 0 ≤ |m| ≤ n を満たす整数である.

• Legendre 函数・陪函数の従う漸化式 P_{n}^{m}(μ) は, 次の漸化式に従う.

(n − |m| + 1)\sqrt{\frac{1}{2n + 3} \frac{(n + 1 + |m|)!}{(n + 1 − |m|)!}}P_{n+1}^{m} − (2n + 1)\sqrt{\frac{1}{2n + 1} \frac{(n + |m|)!}{(n − |m|)!}}μP_{n}^{m} (A.11)

+(n + |m|)\sqrt{\frac{1}{2n − 1} \frac{(n − 1 + |m|)!}{(n − 1 − |m|)!}}P_{n−1}^{m} = 0  (A.12)
∴ P_{n+1}^{m} = \sqrt{\frac{(2n + 1)(2n + 3)}{(n − |m| + 1)(n + |m| + 1)}}μP_{n}^{m}  (A.13)
−\sqrt{\frac{(2n + 1)(2n + 3)}{(n − |m| + 1)(n + |m| + 1)}}\sqrt{\frac{(n − |m|)(n + |m|)}{(2n + 1)(2n − 1)}}P_{n−1}^{m} (A.14)

ただし, m, n は 0 ≤ |m| ≤ n − 1, または m = n = 0 を満たす整数である.

さらに次の関係式が成り立つ.

(1 − μ^{2})\frac{d}{dμ} P_{n}^{m} = (n + |m|)\sqrt{\frac{(n − |m|)(2n + 1)}{(n + |m|)(2n − 1)}}P_{n−1}^{m} − nμP_{n}^{m}  (A.15)

ただし, m, n は 0 ≤ |m| ≤ n − 1 を満たす整数である.

• 完全規格直交性 P_{n}^{m}(μ) (n = |m|, |m + 1, ⋯) は次の直交関係を満たす.

∫_{−1}^{1} P_{n}^{m}(μ)P_{n′}^{m}(μ)dμ = 2δ_{nn′}  (A.16)

ただし, m, n, n′ は 0 ≤ |m| ≤ n, n′ を満たす整数である.

−1 ≤ μ ≤ 1 で定義される連続関数 A(μ) は {P_{n}^{m}|n = |m|, |m + 1|, ⋯} を用いて

A(μ) = ∑_{n=|m|}^{∞} Ã_{n}^{m}P_{n}^{m}(μ),  (A.17)
Ã_{n}^{m} = \frac{1}{2} ∫_{−1}^{1} A(μ)P_{n}^{m}(μ)dμ  (A.18)

と表される.

球面調和函数 Y_{n}^{m}

• 定義球面調和函数 Y_{n}^{m}(λ, ϕ) は Legendre 函数 P_{n}^{m}(sin ϕ), 三角関数^{2}exp(imλ) を用いて次のように定義される.

Y_{n}^{m}(λ, ϕ) ≡ P_{n}^{m}(sin ϕ) exp(imλ)  (A.19)

ただし, m, n は 0 ≤ |m| ≤ n を満たす整数である.

• 球面調和函数の満たす方程式

Y_{n}^{m}(λ, ϕ) は次の方程式を満たす.

[\frac{1}{cos ϕ} \frac{∂}{∂ϕ} (cos ϕ\frac{∂}{∂ϕ} ) + \frac{1}{cos^{2} ϕ} \frac{∂^{2}}{∂λ^{2}} + n(n + 1)] Y_{n}^{m} = 0  (A.20)

すなわち,

[\frac{∂}{∂μ} ((1 − μ^{2})\frac{∂}{∂μ} ) + \frac{1}{1 − μ^{2}} \frac{∂^{2}}{∂λ^{2}} + n(n + 1)] Y_{n}^{m} = 0  (A.21)

(A.22)

の解である. ただし, m, n は 0 ≤ |m| ≤ n を満たす整数である.

• 完全規格直交性

Y_{n}^{m} は次の直交関係を満たす.

∫_{−1}^{1} Y_{n}^{m}(λ, ϕ)Y_{n′}^{m′∗}(λ, ϕ)d(sin ϕ)dλ = 4πδ_{mm′}δ_{nn′}  (A.23)

ただし, m, m′, n, n′ は 0 ≤ |m| ≤ n と 0 ≤ |m′| ≤ n′ とを満たす整数である.

球面上で定義される連続関数 A(λ, ϕ) は {Y_{n}^{m}|m = 0, 1, 2, ⋯, n = |m|, |m + 1|, ⋯} を用いて

A(λ, ϕ) = ∑_{m=0}^{∞} ∑_{n=|m|}^{∞} Ã_{n}^{m}Y_{n}^{m}(λ, ϕ),  (A.24)
Ã_{n}^{m} = \frac{1}{4π} ∫_{−1}^{1} d(sin ϕ) ∫_{0}^{2π} dλA(λ, ϕ)Y_{n}^{m∗}(λ, ϕ)  (A.25)

と表される.

A.1.2 球面調和函数の空間微分

ここでは, 球面調和函数 Y_{n}^{m}(ϕ, λ) の • x 微分

• y 微分

• 2 次元ラプラシアンの計算をする.

2exp(imλ) は ∫_{0}^{2π} exp(imλ) exp(−im′λ)dλ = 2πδ_{mm′} を満たす. ただし, m, m′ は整数である.

2008/06/10(地球流体電脳倶楽部)  zahyou/kyuumen.tex

\underline{DCPAM4 第 2 部離散化  A.1 球面調和函数  5}

x 微分

\frac{1}{r cos ϕ} \frac{∂Y_{n}^{m}}{∂λ} = \frac{1}{r cos ϕ} \frac{∂}{∂λ} (P_{n}^{m}(sin ϕ) exp(imλ)) = \frac{im}{r cos ϕ}P_{n}^{m}(sin ϕ) exp(imλ)  (A.26)

y 微分

\frac{1}{r} \frac{∂Y_{n}^{m}}{∂ϕ} = \frac{1}{r} \frac{∂}{∂ϕ} (P_{n}^{m}(sin ϕ) exp(imλ)) = \frac{\sqrt{1 − μ^{2}}}{r} \frac{d}{dμ} P_{n}^{m}(μ) exp(imλ)  (A.27)

2 次元ラプラシアン

∇_{H}^{2}Y_{n}^{m} ≡ \frac{1}{r^{2}} [\frac{∂}{∂μ} ((1 − μ^{2})\frac{∂}{∂μ} ) + \frac{1}{1 − μ^{2}} \frac{∂^{2}}{∂λ^{2}} ] Y_{n}^{m}  (A.28)

= \frac{1}{r^{2}} [\frac{1}{cos ϕ} \frac{∂}{∂ϕ} (cos ϕ\frac{∂}{∂ϕ} ) + \frac{1}{cos^{2} ϕ} \frac{∂^{2}}{∂λ^{2}} ] Y_{n}^{m}  (A.29)

= −\frac{n(n + 1)}{r^{2}}Y_{n}^{m}  (A.30)

A.1.3 コメント — 全波数について

球面調和函数 Y_{n}^{m}(λ, ϕ) において n のことを全波数と呼ぶ.

全波数には, 座標系の回転に関して不変である, という特徴がある. すなわち, 任意の Y_{n}^{m}(λ, ϕ) は回転して得られる座標系 (λ′, ϕ′) における全波数 n の球面調和函数 {Y_{n}^{m}(λ′, ϕ′)|m = −n, −n+1, ⋯, n} の和で表現できる：

Y_{n}^{m}(λ, ϕ) = ∑_{m′=−n}^{n} A_{n}^{m′}Y_{n}^{m′∗}(λ′, ϕ′)  (A.31)

のである^{3}. この特徴は, 球面調和函数が 2 次元ラプラシアンの固有値であることによっている^{4}.

3この特徴を言い替えれば, 全波数 n の球面調和函数の重ね合わせで表現できる分布関数は座標系を回転させた系においても全波数 n の球面調和函数の重ね合わせで表現できることになる.

4∇_{H}^{2} ≡ \frac{1}{r^{2}} [\frac{∂}{∂ϕ} (cos ϕ \frac{∂}{∂ϕ} ) + \frac{1}{cos^{2} ϕ} \frac{∂^{2}}{∂λ^{2}} ] の, 固有値を − \frac{n(n+1)}{r^{2}} とする固有関数であることと, スカラー演算子 ∇_{H}^{2}

が座標系の回転に関して不変な演算子であることとに起因する.

すなわち, ∇_{H}^{2}Y_{n}^{m}(λ, ϕ) = − \frac{n(n+1)}{r^{2}} Y_{n}^{m}(λ, ϕ) より, 球面調和函数 Y_{n}^{m} exp(imλ) は固有値を − \frac{n(n+1)}{r^{2}} とする ∇_{H}^{2} の固有関数である. {Y_{n}^{m}|n = 0, 1, 2, ⋯, m = −n, −n + 1, ⋯, n} の完全直交性より, {Y_{n}^{m}|m = −n, −n + 1, ⋯, n} は ∇_{H}^{2}f = − \frac{n(n+1)}{r^{2}} f の解空間を張っている基底である.

座標系を回転させて, 新たな座標系での球面調和函数 Y_{n}^{m}(λ′, ϕ′) の和の形で前の座標系での球面調和函数 Y_{n}^{m}(λ, ϕ) を表現することを考えよう.

絶対系で見て同じ位置の値を比べると, 2 次元ラプラシアンを演算した値は不変なので, 前の座標系での球面調和函数 Y_{n}^{m}(λ′, ϕ′) は新たな座標系においても ∇_{H}^{′2}Y_{n}^{m} = − \frac{n(n+1)}{r^{2}} Y_{n}^{m} の解である. 新たな座標系の球面調和函数の集合

\{Y_{n}^{m}(λ′, ϕ′)|m = −n, −n + 1, ⋯, n\} も ∇_{H}^{′2}Y_{n}^{m} = − \frac{n(n+1)}{r^{2}} Y_{n}^{m} の解空間の基底である. したがって, 前の座標系の球面調和函数は新たな座標系の球面調和函数の和の形で書ける.

zahyou/kyuumen.tex  2008/06/10(地球流体電脳倶楽部)

\underline{6 DCPAM4 第 2 部離散化付録 A 座標系・変換公式}

A.1.4 グラフ

P_{n}^{m}(μ) の概形をつかむために, 2 で規格化した P_{n}, P_{n}^{1}, P_{n}^{2} ^{5} のグラフを示す.

岩波公式集の Legendre 函数 P̃_{n} のグラフ（森口, 宇田川, 一松, 1960）

Legendre 函数\overline{P_{n}^{1}} = P_{n}^{1}/\sqrt{2}, \overline{P_{n}^{2}} = P_{n}^{2}/\sqrt{2} のグラフ（森口, 宇田川, 一松, 1960）

5(2005/4/4 石渡) 関数形も書いておきたい. グラフは自分で描きたい.

2008/06/10(地球流体電脳倶楽部)  zahyou/kyuumen.tex

\underline{DCPAM4 第 2 部離散化  A.2 微分公式, ＧＣＭの変数の微分関係式  7}

A.2 微分公式, ＧＣＭの変数の微分関係式

6

ここでは, スカラー量, ベクトルの微分を計算する. さらにそれらを元に, 発散 D, 渦度 ζ, 速度ポテンシャル χ, 流線関数 ψ と (U, V) との関係を付ける.

A.2.1 スカラー量の微分

スカラー量 f(λ, ϕ) の x 微分は\frac{1}{r cos ϕ} \frac{∂f}{∂λ} で与えられる.

f の y 微分は\frac{1}{r} \frac{∂f}{∂ϕ} (= \frac{cos ϕ}{r} \frac{∂f}{∂μ} ) で与えられる.

f の 2 次元ラプラシアンは

∇_{H}^{2}f ≡ \frac{1}{r^{2}} [\frac{1}{cos ϕ} \frac{∂}{∂ϕ} (cos ϕ\frac{∂}{∂ϕ} ) + \frac{1}{cos^{2} ϕ} \frac{∂^{2}}{∂λ^{2}} ] f  (A.32)
= \frac{1}{r^{2}} [\frac{∂}{∂μ} \{(1 − μ^{2})\frac{∂}{∂μ} \} + \frac{1}{1 − μ^{2}} \frac{∂^{2}}{∂λ^{2}} ] f  (A.33)

で与えられる.

A.2.2 ベクトル量の微分

2 次元ベクトル場 v = (v_{1}, v_{2}) の水平発散は

div_{H}v ≡ \frac{1}{r cos ϕ} \frac{∂v_{1}}{∂λ} + \frac{1}{r cos ϕ} \frac{∂}{∂ϕ} (v_{2} cos ϕ)  (A.34)
= \frac{1}{r\sqrt{1 − μ^{2}}} \frac{∂v_{1}}{∂λ} + \frac{1}{r} \frac{∂}{∂μ} (\sqrt{1 − μ^{2}}v_{2})  (A.35)

で与えられる. v の回転の r 成分は,

(rotv)_{r} ≡ \frac{1}{r cos ϕ} \frac{∂v_{2}}{∂λ} − \frac{1}{r cos ϕ} \frac{∂}{∂ϕ} (v_{1} cos ϕ)  (A.36)
= \frac{1}{r\sqrt{1 − μ^{2}}} \frac{∂v_{2}}{∂λ} − \frac{1}{r} \frac{∂}{∂μ} (\sqrt{1 − μ^{2}}v_{1})  (A.37)

で与えられる.

以上で得られた微分公式を元に, 以下に実際に GCM で使用する便利な微分の公式を並べておく.

91/12/04 保坂征宏

6(2004/5/27 石渡) ここにあるのは変かも.

zahyou/spbibun.tex  2008/06/10(地球流体電脳倶楽部)

A.2.3 発散

水平分布する速度場

(u, v) ≡ (\frac{U}{cos ϕ}, \frac{V}{cos ϕ})  (A.38)

の水平発散 D を, U, V を用いて表す.

D = \frac{1}{r cos ϕ} \frac{∂u}{∂λ} + \frac{1}{r cos ϕ} \frac{∂}{∂ϕ} (v cos ϕ)  (A.39)
= \frac{1}{r cos^{2} ϕ} \frac{∂U}{∂λ} + \frac{1}{r cos ϕ} \frac{∂V}{∂ϕ}  (A.40)
= \frac{1}{r(1 − μ^{2})} \frac{∂U}{∂λ} + \frac{1}{r} \frac{∂V}{∂μ}  (A.41)

A.2.4 渦度

水平分布する速度場

(u, v) = (\frac{U}{cos ϕ}, \frac{V}{cos ϕ})  (A.42)

の渦度 ζ を, U, V を用いて表す.

ζ = \frac{1}{r cos ϕ} \frac{∂v}{∂λ} − \frac{1}{r cos ϕ} \frac{∂}{∂ϕ} (u cos ϕ)  (A.43)
= \frac{1}{r cos^{2} ϕ} \frac{∂V}{∂λ} − \frac{1}{r cos ϕ} \frac{∂U}{∂ϕ}  (A.44)
= \frac{1}{r(1 − μ^{2})} \frac{∂V}{∂λ} − \frac{1}{r} \frac{∂U}{∂μ}  (A.45)

A.2.5 速度ポテンシャル, 流線関数と (U, V)

速度ポテンシャル χ, 流線関数 ψ は

D ≡ ∇_{H}^{2}χ  (A.46)
ζ ≡ ∇_{H}^{2}ψ  (A.47)

で定義される. (U, V) を χ, ψ で表す.

U = −\frac{1 − μ^{2}}{r} \frac{∂ψ}{∂μ} + \frac{1}{r} \frac{∂χ}{∂λ}  (A.48)
V = \frac{1}{r} \frac{∂ψ}{∂λ} + \frac{1 − μ^{2}}{r} \frac{∂χ}{∂μ}  (A.49)

となる.

2008/06/10(地球流体電脳倶楽部)  zahyou/spbibun.tex

DCPAM4 第 2 部離散化  A.3 Legendre 函数 P_{n} の性質  9

A.3 Legendre 函数 P_{n} の性質

ここでは Legendre 函数 P_{n} の性質である

1. n − 1 次以下の多項式との積を −1 ≤ μ ≤ 1 まで積分すると零になること

2. P_{n}(μ) が −1 < μ < 1 に n 個の零点を持つこと,

を記す. 1 より Gauss 格子を定義することが保証される. また, 1, 2 は共に Gauss-Legendre の公式の証明に用いられる.

A.3.1 多項式と Legendre 函数の積の積分

P_{n}(μ) は, μ の n 次多項式である. n − 1 次以下の任意の多項式は P_{0} ∼ P_{n−1} の和で表されること, P_{n} の直交性から明らかに, n − 1 次以下の任意の多項式 f(μ) との積を積分すると

∫_{−1}^{1} f(μ)P_{n}(μ)dμ = 0  (A.50)

が成り立つことがわかる.

A.3.2 Legendre 函数の零点

P_{n} は −1 < μ < 1 に n 個の互いに異なる零点を持っている. このことについて, 以下に証明しておく. （寺沢, 1983 の 10.7 節より）

1. f(x) = (x − 1)^{n}(x + 1)^{n} を導入する.

2. f = 0 の解は, x = −1, 1 である. ゆえに, Rolle の定理により, f′ はある α （−1 < α < 1 ）で f′(α) = 0 となる. f′ = 2nx(x^{2} − 1)^{n−1} より, f′ = 0 の解は x = −1, α, 1 のみである.

3. 同様に, f′′ = 0 の解は x = −1, β_{1}, β_{2}, 1 （−1 < β_{1} < β_{2} < 1）のみ.

4. 以上を繰り返すと, f^{(n)} = 0 の解は −1 と 1 の間で互いに異なる n 個の解を持つ. （x = −1, 1 は解でないことに注意せよ. ）

5. したがって, P_{n} = \frac{1}{2^{n}n!} \frac{d^{n}}{dμ^{n}} (μ^{2} − 1)^{n} は −1 と 1 の間で互いに異なる n 個の解を持つ. （証明終り）

この零点の求め方としては, x_{j} = cos \frac{j − 1/2}{n}π を近似解として Newton 法を用いるという方法がある.

91/12/04 保坂征宏

2005/04/04 石渡正樹 zahyou/leg.tex  2008/06/10(地球流体電脳倶楽部)

\underline{10 DCPAM4 第 2 部離散化付録 A 座標系・変換公式}

A.4 積分評価

A.4.1 Gauss の台形公式

ここでは Gauss の台形公式を示す. 波数 M 以下の三角函数で表現される g(λ) （0 ≤ λ < 2π）

g(λ) = ∑_{m=−M}^{m=M} g_{m} exp(imλ)  (A.51)

について M < I を満たすように I をとると,

\frac{1}{2π} ∫_{0}^{2π} g(λ)dλ = \frac{1}{I} ∑_{n=1}^{I} g(λ_{n})  (A.52)
λ_{n} = \frac{2π(n − 1)}{I} (n = 1, 2, ⋯, I)  (A.53)

が成り立つ. これを Gauss の台形公式という. より実用的な公式は,

∑_{n=1}^{I} exp(imλ_{n}) = \{ \begin{matrix}I & (m = 0) \\ 0 & (0 < |m| < I) \end{matrix}  (A.54)
λ_{n} = \frac{2π(n − 1)}{I} (n = 1, 2, ⋯, I)  (A.55)

である. この証明は, I > M (|m| の最大値) より m ≠ 0 の時には exp(imλ_{n}) = exp (\frac{2πim(n − 1)}{I}) において, 全ての n について m(n − 1) が I の整数倍になることがないことを考慮すると明らかである ( m, n − 1 はともに I よりも小さい整数なので, m(n − 1) は I の整数倍にならない) ^{7}. 以下に Gauss の台形公式の証明を記す. まず, 左辺を計算すると,

\frac{1}{2π} ∫_{0}^{2π} g(λ)dλ = ∑_{m=−M}^{M} \frac{1}{2π} g_{m} ∫_{0}^{2π} exp(imλ)dλ = g_{0}  (A.57)

である. ここで, ∫_{0}^{2π} exp(imλ)dλ は m = 0 の項しか残らないことを使った. 一方右辺は

\frac{1}{I} ∑_{n=1}^{I} g(λ_{n}) = \frac{1}{I} ∑_{n=1}^{I} ∑_{m=−M}^{M} g_{m} exp(imλ_{n})  (A.58)
= g_{0} + ∑_{m=−M,m≠0}^{M} \frac{g_{m}}{I} ∑_{n=1}^{I} (exp(\frac{2πim}{I}))^{n−1}  (A.59)
∑_{n=1}^{I} exp \{im \frac{2π(n − 1)}{I}\} = \frac{1 − (e^{\frac{im2π}{I}} )^{I}}{1 − e^{\frac{im2π}{I}}} = \frac{1 − e^{im2π}}{1 − e^{\frac{im2π}{I}}} = 0  (A.56)

91/12/04 保坂征宏

2005/04/04 石渡正樹

7等比級数の和を直接計算しても良い.

2008/06/10(地球流体電脳倶楽部)  zahyou/gaussgr.tex

\underline{DCPAM4 第 2 部離散化  A.4 積分評価  11}

ここで, 上に示した「より実用的な公式」により

∑_{n=1}^{I} (exp(\frac{2πim}{I}))^{n−1} = 0 (m ≠ 0)  (A.60)

が成り立つ. したがって,

\frac{1}{2π} ∫_{0}^{2π} g(λ)dλ = \frac{1}{I} ∑_{n=1}^{I} g(λ_{n})  (A.61)

となる.

A.4.2 Gauss-Legendre の公式

f(μ) を 2J − 1 次以下の多項式とする. P_{n} を 2 で規格化した n 次の Legendre 函数とする. このとき, ∫_{−1}^{1} fdμ は P_{J} の零点である Gauss 格子 μ_{j}（j = 1, 2, ⋯, J）における f の値 f(μ_{j}) のみを用いて, 次式にもとづいて正確に評価することができる.

∫_{−1}^{1} f(μ)dμ = 2 ∑_{j=1}^{J} f(μ_{j})w_{j}  (A.62)
w_{j} = \frac{1}{2} ∫_{−1}^{1} \frac{P_{J}(μ)}{(μ − μ_{j})PJ^{′}(μ_{i})}dμ = \frac{(2J − 1)(1 − μ_{j}^{2})}{(JP_{J−1}(μ_{j}))^{2}}.  (A.63)

ここで, w_{j} は Gauss 荷重と呼ばれる.

以下では上の式を証明する. ただし, Legendre 函数としては, 最初は岩波公式集の Legendre 函数 P̃_{n} を用い, 最後に 2 で規格化した Legendre 函数 P_{n} に直すことにする^{8}.

\underline{STEP 1} Lagrange 補間の導入 f(μ) を K 次多項式（0 ≤ K ≤ 2J − 1）とする. P̃_{n} を岩波公式集の Legendre 函数（Rodrigues の公式）とする.

∫_{−1}^{1} P̃_{n}(μ)P̃_{n′}(μ)dμ = \frac{2}{2n + 1}δ_{nn′}  (A.64)

L(μ) を, f(μ_{j}) を Lagrange 補間公式にしたがって補間した多項式として定義する.

L(μ) ≡ ∑_{j=1}^{J} f(μ_{j}) ∏_{k=1,k≠j}^{J} \frac{μ − μ_{k}}{μ_{j} − μ_{k}}  (A.65)

このとき, 各 j について L(μ_{j}) = f(μ_{j}) である. ここで L は, 0 ≤ K ≤ J − 1 の時（f が J − 1 次以下の多項式）のときは厳密に L = f になる^{9} ことに注意せよ.

したがって, 関数 f(μ) − L(μ) は • 0 ≤ K ≤ J − 1 の時, 0 である.

8私が混乱しないためにこのような手続きを踏む. 実際, 公式集を含む他の文献には P̃_{n}^{m} の公式が書かれていることが多いので, このように書く方が他と参照しやすいであろう.

9このことは L − f が J − 1 次以下の多項式であること, J 個の零点 μ_{j} を持つことから明らか.

zahyou/gaussgr.tex  2008/06/10(地球流体電脳倶楽部)

• J ≤ K ≤ 2J − 1 の時, μ = μ_{j} を零点とする K 次多項式である. μ_{j} は J 次多項式 P̃_{J}(μ) の零点であることを思い出すと, f − L は P̃_{J}(μ) で割り切れるので, ある K − J 次多項式 S(μ) を用いて,

f(μ) − L(μ) = P̃_{J}(μ)S(μ)  (A.66)

と書くことができる.

f(μ) − L(μ) を μ について −1 から 1 まで積分する. J ≤ K ≤ 2J − 1 の時については Legendre 函数の直交性より, P̃_{J}(μ)S(μ) の積分は零である. したがって,

∫_{−1}^{1} f(μ)dμ = ∫_{−1}^{1} L(μ)dμ  (A.67)

= ∑_{j=1}^{J} f(μ_{j}) ∫_{−1}^{1} \frac{∏_{k=1}^{J}(μ − μ_{k})}{(μ − μ_{j}) ∏_{k=1,k≠j}^{J}(μ_{j} − μ_{k})}dμ  (A.68)

= ∑_{j=1}^{J} f(μ_{j}) ∫_{−1}^{1} \frac{P̃_{J}(μ)}{(μ − μ_{j})P̃J^{′}(μ_{j})}dμ  (A.69)
= 2 ∑_{j=1}^{J} f(μ_{j})w_{j}  (A.70)

ここで, 証明すべき式の P_{J} は規格化されていて, 上の式の P̃_{J} は規格化されていないのにもかかわらず同じ w_{j} が使われているが, P̃_{J} と P_{J} の規格化定数は同じなので consistent である.

\underline{STEP 2} w_{j} = \frac{1}{2} ∫_{−1}^{1} \frac{P̃_{J}(μ)}{(μ − μ_{j})P̃J^{′}(μ_{j})}dμ の漸化式を用いた変形

漸化式 (岩波の Lgendre 関数・陪関数の従う漸化式) において m = 0 とした式

(n + 1)P̃_{n+1}(μ) = (2n + 1)μP̃_{n}(μ) − nP̃_{n−1}(μ) (n = 0, 1, 2, ⋯)  (A.71)

より,

(n + 1) | \begin{matrix}P̃_{n+1}(x) & P̃_{n}(x) \\ P̃_{n+1}(y) & P̃_{n}(y) \end{matrix} | = | \begin{matrix}(2n + 1)xP̃_{n}(x) − nP̃_{n−1}(x) & P̃_{n}(x) \\ (2n + 1)yP̃_{n}(y) − nP̃_{n−1}(y) & P̃_{n}(y) \end{matrix} |  (A.72)
= (2n + 1)(x − y)P̃_{n}(x)P̃_{n}(y)  (A.73)
+n(−P̃_{n−1}(x)P̃_{n}(y) + P̃_{n−1}(y)P̃_{n}(x))  (A.74)
= (2n + 1)(x − y)P̃_{n}(x)P̃_{n}(y) + n | \begin{matrix}P̃_{n}(x) & P̃_{n−1}(x) \\ P̃_{n}(y) & P̃_{n−1}(y) \end{matrix} | (A.75)

となる. この式を n = 0, 1, ⋯, n − 1 について加えると,

n | \begin{matrix}P̃_{n}(x) & P̃_{n−1}(x) \\ P̃_{n}(y) & P̃_{n−1}(y) \end{matrix} | = ∑_{k=0}^{n−1}(2k + 1)(x − y)P̃_{k}(x)P̃_{k}(y)  (A.76)

が成り立つ. ここで n = J, x = μ, y = μ_{j} とすると P̃_{J}(μ_{j}) = 0 より,

JP̃_{J}(μ)P̃_{J−1}(μ_{j}) = ∑_{k=0}^{J−1}(2k + 1)(μ − μ_{j})P̃_{k}(μ)P̃_{k}(μ_{j}).  (A.77)

よって,

\frac{P̃_{J}(μ)}{μ − μ_{j}} = \frac{∑_{k=0}^{J−1}(2k + 1)P̃_{k}(μ)P̃_{k}(μ_{j})}{JP̃_{J−1}(μ_{j})}  (A.78)

である. したがって,

w_{j} = \frac{1}{2} ∫_{−1}^{1} \frac{P̃_{J}(μ)}{(μ − μ_{j})P̃J^{′}(μ_{j})}dμ  (A.79)
= \frac{1}{2JP̃_{J−1}(μ_{j})P̃J^{′}(μ_{j})} ∑_{k=0}^{J−1}(2k + 1)P̃_{k}(μ_{j}) ∫_{−1}^{1} P̃_{k}(μ)dμ  (A.80)
= \frac{1}{JP̃_{J−1}(μ_{j})P̃J^{′}(μ_{j})}  (A.81)

(A.82)

である. ただし, (A.80) における積分は, k=0 の時のみ 0 でない値を持つこと, および P̃_{0} = 1 を使った. さらに, 漸化式

(1 − μ^{2})\frac{∂P̃_{n}}{∂μ} = nP̃_{n−1}(μ) − nμP̃_{n}(μ)  (A.83)

で n = J, μ = μ_{j} とする. P̃_{J}(μ_{j}) = 0 より,

w_{j} = \frac{1 − μ_{j}^{2}}{(JP̃_{J−1}(μ_{j}))^{2}}  (A.84)

となる.

\underline{STEP3} P̃_{n} の規格化 P_{n} を

∫_{−1}^{1} P_{n}(μ)P′_{n}(μ)dμ = 2  (A.85)

になるように規格化する. P̃_{J−1} = \sqrt{\frac{1}{2(J − 1) + 1}}P_{J−1} より,

w_{j} = \frac{1 − μ_{j}^{2}}{(J\sqrt{\frac{1}{2J−1}}P_{J−1}(μ_{j}))^{2}} = \frac{(2J − 1)(1 − μ_{j}^{2})}{(JP_{J−1}(μ_{j}))^{2}}  (A.86)

となる.

\underline{まとめ} \underline{14 DCPAM4 第 2 部離散化付録 A 座標系・変換公式} 以上より

zahyou/gaussgr.tex  2008/06/10(地球流体電脳倶楽部)

∫_{−1}^{1} f(μ)dμ = 2 ∑_{j=1}^{J} f(μ_{j})w_{j}  (A.87)
w_{j} = \frac{(2J − 1)(1 − μ_{j}^{2})}{(JP_{J−1}(μ_{j}))^{2}}  (A.88)

\underline{DCPAM4 第 2 部離散化  A.5 球面調和函数の離散的直交関係  15}

A.5 球面調和函数の離散的直交関係

ここでは球面直交関数の離散的直交関係である選点直交性を示す.

∑_{j=1}^{J} ∑_{i=1}^{I} P_{n}^{m}(μ_{j})P_{n′}^{m′}(μ_{j}) exp(imλ_{i}) exp(−im′λ_{i})w_{j} = Iδ_{nn′}δ_{mm′}  (A.89)

ここで, i, j, m, m′, n, n′, I, J, M, N(m) は整数で, 1 ≤ i ≤ I, 1 ≤ j ≤ J, 0 ≤ |m|, |m′| ≤ M, |m| ≤

n ≤ N, |m′| ≤ n′ ≤ N であり, M ≤ [\frac{I}{2}] , N(m) ≤ J − 1 を満たす. また, w_{j} は Gauss 荷重,

λ_{i} = \frac{2π(i − 1)}{I}, μ_{j} は P_{J}(μ) の零点である. [ ] はそれを越えない最大の整数を表す. これは, 有限な直交多項式系において成り立つ選点直交性と呼ばれる性質である^{10}.

この式を証明する. Legendre 函数・陪函数の定義・（連続系での）直交性, Gauss の台形公式, Legendre 函数の零点を用いた多項式の積分評価を既知とすると,

∑_{j=1}^{J} ∑_{i=1}^{I} P_{n}^{m}(μ_{j})P_{n′}^{m′}(μ_{j}) exp(imλ_{i}) exp(−im′λ_{i})w_{j}  (A.90)
= I ∑_{j=1}^{J} P_{n}^{m}(μ_{j})P_{n′}^{m′}(μ_{j})w_{j}δ_{mm′}  (A.91)

ここで Gauss の台形公式を使った. 更に変形すると

∑_{j=1}^{J} ∑_{i=1}^{I} P_{n}^{m}(μ_{j})P_{n′}^{m′}(μ_{j}) exp(imλ_{i}) exp(−im′λ_{i})w_{j}  (A.92)
= I ∑_{j=1}^{J} P_{n}^{m}(μ_{j})P_{n′}^{m}(μ_{j})w_{j}  (A.93)
= \frac{I}{2} ∫_{−1}^{1} P_{n}^{m}(μ)P_{n′}^{m}(μ)dμ  (A.94)

ここで, Gauss-Legendre の公式を使った. 更に, 連続系の Legendre 函数・陪函数の直交性より

∑_{j=1}^{J} ∑_{i=1}^{I} P_{n}^{m}(μ_{j})P_{n′}^{m′}(μ_{j}) exp(imλ_{i}) exp(−im′λ_{i})w_{j}  (A.95)
= Iδ_{nn′}δ_{mm′}  (A.96)

が得られる. 以上により, 離散化した球面調和関数の選点直交性が示された.

余談ではあるが, 直交多項式系においては離散的な直交関係としては選点直交性のほかに次のような直交関係も知られている^{11}. {f_{k}(μ)}(k = 0, 1, 2, ⋯) を [a, b] で定義された重み w(μ), 規格化定数 λ_{k} の直交多項式 (∫_{a}^{b} f_{k}(μ)f_{k′}(μ)w(μ)dμ = λ_{k}δ_{kk′}) とする. μ_{j}, μ_{j′}(1 ≤ j, j′ ≤ J) を f_{J}(μ)

91/11/12 保坂征宏

92/01/07 保坂征宏

2005/04/04 石渡正樹

10別の離散的直交関係については後で述べる.

11以下については, 森, 1984 「数値解析法」が詳しい.

zahyou/risan.tex  2008/06/10(地球流体電脳倶楽部)

の零点, w_{j} = w(μ_{j}) とすれば, 選点直交性

∑_{j=0}^{J−1} f_{k}(μ_{j})f_{k′}(μ_{j})w_{j} = λ_{k}δ_{kk′}  (A.97)

のほかに,

∑_{k=0}^{J−1} \frac{f_{k}(μ_{j})f_{k}(μ_{j′})}{λ_{k}} = \frac{1}{w_{j}} δ_{jj′}  (A.98)

が成り立つ.

実際, Legendre 函数 {P_{n}}(n = 0, 1, 2, ⋯, J − 1) についてはこの関係が成り立つ. すなわち, w_{j} を GCM で用いている Gauss 荷重として,

∑_{n=0}^{J−1} P_{n}(μ_{j})P_{n}(μ_{j′}) = \frac{1}{w_{j}} δ_{jj′}  (A.99)

である. しかし, GCM では Legendre 函数 P_{J} の零点でのみ値を計算することと, 波数切断の関係とから, Legendre 陪函数 {P_{n}^{m}}(n = |m|, |m| + 1, |m| + 2, ⋯, N) の離散的直交関係は意味がない^{12}. Legendre 函数の直交関係についても, 波数切断により P_{n} は n = 0, 1, 2, ⋯, N < J − 1 しか扱わないので^{13} 実際には意味がない.

三角関数についても同様な離散的直交関係がある. 選点直交性

∑_{i=0}^{I−1} exp(imλ_{i}) exp(−im′λ_{i}) = Iδ_{mm′}  (A.100)

のほかに,

∑_{m=− \frac{I}{2} +1}^{\frac{I}{2}} exp(imλ_{i}) exp(−imλ_{i′}) = Iδ_{ii′}  (A.101)

も成り立つ. （ただし, I は偶数で I = 2M. I が奇数の場合には, I = 2M + 1 として, m についての和は −\frac{I − 1}{2} ∼ \frac{I − 1}{2} でとる. ）しかし GCM では, 波数切断により |m| の最大値 M は\frac{I}{3} 以下の値なのでやはり意味がない^{14}.

12そもそも, ここで述べている直交関係は f_{k}（k = 0, 1, 2, ⋯, K − 1）が k 次多項式であるような直交多項式系において成り立つものである. Legendre 陪函数は m が奇数のときは多項式でないし, m が偶数であっても P_{n}^{m} は n 次多項式であって, n − m 次多項式ではない. その場合にも直交多項式の議論を拡張してここで述べている直交関係を使えるのか, については未調査である.

13T42 ならば, m = 0 で J = 63, N = 42, R21 ならば, m = 0 で J = 63, N = 21, である.

14T42 ならば I = 128 に対して M = 42 , R21 ならば I = 64 に対して M = 21 である.

2008/06/10(地球流体電脳倶楽部)  zahyou/risan.tex

\underline{DCPAM4 第 2 部離散化  A.6 スペクトルの係数と格子点値とのやり取り  17}

A.6 スペクトルの係数と格子点値とのやり取り

ここではスペクトルの係数と格子点値との変換法について述べる. 実際の GCM 計算において必要になるのは

• スペクトルの係数と格子点値との値のやり取り

• 速度の格子点値の発散 D ・渦度 ζ のスペクトルの係数への変換

• 速度ポテンシャル χ, 流線関数 ψ （もとは発散, 渦度）のスペクトルの係数から速度の格子点値の作成

である.

A.6.1 スペクトルの係数と格子点値との値のやり取り

スカラー関数 A(λ, ϕ) の格子点値とスペクトルの係数とのやり取りは以下のとおりである. ただし, 格子点値は A_{ij} (i = 1, 2, ⋯, I, j = 1, 2, ⋯, J) , スペクトルの係数は Ã_{n}^{m} (m = −M, −M + 1, ⋯, M, n = |m|, |m| + 1, ⋯, N(m)) とする.

A_{ij} ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} Ã_{n}^{m}Y_{n}^{m}(λ_{i}, ϕ_{j})  (A.102)
Ã_{n}^{m} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} A_{ij}Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j}  (A.103)
w_{j} = \frac{(2J − 1)(1 − sin^{2} ϕ_{j})}{(JP_{J−1}(sin ϕ_{j}))^{2}}  (A.104)

以後この文書では簡単のために, ∑_{m=−M}^{M} ∑_{n=|m|}^{N} を ∑_{m,n} と, ∑_{i=1}^{I} ∑_{j=1}^{J} を ∑_{i,j} と表記する.

A.6.2 スペクトルの係数と格子点値との値のやり取り〜東西微分編

まず,

g ≡ \frac{∂f}{∂λ}

を考える.

東西微分（λ 微分）は次式で評価する.

g_{ij} ≡ [\frac{∂}{∂λ} (∑_{m,n} f̃_{n}^{m}Y_{n}^{m}(λ, ϕ))]_{ij}  (A.105)

すなわち,

g_{ij} = ∑_{m,n} imf̃_{n}^{m}Y_{n}^{m}(λ_{i}, ϕ_{j})  (A.106)

91/12/04 保坂征宏 zahyou/henkan.tex  2008/06/10(地球流体電脳倶楽部)

\underline{18 DCPAM4 第 2 部離散化付録 A 座標系・変換公式} である. 変換公式 (A.103) で A を g とみなしたものと (A.106) とを比較すれば明らかに^{15},

g̃_{n}^{m} = imf̃_{n}^{m}  (A.107)

よって,

g̃_{n}^{m} = \frac{1}{I} ∑_{i,j} imf_{ij}Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j}  (A.108)

である. 次に,

h ≡ \frac{g}{r cos^{2} ϕ} = \frac{1}{r cos^{2} ϕ} \frac{∂f}{∂λ} [= \frac{∂}{∂x} (\frac{f}{cos ϕ})]

とする. f と h とのやり取りを考える. (A.105) より明らかに,

h_{ij} = \frac{1}{r cos^{2} ϕ_{i}} g_{ij}
∴ h_{ij} = \frac{1}{r cos^{2} ϕ_{j}} ∑_{m,n} imf̃_{n}^{m}Y_{n}^{m}(λ_{i}, ϕ_{j})

一方, (A.107) より

h̃_{n}^{m} = ˜ [\frac{∂}{∂λ} (\frac{f}{r cos^{2} ϕ})]_{n}^{m} = im ˜ (\frac{f}{r cos^{2} ϕ})_{n}^{m}
= \frac{1}{I} ∑_{i,j} im (\frac{f}{r cos^{2} ϕ})_{ij} Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j}
= \frac{1}{I} ∑_{i,j} imf_{ij}Y_{n}^{m∗}(λ_{i}ϕ_{j})\frac{w_{j}}{r cos^{2} ϕ_{j}}  (A.109)

A.6.3 スペクトルの係数と格子点値との値のやり取り〜南北微分編

まず,

p ≡ \frac{∂f}{∂μ}

を考える. 南北微分（μ 微分）は次式で評価する.

p_{ij} ≡ [\frac{∂}{∂μ} (∑_{m,n} f̃_{n}^{m}Y_{n}^{m})]_{ij}  (A.110)

15より正確には, (g_{ij} =) ∑_{m,n} imf̃_{n}^{m}Y_{n}^{m} = ∑_{m,n} g̃_{n}^{m}Y_{n}^{m} の両辺に左から ∑_{i,j} Y_{n}^{m∗}(λ_{i}, ϕ_{j})w_{j} を演算すれば, im′f̃_{n′}^{m′} =

g̃_{n′}^{m′} として得られる.

2008/06/10(地球流体電脳倶楽部)  zahyou/henkan.tex

\underline{DCPAM4 第 2 部離散化  A.6 スペクトルの係数と格子点値とのやり取り  19}

すなわち,

p_{ij} = ∑_{m,n} f̃_{n}^{m} \frac{dP_{n}^{m}}{dμ} |_{j} exp(imλ_{i})  (A.111)

である.

p_{n}^{m} = \frac{1}{I} ∑_{i,j} p_{ij}Y_{n}^{m∗}w_{j}
= \frac{1}{I} ∑_{i,j} ( ∑_{m′,n′} f̃_{n′}^{m′} \frac{dP_{n′}^{m′}}{dμ} |_{j} exp(im′λ_{i})) P_{n}^{m}(μ_{j}) exp(−imλ_{i})w_{j}
= −\frac{1}{I} ∑_{i,j} ( ∑_{m′,n′} f̃_{n′}^{m′}P_{n′}^{m′}(μ_{j}) exp(im′λ_{i})) \frac{dP_{n}^{m}}{dμ} |_{j} exp(−imλ_{i})w_{j}
= −\frac{1}{I} ∑_{i,j} f_{ij} \frac{dP_{n}^{m}}{dμ} |_{j} exp(−imλ_{i})w_{j}

ここで, 2 行目から 3 行目の等号では,

∑_{i=1}^{I} ∑_{j=1}^{J} f_{n′}^{m′}P_{n}^{m}(μ_{j}) exp(imλ_{i}) \frac{dP_{n′}^{m′}}{dμ} |_{j} exp(−im′λ_{i})w_{j}

= − ∑_{i=1}^{I} ∑_{j=1}^{J} f_{n′}^{m′} \frac{dP_{n}^{m}}{dμ} |_{j} exp(−imλ_{i})P_{n′}^{m′}(μ_{j}) exp(im′λ_{i})w_{j}  (A.112)

を用いた^{16}. 次に,

q ≡ (1 − μ^{2})\frac{∂f}{∂μ} = (1 − μ^{2})p

∑_{i} ∑_{j} f_{n′}^{m′}P_{n}^{m}(μ_{j}) exp(imλ_{i}) \frac{dP_{n′}^{m′}}{dμ} |_{j} exp(−im′λ_{i})w_{j}

= I ∑_{j} f_{n′}^{m′}P_{n}^{m}(μ_{j}) \frac{dP_{n′}^{m′}}{dμ} |_{j} w_{j}δ_{mm′} = I ∑_{j} f_{n′}^{m}P_{n}^{m}(μ_{j}) \frac{dP_{n′}^{m}}{dμ} |_{j} w_{j}δ_{mm′}
= \frac{I}{2} ∫_{−1}^{1} f_{n′}^{m}P_{n}^{m}(μ) \frac{dP_{n′}^{m}}{dμ} dμδ_{mm′}

∑_{i} ∑_{j} f_{n′}^{m′}P_{n}^{m}(μ_{j}) exp(imλ_{i}) \frac{dP_{n′}^{m′}}{dμ} |_{j} exp(−im′λ_{i})w_{j}

= − \frac{I}{2} ∫_{−1}^{1} f_{n′}^{m}P_{n′}^{m}(μ) \frac{dP_{n}^{m}}{dμ} dμδ_{mm′}
= −I ∑_{j} f_{n′}^{m}P_{n′}^{m}(μ_{j}) \frac{dP_{n}^{m}}{dμ} |_{j} w_{j}δ_{mm′}
= − ∑_{i} ∑_{j} f_{n′}^{m′}P_{n′}^{m′}(μ_{j}) exp(im′λ_{i}) \frac{dP_{n}^{m}}{dμ} |_{j} exp(−imλ_{i})w_{j}

16この証明は以下のとおりである.

ここで, 部分積分すると zahyou/henkan.tex  2008/06/10(地球流体電脳倶楽部)

\underline{20 DCPAM4 第 2 部離散化付録 A 座標系・変換公式} とする. (A.110) より明らかに,

q_{ij} = (1 − μ_{j}^{2}) ∑_{m,n} f̃_{n}^{m} \frac{dP_{n}^{m}}{dμ} |_{j} exp(imλ_{i})

である. 一方,

q̃_{n}^{m} = \frac{1}{I} ∑_{i,j} q_{ij}Y_{n}^{m∗}w_{j}
= \frac{1}{I} ∑_{i,j} ((1 − μ_{j}^{2}) ∑_{m′,n′} f̃_{n′}^{m′} \frac{dP_{n′}^{m′}}{dμ} |_{j} exp(im′λ_{i})) P_{n}^{m}(μ_{j}) exp(−imλ_{i})w_{j}
= −\frac{1}{I} ∑_{i,j} ( ∑_{m′,n′} f̃_{n′}^{m′}P_{n′}^{m′}(μ_{j}) exp(im′λ_{i}))

× \frac{d}{dμ} ((1 − μ^{2})P_{n}^{m})|_{j} exp(−imλ_{i})w_{j}

= −\frac{1}{I} ∑_{i,j} f_{ij} \frac{d}{dμ} ((1 − μ^{2})P_{n}^{m})|_{j} exp(−imλ_{i})w_{j}

が成り立つ. ここで, ２行目から３行目において,

∑_{i=1}^{I} ∑_{j=1}^{J} f_{n′}^{m′}(1 − μ_{j}^{2})P_{n}^{m}(μ_{j}) exp(imλ_{i}) \frac{dP_{n′}^{m′}}{dμ} |_{j} exp(−im′λ_{i})w_{j}

= − ∑_{i=1}^{I} ∑_{j=1}^{J} f_{n′}^{m′} \frac{d}{dμ} ((1 − μ^{2})P_{n}^{m})|_{j} exp(−imλ_{i})P_{n′}^{m′}(μ_{j}) exp(im′λ_{i})w_{j}

を用いた^{17}.

17この証明は (A.112) の証明と同様である.

2008/06/10(地球流体電脳倶楽部)  zahyou/henkan.tex

A.6.4 速度の格子点値から発散・渦度のスペクトルの係数への変換

速度場を

(u, v) = (\frac{U}{cos ϕ}, \frac{V}{cos ϕ})  (A.113)

とする. ここでは, (U_{ij}, V_{ij}) から D̃_{n}^{m}, ζ̃_{n}^{m} を求める^{18}.

まず,

D = \frac{1}{r(1 − μ^{2})} \frac{∂U}{∂λ} + \frac{1}{r} \frac{∂V}{∂μ}  (A.114)

より,

D̃_{n}^{m} = \frac{1}{I} ∑_{ij} (\frac{im}{r cos^{2} ϕ_{j}} U_{ij}P_{n}^{m} − \frac{1}{r}V_{ij} \frac{dP_{n}^{m}}{dμ} |_{j}) exp(−imλ_{i})w_{j}
= \frac{1}{I} ∑_{ij} (imU_{ij}P_{n}^{m} − V_{ij} cos^{2} ϕ_{j} \frac{dP_{n}^{m}}{dμ} |_{j}) exp(−imλ_{i})\frac{w_{j}}{r cos^{2} ϕ_{j}}  (A.115)

同様に,

ζ = \frac{1}{r(1 − μ^{2})} \frac{∂V}{∂λ} − \frac{1}{r} \frac{∂U}{∂μ}  (A.116)

より,

ζ̃_{n}^{m} = \frac{1}{I} ∑_{ij} (\frac{im}{cos^{2} ϕ}V_{ij}P_{n}^{m} + \frac{1}{r}U_{ij} \frac{dP_{n}^{m}}{dμ} |_{j}) exp(−imλ_{i})w_{j}
= \frac{1}{I} ∑_{ij} (imV_{ij}P_{n}^{m} + U_{ij} cos^{2} ϕ_{j} \frac{dP_{n}^{m}}{dμ} |_{j}) exp(−imλ_{i})\frac{w_{j}}{r cos^{2} ϕ_{j}}  (A.117)

A.6.5 χ, ψ のスペクトルの係数から速度の格子点値への変換

ここでは χ_{n}^{m}, ψ_{n}^{m} から U_{ij}, V_{ij} を求める方法を記す.

まず,

U = −\frac{1 − μ^{2}}{r} \frac{∂ψ}{∂μ} + \frac{1}{r} \frac{∂χ}{∂λ}  (A.118)

より,

U_{ij} = ∑_{m,n} (−\frac{cos^{2} ϕ_{j}}{r} ψ̃_{n}^{m} \frac{dP_{n}^{m}}{dμ} |_{j} + \frac{1}{r}imχ̃_{n}^{m}P_{n}^{m}(sin ϕ_{j})) exp(imλ_{i})
= ∑_{m,n} (−\frac{1}{r}ψ̃_{n}^{m} cos^{2} ϕ_{j} \frac{dP_{n}^{m}}{dμ} |_{j} + \frac{1}{r}imχ̃_{n}^{m}P_{n}^{m}(sin ϕ_{j})) exp(imλ_{i})  (A.119)

18この項の計算については後ろの補遺「μ 微分についての補助計算」を参照せよ.

zahyou/henkan.tex  2008/06/10(地球流体電脳倶楽部)

\underline{22 DCPAM4 第 2 部離散化付録 A 座標系・変換公式} である. 同様に,

V = \frac{1}{r} \frac{∂ψ}{∂λ} + \frac{1 − μ^{2}}{r} \frac{∂χ}{∂μ}  (A.120)

より,

V_{ij} = ∑_{m,n} (\frac{1}{r}imψ̃_{n}^{m}P_{n}^{m}(sin ϕ_{j}) + \frac{cos^{2} ϕ_{j}}{r} χ̃_{n}^{m} \frac{dP_{n}^{m}}{dμ} |_{j}) exp(imλ_{i})
= ∑_{m,n} (\frac{1}{r}imψ̃_{n}^{m}P_{n}^{m}(sin ϕ_{j}) + \frac{1}{r}χ̃_{n}^{m} cos^{2} ϕ_{j} \frac{dP_{n}^{m}}{dμ} |_{j}) exp(imλ_{i})  (A.121)

である.

2008/06/10(地球流体電脳倶楽部)  zahyou/henkan.tex

\underline{DCPAM4 第 2 部離散化  A.7 スペクトルの係数同士の関係  23}

A.7 スペクトルの係数同士の関係

ここではスペクトルの係数同士の便利な公式を挙げておく. g = \frac{∂ f}{∂λ} の時

g̃_{n}^{m} = imf̃_{n}^{m}  (A.122)

h = ∇_{H}^{2}f の時

h̃_{n}^{m} = −\frac{n(n + 1)}{r^{2}} f̃_{n}^{m}  (A.123)

(A.122) については「スペクトルの係数と格子点値とのやり取り」に証明を示した. ここでは, (A.123) について証明しておく. 微分評価の定義より,

h_{ij} = (∇_{H}^{2} ∑_{m,n} f̃_{n}^{m}Y_{n}^{m})|_{ij} = − ∑_{m,n} \frac{n(n + 1)}{r^{2}} f̃_{n}^{m}Y_{n}^{m}|_{ij}

である. ところで,

h_{ij} = ∑_{m,n} h̃_{n}^{m}Y_{n}^{m}|_{ij}

である. この２つの式の右辺に左から ∑_{i,j} Y_{n′}^{m′∗}|_{ij} を演算して比較すると,

f̃_{n′}^{m′} = −\frac{n(n + 1)}{r^{2}} h̃_{n′}^{m′}

を得る.

91/12/23 保坂征宏 zahyou/spkeisuu.tex  2008/06/10(地球流体電脳倶楽部)

A.8 波数切断

GCM では, 物理量を球面調和函数 P_{n}^{m}(sin ϕ) exp(imλ) で展開したり波数空間で計算するときに, 計算資源の都合上, ある一定波数以下の波数のみを考慮して計算する. そのことを波数切断するという^{19}. 以下ではまず, 切断の基礎知識として切断の仕方・流儀を述べ, ついで, 切断における事情を述べた上で切断波数の決め方を記す.

A.8.1 波数切断の仕方

波数切断の仕方については, 東西波数（m） , 南北波数（n−m）のそれぞれの切断の方法にいくつかの流儀がある. 一般によく用いられるものは三角形切断（Triangle） , 平行四辺形切断（Rhomboidal : 偏菱形）と呼ばれるものである. 三角形切断の場合について計算する波数領域を波数平面上に書くと (A.1) のようになる. 平方四辺形切断の場合は, (A.2) である.

[figure: n − m; m; N_{tr}; N_{tr}; n; m; N_{tr}; N_{tr}] 図 A.1: 三角形切断の場合の波数領域

[figure: n − m; m; N_{tr}; N_{tr}; n; m; N_{tr}; 2N_{tr}; N_{tr}] 図 A.2: 平方四辺形切断の場合の波数領域

三角形切断, 平行四辺形切断, という名称は波数平面上（(n, m) 平面）での形状による^{20}.

91/11/12 保坂征宏

19後述するように, 現実的には波数切断を決めると同時に格子点数が決まる. すなわち, 以上の理由は格子点数を大きくとれないことの理由でもある.

20平方四辺形切断には, n の最大値を m の最大値の 2 倍にしないようなとり方もある. 詳しくは五角形切断に関する脚注参照.

2008/06/10(地球流体電脳倶楽部)  zahyou/setudan.tex

\underline{DCPAM4 第 2 部離散化  A.8 波数切断  25}

より一般的な切断方法は五角形切断 ((A.3)) である.

[figure: n − m; m; J; K − M; M; K − J; n; m; J; K; M; K − J] 図 A.3: 五角形切断の場合の波数領域

三角形切断, 平行四辺形切断はそれぞれ, 五角形切断において

• 三角形切断 J = K = M = N_{tr}

• 平行四辺形切断 K = 2N_{tr}, J = M = N_{tr}

であるような特別な場合である^{21}.

三角形切断と平行四辺形切断の違いについて, 世の中では次のように言われている^{22}.

• 三角形切断の水平分解能は, 経度方向のみならず緯度方向にも一定である^{23}. 分解能を上げてスケールの細かい波を表現できるようになった場合を考える. 物理的にスケールの小さい波には指向性がないことと, 水平分解能に方向依存性がないこととは調和的である. また, このことは, ある三角形波数切断した球面調和函数により表現される球面上の分布は極の位置を変えても同じ三角形波数切断した球面調和函数により正確に表現されることの言い替えでもある.

• 平行四辺形切断の場合, 各東西波数について同じだけの南北波数をとれる.

A.8.2 切断波数の決め方

ここでは切断波数と南北格子点数の決め方について記す. これらは切断の仕方を決めた後に, 使用する計算資源がネックになって決まる. その際, FFT の仕様, aliasing の回避, という２つの数値的な事情を考慮した上で決める必要がある.

21単に K = J +M であるものも平方四辺形切断と呼ばれる. だが, 例えば R21 と呼ばれるものは, K = 42, J = M = 21 のものである.

22気象庁予報部, 1982 の p.47 より.

23分解能が緯度方向に変化することについては, 平行四辺形切断に限らず, 三角形切断以外のどれでも起こる.

zahyou/setudan.tex  2008/06/10(地球流体電脳倶楽部)

FFT の仕様の事情というのは, 話は簡単で, 東西方向に「格子 ⇔ スペクトル」変換するために用いる FFT が効率よく動くための格子点数・波数がある^{24}ことである.

一方, aliasing に関する事情は複雑である. ここで扱っているスペクトルモデルでは, 格子点でのみ値を計算している. いわゆるスペクトルを使うのは, 単に格子点上での水平微分項の評価をする時のみである. その意味で, 「微分の評価にのみスペクトルを用いるグリッドモデル」と言ってもよい. そのように受け止めると, 格子点値を”正しく”計算することを目指し, また, 考慮する波数は厳密にスペクトルの係数と格子との変換を行なうことのできる波数, すなわち変換において情報の落ちないだけの波数をとらねばならないように思える. ところが実際には, スペクトルモデル的な配慮 — ある波数以下についてのみ正しく計算し, それ以上の波数については計算しない — により切断波数・格子点数が決められている. また, 後述する理由により情報は（非線形 aliasing のことを考えずとも）必ず落ちてしまうのである^{25}.

さて, 以下では aliasing に関する事情を具体的に述べながら, 切断波数に対する格子点数の決め方を記そう. 球面上に連続分布している物理量を球面調和函数で展開する. ある波数 M, N(m) 以下（例えば, T42 ならば M = 42, N = 42 ）については線形項・非線形項の両方について厳密に計算できるように I, J を決めることを目指す.

M, N を仮に固定したとして, まずは線形項について切断波数以下のスペクトルの係数のわかっている物理量 A を格子点値に変換しさらにスペクトルの係数に正しくもどすことを考える. A は −M ≤ m ≤ M, |m| ≤ n ≤ N(m) の m, n については Ã_{n}^{m} がわかっているとする. 格子点値は, 1 ≤ i ≤ I, 1 ≤ j ≤ J について

A_{ij} ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} Ã_{n}^{m}P_{n}^{m}(sin ϕ_{j}) exp(imλ_{i})  (A.124)

で与えられる. これらの格子点値から逆に Ã_{n}^{m}（ − M ≤ m ≤ M, |m| ≤ n ≤ N）を計算する. 離散化した系での積分を Gauss の公式, Gauss-Legendre の公式で評価すれば,

Ã_{n}^{m} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} A_{ij}P_{n}^{m}(sin ϕ_{j}) exp(−imλ_{i})w_{j}  (A.125)

である. ここで, w_{j} は ϕ_{j} における重みである. A_{ij} の定義を代入すれば,

Ã_{n}^{m} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} ( ∑_{m′=−M}^{M} ∑_{n′=|m′|}^{N} Ã_{n′}^{m′}P_{n′}^{m′}(sin ϕ_{j}) exp(im′λ)) P_{n}^{m}(sin ϕ_{j}) exp(−imλ_{i})w_{j}

= \frac{1}{I} ∑_{m′=−M}^{M} ∑_{n′=|m′|}^{N} Ã_{n′}^{m′} ∑_{i=1}^{I} exp(i(m′ − m)λ) ∑_{j=1}^{J} P_{n}^{m}(sin ϕ_{j})P_{n′}^{m′}(sin ϕ_{j})w_{j}  (A.126)

24コード依存性がある. 通常, 2 のべき乗が好ましいとされる. コードによっては, 2,3,5 のべき乗の積でもよいものもある.

25実際の GCM では格子点値からスペクトルに変換する際に情報は落ちている. したがって, 格子 - スペクトル - 格子という変換を行なうと元にはもどらない.

例えば T42 の場合, 自由度は 1 + (2 × 1 + 1) + ⋯ + (2 × 42 + 1) = 43^{2} = 1849 に対して格子点数は 128 × 64 = 8192 である. R21 の場合も, 自由度は (2 × 21 + 1) × (21 + 1) = 946 に対して, 格子点数は 64 × 64 = 4096 である. すなわち, 3/4 以上の情報は格子点値からスペクトルに変換するときに落ちている.

工夫すれば情報が落ちないうまい方法があるかも知れないが, 今のところ見つけていないし多分見つからない.

もちろん, スペクトル - 格子 - スペクトルという変換では元にもどる（ように決めている） .

2008/06/10(地球流体電脳倶楽部)  zahyou/setudan.tex

となる. この計算が Ã_{n}^{m} を正しく評価している（すなわち元にもどる）ための I, J の条件は, −M ≤ m ≤ M, |m| ≤ n ≤ N を満たす m, n について

∑_{i=1}^{I} exp(i(m′ − m)λ) = Iδ_{mm′}  (A.127)
∑_{j=1}^{J} P_{n}^{m}(sin ϕ_{j})P_{n′}^{m}(sin ϕ_{j})w_{j} = δ_{nn′}  (A.128)

が成り立つことである. 三角関数の和による評価が正しいための条件は, ここに登場する波数 |m′−m| が最大で 2M の値をとるので, Gauss の公式の適用条件より, 格子点数 I が I ≥ 2M + 1 を満たすことである. Legendre 函数の積の和による評価が正しいための条件は, ここに登場する計算が n + n′ 次の多項式^{26}の評価であることから, Gauss - Legendre の公式の適用条件より, 格子点数 J が 2J − 1 ≥ max[n + n′] = 2max[N] を満たすことである. ここで, max[n + n′] は n + n′ の最大値を, max[N] は N の最大値を表す.

ちなみに, 格子点値からスペクトルの係数に変換し格子点値にもどすという立場からすれば, この Gauss-Legendre の公式の適用条件というのが情報を落とさずには済まない理由である^{27}. このことを以下に述べる. 情報を落とさずに格子点値をスペクトルの係数に変換し格子点値にもどすには, あらゆる東西波数について南北方向の格子点数 J と同じだけの個数の Legendre 函数が必要である. 東西波数 m の場合, 登場する Legendre 陪函数の n は n = |m|, |m| + 1, ⋯, |m| + J − 1 である. P_{n}^{m}P_{n′}^{m} の次数は n + n′ であるから, 最大で 2J + 2|m| − 2 である. これが 2J − 1 以下になるのは m = 0 の時のみである. m ≠ 0 の場合は高次の Legendre 函数は計算してはならない. つまり情報を落とさざるをえない^{28}.

改めて M, N を固定するという立場にもどって, 切断波数以下のスペクトルの係数のわかっている物理量 B, C の積からそれらの格子点値を用いて B と C との積（非線形項） A のスペクトルの係数を正しく求めるための I, J の条件を考える.

A = BC  (A.129)
B = ∑_{m=−M}^{M} ∑_{n=|m|}^{N} (B̃_{n}^{m} exp(imλ)) P_{n}^{m}(sin ϕ)  (A.130)
C = ∑_{m=−M}^{M} ∑_{n=|m|}^{N} (C̃_{n}^{m} exp(imλ)) P_{n}^{m}(sin ϕ)  (A.131)

なる物理量 A, B, C があるとする^{29}. B,C の −M ≤ m ≤ M, |m| ≤ n ≤ N におけるスペクトルの係数 B̃_{n}^{m}, C̃_{n}^{m} を用いて A のスペクトルの係数 Ã_{n}^{m} を 0 ≤ m ≤ M, |m| ≤ n ≤ N については正しく計算することを考える.

26ここで, 三角関数の和が Iδ_{mm′} となることを用いた. 一般には（m, m′ の偶奇が一致しない場合には） P_{n}^{m}P_{n′}^{m′} は多項式にならない.

27Gauss の公式の適用条件と情報欠落との関係についてコメントしておく. 格子点数 I が奇数の場合には, スペクトルで同じ情報量を持つためには波数\frac{I−1}{2} までを考慮すればよいので, 情報は欠落しないことは明らかである. 一方, I が偶数の場合には, 情報は欠落させないためには波数\frac{I}{2} が必要であるが, この波数は Gauss の公式の適用条件を満たさない. しかしこの場合にも, （私は根拠を調べていないが, 少なくとも）経験的には FFT および逆 FFT によって格子 - スペクトル - 格子変換によって情報が落ちないことが知られている.

28この事情により, 非線形項の場合を考えてさらに著しく落とすことが必要になることが次節からわかる.

29A, B, C とも実数である. すなわち, B̃_{n}^{m} = B̃_{n}^{m∗}, etc. となっている.

zahyou/setudan.tex  2008/06/10(地球流体電脳倶楽部)

Ã_{n}^{m} ≡ ˜ (BC)_{n}^{m}  (A.132)
= \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} B_{ij}C_{ij}P_{n}^{m}(sin ϕ_{j}) exp(−imλ_{i})w_{j}  (A.133)

= \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} ( ∑_{m′=−M}^{M} ∑_{n′=|m′|}^{N} B̃_{n′}^{m′} exp(im′λ_{i})P_{n′}^{m′}(sin ϕ_{j}))  (A.134)

× ( ∑_{m′′=−M}^{M} ∑_{n′′=|m′′|}^{N} C̃_{n′′}^{m′′} exp(im′′λ_{i})P_{n′′}^{m′′}(sin ϕ_{j})) P_{n}^{m}(sin ϕ_{j}) exp(−imλ_{i})w_{j} (A.135)

= \frac{1}{I} ∑_{m′=−M}^{M} ∑_{n′=|m′|}^{N} ∑_{m′′=−M}^{M} ∑_{n′′=|m′′|}^{N} B̃_{n′}^{m′} C̃_{n′′}^{m′′}  (A.136)
× ∑_{i=1}^{I} exp(i(m′ + m′′ − m)λ_{i}) ∑_{j=1}^{J} P_{n′}^{m′}(sin ϕ_{j})P_{n′′}^{m′′}(sin ϕ_{j})P_{n}^{m}(sin ϕ_{j})w_{j} (A.137)

この計算が Ã_{n}^{m} を 0 ≤ m ≤ M, |m| ≤ n ≤ N について正しく評価しているための, I, J の条件を線形項の場合と同様に考えると, 格子点数 I が I ≥ 3M + 1 を, 格子点数 J が 2J − 1 ≥ max[n + n′ + n′′] = 3max[N] を満たすことである. ここで, max[n + n′ + n′′] は n + n′ + n′′ の最大値を, max[N] は N の最大値を表す.

再び格子点値からスペクトルの係数に変換し格子点値にもどすという立場からすれば, これらの I, J に関する条件から, 南北成分のみならず, 東西成分についても変換によって情報が落ちてしまうことがわかる.

これまでに述べた M, N を固定したときに格子点数 I, J がとらねばならない個数について, 線形項・非線形項の２つの場合のうち条件が厳しいのは, 明らかに非線形項の場合である. この条件以下の格子点数しかとらない場合には, aliasing をおこすことになる.

以上, FFT, aliasing という２つの事情を考えて格子点数と切断波数とは同時に決められる. 具体的手順は以下のとおりである.

1. 波数切断の仕方を決める.

2. FFT のかけやすい数を選ぶ. それを東西格子点数 I とする.

3. 東西方向の波数の最大値 M を M = [\frac{I − 1}{3}] にする. ただし [ ] はそれを越えない最大の整数を表す記号である.

4. 最大全波数 N_{max} を決める. 三角形切断ならば N_{max} = M, 平行四辺形切断ならば N_{max} = 2M である.

5. 南北方向の格子点数 J を J ≥ \frac{3N_{max}+1}{2} を満たす数に選ぶ. （今の GCM では偶数でなくてはならない. ）

2008/06/10(地球流体電脳倶楽部)  zahyou/setudan.tex

\underline{DCPAM4 第 2 部離散化  A.8 波数切断  29}

例えば, T42 の場合には M = 42, N = 42, 東西格子点数 I が 128, 南北格子点数 J が 64 である. R21 の場合には M = 21, N = 42, 東西格子点数 I が 64, 南北格子点数 J が 64 である. 参考までに, 線形モデルの場合について決め方を示しておく.

1. 波数切断の仕方を決める.

2. FFT のかけやすい数を選ぶ. それを東西格子点数 I とする.

3. 東西方向の波数の最大値 M を M = [\frac{I}{2}] にする. ただし [ ] はそれを越えない最大の整数を表す記号である^{30}.

4. 最大全波数 N_{max} を決める. 三角形切断ならば N_{max} = M, 平行四辺形切断ならば N_{max} = 2M である.

5. 南北方向の格子点数 J を J ≥ \frac{2N_{max}+1}{2} を満たす数に選ぶ. 例えば, 三角形切断の場合には, I = 128 とすると, M = 64, N = 64, J = 65 となる. つまり T64 では I = 128, J = 65 である. 平方四辺形切断の場合には, I = 64 とすると, M = 32, N = 64, J ≥ 65 となる. つまり R32 では I = 64, J = 65 でよい^{31}.

30ここで, I が偶数のときについては Gauss の公式の適用条件を越えて最大波数\frac{I}{2} まで計算できるという知識を用いた.

31これらの場合でも, 南北方向の細かい情報は格子 - スペクトル - 格子変換によって落ちていることに注意せよ.

zahyou/setudan.tex  2008/06/10(地球流体電脳倶楽部)

A.9 スペクトルモデルと差分モデル

世の中の多くの GCM の離散化の方法としては, 鉛直方向については必ずレベルと称する差分による離散化を行なうが, 水平方向については, 差分する方法（この方法を用いるモデルをグリッドモデルという）と球面調和函数で展開してその係数の時間変化を計算する方法（力学過程において^{32} この方法を用いるモデルをスペクトルモデルという）とが用いられる. その二つの方法については一長一短がある. ここでは双方の特徴について列挙しておく^{33}.

• スペクトルモデルには水平空間差分の誤差がない. これが位相の遅れがないことに通じる（らしい） .

• もっとも, グリッド間隔 1.875 度（波数 63 相当）以上では, 格子点モデルでの差分誤差も十分小さくなり, ほぼ等しい性能といえる.

• 極は特異点であり, 単純には扱えない^{34}. スペクトルモデルではうまく関数系を選ぶことで困難を回避できる. 格子点法では数値的な技巧が必要である（らしい） .

• 保存量を作ることは出力結果の解釈に使いやすいという物理的な理由と, 数値的な発散をおさえやすいという数値的な理由とにより奨励される. 格子点モデルの場合, 技巧を用いることで保存を維持できる. スペクトルモデルの場合, さほどの技巧を用いることなく保存を維持できる.

• 格子点モデルには非線形不安定がある（aliasing） .

• スペクトルモデルの方が, 空間微分を含まないだけプログラムが簡単になる.

• スペクトル法はグリッド法よりも境界条件の点で柔軟でない.

• スペクトルモデルはグリッドモデルに比べて水蒸気等の局地的な現象の表現には適さないといわれる. もっとも, グリッドのあらい格子点モデルではスペクトルモデルに比べてさして優れているとはいえない.

• スペクトルモデルでは一点の影響が（本来は影響が及ばない）遠く離れた点にも与えられてしまう.

• FFT を用いると, 少なくともある程度の解像度までは, スペクトルモデルの方が格子点モデルよりも速い（らしい） .

ちなみに, dcpam はスペクトルモデルに分類される.

91/12/23 保坂征宏 32adjustment 等の意味をなど考えると, 特に物理過程においては, 格子点で考える方が物理的に当然であるように思う.

そのためであろうか, スペクトルモデルである東大版 GCM でも物理過程を格子点で計算している. 他のスペクトルモデルについてもそうであるかどうかは未調査.

33出典は, スペクトル法による数値予報（その原理と実際）（1.6） 34問題点その１. グリッドモデルでは緯度経度図で等間隔に格子点をとると, 極でも CFL を満たすようにするために, 時間差分を細かくしなければならない. 他は未調査.

2008/06/10(地球流体電脳倶楽部)  zahyou/zagai.tex

第B章謝辞

第B章謝辞

B.1 ライセンス規定

COPYRIGHT を参照ください.

B.2 使用上の注意

dcpam4 は研究・教育の場で用いられることを前提としております. 教育現場においては自由に使用・改変・再配布していただいて結構です. 利用する場合には正式なライセンス規定に従って頂くようお願いします. dcpam4 を利用して得られた科学技術的成果を論文や Web 等にて発表する際には, その旨を記し, リファレンスに挙げて頂きますようお願いします.

引用例 (和文)

森川靖大, 石渡正樹, 土屋貴志, 山田由貴子, 高橋芳幸, 小高正嗣, 堀之内武, 林祥介, DC-PAM 開発グループ, 2007: 惑星大気モデル DCPAM, http://www.gfd-dennou.org/library/dcpam/, 地球流体電脳倶楽部.

引用例 (英文)

Morikawa,Y., Ishiwatari,M., Tsuchiya,T., Yamada,Y., Takahashi,O.Y., Odaka,M., Horinouchi,T., Hayashi,Y.-Y., DCPAM Devlopment Group, 2007: DCPAM: planetary atmosphere model, http://www.gfd-dennou.org/library/dcpam/, GFD Dennou Club.

B.3 開発グループメンバー

B.3.1 2007 年度

プログラム製作森川靖大, 石渡正樹

プログラム製作協力土屋貴志, 山田由貴子, 高橋芳幸, 小高正嗣, 堀之内武, 林祥介 \underline{32 DCPAM4 第 2 部離散化付録 B 謝辞} 解説文書作成石渡正樹, 森川靖大

2008/06/10(地球流体電脳倶楽部)  zahyou/zagai.tex(../../CREDITS.tex)