計算流体工学
大気環境シミュレーションの現状
佐藤正樹
埼玉工業大学/地球フロンティア研究セン
ター
[email protected]
2004年7月6,13日 東京工業大学・総合理工学研究科
内容
•
•
•
•
大気環境シミュレーションの概観
数値計算手法
非静力学モデル:領域モデル
大気大循環モデル:全球モデル
気候シミュレーションの現状
• 気候モデリングにおけるコンピュータ環境
• 次世代大気大循環モデル
大気環境シミュレーションの概観
大気環境系の概観
• 気候系
大気、海洋、陸面(雪氷、土壌、植生)、海氷、生
態系
• 空間スケール
地球 10000km ~ 積雲 1km~乱流 10m
• 時間スケール
– 地球の歴史 数10億年~第四紀 数100万年
– 人為的気候変化 100年
– 10年スケール気候変化、エルニーニョ 3年
– 年周期、季節変化、日変化、短期予報(数時間)
気候システム
数値シミュレーションによる気候研究
大気の空間スケール
積雲クラスタ~ 100km
積雲
10km
12740km
10km
1km
気候変動の時間スケール
IPCC2001: http://www.ipcc.ch
大気の数値モデル
• メソモデル:雲解像モデル
100km×100km:非静力学モデル
• 領域モデル
1000km×1000km 領域:静力/非静力学モデル
• 大気大循環モデル
全球3次元: 静力学モデル
数値モデルの格子数と分解能
• 格子数
数100×数100×数10=数105
⇒ 1000×1000×100=108 (次世代:ES)
• 領域幅と格子間隔:
大循環 10000km:Δx=100km ⇒ 10km
領域 1000km:Δx=10km ⇒ 1km
メソ
100km:Δx=1km
•
–
–
–
•
–
–
–
–
大気モデルの構成
(領域・メソモデル・大気大循環モデ
ル)
力学過程
方程式系
空間差分法,時間差分法
数値拡散
物理過程
水循環過程(雲物理,積雲パラメタリゼーション)
放射過程(太陽放射,赤外放射)
乱流過程(接地境界層,subgrid乱流)
化学変化・物質輸送(オゾン,エアロゾル)
数値計算手法
数値計算手法
•方程式系
•大気の基本バランスと波動
•CFL条件
•音波・重力波の扱い
•移流計算
•擬スペクトル法・変換法
流体を記述する方程式系
•非圧縮ブシネスク系:工学or理論
•完全圧縮性成層流体
非静力学方程式
•プリミティブ方程式系
静力学近似
•浅水波方程式系
•準地衡風方程式系
大気の基本バランス
• 静力学平衡
1 
0
pg
 z
• 地衡風平衡
1 
fv  
pg
 x
1 
 fu  
pg
 y
f  2 sin 
工学における流体の方程式系(1)
非圧縮の式
v
1
 v  v  
p  F
t
0
運動方程式
v  0
連続の式
工学における流体の方程式系(2)
ブシネスク方程式系
v
1
'
 v  v  
p 
g F
t
0
0
運動方程式
v  0
連続の式
T
Q
 v  T 
t
Cp
熱エネルギーの式
    T  T0 
状態方程式
大気の方程式系
圧縮性成層流体の方程式系
v
1
 v  v   p  g  F
t


   v   0
t
T
1 dp
Q
 v  T 

t
C p  dt C p 
p  Rd T
運動方程式
連続の式
熱エネルギーの式
状態方程式
大気の方程式系
回転系の圧縮性成層流体の方程式系
v
1
 v  v  2Ω v   p  g  F 運動方程式
t


連続の式
   v   0
t
T
1 dp
Q
熱エネルギーの式
 v  T 

t
C p  dt C p 
p  Rd T
状態方程式
波動
• 音波:圧縮性
音速
cs  Rd T  330m/s
• 重力波:成層効果:浮力
Brunt Vaisala 振動数
g  T g 
1
N

 10 min
T  z C p 
• ロスビー波:緯度方向の回転角速度の変化
df


 2 cos ; c  2
 U  10m/s
2
dy
kx  k y
音波
線形化、1次元、断熱
u
1 


p 
t
 0 x

 

   0 u   0
t x


1 p

 2

t cs t

1  p  p
 2 0
2
2
x
cs t
2
2
p  p0 expikx  it 
   cs k
2
2
2
分散関係
重力波
ブシネスク、線形、2次元、断熱
u
1 

p
t
 0 x

 2  2 2 
dTs  2
w0
 2  2  2 w  g
2
dz x
w
1 
 '  t  x z 

p
g
t
 0 z
0 

v  0
 w  w expikx  im z  it 
0

dTs
T
2

w
0
dTs
k
2
t
dz
g
   2
2
k m
dz


  T  T0 

分散関係
0

CFL(Courant-Friedrichs-Lewy)条件
移流方程式
q

u
q0
t
x
u  0
差分化, von Neumann 解析
q
q
n
j
n 1
q
t
j
n
j
q j q
u
x
n
n
j 1
 0 (upstream)
  expij ,  kx
n

   1  1  e i

t
 0    1 for   1  ;  u
(Courantnumber)
x
大気モデルの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)
• ロスビー波:移流と同等
音波の差分解法
 u n 1 j  u n j
1 p n j  p n j 1


t
0
x

 n 1
n
n
n
p

p
u

u
j
j
j 1
j
2

   0 cs

t
x

n: Forward (explicit, Euler)
n+1: Backward (implicit)
pn+1 ,un :Forward-Backward (semi-implicit)
固有値λ解析:
u n j  U n expij ,
pn j
 P n expi j  1 / 2 
0
AX n1  BX n A  B  0   2cs t sin  : Coulant number
Forward
x
2
 i / cs  U n 
U n 1   1
 n 1   
 n    1 i ,   1: unstable

1  P 
 P   ics
Backward
n 1
n
1
1
i

/
c




U 
s U

,   1: stable

ic
  n 1   n 
1  i
1
P
P
 s

  
Forward-Backward
2
2
  1   i 1  ,
n 1
n
 1 0 U  1  i / cs  U 
2
4
 n
ic 1  n1   0

1   P    1 : st able
 s
 P  
非静力学モデル
(メソモデル・積雲解像モデル)
代表的な非静力学モデル
•
•
•
•
•
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法)、非弾性系
 v n 1  v n
1
n
n 1



v


v



p
 g F

s
 t

n 1
n
1
p

p
n 1

    s v
2
 cs
t



cs→∞ とすれば非弾性系. ρsは z だけの関数とする.
n
 1

p
1
2
n 1
n
 2
p  2






v
s
2
 c t 2


t
cs t
 s


   v  v    g   F 
n
s
s
s
ポアソン方程式

ポアソン方程式の解法
 2( x, y, z, t )  S ( x, y, z, t )
  
 2  2 S
2
x
y
z
 i 1, j ,k  2 i , j ,k   i 1, j ,k  i. j 1.k  2 i , j ,k   i , j 1,k

2
x
y 2
⇒
 i , j ,k 1  2 i , j ,k   i , j ,k 1

 Si , j ,k
2
z
2
2
•
•
•
•
2
緩和法:Jacobi法, Gauss-Seidel法, SOR法
共役勾配法
スペクトル法
マルチグリッド法:並列計算
ポアソン方程式の解法
•
•
•
•
緩和法:Jacobi法, Gauss-Seidel法, SOR法
共役勾配法
スペクトル法
マルチグリッド法
並列計算
水平陽鉛直陰解法(HEVI法)
 v H n 1  v H n
1
n
n
n 1
n 1



v


v



p

F
:

u
,
v

H
H
H

t

s

 wn 1  wn

1  n 1
n
 v  w  
p  g  Fz 

 s z

n 1
n 1
 t
:

w
,
p

 1 p n 1  p n

n 1
 2
  H   s v H

 s wn1 

t
z
 cs







 1
 2  n 1
pn
1
1 
n 1
n
 2


p





v


w
H
s H
s
2
2
 c t 2 z 2 

t

t

z
cs t
 s

   v  w   g   F 
n
s
s
s
z
鉛直方向の3重対角行列


Splitting method
• 音波・重力波に関係する項を短い時間間隔Δτで計算
• 他の項を長い時間隔Δtで計算
• divergence damping 等数値拡散が必要
 v H     v H 
1
t
t




v


v



p

F

H
H
H



s

 w    w
1    
t
t



v


w


p
 g  Fz


 s z


  

1
p

p

  
  





v


w
H
s H
s
 cs 2

z






非静力学モデルの時間ステップ
• 音波の扱いが問題
• 陰解法(HIVI法):時間ステップの制約がない
• HEVI法:水平伝播する音波により時間ステップが
制約される
⇒ splitting 法により、他の項を長い時間ステップで
計算
• 次に問題となるのは移流
音波 cs=300m/s, 風速 U=数10m/s
移流スキーム
• Euler法,フラックス型
有限体積法 ⇒ 保存系
• semi-Lagrange 法
一般には保存は満たされない
CIP
• 保存型 semi-Lagrange 法
Godnov, van-Leer,PPM,Lin and Rood(1996)
移流スキームの性質
•
•
•
•
diffusion:高精度スキーム
dispersion: monotonicity(shape-preserving)
一般に高精度であればより振動的
高精度(2次、3次)かつ limiter により単調
性をできるだけ確保
移流方程式
dq
0
保存則
dt
 
 u   0 連続の式
t x

 q
:移流型
 t  u x q  0

 q    qu  0 :フラックス型
 t
x
Euler 法
•ρq の領域積分が保存する
•セル境界におけるq の評価法により精度,単調性が依存
•Courant数<1 を要求する
  q  
  qu   0
t
x
q 
 q  j

t
f  qu
n 1
n
j
f n j 1  f n j1
2
2
x
0
Semi-Lagrange 法
• q が Lagrange 的に保存
•セル内のqの補間方法により単調性,精度が依存
•Courant数に制限がない (but Δt・du/dx<1:Lipschitz条件)
• ρq の領域積分は保存しない
dq
0
dt
q( x, t )  q( x  ut , t  t )
q n 1 j  q n d . p.
保存型 Semi-Lagrange 法
•ρq の領域積分が保存する
• q が近似的に Lagrange 的に保存
•セル内のqの補間方法により単調性,精度が依存
•Coulant数に制限がない (but Δt・du/dx<1:Lipschitz 条件)
q 
 q  j
2

t
x
f  q( x  ut , t  t )u
n 1
n
j
f n j 1  f n j1
2
0
物理過程
•
–
–
–
–
–
物理過程
水循環過程(雲物理,積雲パラメタリゼーション)
放射過程(太陽放射,赤外放射)
乱流過程(接地境界層,subgrid乱流)
重力波抵抗
化学変化・物質輸送(オゾン,エアロゾル)
物理過程のスケール依存
•
–
–
•
–
–
雲物理過程
メソモデル:微物理過程:雨滴・雪霰の生成・成長過程
大循環モデル:積雲パラメタリゼーション
乱流過程
分解能数10mまで(1km以下でも使う): LES
数km以上:統計乱流モデル
Moller and Yamada,1次,2次,3次クロージャ
例:積雲のパラメタリゼーション
• 積雲のスケール1km vs 大循環モデルの分解能100km
• 分解能以下の現象をパラメタリゼーション
Arakawa and Schubert:格子内で準平衡を過程
• 熱の式、水の変化式にソース項として考慮
• パラメタリゼーションの方式は分解能に応じて変更
• パラメタリゼーションは準理論、経験に基づくものであり、
不確定性が大きい
1km
100km
大気大循環モデル
大気大循環モデル(AGCM)
• 静力学近似に従うプリミティブ方程式系
• 球面調和関数による擬スペクトル変換
– 毎ステップに格子空間とスペクトル空間の
ルジャンドル変換を行う
– 微分の計算を格子空間で行う
• 重力波・音波をsemi-implicit法で計算
• 渦度,発散を予報し,ポアソン方程式を解
いて速度を求める
AGCMの方程式系
経度方向運動方程式
緯度方向運動方程式
静力学平衡
連続の式
熱エネルギーの式
u
1 1 p
  v  u  fv 
 Fx
t
 a 
v
1 1 p
  v  v  fu 
 Fy
t
 a cos 
1 p
0
g
 z

    v   0
t
T
1 dp Q
  v  T 

t
C p  dt C p
状態方程式
p  Rd T
水蒸気の式
q
  v  q  S q
t
方程式の変形
• 鉛直座標: σ=p/ps =圧力/地表面圧力
ジオポテンシャル Φ=gz
Rd T


⇒静力学平衡
z

• 渦度・発散
v u
1
   
x y
a cos
u v
1
  
x y
a cos
v
1 u cos 

 a cos

u
1  v cos 

 a cos

渦度・発散方程式
流線関数ψ
速度ポテンシャルχ


v u

2
u   y , v  x ,   x  y   H 


u   , v   ,   u  v   2 
H


x

y

x

y

1 p
 u


v


u

fv


 0 x
水平方向運動方程式  t

 v   v  v  fu  1 p
 t
 0 y
 
  H    f v

渦度方程式
 t
⇒

2




p
v
2
H

発散方程式

 k   H    f v   H  

 t

2
 0

ζ,δを予報⇒ψ,χについて解く⇒u,vを求める
スペクトル法とルジャンドル変換
球面調和関数Y,ルジャンドル陪関数P
 1
2
1
 
  m
n(n  1) m


 H Yn   2 2
 2
cos Yn   2 Yn
2

 
a
 a cos   a cos  
2
m
Yn  ,    Pn sin  eim
m
m
球面上のポアソン方程式
H X  s
m
m
m
m
X   Xˆ n Yn , s   sˆn Yn
2
n ,m
n ,m
n(n  1) m
m
ˆ
 Xn  
sˆn
2
a
スペクトル法の切断波数
三角形切断: 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
σ座標渦度発散系の方程式系


 v H       v H 
t

Rd Tv





1 Av
1  ( Au cos )


 D
t a cos  a cos

2



1 Au
1  ( Av cos )
v
2
H 



  H    Rd T 
 D

t a cos  a cos

2 

T
1 uT '
1  (vT ' cos )
T
d Q


 T '   
 T

 DT
t
a cos  a cos


dt C p
q
1 uq
1  (vq cos )
q


 q  
 S q  Dq
t
a cos  a cos


AGCM実験例:標準比較実験
• Dynamical Core: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
• NCAR
semi-Lagrange
T42L49
Aqua Planet 実験
Neale and Hoskins
海面水温分布
降水量
Model Version
AMIP
Atmospheric
Model
Intercomparison
Project
BMRC BMRC2.3 (R31 L9) 1990
BMRC BMRC3.7 (R31 L17) 1995
BMRC BMRC3.7.1 (R31 L17) 1995
CCC GCMII (T32 L10) 1990
CCSR/NIES AGCM (T21 L20) 1995
CNRM EMERAUDE (T42 L30) 1992
CNRM ARPEGE Cy11 (T42 L30) 1995
COLA COLA1.1 (R40 L18) 1993
CSIRO CSIRO9 Mark1 (R21 L9) 1992
CSU CSU91 (4x5 L17) 1991
CSU CSU95 (4x5 L17) 1995
DERF GFDLSM392.2 (T42 L18) 1993
DERF GFDLSM195(T42 L18) 1995
DNM A5407.V1 (4x5 L7) 1991
DNM A5407.V2 (4x5 L7) 1995
ECMWF ECMWF Cy36 (T42 L19) 1990
GFDL CDG1 (R30 L14) 1992
GISS Model II Prime (4x5 L9) 1994
GLA GCM-01.0 AMIP-01 (4x5 L17) 1992
GSFC GEOS-1 (4x5 L20) 1993
IAP IAP-2L (4x5 L2) 1993
JMA GSM8911 (T42 L21) 1993
LLNL LLNL/UCLA MPP1 (4x5 L15) 1995
LMD LMD5 (3.6x5.6 L11) 1991
LMD LMD6b (3.6x5.6 L11) 1995
LMD LMD6s (3.6x5.6 L11) 1995
MGO AMIP92 (T30 L14) 1992
MPI ECHAM3 (T42 L19) 1992
MPI ECHAM4 (T42 L19) 1995
MRI GCM-II (4x5 L15) 1993
MRI GCM-IIb (4x5 L15) 1995
NCAR CCM2 (T42 L18) 1992
NMC MRF (T40 L18) 1992
NRL NOGAPS3.2 (T47 L18) 1993
NRL NOGAPS3.4 (T47 L18) 1995
NTU GCM (T42 L13) 1995
RPN NWP-D40P29 (T63 L23) 1993
SUNYA CCM1-TG (R15 L12) 1990
SUNYA/NCAR GENESIS1.5 (T31 L18) 1994
SUNYA/NCAR GENESIS1.5A (T31 L18) 1995
UCLA AGCM6.4 (4x5 L15) 1992
UGAMP UGCM1.3 (T42 L19) 1993
UIUC MLAM-AMIP (4x5 L7) 1993
UKMO HADAM1 (2.5x3.75 L19) 1993
YONU Tr5.1 (4x5 L5) 1994
YONU Tr7.1 (4x5 L7) 1995
Horizontal
Representation
Resolution
spectral
rhomboidal 31
spectral
rhomboidal 31
spectral
rhomboidal 31
spectral
T32
spectral
T21
spectral
T42
spectral
T42
spectral
rhomboidal 40
spectral
rhomboidal 21
finite difference
4 x 5 degrees
finite difference
4 x 5 degrees
spectral
T42
spectral
T42
finite difference
4 x 5 degrees
finite difference
4 x 5 degrees
spectral
T42
spectral
rhomboidal 30
finite difference
4 x 5 degrees
finite difference
4 x 5 degrees
finite difference
4 x 5 degrees
finite difference
4 x 5 degrees
spectral
T42
finite difference
4 x 5 degrees
finite difference
50 sinlat x 64 lon
finite difference
50 sinlat x 64 lon
finite difference
50 sinlat x 64 lon
spectral
T30
spectral
T42
spectral
T42
finite difference
4 x 5 degrees
finite difference
4 x 5 degrees
spectral
T42
spectral
T40
spectral
T47
spectral
T47
spectral
T42
spectral semi-Lagrangian T63
spectral
rhomboidal 15
spectral
T31
spectral
T31
finite difference
4 x 5 degrees
spectral
T42
finite difference
4 x 5 degrees
finite difference
2.5 x 3.75 degrees
finite difference
4 x 5 degrees
finite difference
4 x 5 degrees
Coordinates
sigma
sigma
sigma
hybrid
sigma
hybrid
hybrid
sigma
sigma
modified sigma
modified sigma
sigma
sigma
sigma
sigma
hybrid
sigma
sigma
sigma
sigma
modified sigma
hybrid
modified sigma
sigma
sigma
sigma
sigma
hybrid
hybrid
hybrid
hybrid
hybrid
sigma
hybrid
hybrid
sigma
sigma
sigma
hybrid/sigma
hybrid/sigma
modified sigma
hybrid
sigma
hybrid
modified sigma
modified sigma
Vertical
No. Levels
9 (3, 3)
17 (5, 5)
17 (5, 5)
10 (3, 4)
20 (5, 8)
30 (4, 20)
30 (4, 20)
18 (5, 4)
9 (3, 3)
17 (2, 6)
17 (2, 6)
18 (5, 5)
18 (5, 5)
7 (1, 1)
7 (1, 1)
19 (5, 7)
14 (4, 4)
9 (2, 2)
17 (5, 4)
20 (5, 7)
2 (0, 0)
21 (6, 7)
15 (2, 9)
11 (3, 2)
11 (3, 2)
11 (3, 2)
14 (5, 4)
19 (5, 7)
19 (5, 7)
15 (1, 9)
15 (1, 9)
18 (4, 7)
18 (5, 4)
18 (5, 5)
18 (5, 5)
13 (3, 4)
23 (7, 7)
12 (3, 5)
18 (4, 7)
18 (4, 7)
15 (2, 9)
19 (5, 7)
7 (3, 0)
19 (4, 7)
5 (1, 1)
7 (3, 1)
Bottom, Top
991, 9 hPa
991, 9 hPa
991, 9 hPa
980, 5 hPa
995 hPa
995, 0.01 hPa
995, 0.01 hPa
995, 10 hPa
979, 21 hPa
variable, 51 hPa
variable, 51 hPa
998, 2 hPa
998, 2 hPa
929, 71 hPa
929, 71 hPa
996, 10 hPa
997, 15 hPa
975, 10 hPa
994, 12 hPa
994, 10 hPa
800, 200 hPa
995, 10 hPa
variable, 1 hPa
979, 4 hPa
979, 4 hPa
979, 4 hPa
992, 13 hPa
996, 10 hPa
996, 10 hPa
variable, 1 hPa
variable, 1 hPa
992, 3 hPa
995, 21 hPa
995, 1 hPa
995, 1 hPa
962, 1 hPa
1000, 10 hPa
991, 9 hPa
993, 5 hPa
993, 5 hPa
variable, 1 hPa
996, 10 hPa
990, 200 hPa
997, 5 hPa
900, 100 hPa
990, 100 hPa
気候モデリングにおける
コンピュータ環境
計算機環境
多様な計算機環境
• スカラー vs ベクトル
• 単一PE vs 並列
• PC,WS vs 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)
• 領域幅と格子間隔:
大循環 10000km:Δx=100km ⇒ 10km
領域 1000km:Δx=10km ⇒ 1km
メソ
100km:Δx=1km
• 次世代(今後10年以内)に地球シミュレー
タ(ES)など超分散並列コンピュータによっ
て従来より一桁上の分解能の計算が可能
になる。
• 大気モデルは、統合化に向かう。
– 大循環モデルは格子間隔10km
– 領域モデルは格子間隔1km
• 積雲を分解するのに必要な格子間隔:1km
• 次世代大循環モデルでも分解能は不十分
地球フロンティアにおける
次世代大気大循環モデル
次世代大気大循環モデル
• 全球雲解像モデル:NICAM:Nonhydrostatic Icosahedral
Atmospheric Model (非静力学正20面体大気モデル)
• 地球フロンティア研究センター・地球環境モデリング研究
プログラムで開発中
地球シミュレータの性能を最大限生かしたモデル
• 高分解能全球3次元気候モデル:長期気候予測
• 準一様な格子モデル:格子間隔~10km
• 非静力学方程式系
• 積雲をできるだけ顕わに表現し,統計平衡を仮定した積
雲パラメタリゼーションを用いない
次世代大気大循環モデルの開発
•
–
–
–
–
開発の手順
準一様格子の浅水波モデル(1層モデル)
保存則を満たす非静力学方程式モデル
物理過程の検討(雲物理,放射,乱流)
計算性能(ベクトル化,並列化,データ処
理,画像解析)
格子構造
正20面体モデル
等角立方格子モデル
正20面体格子の生成
glevel-1
glevel-3
glevel-2
glevel-4
Number
of grid
points of
n-level:
  2
10 2
n 2
level 7
level 8
(~56km)
(~28km)
level 9
level 10
(~14km)
(~7km)
等角立方格子の生成
Number of
grid points:
6N 2
10x10
20x20
40x40
80x80
格子間隔の比較
浅水波モデルの結果
Williamson Test Case 5
Conformal cubic grid Spectrum model
(T213)
Icosahedral grid
非静力学コア
質量・エネルギー保存スキーム
Cold bubble 実験
エネルギーの時間変化
Squall line
(GCSS CASE1)
降水量の時間変化
降水・質量変化
雲水量と雨
Life Cycle of Extratropical Cyclone Exp.(1)
Polvani and Scott(2002)
glevel 6
Δx=112km
glevel 8
Δx=28km
glevel 10
Δx=7km
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
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

AFES ( green line)
Elapse time increases
in the sense O(n3).
 Legendre transformation

NICAM ( red line)
Elapse time increases
in the sense O(n2).

In all resolutions,
NICAM is faster than
AFES!
Caution: what is the resolvable scale? 2Δx or 4Δx?
スペクトルモデルとの比較
 Maximum time step and one-day simulation time
If L = 2 πR / N = 4 Δx
NICAM
t
1day
time
AFES
t
1day
time
gl7
450
6.70
gl8
225
32.1
gl9
113
210
gl10
57
1519
gl11
29
12200
T159
400
8.02
T319
200
27.9
T639
100
184
T1279
50
1884
T2559
25
24930
• Time step of NICAM can be larger than AFES
At T1279 & glevel-10, NICAM is faster than AFES.
開発の現状
• 現状
– 保存型方程式に基づく正二十面体格子3次元球面
非静力学モデルの力学コアを完成
– 全球を水平分解能 3.5km で分割した実験を実施
– 水惑星実験
– 海陸配置の導入、地形
– 雲物理過程に対する依存性
– エアロゾル
まとめ
• 気候モデルは拡張しつづけている.
• 大気モデルは気候モデルのごく一部分である.
• 代表的な大気モデルは,領域モデル, 大気大循環モデル
である.領域モデルはメソモデル、非静力学モデル,雲解
像モデルを包含する総称となりつつある.
• 大気モデルは,力学過程,物理過程からなる.
• 領域モデルはフリーで利用できるものが数多くある.
• 大気大循環モデルも利用可能なものはあるが、物理過程
などの不確かさを十分理解して使う必要がある.
• 数km格子になるまで,積雲パラメタリゼーションに伴う不
確定性が避けられない.
• 複雑なモデルの評価のために標準実験が提唱さ
れている.
• コンピュータの進歩により大気モデルの領域依
存性は統合化されるであろう.
• 地球フロンティアで開発中の次世代大気大循環
モデルは,3.5km格子で全球を覆う非静力学大
循環モデルである.
• 近年のモデル開発は,数理解析(差分スキーム),
理論(大気科学),素過程(雲物理,放射,乱流),
計算科学,システム開発などの総合的な分野の
協力が必要である。
ダウンロード

東京工業大学大学院総合理工学研究科「計算流体工学」