ppt ファイル

SPMODEL の一環としての
球面浅水モデル
小高正嗣
石渡正樹
竹広真一
石岡圭一
林祥介
（北大・理）
（北大・地球環境）
（京大・数解研）
（京大・理）
（北大・理）
はじめに: 我々の目標

理解のための数値モデルがほしい



数式を理解するように数値モデルとその結果を
理解したい
簡単モデルと複雑モデルの結果をつなげて理解
したい
そのような数値モデルの要件



可読性が高い
拡張性が高い
階層化されたモデル群の一つ
これまでの取り組み: SPMODEL

SPMODEL(竹広他, 2002)



GFD スペクトルモデルを階層的に整備
ソースコードの可読性と拡張性を重視
特徴

スペクトル計算のための配列関数



ISPACK(石岡, 2002)の F90 インターフェース
数式との対応のよいソースコード作成を支援
データ入出力に gt4f90io ライブラリ

データは gtool4 netcdf 形式で簡単に出力
SPMODEL のこれまでの成果と課題

簡単な方程式系の場合, ソースコードの可
読性を向上

より複雑なモデルのベースとなりうるか?



GCM (大気大循環モデル)
3 次元球殻対流/MHDモデル
…
本研究の目的


SPMODEL の一環として
球面浅水方程式モデルを開発
SPMODEL の可能性を検討

数値モデルコードの可読性と拡張性
簡単モデルから複雑モデルへ発展すると?
 GCM のベースに使えるか?


数値モデルの実用性について
計算速度は実用に耐えるか?
 スパコンで使えるか?

球面浅水方程式

渦度発散型で表現した場合

1
      2 J  ,   0
t
R
D
1
      2 J  ,     2 g h  hs   E 
t
R
h
1
   h   2 J h,   0
t
R
SPMODEL の書法と配列関数

変数
xy_data
w_data

スペクトル変換・逆変換
w_xy(xy_data)
xy_w(w_data)

格子点データ
スペクトルデータ
スペクトル変換
逆変換
微分演算
xy_GradLon_w(w_data)
xy_GradLat_w(w_data)
w_Div_xy_xy(xy_data,xy_data)
w_Jacobian_w_w(w_data,w_data)
勾配型経度微分
勾配型緯度微分
発散
ヤコビアン
w_Lapla_w(w_data)
ラプラシアン
ソースコード例

リーフロッグスキーム用いた場合

1
      2 J  ,   0
t
R
w_EtaA=w_EtaB + 2*dt* &
(- w_Div_xy_xy(xy_w(w_Eta)*xy_GradLon_w(w_Chi)/R,
xy_w(w_Eta)*xy_GradLat_w(w_Chi)/R)/R
+ w_Jacobian_w_w(w_Eta, w_Psi)/R**2)
コードチェックのためのテスト計算
Williamson et al. (1992) のテスト

Case5:孤立した山を越える流れ




中緯度に山を置く
初期条件は剛体回転の地衡風と
パラメータは Williamson et al. (1992) に同じ
数値粘性なし
初期値の分布(地形, 高さ, 東西風)
その他の設定, 計算環境

分解能と時間積分法



T21, T42, T63, T106
semi-implicit リープフロッグ
+ Asselin フィルター
使用計算機環境


Intel Xeon 2.4GHz + Fujitsu F90
Fujitsu VPP800
計算結果: 相対渦度分布(T106)
計算時の CPU 時間

15 日モデル時間の計算


表示は hh:mm:ss
実行速度はまずまず
Xeon
2.4GHz
T42
T63
T106
00:00:33
00:01:09
00:13:52
00:00:09
----------
00:01:53
VPP800
まとめ
SPMODEL による球面浅水方程式モデルを開発

数値コードの可読性, 拡張性は高い



semi-implicit でも読みやすいコードが書ける
浅水多層モデルへの拡張, 強制項の導入は容易
実用に耐えうる数値モデル


計算時間はあまり犠牲になっていない
パソコンでもスパコンでも計算可能
GCM のベースとなることが期待できる
今後の課題

実用性の向上


他のスパコン環境 (SX6, SR8000) への移植
異なる計算設定を与えた場合の性能チェック


数値粘性がある場合, 強制項がある場合,…
GCM 開発へ向けた検討: 鉛直微分の扱い


差分: 可読性のよいコードが書けるか?
スペクトル: 適用可能か? 実用になるか?
参考 URL

SPMODEL


ISPACK


http://www.gfddennou.org/arch/spmodel/
http://www.gfd-dennou.org/arch/ispack/
gt4f90io

http://www.gfd-dennou.org/arch/gtool4/
付録図(1): 高さ(T106, 0, 5, 10, 15 日後)
付録図(2):初期値からの相対誤差
(左)全エネルギー (右)全エンストロフィー
付録図(3):
(左)渦度の全球積分値 (右)発散の全球積分値
ソースコード例: 時間積分(2)

リープフロッグ+台形semi-implicit の場合
1  g ht  h'
2
2
n 1


 1  g ht 2 2 h'n 1 2tF
n
1


F     h'    2 J (h' , )  hD n 1
R


1

2 
 ht      2 J  ,    E 
R


n
ソースコード例: 時間積分(2)
call get_dHsfc
w_HsfcA=((1 + dt**2*Grav*H0*n(:.1))w_HsfcB
+ 2*dt*w_dHsfc)/(1 - dt**2*Grav*H0*n(:.1))
…
end
contains
soubroutin get_dHsfc
…
w_dHsfc=-(w_Div_xy_xy(xy_Hsfc*xy_GradLon_w(w_Chi)/R,
xy_Hsfc*xy_GradLat_w(w_Chi)/R)/R
- w_Jacobian_w_w(w_Hsfc,w_Psi)/R**2)
- H0*w_DivB
-H0*dt*(w_Div_xy_xy(xy_Eta*xy_GradLon_w(w_Psi)/R,
xy_Eta*xy_GradLat_w(w_Psi)/R)/R
+ w_Jacobian_w_w(w_Eta,w_Chi)/R**2
+ w_Lapla_w(w_E)/R**2)

Download Report