GFD セミナー特別編 2024 2024/03/16

ISPACK の開発とその大気力学研究への応用
石岡圭一 (京大・理・地惑)E-mail: ishioka@gfd-dennou.org

ISPACK 開発の契機と, 林先生との御縁のふりかえり
4 回生で余田先生の下で卒論研究をスタートし (1989 年), そこで, 球面のス
ペクトル法の勉強を, 気象庁予報部電子計算機室報告・別冊第 28 号「スペクトル法による数値予報 (その原理と実際)」(金光正郎・佐藤信夫・伊佐真
好・工藤達也, 1982) で勉強することから始めた. 当時, 山田先生も防災研にいらっしゃって, 御助言をいただいた. その後, 余田先生と山田先生が, 球
面のスペクトル法を使った減衰性 2 次元乱流の論文, Yoden and Yamada(1993), を出版される.
私自身は, 4 回生の頃から, 自分で FFT から書き始め, M1 の頃にはそれなりに使える球面のスペクトル法の変換プログラムを作って, それを用いて, 帯状流の安定性に関する修論研究, D 論研究をやり始めていた. 当時は, 京大のスパコンは富士通のベクトル計算機 VP400E (1.4GFLOPS) とか, VP2600(5GFLOPS) とかだった. 私のプログラムも, ひたすら「ベクトル長を長くする」ということだけ考えたプログラムだった.

その当時, 気象学会の力学のセッションで発表するたび, 林先生から常に建設的なコメントをいただいていた. また院生の頃に, 気象研究ノート 156
号「ロスビー波」の中の林先生の名著「二次元ロスビー波の線形論」(林,1987) を読んで, 感銘を受けていた.
その後, 京大で博士をとって (1995 年), 東大数理にいらっしゃった山田先生に学振 PD で受け入れてもらう. そこに林先生もいらっしゃって, 計算機環
境も含め, 大変お世話になり, GFD 電脳倶楽部の活動にもそれ以来, 巻き込
んでいただいた.
それまで開発してきていた球面調和関数変換のライブラリを OSS として整備・公開することを林先生に勧められ, それが ISPACK の整備につながった.

ispack-0.0 のリリースが 1997 年 4 月 15 日. また, 当時, 林先生の下に豊田英司さん (現・気象庁) がいらっしゃって, ライセンスのありかたについて高みから解説して下さっていたので, その影響があり, LGPL を採用した次第.

ということで, ISPACK のある意味, 「プロデューサー」は林先生と言っても過言ではない.
また, 東大数理には, 林さんが組織した nettomo/uttomo という手練の猛者
達 (例えば, 現・京大数学の教授の加藤周さん, 現・東大情報システム本部
講師の中村誠さん, 現・筑波大学計算科学研究センター教授の額田彰さん,
等々) が (当時は学生として) いて, そこでいろいろ教わった. 私に Debianを紹介してくれたのは, 確か中村誠さん. それで電脳サーバー達は, 私がDebian に改宗したことによって, 今でも Debian で動いている (それまでは,
私は Slackware 使っていた). ちなみに, Debian な電脳サーバーは, 1999 年に Pentium3 550MHz な CPU のマシンでスタートした (UMAX の筐体だった). ちなみに, 当時の Debian は Slink. 当時, nettomo/uttomo の人々に出会わなければ, サーバー管理とかとは無縁の人生を歩んでいたかも. Linux は使い続けていたと思うけど.

というわけで, ISPACK は, 林先生なしでは存在しえなかったということで,今日は, 林先生への感謝も込めて, ISPACK の高速化の核心となっている,

Ishioka (2018) で開発した漸化式の話や, ISPACK の応用例の一つとしての3 次元スペクトル法の大気力学コアの話 (Ishioka, et al., 2022), および, それから派生した, Pekeris 波の等価深度の見積りの話 (Ishioka, 2023; Ishizaki,

et al., 2023) などについて紹介する.





Ishioka, K., 2018:
A new recurrence formula for efﬁcient computation of sphericalharmonic transform. JMSJ, 96, 241–249.
に関する話題.

球面のスペクトル法のざっくりとしたレビュー
従属変数を球面調和関数で展開.
MX

f (λ, µ, t) =n
Xn=0 m=−nm(λ, µ).
(t)Y
am
n
n

λ: 経度, µ = sin θ, θ: 緯度, t: 時刻, M : 切断波数.
Ynm(λ, µ) = Pn|m|(µ)eimλ.
Pnm(µ): Legendre 陪関数.

Pnm(µ) ≡

v
u
u
u
u
t

(n − m)! 1(2n + 1)(n + m)! 2nn!
n+md×(1 − µ2)m/2 n+m (µ2 − 1)ndµ(0 ≤ m ≤ n).

Legendre 陪関数による展開
球面調和関数展開では, 以下の Legendre 陪関数による展開に最も計算コストがかかる.
逆変換:

gjm =

M
X
n=m

m(µ ) (m = 0, · · · , M ; j = 1, · · · , J).
sm
P
j
n n

正変換:

sm
n =

J
X
j=1

wj gjmPnm(µj ) (m = 0, · · · , M ; n = m, · · · , M ).

J: ガウス緯度の個数, µj : ガウス緯度, wj : ガウス重み.

必要な計算量

Legendre 陪関数の m ≥ 1 の成分に対しては実部・虚部を扱わねばならないことと, ルジャンドル陪関数の対称性から添字 j に関する計算は半分のJ/2 でよいことを考慮すると, 逆変換・正変換に必要な計算量は N はルジャ
ンドル陪関数をすべて事前に計算しておくとしても (乗算と加算が必要なことを考慮して)J
N = · (M + 1)2 · 2 = J(M + 1)2.
2
さらに, 正変換時の数値計算で誤差が出ないことを保証するためには3 M としておく必要があるので,
J>2

3 3
N ∼ M .
2
M が大のときは, この変換計算をいかに効率化するかが全体の計算スピー
ドを決める.

x64 な CPU 上 (というか, ベクトル機でない普通の CPU) での最適化
今日, x64 な CPU が広く普及しているので, コードをそれに合わせてチューンしておくことは有意義.

x64 な CPU はメモリバンド幅が狭いので, キャッシュの活用が重要.
しかし, Legendre 陪関数変換を素朴に実装すると, 行列 × ベクトル型の演算になるので, キャッシュに載りにくい.

→ 解決策: Legendre 陪関数を変換と同時に計算する.

ただし, Legendre 陪関数を計算する計算量は変換そのものの計算量と同じオーダー. なので, Legendre 陪関数の計算コストを下げられる工夫が可能
ならやるべき.

Legendre 陪関数は, 普通, 以下の漸化式で計算される.
m (µ) = (µP m(µ) − ϵmP m (µ))/ϵm ,
Pn+1
n
n n−1
n+1
r

2 − m2)/(4n2 − 1).
ϵm
(n
=
n
m を事前に計算しておくと
この漸化式を 1 回計算するには, ϵm
および
1/ϵ
n
n

すれば, 乗算 3 回と減算 1 回 (積和算としてカウントすると 3 回).

Ishioka (2018) では 変換の式および漸化式を工夫することによって, 漸化式部分のコストを削減した. 積和算命令にも馴染むものにした (1 段あたり
実質積和算 1 回).

着想に至った経緯
最初は, Clenshaw の漸化式みたいのが使えないかと思っていた. すなわち,
例えば逆変換の方の計算:

gjm =

M
X
n=m

m(µ )
sm
P
j
n n

は, 漸化式m (µ) = (µP m(µ) − ϵmP m (µ))/ϵm
Pn+1
n
n n−1
n+1

と組み合わせると, 以下のように計算できる (上付き添字の m は省略).

yM = sM ,yM −1 = (µj /ϵM )yM + sM −1,
yn = (µj /ϵn+1)yn+1 − (ϵn+1/ϵn+2)yn+2 + sn
(n = M − 2, M − 3, . . . , m + 1)

を計算し,

gj = ym+1Pm+1 + (sm − (ϵm+1/ϵm+2)ym+2)Pmとすれば, gj が求まる. このアルゴリズムだと, 漸化式の計算と逆変換が, n
一段あたり, 乗算 3 回, 加算 2 回で済むので, 漸化式と逆変換を素朴にやるのに比べると, 乗算が 1 回省略できている勘定にはなる. 積和算という意味でも 3 回でできるので, 積和算 1 回分は得. ただし, 展開係数の実部虚部を考えると, Pn の漸化式と逆変換を Clenshaw の漸化式でまとめても得にならない (むしろ Pn が陽に求まらないので計算量増える).
ということで, 漸化式自体の計算量を減らせないかということを考えて, IS-

PACK1 の頃は, 次にあげる手法の前身となるアルゴリズムで, 漸化式 1 段あたり乗算 2 回, 加算 1 回には減らせていた (積和算でカウントすると, 3 回にだったものが 2 回). ただ, 今一つな感じだったが, 次のページ以降のアルゴ
リズムを「発見」して効率化された.

元々の漸化式 (以下, ごちゃごちゃするので, 上付き添字 m は省略)

µPn(µ) = ϵn+1Pn+1(µ) + ϵnPn−1(µ)
を 2 段階組み合わせると,µ2Pn(µ) = ϵn+1ϵn+2Pn+2(µ) + {(ϵn+1)2 + (ϵn)2}Pn(µ)+ ϵnϵn−1Pn−2(µ)となる. 赤道反対称の Legendre 陪関数について考えて,

n = m + 2l + 1 (l = 0, 1, . . .) とすると, 上の漸化式は,µ2Pm+2l+1(µ) =ϵm+2l+2ϵm+2l+3Pm+2l+3(µ)+ {(ϵm+2l+2)2 + (ϵm+2l+1)2}Pm+2l+1(µ)+ ϵm+2l+1ϵm+2lPm+2l−1(µ),と表せる. ここで,Pm+2l+1(µ) = µαlpl(µ),

と pl を導入すれば,µ2αlpl(µ) =ϵm+2l+2ϵm+2l+3αl+1pl+1(µ)+ {(ϵm+2l+2)2 + (ϵm+2l+1)2}αlpl(µ)+ ϵm+2l+1ϵm+2lαl−1pl−1(µ),となる. この漸化式を簡単にするために, αl に対して

ϵm+2l+2ϵm+2l+3αl+1 = −ϵm+2l+1ϵm+2lαl−1,を満すことを課す. この式の両辺に αl を掛けて変形すると,
ϵm+2l+1ϵm+2l
αl+1αl = −αlαl−1,ϵm+2l+3ϵm+2l+2を得る. この式を繰り返し適用すると,

αl+1αl = (−1)lϵm+3ϵm+2
α1 α0 ,
ϵm+2l+3ϵm+2l+2を得る. 従って, あとは α0 と α1 を定めさえすれば, 順次, αl (l = 2, 3, . . .)を定めることができる.

α0 と α1 の定めかたには任意性があるが, 後々の便利のため, 以下のように定めることにする.

1

ϵm+1
α0 =
, α1 =.
ϵm+1
ϵm+2ϵm+3
こう定めると, 結局,(−1)lαl+1αl =.
ϵm+2l+3ϵm+2l+2となる. 先の pl(µ) に対する漸化式の両辺に (−1)lαl を掛けて, 上の式を用いると,pl+1(µ)= (−1)l(αl)2µ2 − (−1)l(αl)2 (ϵm+2l+2)2 + (ϵm+2l+1)2+ pl−1(µ).
となる.



pl(µ)

さらに, al と bl を

al = (−1)l(αl)2,

bl = −al (ϵm+2l+2)2 + (ϵm+2l+1)2 ,と導入すれば, 結局, pl(µ) に対する漸化式

pl+1(µ) = (alµ2 + bl)pl(µ) + pl−1(µ).
が得られる. この漸化式においては, al と bl および µ2 を前もって計算しておけば, 漸化式一段あたり, 2 回の積和算 (乗算 2 回, 加算 2 回) で計算できる. この漸化式は元の Pn で考えれば, n を 2 つ増す漸化式に相当するので,

n の一段あたり, 1 回の積和算で済んでいて, 元々の 3 回の積和算の必要性に比べて, 計算量が 1/3 になることになる.

ただし, この n について 2 つずつ増える漸化式を Pn(µ) の偶奇モード両方について計算していては全く意味が無い (むしろ以前のアルゴリズムの方が速いことになる) ので, どちらか一方についてのみ計算することになる. 実
際には, µ = 0 付近の精度を考えて, 奇モードのみについて計算する.
隅モードについては,µPn(µ) = ϵn+1Pn+1(µ) + ϵnPn−1(µ)
を用いて, 2 つの奇モードの重ね合わせ (の 1/µ) として計算できる.
このアルゴリズムを利用することの副産物として, 乗算の削減だけでなく,Legendre 陪関数に関するロードとストアが半分で済むので, その分のコスト削減にもつながる. 実際には, 例えば, 逆変換は以下のように計算される.



gj =

M −m
2
X

l=0



(tlαl)pml (µj ) + µjM −m−12
X

l=0

(sm+2l+1αl)pml (µj ).

ここに,





sm+2lϵm+2l+1 + sm+2l+2ϵm+2l+2 (l = 0, 1, . . . , [ M −m
] − 1)
2
tl = 
M −m ]).

 s
ϵ
(l
=
[
m+2l m+2l+12

その他の工夫:
論文では, m が大きい場合に生じるアンダーフローの回避の手法についても提案している. また, ISPACK3 の実装では, 変換のホットスポット部分はアセンブリ言語で書くなどもして高速化している.

他の実装との速度比較
比較対象: SHTns
作者: Nathanaël Schaeffer (ISTerre, Université de Grenoble)

URL: http://users.isterre.fr/nschaeff/SHTns/
論文: Schaeffer(2013): Efﬁcient spherical harmonic transforms aimed at

pseudospectral numerical simulations, Geochemistry, Geophysics,
Geosystems, vol. 14, pp.751–758.
開発履歴: Changelog を見ると, 2010 年に v1.0 がリリースされていたようだ. 現在は v3.6.6 (2024 年 2 月 29 日リリース) (ただし, 以下にでてくる
ベンチマーク結果は, v3.4.6 (2021 年 5 月 15 日リリース) を用いたもの)※実装の方向性は私のものに非常に近い (Legendre 陪関数を on-the-ﬂy で計算するとか, SIMD 命令の活用とか).

SHTns も, v3.2-r660 以降では, “New recurrence formula of Ishioka (2018)leading to faster transforms, especially for large sizes. see https://doi.org/10.
019” として私の手法が取り入れられているので, それとの比較も入れる.
以前は, conﬁgure 時に –enable-ishioka というスィッチをオプションを付けてから make する必要があったが, 今は (v3.4 以降) そちらがデフォルトになっている. 従って, あえて –disable-ishioka というスィッチを有効にした場合とも比較しておく.

ベンチマーク結果 (AVX2)
プラットホーム:

CPU: Broadwell Xeon E5-2699v4 (22 コア)×2
コンパイラ: ifort 19.0.0.117 -xHost -qopenmp -align array64byte (SHTns
は, gcc 6.3.0, option -O2 -march=native -ffast-math)

44 スレッドの場合のみを示す.

スペクトル⇔格子点値 1 回あたりの変換にかかる時間 (sec)(SHTns’ は, SHTns を –disable-ishioka で make した場合).
M = 1023
M = 2047
bwd
fwd
bwd
fwd
SHTns’ 0.0021 0.0022 0.014 0.015
SHTns
0.0015 0.0015 0.0097 0.0098
ISPACK3 0.0017 0.0014 0.011 0.0094
M = 4095 M = 8191
bwd fwd
bwd fwd
SHTns’ 0.11 0.11 0.86 0.86
SHTns
0.070 0.071 0.53 0.52
ISPACK3 0.075 0.064 0.53 0.45

M = 16383
bwd fwd
7.4 7.3
4.5 4.3
4.0 3.4

ベンチマーク結果その 2(AVX512)
プラットホーム:

CPU: Skylake-W Xeon W-2135 (3.7GHz, 6 コア)
コンパイラ: ifort 19.1.3.304 -xHost -qopenmp -align array64byte (SHTns
は, gcc 8.3.0, option -O2 -march=native -ffast-math)

6 スレッドの場合のみを示す.
スペクトル⇔格子点値 1 回あたりの変換にかかる時間 (sec)

M = 1023
M = 2047
bwd
fwd
bwd fwd
SHTns
0.0044 0.0044 0.033 0.033
ISPACK3 0.0048 0.0045 0.029 0.027
M = 4095 M = 8191
bwd fwd
bwd fwd
SHTns
0.20 0.19 1.6 1.5
ISPACK3 0.18 0.18 1.3 1.3

M = 16383
bwd fwd
11.0 10.9
9.9 9.6

インパクト:Ishioka (2018) で提案した漸化式は, 上で紹介した SHTns だけでなく, 他の開発者の球面調和関数変換ライブラリにも既に利用されている. cf. lib-

sharp2 by Martin Reinecke ( https://gitlab.mpcdf.mpg.de/mtr/libsharp/ ).
そこでは, “This update features signiﬁcant speedups thanks to importantalgorithmic discoveries by Keiichi Ishioka と謝辞を書いてもらっている.
ということで, GFD 電脳倶楽部の方向性としても, それなりにインパクトのある仕事ができたかなと思っている. また, SHTns は, README の中で,

“If you use Ishioka’s recurrence (the default since SHTns v3.4), you may
also want to cite his paper:” として, Ishioka (2018) の文献情報を入れてくれているので, SHTns を利用して成果を出した論文が出ると, ほぼ自動的にIshioka (2018) も引用してくれるのでありがたい. 特に, SHTns は, MagIC(https://github.com/magic-sph/magic) というダイナモコードで使われているので, ダイナモ業界の論文でも Ishioka (2018) を引用してもらっている.

Ishioka, K., N. Yamamoto, and M. Fujita, 2022:
A Formulation of a three-dimensional spectral model for the primitive
equation. JMSJ, 100, 445–469.
に関する話題.

プリミティブ方程式系に対する 3 次元スペクトル法の定式化
プリミティブ方程式系にもとづく静力学全球大気モデルは, 水平方向の離散化に球面調和関数展開のスペクトル法を用いているものが多いが, 鉛直
方向は主に差分法が使われている.
鉛直方向にもスペクトル法を使う, 3 次元スペクトル法のモデルの定式化について, 山本直人君 (2016 年度修士論文), 藤田雅人君 (2019 年度修士論文)で取り組んでもらったが, それらの結果と, 私が追加した計算結果や, 傾圧
トイモデルの定式化などの内容も含めて, Ishioka, et al. (2022) として論文にまとめたので, それに沿って話をしていく.

開発の動機: 修士課程以降, ISPACK(およびその前身) を使って, 球面上の 2次元流体 (順圧/浅水) の不安定や乱流の研究を自分自身, および学生さん達といろいろやってきていたが, さすがに 3 次元の話がやりたくなっていた.
修論研究の学生さん (永井君) の修論 (2018 年度) として, 独自に球殻対流の3 次元スペクトル法の実装 (動径方向に普通と違って Legendre 多項式を用いる) はやってみてもらったことはあったが, 地球大気を念頭に置いた, プ
リミティブ方程式系のものを独自にちゃんと実装したくなった.
球面上のプリミティブ方程式系の実装として, 水平方向にスペクトル法, 鉛
直方向には差分法というものが世の中のスタンダードで, DCPAM も含めて, ソースが公開されているものも昔から多数存在している.
ということで, 折角新しいものをスクラッチから作るなら, 普通と違うものを作ってみたかった. そこで, プリミティブ方程式系で, 鉛直方向にもスペクトル法を使う実装ができないかということでスタートした.

このような試みは, 沢山ありそうな気がするのだが, 調べてみると, 意外とあまり無くて, 鉛直方向に Legendre 多項式展開を用いた Machenhauerand Daley (1974) と, Chebyshev 多項式展開を用いた黒木・村上 (2015) く
らいしか見つからなかった.
どちらの先行研究も σ = p/ps 座標のプリミティブ方程式系を対象としているが (p は圧力, ps は地表面圧力), Machenhauer and Daley (1974) では

σ → 0 での特異性 (後述) の扱いは提案はあるもののアドホックであり, また, 黒木・村上 (2015) ではその考察が無い. また, 時代的に, 当然ながら現代的なベンチマーク計算などはやらておらず, 水平切断波数 9(平行四辺形
切断), 鉛直の切断の次数 5 の計算のみ. 黒木・村上 (2015) では, 現代的なベンチマーク計算も行われているが, 特異性の回避などに関する考察は無い.

何が問題?
支配方程式の中に現れる, 静力学平衡の式を鉛直積分した形の式:′(λ, µ, σ ′, t)
τ
σ
′.
dσ
Φ′ = Φ′s − 1σ′Zにおいて, τ ′ を普通に σ 方向に Legendre 多項式展開していると, σ → 0 で
右辺が対数発散する. さらに, 右辺の積分結果に対数関数を含むようになるため, 変換法を使った場合にエイリアスエラーを除去できなくなる. ここで, λ は経度, µ はサイン緯度, Φ′ はジオポテンシャルの擾乱成分, Φ′s は地
表のジオポテンシャルの擾乱成分, τ ′ は温度場の擾乱成分. どちらも無次元化されている.
問題解決へのアイディア:

τ ′(λ, µ, σ ′, t) を, σ × (σ の多項式) の形で展開する. また, 多項式展開の部分は, Legendre 多項式展開を用い, 展開係数の時間発展の定式化には, ガ
ラーキン法を用いる. すなわち, 残差にかける重みとして, σ × Pn(η) を用
いる. ここに, η = 1 − 2σ で, Pn(η) は n 次の Legendre 多項式.

また, 発散場・渦度場は, 単に Pn(η) で展開する. このような定式化によって, τ ′ の時間発展部分だけは余分に 5 重対角行列を係数とする連立一次方程式を解く手間が増えるが, コスト的には大したことはない. ちゃんとガラー
キン法で定式化することによって, Semi-implicit 法を適用する際に現れる行列が対称行列になるなど, 諸々, とても美しく定式化される.

Held and Suarez (1994) のベンチマーク設定の計算の例 (水平切断波数 T85,
鉛直方向切断次数 13, 鉛直方向グリッド数 20) (Ishioka, et al. (2022) Fig. 2)1000 日間平均の帯状平均温度場 (上パネル) および風速場 (下パネル).

Jablonowski and Williamson (2006) のテストケース設定 (ただし, 初期擾
乱を 1/1000 にしている) の計算の例 (水平切断波数 T85, 鉛直方向切断次数170, 鉛直方向グリッド数 256) (Ishioka, et al. (2022) Fig. 5)

時間発展の 19 日目の地表面気圧場.

Jablonowski and Williamson (2006) のテストケース設定 (ただし, 初期擾
乱を 1/1000 にしている) の計算における, 鉛直方向切断次数 170 の結果を基準とした誤差の鉛直方向切断次数依存性 (水平切断波数 T85) (Ishioka, etal. (2022) Fig. 6)

縦軸が地表面気圧の誤差の l2 ノルム. 横軸が鉛直方向切断次数.

ということで, Held and Suarez (1994) のテストケースでは, 鉛直グリッド
数は 20 だが, 鉛直の自由度は 14 であるので, 差分法より少ない自由度で精度良く計算できていることが分かる. その他, Polvani, et al. (2004) のテス
トケースでもチェックを行っている.

Jablonowski and Williamson (2006) のテストケース設定 (ただし, 初期擾
乱を 1/1000 にしている) では, 鉛直の切断次数を上げていくと, 傾圧不安定
擾乱が発達した段階での地表面気圧の気圧の誤差が, 急速に小さくなっていくことが確認できた.
その他, 固有値解析や, スポンジ層の設定, 鉛直自由度を落したトイモデルの構築など, いろいろ副産物もある. その一例として, Lamb 波の計算も行っている.

等温静止大気で回転が無い場合を考えると, 3 次元スペクトルモデルを線形化したものから得られる一般化固有値問題の式は, 以下のようになる. ここに, c は無次元化された位相速度, A は三重対角行列, B は五重対角対称行列.
 √κδ̂k,0 
0









.
.




.
I
O
. 



 √




κ
δ̂



0
k,L






′




kτ̂
− ic 
0  
k,0 



.
.




B
O
. 
.








′

  kτ̂

0



k,L−1



√
0 ··· 0 0 ··· 0 1
k κŝk


 √
(1)
κ
δ̂
1


k,0

 √







κ
δ̂


0
k,1




√



.
.


.

O
κA
.




 √






κ
δ̂
0


k,L


.

= 

′


kτ̂


0
k,0 


√



T
.


.

−
κA
O
.
.








′

0   kτ̂k,L−1 



√

−1 0
···
0 0 ··· 0 0
k κŝk3 次元スペクトル法の定式化に基いて計算した Lamb 波の鉛直構造 (Ishioka,
et al. (2022) Fig. 9)

横軸は振幅, 縦軸は σ. 点線は理論解 σ −κ. 実線が数値解. また, κ = 2/7.

√
計算された Lamb モードの ( RT0 で無次元化された) 位相速度の鉛直切断q

次数 L 依存性. 真値は 7/5 ≈ 1.183216.
L = 10 L = 20 L = 40 L = 80
1.170342 1.176177 1.179378 1.181121

3 次元スペクトル法の定式化で計算した Lamb 波の位相速度の誤差の切断次数依存性 (Ishioka, et al. (2022) Fig. 10)Ishioka, K., 2023:
What is the equivalent depth of the Pekeris mode?
JMSJ, 101, 139–148.
およびIshizaki, H., T. Sakazaki, and K. Ishioka, 2023:
Estimation of the equivalent depth of the Pekeris mode using reanalysis
data. JMSJ, 101, 461-469.
に関する話題.

Ishioka, et al. (2022) が accept されたのが 2021 年 12 月 17 日. その約 1ヶ月
後, 2022 年 1 月 15 日に Hunga Tonga-Hunga Ha’apai 火山の噴火があった.

The HTHH eruption as viewed from the Himawari-8 weather satellite(Proud, et al. (2022) Fig. 1)
04:20UTC

04:50UTC

05:40UTCその後, この噴火で, Lamb 波だけでなく, Pekeris 波が検出されたことを示す論文 (Watanabe et al., 2022) が投稿されたことを知る (2022 年 4 月).

The 10-minute difference of the 9.6µm
brightness
temperatureobserved
by Himawari-8 (08:40UT−08:30UT)(Watanabe, et al. (2022) Fig. 1)

Watanabe, et al. (2022) において, 擾乱伝播の位相速度の方位角平均から,Lamb 波の等価深度は 10.1km, Pekeris 波の等価深度は 6.1km と見積られ,それらの値は, Salby (1979) が, 大気自由振動の鉛直構造方程式を U.S.

STANDARD ATMOSPHERE, 1976 (以下, USSA76 として引用) の大気温度構造に対して計算して得られた値の 9.6 km と 5.8 km と対比されていた.
これに興味を持ち, Ishioka, et al. (2022) の定式化で, USSA76 の大気温度
構造を与えた場合, 等価深度がいくらになるのか計算してみた.

3 次元スペクトルモデルでの固有値計算の結果 (等価深度 (km))
鉛直格子点数 鉛直切断波数 Lamb モード Pekeris モード
1024
1023
9.90
6.32
1024
682
9.90
6.30
512
341
9.90
6.24
256
170
9.89
6.19
128
85
9.88
6.03

ということで, 鉛直の切断次数をかなり大きくしても, Salby (1979) の結果と合わない.
何か Salby (1979) が考慮していてこちらで見落していることもあるのかと思い, Salby (1979) がやったことを徹底的に見直してみた. 結果, Salby(1979) の計算に問題があったことが分かり, それを指摘してちゃんとした精度で等価深度を計算しなおした論文 (Ishioka, 2023) を出版した.
以下, 背景となる大気の自由振動の話からおさらい.

大気の自由振動 (その 1)
プリミティブ方程式のもと, 静止大気を基本場とした場合, 擾乱に関して線形化したされた大気運動の基礎方程式は, 水平方向と鉛直方向に変数分離
され, 水平方向は全球では Laplace の潮汐方程式の形になり, 「等価深度」

h の深さの線形浅水方程式と同じになる. このとき, 自由振動の等価深度は,
鉛直構造方程式の固有値として定まり, 鉛直温度プロファイルによって値が決まる.

大気の自由振動 (その 2)
等温大気のとき, その場合のスケールハイトを H とすると, h は唯一しか解がなく = γH となる (γ は比熱比) (Lamb, 1911)
現実大気の等価深度がどうなるかは, 大気全体の鉛直温度構造が分からな
いと鉛直構造方程式が解けないため進展がなかったが, Taylor (1929) が 1883
年の Krakatoa の噴火の際の気圧擾乱の伝播速度から, 等価深度が 10.4km程度であると見積った. (現在では, この等価深度のモードは Lamb モードと呼ばれている).

大気の自由振動 (その 3)
この 約 10km というのが大気の自由振動の唯一の等価深度と思われたが,熱帯域の大気の熱潮汐の振幅が 1 日周期より半日周期の方が大きいことの原因として, 熱強制と自由振動との共鳴が疑われていた (これは違うということが現在では判明しているが). 8km 程度の等価深度のモードが存在するかどうかが Kelvin (1882) 以来, 大きな関心事であった. そのような時代背景のもと, Pekeris (1937) が, 成層圏界面の高温域の温度によっては, Lamb
モードと別個に, 8km 程度の等価深度のモードも存在しうることを示した.

Pekeris (1937) で使われた鉛直温度プロファイルと, 計算された等価深度(Pekeris (1937) の Fig. 2 と Table 1)

大気の自由振動 (その 4)
より現実的な鉛直温度分布を与えた等価深度の計算は Salby (1979) によって USSA76 に対して行われたが, そこでは, Lamb モードに相当する等価深度として 9.6km が得られただけでなく, 5.8km の等価深度も得られた. 後者の方が Pekeris が存在可能性を指摘したモードに対応するものだが (Salbyは ducted mode と呼んだ), そのようなモードの実在性についてはこれまで検出されていなかった.

大気の自由振動 (その 5)

Watanabe et al. (2022) では, Himawari-8 の輝度温度データの時間差分をとることにより, 位相速度 315m/s 程度の波 (等価深度 10.1km に対応) と,245m/s 程度 (等価深度 6.1km に対応) の波が検出された.
前者は Lamb モードに相当するものだが, 後者は Pekeris が存在可能性を指摘したモードに相当するものと考えられるため, これを Pekeris モードと呼ぶことを提案した. このように Pekeris モードが検出されたのはこの論文が初めてのことである. この論文では, GCM を用いた数値実験や, ERA5
再解析データでも Pekeris モードが見出されることを示している.

Power spectrum for eastward-propagating wavenumber 1 and 2 components for equatorially symmetric components of the surface pressure
calculated with ERA5 data for 20◦S–20◦N during 1950–2016 (Watanabe, et al. (2022) Fig. 9d)Gray solid (denoted by “P”) and dashed lines (denoted by “L”) are for resonant frequencies for h=6.1 km
and h=10.1 km, respectively.

大気の自由振動 (その 6)以上の背景のもと, 私も Ishioka, et al. (2022) で提案した鉛直スペクトル法の GCM のコードを用いて, USSA76 の温度分布に関して, 等価深度の計算をしてみたが, 既に述べたように, Salby (1979) が計算したような 5.8km程度というほどの等価深度にはならなかった. さらに, Lamb モードの等価
深度も 9.6km よりは明らかに大きな値になっていた.
そこで, Salby (1979) の計算手法を見直してみたところ, 計算手法に問題があって大幅な精度低下を引き起していることが分かったので, その問題点と, それを修正した場合どうなるかについて報告する.

Salby (1979) はどうやって等価深度を計算したのか?
鉛直構造方程式と下端境界条件は, 線形化されたプリミティブ方程式系から以下のように導かれる:






1
d
1
d
−  ′(H̃Z) + αZ = 0,dζ H̃ (H̃ + κ) dζ
H̃Z ′ − κZ = 0 (ζ = 0).



(2)
(3)

Z(ζ) は圧力擾乱の鉛直構造を与える関数, ζ = z/H, z は幾何高度 H は定められたスケールハイト, κ = (γ − 1)/γ で, γ は比熱比 (7/5), α = H/h
で, h は等価深度.

R0T (ζ)H̃(ζ) =
g0 H

(局所的なスケールハイトを H で無次元化したもの).
(4)

R0: 乾燥大気の気体定数, g0: 重力加速度, T (ζ): 鉛直温度プロファイル.

Salby (1979) は, 鉛直構造方程式 (2) と下端境界条件 (3) をそのままの形では扱わず, 以下のように変換してから使っていた.

Z(ζ) = eξ/2H̃ −1[H̃ ′ + κ]1/2v(ζ)(5)

ここで, ξ は以下のように定義される.
Z

ζ dη
.
ξ(ζ) = 0
H̃(η)

(6)

(2) と (3) は, 以下のように書ける.
v ′′ + k2(ζ; α)v = 0,


′′
H̃ H̃
1

′
′
v = 0H̃v +  − (H̃ + κ) +2
2(H̃ ′ + κ)(7)(ζ = 0).

(8)

ここに, k2(ζ; α) は「屈折率」で, 以下のように定義される.
k2 = −

1
H̃ ′′′
3(H̃ ′′)2H̃ ′′α(H̃ ′ + κ)
+
−
−
+
.
2
′
′
2
′
4H̃
2(H̃ + κ) 4(H̃ + κ)
2H̃(H̃ + κ)
H̃(9)USSA76 の温度プロファイルから計算される H̃ と屈折率 k2 の分布(Salby (1979) Fig. 1 and 3)

ここに, H = R0T∗/g0 で, T∗ = 250 K としているので, H ≈ 7.32 km.

Salby (1979) は (7) を十分な高高度で放射境界条件を定めて, そこからから下向きに数値的に解いていき, 下端境界条件 (8) が (近似的に) 満たされる αの値を探した.

Salby (1979) Fig. 2 (α = H/h)

2 つの明瞭な凹みが α = 0.764 と α = 1.25 に見える. これらは, それぞれ, 等価深度 9.6 km と 5.8 km に対応する. 前者は Lamb モードに対応し,後者は, Pekeris モードに対応する. (Salby (1979) は「ducted mode」と呼んでいるが).

Salby (1979) の手法の問題点
USSA76 のように, 鉛直温度勾配に不連続な点があるような分布を用いる
場合, H̃ ′′ は不連続点で δ 関数的になり, H̃ ′′′ は δ 関数の微分になってしまう.
これは, 「屈折率」k2 を定める際に問題となる.

USSA76 の温度分布 (Ishioka, 2023 のFig. 1b)

Salby (1979) の手法の問題点を克服するには?
以下のような単純な変換だけを考える.

Z(ζ) = eξ/2Z̃(ζ).

(10)

この Z̃ を使うと, 鉛直構造方程式 (2) と下端境界条件 (3) は, 以下のように書
ける.










d
1  1
1 
 d



 (H̃ Z̃) + αZ̃ = 0,−+dζ 2H̃ H̃ ′ + κ dζ2H̃H̃  d1  (H̃ Z̃) = 0
H̃ Z̃ − ′+(ζ = 0).
H̃ + κ dζ
2H̃



(11)
(12)

式 (11) を十分高高度で放射境界条件を与えて下向きに数値計算していき,
下端境界条件 (12) が (近似的に) 満される α を探す.

下端境界条件の誤差の α 依存性 (Ishioka (2023) Fig. 2)

明確な凹みが α = 0.739 と α = 1.107 で見られる. これらは, それぞれ,
等価深度

h = 9.90 km と h = 6.61 km に対応する.

すなわち, Salby (1979) の計算結果は, Pekeris モードについて特に不正確.

より凝った計算
以上のように, Salby (1979) の結果は不正確であるというのが結論だが, Salby(1979) がもしかすると論文に書かれていない細かいことを考慮してその結果を出してたのかとも思い, 念のため, 重力加速度, 平均分子量の高度依存性まで考慮した計算もやってみた.
温度

重力加速度

(Ishioka (2023) の Fig. 1 と Fig. 4a, 4b)

平均分子量このように凝ったことをしても, 結局誤差の α 依存性は殆んど変わらず, 求
まる等価深度は, Lamb モードが h = 9.90 km, Pekeris モードが h = 6.63

km であった.
結局, 有効数字 2 桁の範囲では, USSA76 の温度分布を用いる場合, Lamb
モードの等価深度は, h = 9.9 km, Pekeris モードの等価深度は, h = 6.6

km というのが結論.
※ Watanabe, et al. (2023) の結果と比較すると, Lamb モードの等価深度が約 10 km というのは整合的だが, Pekeris モードの等価深度が 6.6 km とい
うのは, 6.1 km という見積りに比べて大きい. ただし, Ishioka (2023) では,

USSA76 の温度分布を使っているので, Watanabe, et al. (2023) で Pekeris
モードを検出した緯度帯や季節等を考慮すると, 温度構造が違うので, 等価深度は違っている可能性もある.

ということで, ERA5 のデータを使って, HTHH 噴火時の熱帯域の大気の鉛直温度プロファイルを用いて, 等価深度の計算をやってみた (Ishizaki, et al.,2023).
結論として, そのように現実的な温度分布を用いても, Lamb モードの等価深度は, h = 10.1 km, Pekeris モードの等価深度は, h = 6.5 km という値であり, USSA76 を用いた場合と大きくは違わない.

さらに, ERA5 の 67 年分の熱帯の地表面気圧データのスペクトル解析から,

Pekris モードに相当するピークには, 等価深度 6.1 km よりも 6.5 km とした場合の方がより合致しているように見える (Ishizaki, et al. (2023) Fig. 5).

この部分のまとめと考察

Pekeris モードの等価深度は, USSA76 の鉛直温度分布では 6.6 km, HTHH
噴火時の熱帯の温度分布を考えても 6.5 km となる. この 6.5 km というのは, ERA5 の 67 年分の熱帯の地表面気圧データのスペクトル解析とは整合的だが, Watanabe, et al. (2022) の見積もりの 6.1 km とはずれている.
このずれの要因として, Watanabe, et al. (2022) で検出された遅い方の波動はまだ完全に Pekeris モードだけが卓越しているわけではない transient
な状態を見ているという可能性や, また, 平均流が影響している可能性もある.

全体のまとめ
以上, ISPACK の開発の経緯とそれに関連する最近の私の研究として, Ish-ioka (2018), Ishioka, et al. (2022), Ishioka (2023), および Ishizaki, etal. (2023) について紹介した.
ISPACK も含めて, 力学コアを自前で作るとか, いろいろ自分たちでやっていると, 温故知新な感じでいろいろ新たな発見があることに改めて気づかされる (「車輪の再発明」と言われかねないが).
自分自身でいろいろコードを書いてきたことについては, GFD 電脳倶楽部の先達の皆さんから学ばせてもらったことが非常に大きいので, その活動を常にリードされてきた林先生に改めて感謝したい.

参考文献

Held, I. M., and M. J. Suarez, 1994: A proposal for the intercomparison of the dynamical cores of atmospheric general circulation models.
Bull. Amer. Meteor. Soc., 75, 1825–1830.
Ishioka, K., 2018: A new recurrence formula for efﬁcient computation of
spherical harmonic transform. J. Meteor. Soc. Japan, 96, 241–249.
Ishioka, K., 2023: What is the equivalent depth of the Pekeris mode? J.
Meteor. Soc. Japan, 101, 139–148.
Ishioka, K., N. Yamamoto, and M. Fujita, 2022: A formulation of a threedimensional spectral model for the primitive equation. J. Meteor. Soc.
Japan, 100, 445–469.
Ishizaki, H., T. Sakazaki, and K. Ishioka, 2023: Estimation of the equivalent depth of the Pekeris mode using reanalysis data. J. Meteor. Soc.
Japan, 101, 461-469.

Jablonowski, C., and D. L. Williamson, 2006: A baroclinic instability test
case for atmospheric model dynamical cores. Quart. J. Roy. Meteor.
Soc., 132, 2943–2975.
Kelvin, Lord (W. Thomson), 1882: 2. On the thermodynamic acceleration of the Earth’s rotation. Proc. Roy. Soc. Edinburgh, 11, 396–405.
Lamb, H., 1911: On atmospheric oscillations. Proc. Roy. Soc. Lond.,
A84, 551–572.
Machenhauer, B., and R. Daley, 1974: Hemispheric spectral model.
GARP Pub. Ser., No. 14, 226–251.
Pekeris, C. L., 1937: Atmospheric oscillations. Proc. Roy. Soc. Lond.,
A158, 650–671.
Polvani, L. M., R. K. Scott, and S. J. Thomas, 2004: Numerically converged solutions of the global primitive equations for testing the dynamical core of atmospheric GCMs. Mon. Wea. Rev., 132, 2539–
2552.

Proud, S. R., A. T. Prata, and S. Schmauß, 2022: The January 2022
eruption of Hunga Tonga-Hunga Ha’apai volcano reached the mesosphere. Science, 378, 554–557.
Salby, M. L., 1979: On the solution of the homogeneous vertical structure problem for long-period oscillations. J. Atmos. Sci., 36, 2350–
2359.
Schaeffer, N., 2013: Efﬁcient spherical harmonic transforms aimed at

pseudospectral numerical simulations, Geochemistry, Geophysics, Geosy
tems, 14, 751–758.
Taylor, G. I., 1929: Waves and tides in the atmosphere. Proc. Roy. Soc.
Lond., A126, 169–183.
Watanabe, S., K. Hamilton, T. Sakazaki, and M. Nakano, 2022: First
detection of the Pekeris internal global atmospheric resonance: evidence from the 2022 Tonga eruption and from global reanalysis data.
Journal of the Atmospheric Sciences, 79, 3027–3043.

Yoden, S. and M. Yamada, 1993: A numerical experiment on
two-dimensional decaying turbulence on a rotating sphere. J. Atmos.
Sci., 50, 631–644.
気象庁, 1982: 気象庁予報部電子計算機室報告・別冊第 28 号「スペクトル法による数値予報 (その原理と実際)」, (執筆者: 金光正郎・佐藤信夫・伊佐真好・工藤達也)
黒木志洸・村上茂教, 2015: 鉛直方向に Chebyshev 多項式展開を用いた

3 次元スペクトルモデルの構築. 日本気象学会 2015 年春季大会予稿集,
P434.
林祥介, 1987: 2 次元定常ロスビー波の線形論. 気象研究ノート, 156, 235–

254.