レッスン
理論 — 格子上の気象
予報モデルは、大気を格子上の数値として保持し、それを時間ステップごとに進めていきます。モデルが行う最も単純な処理は移流です。これは風がパターンを運ぶ現象であり、 ∂u/∂t + c ∂u/∂x = 0 で表されます。このページでは、円周 160 km のリング上でこれを解きますが、この設定には重要な意味があります。移流には厳密解が存在するからです。時間 t 後の場は、初期の場を c·t だけシフトさせたものにほかなりません。したがって、ここで報告されるすべての誤差は、より細かい計算結果に対してではなく、真理に対して測定されています。クーラン数 ν = c·Δt/Δx は、風が1ステップでいくつの格子をまたぐかを示しています。そして CFL条件 ν ≤ 1 は、ここにある陽解法スキームが従うか、あるいは破綻するかの法則です。クーラン、フリードリヒス、そしてレーヴィがこれを発表したのは 1928 年のことであり、誰かが数値予報を実行する 18 年も前のことでした。
各記号の意味
c- 風速 (m/s)。ここでは定数ですが、実際のモデルでは設定値ではなく解の一部となります
Δx- 格子間隔。 160 km を格子数で割った値です
Δt- 時間ステップ (秒)
ν- クーラン数、 c·Δt/Δx 。このページのすべては、この値が 1 より大きいか小さいかにかかっています
N- リング上の格子数。 10 から 400 まで
- 前提
- 一つの方程式、一次元、単一の風速、そして境界のないリング。厳密解がわかっているため、報告される誤差はスキームそのものに起因する純粋なものです。またここでは、各々の新しい値が古い値のみから構築される陽的なスキームを前提としています。これこそがCFL条件が真に対象とする前提であり、実際のモデルが部分的に回避しているものでもあります。
- 成り立たない場合
- 予報センターが扱うほとんどの場面において、この前提は成り立ちません。実際のモデルは複数の波の速度を同時に扱い、その中で最も速いのは約 340 m/s の音速です。誰も予報しようとしていない音波のためにCFL条件に従おうとすれば、気象現象が必要とするよりも 30倍 も小さな時間ステップを強いられます。運用モデルが半陰的解法を組み合わせてこれらの波を別個に扱うのはこのためです。風も一定ではないため、 ν は平均的な風ではなく、格子上のどこ吹く最も速い風に対しても条件を満たす必要があります。また、領域限定モデルには境界が存在し、そこでどう振る舞うべきか自体が独立した研究テーマとなっています。陰的スキームであれば、配列を走査する代わりに毎ステップ連立方程式を解くという代償を払うことで、CFL限界を完全に回避することができます。
全プロセスの詳細解説
-
4 km のモデルには 40 分かかる。では 500 m にはどれだけかかるか？ 7 ステップ
ある 4 km 予報モデルは計算機の 40 分間で対象領域をカバーします。では、同じ予報を 500 m で行うとどれだけのコストがかかるでしょうか。 先を読み進める前に計算し、三つの要因のうちどれが想定外であったかに注目してください。安定限界のぎりぎりを押して、ツールの表を 1 km の行に合わせておきましょう。
-
まず格子数から始めます。間隔を半分にすると各辺の数は 2倍 になります。予報領域には2つの水平次元があるため、1回の半減で格子数は 4倍 になります。 4 km で 40 × 40 なら、 2 km では 80 × 80 です。
-
しかし、時間ステップは自由に保てるものではありません。これが多くの人の陥る罠です。CFL条件がそれを格子間隔に結びつけているため、 Δx を半分にすれば許容される最大の Δt も半分になります。風速 10 m/s 、 ν = 0.9 のとき、 4 km の格子では 360 s が許されますが、 2 km の格子では 180 s になります。シミュレーション1日あたりのステップ数は 240 から 480 に増加します。
-
これらを掛け合わせます。格子数が 4倍 でステップ数が 2倍 なので、1回半減させるごとの計算量は 8倍 となります。 4 km での1日あたり 384,000 回の格子更新に対し、 2 km では 3,072,000 回です。表の最後の列にはどの行にも 8.00 と記載されていますが、これは見積もりではなく実際の倍率です。
-
4 km から 500 m へは3回の半減となるため、倍率は 8³ = 512 になります。表は端から端までこれと一致しています。 4 km での1日あたり 384,000 回の更新が、 500 m では 196,608,000 回になります。
-
40 分は 40 × 512 = 20,480 分になります。つまり、たった1日分の気象に対して2週間かかるということです。これこそが、予報格子が誰もが望むほど細かくはならず現在のサイズにとどまっている理由であり、解像度を 2倍 にしようとするセンターが 8倍 の性能の計算機を買わなければならない理由に対する答えなのです。
-
この見積もりはその単純な構造に対して正直です。格子更新1回あたりの計算量が変わらないことを前提としており、これは一つの移流項についてはほぼ正当ですが、物理モデルが複雑になるにつれて過小評価になりがちです。また、 ν が 0.9 にとどまることを仮定していますが、実際のモデルでそれを約束することはできません。 ν は格子上で最も速い風に対して条件を満たさなければならず、ジェット気流がモデル開発者に都合よく配慮してくれることなどないからです。
-
このページが触れない点について補足します。この算術的現実に直面したセンターは、ただ大人しくコストを払うわけではありません。増大する係数の一部は、毎ステップ連立方程式を解くコストと引き換えに速い波のCFL条件を回避する半陰的スキームによって取り戻されます。また、ネスティングや可変解像度格子を用いて、必要な場所にだけ解像度を費やすことでも相殺されます。どちらも実在する手法ですが、どちらも無料（タダ）ではなく、このページではそのどちらも計算に入れていません。
解答
計算量は 512 倍。したがって、たった1日の予報におよそ2週間かかります。 1回の半減につき格子数が 4倍、ステップ数が 2倍。これが3回繰り返されるため 8³ = 512 となります。誰も予想しないのが第三の要因です。すなわち、クーラン条件があなたに代わって時間ステップを消費してしまうのです。
-
学習の道すじ
天気予報の仕組み
参考文献 (4)
- The condition itself, published in 1928 — eighteen years before anyone integrated a forecast, and about hyperbolic equations rather than about weather: R. Courant, K. Friedrichs and H. Lewy, "Über die partiellen Differenzengleichungen der mathematischen Physik." Mathematische Annalen 100(1), 1928, pp. 32-74. DOI 10.1007/BF01448839
- The first numerical forecast that worked, and the paper that had to choose a grid and a time step to make it work: J. G. Charney, R. Fjørtoft and J. von Neumann, "Numerical Integration of the Barotropic Vorticity Equation." Tellus 2(4), 1950, pp. 237-254. DOI 10.3402/tellusa.v2i4.8607
- Richardson's 1922 attempt, which is where the subject starts and where the arithmetic first bit — his six-hour forecast of pressure change came out 145 hPa, against a real change of almost nothing: L. F. Richardson, Weather Prediction by Numerical Process. Cambridge University Press, 1922.
- Where the stability analysis of these two schemes is done properly, including why a second-order centred difference is unconditionally unstable for pure advection: D. R. Durran, Numerical Methods for Fluid Dynamics: With Applications to Geophysics, 2nd edition. Springer, 2010.