dcpam4 支配方程式系とその離散化

地球流体電脳倶楽部

2008 年 06 月 26 日 (dcpam4-20080626)

\underline{ii dcpam4 支配方程式系とその離散化}

目次

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

第 2 章座標系・変換公式 2 2.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 座標系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3 水平格子点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.4 鉛直レベル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.5 水平スペクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.5.1 水平スペクトルの基底の導入 . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.5.2 波数切断 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.5.3 離散化したスペクトルの基底の直交性 . . . . . . . . . . . . . . . . . . . . 5 2.5.4 格子点値とスペクトルの係数との変換法 . . . . . . . . . . . . . . . . . . . 5 2.5.5 内挿公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.5.6 空間微分の評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.6 参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

第 3 章力学過程 8 3.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2 支配方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.2.1 連続の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2.2 静水圧の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2.3 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2.4 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2.5 水蒸気の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2.6 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2.7 波数依存型の水平拡散項 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.3 鉛直離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.1 連続の式, 鉛直速度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.2 静水圧の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3.3 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3.4 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3.5 水蒸気の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.4 水平離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.4.1 連続の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.4.2 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.4.3 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.4.4 水蒸気の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.5 時間積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.5.1 leap frog による時間積分と時間フィルター . . . . . . . . . . . . . . . . . . 16 3.5.2 semi-implicit 時間積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.6 参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

第 4 章積雲パラメタリゼーション 19 4.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 湿潤対流調節 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

4.2.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2.2 水蒸気が少ないという近似を行う場合 . . . . . . . . . . . . . . . . . . . . 20 4.2.3 温度と比湿の調節量の計算方法 . . . . . . . . . . . . . . . . . . . . . . . . 20

4.3 Kuo スキーム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.4 浅い積雲 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.5 荒川シューバートスキーム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.6 参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

第 5 章大規模凝結 26 5.1 大規模凝結 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

第 6 章放射 27 6.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.2 入射放射 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

6.2.1 入射フラックス分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.2.2 年平均日射の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.2.3 同期回転惑星の場合の入射フラックス . . . . . . . . . . . . . . . . . . . . 28

6.3 短波放射 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

6.3.1 短波放射フラックス . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 6.4 長波放射 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

第 7 章地表面過程 30 7.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 7.2 SWANP (未作成) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 7.3 バルク法 (未作成) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

第 8 章蒸発・凝結による地表面気圧変化 31 8.1 基本的な考え方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 8.2 表面気圧変化の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 8.3 比湿の補正 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

第 9 章鉛直フィルター 33 9.1 背景と目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 9.2 基本的な手続き . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 9.3 調節のための基本温度場 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 \underline{iv dcpam4 支配方程式系とその離散化}

9.4 調節する部分の決定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 9.5 調節および誤差の補正 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

第 10 章飽和比湿 36 10.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 10.2 Nakajima et al. (1992) で用いられた式 . . . . . . . . . . . . . . . . . . . . . . . . 36 10.3 参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

第 11 章惑星大気の物理定数 37 11.1 地球大気の物理定数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 11.2 木星大気の物理定数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

付録 A 座標系・変換公式に関する解説 39

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

A.2 微分公式, GCM の変数の微分関係式 . . . . . . . . . . . . . . . . . . . . . . . . . 45 A.2.1 スカラー量の微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 A.2.2 ベクトル量の微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 A.2.3 発散 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 A.2.4 渦度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 A.2.5 速度ポテンシャル, 流線関数と (u, v) . . . . . . . . . . . . . . . . . . . . . 46

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

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

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

A.6.1 スペクトルの係数と格子点値との値のやり取り . . . . . . . . . . . . . . . . 54 A.6.2 スペクトルの係数と格子点値との値のやり取り〜東西微分編 . . . . . . . . 54 A.6.3 スペクトルの係数と格子点値との値のやり取り〜南北微分編 . . . . . . . . 55 A.6.4 χ, ψ のスペクトルの係数から速度の格子点値への変換 . . . . . . . . . . . . 57

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

A.8.1 波数切断の仕方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 A.8.2 切断波数の決め方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

A.9 スペクトルモデルと差分モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 A.10 参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

付録 B 使用上の注意とライセンス規定 66 \underline{dcpam4 支配方程式系とその離散化 v}

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

B.3.1 2007 年度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

\underline{dcpam4 支配方程式系とその離散化  1.1 この文書について  1}

第1章この文書について

1.1 この文書について

この文書は, 地球流体電脳倶楽部で開発中の惑星大気モデル (Dennou-Club Planetary Atmospheric Model) のバージョン 4 である dcpam4 の支配方程式系およびその離散化手法を解説したものである.

dcpam4 は開発中であり, 本文書の内容とソースコードとは必ずしも一致しない.

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

\underline{2 dcpam4 支配方程式系とその離散化第 2 章座標系・変換公式}

第2章座標系・変換公式

2.1 はじめに

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

2.2 座標系

座標系は, 水平方向には緯度 φ, 経度 λ を, 鉛直方向には σ ≡ \frac{p}{p_{s}} をとる. ここで p は気圧, p_{s} は地表面気圧である.

座標の取り方に関する詳細は別紙『支配方程式系の導出に関する参考資料^{1}』の『座標系の取り方』を参照せよ.

2.3 水平格子点

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

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

• 経度方向の格子点

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

1http://www.gfd-dennou.org/library/dcpam/dcpam4/dcpam4 current/doc/physical model/htm/physical model.htm

2以下, J は偶数とする. dcpam4 では, (Gauss 緯度としてとる場合には) J は偶数でなければならない.

3J 次の Legendre 函数 P_{J} (μ) はを満たす J 次多項式であり, P_{J} (μ) の零点は全て −1 < μ < 1 にある. なお, Gauss 緯度は近似的には sin^{−1} (cos \frac{j − 1/2}{J} π) で与えられる.

2008/06/26(地球流体電脳倶楽部)  coordinates/coordinates.tex(coordinates/crd-hdiscrete-gauss.tex)

\underline{dcpam4 支配方程式系とその離散化  2.4 鉛直レベル  3}

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

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

ととる.

2.4 鉛直レベル

Arakawa and Suarez (1983) のスキームを用いる. とり方は以下のとおりである^{4}. 下の層から上へと層の番号をつける. 整数レベルと半整数レベルを定義する^{5}. 半整数レベルでの σ の値 σ_{k−1/2} (k = 1, 2, ⋯, K) を定義する. ただし, レベル\frac{1}{2} は下端 (σ = 1), レベル K + \frac{1}{2} は上端 (σ = 0) とする. 整数レベルの σ の値 σ_{k} (k = 1, 2, … K) は次の式から求める.

σ_{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} は乾燥空気の定圧比熱である^{6}. ま

た, レベル加重 Δσ は以下のように定義される.

\begin{aligned}Δσ_{k} & ≡ σ_{k−1/2} − σ_{k+1/2},  (1 < k < K) \\ Δσ_{1/2} & ≡ σ_{1/2} − σ_{1} = 1 − σ_{1}, \\ Δσ_{K+1/2} & ≡ σ_{K} − σ_{K+1/2} = σ_{K}. \end{aligned}  (2.4)

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

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

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

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

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

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

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

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

• 等温位大気はいつまでも等温位に留まる

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

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

coordinates/coordinates.tex(coordinates/crd-spectral.tex)  2008/06/26(地球流体電脳倶楽部)

2.5 水平スペクトル

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

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

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

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

Y_{n}^{m}(λ, φ) ≡ P_{n}^{m}(sin φ) exp(imλ),  (2.5)

ただし, 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.6)
∫_{−1}^{1} P_{n}^{m}(μ)P_{n′}^{m}(μ)dμ = 2δ_{nn′}  (2.7)

である. なお, P_{n}^{0} を P_{n} とも書く. また sin φ = μ であることを再掲しておく.

2.5.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) である.

2008/06/26(地球流体電脳倶楽部)  coordinates/coordinates.tex(coordinates/crd-spectral.tex)

\underline{dcpam4 支配方程式系とその離散化  2.5 水平スペクトル  5}

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

球面調和函数と波数切断に関する詳細は, 第 A.1 節および第 A.8 節を参照せよ.

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

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

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

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

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

2.5.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} とは次の変換則に従う^{8}.

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

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

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

7詳しくは第 A.5 節を参照せよ.

8正変換, 逆変換時の係数は consistent に与えてさえいれば問題がない.

9さらに, 実際の計算手続きとしては, 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 まで加えればよい.

coordinates/coordinates.tex(coordinates/crd-spectral.tex)  2008/06/26(地球流体電脳倶楽部)

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

Ã_{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.13)

2.5.5 内挿公式

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

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

として得る.

2.5.6 空間微分の評価

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

• λ 微分

(\frac{∂f}{∂λ} )_{ij} ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} imf̃_{n}^{m}Y_{n}^{m}(λ_{i}, φ_{j}),  (2.15)

\begin{matrix} \\ ( & \frac{∂f}{∂λ} & ) \end{matrix}_{n}^{m} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} imf_{ij}Y_{n}^{m∗}(λ_{i}, φ_{j})w_{j}.  (2.16)

• φ 微分

(\frac{∂f}{∂φ} )_{ij} ≡ ∑_{m=−M}^{M} ∑_{n=|m|}^{N} f̃_{n}^{m} \frac{dP_{n}^{m}}{dφ} |_{j} exp(imλ_{i}),  (2.17)

\begin{matrix} \\ ( & \frac{∂f}{∂φ} & ) \end{matrix}_{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.18)

2.6 参考文献

Arakawa, A., Suarez, M. J., 1983: Vertical differencing of the primitive equations in sigma coordinates. Mon. Wea. Rev., 111, 34–35.

2008/06/26(地球流体電脳倶楽部)  coordinates/coordinates.tex(coordinates/crd-references.tex)

\underline{dcpam4 支配方程式系とその離散化  2.6 参考文献  7}

気象庁予報部, 1982 : スペクトル法による数値予報（その原理と実際） . 気象庁, 111pp. Haltiner, G.J., Williams, R.T., 1980: Numerical Prediction and Dynamic Meteorology (2nd ed.).

John Wiley & Sons, 477pp. 森口, 宇田川, 一松編 ,1956 : 岩波数学公式 I . 岩波書店, 318pp. 森口, 宇田川, 一松編 ,1960 : 岩波数学公式 III . 岩波書店, 310pp. 一松信, 1982 : 数値解析. 朝倉書店, 163pp. 森正武, 1984 : 数値解析法. 朝倉書店, 202pp. 寺沢寛一, 1983 : 自然科学者のための数学概論（増訂版） . 岩波書店, 711pp.

coordinates/coordinates.tex(coordinates/crd-references.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{8 dcpam4 支配方程式系とその離散化第 3 章力学過程}

第3章力学過程

3.1 はじめに

この章では力学過程の支配方程式を記し, その支配方程式の離散化を行う.

ここで述べる力学過程とは, 流体の支配方程式における外力項を除いた部分を指す. 外力項である放射や鉛直乱流拡散や雲などに関する過程については別紙を参照のこと.

離散化については, 空間に関する離散化である鉛直離散化と, 水平離散化の方法ならびに時間に関する離散化を行う.

3.2 支配方程式

ここでは力学過程の支配方程式系を示す. この方程式系の詳細に関しては, Haltiner and Williams (1980) もしくは別紙『支配方程式系の導出に関する参考資料^{1}』の『力学過程の支配方程式系の導出』を参照せよ.

3.2.1 連続の式

\frac{∂π}{∂t} + \frac{u}{a cos φ} \frac{∂π}{∂λ} + \frac{v}{a} \frac{∂π}{∂φ} = −D − \frac{∂σ̇}{∂σ} .  (3.1)

3.2.2 静水圧の式

\frac{∂Φ}{∂σ} = −\frac{RT_{v}}{σ} .  (3.2)

1http://www.gfd-dennou.org/library/dcpam/dcpam4/dcpam4 current/doc/physical model/htm/physical model.htm 2008/06/26(地球流体電脳倶楽部)  dynamics/dynamics.tex(dynamics/dyn-equations.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  3.2 支配方程式  9}

3.2.3 運動方程式

\frac{∂ζ}{∂t} = \frac{1}{a cos φ} \{\frac{∂v_{A}}{∂λ} − \frac{∂(u_{A} cos φ)}{∂φ}\} + D(ζ),  (3.3)
\frac{∂D}{∂t} = \frac{1}{a cos φ} \{\frac{∂u_{A}}{∂λ} + \frac{∂(v_{A} cos φ)}{∂φ}\} − ∇_{σ}^{2}(Φ + R\overline{T} π + KE) + D(D).  (3.4)

3.2.4 熱力学の式

\frac{∂T}{∂t} = −\frac{1}{a cos φ} \{\frac{∂(uT′)}{∂λ} + \frac{∂(vT′ cos φ)}{∂φ}\} + T′D

− σ̇ \frac{∂T}{∂σ} + κT_{v} (\frac{∂π}{∂t} + \frac{u}{a cos φ} \frac{∂π}{∂λ} + \frac{v}{a} \frac{∂π}{∂φ} + \frac{σ̇}{σ} ) + \frac{Q}{C_{p}} + D(T) + D′(v).

(3.5)

3.2.5 水蒸気の式

\begin{aligned}\frac{∂q}{∂t} & = −\frac{1}{a cos φ} \{\frac{∂(uq)}{∂λ} + \frac{∂(vq cos φ)}{∂φ}\} + qD \\ & − σ̇ \frac{∂q}{∂σ} + S_{q} + D(q). \end{aligned}  (3.6)

ここで, 独立変数は以下の通りである.

φ : 緯度 [deg.],  (3.7)
λ : 経度 [deg.],  (3.8)
σ ≡ p/p_{s},  (3.9)
t : 時間 [s].  (3.10)

ここで, p は気圧, p_{s} は地表面気圧である. モデルで時間発展を計算することとなる予報変数は以下の通りである.

π (φ, λ) ≡ ln p_{s},  (3.11)
T (φ, λ, σ) : 気温 [K],  (3.12)
q (φ, λ, σ) : 比湿 [kg kg^{−1}],  (3.13)
ζ (φ, λ, σ) ≡ \frac{1}{a cos φ} \{\frac{∂v}{∂λ} − \frac{∂(u cos φ)}{∂φ}\} : 渦度 [s^{−1}],  (3.14)
D (φ, λ, σ) ≡ \frac{1}{a cos φ} \{\frac{∂u}{∂λ} + \frac{∂(v cos φ)}{∂φ}\} : 発散 [s^{−1}].  (3.15)

ここで u は東西風速, v は南北風速であり, それぞれ (φ, λ, σ) の関数である. 各時間ステップで診

dynamics/dynamics.tex(dynamics/dyn-equations.tex.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{10 dcpam4 支配方程式系とその離散化第 3 章力学過程} 断的に求められる変数は以下の通りである.

Φ ≡ gz : ジオポテンシャル高度 [m^{2} s^{−2}],  (3.16)
σ̇ ≡ \frac{dσ}{dt} ≡ \frac{∂σ}{∂t} + \frac{u}{a cos φ} \frac{∂σ}{∂λ} + \frac{v}{a} \frac{∂σ}{∂φ} + \frac{∂σ}{∂σ} ,  (3.17)
\overline{T} (σ) : 水平平均気温 [K],  (3.18)
T′ (φ, λ, σ) ≡ T − \overline{T} ,  (3.19)
T_{v} (φ, λ, σ) ≡ T {1 + (ϵ_{v}^{−1} − 1) q} ,  (3.20)
T′_{v} (φ, λ, σ) ≡ T_{v} − \overline{T} ,  (3.21)
u_{A} (φ, λ, σ) ≡ (ζ + f)v − σ̇ \frac{∂u}{∂σ} − \frac{RT′_{v}}{a cos φ} \frac{∂π}{∂λ} + F_{λ},  (3.22)
v_{A} (φ, λ, σ) ≡ −(ζ + f)u − σ̇ \frac{∂v}{∂σ} − \frac{RT′_{v}}{a} \frac{∂π}{∂φ} + F_{φ},  (3.23)
∇_{σ}^{2} ≡ \frac{1}{a^{2} cos^{2} φ} \frac{∂^{2}}{∂λ^{2}} + \frac{1}{a^{2} cos φ} \frac{∂}{∂φ} (cos φ\frac{∂}{∂φ} ) ,  (3.24)
KE (φ, λ, σ) ≡ \frac{u^{2} + v^{2}}{2}  (3.25)
D(ζ) :  渦度水平拡散,  (3.26)
D(D) :  発散水平拡散,  (3.27)
D(T) :  温度水平拡散,  (3.28)
D(q) :  比湿水平拡散,  (3.29)
F_{λ} (φ, λ, σ) :  小規模運動過程 (経度方向),  (3.30)
F_{φ} (φ, λ, σ) :  小規模運動過程 (緯度方向),  (3.31)
Q (φ, λ, σ) :  放射, 凝結, 小規模運動過程等による加熱・温度変化, (3.32)
S_{q} (φ, λ, σ) :  凝結, 小規模運動過程等による水蒸気ソース,  (3.33)
D′ (v) :  摩擦熱.  (3.34)

各水平拡散 (3.26)〜(3.29) に関しては 3.2.7 節で説明される. 定数は以下の通りである.

a : 惑星半径 [m],  (3.35)
R : 乾燥大気の気体定数 [J kg^{−1} K^{−1}],  (3.36)
C_{p} : 乾燥大気の大気定圧比熱 [J kg^{−1} K^{−1}],  (3.37)
f : コリオリパラメータ [s^{−1}],  (3.38)
κ ≡ R/C_{p},  (3.39)
ϵ_{v} : 水蒸気分子量比.  (3.40)

3.2.6 境界条件

鉛直流に関する境界条件は

σ̇ = 0 at σ = 0, 1.  (3.41)
2008/06/26(地球流体電脳倶楽部)  dynamics/dynamics.tex(dynamics/dyn-equations.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  3.3 鉛直離散化  11}

である. よって (3.1) から, 地表気圧の時間変化式と σ 系での鉛直速度 σ̇ を求める診断式

\frac{∂π}{∂t} = − ∫_{0}^{1} v_{H} · ∇_{σ}πdσ − ∫_{0}^{1} Ddσ,  (3.42)
σ̇ = −σ \frac{∂π}{∂t} − ∫_{0}^{σ} Ddσ − ∫_{0}^{σ} v_{H} · ∇_{σ}πdσ,  (3.43)

が導かれる. ただし熱的境界条件については第 7 章において記述する.

3.2.7 波数依存型の水平拡散項

水平拡散項は, 次のように ∇^{N_{D}} の形で計算する.

D(ζ) = −K_{HD} [(−1)^{N_{D}/2}∇^{N_{D}} − (\frac{2}{a^{2}} )^{N_{D}/2}] ζ,  (3.44)
D(D) = −K_{HD} [(−1)^{N_{D}/2}∇^{N_{D}} − (\frac{2}{a^{2}} )^{N_{D}/2}] D,  (3.45)
D(T) = −(−1)^{N_{D}/2}K_{HD}∇^{N_{D}}T,  (3.46)
D(q) = −(−1)^{N_{D}/2}K_{HD}∇^{N_{D}}q.  (3.47)

この水平拡散項は計算の安定化のための意味合いが強い. 小さなスケールに選択的な水平拡散を表すため, N_{D} としては, 4∼16 を用いる.

3.3 鉛直離散化

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

3.3.1 連続の式, 鉛直速度

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

ここで,

v · ∇π = \frac{u}{a cos φ} \frac{∂π}{∂λ} + \frac{v}{a} \frac{∂π}{∂φ} .  (3.51)
dynamics/dynamics.tex(dynamics/dyn-vdiscrete-as83.tex.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{12 dcpam4 支配方程式系とその離散化第 3 章力学過程}

3.3.2 静水圧の式

\begin{aligned}Φ_{1} & = Φ_{s} + C_{p}(σ_{1}^{−κ} − 1)T_{v,1} \\ & = Φ_{s} + C_{p}α_{1}T_{v,1}. \end{aligned}  (3.52)
\begin{aligned}Φ_{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} \\ & = C_{p}α_{k}T_{v,k} + C_{p}β_{k−1}T_{v,k−1}. \end{aligned}  (3.53)

ここで,

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

3.3.3 運動方程式

\frac{∂ζ_{k}}{∂t} = \frac{1}{a cos φ} \{\frac{∂v_{A}_{,k}}{∂λ} − \frac{∂(u_{A}_{,k} cos φ)}{∂φ}\} + D(ζ_{k}),  (3.56)
\frac{∂D}{∂t} = \frac{1}{a cos φ} \{\frac{∂u_{A}_{,k}}{∂λ} + \frac{∂(v_{A}_{,k} cos φ)}{∂φ}\} − ∇_{σ}^{2}(Φ_{k} + C_{p}κ̂_{k}\overline{T}_{k}π + (KE)_{k}) + D(D_{k}). (3.57)

ここで,

\begin{aligned}u_{A}_{,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 cos φ} \frac{∂π}{∂λ} + F_{λ,k}, \end{aligned}  (3.58)
\begin{aligned}v_{A}_{,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} \frac{∂π}{∂μ} + F_{φ,k}, \end{aligned}  (3.59)
\begin{aligned}κ̂_{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}} , \end{aligned}  (3.60)
T_{v,k} = T_{k} − \overline{T}_{k},  (3.61)
(KE)_{k} = \frac{v_{k}^{2} + v_{k}^{2}}{2} .  (3.62)

2008/06/26(地球流体電脳倶楽部)  dynamics/dynamics.tex(dynamics/dyn-vdiscrete-as83.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  3.3 鉛直離散化  13}

3.3.4 熱力学の式

\begin{aligned}\frac{∂T_{k}}{∂t} & = −\frac{1}{a cos φ} \{\frac{∂(u_{k}T′_{k})}{∂λ} + \frac{∂(v_{k}T′_{k} cos φ)}{∂φ}\} + Ĥ_{k} \\ & + \frac{Q_{k}}{C_{p}} + D(T_{k}) + D′(v). \end{aligned}  (3.63)

ここで,

Ĥ_{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.64)

\begin{aligned}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}^{κ}} \\ & = a_{k}T_{k} + b_{k−1}T_{k−1}, \end{aligned}  (3.65)

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

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

3.3.5 水蒸気の式

\frac{∂q_{k}}{∂t} = −\frac{1}{a cos φ} \{\frac{∂(u_{k}q_{k})}{∂λ} + \frac{∂(v_{k}q_{k} cos φ)}{∂φ}\} + R_{k} + S_{q,k} + D(q_{k}).  (3.68)

ここで,

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.69)

dynamics/dynamics.tex(dynamics/dyn-hdiscrete-spectral.tex.tex)  2008/06/26(地球流体電脳倶楽部)

3.4 水平離散化

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

3.4.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.70)

ここで,

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

3.4.2 運動方程式

\frac{∂ζ̃_{n}^{m}}{∂t} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} im(v_{A})_{ij}Y_{n}^{m∗}(λ_{i}, φ_{j})\frac{w_{j}}{a cos φ_{j}}

+ \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J}(u_{A})_{ij} cos φ_{j} \frac{∂}{∂φ} Y_{n}^{m∗}(λ_{i}, φ_{j})\frac{w_{j}}{a cos φ_{j}}

+ D̃_{M,n}^{m}ζ̃_{n}^{m},

(3.72)

\frac{∂D̃_{n}^{m}}{∂t} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} im(u_{A})_{ij}Y_{n}^{m∗}(λ_{i}, φ_{j})\frac{w_{j}}{a cos φ_{j}}

− \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J}(v_{A})_{ij} cos φ_{j} \frac{∂}{∂φ} Y_{n}^{m∗}(λ_{i}, φ_{j})\frac{w_{j}}{a cos φ_{j}}

− \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.73)

ここで,

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

2008/06/26(地球流体電脳倶楽部)  dynamics/dynamics.tex(dynamics/dyn-hdiscrete-spectral.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  3.5 時間積分  15}

3.4.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 cos φ_{j}}

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

+ \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.75)

ここで,

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

3.4.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 cos φ_{j}}

+ \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} v_{ij}q_{ij} cos φ_{j} \frac{∂}{∂φ} Y_{n}^{m∗}(λ_{i}, φ_{j})\frac{w_{j}}{a cos φ_{j}}

+ \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.77)

3.5 時間積分

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

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

dynamics/dynamics.tex(dynamics/dyn-tintegral-lf-si.tex.tex)  2008/06/26(地球流体電脳倶楽部)

3.5.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.78)

Ȧ_{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.79)

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

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

と Ā を求める. ϵ_{f} はフィルターの係数であり, 標準的な値である 0.05 を使用する.

3.5.2 semi-implicit 時間積分

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

\frac{∂π}{∂t} = (\frac{∂π}{∂t} )_{NG} − C · D,  (3.81)

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

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

ここで, 非重力波項は,

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

2008/06/26(地球流体電脳倶楽部)  dynamics/dynamics.tex(dynamics/dyn-tintegral-lf-si.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  3.5 時間積分  17}
σ̇_{k−1/2}^{NG} = −σ_{k−1/2} (\frac{∂π}{∂t} )^{NG} − ∑_{l=k}^{K} v_{l} · ∇πΔσ_{l},  (3.85)

(\frac{∂D}{∂t} )^{NG} = \frac{1}{a cos φ} \{\frac{∂u_{A}_{,k}}{∂λ} + \frac{∂(v_{A}_{,k} cos φ)}{∂φ}\} − ∇_{σ}^{2}(KE)_{k} − ∇_{σ}^{2} ∑_{k=1}^{K} W_{kl}(T_{v,l} − T_{l}) + D(D_{k}),

(3.86)

(\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.87)
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}(ˆ\overline{T}_{k−1/2} − \overline{T}_{k}) + σ̇_{k+1/2}^{NG}(\overline{T}_{k} − ˆ \overline{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.88)

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

C_{k} = Δσ_{k},  (3.89)
W_{kl} = C_{p}α_{l}δ_{k≥l} + C_{p}β_{l}δ_{k−1≥l},  (3.90)
G_{k} = κ̂_{k}C_{p}\overline{T}_{k},  (3.91)
\underline{h} = \underline{Q}\underline{S} − \underline{R},  (3.92)
Q_{kl} = \frac{1}{Δσ_{k}} (ˆ\overline{T}_{k−1/2} − \overline{T}_{k})δ_{k=l} + \frac{1}{Δσ_{k}} (\overline{T}_{k} − ˆ \overline{T}_{k+1/2})δ_{k+1=l},  (3.93)
S_{kl} = σ_{k−1/2}Δσ_{l} − Δσ_{l}δ_{k≤l},  (3.94)
R_{kl} = − (\frac{α_{k}}{Δσ_{k}} Δσ_{l}δ_{k≤l} + \frac{β_{k}}{Δσ_{k}} Δσ_{l}δ_{k+1≤l}) \overline{T}_{k}.  (3.95)

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

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

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

δ_{t}π = (\frac{∂π}{∂t} )_{NG} − C · \overline{D}^{t},  (3.98)
dynamics/dynamics.tex(dynamics/dyn-tintegral-lf-si.tex.tex)  2008/06/26(地球流体電脳倶楽部)

δ_{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.99)

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

(3.98), (3.99), (3.100) より,

\{(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.101)

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

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

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

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

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

3.6 参考文献

Arakawa, A., Suarez, M. J., 1983: Vertical differencing of the primitive equations in sigma coordinates. Mon. Wea. Rev., 111, 34–35.

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

Bourke, W.P., 1988: Spectral methods in global climate and weather prediction models. Physically-Based Modelling and Simulation of Climates and Climatic Change. Part I., M.E. Schlesinger (ed.), Kluwer Academic Publishers, Dordrecht, 169–220.

Haltiner, G.J., Williams, R.T., 1980: Numerical Prediction and Dynamic Meteorology (2nd ed.). John Wiley & Sons, 477pp.

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

2008/06/26(地球流体電脳倶楽部)  dynamics/dynamics.tex(dynamics/dyn-references.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  4.2 湿潤対流調節  19}

第4章積雲パラメタリゼーション

4.1 はじめに

大気大循環モデルにおいては積雲を様に表現するだけの分解能を持たないので, 雲の発生する条件並びに雲が大気大循環に与える影響については何らかの方法で評価せざるを得ない. 雲が発生する条件および雲が大気大循環に与える影響のうちの熱・運動量輸送効果については^{1}, 大規模場の速度や熱力学的諸量から評価することが多い. この評価方法は一般に積雲パラメタリゼーションと呼ばれ, 特に以下の型のものが良く用いられる.

• 湿潤対流調節 (Manabe et al., 1965)

• クオスキーム (Kuo, 1974) ^{2}

• 浅い積雲^{3}

• 荒川シューバートスキーム^{4}

また, そもそも大気が過飽和状態にあれば降水が起こる. これを大規模凝結という. これについては別紙『大規模凝結』を参照のこと.

以下では dcpam4 に実装されている各種パラメタリゼーションについて解説する.

4.2 湿潤対流調節

4.2.1 はじめに

連続した 2 つのレベルの間の層において, 次の条件が満たされる場合に調節を行う.

1. 温度減率が湿潤断熱減率よりも大きい.

2. 飽和もしくは過飽和.

1雲が大気大循環に与える他の効果として放射が知られる.

2dcpam4 には今後実装の予定.

3dcpam4 には今後実装の予定.

4dcpam4 には実装されていない.

cumulus/cumulus.tex(cumulus/cumulus-adjust.tex.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{20 dcpam4 支配方程式系とその離散化第 4 章積雲パラメタリゼーション}

4.2.2 水蒸気が少ないという近似を行う場合

上記の条件 (1) に関して, エントロピー S の高度変化に関する条件

\frac{dS}{dz} < 0  (4.1)

を, 「水蒸気が少ない」という近似を用いて書きかえると

\frac{∂T}{∂z} − \frac{RT}{C_{p}p} \frac{∂p}{∂z} + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂z} < 0  (4.2)

となる. ここで, C_{p} は乾燥大気の定圧比熱, L は潜熱, q^{∗} は飽和比湿である. 上記の条件 (2) はそのまま用いる. これらを用いて温度と比湿を調節するのがこのモデルでのデフォルトの湿潤対流調節スキームである. 以下, スキームの定式化の説明を行う^{5}.

4.2.3 温度と比湿の調節量の計算方法

比湿と温度を, (q̂, T̂) から (q, T) へ調節するものとする. 条件式は以下の通りである.

q_{k−1} = q^{∗}(T_{k−1}, p_{k−1}),  (4.3)
q_{k} = q^{∗}(T_{k}, p_{k}),  (4.4)
T_{k−1} − T_{k} + \frac{L}{C_{p}} \{q^{∗}(T_{k−1}, p_{k−1}) − q^{∗}(T_{k}, p_{k})\} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{p_{k−1/2}} \frac{T_{k−1} + T_{k}}{2} = 0,  (4.5)

(C_{p}T_{k} + Lq_{k})Δp_{k} + (C_{p}T_{k−1} + Lq_{k−1})Δp_{k−1} = (C_{p}T̂_{k} + Lq̂_{k})Δp_{k} + (C_{p}T̂_{k−1} + Lq̂_{k−1})Δp_{k−1}. (4.6) 解は以下のようになる.

• q_{k}, q_{k−1} (q ≪ 1 は使っていない)

比湿に関しては, 過飽和の場合には以下のように飽和状態に調節する.

q_{k} = q^{∗}(T̂_{k}, p_{k}) + \frac{∂q^{∗}}{∂T} |T̂_{k} ΔT_{k},  (4.7)
q_{k−1} = q^{∗}(T̂_{k−1}, p_{k−1}) + \frac{∂q^{∗}}{∂T} |T̂_{k−1} ΔT_{k−1}.  (4.8)

ここで,

ΔT_{k} = T_{k} − T̂_{k},  (4.9)
ΔT_{k−1} = T_{k−1} − T̂_{k−1}  (4.10)

である. これらの導出に関しては以下を参照のこと.

5以下は差分法と混ざった話になってしまっているので, あとで連続系の話と離散系の話とに分ける必要がある.

2008/06/26(地球流体電脳倶楽部)  cumulus/cumulus.tex(cumulus/cumulus-adjust.tex.tex)

• ΔT_{k−1} (q ≪ 1 は使っていない) (4.6), (4.7), (4.8) より,

(C_{p}T_{k−1} + Lq_{k−1})Δp_{k−1} − (C_{p}T̂_{k−1} + Lq̂_{k−1})Δp_{k−1} = −(C_{p}T_{k} + Lq_{k})Δp_{k} + (C_{p}T̂_{k} + Lq̂_{k})Δp_{k},
\{ΔT_{k−1} + \frac{L}{C_{p}} (q_{k−1} − q̂_{k−1})\} Δp_{k−1} = − \{ΔT_{k} + \frac{L}{C_{p}} (q_{k} − q̂_{k})\} Δp_{k},

[ΔT_{k−1} + \frac{L}{C_{p}} \{q^{∗}(T̂_{k−1}) + \frac{∂q^{∗}}{∂T} |_{k−1} ΔT_{k−1} − q̂_{k−1}\}] Δp_{k−1}

= − [ΔT_{k} + \frac{L}{C_{p}} \{q^{∗}(T_{k}) + \frac{∂q^{∗}}{∂T} |_{k} ΔT_{k} − q̂_{k}\}] Δp_{k},

\{1 + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T} |_{k−1}\} ΔT_{k−1}Δp_{k−1} + \frac{L}{C_{p}} \{q^{∗}(T̂_{k−1}) − q̂_{k−1}\} Δp_{k−1}

= − [1 + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T} |_{k}] ΔT_{k}Δp_{k} − \frac{L}{C_{p}} \{q^{∗}(T̂_{k}) − q̂_{k}\} Δp_{k} (4.11)

である. ここで

γ_{k} ≡ \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T} |_{k}  (4.12)

とおくと,

{1 + γ_{k−1}} ΔT_{k−1}Δp_{k−1}

= − (1 + γ_{k}) ΔT_{k}Δp_{k} + \frac{L}{C_{p}} [\{q̂_{k−1} − q^{∗}(T̂_{k−1})\} Δp_{k−1} + \{q̂_{k} − q^{∗}(T̂_{k})\} Δp_{k}] (4.13)

となる. ここで,

ΔQ̂ ≡ {q̂_{k−1} − q^{∗}(T̂_{k−1})} Δp_{k−1} + {q̂_{k} − q^{∗}(T̂_{k})} Δp_{k}  (4.14)

とおき, ΔT_{k−1} について解けば,

\{1 + γ_{k−1}\} ΔT_{k−1}Δp_{k−1} = − (1 + γ_{k}) ΔT_{k}Δp_{k} + \frac{L}{C_{p}} ΔQ̂,
ΔT_{k−1} = −\frac{1 + γ_{k}}{1 + γ_{k−1}} \frac{Δp_{k}}{Δp_{k−1}} ΔT_{k} + \frac{1}{1 + γ_{k−1}} \frac{L}{C_{p}} ΔQ̂\frac{1}{Δp_{k−1}} (4.15)

となる. ΔT_{k} の導出については以下を参照のこと.

• ΔT_{k}

cumulus/cumulus.tex(cumulus/cumulus-adjust.tex.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{22 dcpam4 支配方程式系とその離散化第 4 章積雲パラメタリゼーション}

(4.5), (4.7), (4.8), (4.9), (4.10) より,

T_{k−1} − T_{k} + \frac{L}{C_{p}} \{q^{∗}(T_{k−1}, p_{k−1}) − q^{∗}(T_{k}, p_{k})\} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{p_{k−1/2}} \frac{T_{k−1} + T_{k}}{2} = 0,

T̂_{k−1} + ΔT_{k−1} − (T̂_{k} + ΔT_{k}) + \frac{L}{C_{p}} \{q^{∗}(T̂_{k−1}) + \frac{∂q^{∗}}{∂T} |_{k−1} ΔT_{k−1} − q^{∗}(T̂_{k}) − \frac{∂q^{∗}}{∂T} |_{k} ΔT_{k}\}

− \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{p_{k−1/2}} \frac{T̂_{k−1} + ΔT_{k−1} + T̂_{k} + ΔT_{k}}{2} = 0,

T̂_{k−1} − T̂_{k} + \frac{L}{C_{p}} \{q^{∗}(T̂_{k−1}) − q^{∗}(T̂_{k})\} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{p_{k−1/2}} \frac{T̂_{k−1} + T̂_{k}}{2}

+ ΔT_{k−1} − ΔT_{k} + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T} |_{k−1} ΔT_{k−1} − \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T} |_{k} ΔT_{k} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{p_{k−1/2}} \frac{ΔT_{k−1} + ΔT_{k}}{2} = 0. (4.16)

ここで,

S_{k−1/2} ≡ T̂_{k−1} − T̂_{k} + \frac{L}{C_{p}} \{q^{∗}(T̂_{k−1}) − q^{∗}(T̂_{k})\} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{p_{k−1/2}} \frac{T̂_{k−1} + T̂_{k}}{2}  (4.17)

とおく. すると湿潤不安定の条件は

S_{k−1/2} > 0  (4.18)

で表される. (4.17) を用いて (4.16) を変形すると,

ΔT_{k−1} − ΔT_{k} + γ_{k−1}ΔT_{k−1} − γ_{k}ΔT_{k} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{2p_{k−1/2}} ΔT_{k−1} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{2p_{k−1/2}} ΔT_{k} = −S_{k−1/2},

(−1 − γ_{k} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{2p_{k−1/2}} ) ΔT_{k} + (1 + γ_{k−1} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{2p_{k−1/2}} ) ΔT_{k−1} = −S_{k−1/2} (4.19)

となる. ここで (4.15) を用いると,

(−1 − γ_{k} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{2p_{k−1/2}} ) ΔT_{k}

+ (1 + γ_{k−1} − \frac{R}{C_{p}} \frac{Δp_{k−1/2}}{2p_{k−1/2}} ) (−\frac{1 + γ_{k}}{1 + γ_{k−1}} \frac{Δp_{k}}{Δp_{k−1}} ΔT_{k} + \frac{1}{1 + γ_{k−1}} \frac{L}{C_{p}} ΔQ̂\frac{1}{Δp_{k−1}} )

= −S_{k−1/2}  (4.20)

κ = R/C_{p} を使うと^{6},

(−1 − γ_{k} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) ΔT_{k}

+ (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) (−\frac{1 + γ_{k}}{1 + γ_{k−1}} \frac{Δp_{k}}{Δp_{k−1}} ΔT_{k} + \frac{1}{1 + γ_{k−1}} \frac{L}{C_{p}} ΔQ̂\frac{1}{Δp_{k−1}} )

= −S_{k−1/2}  (4.21)

となる. ΔT_{k} に関する項を左辺にまとめると,

[(−1 − γ_{k}) − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} + (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) (−\frac{1 + γ_{k}}{1 + γ_{k−1}} \frac{Δp_{k}}{Δp_{k−1}} )] ΔT_{k}

= −S_{k−1/2} − (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) \frac{1}{1 + γ_{k−1}} \frac{L}{C_{p}} ΔQ̂\frac{1}{Δp_{k−1}}  (4.22)

6もっと前から使うべき?? (コードとの対応によっては κ を使わない方がよいのかも)

2008/06/26(地球流体電脳倶楽部)  cumulus/cumulus.tex(cumulus/cumulus-adjust.tex.tex)

dcpam4 支配方程式系とその離散化  4.3 Kuo スキーム  23

となる. ΔT_{k} について解くと, ΔT_{k}

= \frac{−S_{k−1/2} − (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) \frac{1}{1 + γ_{k−1}} \frac{L}{C_{p}} ΔQ̂\frac{1}{Δp_{k−1}}}{(−1 − γ_{k}) − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} + (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) (−\frac{1 + γ_{k}}{1 + γ_{k−1}} \frac{Δp_{k}}{Δp_{k−1}} )}

= \frac{S_{k−1/2} + (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) \frac{1}{1 + γ_{k−1}} \frac{L}{C_{p}} ΔQ̂\frac{1}{Δp_{k−1}}}{(1 + γ_{k}) + κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} + (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) \frac{1 + γ_{k}}{1 + γ_{k−1}} \frac{Δp_{k}}{Δp_{k−1}}}

= \frac{(1 + γ_{k−1})Δp_{k−1}S_{k−1/2} + (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) \frac{L}{C_{p}} ΔQ̂}{(1 + γ_{k}) (1 + γ_{k−1})Δp_{k−1} + κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} (1 + γ_{k−1})Δp_{k−1} + (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) (1 + γ_{k})Δp_{k}}

= \frac{(1 + γ_{k−1})Δp_{k−1}S_{k−1/2} + (1 + γ_{k−1} − κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} ) \frac{L}{C_{p}} ΔQ̂}{(1 + γ_{k}) (1 + γ_{k−1})(Δp_{k−1} + Δp_{k}) + κ\frac{Δp_{k−1/2}}{2p_{k−1/2}} \{(1 + γ_{k−1})Δp_{k−1} − (1 + γ_{k})Δp_{k}\}}

(4.23)

が得られる.

4.3 Kuo スキーム

Kuo スキームにおいては, 雲底 k_{B} と雲頂 k_{T} の間で,

T_{k} = T̂_{k} + \frac{L}{C_{p}} (1 − b)I \frac{(T_{k}^{c} − T_{k})}{∑_{k=k_{B}}^{k_{T}}(T_{k}^{c} − T_{k})\frac{Δp_{k}}{g}}  (4.24)
q_{k} = q̄_{k}^{t−Δt} + bI \frac{(q_{k}^{c} − q_{k})}{∑_{k=k_{B}}^{k_{T}}(q_{k}^{c} − q_{k})\frac{Δp_{k}}{g}}  (4.25)

ここで, I は, t − Δt から t + Δt の間の雲底雲頂間の水蒸気収束量であり,

I = ∑_{k=k_{B}}^{k_{T}}(q̂_{k} − q̄_{k}^{t−Δt})\frac{Δp_{k}}{g} .  (4.26)

T_{k}^{c}, q_{k}^{c} は雲内の温度, 比湿であり後述する. b は水蒸気収束を凝結加熱と加湿に分配するパラメータであり, Anthes(1977) に従って,

b = \{ \begin{matrix}(\frac{1−⟨RH⟩}{1−RH_{c}})n & ⟨RH⟩ ≥ RH_{c} \\ 1 & ⟨RH⟩ < RH_{c} \end{matrix}  (4.27)
cumulus/cumulus.tex(cumulus/cumulus-kuo.tex.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{24 dcpam4 支配方程式系とその離散化第 4 章積雲パラメタリゼーション}

と評価する. RH ≡ q̂/q^{∗}(T̂, p) は相対湿度であり,

⟨RH⟩ = \frac{∑_{k=k_{B}}^{k_{T}} q̂\frac{Δp_{k}}{g}}{∑_{k=k_{B}}^{k_{T}} q^{∗}(T̂, p)\frac{Δp_{k}}{g}} .  (4.28)

雲底 k_{B} は, 地表の気塊を断熱的に持ち上げたときに凝結しはじめる高度 (LCL; 持ち上げ凝結高度) より上で成層が湿潤不安定であるもっとも下のレベルであるとする. レベル k で (T_{k}^{a}, q_{k}^{a}) である気塊の k + 1 への断熱的な気塊の持ち上げによる (T_{k+1}^{a}, q_{k+1}^{a}) への変化は,

T_{k+1}^{a} = T_{k}^{a} (\frac{p_{k+1}}{p_{k}} )^{κ}  (4.29)
q_{k+1}^{a} = q_{k}^{a}  (4.30)

となる. このとき,

q_{k+1}^{a} > q^{∗}(T_{k+1}^{a}, p_{k+1})  (4.31)

であれば, 過飽和凝結が起こり,

T_{k+1}^{a} ← T_{k+1}^{a} + \frac{L}{C_{p}} \frac{q_{k+1}^{a} − q^{∗}(T_{k+1}^{a}, p_{k+1})}{1 + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T}}  (4.32)

q_{k+1}^{a} ← q_{k+1}^{a} − \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T} \frac{q_{k+1}^{a} − q^{∗}(T_{k+1}^{a}, p_{k+1})}{1 + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T}},  (4.33)

となる.

まず, T_{1}^{a} = T̂_{1}, q_{1}^{a} = q̂_{1} から出発して逐次このような持ち上げ操作を行ない, はじめて凝結が起きたレベルを LCL とする. 成層が湿潤不安定であるかどうかは, T_{k}^{a} = T̂_{k}, q_{k}^{a} = q^{∗}(T̂_{k}, p_{k}) から出発して, T_{k+1}^{a}, q_{k+1}^{a} を求めたとき,

T_{k+1}^{a} > T̂_{k+1}  (4.34)

となるかどうかによって判断する.

雲内の温度, 比湿 T_{k}^{c}, q_{k}^{c} は雲底 k_{B} で T_{k_{B}}^{a} = T̂_{k_{B}}, q_{k_{B}}^{a} = q^{∗}(T̂_{k_{B}}, p_{k_{B}}) から出発して断熱的に持ち上げた気塊の温度, 比湿 T_{k}^{c} = T_{k}^{a}, q_{k}^{c} = q_{k}^{a} とする. q_{k}^{c} = q^{∗}(T_{k}^{c}, p_{k}) である. 雲頂 k_{T} は, 雲底よりも上でこのようにきめた雲内の温度が周辺よりも低くなる, すなわち,

T_{k}^{c} < T̂_{k}  (4.35)

となるような一番下のレベルである.

降水量 P は,

P = (1 − b)I.  (4.36)

である.

2008/06/26(地球流体電脳倶楽部)  cumulus/cumulus.tex(cumulus/cumulus-references.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  4.6 参考文献  25}

4.4 浅い積雲 4.5 荒川シューバートスキーム 4.6 参考文献

Kuo, H.L., 1974: Further studies of the parameterization of the influence of cumulus convection

of large-scale flow. J. Atmos. Sci., 34, 1232–1240. Manabe, S., Smagorinsky, J., Strickler, R.F., 1965: Simulated climatology of a general circula-

tion model with a hydrologic cycle. Mon. Weather Rev., 93, 769–798.

cumulus/cumulus.tex(cumulus/cumulus-references.tex.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{26 dcpam4 支配方程式系とその離散化第 5 章大規模凝結}

第5章大規模凝結

5.1 大規模凝結

各格子点において, 比湿が飽和比湿よりも大きい場合には, 大規模凝結が起こって飽和状態に調節される. 凝結した水は速やかに降水となって落下し, 雨水の蒸発は考えない. この大規模凝結は, 時間積分の後, 調節として積雲対流のパラメタリゼーションに引き続いて行なわれる. 大規模凝結の起こる条件は,

q̂ > q^{∗}(T̂, p)  (5.1)

である. ここで, T̂, q̂ は時間積分してまだ調節の完全に済んでいない温度, 比湿である. 調節は,

q = q^{∗}(T, p)  (5.2)

を満たし, かつ

C_{p}T + Lq = C_{p}T̂ + Lq̂  (5.3)

となるように行なう. Newton 法によって,

T = T̂ + \frac{L}{C_{p}} \frac{q̂ − q^{∗}(T̂, p)}{1 + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T}}  (5.4)
q = q̂ − \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T} \frac{q̂ − q^{∗}(T̂, p)}{1 + \frac{L}{C_{p}} \frac{∂q^{∗}}{∂T}}.  (5.5)

また, このときの降水量 P は,

P = \frac{Δp}{g}(q̂ − q).  (5.6)

標準的には, これで求められた T, q を T̂, q̂ と読み換え, さらに 2 回の逐次近似を行なう.

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

\underline{dcpam4 支配方程式系とその離散化  6.2 入射放射  27}

第6章放射

6.1 はじめに

放射過程としては太陽から射出された短波放射と地球において射出された長波放射とに分けてとり扱う.

6.2 入射放射

この節では大気上端における中心星 (太陽系惑星の場合は太陽) からの入射放射を与える式についての解説をおこなう.

6.2.1 入射フラックス分布

大気上端における入射放射フラックスの分布の式を書きくだす.

入射フラックス F_{S}^{I} は, 太陽定数を S_{0}, 太陽地球間の距離の, その時間平均値との比を r_{S}, 入射角を ζ とすると,

F_{S}^{I}(λ, φ) = −S_{0}r_{S}^{−2} cos ζ.  (6.1)

ζ は次の式で与えられる.

cos ζ = cos φ cos δ_{S} cos H + sin φ sin δ_{S}  (6.2)

δ_{S} は太陽の赤経, H は時角 (地方時から π を引いたもの) である.

6.2.2 年平均日射の場合

年平均入射量および年平均入射角は, 近似的に, 次のようになる.

\overline{F_{S}^{I}}(φ) ≃ −S_{0}(A_{ins} + B_{ins} cos^{2} φ),  (6.3)
\overline{cos ζ} ≃ A_{ζ} + B_{ζ} cos^{2} φ.  (6.4)
radiation/radiation.tex(radiation/radiation-incoming.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{28 dcpam4 支配方程式系とその離散化第 6 章放射}

\begin{matrix}惑星名 & A_{ins} & B_{ins} & A_{ζ} & B_{ζ} \\ 地球 & 0.127 & 0.183 & 0.410 & 0.590 \end{matrix} 表 6.1: 各惑星における A^{ins}, B^{ins}, A^{ζ}, B^{ζ} の値

大気上端におけるアルベド A を考慮すると

\overline{F_{S}^{I}}(φ) ≃ −S_{0}(1 − A)(A_{ins} + B_{ins} cos^{2} φ)  (6.5)

となる. また, モデルで使用する場合\overline{cos ζ} よりも\overline{sec ζ} の方が便利である. ^{1} \overline{sec ζ} の式は

\overline{sec ζ} ≃ \frac{1}{A_{ζ} + B_{ζ} cos^{2} φ}.  (6.6)

A_{ins}, B_{ins}, A_{ζ}, B_{ζ} の値を表 6.1 に示す.

6.2.3 同期回転惑星の場合の入射フラックス

δ を外から与える. これは太陽直下点の緯度. 次に太陽直下点の経度 (degree) λ_{subsolar} を与える. これにより時角は

H = λ − λ_{subsolar} ∗ \frac{180}{π}  (6.7)

となる.

cos ζ = sin φ sin δ + cos φ · \sqrt{1 − sin^{2} δ} · cos H  (6.8)

入射フラックス分布は

F_{S}^{I}(λ, φ) = −S_{0}(1 − A) cos ζ  (6.9)

6.3 短波放射

6.3.1 短波放射フラックス

短波放射過程においては, 水蒸気とそれ以外の大気による吸収のみを考慮し多重散乱は考慮しない. 吸収係数の異なった N_{S} 個の波長帯を考える (k-distribution method). F_{S} は,

F_{S}(z) = ∑_{i}^{N_{S}} a_{i} [(1 − α_{A})F_{S}^{I} exp (−τ_{S,i}(z) sec ζ)
−α_{g}(1 − α_{A})F_{S}^{I} exp (−τ_{S,i}(0) sec ζ) exp (−(τ_{S,i}(0) − τ_{i}(z)) sec ζ_{0})] (6.10)

1(2007-05-23 石渡) なんでだっけ?

2008/06/26(地球流体電脳倶楽部)  radiation/radiation.tex(radiation/radiation-short.tex)

\underline{dcpam4 支配方程式系とその離散化  6.4 長波放射  29}

ここで, F_{S}^{I} は大気上端からの入射, ζ は入射角, ζ_{0} は散乱光の相当入射角で, sec ζ_{0} = 1.66 とする. α_{A} は大気の散乱によるアルベドであり, 一定値を与える. α_{g} は地表面のアルベドである.

τ_{S,i}(z) は, 大気上端を 0 とした光学的厚さであり,

τ_{S,i}(z) = ∫_{z}^{∞} k_{S,i}ρqdz + ∫_{z}^{∞} k̄_{S,i}ρdz  (6.11)

k_{S,i} は波長帯 i の水蒸気に対する吸収係数, k̄_{S,i} は波長帯 i の水蒸気以外の大気に対する吸収係数である. これら吸収係数は z 等に依存しない一定値を与える. a_{i} は波長帯 i の放射エネルギーの全体に対する割合である.

地表面での吸収は,

F_{S}(0) = ∑_{i}^{N_{S}} a_{i}(1 − α_{g})(1 − α_{A})F_{S}^{I} exp(−τ_{S,i}(0) sec ζ),  (6.12)

で与えられる.

6.4 長波放射

長波放射過程においては, 水蒸気とそれ以外の大気による吸収と射出のみを考慮する. 吸収係数の異なった N_{R} 個の波長帯を考える (k-distribution method). F_{R} は,

F_{R}(z) = (πB(T_{g}) − πB(T_{s})) T^{f}(z, 0) + πB(T (z_{T}))T^{f}(z, z_{T}) − ∫_{0}^{z_{T}} \frac{dπB}{dξ} T^{f}(z, ξ)dξ (6.13)

ここで, T^{f}(z_{1}, z_{2}) は, z = z_{1}, z_{2} 間のフラックス透過関数, πB ≡ σ_{SB}T^{4} は放射源関数である.

フラックス透過関数,T^{f}(z_{1}, z_{2}) は,

T^{f}(z_{1}, z_{2}) = ∑_{i=1}^{N_{R}} b_{i} exp (−δ_{R}|τ_{R,i}(z_{1}) − τ_{R,i}(z_{2})|)  (6.14)

τ_{i}(z) は, 大気上端を 0 とした光学的厚さであり,

τ_{R,i}(z) = ∫_{z}^{∞} k_{R,i}ρqdz + ∫_{z}^{∞} k̄_{R,i}ρdz  (6.15)

k_{R,i} は波長帯 i の水蒸気に対する吸収係数, k̄_{R,i} は波長帯 i の水蒸気以外の大気に対する吸収係数である. これら吸収係数は z 等に依存しない一定値を与える. b_{i} は波長帯 i の放射エネルギーの全体に対する割合であり, 一定値をとると近似する. また, δ_{R} = 1.5 を用いる.

radiation/radiation.tex(radiation/radiation-long.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{30 dcpam4 支配方程式系とその離散化第 7 章地表面過程}

第7章地表面過程

7.1 はじめに

ここでは地表面モデルについてまとめる.

大気大循環モデルにおける地表面のモデルにはいろいろな種類がある. しかし土壌モデル等を含まない現 dcpam4 において, 地表面での過程と言えるものはバルク法による熱や運動量等のやり取りくらいしかない. それでも, 地表面における諸量は幾つもあり, 放射や重力波抵抗等の大気中の現象にも影響を与える. 地表面における諸量には温度, アルベド, 粗度, 熱容量, 湿潤度 (蒸発のしやすさ), 高度分散 (重力波抵抗パラメタリゼーション用) などがある. これらの多くはファイルなどの形でモデルの外から与える必要がある. 地表面の温度は, 海面の場合外から与えられるが, 陸面の場合予報変数である.

7.2 SWANP (未作成)

7.3 バルク法 (未作成)

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

\underline{dcpam4 支配方程式系とその離散化  8.2 表面気圧変化の計算  31}

第8章蒸発・凝結による地表面気圧変化

8.1 基本的な考え方

大気の質量は, 水蒸気の蒸発と凝結 (降水過程) によって変化する. 現在の地球大気の場合には通常この効果は無視することができるけれども, 大気中の水蒸気量が増大した場合には蒸発と凝結による大気質量変化を考慮する必要がある.

蒸発と凝結による大気質量変化をまじめに考慮する場合には, 各高度レベルにおいて質量変化を計算することになる. しかし, そのような取り扱いは σ 座標系においては非常に繁雑になる. dcpam の蒸発・凝結による大気質量変化を計算するモジュールにおいては, 水蒸気の凝結・蒸発による大気質量の変化については, それぞれの鉛直コラムにおけるトータルの質量変化, すなわち蒸発量と凝結量の差の鉛直積分だけを考慮し, 各時間ステップの最後で表面気圧 p_{s} と比湿 q を補正するということをおこなっている.

8.2 表面気圧変化の計算

各時間ステップにおける表面気圧 p_{s} の補正量を Δp_{s}, 凝結と蒸発による比湿変化を Δq とおくと

Δp_{s} = g ∫ ρ̂Δqdẑ = −g ∫ ρ̂Δq \frac{dp̂}{−ρ̂g} = ∫ Δqdp̂ = p̂_{s} ∫ Δqdσ  (8.1)

ただし, p̂_{s} は補正前の表面気圧, ρ̂ は補正前の大気密度である. Δq は蒸発による比湿の増加項と凝結による比湿の減少の項から成る. 順番に考える.

1. 地表面からの蒸発による水蒸気量増加地表面における蒸発フラックスを F_{evap} ∼ LC(q^{∗} − q_{s})ρ_{s}, Δt をタイムステップとすると, 蒸発による最下層の比湿変化 Δq_{evap} は

\frac{Δq_{evap}}{Δt} = \frac{F_{evap}L^{−1}ρ_{s}^{−1}}{Δz} = \frac{F_{evap}L^{−1}ρ_{s}^{−1}}{Δp/(ρ_{s}g)} = F_{evap} \frac{g}{LΔp}  (8.2)

となる. これより

Δq_{evap} = F_{evap} \frac{g}{LΔp}Δt  (8.3)

ps-correct/ps-correct.tex  2008/06/26(地球流体電脳倶楽部)

\underline{32 dcpam4 支配方程式系とその離散化第 8 章蒸発・凝結による地表面気圧変化}

2. 凝結による水蒸気減少

Δq_{cond} = S_{q}^{cond}Δt  (8.4)

ただし, S_{q}^{cond} は凝結による比湿の変化率を表す.

よって,

Δp_{s} = F_{evap} \frac{g}{L}Δt + p̂_{s} ∫ S_{q}^{cond}Δtdσ  (8.5)

となる.

8.3 比湿の補正

上のように, 表面気圧 p_{s} を補正したならば各鉛直レベルにおける比湿の値も補正しなければならない. 補正の前後で水蒸気の質量 (絶対湿度) が変わっているわけではないからである. したがって, 水蒸気の質量が補正の前後で等しくなるように比湿 q を補正する. 補正の前後で水蒸気の質量が等しくなるという条件は以下の式で表される.

ρ̂q̂ = ρq,  (8.6)

ただし, ρ̂, q̂ は補正前の大気密度・比湿の値を, ρ, q は補正後の密度・比湿の値を表す. この式を変形する.

\frac{p̂_{s}σ}{RT} q̂ = \frac{p_{s}σ}{RT} q  (8.7)

p_{s} = p̂_{s} + Δp_{s} を代入すると

\frac{p̂_{s}σ}{RT} q̂ = \frac{(p̂_{s} + Δp_{s})σ}{RT} q,
p̂_{s}q̂ = (p̂_{s} + Δp_{s})q,

(8.8)

よって, 最終的に

q = \frac{p̂_{s}}{p̂_{s} + Δp_{s}} q̂  (8.9)

となる.

2008/06/26(地球流体電脳倶楽部)  ps-correct/ps-correct.tex

\underline{dcpam4 支配方程式系とその離散化  9.4 調節する部分の決定  33}

第9章鉛直フィルター

9.1 背景と目的

太陽定数を増加させた計算を行ったところ 2-grid noise が生じた. このノイズの振幅は非常に大きくなり計算が破綻する.

このノイズを消すためにフィルターをかける. 物理的実体は無い. 強いて言えば, 重力波でつぶれるかもね, 程度である.

9.2 基本的な手続き

1. 調節のための基本温度場を決める.

2. 調節を行う部分を決める. 具体的には, 全層か一部か. 一部としたらどの領域か

3. 温度構造を基本場を調節する. その後, 誤差の補正を行う.

9.3 調節のための基本温度場

ギザギザ成分を取り除いた分布を基本場としたい. そこで「基本場」は次で与えられるとする.

T_{Bk} = \frac{T_{k+1/2} + T_{k−1/2}}{2}  (9.1)

ここで, 1 ≤ k ≤ KMAX

9.4 調節する部分の決定

1. 全層の場合何も考えずに全部やる

2. ギザギザ部分のみ調節する場合以下の式を満たす level k はギザギザな点と判定する.

(T_{k} − T_{k−1}) · (T_{k+1} − T_{k}) < 0  (9.2)

\underline{34 dcpam4 支配方程式系とその離散化第 9 章鉛直フィルター}

実際には, こんな凝ったやってないと思う.

9.5 調節および誤差の補正

節のタイトルが正しくないような気がする. 「誤差の補正」ではないなあ.

1. 全層で調節を行い, 全層に誤差をばらまく場合.

これについては未実装である. これはもっとも安直な調節方法である. 調節前の温度の値を T_{k}, 調節後を T̂_{k} とする. 何も考えずに調節するなら

T̂_{k} = T_{k} + S_{grst}(T_{Bk} − T_{k})  (9.3)

温度変化量は

ΔT_{k} ≡ T̂_{k} − T_{k} = S_{grst}(T_{Bk} − T_{k})  (9.4)

この場合, 全体で

∑_{k=1}^{KMAX} c_{v}ΔT_{k}Δp_{k} = ∑_{k=1}^{KMAX} c_{v}S_{grst}(T_{Bk} − T_{k})(p_{s}Δσ_{k})  (9.5)

の内部エネルギーをコラムに与えてしまうことになる. これが誤差になる. 調節の前後で全体の内部エネルギー量に変化がないようにするためには

∑_{k=1}^{KMAX} ΔT_{k}Δp_{k} = 0  (9.6)

としてやらなければいけない. ただし, ここで比熱が一定の場合を仮定した. これを解決するもっとも安直な方法は上の誤差を全層にばらまくこと. この場合, 温度の調節量は以下の式で与えられる.

T̂_{k} = T_{k} + S_{grst}(T_{Bk} − T_{k}) − \frac{∑_{k=1}^{KMAX} S_{grst}(T_{Bk} − T_{k})Δp_{k}}{∑_{k=1}^{KMAX} Δp_{k}}  (9.7)

こうすると

∑_{k=1}^{KMAX} T̂_{k}Δp_{k} = ∑_{k=1}^{KMAX} \{S_{grst}(T_{Bk} − T_{k})Δp_{k} − \frac{∑_{k=1}^{KMAX} S_{grst}(T_{Bk} − T_{k})Δp_{k}}{∑_{k=1}^{KMAX} Δp_{k}} Δp_{k}\} (9.8)

となる.

2. 全層で調節し, 誤差を局所的に解消していく場合 AGCM5 の場合だと, p2grstA.F

2008/06/26(地球流体電脳倶楽部)  vfilter-conserve/vfilter-conserve.tex

\underline{dcpam4 支配方程式系とその離散化  9.5 調節および誤差の補正  35}

3 点トリオで考え, そこで生じた誤差をその 3 層にばらまくことにする. 下から順に, 以下の式で温度の調節量を計算していく.

T̂_{k−1} = T_{k−1} + S_{grst}(T_{Bk−1} − T_{k−1}) − \frac{∑_{k=k−1}^{k+1} S_{grst}(T_{Bk} − T_{k})Δp_{k}}{∑_{k=k−1}^{k+1} Δp_{k}} , (9.9)
T̂_{k} = T_{k} + S_{grst}(T_{Bk} − T_{k}) − \frac{∑_{k=k−1}^{k+1} S_{grst}(T_{Bk} − T_{k})Δp_{k}}{∑_{k=k−1}^{k+1} Δp_{k}} ,  (9.10)
T̂_{k+1} = T_{k+1} + S_{grst}(T_{Bk+1} − T_{k+1}) − \frac{∑_{k=k−1}^{k+1} S_{grst}(T_{Bk} − T_{k})Δp_{k}}{∑_{k=k−1}^{k+1} Δp_{k}} (9.11)

この T̂ を用いて, 1 つ上に上がり,

ˆT̂_{k} = T̂_{k} + S_{grst}(T_{Bk} − T̂_{k}) − \frac{∑_{k=k}^{k+2} S_{grst}(T_{Bk} − T̂_{k})Δp_{k}}{∑_{k=k}^{k+2} Δp_{k}} ,  (9.12)
ˆT̂_{k+1} = T̂_{k+1} + S_{grst}(T_{Bk+1} − T̂_{k+1}) − \frac{∑_{k=k}^{k+2} S_{grst}(T_{Bk} − T̂_{k})Δp_{k}}{∑_{k=k}^{k+2} Δp_{k}} , (9.13)
T̂_{k+2} = T_{k+1} + S_{grst}(T_{Bk+2} − T_{k+2}) − \frac{∑_{k=k}^{k+2} S_{grst}(T_{Bk} − T̂_{k})Δp_{k}}{∑_{k=k}^{k+2} Δp_{k}}  (9.14)

ただし, 右辺 3 項目の分子の和において k が k + 2 の場合は, T̂_{k+2} ではなく T_{k+2} である (まだ調節を行っていないので). 以上を level を 1 つづつ上がりながら順番に行う.

3. ギザギザ部分だけを調節し誤差を 3 点ごとにばらまく下から順に

(T_{k} − T_{k−1}) · (T_{k+1} − T_{k}) < 0  (9.15)

となる部分の 3 点トリオ出 (2) と同じことを行う.

4. ギザギザ部分だけを調節し, 誤差をギザギザ部分全体に均等にばらまく場合連続してギザギザ部分になっているところを判定し (1) と同じことをする. 下記のように判定する.

(a) 下から順に登っていって初めて

(T_{k} − T_{k−1}) · (T_{k+1} − T_{k}) < 0  (9.16)

となったところで, (k − 1) level がギザギザ部分の底と判定する.

(b) ギザギザ部分で登っていって

(T_{k} − T_{k−1}) · (T_{k+1} − T_{k}) ≥ 0  (9.17)

となったら, k level がギザギザ部分の天井. このように判定されたギザギザ部分全体に誤差をばらまく.

\underline{36 dcpam4 支配方程式系とその離散化第 10 章飽和比湿}

第10章飽和比湿

10.1 はじめに

飽和比湿とは, 与えられた温度圧力のもとで飽和蒸気圧を与える比湿の値である.

飽和比湿を計算するための式として, dcpam4 では Nakajima et al. (1992) で使われた式を用意している. 以下に説明を記す.

10.2 Nakajima et al. (1992) で用いられた式

Nakajima et al. (1992) では, Eisenberg and Kauzmann (1961; この訳書が「水の構造と物性」だと思われる) で与えられている水蒸気の飽和蒸気圧の表と式を近似的に表現する

e^{∗}(T) = p_{0}^{∗} exp (−\frac{L}{RT} )  (10.1)

を用いている. ここで, p_{0}^{∗} = 1.4 × 10^{11} Pa である. この式から飽和比湿の式を書き下すと

q^{∗}(T, p) = p_{0}^{∗} exp (−\frac{L}{RT} ) \frac{1}{p}  (10.2)

である.

10.3 参考文献

Nakajima, S., Hayashi, Y.-Y., Abe, Y., 1992: A study on the “runaway greenhouse effect” with a one dimensional radiative convective equilibrium model. J. Atmos. Sci., 49, 2256–2266.

カウズマン・アイゼンバーグ著, 関集三・松尾隆祐訳, 1975: 水の構造と物性, みすず書房, pp.302.

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

\underline{dcpam4 支配方程式系とその離散化  11.2 木星大気の物理定数  37}

第11章惑星大気の物理定数

11.1 地球大気の物理定数

地球大気の基本的な物理定数を以下に示す.

惑星半径  a  m  6.37 ×10^{6}
重力加速度  g  m s^{−2}  9.8
乾燥大気の定圧比熱  C_{p}  J kg^{−1} K^{−1} 1004.6
乾燥大気の気体定数  R  J kg^{−1} K^{−1} 287.04
蒸発潜熱  L  J kg^{−1}  2.5 ×10^{6}
水蒸気定圧比熱  C_{v}  J kg^{−1} K^{−1} 1810.
水蒸気気体定数  R_{v}  J kg^{−1} K^{−1} 461.
液体水の密度  d_{H_{2}O}  J kg^{−1} K^{−1} 1000.
水蒸気分子量比  ϵ_{v}  0.622
仮温度の係数  δ_{v} = ϵ_{v}^{−1} − 1  0.606
乾燥大気の定圧比熱と気体定数の比 κ = R/C_{p}  0.286
Stefan Bolzman 定数  σ_{SB}  W m^{−2} K^{−4} 5.67 ×10^{−8}
Kálman 定数  k  0.4

11.2 木星大気の物理定数

木星大気の基本的な物理定数を以下に示す. ここでは, まだ凝結成分として H_{2}O しか考慮されていない.

constants/constants.tex(constants/cst-jupiter.tex.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{38 dcpam4 支配方程式系とその離散化第 11 章惑星大気の物理定数}

惑星半径  a  m  1.43 ×10^{10}
重力加速度  g  m s^{−2}  23.1
乾燥大気の定圧比熱  C_{p}  J kg^{−1} K^{−1} 11900.9264
乾燥大気の気体定数  R  J kg^{−1} K^{−1} 3611.44466
蒸発潜熱  L  J kg^{−1}  2.5 ×10^{6}
水蒸気定圧比熱  C_{v}  J kg^{−1} K^{−1} 1810.
水蒸気気体定数  R_{v}  J kg^{−1} K^{−1} 461.
液体水の密度  d_{H_{2}O}  J kg^{−1} K^{−1} 1000.
水蒸気分子量比  ϵ_{v}  7.8250532
仮温度の係数  δ_{v} = ϵ_{v}^{−1} − 1  −0.8722053
乾燥大気の定圧比熱と気体定数の比 κ = R/C_{p}  0.30345912
Stefan Bolzman 定数  σ_{SB}  W m^{−2} K^{−4} 5.67 ×10^{−8}
Kálman 定数  k  0.4

2008/06/26(地球流体電脳倶楽部)  constants/constants.tex(constants/cst-jupiter.tex.tex)

\underline{dcpam4 支配方程式系とその離散化  A.1 球面調和函数  39}

付録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}

\underline{40 dcpam4 支配方程式系とその離散化付録 A 座標系・変換公式に関する解説}

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

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

spectral/spectral.tex(spectral/spl-spherical-harmonics.tex)  2008/06/26(地球流体電脳倶楽部)

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} とも書く.

• 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)

と表される.

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-spherical-harmonics.tex)

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)

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

• 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.9)

ただし, 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}

+ (n + |m|)\sqrt{\frac{1}{2n − 1} \frac{(n − 1 + |m|)!}{(n − 1 − |m|)!}}P_{n−1}^{m} = 0,  (A.10)

∴ P_{n+1}^{m} = \sqrt{\frac{(2n + 1)(2n + 3)}{(n − |m| + 1)(n + |m| + 1)}}μP_{n}^{m}

− \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.11)

ただし, 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.12)

ただし, 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.13)

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

spectral/spectral.tex(spectral/spl-spherical-harmonics.tex)  2008/06/26(地球流体電脳倶楽部)

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

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

と表される.

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

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

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

ただし, 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.17)

すなわち,

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

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

• 完全規格直交性 Y_{n}^{m} は次の直交関係を満たす.

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

ただし, 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.20)
Ã_{n}^{m} = \frac{1}{4π} ∫_{−1}^{1} d(sin φ) ∫_{0}^{2π} dλA(λ, φ)Y_{n}^{m∗}(λ, φ)  (A.21)

と表される.

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-spherical-harmonics.tex)

\underline{dcpam4 支配方程式系とその離散化  A.1 球面調和函数  43}

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

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

• y 微分

• 2 次元ラプラシアンの計算をする. 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.22)

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.23)

2 次元ラプラシアン

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

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

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

(A.24)

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.25)

のである^{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} spectral/spectral.tex(spectral/spl-spherical-harmonics.tex)  2008/06/26(地球流体電脳倶楽部)

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)

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

すなわち, ∇_{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} の解空間の基底である. したがって, 前の座標系の球面調和函数は新たな座標系の球面調和函数の和の形で書ける.

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-spectral-differentiation.tex)

dcpam4 支配方程式系とその離散化  A.2 微分公式, GCM の変数の微分関係式  45

A.2 微分公式, GCM の変数の微分関係式

ここでは, スカラー量, ベクトルの微分を計算する. さらにそれらを元に, 発散 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

= \frac{1}{r^{2}} [\frac{∂}{∂μ} \{(1 − μ^{2})\frac{∂}{∂μ} \} + \frac{1}{1 − μ^{2}} \frac{∂^{2}}{∂λ^{2}} ] f

(A.26)

で与えられる.

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 φ)

= \frac{1}{r\sqrt{1 − μ^{2}}} \frac{∂v_{1}}{∂λ} + \frac{1}{r} \frac{∂}{∂μ} (\sqrt{1 − μ^{2}}v_{2})

(A.27)

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

(rotv)_{r} ≡ \frac{1}{r cos φ} \frac{∂v_{2}}{∂λ} − \frac{1}{r cos φ} \frac{∂}{∂φ} (v_{1} cos φ)

= \frac{1}{r\sqrt{1 − μ^{2}}} \frac{∂v_{2}}{∂λ} − \frac{1}{r} \frac{∂}{∂μ} (\sqrt{1 − μ^{2}}v_{1})

(A.28)

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

spectral/spectral.tex(spectral/spl-spectral-differentiation.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{46 dcpam4 支配方程式系とその離散化付録 A 座標系・変換公式に関する解説}

A.2.3 発散

水平分布する速度場の水平発散 D を u, v を用いて表す

D = \frac{1}{r cos φ} \frac{∂u}{∂λ} + \frac{1}{r cos φ} \frac{∂}{∂φ} (v cos φ).  (A.29)

A.2.4 渦度

水平分布する速度場の渦度 ζ を u, v を用いて表す

ζ = \frac{1}{r cos φ} \frac{∂v}{∂λ} − \frac{1}{r cos φ} \frac{∂}{∂φ} (u cos φ).  (A.30)

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

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

D ≡ ∇_{H}^{2}χ,  (A.31)
ζ ≡ ∇_{H}^{2}ψ  (A.32)

で定義される. (u, v) を χ, ψ で表す.

u = −\frac{1}{r} \frac{∂ψ}{∂φ} + \frac{1}{r cos φ} \frac{∂χ}{∂λ},  (A.33)
v = \frac{1}{r cos φ} \frac{∂ψ}{∂λ} + \frac{1}{r} \frac{∂χ}{∂φ}  (A.34)

となる.

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 の公式の証明に用いられる.

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-legendre.tex)

\underline{dcpam4 支配方程式系とその離散化  A.4 積分評価  47}

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.35)

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

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 法を用いるという方法がある.

A.4 積分評価

A.4.1 Gauss の台形公式

ここでは Gauss の台形公式を示す.

波数 M 以下の三角函数で表現される g(λ) （0 ≤ λ < 2π）

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

spectral/spectral.tex(spectral/spl-gaussgr.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{48 dcpam4 支配方程式系とその離散化付録 A 座標系・変換公式に関する解説} について M < I を満たすように I をとると,

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

λ_{n} = \frac{2π(n − 1)}{I} (n = 1, 2, ⋯, I)

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

∑_{n=1}^{I} exp(imλ_{n}) = \{ \begin{matrix}I & (m = 0), \\ 0 & (0 < |m| < I), \end{matrix}  (A.38)

λ_{n} = \frac{2π(n − 1)}{I} (n = 1, 2, ⋯, I)

である. この証明は, 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 の整数倍にならない) ^{6}. 以下に Gauss の台形公式の証明を記す. まず, 左辺を計算すると,

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

である. ここで, ∫_{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})
= g_{0} + ∑_{m=−M,m≠0}^{M} \frac{g_{m}}{I} ∑_{n=1}^{I} (exp(\frac{2πim}{I}))^{n−1} .  (A.41)

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

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

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

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

となる.

∑_{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.39)

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-gaussgr.tex)

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.44)
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.45)

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

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

\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.46)

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

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

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

したがって, 関数 f(μ) − L(μ) は

• 0 ≤ K ≤ J − 1 の時, 0 である.

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

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

と書くことができる.

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

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

spectral/spectral.tex(spectral/spl-gaussgr.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{50 dcpam4 支配方程式系とその離散化付録 A 座標系・変換公式に関する解説} f(μ) − L(μ) を μ について −1 から 1 まで積分する. J ≤ K ≤ 2J − 1 の時については Legendre 函数の直交性より, P̃_{J}(μ)S(μ) の積分は零である. したがって,

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

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

= ∑_{j=1}^{J} f(μ_{j}) ∫_{−1}^{1} \frac{P̃_{J}(μ)}{(μ − μ_{j})P̃J^{′}(μ_{j})}dμ

= 2 ∑_{j=1}^{J} f(μ_{j})w_{j}  (A.49)

ここで, 証明すべき式の 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.50)

より,

(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} |
= (2n + 1)(x − y)P̃_{n}(x)P̃_{n}(y)
+ n(−P̃_{n−1}(x)P̃_{n}(y) + P̃_{n−1}(y)P̃_{n}(x))
= (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.51)

となる. この式を 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.52)

が成り立つ. ここで 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.53)

よって,

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-gaussgr.tex)

\underline{dcpam4 支配方程式系とその離散化  A.4 積分評価  51}

である. したがって,

w_{j} = \frac{1}{2} ∫_{−1}^{1} \frac{P̃_{J}(μ)}{(μ − μ_{j})P̃J^{′}(μ_{j})}dμ

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

= \frac{1}{JP̃_{J−1}(μ_{j})P̃J^{′}(μ_{j})}  (A.55)

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

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

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

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

となる. \underline{STEP3} P̃_{n} の規格化 P_{n} を

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

になるように規格化する. 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.59)

となる. \underline{まとめ} 以上より

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

spectral/spectral.tex(spectral/spl-spherical-orthogonal.tex)  2008/06/26(地球流体電脳倶楽部)

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.62)

ここで, 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}(μ) の零点である. [ ] はそれを越えない最大の整数を表す. これは, 有限な直交多項式系において成り立つ選点直交性と呼ばれる性質である^{9} .

この式を証明する. 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}

= I ∑_{j=1}^{J} P_{n}^{m}(μ_{j})P_{n′}^{m′}(μ_{j})w_{j}δ_{mm′}.  (A.63)

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

∑_{j=1}^{J} ∑_{i=1}^{I} P_{n}^{m}(μ_{j})P_{n′}^{m′}(μ_{j}) exp(imλ_{i}) exp(−im′λ_{i})w_{j}

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

= \frac{I}{2} ∫_{−1}^{1} P_{n}^{m}(μ)P_{n′}^{m}(μ)dμ.  (A.64)

ここで, 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}

= Iδ_{nn′}δ_{mm′}  (A.65)

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

余談ではあるが, 直交多項式系においては離散的な直交関係としては選点直交性のほかに次のような直交関係も知られている^{10} . {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}(μ)

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

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-spherical-orthogonal.tex)

\underline{dcpam4 支配方程式系とその離散化 A.6 スペクトルの係数と格子点値とのやり取り 53}

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

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

のほかに,

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

が成り立つ.

実際, 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.68)

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

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

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

のほかに,

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

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

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

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

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

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

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

spectral/spectral.tex(spectral/spl-spectral-transform.tex)  2008/06/26(地球流体電脳倶楽部)

• 速度の格子点値の発散 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.71)
Ã_{n}^{m} = \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} A_{ij}Y_{n}^{m∗}(λ_{i}, φ_{j})w_{j},  (A.72)
w_{j} = \frac{(2J − 1)(1 − sin^{2} φ_{j})}{(JP_{J−1}(sin φ_{j}))^{2}}.  (A.73)

以後この文書では簡単のために, ∑_{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.74)

すなわち,

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-spectral-transform.tex)

\underline{dcpam4 支配方程式系とその離散化 A.6 スペクトルの係数と格子点値とのやり取り 55} である. 変換公式 (A.72) で A を g とみなしたものと (A.75) とを比較すれば明らかに^{14} ,

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

よって,

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

である. 次に,

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

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

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.76) より

h̃_{n}^{m} = \begin{matrix} \\ [\frac{∂}{∂λ} (\frac{f}{r cos^{2} φ})] & {m \atop n} \end{matrix} = im\begin{matrix} \\ (\frac{f}{r cos^{2} φ}) & {m \atop n} \end{matrix}

= \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.78)

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

まず,

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

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

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

14より正確には, (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′} として得られる.

spectral/spectral.tex(spectral/spl-spectral-transform.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{56 dcpam4 支配方程式系とその離散化付録 A 座標系・変換公式に関する解説} すなわち,

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

である. よって,

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.81)

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

q ≡ cos^{2} φ\frac{∂f}{∂φ} = cos^{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}.

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

ここで, 部分積分すると 2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-spectral-transform.tex)

\underline{dcpam4 支配方程式系とその離散化 A.6 スペクトルの係数と格子点値とのやり取り 57} とする. (A.79) より明らかに,

q_{ij} = cos^{2} φ_{j} ∑_{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} (cos^{2} φ_{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{d}{dφ} (cos^{2} φP_{n}^{m})|_{j} exp(−imλ_{i})w_{j}

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

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

∑_{i=1}^{I} ∑_{j=1}^{J} f_{n′}^{m′} cos^{2} φ_{j}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φ} (cos^{2} φP_{n}^{m})|_{j} exp(−imλ_{i})P_{n′}^{m′}(φ_{j}) exp(im′λ_{i})w_{j}

を用いた^{16} .

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

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

u = −\frac{1}{r} \frac{∂ψ}{∂φ} + \frac{1}{r cos φ} \frac{∂χ}{∂λ}  (A.82)

より,

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

である. 同様に,

v = \frac{1}{r cos φ} \frac{∂ψ}{∂λ} + \frac{1}{r} \frac{∂χ}{∂φ}  (A.84)

16この証明は (A.81) の証明と同様である.

spectral/spectral.tex(spectral/spl-spectral-transform.tex)  2008/06/26(地球流体電脳倶楽部)

より,

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

である.

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

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

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

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

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

(A.86) については「スペクトルの係数と格子点値とのやり取り」に証明を示した. ここでは, (A.87) について証明しておく.

微分評価の定義より,

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′}

を得る.

A.8 波数切断

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

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-truncation.tex)

\underline{dcpam4 支配方程式系とその離散化  A.8 波数切断  59}

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) 平面) での形状による^{18}.

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

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

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

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

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

spectral/spectral.tex(spectral/spl-truncation.tex)  2008/06/26(地球流体電脳倶楽部)

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

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

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

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

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

A.8.2 切断波数の決め方

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

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

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

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

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

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

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-truncation.tex)

さて, 以下では 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.88)

で与えられる. これらの格子点値から逆に Ã_{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.89)

である. ここで, 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.90)

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

23実際の 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 以上の情報は格子点値からスペクトルに変換するときに落ちている.

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

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

spectral/spectral.tex(spectral/spl-truncation.tex)  2008/06/26(地球流体電脳倶楽部)

∑_{i=1}^{I} exp(i(m′ − m)λ) = Iδ_{mm′},  (A.91)
∑_{j=1}^{J} P_{n}^{m}(sin φ_{j})P_{n′}^{m}(sin φ_{j})w_{j} = δ_{nn′}  (A.92)

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

ちなみに, 格子点値からスペクトルの係数に変換し格子点値にもどすという立場からすれば, この Gauss-Legendre の公式の適用条件というのが情報を落とさずには済まない理由である^{25} . このことを以下に述べる. 情報を落とさずに格子点値をスペクトルの係数に変換し格子点値にもどすには, あらゆる東西波数について南北方向の格子点数 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 函数は計算してはならない. つまり情報を落とさざるをえない^{26} .

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

A = BC,  (A.93)
B = ∑_{m=−M}^{M} ∑_{n=|m|}^{N} (B̃_{n}^{m} exp(imλ)) P_{n}^{m}(sin φ),  (A.94)
C = ∑_{m=−M}^{M} ∑_{n=|m|}^{N} (C̃_{n}^{m} exp(imλ)) P_{n}^{m}(sin φ)  (A.95)

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

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

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

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

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

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-truncation.tex)

Ã_{n}^{m} ≡ { \atop (BC)_{n}^{m}}

= \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} B_{ij}C_{ij}P_{n}^{m}(sin φ_{j}) exp(−imλ_{i})w_{j}

= \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}))

× ( ∑_{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}

= \frac{1}{I} ∑_{m′=−M}^{M} ∑_{n′=|m′|}^{N} ∑_{m′′=−M}^{M} ∑_{n′′=|m′′|}^{N} B̃_{n′}^{m′} C̃_{n′′}^{m′′}

× ∑_{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.96)

この計算が Ã_{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} を満たす数に選ぶ. (dcpam4 では偶数でなくてはならない. )

spectral/spectral.tex(spectral/spl-truncation.tex)  2008/06/26(地球流体電脳倶楽部)

例えば, 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}] にする. ただし [ ] はそれを越えない最大の整数を表す記号である^{28} .

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 でよい^{29} .

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

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

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

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

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

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

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

30adjustment 等の意味をなど考えると, 特に物理過程においては, 格子点で考える方が物理的に当然であるように思う.

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

31出典は, スペクトル法による数値予報（その原理と実際）（1.6）

32問題点その１. グリッドモデルでは緯度経度図で等間隔に格子点をとると, 極でも CFL を満たすようにするために, 時間差分を細かくしなければならない. 他は未調査.

2008/06/26(地球流体電脳倶楽部)  spectral/spectral.tex(spectral/spl-spectral-vs-difference.tex)

\underline{dcpam4 支配方程式系とその離散化  A.10 参考文献  65}

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

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

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

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

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

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

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

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

A.10 参考文献

気象庁予報部, 1982 : スペクトル法による数値予報（その原理と実際） . 気象庁, 111pp. 森口, 宇田川, 一松編 ,1956 : 岩波数学公式 I . 岩波書店, 318pp. 森口, 宇田川, 一松編 ,1960 : 岩波数学公式 III . 岩波書店, 310pp. 一松信, 1982 : 数値解析. 朝倉書店, 163pp. 森正武, 1984 : 数値解析法. 朝倉書店, 202pp. 寺沢寛一, 1983 : 自然科学者のための数学概論（増訂版） . 岩波書店, 711pp.

spectral/spectral.tex(spectral/spl-references.tex)  2008/06/26(地球流体電脳倶楽部)

\underline{66 dcpam4 支配方程式系とその離散化付録 B 使用上の注意とライセンス規定}

付録B 使用上の注意とライセンス規定

B.1 ライセンス規定

COPYRIGHT^{1} を参照ください.

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.

1http://www.gfd-dennou.org/library/dcpam/dcpam4/dcpam4 current/COPYRIGHT 2008/06/26(地球流体電脳倶楽部)  CREDITS.tex

\underline{dcpam4 支配方程式系とその離散化  B.3 開発グループメンバー  67}

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

B.3.1 2007 年度

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

プログラム製作協力土屋貴志, 山田由貴子, 高橋芳幸, 小高正嗣, 堀之内武, 林祥介

解説文書作成石渡正樹, 森川靖大

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