計算流体工学 大気環境シミュレーションの現状 佐藤正樹 埼玉工業大学／地球フロンティア研究センター satoh@sit.ac.jp ２００４年７月６，１３日 東京工業大学・総合理工学研究科

内容 大気環境シミュレーションの概観 数値計算手法 非静力学モデル：領域モデル 大気大循環モデル：全球モデル 気候シミュレーションの現状 　　気候シミュレーションの現状 気候モデリングにおけるコンピュータ環境 次世代大気大循環モデル

大気環境シミュレーションの概観

大気環境系の概観 気候系 大気、海洋、陸面（雪氷、土壌、植生）、海氷、生態系 空間スケール 　大気、海洋、陸面（雪氷、土壌、植生）、海氷、生態系 空間スケール 　地球 10000km ～ 積雲 1km～乱流 10m 時間スケール 地球の歴史 数10億年～第四紀　数100万年 人為的気候変化 100年 10年スケール気候変化、エルニーニョ ３年 年周期、季節変化、日変化、短期予報（数時間）

気候システム

数値シミュレーションによる気候研究

大気の空間スケール 積雲クラスタ～ 100km 積雲 １0km １km １0km 12740km

気候変動の時間スケール IPCC2001:　http://www.ipcc.ch

大気の数値モデル メソモデル：雲解像モデル 100km×100km：非静力学モデル 領域モデル 大気大循環モデル 　全球3次元： 静力学モデル

気象庁数値予報モデル　ｂｙ　室井

数値モデルの格子数と分解能 格子数 数100×数100×数10=数105 ⇒ 1000×1000×100=108 （次世代：ES） 　　数100×数100×数10=数105 　　⇒ 1000×1000×100=108 （次世代：ES） 領域幅と格子間隔： 　大循環 10000km：Δｘ＝100km ⇒ 10km 領域　　1000km：Δｘ＝10km ⇒ 1km メソ　　　100km：Δｘ＝1km

大気モデルの構成 （領域･メソモデル・大気大循環モデル） 力学過程 方程式系 空間差分法，時間差分法 数値拡散 物理過程 水循環過程（雲物理，積雲パラメタリゼーション） 放射過程（太陽放射，赤外放射） 乱流過程（接地境界層，subgrid乱流） 化学変化･物質輸送（オゾン，エアロゾル）

数値計算手法

方程式系 大気の基本バランスと波動 ＣＦＬ条件 音波・重力波の扱い 移流計算 擬スペクトル法・変換法 数値計算手法 方程式系 大気の基本バランスと波動 ＣＦＬ条件 音波・重力波の扱い 移流計算 擬スペクトル法・変換法

非圧縮ブシネスク系：工学or理論 完全圧縮性成層流体 非静力学方程式 プリミティブ方程式系 静力学近似 浅水波方程式系 準地衡風方程式系 流体を記述する方程式系 非圧縮ブシネスク系：工学or理論 完全圧縮性成層流体 非静力学方程式 プリミティブ方程式系 静力学近似 浅水波方程式系 準地衡風方程式系

大気の基本バランス 静力学平衡 地衡風平衡

工学における流体の方程式系（１） 非圧縮の式 運動方程式 連続の式

工学における流体の方程式系（２） ブシネスク方程式系 運動方程式 連続の式 熱エネルギーの式 状態方程式

大気の方程式系 圧縮性成層流体の方程式系 運動方程式 連続の式 熱エネルギーの式 状態方程式

大気の方程式系 回転系の圧縮性成層流体の方程式系 運動方程式 連続の式 熱エネルギーの式 状態方程式

波動 音波：圧縮性　 重力波：成層効果：浮力 ロスビー波：緯度方向の回転角速度の変化 音速 Brunt　Vaisala 振動数

音波 線形化、１次元、断熱 分散関係

重力波 ブシネスク、線形、２次元、断熱 分散関係

CFL（Courant-Friedrichs-Lewy）条件 移流方程式 差分化,　von　Neumann　解析

大気モデルのCFL条件 Δx=10～100km，Δz=1km 音波 　Δt<Δx/cs=30s～5min：Δt<Δz/cs=3～30s 重力波 　Δt< 1/N = 10min 移流：U=30m/s，W=数cm/s (メソ：W=10m/s) Δt< Δx/U=3～30min 　　(Δt<Δz/W=100s) ロスビー波：移流と同等

音波の差分解法 n: Forward （explicit, Euler） n+1： Backward (implicit) ｐn+1 ,un ：Forward-Backward (semi-implicit)

固有値λ解析： Forward Backward Forward-Backward

非静力学モデル （メソモデル・積雲解像モデル） 非静力学モデル （メソモデル・積雲解像モデル）

代表的な非静力学モデル RAMS http://www.aster.com MM5 http://www.mmm.ucar.edu/mm5 ARPS http://www.caps.ou.edu WRF http://wrf-model.org JMA-NHM (気象庁) 　　http://www.mri-jma.go.jp/Project.mrinpd/INDEXJ.htm http://pfi.kishou.go.jp (気象庁プラットフォーム) CReSS (名大) http://www.rain.ihas.nagoya-u.ac.jp/CReSS　

JMA-NHM：気象庁非静力学モデル 衛星可視画像 シユレーション結果（雲水量） 1997年1月22日

非静力学モデルの力学フレーム 音波の扱いによって方程式系が分類される 非弾性系 弾性系 or 完全圧縮系 陽解法（HEVE法） 完全陰解法（HIVI法） 水平陽解法、鉛直陰解法（HEVI法） 水平伝播する音波、重力波を分離 (split-explicit)

陽解法（HEVE法） 大気モデルは通常水平格子間隔の方が鉛直格子間隔より一桁大きい 水平格子間隔 Δx = 10km 　CFL条件⇒Δt = 10km/300m/s = 30s 鉛直格子間隔 Δz = 500m 　CFL条件⇒Δt = 500m/300m/s = 1.7s ⇒陽解法（HEVE法）は実用的でない。

陰解法（HIVI法）、非弾性系 cs→∞　とすれば非弾性系. ρsは z だけの関数とする. ポアソン方程式

ポアソン方程式の解法 ⇒ 緩和法：Jacobi法, Gauss-Seidel法, SOR法 共役勾配法 スペクトル法 マルチグリッド法：並列計算

ポアソン方程式の解法 緩和法：Jacobi法, Gauss-Seidel法, SOR法 共役勾配法 スペクトル法 マルチグリッド法 　　並列計算

水平陽鉛直陰解法（HEVI法） 鉛直方向の3重対角行列

Splitting method 音波・重力波に関係する項を短い時間間隔Δτで計算 他の項を長い時間隔Δtで計算 divergence damping 等数値拡散が必要

非静力学モデルの時間ステップ 音波の扱いが問題 陰解法（HIVI法）：時間ステップの制約がない HEVI法：水平伝播する音波により時間ステップが制約される ⇒ splitting 法により、他の項を長い時間ステップで計算 次に問題となるのは移流 音波 cs=300m/s, 風速 U=数10m/s

移流スキーム Euler法，フラックス型 有限体積法 ⇒ 保存系 semi-Lagrange 法 一般には保存は満たされない CIP Godnov, van-Leer，PPM，Lin and Rood（1996）

移流スキームの性質 diffusion：高精度スキーム dispersion： monotonicity(shape-preserving) 一般に高精度であればより振動的 高精度（2次、３次）かつ limiter により単調性をできるだけ確保

移流方程式 保存則 連続の式 ：移流型 ：フラックス型

Euler　法 ρq の領域積分が保存する セル境界におけるq の評価法により精度，単調性が依存 Courant数<1　を要求する

Semi-Lagrange 法 q が Lagrange 的に保存 セル内のqの補間方法により単調性，精度が依存 Courant数に制限がない (but Δt・du/dx<1:Lipschitz条件) ρq の領域積分は保存しない

保存型 Semi-Lagrange 法 ρq の領域積分が保存する q が近似的に Lagrange 的に保存 Coulant数に制限がない (but Δt・du/dx<1:Lipschitz　条件)

物理過程 物理過程 水循環過程（雲物理，積雲パラメタリゼーション） 放射過程（太陽放射，赤外放射） 乱流過程（接地境界層，subgrid乱流） 重力波抵抗 化学変化･物質輸送（オゾン，エアロゾル）

物理過程のスケール依存 雲物理過程 メソモデル：微物理過程：雨滴・雪霰の生成・成長過程 大循環モデル：積雲パラメタリゼーション 乱流過程 分解能数10mまで（1km以下でも使う）： LES 数km以上：統計乱流モデル 　　Moller and Yamada，１次，２次，３次クロージャ

例：積雲のパラメタリゼーション 1km 100km 積雲のスケール1km vs 大循環モデルの分解能100km 分解能以下の現象をパラメタリゼーション Arakawa and　Schubert：格子内で準平衡を過程 熱の式、水の変化式にソース項として考慮 パラメタリゼーションの方式は分解能に応じて変更 パラメタリゼーションは準理論、経験に基づくものであり、不確定性が大きい 1km 100km

大気大循環モデル

大気大循環モデル(AGCM) 静力学近似に従うプリミティブ方程式系 球面調和関数による擬スペクトル変換 毎ステップに格子空間とスペクトル空間のルジャンドル変換を行う 　微分の計算を格子空間で行う 重力波・音波をsemi-implicit法で計算 渦度，発散を予報し，ポアソン方程式を解いて速度を求める

AGCMの方程式系 経度方向運動方程式 緯度方向運動方程式 静力学平衡 連続の式 熱エネルギーの式 状態方程式 水蒸気の式

方程式の変形 鉛直座標： σ=p/ps =圧力/地表面圧力 ジオポテンシャル Φ=gz 　　⇒静力学平衡 渦度・発散

渦度・発散方程式 ⇒ ζ，δを予報⇒ψ，χについて解く⇒ｕ，ｖを求める 流線関数ψ 速度ポテンシャルχ 水平方向運動方程式 渦度方程式

スペクトル法とルジャンドル変換 球面調和関数Ｙ，ルジャンドル陪関数Ｐ 球面上のポアソン方程式

スペクトル法の切断波数 代表的な切断波数 T21 → J=64 Δx=625km T42 → J=128 Δx=312km 三角形切断： N=M，m:東西波数，n:南北波数（節の数） 非線形項の計算安定性（変換法）： 経度方向の格子点数 I=3M+1 以上：FFTとの親和性 緯度方向の格子点数 J=2N 以上 代表的な切断波数 T21　→ J=64　Δx=625km T42　→ J=128 Δx=312km T106 → J=320 Δx=125km T213 → J=640 Δx=62.5km

σ座標渦度発散系の方程式系

ＡＧＣＭ実験例：標準比較実験 海面温度の気候値を与え，長期間積分． Ｄｙｎａｍｉｃａｌ Ｃｏｒｅ：Held and Suarez(1994) 物理過程を簡略化. 力学過程の効果を調べる． 水惑星(Aqua　Planet)：Neale and Hoskins(2000) 地表面条件を簡略化．全球が海洋．水循環を調べる． AMIP（Atmospheric Model Intercomparison Project) http://www-pcmdi.llnl.gov　 海面温度の気候値を与え，長期間積分．

Dynamical Core：Held and Suarez DWD 　 T106L19 NCAR 　 Eulerian 　T42L49 　semi-Lagrange

Aqua Planet 実験 Neale and Hoskins 降水量 Aqua　Planet　実験 Neale　and　Hoskins 海面水温分布

AMIP Atmospheric Model Intercomparison Project

気候モデリングにおける コンピュータ環境

計算機環境 多様な計算機環境 スカラー ｖｓ ベクトル 単一PE ｖｓ 並列 PC，WS ｖｓ Super computer PC　or　PC　unix, Work　Station PC　cluster, WS cluster scalar parallel: massive parallel vector computer：　super computer vector parallel vector massive parallel: Earth　Simulator（640node×8PE）

ES (地球シミュレータ） 超並列ベクトル機：640node×8PE ピーク性能40TFLOPS 主記憶10TB

気候モデルの計算機環境の現状と今後 従来のモデルは単一PE，ベクトル機用にコーディングされている 今後は並列化対応が避けられない：アルゴリズムの選択に影響を与える：スペクトル法，ポアソン解法 スカラー機とベクトル機ではコーディング方法が異なる 最近の米国のモデルはスカラーパラレル用に作られている スカラーパラレル機の限界が指摘されている（USGCRP2000） ベクトルパラレルは速いが高い ベクトル機の今後の動向は不透明

数値モデルの格子数と分解能 格子数 数100×数100×数10=数105 ⇒ 1000×1000×100=108 （次世代：ES） 　　数100×数100×数10=数105 　　⇒ 1000×1000×100=108 （次世代：ES） 領域幅と格子間隔： 　大循環 10000km：Δｘ＝100km ⇒ 10km 領域　　1000km：Δｘ＝10km ⇒ 1km メソ　　　100km：Δｘ＝1km

次世代（今後１０年以内）に地球シミュレータ（ES)など超分散並列コンピュータによって従来より一桁上の分解能の計算が可能になる。 大気モデルは、統合化に向かう。 大循環モデルは格子間隔10km 領域モデルは格子間隔１ｋｍ 積雲を分解するのに必要な格子間隔：1km 次世代大循環モデルでも分解能は不十分

地球フロンティアにおける 次世代大気大循環モデル

次世代大気大循環モデル 全球雲解像モデル：NICAM：Nonhydrostatic　Icosahedral　Atmospheric　Model　（非静力学正20面体大気モデル） 地球フロンティア研究センター・地球環境モデリング研究プログラムで開発中 　地球シミュレータの性能を最大限生かしたモデル 高分解能全球３次元気候モデル：長期気候予測 準一様な格子モデル：格子間隔～10km 非静力学方程式系 積雲をできるだけ顕わに表現し，統計平衡を仮定した積雲パラメタリゼーションを用いない

次世代大気大循環モデルの開発 開発の手順 準一様格子の浅水波モデル（１層モデル） 保存則を満たす非静力学方程式モデル 物理過程の検討（雲物理，放射，乱流） 計算性能（ベクトル化，並列化，データ処理，画像解析）

格子構造 正20面体モデル 等角立方格子モデル

正20面体格子の生成 Number of grid points of n-level: glevel-1 glevel-2

level 7 (~56km) level 8 (~28km) level 9 (~14km) level 10 (~7km)

等角立方格子の生成 Number of grid points: 20x20 10x10 40x40 80x80

格子間隔の比較

浅水波モデルの結果 Williamson Test Case 5 Icosahedral grid Conformal cubic grid Spectrum model (T213)

非静力学コア 質量・エネルギー保存スキーム Cold bubble 実験 エネルギーの時間変化

Squall line （GCSS CASE1) 雲水量と雨 降水量の時間変化 降水・質量変化

Life Cycle of Extratropical Cyclone Exp.(1) Polvani and Scott(2002) glevel 10 Δx=7km glevel 6 Δx=112km glevel 8 Δx=28km

Life Cycle of Extratropical Cyclone Exp.(2) glevel 6 glevel 8 glevel 10

計算効率 Configuration Horizontal resolution : glevel-8 Vertical layers : 100 Fixed 　The used computer nodes 　 increases from 10 to 80. Results Red 　：　actual speed-up line Green　：　ideal speed-up line For the parallel performance check, I performed the model, changing the number of nodes from 10 to 80. This figure shows the parallel performance of our model. A horizontal axis is the number of PNs. PNs stand for Processor Nodes. One processor node have 8 vector processor. A vertical axis is a GFLOPS value. Flops stand for floating point number operation per second. When the amount of calculation of a model is divided equally to CPU, a gflops value is proportional to the number of CPUs, it becomes like a green line. In our model, gflops value increases like a red line. Ideally, the gflops value in 80 nodes is 8 times in 10 nodes. In our model, the gflops value in 80 nodes is 7.48 times in 10 nodes. The GFLOPS value at 80 nodes is 7.48 times faster than that at 10 nodes. Sustained performance at 80 nodes is 43.5%

スペクトルモデルとの比較 Elapse time of 1step for NICAM and AFES ＡＦＥＳ ( green line) Elapse time increases in the sense Ｏ(n3). Legendre transformation ＮＩＣＡＭ ( red line) in the sense Ｏ(n2). In all resolutions, NICAM is faster than AFES! Caution: what is the resolvable scale? 2Δx or 4Δx?

スペクトルモデルとの比較 NICAM gl7 gl8 gl9 gl10 gl11 450 225 113 57 29 1day time 　Maximum time step and one-day simulation time If L = 2 πR / N = 4 Δx NICAM gl7 gl8 gl9 gl10 gl11 450 225 113 57 29 1day time 6.70 32.1 210 1519 12200 AFES T159 T319 T639 T1279 T2559 400 200 100 50 25 1day time 8.02 27.9 184 1884 24930 時間刻み幅は、NICAM,AFESともT159でのHeld & Suarezテストでとり得る最大値をもとに、 見積もった値です。 1日積分は、1ステップあたりの計算時間とステップ数から計算した値です。 　Time step of NICAM can be larger than AFES At T1279 & glevel-10, NICAM is faster than AFES.

開発の現状 現状 保存型方程式に基づく正二十面体格子３次元球面非静力学モデルの力学コアを完成 全球を水平分解能 3.5km で分割した実験を実施 水惑星実験 海陸配置の導入、地形 雲物理過程に対する依存性 エアロゾル

まとめ 気候モデルは拡張しつづけている． 大気モデルは気候モデルのごく一部分である． 代表的な大気モデルは，領域モデル, 大気大循環モデルである．領域モデルはメソモデル、非静力学モデル，雲解像モデルを包含する総称となりつつある． 大気モデルは，力学過程，物理過程からなる． 領域モデルはフリーで利用できるものが数多くある． 大気大循環モデルも利用可能なものはあるが、物理過程などの不確かさを十分理解して使う必要がある． 数km格子になるまで，積雲パラメタリゼーションに伴う不確定性が避けられない．

複雑なモデルの評価のために標準実験が提唱されている． コンピュータの進歩により大気モデルの領域依存性は統合化されるであろう． 地球フロンティアで開発中の次世代大気大循環モデルは，3.5km格子で全球を覆う非静力学大循環モデルである． 近年のモデル開発は，数理解析（差分スキーム)，理論（大気科学），素過程（雲物理，放射，乱流），計算科学，システム開発などの総合的な分野の協力が必要である。