AGCM5 第3部コード解説

地球流体電脳倶楽部

2002 年 06 月 20 日 (AGCM5-3.4)

\underline{ii AGCM5 第 3 部コード解説}

目次

1  コードの概要  1

1.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 変数の命名法について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 メインプログラムの流れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2  力学過程  6

2.1 力学過程の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1.1 セミインプリシット時間積分の概要 . . . . . . . . . . . . . . . . . . . . . . 6 2.2 力学過程メイン DYNMCS の流れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 セミインプリシット時間積分の詳細 . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3.1 非重力波項の計算 (サブルーチン GRDDYN) . . . . . . . . . . . . . . . . . . 11 2.3.2 セミインプリシット時間積分 (サブルーチン TINTGR) . . . . . . . . . . . . 14

3  物理過程  15

3.1 物理過程の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 物理過程メイン PHYSCS の流れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.3 調節的物理過程 PADJST の流れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

A  プログラムの構造  19

A.1 サブルーチン呼び出しの概略 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/dai3bu.tex

\underline{AGCM5 第 3 部コード解説 iii}

A.2 変数の引渡し図 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 B  謝辞  30
/model/agcm/agcm5/doc/dai3bu/dai3bu.tex  2002/006/20(地球流体電脳倶楽部)

第1章コードの概要

1.1 はじめに

本稿は AGCM5 のプログラムの構造, 計算の手順を解説する.

1.1.1 変数の命名法について

モデル中の主要な変数の命名法について述べる. これはコードと本稿の理解に必要である.

[例] 変数 GTVOR は渦度 (VOR) の時間変化 (2 文字目が T) のグリッド値 (1 文字目が G) である. 同様に, GAVOR は渦度そのもの (2 文字目が A) のグリッド値である.

[命名規則] (予報変数と主要な時間変化項)

• 第 1 文字目: G の場合グリッド値, W の場合スペクトル (球関数展開係数).

• 第 2 文字目: T の場合時間変化率, A か B か D の場合変化率でない値そのもの. A, B はリープフロッグの早い方 (B) と遅い方 (A) の時刻の値を区別するために用いる. 区別の必要がないときは D が用いられることが多い.

• 3 文字目以降: 物理量を表す.

– U: 東西風速 u (グリッドの場合) または U ≡ u cos φ (スペクトルの場合)

– V: 南北風速 v (グリッドの場合) または V ≡ v cos φ (スペクトルの場合)

– VOR: 渦度 ζ

– DIV: 発散 D

– T: 気温 T

– TV: 仮温度 T_{v}

– PS: 地表面気圧 p_{s}

– PI: ln p_{s}

– Q: 比湿 q

97/03/06 堀之内武

\underline{2 AGCM5 第 3 部コード解説第 1 章コードの概要}

     但し, 時間変化率 GTUA, GTVA, GTH, GTR は物理過程 PHYSCS と力学過程 DYNMCS 中で徐々

に計算されるので, 3 文字目以降が物理量 UA, V A, H, R (後述. 第 1 部, 第 2 部も参照のこと) と一致するのは時間積分ループの終わり頃である.

• 半整数レベルに補間した変数の名前の最後には M をつける.

その他の変数について: 第 2 文字目の D 等が省略されたり (例: GSIGD: 鉛直速度 σ̇), さらに第 1 文字目も省略されたりするが, 比較的わかりやすい名前がつけられている. なお, 鉛直速度はもともと半整数レベルで定義されるため, 接尾辞 M を伴わない.

2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/gaiyou/intro.tex

\underline{AGCM5 第 3 部コード解説  1.2 メインプログラムの流れ  3}

1.2 メインプログラムの流れ

図 1.1 にメインプログラムとその直接の下請ルーチンから見たモデルの構造の概要を示す. さらに入出力関連の部分も書き加えた. AGCM5 が必要とする入力ファイルは初期値ファイルと地表面に関する境界値, 地形ファイルである. 但し, 初期値ファイル以外は省略して, デフォルトまたは陽に指定した一定値を用いることが出来る. 出力は結果を残すためのヒストリー出力と時間積分を延長するための再スタートファイルがある. 後者は延長積分時に初期値ファイルとして使われる.

モデルは初期化部分と時間積分のメインループの二つの部分に分けることが出来る. 前者は初期値の読み込みと若干の定数の設定を行う. 後者が言うまでもなくモデルの主要部分である. 境界値, 地形ファイルの読み込みはループに組み込まれてはいるが, 現在のところ下請のサブルーチン RDCON の仕様により, 最初の一回しか読み込まれないので時間変化しない. 改良が望まれる.

時間積分にはリープフロッグのセミインプリシット法が用いられているので, 時間の流れは単純ではない. 図 1.2 は時間積分ループの時間とデータの流れの概要を示す. 図に「時刻を更新するとき」とあるが, 通常のステップでは時刻 t − Δt の値を t での値に書き換えて次のステップにのぞむ. 更新されないのは, 単一時刻の初期値から出発するため時間刻を徐々に伸ばしていくときだけである. リープフロッグの計算モードの除去は時々の再スタートでなく, 時間フィルターにより行われる.

モデルの時間積分はリープフロッグスキームのセミインプリシット法により行われる. 物理過程で計算された時間変化項は力学過程に引き渡され, 時間積分は力学過程の中で行われる. ただし, 物理過程のうち, 降水過程などの調節過程は時間積分後に適用される.

97/03/06 堀之内武 /model/agcm/agcm5/doc/dai3bu/gaiyou/main.tex  2002/006/20(地球流体電脳倶楽部)

図 1.1: AGCM5 の構造の概要. アルファベット太字で書かれているのはメインプログラム GCM5 (aadmn.F) 中で呼ばれているサブルーチンである. さらに, 入出力に関連する部分を加えた. HISTIN により出力のためのデータを各所で HISTOU 下の作業配列にためておき, HISTOU のコールで一気に出力する. なお, メインで呼ばれているサブルーチンのうち, 骨格を把握するためは必要ないものは割愛した.

2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/gaiyou/main.tex

図 1.2: 時間積分ループの時間とデータの流れ. GA, GB はそれぞれ現在とその前のステップのグリッド値を表す. 時間変化率 GT のうち表面にでるのは物理過程によるもののみである. そこでは t − Δt での時間変化率を基本とするが, 拡散などはインプリシットで扱い, t + Δt での値を推定値する. ∘→ はデータの流れを, •→ は制御パラメターを表す.

/model/agcm/agcm5/doc/dai3bu/gaiyou/main.tex  2002/006/20(地球流体電脳倶楽部)

第2章力学過程

2.1 力学過程の概要

本章では力学過程のコードの構造と計算手順を説明する. ここでいう力学過程とはスペクトルで計算される部分のことを指し, 移流項や水平拡散項, 圧力勾配項などである. あるいは物理過程以外と言っても良いかもしれない. スペクトルとは球面調和関数展開の係数のことであり, スペクトルへの変換とは球関数展開に他ならない.

実際の時間積分は「力学過程」の中で行われる. 物理過程で計算された時間変化項は力学過程に引き渡され, 「力学過程」の時間変化項に加えられ, セミインプリシット法による時間積分に組み込まれる. 但し, 降水過程などの調節過程は時間積分後に適用される.

2.1.1 セミインプリシット時間積分の概要

セミインプリシット法についてはドキュメント第 2 部の 2.4 節に解説があるので参照のこと. 時刻 t − Δt から t + Δt へのセミインプリシットの計算手順は以下のようにまとめられる.

1. 発散項 D のみに関する単一のセミインプリシット方程式,

\underline{M} \overline{D}^{t} = f  (2.1)

を解く. ここで\overline{D}^{t} ≡ \frac{1}{2}(D^{t+Δt} + D^{t−Δt}) である.

2. \overline{D}^{t} と時刻 t における時間変化の非重力波 (NG) 項から各変数の t + Δt での値を求める.

以上の手順の説明は 2.3 節で行う.

97/03/06 堀之内武

2.2 力学過程メイン DYNMCS の流れ

図 2.1 に力学過程メインルーチン DYNMCS の流れを示す. 太文字はサブルーチン名, タイプライター文字は変数名, 変数名等のうしろの括弧は中身や時刻に関する補足説明である. 変数の命名法は第 1 章に述べた. 内部変数の一部は変数名でなく物理量名をそのまま記してある. このうち\underline{M} はセミインプリシット用行列 (2.1 参照), その他の量についてはは第 1, 2 部を参照のこと.

DYNMCS で呼ばれている GRDDYN は時間変化率の非重力波 (NG) 項を計算する. 物理過程 PHYSCS から引き渡された時間変化項 GT* (* には適当な物理量名が入る) に, 時刻 t のデータから計算した非線型項等を加えたものが NG 項である. TINTGR はセミインプリシット時間積分 (2.1.1 の 1,2 とも) を行う. 以上の 2 つが DYNMCS の主な下請ルーチンであるので, 節を改めて詳しく説明する.

GRDDYN 以前に呼ばれるサブルーチンは必要な定数等の準備にあてられる. DMTRIX までは基本的に時間に依存しない. SETPHS は外部のファイルから地表面のジオポテンシャルを読み込む. 但し, 現在のところ下請のサブルーチン RDCON の仕様により, 最初の一回しか読み込まれないようになっている. なお, 地表面高度を一定とする場合はファイルを用意しなくて良い. TINTGR の前後, 計 4 つのサブルーチンは主にグリッド値とスペクトルの間の変換を行う. その後, 水平粘性に失われたなわれた運動エネルギーを補償する摩擦熱を行い, DYNMCS 外に出ない量の出力 (HISTIN) を行う. 但し, ここでは配列にためておくだけで, 実際の出力はメインループの最後の方でまとめて行われる (図 1.1). 最後の「診断」はモニターのためであり, 時間積分には影響しない.

97/03/06 堀之内武 /model/agcm/agcm5/doc/dai3bu/dynamics/dynmcs.tex  2002/006/20(地球流体電脳倶楽部)

\underline{8 AGCM5 第 3 部コード解説第 2 章力学過程}

図 2.1: 力学過程メイン DYNMCS の流れ. 太文字はサブルーチン名, タイプライター文字は変数名, 変数名等のうしろの括弧は中身や時刻に関する補足説明である. DYNMCS が呼ぶサブルーチンと引数を示しているが, 座標変数などの定数の入力はほとんど省略した. DYNMCS の内部変数の一部は変数名でなく物理量名をそのまま記した. 呼ばれているサブルーチンの中で実質的な計算をしてい

るのは GRDDYN と TINTGR である.

2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/dynamics/dynmcs.tex

\underline{AGCM5 第 3 部コード解説  2.3 セミインプリシット時間積分の詳細  9}

2.3 セミインプリシット時間積分の詳細

2.1.1 節で述べたように, セミインプリシット時間積分は 2 つのステップに分けられる. これを詳しく書くと以下のようになる.

●第 1 段階: \overline{D}^{t} の計算

\underline{M} \overline{D}^{t} = f  (2.2)

という形の式（ドキュメント第 2 部の (2.55) 式: 間違い有り）を解くのが第 1 段階である. ここでは第 2 部の 2.4 節と同様, 太字で鉛直に離散化したベクトルを表し, 下線で行列を表す. 水平方向に関しては, 各スペクトル係数が独立に計算されるので離散化したことが表向きは現れない. D は発散で,

\overline{D}^{t} ≡ \frac{1}{2}(D^{t+Δt} + D^{t−Δt})  (2.3)

である. 行列\underline{M} は

\underline{M} ≡ (1 − 2ΔtD_{H})(1 − 2ΔtD_{M})\underline{I} − (Δt)^{2}[\underline{W} \underline{h} + (1 − 2ΔtD_{H})GC^{T}]∇_{σ}^{2}  (2.4)

であり (第 2 部の (2.55) 式と若干違うことに注意), ベクトル f は,

f ≡ (1 − 2ΔtD_{H})(1 − ΔtD_{M})D^{t−Δt} + (1 − 2ΔtD_{H})Δt (\frac{∂D}{∂t} )_{NG}

−Δt∇_{σ}^{2} \{(1 − 2ΔtD_{H})Φ_{S} + \underline{W} [(1 − ΔtD_{H})T^{t−Δt} + Δt (\frac{∂T}{∂t} )_{NG}]

+(1 − 2ΔtD_{H})G [π^{t−Δt} + Δt (\frac{∂π}{∂t} )_{NG}]\} .  (2.5)

である. ここでは第 2 部の (2.55) 式では (1 − 2ΔtD_{H}) となっている T^{t−Δt} の係数を (1 − ΔtD_{H}) に直していることに注意.

\underline{M} はサブルーチン DMTRIX で計算され, LU 分解される. \underline{M} は時間刻みの大きさが変らない限り設定し直す必要がない. f の計算と (2.2) を解く作業はサブルーチン TINTGR で行われる.

●第 2 段階: 時間積分時間積分も TINTGR が行う. 第 1 段階で\overline{D}^{t} が求まったので発散の時間積分は容易である:

D^{t+Δt} = 2\overline{D}^{t} − D^{t−Δt}.  (2.6)

他の物理量についても, \overline{D}^{t} と時刻 t における時間変化率の非重力波 (NG) 項から, t + Δt での値が求まる:

X^{t+Δt} = γ_{X}^{−1} \{X^{t−Δt} + 2Δt [(\frac{∂X}{∂t} )_{NG} + G_{X} \overline{D}^{t}]\} ,  (2.7)

97/03/06 堀之内武

97/05/19 堀之内武 /model/agcm/agcm5/doc/dai3bu/dynamics/semiimp.tex  2002/006/20(地球流体電脳倶楽部)

\underline{10 AGCM5 第 3 部コード解説第 2 章力学過程} ここで,

\begin{aligned}X & = & π, & ζ, & T , & q, \\ (\frac{∂X}{∂t} )_{NG} & = & (\frac{∂π}{∂t} )_{NG} , & (\frac{∂ζ}{∂t} )_{NG} , & (\frac{∂T}{∂t} )_{NG} , & (\frac{∂q}{∂t} )_{NG} , \\ γ_{X} & = & 1, & (1 − 2ΔtD_{M}), & (1 − 2ΔtD_{H}), & (1 − 2ΔtD_{H}), \\ G_{X} & = & −C^{T} , & \underline{0}, & −\underline{h}, & \underline{0}. \end{aligned}  (2.8)

以下の各小節で, 具体的な計算手順とプログラムソースとの対応を述べる. ソースを見る際にはタグファイルを用いた支援システム等を利用せよ (cf. 「極楽 AGCM 体験」 ).

2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/dynamics/semiimp.tex

当該サブルーチンからの出力変数 (入力に上書きする場合を含む) の出力時 GTU•
当該サブルーチン内の作業変数または, 入/出力変数が入/出力時の値と違う場合GTU
当該サブルーチンへの入力変数  GTU
当該サブルーチンへの入力に上書きして出力される変数の入力時  GTU•

表 2.1: 数式と対応させて書く場合のモデル変数の表記法. 変数名 GTU の場合を例にして書く. 表の上二行 (黒) は当該サブルーチン内で設定/変更された変数を表し, 下二行 (灰) は引数として入力された変数を表す. 入力された変数もサブルーチン内で値が書き変わると色が黒くなる. 本文中で入力変数の下に付した文字列は, 値を設定したサブルーチン名を示す. 作業変数等の下の文字列は, 値を設定する下請サブルーチンを示す.

2.3.1 非重力波項の計算 (サブルーチン GRDDYN)

非重力波 (NG) 項はサブルーチン GRDDYN で計算される.

NG 項はドキュメント第 2 部の式 (2.38)∼(2.41), (2.9), (2.22) で与えられる. これらを再掲し, 各項の下にモデル中の変数名を記す. 以下で, 数式と対応させて書く場合のモデル変数の表記法は表 2.1 に従って統一する.

(\frac{∂π}{∂t} )_{NG} = \underbrace{Z_{k}}_{GTPI •}  (2.9)

(\frac{∂D}{∂t} )^{NG} = \frac{∂}{a(1 − μ^{2})∂λ} \underbrace{(UA)_{k}}_{GTUA •} +\frac{∂}{a∂μ} \underbrace{(VA)_{k}}_{GTVA •} −∇_{σ}^{2} \underbrace{[(KE)_{k} + ∑_{k=1}^{K} W_{kl}(T_{v,l} − T_{l})]}_{GTKE •} (2.10)

(\frac{∂ζ}{∂t} )^{NG} = \frac{∂}{a(1 − μ^{2})∂λ} \underbrace{(VA)_{k}}_{GTVA •} −\frac{∂}{a∂μ} \underbrace{(UA)_{k}}_{GTUA •}  (2.11)

(\frac{∂T_{k}}{∂t} )^{NG} ≡ −\frac{∂}{a(1 − μ^{2})∂λ} \underbrace{(U_{k}T′_{k})}_{GTUT •} −\frac{∂}{a∂μ} \underbrace{(V_{k}T′_{k})}_{GTVT •} + \underbrace{H_{k}}_{GTH •}  (2.12)

(\frac{∂q_{k}}{∂t} )^{NG} ≡ −\frac{∂}{a(1 − μ^{2})∂λ} \underbrace{(U_{k}q′_{k})}_{GTUQ •} −\frac{∂}{a∂μ} \underbrace{(V_{k}q′_{k})}_{GTVQ •} + \underbrace{R_{k}}_{GTR•}  (2.13)

97/03/06 堀之内武 /model/agcm/agcm5/doc/dai3bu/dynamics/grddyn.tex  2002/006/20(地球流体電脳倶楽部)

GRDDYN の出力変数は, 以上の 10 個の時間変化項と, 2 つの出力変数

GSIGD •_{PSDOT} ≡ σ̇,
GPID • ≡ \frac{∂π}{∂t} + ∑_{l=1}^{K} v_{l} · ∇πΔσ_{l} = − ∑_{l=1}^{K} D_{l}Δσ_{l},

である^{1}.

GRDDYN の入力変数は時刻 t でのグリッド値 GA * (∗ = 変数名) と, 物理過程 PHYSCS で計算された u, v, T, q の時間変化項 (それぞれ GTUA •, GTVA •, GTH •, GTR• に格納; これらは PHYSCS の中では GTU, GTV, GTT, GTQ という名前である) と, 水平・鉛直離散化や熱力学に関する定数配列である.

次に時間変化項の具体的な中身とモデル変数の対応を見ることにする. 但し 1.1.1 節の命名規則から明らかなものは省略する. 第 2 部の式 (2.10), (2.12), (2.42), (2.23) より,

\underbrace{Z_{k}}_{GTPI •} = − \underbrace{∑_{k=1}^{K} v_{k} · \overbrace{∇π}^{(2.15),(2.16) 参照} Δσ_{k}}_{GSUMPA_{k=1}_{PSDOT}}  (2.14)

\underbrace{(UA)_{k}}_{GTUA •} = (ζ_{k} + \underbrace{f}_{CORIOL_{DSETC}} )V_{k} − \frac{1}{2Δσ_{k}} [σ̇_{k−1/2}(U_{k−1} − U_{k}) + σ̇_{k+1/2}(U_{k} − U_{k+1})]

−C_{p} \underbrace{κ̂_{k}}_{SKAPPA_{DSETV}} \overbrace{T′_{v,k}}^{GATVED}\frac{1}{a}\underbrace{\frac{∂π}{∂λ}}_{GAPIX_{PIXPIY}}+\underbrace{F_{λk}}_{GTUA •_{PHYSCS}} cos φ ,  (2.15)

\underbrace{(VA)_{k}}_{GTVA •} = −(ζ_{k} + \underbrace{f}_{CORIOL_{DSETC}} )U_{k} − \frac{1}{2Δσ_{k}} [σ̇_{k−1/2}(V_{k−1} − V_{k}) + σ̇_{k+1/2}(V_{k} − V_{k+1})]

−C_{p} \underbrace{κ̂_{k}}_{SKAPPA_{DSETV}} \overbrace{T′_{v,k}}^{GATVED}\frac{1}{a} \underbrace{(1 − μ^{2})\frac{∂π}{∂μ}}_{GAPIY_{PIXPIY}} +\underbrace{F_{φk}}_{GTVA •_{PHYSCS}} cos φ ,

\underbrace{(KE)_{k} + ∑_{k=1}^{K} W_{kl}(T_{v,l} − T_{l})}_{GTKE •} = \frac{U_{k}^{2} + V_{k}^{2}}{2(1 − μ^{2})} + \underbrace{∑_{k=1}^{K} W_{kl}(T_{v,l} − T_{l})}_{GHYDRO を呼ぶと W_{kl} がかかる} ,  (2.16)

1GSIGD は当然だろうが GPID を出力する意味は良くわからない. \frac{∂ π}{∂t} ではないことに注意.

2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/dynamics/grddyn.tex

\underbrace{H_{k}}_{GTH •} =\underbrace{T′_{k}}_{GDTED}D_{k} + \underbrace{\frac{Q_{k}}{c_{p}}}_{GTH •_{PHYSCS}}

− \frac{1}{Δσ_{k}} [σ̇_{k−1/2}( \underbrace{T̂′_{k−1/2}}_{GATEDM} −\underbrace{T̂′_{k}}_{GATED} ) + σ̇_{k+1/2}(T′_{k} − T̂′_{k+1/2})]

− \frac{1}{Δσ_{k}} [ \underbrace{σ̇_{k−1/2}^{NG}}_{GSIGDA_{PSDOT}} (\underbrace{ˆT̄_{k−1/2}}_{TBARM_{DSETV}}−\underbrace{T̄_{k}}_{TBAR_{DSETV}} ) + σ̇_{k+1/2}^{NG}(T̄_{k} − ˆ T̄_{k+1/2})]

+\underbrace{κ̂_{k}}_{SKAPPA_{DSETV}} T_{v,k} \underbrace{v_{k} · ∇π}_{GAPADV_{PSDOT}}

−\frac{1}{Δσ_{k}} \underbrace{α_{k}}_{SALPHA_{DSETV}} T_{v,k} \underbrace{∑_{l=k}^{K} v_{l} · ∇πΔσ_{l}}_{GSUMPA_{PSDOT}} −\frac{1}{Δσ_{k}} \underbrace{β_{k}}_{SBETA_{DSETV}} T_{v,k} \underbrace{∑_{l=k+1}^{K} v_{l} · ∇πΔσ_{l}}_{GSUMPA_{PSDOT}}

−\frac{1}{Δσ_{k}} α_{k} \underbrace{T′_{v,k}}_{GATVED} \underbrace{∑_{l=k}^{K} D_{l}Δσ_{l}}_{GSUMPD_{PSDOT}} −\frac{1}{Δσ_{k}} β_{k} \underbrace{T′_{v,k}}_{GATVED} \underbrace{∑_{l=k+1}^{K} D_{l}Δσ_{l}}_{GSUMPD_{PSDOT}} .  (2.17)

ここで,

\underbrace{σ̇_{k−1/2}^{NG}}_{GSIGDA_{PSDOT}} = −σ_{k−1/2} \underbrace{(\frac{∂π}{∂t} )^{NG}}_{(−1)×GSUMPA_{k=1}_{PSDOT}} − \underbrace{∑_{l=k}^{K} v_{l} · ∇πΔσ_{l}}_{GSUMPA_{PSDOT}} .  (2.18)

そして,

\underbrace{R_{k}}_{GTR•} = q_{k}D_{k} + \underbrace{R_{PHYSCS,k}}_{GTR•_{PHYSCS}}

−\frac{1}{2Δσ_{k}} [σ̇_{k−1/2}(q_{k−1} − q_{k}) + σ̇_{k+1/2}(q_{k} − q_{k+1})] .  (2.19)

/model/agcm/agcm5/doc/dai3bu/dynamics/grddyn.tex  2002/006/20(地球流体電脳倶楽部)

2.3.2 セミインプリシット時間積分 (サブルーチン TINTGR)

9 ページで述べた時間積分の 2 つの手続きは TINTGR で行われる. ここでは TINTGR 内の計算の流れを, 物理量とモデル変数の対応を示しながら説明する.

＊ (2.2) の f を求める:

時間積分の第 1 段階は (2.2) を解くことである. 時間に依存しない行列\underline{M} は DMTRIX で計算され, LU 分解されて DMLU として TINTGR に引き渡されるので, 最初にするべきことは f を求めることである. (2.5) を再掲し, 各項の下にモデル内の変数名を記す.

\underbrace{f}_{WBRDIV} ≡ (1 − 2Δt\underbrace{D_{H}}_{DIFT_{DSETC}})(1 − Δt\underbrace{D_{M}}_{DIFV_{DSETC}} )\underbrace{D^{t−Δt}}_{WDDIV •_{GD2WD}} +(1 − 2ΔtD_{H})Δt \underbrace{(\frac{∂D}{∂t} )_{NG}}_{WTDIV_{TENG2W}}

−Δt \underbrace{∇_{σ}^{2}}_{FLAPLA_{DSETC}}\{(1 − 2ΔtD_{H}) \underbrace{Φ_{S}}_{WPHIS_{SETPHS}}+ \underbrace{\underline{W} \overbrace{[(1 − ΔtD_{H})\underbrace{T^{t−Δt}}_{WDT •_{GD2WD}}+Δt \underbrace{(\frac{∂T}{∂t} )_{NG}}_{WTTTENG2W}]}^{WDTX}}_{WDPHIX_{WHYDRO を呼ぶと\underline{W} がかかる}}

+(1 − 2ΔtD_{H})\underbrace{G}_{GG_{DSETV}}\underbrace{[ \overbrace{π^{t−Δt}}^{WDPIGD2WD} +Δt \overbrace{(\frac{∂π}{∂t} )_{NG}}^{WTPITENG2W} ]}_{WDPTX}\} .  (2.20)

上式の WD * • はサブルーチン GD2WG において時刻 t − Δt のグリッド値 GB∗ からスペクトルに変換され, WT * は TENG2W で GT∗ からスペクトルに変換される. 後者は変換時に (2.9)∼(2.13) の微分操作や加算を行う.

＊ (2.2) を解く:

サブルーチン LUSOLV により解かれる. 結果は WBRDIV に上書きされる. その他の LUSOLV の引数は\underline{M} とそのピボット (DMLU と KPIV ) と配列の大きさに関する定数である.

＊ (2.6) と (2.7) を解く:

モデル中の変数との対応は, C (≡ \{Δσ_{k}\}) が DSIG , \underline{h} が TMTX である. 非重力波 (NG) 項は WT *, スペクトル値は WD * • → WD * • となる.

97/03/06 堀之内武 2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/dynamics/tintgr.tex

第3章物理過程

3.1 物理過程の概要

現在の AGCM5 に存在する物理過程は以下の通り.

• 放射過程 (RADFLX)

標準のバンド数は, 長波が 4, 短波が 1 (変更可).

• 鉛直拡散 (VDFUVT)

標準は Yamada and Mellor のレベル 2 クロージャー. オプションで係数を一定にできる.

• 地表フラックス (SFCUVT)

標準はバルク法. オプションで係数を一定にできる.

• 重力波抵抗 (GRAVTY)

McFarlane のパラメタリゼーション

• 積雲パラメタリゼーション (CUMADJ) 標準は Kuo スキーム. オプションで湿潤対流調節

• 大規模凝結 (LSCOND)

• 乾燥対流調節 (DADJST)

• 負の水蒸気除去 (RMNEGQ)

以上のうち, 重力波抵抗まではサブルーチン PHYSCS にまとめられており, 積雲パラメタリゼーション以降はサブルーチン PADJST にまとめられている. 前者は次のステップを予報するための時間変化率を計算し, 後者は予報の結果を「調節」する.

NAMELIST 入力や非標準ルーチンのリンクにより, 各過程の削除や様々な設定変更ができる. また, 利用者が定義した過程を導入することもできる. 詳しくは「利用者の手引き」とドキュメント第 4 部を参照のこと.

97/03/06 堀之内武

\underline{16 AGCM5 第 3 部コード解説第 3 章物理過程}

3.2 物理過程メイン PHYSCS の流れ

降水などの「調節」過程以外の物理過程は, PHYSCS にまとめられている. その内容は,

   • 放射過程 RADFLX
   • 鉛直拡散 VDFUVT

• 地表フラックス SFCUVT

   • 重力波抵抗 GRAVTY

である.

PHYSCS では時刻 t − Δt でのグリッド値から各予報変数の時間変化率を求める. 従って PHYSCS で計算される時間変化率は基本的に時刻 t − Δt でのそれであるが, 鉛直拡散と地表面過程 (バルク法など) については, 陰解法的な扱いにより時刻 t − Δt における時間変化率を推定する. これは, 陰解法的な扱いをする過程による時間変化だけを取り込む擬似的なものである.

PHYSCS は時間変化率を計算するだけであり, 実際の時間積分は「力学過程」 DYNMCS で行われる. 第 1 章を参照のこと.

図 3.1 に力学過程メインルーチン PHYSCS の流れを示す. 太文字はサブルーチン名, タイプライター文字は変数名を表す. 変数名中の GT[UV] などは “GTU と GTV” を表す. GT[U|V] も同様である. [|U] は “” または “U” である.

第 1 章で述べたように, GD から始る変数は時刻 t − Δt における物理量のグリッド値を表す. GT から始る変数は時間変化率である.

GR から始る変数は地表面に関する量で, HC は熱容量, , ALB はアルベド, BET は湿潤度 (蒸発のしやすさ), Z0 は粗度長, ZSD は高度分散 (重力波抵抗用) である.

FLUX で終わるのは鉛直フラックスで, 半整数レベルで定義される. 鉛直フラックスの最初の文字 U, V, T, Q は, それぞれ u, v, T, q の時間変化率の式に現れるフラックスであることを示す. R, S はそれぞれ長波, 短波のフラックスであるから, これも T の時間変化に関係する. DRFLUX は長波の補正のために使われるものである.

MTX で終わるのは陰解法のための行列である. MTXS 終わるのはその中でも地表面部分だけに関するもので, PHYTND 中で MTX に組み込まれる.

3.3 調節的物理過程 PADJST の流れ

PADJST は最新の時刻のグリッド値を「調節」する. PADJST 内の各過程は以下の順序で適用される.

図 3.1: 物理過程メイン PHYSCS の流れ. 太文字はサブルーチン名, タイプライター文字は変数名を示す. GT[UV] や GT[U|V] は GTU と GTV を表す. 個々の変数名の見方については本文を参照のこと.

/model/agcm/agcm5/doc/dai3bu/physics/physcs.tex  2002/006/20(地球流体電脳倶楽部)

\underline{18 AGCM5 第 3 部コード解説第 3 章物理過程}

1. 負の水蒸気除去 RMNEGQ

2. 積雲パラメタリゼーション CUMADJ

3. 大規模凝結 LSCOND

4. 乾燥対流調節 DADJST

5. 負の水蒸気除去 RMNEGQ (2 回目)

第A章プログラムの構造

A.1 サブルーチン呼び出しの概略

★セットアップ

┏ YPREP : システム前処理 ┣ SETPAR* : 時間積分・時間制御パラメタ読み込み設定

┃┗┳━ ASUTIM ・ ASDELT  :  Atime routine への設定
┃┃  start ・ restart ファイルの指定

┃┗━ PCONST : 物理定数の設定 ┃ ┣ SETCOR : モデル座標の設定 ┃┗┳━ SETLON : 経度の計算 ┃┣━ SETLAT : 緯度の計算 ┃┃┗━━━ GAUSS : ガウス緯度の計算

┃┣━ SETSIG         : σ 1 の設定
┃┃┗━━━ SETSGM : σ 1/2 の指定

┃┃┗━━━ SETSIX* : σ 1/2 の読み込み

┃┗━ HSCORD         :  Icord routine への設定

\underline{20 AGCM5 第 3 部コード解説付録 A プログラムの構造}

┃ ┃ ┃ ┃★初期値読み込み生成 ┃ ┣ RDSTRT : 初期値の読み込みと再生成 ┃┗━┳━ IFLOPN : ファイルのオープン ┃┣━ RDGDST : 初期値の読み込み

90/11/30 山中康裕

94/02/08 保坂征宏 TeX 化

97/03/08 堀之内武 update

┃┃┗━━━ GDREAD :  GTOOL3 入力

┃┣━ RDGDSV : 渦度発散の読み込み ┃┃┗━━━ GDREAD : 力学過程のセットアップ ┃┗━ DSTRT : 初期値生成 ┃┗━┳━ DSETC : 力学過程のセットアップ ┃┃┗━┳━ DYNPAR*: 力学パラメタ（水平拡散 etc）の設定 ┃┃┣━ GAUSS : ガウス緯度荷重の計算 ┃┃┣━ SPSTUP : スペクトル変換の準備

┃┃┃┗━━━ RFFTIM ・ SETPNM

┃┃┣━ SETCOF : スペクトル計算の係数の設定 ┃┃┗━ SETDIF : 水平拡散の係数のセット

┃┗━ STRTV          :  n-1/n step の Grid → Wave → Grid 生成
┃┗━┳━ GUV2V  :  U,V(Grid) → D,phi(Wave)
┃┣━ G2W    :  T/Ps/Q(Grid) → T/Ps/Q(Wave)

┃┃┗━━━ SPG2W

┃┗━ GENGD  :  D,Phi,T,Ps,Q,(Wave) → U,V,T,Ps,Q(Grid)
┃┗━━━ W2G ／ SPG2W

┃ ┣ RDSTRG : 地表初期値の読み込み ┃┗━┳━ IFLOPN : ファイルのオープン

┃┗━━━ GDREAD     :  GTOOL3 入力

┃

┣ SETTIM                 :  A 時刻を記憶 (step routine TIMSTP への設定)
┣ AHSTRG                 :  U,V,T,Ps,Q,Vor,Div,Tg の Ihist routine への登録

┃┗━━━ HISTRG* : hist 管理への登録 (個別変更読み込み) ┃┗━━━ IFLOPN : ファイルのオープン ┃ ┃━━━━━━━━メインループ━━━━━━━━━━━━━━━━━━━━━━━━━━ ━━ ┃ ┃★時間制御 ┃

┣ TIMSTP                 : フラグセット・ Astep routine 内部時間更新準備

┃ ┃★ (ｔ−Δｔ) のセーブ ┃

┣ SAVEGB                 :  Asavb routine へ U,V,T,Ps,Q,Vor,Div をセーブ

\underline{AGCM5 第 3 部コード解説  A.1 サブルーチン呼び出しの概略  21}

┃┃ ┃ ┃★物理過程 ┃

┣ PHYSCS
┃┗━┳━ PHSTRG         :  Evap,Sens,Rain,OLR,OSR,
┃┃┃  Qvdf,Qcnd,QradL,QradS の Ihist routine への登録

┃┃┗━━━ HISTRG*: hist 管理への登録 (個別変更読み込み) ┃┃┗━━━ IFLOPN : ファイルのオープン ┃┣━ TINTP : σ 1/2 の温度補間を診断 ┃┣━ PZLEVS : 気圧・高度を診断 ┃┣━ PGRNDM : 地表の諸量（地表面のファイル／標準値読み込み） ┃┃┗━━━ RDCON : SST ／熱容量／α／β／Ｚ 0 ／高度分散／地表状態の読み込み・設定 ┃┃ : 地表モデルからの境界条件 (現在地表モデルと言えるようなものはない) ┃┃ ┃┣━ PHMTX0 : インプリシット行列の初期化 (時間変化項) ・フラックスの初期化 ┃┣━ RADFLX : 放射フラックスの計算 ┃┃┗━┳━ LNGFLX*: 長波フラックスの計算 ┃┃┣━ SHTINS : 短波入射・入射角 ┃┃┗━ SHTFLX : 短波フラックスの計算 ┃┣━ VDFUVT* : U,V,T,Q の鉛直拡散フラックスの計算

┃┃┣━━━ HISTRG*:  hist 管理への登録 (拡散係数 DF[MHE])

┃┃┃┗━━━ IFLOPN : ファイルのオープン

┃┃┣━━━ VDFCOF*: 拡散係数の計算 (Fix/K-theory/Meller-Yamada)
┃┃┣━━━ SHLCOF : 浅い積雲対流：拡散型 (Tiedtke による parametarization)
┃┃┗━━━ HISTIN :  Ihist routine へ出力 (拡散係数 DF[MHE])
┃┣━ SFCUVT*        :  U,V,T,Q の地表フラックスの計算
┃┃┗━━━ BLKCOF*:  U,T,Q のバルク係数 (Kondo/Louis)

┃┣━ PHYTND : U,V,T,Q の物理過程による温度変化率の計算 ┃┃┗━━━ LUMAK3 ・ LUSOL3 : 行列のＬＵ分解・ＬＵ分解による解の計算 ┃┣━ FLXCOR : U,V,T,Q の鉛直拡散フラックスの (新時刻での) 補正 ┃┣━ RADCOR : 地表面温度・地表気温の変化による放射フラックスの補正 ┃┣━ RADTND : 放射による温度変化率

┃┣━ GRAVTY         : 重力波抵抗 (McFarlane)
┃┣━ HISTIN         :  Evap,Sens,OLR,OSR,Qvdf,QradL,QradS 等の Ihist routine

へ出力 ┃┃ : 拡散＋放射による時間変化率の合計

┃┣━ PDIAG          : 物理過程の診断 (Suface,Top,Water)

┃┃┗━━━ QDIAGP*: 診断の制御（パラメタの読み込み） ┃┗ TGINTG : 地表面温度の積分 ┃ ┃★力学過程・時間積分 ┃

┣ DYNMCS

┃┗━┳━ DSETC : 力学過程のセットアップ（前述. 通常ここでは実行せず） ┃┣━ DSETV : 力学過程のセットアップ ┃┃┗━┳━ SETVDF : 鉛直差分ファクタの設定 ┃┃┃┗━━━ HYDROS ・ HYDROQ : 静水圧の係数設定 ┃┃┗━ SETMTX : セミインプリシット行列用の係数設定

┃┣━ DHSTRG         : σ,lnPs の時間変化率の Ihist routine への登録

┃┃┗━━━ HISTRG*: hist 管理への登録 (個別変更読み込み) ┃┃┗━━━ IFLOPN : ファイルのオープン ┃┣━ DMTRIX : セミインプリシット渦度方程式用の LU 行列 ┃┃┗━━━ LUMAKE : 行列のＬＵ分解 ┃┣━ SETPHIS : 地形 (Φ s) の読み込みとスペクトルへの変換 ┃┃┣━━━ RDCON : 定数ファイルの読み込み

┃┃┃┗━━━ IFLOPN/GDREAD

┃┃┗━━━ G2W : 格子→スペクトル ┃┃ ┃┣━ GRDDYN : 格子点上の力学項の計算 ┃┃┗━━━ PSDOT : 地表気圧・σ速度の計算 ┃┣━ TENG2W : 時間変化項の計算（スペクトルへ変換） ┃┃┗━━━ G2W : 格子→スペクトル

┃┣━ GD2WD          : 予報変数 (D,Phi,T,Ps,Q) のスペクトルへの変換

┃┃┗━━━ G2W : 格子→スペクトル ┃┣━ TINTGR : スペクトル上での時間積分 ┃┃┗━┳━ WHYDRO : ジオポテンシャルのスペクトルを計算 ┃┃┗━ LUSOLV : 発散の時間平均スペクトルを解く

┃┃  :  Vor,Div,T,Q,lnPs の時間変化項を求め積分

┃┣━ MOVGD : スペクトルの更新

┃┣━ GENGD          :  D,Phi,T,Ps,Q(Wave) → U,V,T,Ps,Q(Grid)
┃┃┗━━━ W2G ／ SPG2W

┃┣━ TENDIF : 運動量の拡散項の計算 ┃┣━ CORDIF : 運動量の拡散項による熱補正 ┃┃┗━━━ W2G : スペクトル→格子

┃┣━ HISTIN         : σ,lnPs の時間変化率の Ihist routine へ出力
┃┣━ DDIAG          :  (total energy,enstrophy,KE,IE,PE,LE,total mass) の出

力

┃┃┗━┳━ DDIHRG :  TE etc の Ihist routine への登録

┃┃┃┗━━━ HISTRG*: hist 管理への登録 (個別変更読み込み) ┃┃┃┗━━━ IFLOPN : ファイルのオープン

┃┃┣━ GHYDRO : 温度→ジオポテンシャル (potential energy のため)
┃┃┗━ HISTIN :  TE etc の Ihist routine への出力
┃┗━ SPCTLD         :  D,Phi,T,Ps,Q(Wave) の診断

\underline{24 AGCM5 第 3 部コード解説付録 A プログラムの構造} ┃┗━┳━ QDIAGP*: 診断の制御（パラメタの読み込み） ┃┗━ PRINTD : 診断の印字出力 ┃ ┃★調節 ┃

┣ PADJST                 : 物理過程２

┃┗━┳━ TINTP : σ 1/2 の温度補間を診断 ┃┣━ PZLEVS : 気圧・高度を診断 ┃┣━ RMNEGQ : 負の水蒸気量除去 ┃┣━ CUMADJ* : 積雲パラメタリゼーション (Kuo スキーム／ Manabe の湿潤対流調節) ┃┃┗ CCLLEV : 雲層のチェック (Kuo スキーム)

┃┣━ LSCOND*        : 大規模凝結

┃┣━ DADJST : 乾燥対流調節 ┃┣━ RMNEGQ : 負の水蒸気量除去 (2 回目)

┃┣━ HISTIN         :  Rain,Qcnd,Scnd 等の Ihist routine へ出力
┃┗━ PDIAG2         : 物理過程の診断 (Rain)

┃┗━━━ QDIAGP*: 診断の制御（パラメタの読み込み）

┣ TFILT
┃┗━━━ LOADGB         :  Asavb routine から U,V,T,Ps,Q,Vor,Div をロード

┃ ┃★ヒストリ出力 ┃

┣ HISTIN                 :  U,V,T,Ps,Q,Vor,Div,Tg の Ihist routine への出力
┣ HISTRP                 :  Ihist routine に登録された物理量名の一覧印字

┣ HISTOU : Ihist routine に登録された物理量の印字・ファイル出力 ┃┣━━━ HHEDIT ・ GZDBWX : ヘッダー編集・精度変換書き込み

┃┗━━━ HMONIT

┃┗━ PRINTD : 診断の印字出力 ┃ ┃★時刻変更 ┃

┣ ADVSTP                 :  Astep routine の時間更新

┃ ┃★リスタートファイルの書き出し ┃ ┣ WRRSTR : リスタートファイルの書き出し ┃┗━┳━ IFLOPN : ファイルのオープン

┃┣━ WRRSTX         : 基本リスタートファイル U/V/T/Ps/Q (steps n-1&n) の書き出

し

┃┃┗━ GDWRIT     :  GTOOL3 出力

┃┗━ WRRSTY : 渦度発散リスタートファイル (steps n-1&n) の書き出し ┣ WRRSTG : 地表リスタートファイルの書き出し ┃┗━┳━ IFLOPN : ファイルのオープン

┃┗━ GDWRIT         :  GTOOL3 出力

\underline{26 AGCM5 第 3 部コード解説付録 A プログラムの構造} ┗ YFINE : システム後処理

\underline{AGCM5 第 3 部コード解説  A.2 変数の引渡し図  27}

A.2 変数の引渡し図

この図は現在の AGCM5 には対応してないので改定が必要である.

92/??/?? 沼口敦 94/02/10 保坂征宏 TeX 化

/model/agcm/agcm5/doc/dai3bu/appndx/senzu.tex  2002/006/20(地球流体電脳倶楽部)

2002/006/20(地球流体電脳倶楽部)  /model/agcm/agcm5/doc/dai3bu/appndx/senzu.tex

\underline{AGCM5 第 3 部コード解説  A.2 変数の引渡し図  29}

/model/agcm/agcm5/doc/dai3bu/appndx/senzu.tex  2002/006/20(地球流体電脳倶楽部)

第B章謝辞

資源 agcm5

プログラム製作沼口敦

プログラム製作協力保坂征宏, 竹広真一, 石渡正樹, 堀之内武

文章保坂征宏, 竹広真一, 石渡正樹, 堀之内武, 余田成男

版権と引用本資源の版権は地球流体電脳倶楽部に属する. 地球流体電脳倶楽部の英語名は’GFD-DENNOU Club’ である. 資源の利用にあたっては地球流体電脳倶楽部の定める規定にしたがっていただきたい. 原則として, 教育的目的の場合には自由に使用・改変して良いものとしている. 詳細は心得を参照されたい. 地球流体電脳倶楽部の定める規定は地球流体電脳倶楽部サーバー (dennou) 上の

    ftp://dennou-k.gaia.h.kyoto-u.ac.jp/saloon/dennou/To.Users

というファイルに収められている. 電脳サーバーは次の 4 箇所におかれている.

    dennou-h.ees.hokudai.ac.jp
    dennou-t.ms.u-tokyo.ac.jp
    dennou-k.gaia.h.kyoto-u.ac.jp
    dennou-q.geo.kyushu-u.ac.jp

agmc5 を用いて作成された著作物等には, 英語では例えば “The calculations were performed by agcm5, GFD-DENNOU Club.” 日本語では例えば「計算は地球流体電脳倶楽部の agcm5 をもちいた.」のように明記していただければ幸いである. 「著者」に対応するものとしては地球流体電脳倶楽部 (GFD-DENNOU Club) をあげていただければ幸いである.

沿革 agcm5 は, 気象庁数値予報課の全球現業予報モデル (1985 年に住明正・増田耕一両氏が気象庁から東大理学部地球物理学科気象研究室に移植した) を参考に, 沼口敦氏が開発した大気大 \underline{AGCM5 第 3 部コード解説 31}

93/07/09 林祥介 97/03/11 余田成男

循環モデルである. agcm5 の 5 は氏の開発の第 5 版を意味する数字であり, 1988 年の作業開始以来, 5 回にわたる大きな変更を経たものである. この間, 地球流体電脳ライブラリ (dcl) との整合性をも考慮し, 大学等での研究・教育用に利用できる形に整理がなされてきた. 残念ながら, agcm5 は現実大気のシミュレーションモデルとしてのチューンがなされていないモデルであり, 若干の水過程と若干の放射過程を含むものの, その実体は単なる球面プリミティブモデルといったほうが正しい.

国立環境研究所, 並びに, 東京大学気候システム研究センターでは, 沼口敦氏が製作の中心となり agcm5 をベースとして気候シミュレーションモデル, CCSR/NIES モデルが作成された.

地球流体電脳倶楽部版 agcm5 は沼口氏によって提供された資源に, マニュアルの整備, インストール道具の整備, バグフィックスを行なったものである.

/model/agcm/agcm5/doc/dai3bu/CREDITS.tex  2002/006/20(地球流体電脳倶楽部)