森川 靖大

北海道大学大学院理学研究科 地球惑星科学専攻

地球流体力学研究室 博士後期課程 1 年

惑星大気大循環モデル
DCPAM の開発：
データ I/O ライブラリの改良および
力学コアの設計と実装

目次

はじめに

目標

モデルの解説

動作試験

はじめに

「大気大循環モデル」とは

別名：GCM ; General Circulation Model

全球規模の大気循環を表現する数値モデル

気象予測 (天気予報)・気候変動予測 (温暖化等) で活躍

モデルの高度化・複雑化

3 次元の力学過程・物理過程の水平・鉛直離散化

再現性・高速化重視による特殊なチューニング

目標

読みやすい・書き換えやすい GCM の開発

結果の考察を容易にする GCM

系を簡素化した縮小モデルを容易に作成

地球に限らず他の惑星でも使えるGCM

物理過程と惑星パラメータを容易に交換

類似する思想のモデル

FMS (GFDL, 2005)

http://www.gfdl.noaa.gov/~fms/

DCPAM
(Dennou Club Planetary Atmospheric Model)

AGCM5 (SWAMP Project, 1998) を参考に

大学等での研究・教育用に整理された、比較的チューニングされていない GCM

自由に利用可能 (Web で公開)

ドキュメントが整備

座標系・支配方程式系・離散化を踏襲

Fortran90/95 の新機能を活用

GCM の可読性と可変性を高める

モデルの解説 (1)

座標系

支配方程式系

離散化

座標系

水平方向：緯度経度座標系

鉛直方向：座標系

= p/p0 (p0：地表面気圧)

=1

等  面

等 p 面

地形

座標と p 座標

緯度・経度座標

支配方程式系

プリミティブ方程式系

5つの予報変数を解く：渦度、発散D、地表面気圧p0 、温度T、比湿q

境界条件： at  = 0, 1

D ： 水平拡散

Sq ： 水蒸気のソース

Q ： 加熱

運動方程式 (渦度・発散形)

連続の式

静水圧の式

水蒸気の式

熱の式

離散化

水平：スペクトル法 (球面調和関数展開)

一般的に「差分法」に比べて精度が高い

「球面を覆う」大気を扱うのに都合が良い

鉛直：Arakawa and Suarez (1983)

全エネルギー・質量・角運動量・温位を保存

時間：陽的 leapfrog スキーム

離散化が単純

計算速度の向上の観点からは、重力波に対して陰的スキームを用いた 「semi-implicit」 の方が良い

モデルの解説 (２)

モデル内部の階層化

変数命名規則

数式のコードの簡素化

ドキュメント生成システム

データ I/O 部分の簡素化

モデル内部の階層化(1)

F90 の「モジュール機能」を活用

「モジュール」の交換・分離を容易に ⇒ 可変性を向上

モデル内部の階層化(2)

初期化処理の統一 (The FMS Manual (Balaji, 2002) より)

各モジュールで、必ず初期化用サブルーチンを用意

サブルーチンの命名法も統一

個々のモジュールが間接的に依存するモジュールの情報を知る必要無し ⇒ 可変性向上へ

ファイル名 ： module.f90,

例： dynamics.f90

モジュール名 ： module_mod,

例： dynamics_mod

サブルーチン名： module_init,

例： dynamics_init

モデルの解説 (２)

モデル内部の階層化

変数命名規則

数式のコードの簡素化

ドキュメント生成システム

データ I/O 部分の簡素化

変数命名規則

変数の情報が容易に分かる命名規則

SPMODEL (竹広 他, 2004) の命名規則の拡張

変数名規則 ： x_var_t

x_ ： 空間次元

var ： 物理的意味

_t ： 時刻情報

具体例

x_Lon ： 空間次元経度方向の座標データ

xyz_Temp_b ： 温度の3次元格子点データ (時刻 t-t)

wz_Vor_a ： 渦度の水平スペクトルデータ

＋ 鉛直格子点データ (時刻 t+t)

}

一目でわかる

⇒ 可読性が向上

モデルの解説 (２)

モデル内部の階層化

変数命名規則

数式のコードの簡素化

ドキュメント生成システム

データ I/O 部分の簡素化

数式のコードの簡素化 (1)

FORmula TRANslater ではあるが・・・

3次元・球座標・離散化によって方程式が十分複雑に

個々の項をサブルーチンで演算すると、プログラム構造が整理される一方、「数式」とはかけ離れた形に ⇒ 可読性が低下

関数を利用 (SPMODEL (竹広 他, 2004) より)

F90 で提供される「配列を返す関数」を活用

SPMODEL で提供されるスペクトル演算用の関数を利用

「数式」に近い形で記述 ⇒ 可読性を向上

数式のコードの簡素化 (2)

SPMODEL と変数命名規則を併用

具体例

渦度の式

wz_Vor_t = & ! 運動量移流項 UA, VA を渦度の

& ! 時間変化のスペクトルデータに

& wa_Div_xya_xya( xyz_VA , - xyz_UA ) / Rplanet

wz_Vor_a = & ! t+Δt の渦度のスペクトルデータを計算

& wa_xya( xyz_Vor_b ) + 2.0d0*DelTime * wz_Vor_t

xyz_Vor_a = & ! 渦度のスペクトルデータを格子点データに

& xya_wa( wz_Vor_a )

 ： 渦度

t ： 時間

a ： 惑星半径

 ： 経度

 ： sin 

 ： 緯度

UA ： 東西運動量輸送

VA ： 南北運動量輸送

数式のコードの簡素化 (3)

F77 での具体例

渦度の式

CALL SMTG2S !" t-Δt の格子点データをスペクトルに

I ( MMAX , IMAX, ...

I GBVOR, !" 渦度の格子点データ (t-Δt)

O WBVOR, ... ) !" 渦度のスペクトルデータ (t-Δt)

CALL SMTV2S !" 運動量移流項 UA, VA を渦度の

!" 時間変化のスペクトルデータに

I ( MMAX , IMAX , ...

I GTUA , GTVA , ... !" 運動量移流項 UA, VA

O WTVOR, ... ) !" 渦度の時間変化

!" t+Δt の渦度のスペクトルデータを計算

DO 7100 K = 1 , KMAX !" 鉛直方向のループ

DO 7100 NM = 1 , NMDIM !" 水平方向のループ

WAVOR( NM,K ) = WBVOR( NM,K ) + WTVOR( NM,K )*DELT2

7100 CONTINUE

CALL SMTS2G !" 時間変化項をスペクトルに

I ( MMAX, IMAX, ...

I WAVOR, !" 渦度のスペクトルデータ (t+Δt)

O GAVOR, ... ) !" 渦度の格子点データ (t+Δt)

モデルの解説 (２)

モデル内部の階層化

変数命名規則

数式のコードの簡素化

ドキュメント生成システム

データ I/O 部分の簡素化

ドキュメント生成システム(1)

「ドキュメント」の種類 (方言あり)

数学的・技術的記述を含む文書

数値モデルの場合、支配方程式系や離散化などが記載

チュートリアル

モデルの使用法を習得するための学習用手引き

リファレンスマニュアル

個々のプログラムの動作の要約や引数の情報、その他プログラムを利用する際に必要となる仕様を記載

プログラムの内部の理解、プログラムの一部の交換・分離を補助 ⇒ 可変性・可読性の向上

ドキュメント生成システム(2)

リファレンスマニュアルの手動作成

面倒くさい

プログラムとマニュアルで同じことを２度書くことが多い

プログラムの更新の度に書き換えねばならない

プログラムとマニュアルの乖離が多発

自動生成の場合

F90 ソースとしての可読性を低下させない

原理的には、HTML でも XML でもなんでも、埋め込んで自動生成できるが…

引数情報など、ソースコードをそのまま出力

２度書きしたくない

ドキュメント生成システム(3)

F90 ソースにRD (Ruby Documentation Project, 2005)を埋め込み

:

contains

!=begin

!== Procedure Interface

!=== Subroutine module_name_init : モジュールの初期化

!NAMELIST を入力し、グローバル変数を allocate する。

!

subroutine module_name_init(inchar, outint, inoutdata, inoutdb)

!==== Dependency

use type_mod , only : REKIND, DBKIND, INTKIND, STRING

!=end

implicit none

!=begin

!==== NAMELIST

character(TOKEN) :: name = '' ! 名前

integer(INTKIND) :: length = 0 ! 長さ

namelist /module_name_nml/ &

& name , & ! 名前

& length , & ! 長さ

& Loop ! ループの初期値

!==== Input

character(*) , intent(in) :: inchar ! Input Character

!==== Output

integer(INTKIND), intent(out) :: outint ! Output Integer

!==== In/Out

real(REKIND), intent(inout) :: inoutdata ! In/Out Data

real(DBKIND), pointer,optional :: inoutdb ! In/Out Data

!=end

continue

:

end subroutine module_name_init

:

end module sample_mod

!=begin

!= Module module_name_mod : Sample module of RD in F90 source code.

!

! * Developers: Morikawa Yasuhiro

! * Version: $Id: module_name.f90,v 1.4 morikawa $

!

!== Overview

!module_name_mod の概要

! :

!== Error Handling

!module_name_mod におけるエラー処理の詳細。

! :

!=end

module module_name_mod

!=begin

!== Dependency

use type_mod, only: STRING, REKIND, DBKIND, INTKIND

!=end

implicit none

!=begin

!== Public Interface

private

public :: module_name_init, module_name_end ! subroutines

public :: module_name_func ! functions

public :: data1, data2 ! variables

!== Public Data

real(DBKIND), save :: & ! follow data is default values.

&

& data1 = 3.141592653589793 , & ! 円周率

& data2 = 6.371d6 ! 球の半径

:

ドキュメント生成システム(4)

HTML の自動生成

モデルの解説 (２)

モデル内部の階層化

変数命名規則

数式のコードの簡素化

ドキュメント生成システム

データ I/O 部分の簡素化

データ I/O 部分の簡素化(1)

データI/O

Input/Output つまりデータの入力・出力

ソースコードを繁雑に or プログラムの分離・交換を困難に ⇒ 可読性・可変性低下

I/O 部分を簡素化するために

I/Oライブラリに機能を押し付ける

演算部分からは簡素化されたサブルーチンを利用

⇒ 可読性・可変性向上

データ I/O 部分の簡素化(2)

データ I/O ライブラリ gt4f90io

gtool4 Fortran90 ツール・ライブラリ (豊田 他, 2002) の I/O 部分の独立・改良

多次元格子点データのための自己記述的表現方法 gtool4 netCDF 規約 (豊田 他, 2000) に基づくデータの I/O

データの参照効率の向上

I/O コード部分の簡素化を実現

最低限必要は4つ

HistoryCreate : 初期設定 (出力ファイル等々)

HistoryAddvariable : 変数定義

HistoryPut : 変数出力

HistoryClose : 終了処理

データ I/O 部分の簡素化(3)

サンプル

program sample

use gt4_history ! モジュールの使用を宣言

[型宣言] ......

call HistoryCreate( & ! ヒストリー作成

file='sample.nc', title='gt4_history', & ! ・ファイル名、タイトル

..., dims=(/'x','t'/), dimsizes=(/30,0/), & ! ・次元変数、次元サイズ

.......)

call HistoryAddVariable( & ! 変数定義

varname='temp', dims=(/'x','t'/), .... ) ! ・変数名、依存次元、..

[時間積分ループ]

:

call HistoryPut(varname= 'temp', value=temp) ! 変数の出力

:

[時間積分ループ 終わり]

call HistoryClose ! 終了の処理

stop

end program sample

動作試験

Held and Suarez

(1994) テスト

計算結果

Held and Suarez (1994) テスト

与える放射平衡温度

(Held and Suarez, 1994 より)

Held and Suarez(1994) の T63L20 のスペクトルモデルから得られた

帯状平均東西風。200～1200日の平均値 (Held and Suarez, 1994 より)

T63L20

GCM の力学過程の長期的性質を評価

気候モデルとしての GCM の評価

計算設定

左下図のニュートン冷却と  >0.7 の風にレイリー摩擦を与える

1200 日積分を行い、その統計量を見る

計算結果

T21L20, T42L20 で実験

T42L20 のパターンは Held and Suarez の結果と一致

左が T21L20, 右が T42L20 のスペクトルモデルから得られた、帯状平均した東西風 (単位 m/s) 。

T21L20 は200～1200日の平均値。 T42L20 は 200～1250日の平均値

T42L20

T21L20

まとめ

GCM の可読性・可変性が向上

モジュール機能を活かしたプログラムの階層化

次元・物理的意味・時間が分かる変数命名規則

SPMODEL による数式のコードの簡素化

リファレンスマニュアルの自動生成

データI/O部分を GCM の演算部分と分離、簡素なサブルーチンが利用できるよう整備

力学コアの長期的性質の検証を実施

お手本とパターンは類似

得られた結果の詳細な検証、及び高解像度での計算が課題

参考文献 (1)

はじめに・目標

GFDL, 2005: FMS: The flexible modeling system. http://www.gfdl.noaa.gov/~fms/

SWAMP Project, 1998: AGCM5. http://www.gfd-dennou.org/arch/agcm5/. 地球流体電脳倶楽部

モデルの解説(1)

Arakawa, A., Suarez, M. J., 1983: Vertical differencing of the primitive equations in sigma coordinates. Mon. Wea. Rev., 111, 34--35.

モデルの解説(2)

Balaji, V., 2002: The FMS Manual: A developer's guide to the GFDL Flexible Modeling System. http://www.gfdl.noaa.gov/~vb/FMSManual/FMSManual.html

Ruby Documentation Project, 2005: http://www.rubyist.net/rubikitch/RDP.cgi?cmd=view;name=top

竹広真一, 石岡圭一, 森川靖大, 小高正嗣, 石渡正樹, 林祥介, SPMODEL 開発グループ, 2004: 階層的地球流体力学スペクトルモデル集 (SPMODEL), http://www.gfd-dennou.org/arch/spmodel/, 地球流体電脳倶楽部.

豊田 英司, 塚原 大輔, 森川 靖大, 小高正嗣, 石渡 正樹, 堀之内 武, 林 祥介, 地球流体電脳倶楽部 Davis プロジェクト, 2004: gtool4 netCDF 規約, http://www.gfd-dennou.org/arch/gtool4/, 地球流体電脳倶楽部.

豊田 英司, 石渡 正樹, 竹広 真一, 林 祥介, gtool4 開発グループ, 2002: gtool4 Fortran90 ツール・ライブラリ, http://www.gfd-dennou.org/arch/gtool4/, 地球流体電脳倶楽部.

参考文献 (2)

動作試験

Held, I. M., and Suarez, M. J., 1994: A proposal for the intercomparison of the dynamical cores of atmospheric general circuation models. Bull. Am. Meteor. Soc., 75, 1825—1830.

付録

Arakawa, A., Suarez, M. J., 1983: Vertical differencing of the primitive equations in sigma coordinates. Mon. Wea. Rev., 111, 34--35.

Balaji, V., 2004: FMS: the GFDL Flexible Modeling System. http://www.gfd-dennou.org/arch/dcmodel/workshop/2004-09-23/FMS/pub/

Haltiner, G. J., and Williams, R. T., 1980: Numerical Prediction and Dynamic Meteorology. John Wiley & Sons, 477pp.

石岡圭一, 2004: スペクトル法による数値計算入門. 東京大学出版会, 232pp.

室井ちあし, 2001: 統一気象予測ソフトウェアの開発 その理想と実体. 地球流体セミナー 2001 年夏, http://www.gfd-dennou.org/seminars/gfdsemi/2001-08-22/muroi/pub-web. 地球流体電脳倶楽部

Satoh, M., 2004: Atmospheric Circulation Dynamics and General Circulation Models. Springer, 643pp.

時岡達志, 山岬正紀, 佐藤信夫, 1993: 気象の数値シミュレーション. 東京大学出版会, 247pp.

付録：FMS について

FMS (Flexible Modeling System)

目的

結果の考察を容易に (ただし地球限定)

仕様

FMS の利用者は「User code」部分のみを編集

Infrastructure で I/O や並列化を行う

Superstructure は land, ocean, sea ice の結合インターフェース提供

Balaji (2004) より

付録：球座標の格子系の例

左：正20面体格子、右：等角立方格子

Satoh (2004), Fig 20.1 より

付録：差分法とスペクトル法

時岡 他 (1993), 図1-1 より

付録：スペクトル法
(球面調和関数展開)

長所

差分近似に伴う誤差が生じない

Aliasing が無い (後ページ参照)

極の取り扱いが容易

短所

境界条件が複雑な場合に適用が難しい

スペクトル変換のコスト(時間)がかかる (特に分散型並列計算機には向かない)

詳しくは 石岡 (2004) や Satoh(2004) 参照

付録：Aliasing

Haltiner and Williams (1980) Fig5.9 より

d

エイリアジングの例。本来は実線で表される波長 4d/3 の波だが、差分系では格子点 (図中黒丸) の値しか見ないため、破線で描いたような波長4dの波に見えてしまう.

付録：Arakawa and Suarez (1983)

Arakawa and Suarez (1983)

座標のプリミティブ方程式系において、連続系のもついくつかの性質を保持する鉛直差分スキーム

全領域積分した質量を保存

全領域積分したエネルギーを保存

全領域積分した角運動量を保存

全領域積分した温位を保存

静水圧の式が local に決まる (下層の高度が上層のに依存しない

付録：時間積分スキーム

時間切れ・・・

GCMというのは全球大気の数値モデルで、おなじみの天気予報や研究では気候変動予測で使われてます。

では、僕が作りたいのはどういう GCM かというと、従来のものに比べて読みやすくて、書き換えが簡単にできるようなモデルなんです。例えば、GCM の結果を見ようとした時、普通の GCM は 3 次元の力学であるとか、物理過程と呼ばれる小さいスケールでの乱流や放射など、いろいろな要素が入ってるんですね。その際、結果を解釈しようとする時に物理過程を外したり、力学の部分を 2 次元や 1次元といった簡単なものにして比較できると簡単ですよね。今のモデルは高速で予報精度もすごい高いんですが、逆に読みやすさや、組み換えなどは犠牲にされてます。こういった結果を理解するために、可変性の高いモデルを目指したものでは、他にも GFDL の FMS ってのもあります。しかし、これはコードの読みやすさなどの点にはあまりこだわっておらず、僕がやろうとすることとはちょっと方向性が違います。

実際には、 Fortran90 の機能を最大限利用してこれらのことを可能にしようとしてます。サンプルとする GCM として AGCM5 を参照してます。

電脳倶楽部をどう説明しよう？

地球流体力学に関する計算機資源などなど、いろいろお役立ちなものの開発と公開を行っている団体で、日本全国の有志の地球流体関係の研究者が参加している。

1) FORTRAN77 を使ってきたこと、

実際には、 Fortran90 の機能を最大限利用してこれらのことを可能にしようとしてます。サンプルとする GCM として AGCM5 を参照してます。

今のモデルの方程式系は AGCM5 をほぼ踏襲した方程式系になっています。これは、もちろん AGCM5 を参照しているから、という理由もありますし、後でモデルの比較をする際に AGCM5 を利用できるため、という理由もあります。

今のモデルの方程式系は AGCM5 をほぼ踏襲した方程式系になっています。これは、もちろん AGCM5 を参照しているから、という理由もありますし、後でモデルの比較をする際に AGCM5 を利用できるため、という理由もあります。

渦度・発散って何？

比湿って何？ q = ρv /ρ (ρv ：水蒸気密度、ρ = ρd + ρv ：大気の密度)

各項の意味は？

今のモデルの方程式系は AGCM5 をほぼ踏襲した方程式系になっています。これは、もちろん AGCM5 を参照しているから、という理由もありますし、後でモデルの比較をする際に AGCM5 を利用できるため、という理由もあります。

惑星大気大循環モデルの開発： データI/Oライブラリと力学コアの設計

森川 靖大

2004-07-17T18:01:54