- 著者：ゆみや / stepney141
- 2016年3月初版、2020年5月改訂
- この文章についての情報：https://stepney141.hatenablog.com/entry/2020/05/27/071742
- 序論
- 本論
2-1. 天気予報の歴史
2-2. 数値予報の原理
2-3. 実際のプログラムの開発 - 結論
3-1. 結果
3-2. 今後の課題・展望 - 参考文献
関数電卓とは、主に自然科学の計算に用いられる電卓のことで、高性能なものは昔のスーパーコンピュータの性能を超えるものもある。 一方、数値予報とは、大気の動きを表す方程式を解いて将来の天気を予測する方法で、世界初のコンピュータENIACで初めて実行された。 関数電卓の一機種には、ENIACの実に1万倍もの性能を持つものがあることから、僕は「関数電卓でも数値予報の実行が可能なのではないか？」との仮説を立てた。 そして、実際に関数電卓上で動作するプログラムを作成・実行し、一定条件下では関数電卓による数値予報は可能であるとの結論に至った。 しかし、計算に時間がかかりすぎてしまうなど、多くの問題が未だ存在している。
- ※2-1節・2-2節は、中学三年生当時の筆者が不十分な理解で書き上げたものであるため、記述に誤りがある可能性が高いです。
- ※2-1節・2-2節は、卒業論文の提出時に文書の構成上必要で付け加えた導入部分であるため、本筋である関数電卓用プログラムとは関係が薄い内容です。数値予報に関する知識が既にあり、「関数電卓を使った天気予報」の概要のみをご覧になりたい方は、これら2つの節は読み飛ばして下さい。
「関数電卓」とは、科学や工学、数学に代表される自然科学に関する計算を行うための電卓の一種である。 このような種類の電卓としては、1972年にアメリカのコンピュータ機器メーカーであるHewlett-Packard社が発売したHP-35が世界初であった。 それから50年以上が経過した現在でも、日本のカシオ計算機などによって製造・販売が行われている。 関数電卓は、自然科学で登場する数式を計算しなければならないため、一般的な四則計算専用の電卓よりも格段に高機能であり、高性能である。 一例を挙げると、2015年現在カシオが販売している関数電卓の一機種であるfx-375ESの場合、394個の機能・関数を搭載している。 関数電卓の中にも上位機種が存在し、中には独自のプログラミング言語を持ち、パソコンのようにユーザーがプログラムを書いて複雑な計算でもボタン1個で行えるようなものも存在する。 このような電卓は、一般に特に区別して「プログラム電卓」と呼ばれる。 プログラム電卓の中には、1940年代から60年代のスーパーコンピュータとほぼ同程度、あるいはそれ以上の性能を持つ機種も存在する。 一例を挙げると、前述したHewlett-Packard社が2015年現在販売しているプログラム電卓HP50g Graphing Calculatorの場合、搭載している機能・関数は2300個を超え、そのプロセッサには携帯電話やゲーム機などにも使われるARM920Tが使用されており、これは世界初の汎用コンピュータであるENIACと比較しても遜色のない性能を持つ。
「数値予報」とは、大気の動きを支配する運動方程式を数値的に解いて、将来の気圧や時間変化を求める手法のことである。 この運動方程式に観測された現在の気圧や気温の値を初期値として代入し、数値的に解けば、将来の気圧・気温の値を予測することができるのである。 だが、大気の運動は非常に複雑であるため、それを表す運動方程式も必然的に複雑なものとなり、数値予報の実行には複雑な方程式でも迅速に処理することができるスーパーコンピュータが必要である。 しかし、大気の動きは、様々な仮定を施してやることで、正確な大気の動きは表現できなくとも、近似した、単純化した大気の動きを簡単化した運動方程式で表現することができるため、比較的性能の低いコンピュータでも簡単なモデルであれば数値予報を実行することが可能である。 実際に、1950年、アメリカの気象学者であるチャーニーらの研究グループが、「順圧準地衡風モデル」と呼ばれる簡単化した方程式を用いることで、現在のパソコンの100万分の1にも満たない性能である上述のENIACにて、世界で初めての数値予報の実行に成功している。
ここで、先ほどプログラム電卓の例として挙げたHP50gと、世界初の汎用コンピュータで世界初の数値予報にも用いられたENIACの性能の比較をしてみたい。 比較には、1秒間に何百万回命令を実行できるかを表すMIPS（Million Instructions Per Second）値を用いる。 まず、ENIACは毎秒5,000回の演算が可能であったとされる。この数字に従うと、ENIACは0.005 MIPSとなる。次に、HP50gのMIPS値であるが、HP50gが搭載しているプロセッサのARM920Tの1秒間あたりの演算実行回数については情報を見つけることができなかったため、以下の式で計算する。 (MIPS値) = (クロック周波数) / (CPI値) ここで、クロック周波数とは、デジタル回路を同期するために発生させる信号の毎秒の周波数を指し、CPIとは、1つの命令を実行するのに必要な平均のクロック数を指す。 HP50gのARM920Tの公称クロック周波数、CPI値はそれぞれ75 MHz、1.5である。これらの値を用いると 75 * 10^6 / 1.5 / 10^6 = 50 となり、HP50gは50 MIPSであると推測される。 つまり、MIPS値のみのごく単純な比較ではあるが、HP50gはENIACの実に10,000倍もの性能を持っているということになる。 僕は、「プログラム電卓にはENIACを上回る性能のものがある」、そして「世界初の数値予報はENIACで行われた」ということから、「プログラム電卓でも数値予報の実行は可能ではないか？」との仮説を立てた。 そこで、実際の数値予報の技術や歴史について調べ、簡単化した数値予報モデルの関数電卓用のコードを実際に作成・実行することにした。また、使用するプログラム電卓は前述のHP50g Graphing Calculatorに決定した。 HP50gのスペックを以下に示す。
- CPU：Samsung S3C2410A (ARM920Tプロセッサ) 75MHz
- RAM：512 KiB
- フラッシュメモリ：2 MiB
- 画面サイズ：131×80 画素
- 端子：Mini-USBポート、IrDA、3.3V TTLレベル非同期シリアルポート
最初に、数値予報がどのようにして誕生したのかを理解するため、天気予報の歴史をまとめることにする。 人類にとって、日々の天気の移り変わりを把握することはとても重要であったため、気象の研究は古代ギリシャの時代から行われてきた。だが、科学的な気象観測が行われるようになったのはようやくルネサンス時代になってからのことである。17世紀以降、気圧計や温度計が開発され、それに伴い気象学も発達し始めた。 天気予報が成立したのは19世紀頃のことである。1820年にドイツのブランデスにより初めて天気図が作られ、雨や曇りの場所が低気圧になっていることが見出された。また、1840年代にはヨーロッパ各国に電信が広まり、気象観測データを迅速に集めることが可能になった。さらに1855年、パリ天文台長のルヴェリエが気象観測網によるデータの収集、天気図解析による天気予報、そしてそれらを行う国家による気象機関の設立を進言したことで、日本を含む世界各国で気象機関が設立された（日本の東京気象台は1875年成立、後に気象庁となる）。 だが、この頃の天気予報は未だ「主観的」、「経験的」と言えるものであった。空の観察、そして予報官の長年の経験とカンがものを言う、いわば職人技であったということである。
そのような主観的予報が一般的な中、20世紀前半から、流体力学の方程式を使って数学的に天気予報を行う「客観的予報」、つまり数値予報の研究が盛んになった。1901年にアッベが数値予報の原理を考案し、その後1904年にはV. ビヤークネスによって大気の運動方程式が導かれた。 その後、イギリスの気象学者リチャードソンによって手計算での数値予報が試みられた。彼が用いた方程式は仮定や近似を施していないものであり、ましてコンピュータの無かった当時では非現実的な試みであった。彼は数か月もの時間をかけて予測計算を行ったものの、結果得られた値は実際の観測値の100倍以上にまで大きくなってしまい、試みは失敗に終わった。 この結果は1922年に出版され、その中でリチャードソンは次のように述べている。 「おそらく漠とした将来のいつの日か、実際の天気の進みよりも、計算の方が先んじることができるだろう。しかしそれは夢である」 この言葉から、現実的な速度・精度で数値予報を行うという彼の構想は「リチャードソンの夢」と呼ばれた。 その後、ラジオゾンデの発明によって高層気象観測が行われるようになり、低気圧や高気圧の研究が進んだ。そして第二次世界大戦後の1946年、フォン・ノイマンによって汎用コンピュータが開発され、これがチャーニーらによる数値予報の成功につながった。 1950年、チャーニーらは「500hPa面」と呼ばれる位置の大気中層における数値予報をENIACで行ったが、計算時間が非常にかかってしまい、24時間先の予報をするのに24時間以上かかってしまうという結果だった。現実的に考えて、これはもはや「予報」として成立していないが、リチャードソンによる計算よりも明らかに高速であり、当時の気象学者たちは数値予報の実現が現実味を帯びてきたことに非常に驚愕したという。 なお、コンピュータの技術は1950年代以降急速に発達し、1952年にチャーニーらは別のコンピュータを用い、24時間予報を1時間程度で処理することに成功している。
数値予報とは、上述したように、大気の動きを表す運動方程式を解いて将来の大気の状態を予測する方法である。 天気は大気の動きによって変化し、大気の動きは物理現象なので、物理学の法則に従って動く。よって、物理学の方程式でその動きを予測できるのである。 大気中では、高気圧や低気圧のような数千キロメートル規模の巨大なスケールの現象から、雷などの数キロメートル規模の小さなスケールの現象まで、様々な現象が起こっている。 このように様々な現象が起こっているがゆえに、大気中で起こっているすべての現象を予報しようとすると、予報方程式はとてつもなく複雑になってしまう。 だが、その複雑な方程式の中でも、天気の移り変わりにあまり大きな影響を与えない小さな規模の現象を省略しても、大気全体の動きの予測にはあまり影響はないはずである。 チャーニーたちは、このような考えに基づいて複雑な大気の支配方程式を簡略化し、それによって当時のあまり性能の良くないコンピュータでも数値予報の実行を可能にした。
ところで、北半球においては、大気の中層に規則正しく西から東へ進む波があり、地上の低気圧や高気圧はこれに引きずられる形で東に進んでゆく。 数値予報では、この波の動きを予報して、そこから低気圧や高気圧の動きを二次的に予測するという手法をとっている。 ここで「渦度」という言葉が登場する。「渦度」とは、空気の流れが回転している有様を表現する量であり、大気の中層では、空気は別の場所に移動しても渦度は変わらない、つまり空気を回転させる力の度合いは変化しない、という法則が成り立つ。これを「渦度保存の法則」と呼ぶ。 この渦度は、大気中の気圧から求めることができ、そしてこれらの渦度と気圧が分かれば、将来の時点での、渦度を持った回転する空気を流す力の場を求めることができる。 この求めた風の場から将来の渦度と気圧を求め、そこからさらに将来の風の場を求め…ということを続ければ、12時間や24時間先の風の場を求めることができ、そこから低気圧や高気圧の動きを予測する…というプロセスを踏んで気象予測をすることが可能になる。このような流れでチャーニーらは数値予報を実行したのである。
また、このプロセスに基づいて渦度の計算をする際、コンピュータで数値を取り扱いやすいように、計算する領域を数十から数百キロメートルほどの間隔の格子で覆ってやる必要がある。チャーニーたちは格子間隔を約300キロメートルほどに設定して計算を行ったが、これは数十キロメートルほどの比較的中規模なスケール以下の現象を全て無視することを意味する。 この格子間隔は、短くすればするほど小さなスケールの現象も計算でき予測精度も上がるが、その分計算量も増加していくので、必要な精度で打ち止めになるよう適当な長さの格子で覆ってやる必要がある。 例を挙げると、300キロメートルほどの格子間隔でアジア全域を覆った場合、格子点の数は600個ほどになるが、実際に予報する際、この格子点1つ1つでそれぞれの渦度変化を計算し、将来の風の場を求め…という計算をしなくてはならない。 数値予報モデルの精度やコンピュータの性能が圧倒的に進歩している現在は、格子間隔に数キロメートルほどの値を取れるようになっているが、この格子点の問題は、関数電卓という限られた計算資源で数値予報を行うにあたっては非常に重要な問題であると言える。
ここでは、筆者が実際に作成したプログラムの概要を示す。 関数電卓は、単純なマシンスペック面では60年前のコンピュータを遥かに凌ぐとはいえ、実際には入力デバイス・入出力方法・電源・表示用ディスプレイ等の多数の面において、単純比較が不能な大きな相違点がある。また、関数電卓と、パソコンやスマートフォン等現代で一般的なコンピュータとを比較しても、同様に多くの相違点が存在する。このため、プログラムの開発にあたっては既存の資料を直接参考にすることがほとんど出来ず、様々な困難があった。
当初、本研究は「関数電卓で1946年のチャーニーらの試みを再現する」ことを目標としていた。しかし当時中学三年生だった筆者は、自身の貧弱な知識と限られた時間の中ではそれが不可能であることをじきに理解した。そのため、目標は「実際の数値実験用のコードを関数電卓に移植・実行する」ことに変化し、最終的に『気象研究ノート第110号 気象力学に用いられる数値計算法』（日本気象学会刊, 1972） に掲載・詳説されていたFORTRANによるコードを関数電卓用に移植・実行することで決着した。
HP50gの定格クロック周波数は、消費電力を抑えるため本来のプロセッサの性能をフルに利用せず75 MHzに留められているが、マニア有志により230MHzへのオーバークロックを可能にするツールが開発・配布されている。これを用いて計算時間の短縮を図り、またHP50gの性能を可能な範囲で最大限引き出すことを目指した。オーバークロックに伴い、消費電力の増加により計算途中で電源を喪失する可能性があることから、HP50g本体のMini-USB端子経由で外部給電しつつ計算を行うことにした。また、実行時間の計測にはHP50g内蔵の時計機能を利用した。
今回移植した数値実験コードの概要を述べる(参考文献：日本気象学会刊, 1972)。詳しい説明は参考文献に譲る。まず、予報モデルは「非発散順圧モデル」と呼ばれる種類のもので、ポアソン方程式として与えられる。 計算の大まかな流れは以下のようになる。
- 等圧面高度の初期値を生成する
- 等圧面高度から渦度・風速場を計算する
- 緩和法（加速リーブマン法）を用いて高度変化量を求める
- 高度変化量から新たな等圧面高度を時間外挿し、②に戻る
- 2から4までの流れを24回繰り返す
格子間隔は300キロメートルに設定し、格子点の数は縦10個・横21個の計210個としている。 計算領域は東西方向に周期的になっており、いわば本来円筒だった領域を切って広げて1つの長方形としたようになっている。初期値に関しては、移植元FORTRANコードに含まれていた初期場生成ルーチンをそのまま移植しており、原文を引用すると「ジェット型の基本流（J=5.5で最大帯域風速Uc=40 m/sec）に大規模な擾乱による風の場が重なった」ものとなっている。
また、使用した関数電卓であるHP50g Graphing Calculatorのスペックを、改めて以下に示す。
- CPU：Samsung S3C2410A (ARM920Tプロセッサ) 75MHz
- RAM：512 KiB
- フラッシュメモリ：2 MiB
- 画面サイズ：131×80 画素
- 端子：Mini-USBポート、IrDA、3.3V TTLレベル非同期シリアルポート
使用したプログラミング言語は、HP50gに内蔵されているUser RPLという言語である。この言語の特徴は、Forth言語のように逆ポーランド記法を採用していること、そしてスタックマシン上で動作するインタプリタ言語であることである。言語仕様の詳細説明はここでは省く（興味のある方は日本語版Wikipediaの”RPL（プログラミング言語）”の項を参照のこと）。
最終的に作成したコードは下記GitHubリポジトリにアップロードした。 https://github.com/stepney141/HP50g-works/blob/master/THOR1.txt
- Init：初期値となる等圧面高度の生成ルーチン（上記①）
- Sub：等圧面高度から渦度・風速場を計算するルーチン（上記②）
- Relax：緩和法で高度変化量を求めるルーチン（上記③）
- TimeD：時間外挿で等圧面高度を求めるルーチン（上記④）
- Cyclic：周期的境界条件を満たすためのサブルーチン その他、細かなサブルーチンの説明は省略する。
予報モデルの数学的取り扱いなどの詳しい説明は上記参考文献を参照されたい。
計算時間の測定値は13675秒、換算して約3時間47分となった。24回の外挿を繰り返したことを考えると、24時間予報としてはおよそ妥当と言える範囲の計算時間に収まったと言える。しかし、これはクロック周波数を約3倍に上昇させ、関数電卓の限界を引き出した上での結果であるので、通常の速度で予報を行った場合はより長い計算時間を要するはずである。 計算結果を評価するため、筆者のノートパソコンにてオリジナルのFORTRANコードを実行し、同様の24時間予報を試みた。結果、関数電卓による計算とほとんど一致した結果が得られ、関数電卓における数値実験が可能であることが確認できた。
関数電卓で数値実験が可能であることが確認された以上、関数電卓で本格的な数値予報モデルを実行できる可能性も十分にあるといえるが、そのために今後解決しなければならない課題も存在する。
今回のコードでは、オーバークロックをした上でなお計算におよそ4時間を要した。User RPLはインタプリタ言語であり、このことが速度を低下させていると思われる。同じくHP50gに搭載されているコンパイラ言語のSystem RPLや、有志が開発したコンパイラを使用してC言語でコードを書くことで、より高速な実行が可能になると思われる。コード自体を細かく最適化することで多少の高速化が図れる可能性もある。 また、HP50gは2006年に発表され、2015年に生産終了した旧世代の機種である。現在では同じHewlett-Packard社のHP Primeや、Texas-Instruments社のTI-Nspire CX CASシリーズなど、より高速なプロセッサ・より高性能なRAMを搭載した新機種を使用すれば、より容易に高速な実行が可能になると思われる。さらに近年はPython（MicroPython）が利用できる機種も増えており、これらを利用した計算も行いたいと考えている。
数値予報にとって、方程式の初期値となる気象観測データは極めて重要な役割を持つ。気象庁では世界的な気象観測網から受信したデータなどの初期値を有料で公開しているが、価格・利用条件等の諸般の理由から今回の研究では使用することができなかった。また、欧州中期予報センター（ECMWF)も初期値として利用可能なデータを公開しているという情報も得たが、研究当時（2015年）の筆者の語学力の不足・利用条件等の問題があり、これの利用も断念した。 今回は数値実験用のコードを移植することで初期値の問題を回避したが、本格的な数値予報を実行する上では実際の気象データは必要不可欠である。また、データを取得できたとしても、それを関数電卓でも利用可能な形に変換する工程が必須となる。初期値データの確保、そして変換手法の開発が今後の課題である。
今回利用した予報モデルは、大きな擾乱以外を無視するなど、様々な近似がされているごく単純なものであった。しかし、より高い精度の予報を行うためには、近似の度合いが弱まったより複雑なモデルを用いる必要がある。今後は他の予報モデルを実装し、関数電卓の限界をさらに追い求めていきたい。
- 数値予報と現代気象学（新田尚・二宮洸三・山岸米二郎著、東京堂出版）
- 数値予報の基礎知識（二宮洸三著、オーム社）
- 人と技術で語る天気予報史 数値予報を開いた〈黄金の鍵〉（古川武彦著、東京大学出版会）
- 嵐の正体にせまった科学者たち 気象予報が現代のかたちになるまで（John D.Cox著、丸善出版）
- 気象研究ノート第110号 気象力学に用いられる数値計算法（日本気象学会刊、1972）
- 数値予報論（岸保勘三郎著、地人書館）
- 数値予報新講（岸保勘三郎著、東京堂出版）