月刊「海洋」特集号原稿
非静力学モデルの力学フレームの現状と課題
斉藤和雄（気象庁数値予報課）
鉛直方向に加速度がない場合、圧力の鉛直勾配は流体に働く重力と釣り合っており、この関係を静力学平衡と呼ぶ。非静力学モデルは静力学平衡を仮定しないモデルで、扱う現象の水平スケールに原理的な制約がなくなるが、気圧の計算などに工夫が必要になる。ここでは研究と現業予報の双方に資する目的で開発が行われている気象庁の気象研究所/数値予報課非静力学モデルを主な例に、非静力学モデルの力学フレームに関する開発課題について論じる。
１．はじめに
気象数値モデルは、大気の状態を数値表現し、物理法則に基づいて計算して将来の状態を予測するものである。乾燥大気の場合、大気の状態は一般に、気温、気圧、三方向の風速、の５つの変数で表わされる。物理法則としては、密度を含めた６つの変数に対して、三方向の運動方程式（風、あるいは運動量の予報式）、状態方程式（密度の診断式）、連続の式（質量保存の式）、熱力学の式、の６つの基礎方程式により系が閉じることになる。現在の気象庁の現業的な数値予報モデルでは、鉛直方向の運動方程式を静力学平衡の関係に置き換える静力学近似が用いられている。静力学近似は、大気の大規模運動のように鉛直方向のスケールが水平方向に比べて十分小さい場合には良い精度で成り立つが、鉛直方向の不安定に起因する対流や小スケールの地形効果については、必ずしも成り立たない。非静力学モデルでは、水平分解能に原理的な制限がなくなり、小スケールの現象をより正確に表現することが可能になるが、３次元運動場の影響を考慮した気圧の計算が必要となるため、その取扱いに工夫が必要になってくる。
非静力学モデルに限ったことではないが、数値モデルで重要な点として
精度、計算効率、安定性、汎用性
などが挙げられる。これらを全て満足させることは必ずしも容易ではなく、また扱う問題や用いる計算機のアーキテクチャによっても事情は変わってくる。ここでは、日本で開発された非静力学モデルの内、現在最も広く研究用に利用されており、また近い将来の気象庁の現業用数値予報モデルとしても開発が行われている気象研究所/数値予報課非静力学モデル
(MRI/NPD-NHM)を例にとって、力学フレームを中心に開発課題について論じることにする。モデルの現状はMRI/NHM-NHMがベースとなるが、議論のうちのかなりの部分は、他のモデルにも共通する課題として、ここに記述する情報価値を持っていると思う。
２. 気象研究所/数値予報課非静力学モデル
気象研究所/数値予報課非静力学モデル
(MRI/NPD-NHM)は、1999年２月から気象庁気象研究所と気象庁数値予報課で研究・現業予報の双方に資するために共同開発を行っている汎用非静力学メソモデルである。「気象研究所予報研究部で開発された非静水圧モデル」
(Ikawa and Saito, 1991) の発展としての「気象研究所非静力学メソスケールモデル
(MRI-NHM)」（斉藤・加藤、1996; 1999）をベースに、数値予報課で開発されたスプリットイクスプリシット時間積分法（Muroi
et al., 2000）を組み込んだもので、基礎方程式としてマップファクターを含む完全圧縮方程式系
(Saito, 1997) を用い、氷相を考慮するバルク法の雲物理過程、乱流クロージャモデル、境界層過程などを含んでいる。パソコン上で動作する「ＮＨＭ統合環境」（上野ほか、2000）が作られているほか、2000年10月からは、気象庁モデル開発戦略の一環として、国内17の大学・国立研究機関にソースコードが配布され、利用・開発コミュニティが拡がりつつある。モデルの詳細は気象研究所技術報告
(Saito et al., 2001)にまとめられており、気象研究所ホームページ
(http://www.mri-jma.go.jp/Dep/fo/mrinpd/INDEXJ.htm)からも閲覧することが出来る。また本特集号でも8章と14章で関連する解説があるので、参考にして頂きたい。表１に、MRI/NPD-NHMの現状と開発課題についてまとめる。中短期開発課題にはすでに取り組まれているものも含まれている。以下の章では、そのうちの主だったものの詳細について述べる。
３．方程式系
１）非弾性方程式系モデル
非弾性モデル（anelastic model）では、場の変数を、静力学の関係が成り立つ水平一様な基本場とそこからの偏差に分けて取り扱い、基礎方程式系中の密度には基本場の値を代用する
(Ogura and Phillips, 1962)。密度の時間変化を考えないことにより、大気の圧縮性が除去され、これによって気象学的にはノイズである音波が解に含まれなくなる。方程式系は３方向の運動方程式と熱力学の式、気圧診断式としての連続の式の５つで完結する。非弾性近似は、通常の気象現象に対して概ね良い精度で成り立ち、再現実験などにも使われる。非静力学モデルの力学コアとして最も歴史があり、MRI-NHMでも多くのシミュレーションに用いられた実績がある。
問題点としては、方程式系を基本場の密度で線形化するため、密度変化が大きくなるような場では誤差が大きくなる。また、連続の式で圧縮性を除去することは、領域内の質量保存を前提としているため、このままでは平均気圧の変化を表現できない。気圧診断方程式が三次元の楕円型方程式となるため、後述するように大規模計算では計算効率を確保するのが難しいという問題もある。このため、最近のメソスケールモデルでは後述する圧縮系の方程式を使うものが増えてきている。ただし、狭い領域に対し超高分解能のシミュレーションを行う研究用LESモデルとしては、依然有用である。また圧縮系モデルの場合、運動場とバランスした気圧場を初期値に与える必要があるが、MRI/NPD-NHMではモデルの立ち上げを非弾性で行って気圧場を診断するオプションを持っており、ＮＨＭ統合環境で用いられている。フランス気象局で試験的に運用されているMeso-NH
(Lafore et al., 1998)ではマップファクタを考慮する非弾性モデルで、側面フラックスから平均気圧の変化を表現する工夫をしている。またDurran
(1989) は連続の式に非断熱膨張の効果を取り入れることで非弾性近似の精度を上げる提案をしている。
２）準圧縮系モデルと完全圧縮系モデル
準圧縮系モデル(quasi-compressible
model)では、連続の式に大気の圧縮性を認め、発散から気圧を予報するが、密度は基本場の値を代用する。密度の時間変化を考慮しないため、非弾性モデルと同様に状態方程式は陽には用いられない。準圧縮系モデルは、気圧を発散から予報するために音波を解に含むが、方程式系に基本場の密度による線形近似を行っているため、精度は非弾性系モデルとほぼ同じである。完全圧縮系モデル(fully
compressible model) では、連続の式として厳密な式を用いる。気圧の予報方程式は準圧縮系モデルと同様な形をとるが、基本場密度による近似は行われない。完全圧縮モデルは、方程式系に近似を含まないという点で優れているが、反面、密度の時間変化を許すため質量の保存が数値計算の誤差に敏感になるため、計算上の注意が必要である。MRI/NPD-NHMの圧縮系バージョンでは完全圧縮系を標準としているが、準圧縮系のオプションも残している。静力学バージョンは現状では、非弾性バージョンのみに用意されている。
４．圧縮系モデルの積分スキーム
圧縮系モデルは音波を解に含む。大気中の音波の伝播速度は300m/s以上と大きいので、そのままでは安定に数値計算を行うためには、時間積分の時間刻みを小さな値にとらなければならなくなる。気象モデルでは、実用的な計算効率を達成するために大別して以下の２種類の取扱いを行っている。
1) HI-VI法
HI-VI法はhorizontally
implicit - vertically implicit schemeの略で、気圧方程式とともに運動方程式の中で音波に関連する項を鉛直・水平方向ともに陰解法（インププリシット解法）して、時間刻みの制限を緩和する計算手法である。音波関連項のみを陰解法するので、セミインプリシット法とも呼ばれる。Tapp
and White (1976)により提案された手法で、気圧を陰解法するために、形式的に非弾性モデルの気圧診断方程式に似た楕円方程式を解く必要がある。MRI-NHMでは非弾性方程式に代わって用いられており(Saito,
1997)、雲解像モデルとしての研究用モデルの計算スキームとして実績がある。
問題点としては、地形に沿った座標系を用いる場合、音波関連項の陰的取り扱いを線形項のみに限る場合、急峻な地形に対して音波の安定化が十分に行えない場合があることがある。また楕円方程式の解法に用いられる直接法の固有関数展開は、広義のスペクトル法で、大規模計算では計算量が増大する。これらの対策としては、気圧傾度力項を全て陰解法して、気圧楕円方程式を逐次近似法で解くことが考えられる。
2) HE-VI法
HE-VI法はhorizontally
explicit - vertically implicit schemeの略で、音波関連項を分割して小さな時間刻みで解くとともに、鉛直方向に対してのみ陰解法して計算効率を高める手法である。音波項を分割（タイムスプリット）して水平方向には陽解法する（イクスプリシットに解く）ため、スプリットイクスプリシット法とも呼ばれる。気圧方程式は鉛直の１次元楕円方程式となる。HE-VI法では水平格子間隔が鉛直方向のそれに近づく場合は音波についての時間刻みは小さな値をとらなければならない。
通常、小さな時間刻みの積分はKlemp
and Wilhelmson (1978)に基づいてforward-backwardで行うが、Satomura
(1989)は小さな時間刻みも中央差分で解く方法を用いて時間差分の精度を高めている。また米国で開発中のWRF
(http://wrf-model.org/)では、大きな時間刻みの積分にRunge-Kutta法を用いて時間差分の精度を高めている。
５．移流スキームの精度
MRI/NPD-NHMを始めとする殆どの非静力学モデルは格子点法を用いている[1]。これに対して、現在の気象庁の現業数値予報モデルは全てスペクトル法を用いている。スペクトル法の最大の利点は、基底関数の空間差分が数学的に厳密に求められるため、移流の表現に位相速度による波数分散性がないことである。MRI/NPD-NHMでは、格子点で通常用いられる中央差分形式に加え、差分表現による波数分散性に伴うにせの極値の出現を抑えるための風上値による補正スキーム
(Kato, 1998) が用意されている。他の非静力学モデルでは、4次精度などより高次の差分形式を用いるものや、セミラグランジアン法を用いる非静力学モデルも作られている。比較的最近の数値計算技術としては、結合コンパクト差分
(Chu and Phan, 1992) や有限体積法 (Carpenter et al.,
1990など)、セミラグランジアン法の発展としてのＣＩＰ法（Cubic
Interpolate Profile＝微分量もセミラグランジアン的に移流させて精度を高める手法、Yabe
et al., 2001など）などが提唱されてきており、これらの非静力学モデルへの応用は今後の課題である。
６．重力波の扱い
重力波は安定成層大気で浮力を復元力とする波で、ブラントバイサラの振動数Nを持つ。重力波は伝播速度が速い上に、陰的取り扱いを行わない場合、安定に積分できる時間刻みは1/Nで制限される。通常の大気では1/Nは1分程度だが、成層圏や逆転層など安定度の大きな場所では局所的には30秒以下になる場合がある。これまでMRI/NPD-NHMでは重力波については特別な扱いをしていなかったが、HE-VI法の小さな時間刻みで重力波を計算するオプションを最近作成している。HI-VI法の場合には、重力波を音波とともにインプリシットに扱い、セミラグランジアン法と組み合わせることにより、時間刻みの制限が格子間隔によらずに決められる（Golding,
1992など）ことが期待出来る。
７．座標系
現在地形を含む非静力学モデルの鉛直座標系として最も多く使われているのは、地形に沿った座標系
(z*座標系)である。この座標は地表で０、モデルのトップで水平な高度面となるので、境界条件の扱いが容易になるという利点を持っている。ただし斜面上では座標面が傾きを持つため、水平傾度力項などをモデル面に沿った差分と鉛直差分の和の形式で表現する必要があり、誤差の原因となる。この問題を解決するため、下部境界条件が複雑になっても傾斜を持たない水平な座標（ステップ座標）で計算を行う非静力学モデルの開発も行われている。また地形に沿った座標系を改良して、モデル中層のある高さ以上で水平な座標系となるようなハイブリッド系を用いることも計算精度の向上に有効と思われる。Satomura
(1989)は、境界適合格子生成により、鉛直にも曲線座標を用いて地形に沿う水平座標との直交性を改善して急傾斜の場合の精度を改善している。いくつかの非静力学モデルでは、従来の静力学モデルとのつながりを考慮して、鉛直座標として気圧座標系を用いている。
水平の座標の取り方は、領域モデルでは、等角投影による直交座標を用いるのが一般的だが、緯経度座標を用いるモデル
(Doms and Schaettler，1997など)もある。Saito
(2001)はMRI/NPD-NHMのマップファクタの扱いを拡張して、等緯経度座標にも用いられるようにする試みを行っている。全球非静力学モデルの開発に関連して、正20面体測地線格子を用いるモデルの開発も始まっている（4章参照）。
８．境界条件
境界条件は領域モデルのシミュレーションを成功させるためには本質的に重要である。非弾性系モデルの場合、内部の気圧場は楕円方程式を解くことによって求められるが、その解は境界条件によって規定される。圧縮系モデルでは気圧は発散を生成項に持つ予報方程式を解くことによって得られるが、HI-VI法では気圧傾向方程式は非弾性と同様な境界値問題としての楕円型方程式となる。HE-VI法では側面境界の影響は音速でモデル内部に伝播する。MRI/NPD-NHMの場合、現状は風、温位、水蒸気のみを一方向ネスティングしている。雲物理変数、特に雲水や雲氷を親モデルから与えてやることは重要な課題である。また乱流エネルギーは内部摩擦を通じてモデルの接地境界層構造に影響するため、ネスティング時の安定性に影響する。これらについての改良の試みが既に始まっている。台風の研究シミュレーションでは、双方向ネスティング技術の開発も重要で、MRI/NPD-NHMについても、気象研究所で取り組み始められている。
上部境界条件は、多くのモデルでは現状では摩擦のない断熱固定壁にレーリーダンピングによる吸収層を設けている。モデル上端に放射境界条件を入れてより現実的な境界条件にしていく必要がある。
下部境界条件は物理過程と結びついており今回のレビューの範囲を逸脱するため、ここでは割愛する。
９．全球非静力学モデル
近年、計算機とモデルの進歩により、非静力学モデルについても、広領域かつ長期間の積分が可能になってきており、雲解像モデルの気候研究への応用も試みられてきている。一方、
全球モデルは高分解能化がすすんでおり、将来は現業数値予報に用いられる全球モデルの水平分解能が静力学近似の限界といわれる10kmに近づくことも予想される。また球面スペクトル法を用いるモデルでは、波と格子の変換に関する計算量が高分解能化により急激に増大するという問題がある。このような事情に関係して、全球非静力学モデルの開発が試み始められている。MRI/NPD-NHMについても、球面直交曲線座標系バージョンをベースとした全球モデルの開発が始まっている
(Saito, 2001)。ただし、極の扱いをはじめ本格的な高分解能全球モデルとしては多くの課題を残している。
全球非静力学については4章に詳しく述べられている。
10．おわりに
紙数の制約もあり、ここでは非静力学モデルの開発課題のうちごく一部のみを簡単に紹介するにとどまった。大規模計算に関しては、地球シミュレータ（13章参照）に代表されるような分散主記憶型並列計算機が用いられるようになってきており、モデルを効率的に動作させるための並列化の問題が重要な要素となってきている。またシミュレーションの多くのケースでは、モデルの計算結果は物理過程に最も強く依存している。さらに現業予報のためには、データ同化と初期値作成の問題は予報精度に決定的に重要であるが、これらについては、研究と現業予報のそれぞれについて、多くの課題が残っている。予報モデルは、力学コア、物理過程、入出力ファイルの処理、など様々なプログラムが組み合わさったソフトウェアであり、さらに現在の数値予報は、データの収集と品質管理、データ同化、可視化と結果の配信など、全体で巨大なシステムを構成しており、ここで述べた課題はそれらの一部についてのものに過ぎないことをお断りしておく。
参考文献
[1]
Carpenter, R.L., K.L. Droegemeier, P.R. Woodward, and C.E. Hane (1990):
Application of the piecewise Parabolic Method (PPM) to meteorological modeling.
Mon. Wea. Rev., 118, 586-612.
[2] Chu P.C. and Fan, C. (1998): A
Three-Point Combined Compact Difference Scheme : J. Comput. Phys., 140, 370-399
[3]
Doms, G. and U. Schaettler (1997): The nonhydrostatic limited-area model LM
(Lokal-Modell) of DWD. Part I: Scientific Documentation. Deutscher Wetterdienst, 155pp.
[4]
Durran, D. R. (1989): Improving the anelastic approximation. J. Atmos. Sci., 46, 1453-1461.
[5]
Golding, B.W. (1992): An efficient non-hydrostatic forecast model. Meteorol.
Atmos. Phys., 50, 89-103.
[6] Ikawa, M. and K. Saito (1991):
Description of a nonhydrostatic model developed at the Forecast Research
Department of the MRI. Technical Reports of the MRI, 28, 238pp
[7]
Juang, H.M. (1992): A spectral fully compressible nonhydrostatic mesoscale
model in hydrostatic sigma coordinates: Formation and preliminary results.
Meteorol. Atmos. Phys., 50, 75-88.
[8]
Kato, T., 1998: Numerical simulation of the band-shaped torrential rain
observed southern Kyushu, Japan on 1 August 1993. J. Meteor. Soc. Japan, 76,
97-128.
[9]
Klemp, J.B. and R.B. Wilhelmson (1978): The simulation of three-dimensional
convective storm dynamics. J. Atmos. Sci.,
35, 1070-1096.
[10]
Lafore, J.P., et al. (1998): The Meso-NH Atmospheric Simulation System. Part 1:
adiabatic formulation and control simulations. Ann. Geophysicae, 16,
90-109.
[11] Muroi, C., K. Saito, T. Kato and H.
Eito (2000): Development of the MRI/NPD nonhydrostatic model. CAS/JSC WGNE
Research Activities in Atmospheric and Oceanic Modelling, 30, 5.25-5.30.
[12]
Ogura, Y. and N.A. Phillips (1962): Scale analysis of deep and shallow water
convection in the atmosphere. J. Atmos. Sci., 19, 173-179.
[13] Saito, K. (1997): Semi-implicit
fully compressible version of the MRI meso-scale nonhydrostatic model.
--Forecast experiment of the 6 August 1993 Kagoshima torrential rain--.
Geophys. Mag. Ser.2, 2, 109-137.
[14] 斉藤和雄・加藤輝之
(1996): 気象研究所非静水圧 ネスティングモデルの改良について.
天気、43,
369-382.
[15] 斉藤和雄・加藤輝之
(1999): 気象研究所非静力学メソスケールモデル、気象研究ノート、196, 169-195.
[16] Saito, K., T. Kato, H. Eito and C.
Muroi (2000): Documentation of the Meteorological Research Institute/Numerical
prediction Division unified nonhydrostatic model. Technical Reports of the MRI,
42, 133pp.
[17] Saito, K. (2001): A global version
of the Meteorological Research Institute/Numerical Prediction Division
Nonhydrostatic Model. CAS/JSC WGNE Research Activities in Atmospheric and
Oceanic Modelling, 31, 6.20-6.21.
[18] Satomura, T.,
1989: Compressible flow simulation on numerically generated grids. J. Meteor.
Soc. Japan, 67, 473-482.
[19]
Tapp, M.C. and P.W. White (1976): A non-hydrostatic mesoscale model. Q. J. R. Meteorol. Soc., 102, 277-296.
[20] 上野幹雄・川畑拓矢・酒井亮太・白川栄一・石田純一・斉藤和雄
(2000): パソコン版気象研究所非静力学モデル. 天気、47,
289-294.
[21] Yabe, T, R. Tanaka, T. Nakamura and
F. Xiao, 2001: An exactly conservative semi-Lagrangian scheme (CIP-CSL) in one
dimension. Mon. Wea. Rev., 129, 332-244.
表1 気象研究所/数値予報課非静力学モデルの現状と開発課題
|
分 類 |
現状 |
中短期開発課題 |
長期課題、その他 |
---------------------------------------------------------------------------------------------------------------------------------
|
基礎方程式系 |
非弾性方程式系 |
基本場密度を用いる、 静力学バージョンあり |
-- |
マップファクターの考慮 連続式に熱膨張を考慮 |
|
|
準圧縮方程式系 |
基本場密度を用いる |
-- |
フラックス形式誤差の補正 |
||
|
完全圧縮方程式系 |
等角投影についてのマップファクターを考慮 |
球面直交座標への対応 |
静力学バージョンの追加 |
||
|
予報変数 |
運動量、気圧、温位、乱流エネルギー |
混合距離の予報変数化 |
全予報変数に密度をかけて完全フラックス形式に |
||
|
鉛直座標系 |
地形に沿った座標系(z*系) |
下層z*、上層zのハイブリッド化 |
境界要素法 ステップ座標系 |
||
|
鉛直格子構造 |
Lorenz型 |
-- |
Charney-Philipsとの比較 |
||
|
水平格子構造 |
Arakawa Cグリッド |
-- |
CDスキームの導入 |
||
|
移流項の計算 |
中央２次フラックス形式 風上値による補正 ４次差分は逐次版のみ |
高次差分の並列化 時間刻みのスプリット |
セミラグランジアン法 CIP法 コンパクト差分法 |
||
|
時間積分法 |
リープフロッグにタイムフィルターを併用 |
-- |
Runge-Kutta法の併用 |
||
|
音波の扱い |
HI-VI法 |
水平方向には基本場成分をインプリシット化 |
-- |
完全インプリシット化 気圧方程式に逐次解法導入 |
|
|
HE-VI法 |
Forward-backward |
-- |
短い時間刻みに中央差分 |
||
|
重力波の扱い |
HI-VI法 |
特別な扱いなし |
-- |
インプリシット化 |
|
|
HE-VI法 |
特別な扱いなし |
スプリットして短い時間刻みで積分 |
インプリシット化 |
||
|
乱流の扱い |
レベル2.5の乱流クロージャモデル |
混合距離の予報変数化 |
レベル３クロージャモデル |
||
|
上部境界条件 |
摩擦のない断熱固定壁＋レーリー摩擦吸収層 |
-- |
開放境界条件 放射条件の導入 |
||
|
側面境界条件 |
1方向放射ネスティング レーリー摩擦吸収層 |
双方向ネスティング |
|
||
|
初期値化 |
RSMの予報値を内挿 非弾性方程式を解く |
デジタルフィルター |
変分法（3次元、4次元）の導入 |
||
|
計算拡散 |
４次の線形拡散 |
-- |
極域フィルター |
||