京を使った高解像度全球大気シミュレーションの成果とこれからの展望
富田浩文（ 理研・計算機構 ）NICAM working team&Team SCALE

内容• 2011年ＳＳ研でお話しした内容– 大気モデリング研究でのこれまでの進展• ES出現前はどうだったか？
• ES出現後はどうなったか？
• 今後は（京では）？

• 今回ＳＳ研でお話しする内容– 京コンピュータで得られた二つのlandmark的な成果• 全球1kmを切る超高解像度実験( Miyamoto et al. GRL2013)• 雲解像することによるＭＪＯの予測可能性へ前進（Miyakawa et al. Nature Comm. 2014）• その先には、、、あくまで、富田の方向性– 全球雲解像から、全球大渦解像へ

気候モデルでのHPCの必要性

計算量が増大する三つの軸：• 高解像度化– サブグリッドスケール現象（乱流や雲）のパラメータを排除• 雲対流系まで解像• e.g. NICAM,NHM,WRF• モデル要素の複合化・精緻化– 精緻化されてきたプロセスモデルの統合• 大気・海洋・陸面モデル＋炭素循環＋生態系＋…
• e.g. MIROC-ESM• アンサンブル数の増大、積分期間の延長– より統計的な評価のため• 非線形系の予測には必須

バランスが重要雲はなぜ重要？

• 大気大循環のエンジンとして :– 熱帯の雲は、熱を中緯度方面へ運ぶ重要な役割をになう。
• 熱帯の雲： 階層構造
– クラウドクラスター、スーパークラウドクラスター、熱帯低気圧、MJO
まずは、積乱雲系のメカニズムや予測可能性お話•Parasol effect
Reflection of
こちらの気候 solar incident
学系のお話
はのちほど。

気候学的にはエネルギー収支に大きく影響:– 二つの効果!
• パラソル効果 :日射を跳ね返す。
• 温室効果 : 地面からの赤外放射をいったん吸収。地面へ再放出する。

– エアロゾルとの関係.
• エアロゾル間接効果: 光学的厚さを変える。

Greehouse effect4
Emision of infrare熱帯の雲活動概観
大気現象の空間スケールと階層構造
スーパークラウドクラスタ～ 数1000km  MJOと台風

クラウドクラスタ～ 数100km

積乱雲
１0km

１0km
地球の直径12740km

１km

全球雲解像モデルＮＩＣＡＭの概要
3. 今回のシミュレーションの特徴と成果
NICAMとは？

Non-hydrostatic ICosahedral Atmospheric Model
非静水圧平衡(≒雲解像)

二十面体

従来

大気モデル

NICAM

従来の全球モデル

格子間隔 数十 ～ 数百 km
雲を大規模場の関数として簡単化

NICAM

格子間隔 数百m ～ 十数 km
雲の対流を陽に計算

ＮＩＣＡＭのアドバンテージ• 従来型AGCM– 解像度100km～数10km– 積雲パラメタリゼーション＋大規模凝結• 統計的/経験的法則に基づく• モデルの不確定要素No.1• 雲を一つ一つ解像することの必要性– これにより、雲の階層構造を表現– パラメタリゼーションの不確定性排除

NICAMのモデリング戦略
NICAM の開発: ~2000

まだまだ、改良が続いている。
開発の基本的ポリシー

雲そのものを露わに表現すること– 高解像度化、物理過程の精緻化は当然の方向• 正20面体格子の採用– 準一様格子点法を用いることは、計算効率上不可避。スペクトル法では、もはやしんどい。

• 非静力学ダイナミカルコア– 一個一個の雲スケールでは、激しい上昇下降流が起こっている。

• 雲微物理による直接表現:– 積雲パラメタリゼーションの不確定性を排除し、雲力学を露わに表現すること（凝結、蒸発、衝突・併合）で現象を理解する。

NICAMの開発歴• NICAM力学過程 ( ~2000)• 水平格子 : 修正型正20面体格子( Tomita et al. 2001, 2002 JCP)
• 力学系 : 全エネルギー保存考慮の非静力学方程( Satoh 2002, 2003 MWR)
 3次元力学コア( Tomita and Satoh 2004 FDR)

– 2004 ： 物理過程を含み最初のバージョン完成• 2004以降• 高次精度＋単調性考慮の移流スキーム(Miura 2007)• Consistency With Continuity( Niwa et al. 2011)– 混合比が移流スキームの拡散などに左右されない。

• 水平格子ファミリー( ストレッチ格子, Tomita 2008, Iga andTomita 2013 )• NICAM-DC は、現在BSD2 license化し、以下で公開。

– http://scale.aics.riken.jp/nicamdc/

Summary of NICAM current status
Ref. Satoh et al. 2014 PEPS in press, / Satoh et al. 2008 J. Comput. Phys.
Tomita & Satoh 2004 Fluid Dyn. Res. /Tomita et al. 2011, ECMWF workshop proceeding Dynamics
Governing equations

Fully compressible non-hydrostatic system

Spatial discretization
Horizontal grid configuration
Vertical grid configuration
Topography

Finite Volume Method
Icosahedral grid with spring dynamics smoothing(Tomita et al. 2001/2002)
Lorenz grid
Terrain-following coordinate

Conservation

Total mass, total energy (Satoh 2002, 2003)

Temporal scheme

Slow mode － explicit scheme （RK2, RK3）
Fast mode － Horizontal Explicit Vertical Implicit scheme Physics
Turbulence/shallow clouds

MYNN 2.0,2.5(Nakanishi and Niino 2004) modified by Noda(2009)

Surface flux

Louis (1979), Uno et al. (1995)

Radiation

MSTRNX (Sekiguchi and Nakajima, 2005)

Cloud microphysics

NSW6 (Tomita 2008) --- 6 caegories of water ( 1moment-bulk)
NDW6(Seiki et al.2013) --- 6 caegories of water ( 2moment-bulk)

Cloud parameterization

NONE

Surface process

MATSIRO(Takata et al.)10熱帯域におけるマッデン・ジュリアン振動の１ヵ月予報が実現可能であることを実証～「京コンピュータ」+ 次世代型超精密気象モデルによる予測 ～3次元可視化ソフト Vapor

Nature Communications 2014年5月6日号 ハイライトペーパーとして掲載
著者： 宮川知己, 佐藤正樹, 三浦裕亮, 富田浩文, 八代尚, 野田暁,
山田洋平, 小玉知央, 木本昌秀, 米山邦夫

MJO の影響による近年の洪水被害
宮川氏（東大）提供

2010年 スリランカ洪水
被災者２０万人

2014年（先月） ソロモン諸島
被災者５万人以上

2011年 タイ洪水
２３０万人に影響

2013年 ジャカルタ
３万人以上避難2012年 オーストラリア
１万人以上孤立ＭＪＯの日本への影響

宮川氏（東大）提供

２０１１年秋の異常高温と厳冬
暖2011年10月下旬～11月 ＭＪＯがインド洋で大きく発達
持続的な加熱の影響がジェット気流を蛇行させて記録的な暑さに。

寒高知大学気象情報頁の画像に加筆
http://weather.is.kochi-u.ac.jp/

気象庁報道発表資料

ＭＪＯの活動がインドネシア付近で活発だった
12月はインド洋の対流は抑制され、一転して寒冷化。

成果をザックリ言うと・ 熱帯の雨の大部分を支配する現象の予測が約１ヶ月先まで可能なことが示された。
・ 日本を含む中緯度の予測精度向上にもつながると期待できる。
Capacity vs Capabilityの観点
これらは、地球シミュレータ時代に可能性にとどまっていた熱帯気象予測において、京を駆使して、中規模ランを大量に流すことによって統計的な知見へ引き上げたこと。

3. 今回のシミュレーションの特徴と成果
今回の計算の位置づけ

２００７年

今回

（NEC イメージ図）
JAMSTECRIKEN/AICS

地球シミュレータ(初代）京コンピュータ

NICAMを使い、世界に先駆けてMJOの再現に成功（Miura et al. 2007, Science)

NICAMによる多数のシミュレーションが可能に。

ただし、多くの事例を行うのは困難!
一発計算で可能性示す。

ＥＳ時代のCapability computingの典型

改良の進んだNICAMで過去１０年の冬季MJOすべての予測実験（19事例 54本のシミュレーション）
によりMJO予測性能を調査

京時代のCapacity computingの典型可能性を示すデモンストレーションから真の科学的成果への昇華

3. 今回のシミュレーションの特徴と成果
顕著な成果(1)宮川氏（東大）提供

NICAMによるMJO予測スキル（論文図に加筆）

有効な予測スコア （0.6以上） を 26 – 28 日程度維持

3. 今回のシミュレーションの特徴と成果
顕著な成果(2)宮川氏（東大）提供

54本の実験におけるMJOの位置（phase）別の平均降水分布（論文図に加筆）（平均16日目）（平均28日目）

前スライドの予測スコアにあまり反映されない降水の水平分布も観測とよく一致「京」を利用した世界初の超高解像度全球大気シミュレーショで積乱雲をリアルに表現～台風や集中豪雨などの発生メカニズムの解明に寄与～

可視化： 吉田風龍画伯 作

Geophysical Research Letter 2013年9月号 米国地球物理連合Research Spotlight
著者： 宮本佳明, 山浦剛、吉田龍二、八代尚、梶川義幸、富田浩文

水平1㎞格子以下のシミュレーション！
SNIPER project (sub km 水平格子間隔) (2012- ) :非公式プロジェクト
理研ＡＩＣＳ＋戦略分野３で共同実施

Grand Challenge like simulation on K Computer!
– 水平860m格子間隔,鉛直100層– K-Computerの1/4~全系を使わしてもらう。
– 目的• 高解像度のリファレンスランを作る。
• 計算科学的には、スケーラビリティチェックを兼ねる。

– 解析：• 全球での積雲対流の統計的描像
––––

収束性
対流の個数
隣の対流との距離などなど• 論文:
Miyamoto et al. (2013), Geophysical Research Letter--- Americal Geophysical Union EOS Research Spotlight

Why challenging? What is a challenge?
High end computerでも、sub-km GCMはデモンストレーションの域（大道芸！）
でも次の時代には、かならず、capacity computing化され、先端的成果に最も近いものになる。

A snapshot of sub-km AGCM (NICAM)
Horizontal: Δ0.87 km: vertical
100levels: integration time 24h

風龍画伯によるNICAM870mアニメーション作品

風龍：本名吉田龍二、理研ＡＩＣＳ複合系気候科学チーム所属、博士（理学）2011年彗星のごとく現れ、京を用いた計算の可視化において、数々の名作を生み出してきた。2014年学位取得後、その技にますますの磨きがかかり、業界（？）でも引く手あまた。今後の活躍（研究も？）が期待される若手のホープである。

対流の個数と対流同士距離
510
number0.25
frequency/total number(a) number of convection(b) distance between convection
D14.0D7.0
D3.5
D1.7D0.8

0.20.154100.10.05

310Δx ≧ 3.5 km:

0.8

1.7

3.5D (km)7– 個数: 解像度により4倍に増える– 対流間距離: 4 格子 => 非物理？

Δx ≦ 1.7 km:– 個数: 増加の傾向弱まる– 距離: >5格子以上にピーク=> だんだん収束始める。

1400

4 5
number of gird

10結論

対流の描像について、
Δ3.5 km  Δ1.7 km の間で変わる- リアリスティックな対流解像には、2.0km以上が、全球雲解像モデルでは必要。

2014年ハイパフォーマンスコンピューティングと計算科学シンポジウム

Efficiency of NICAM on K Computer
H. Yashiro(RIKEN/AICS)

Performance efficiency••••

ESからportingしただけ : ~4%
ステンシル計算のキャッシュ最適化 : ~5%
コードをスリム化、クリーニング :~7%
条件分岐の改良: ~10%

13%
6%

7%
6%

6%
5%
17%

8%

Physics

Dynamics

Weak scaling test•
同じ問題で、同じステップ数•

まあまあ、よいスケーラビリティ

0.9PetaFLOPS
14km

3.5km

800m

400m

3. 今回のシミュレーションの特徴と成果
今回の計算の位置づけ

今回

２００７年

（NEC イメージ図）
JAMSTECRIKEN/AICS

地球シミュレータ時代

京コンピュータ

NICAMでも3.5㎞解像度ランが限界
メモリー量とも演算性能も不足！

NICAMで1発大規模計算を行い、批判に答える！

各研究機関からの指摘：• 3.5km解像度で雲解像と言えるのか？
• 気候値が合わないのでは？
• 解像度はまだ足りない？

2011年度の8月の事例（ＭＪＯ，台風ありの）で、1㎞以下の解像度で鉛直層も100まで上げて、実験
京時代のCapablity computingの典型

ここからが真の勝負！数値的な収束の兆しは見えた！全球対流特性は把握できており、来るべきエクサスケールで、スケール相互作用の解析、衛星データ大量同化実験でいかに？

今後、、、
1. 更なる高解像度化
2. より多くのアンサンブルとデータ同化
3. より高度化された物理過程

富田の興味は上記の中でも、、、１番。

GCRM => GLES?!

なんでLESなのか？
• 大気大循環のエンジンとして :– 熱帯の雲は、熱を中緯度方面へ運ぶ重要な役割をになう。
• 熱帯の雲： 階層構造
– クラウドクラスター、スーパークラウドクラスター、熱帯低気圧、MJO• 気候学的にはエネルギー収支に大きく影
響:– 二つの効果!
• パラソル効果 :日射を跳ね返す。
• 温室効果 : 地面からの赤外放射をいったん吸収。地面へ再放出する。
– エアロゾルとの関係.
• エアロゾル間接効果: 光学的厚さを変える。

Parasol effect
Reflection ofsolar incident

Greehouse effect26
Emision of infrareエアロゾルの雲に対する効果

第一種間接効果
一部反射
寒冷化 ：
直接効果

反射大：
寒冷化度合い大反射小：
寒冷化度合い小エアロゾル：多 短寿命
 分厚い雲エアロゾル：少すかすかの雲

水蒸気凝結

長寿命
エアロゾル

第二種間接効果なんで、LESなのか？
1750年から2005年の間における、各要素の放射強制力(2007年、IPCC第4次評価報告書による)
出展： wikipedia

上層の巻雲も絡むが、境界層での細かい雲も大いに絡む

エアロゾルの雲に対する効果--- 雲の反射率の不確定性---

AGCM milestone from GCRM to GLES?!(roughly estimate)
希望的に、効率10%にしている。
解像度

一日積分計算量（FLOP数）We are here
230P
3.5㎞/L40想定マシン

効率（％）

一日積分経過
時間

一か月
積分経
過時間

何を解像？何がうれしい？

131TFLOPS
(ES2)

15％

3.2時間

4日

いくつかの積雲が集まった対流系同上。日変化などが良くなることが期待される。

1.7km/L40

1840P

10PFLOPS(K computer)10%

800m/L100

36800P

10PFLOPS(K computer)

10%

10時間

12.5日

ひつとひとつの対流のほとんどを解像

400m/L100

295000P

1EFLOPS

10%

50分

24時間

文句なしに対流解像モデル(真の雲解像)

200m/L100

2360E Exa scale era

100m/L100

18880E

50m/L200

302Z

100EFLOPS

対流を精緻に解像

Tentative goal?
10%

50分

全球LESには苦しい。
24時間

全球LES/境界層雲などを直接解像

GLESへ向かおう！
• と、言ったものの、計算量的にかなり大変！！！
当分は、、、
トップダウン型アプローチ:
全球モデリングの高解像度化と物理過程精緻化、、、準一様格子であれば、この先も続くでしょう。。。。

ボトムアップ型アプローチ
LESはほんとに機能するの？
多くの文献で言っていることは本当？
再試、追試が必要。これまでわかっていることの更なる深い理解。

• LESを使って、とりあえず今ある問題に取り組もう！
1. 格子アスペクト比に対する乱流（エネルギースペクトル）へのインパクト2. 高解像度でどの程度、浅い雲が表現できるか？遷移過程はどうか？
3. 湿潤ＬＥＳの理論は？
– などなどSCALE project• We have started the preparation towards the Global
Large Eddy Simulation• Co-design with computer science people
in RIKE/AICS
http://scale.aics.riken.jp/

Team SCALE
Computatio
nal ClimateScienceResearchTeamTL : H.Tomita

SystemSoftwareResearchTeam

TL :Y.Ishikawa

ProgrammingEnvironment ResearchTeam
TL : M.SatoHPC
Programing
Framework
ResearchTeamHPC
Usability
ResearchTeam

TL :N.Maruyama

TL : T.Maeda

Library, and Research for weather/climate researchWith computer science people (SCALE library)格子アスペクト比に対する境界層乱流への影響
S.Nishizawa(RIKEN/AICS): Preparation for submission

Filter length
effectScotti et al. (1993) formulation
dx ＼ dz

Filter length problem:• （Δｘ・Δｙ・Δｚ）^1/3: Usually used.
• This depends on the numerical scheme• In Scale-LES, 4th order advection scheme with 4thorder hyper diffusion to remove 2grid-component.
=> 2X（Δｘ・Δｙ・Δｚ）^1/3

5m

10m

10m

10mAR120m

10mAR230m

10mAR350m

5mAR1060m
100m
150m

30mAR110mAR510mAR6

5mAR2030m

30mAR210mAR1030mAR5Black : dz=10m resolution AR1
Even in AR1, energy piling becomes reduced!

In the case of
Δeq= (Δｘ・Δｙ・Δｚ）^1/3

bule : dz=10m resolution AR10If Aspect ration is large, enegy piling is terrible.

In the case ofΔeq= 2(Δｘ・Δｙ・Δｚ）^1/3

格子アスペクト比に対する境界層乱流への影響

Aspect ratio
effectS.Nishizawa(RIKEN/AICS): Preparation for submission

Scotti et al. (1993) formulation
dx ＼ dz

Aspect ratio problem:• Usually, no consideration : f =1 (left figure)• If f is appropriately set according to grid aspect ratio,
spectrum is drastically improved. ( right figure )bule : dz=10m resolution AR10
If Aspect ration is large, enegy piling is terrible.

10m

10m

10mAR120m

10mAR230m

10mAR350m

5mAR1060m
100m
150m

In the case ofΔeq= (Δｘ・Δｙ・Δｚ）^1/3

5m

30mAR110mAR510mAR65mAR2030m

30mAR210mAR1030mAR5

In the case ofΔeq= 2(Δｘ・Δｙ・Δｚ）^1/3
and f(a1,a2) setting

超高解像度で、めちゃ狭い領域ながら・・・浅い雲の再現?

• Very high resolution run for DYCOMS II RF01
– LES with 5m cubic grid and 2nd moment microphysicalscheme.

Y.Sato (RIKEN/AICS)

浅い雲の形態の遷移過程を見るための実験 Very wide domain, meso-scale LES run
 The transition between open & close cells, and POCS.
aerosol interaction & surface flux strength : key issue?

 5m v-grid & 50m h-gridwith 30km X 700km?

u=5[m/s] at initial2 km(dz =5m)

Gradually decreasing
CCN number••

Gradually increasing
surface flux

247
[km]

768 km (dx=dy=50m) 247

Y.Sato (RIKEN/AICS)
Preparation for sumition[km]m ickd28erio(pRayleigh dumping
to initial value)

••
•

•

•

Wide domain size：784 [km] x 28.8 [km] x 1.7 [km] (x, y, z)
Resolution 50m in the horizontal direction, 5m in the verticaldirection
Eviromnet wind：(u,v) = (5,0) [m/s]Integration time：12hours
Surface flux：–momentum frictional velocity=0.25 m/s–
Latent heat（93 W/m2 →200 W/m2 linearly varied
from left to right )–Sensible heat （16 W/m2 → 40 W/m2 linearly varied
from left to right）
External forcing：–largescale subusidence（1.33e-6）: Corresponding toPeru offshore–
Radiation: same as the DYCOMS-II
Aerosol ：varied from 250 [/cc] to 1 [/cc] with logarithmic

First quick view!

エアロゾルがたくさん
Surface flux 弱い

①②③

④⑤⑥⑦

エアロゾルが少ない
Surface flux：強い⑧⑨おまけ：ちょっと面白いテスト： based on Straska et al.
• 数値粘性を不安定を起こさないまで、うんと小さくしてみる– Dx=dx=1.5mまで• 通常は、200ｍ格子とかそれぐらいのテスト

解像度高めると、渦拡散が陽に見え始める（拡散方程式の形でなく） 渦粘性パラメタリ
ゼーションへ指針

Seiya NISHIZAWA（RIKEN/AICS）

まとめ• 京で行った二つのランドマーク的な研究成果を紹介– MJOの長期予報の可能性が広がる（Miyakawa et al. 2014）• 一つ一つの規模は小さいcapacity computing典型例– 全球1kmを切る大気大循環モデルによる対流の全球特性を把握（Miyamoto et al. 2013）• 大規模実験で、capability computingの典型例• いつの時代も、capability computingは、次の時代のcapacitycomputingとなることが多い。
• 次の(富田の)ターゲット、全球LESまでは、遠い。
– 今から地道に、トップダウン、ボトムアップ研究で基本的なところから入る。
– 特にボトムアップが重要と認識。
– 基礎的、基盤的研究は、なかなか日の目を見ないが、これなくして、明日はない。（あたりまえか、、^^）– 全球ＬＥＳはいつになるんだろうか？10年後？？？