【表計算ソフトウェアによる計算シート例を [別ページ]として追加しました。】
このページは未完成です。今後も内容を大きく書き変える可能性があります。
[温室効果の基本]のページで、 温室効果の強化によって起こる気候変化について、 エネルギー循環系かつフィードバックシステムであるシステムの簡単な例である [ふろおけモデル]と類似していることを述べた。
そして、「ふろおけモデル」のうちでも簡単な、 流入量が一定値であり、流出量がストックに比例する場合について、 [モデルの構成およびRによる計算プログラム]と [表計算ソフトウェアによる計算プログラム]の例を示した。
ここでは、気候モデルのうちでももっとも単純な0次元エネルギー収支モデル ([別ページ]参照)、 ただし定常を前提とせず時間変化がありうるものの簡単な例について、 モデルの構成およびRによる計算プログラムを示す。 [表計算ソフトウェアによる計算プログラム]の例は 別ページに示す。
気候システムがたくわえているエネルギー量を、便宜上、地球の表面積あたりにして、 S とする(SI単位は J/m2)。以下、便宜上「ストック」と呼ぶ。 すると、気候システムのエネルギー収支は次のように書ける。
dS/dt = Fin - Fout ... (1)
ここで、エネルギー収入 Finは次のように書ける。
Fin = (1/4) S0 (1 - αp) ... (2)
ただし、S0は「太陽定数」(1360 W/m2)、 αpは地球の「惑星アルベド」(0.30)で、 ここではひとまず、時間とともに変化することはないパラメータとみなす。 (パラメータ値をを変えてみることはありうる。)
エネルギー支出 Foutは地球の有効放射温度 Teと次のような関係にある (これは有効放射温度の定義)。ここで、σBはシュテファン・ボルツマンの定数(物理定数、 5.67×10-8 W m-2 K-4)。
Fout = σB Te4 ... (3)
また、温室効果の強さを g という文字を使う。同じ文字を使ってしまうが重力加速度ではない。 温室効果の強さの表現は標準がないが、この授業では、 2015年第2回小レポート課題で 「定常状態を仮定して、地表面温度 Tsが有効放射温度のなん倍か」としたのと ほぼ同じ意味で、定常でなくても使える定義として、 「地球から出ていく放射のエネルギーフラックス密度Foutが、 地表面温度の(1/g)倍の温度の黒体放射に等しい」となるような g、ということにする。 つまり、Foutと g とは次のような関係にある。
Fout = σB (Ts/g)4 ... (4)
(地表面積あたりの)エネルギーのストックは、 気候システムの代表温度に(地表面積あたりの)熱容量 C をかけたものであるとみなし、 気候システムの代表温度を地表面温度と同一視することにする。つまり
S = C Ts ... (5)
熱容量をもたらす実体としては海水の一定厚さの層を考えることにする。 地表面積あたりの熱容量は、密度ρ、比熱容量c、厚さhの積である。
C = ρ c h ... (6)
密度 (1.00×103 kg/m3)、比熱容量 (4186 J/(kg K)) は、 それぞれ定数とみなし、(純粋な)水に関する代表的な値(かっこ内)で近似する。 (古い単位の定義で、密度は 1 g/cm3、比熱容量は 1 cal/(g K)であることから換算。) 層の厚さをモデルのパラメータとする。 代表的な値は、季節変化に関与する海洋混合層を想定した場合 100 m、 風成循環(授業ではあとの回で述べる)に関与する海洋表層を想定した場合 500 m、 海洋全層を想定した場合 4000 mである。
以上を、Tsに関する式としてまとめると次のようになる。 (ただしFin、Cはそれぞれ定数とみなし、内わけには立ち入らないことにする)。
dTs/dt = (1/C) {Fin - σB (Ts/g)4} ... (7)
(7)はTsに関する微分方程式だが、右辺がTsの一次式ではないので、 [ふろおけモデルの簡単な例]のように 指数関数の形の簡単な解析解が見つかるわけではない。 しかし、Tsが、ある一定値Ts0からあまり大きく離れない値をとるとすれば、 次の線形近似(一次関数近似)を使うことができ、その解は指数関数を使った形に書ける。
Ts4 ≒ Ts04 + 4 (Ts - Ts0)3 ... (8)
【この節は、式の番号以外は、[ふろおけモデルの簡単な例]と同じだが重ねて示しておく。】
モデルが少し複雑になると、解析解が得られなくなる。 その場合に有効なのは、数値計算によって近似解を求めることである。 有限の時間区間(時間差Δt)の間のSの変化 ΔS = S(t+Δt) - S(t) を考える。 dS/dtは ΔS/Δt の Δt→0 での極限だが、これを、極限をとらないΔS/Δtで近似する。 たとえば (1) の近似は次のようになる。
ΔS/Δt = Fin - Fout ... (9)
これを変形して、S(t)がわかっているときにS(t+Δt)を求める式が得られる。
S(t+Δt) = S(t) + ( Fin - Fout ) Δt ... (10)
いまFinは定数としたが、Foutは時間とともに変化する。 (10)式を使おうとすると、時間区間のうちのどの時刻でのFoutの値を使うかが問題になる。 ここでは、区間の初めの時刻 t での値で代表させることにする。
【微分を有限区間の差分で近似する方法にはいろいろなものがある。 気象学で応用される方法の概略は、別の教材ページ 「数値計算の教材: 時間発展問題の差分計算」で紹介した。 区間の中央で代表させる方法もあり、うまく使えればそのほうが精度がよいのだが、 数値があばれること(「計算不安定」という)をふせぐのにくふうがいるので、 ここでは、区間の初めで代表させる(専門用語で「前方差分」という)方法に限って紹介する。】
Δtが大きすぎると、差分による計算結果は微分方程式の解の近似として不適切になる。 Δtを小さくすると、計算量がふえて、計算機に対する負荷が重くなるほか、誤差の累積が問題になる場合もある。 適切なΔtを選ぶには、試行錯誤を必要とすることもあるが、 まず重要なのは、起こる現象の時間規模よりもやや細かくとることである。 たとえば、exp(-a t)に比例するような指数関数型の変化が起こるならば、Δtは 1/a よりも小さく(たとえばその5分の1くらいに)する。
【この節の前おき部分の文章の大部分は[ふろおけモデルの簡単な例]と同じだが、重ねて示しておく。】
わたしが教材として示すプログラム例は、 百行くらいまでの分量のソースプログラム(人が読み書きする形でのプログラム)を 必要に応じて部分的に書きかえて実行することができる人向けである。 プログラミングをしない人に使ってもらうために配慮したプログラムではない。 また、計算のプロフェッショナルの観点ですぐれたプログラムでもない。
わたしはこれまで、数値計算のプログラムをおもにFortran言語で書いてきた。 「数値計算の教材: 時間発展問題の差分計算」の例文は、Fortran, Pascal, C言語によるものを用意した。
ここでは、R言語によるものを整備したいと考えている。 Rの処理系は、農工大のパソコン教室のWindowsマシンにはすでにはいっているし、 オープンソースソフトウェアとして提供され、 MS Windows, Mac, Linuxをはじめとする多くのOSの計算機上で動くので、 Rに慣れることは多くの進路の人にとって損にならないと思う。 Rは統計的データ解析のための言語であり、 ここで行なうのは数値シミュレーションの基礎となる理論的計算なので、 Rの典型的な使いかたではないのだが、なんとか使うことができる。
なお、Rを含むほとんどのプログラム言語は、数値を扱うことはできるが、 単位のついた物理量を扱う機能をもっていない。 物理量の次元と単位についてはプログラミングの前の段階で人が注意して扱わなければならない。
教材用のプログラムを作るにあたって、次のような方針をとることにする。
プログラム例を、下の「ソースプログラム」という見出しのあとに用意した。 その部分だけの内容をテキストファイルとして、たとえば greenhouse_ex1.R という名まえで保存したとする。 R処理系を起動し、 Windows用の処理系の日本語メニューで言えば「ファイル」の中の「Rコードのソースを読み込み」で そのファイルを指定すれば、実行されるだろう。 しかしむしろ、「ディレクトリの変更」をしたあとで、「R console」のコマンドとして(「>」のプロンプトから) 「source "greenhouse_ex1.R"」のような形で実行するのがよいだろう。 プログラムの編集は、Rの外のテキストエディタ (Windowsの「メモ帳」を含む)ですることもできるが、 メニューの「ファイル」の中の「スクリプトを開く」からRのエディタを使って「保存」するのが便利だろう。
プログラムのどこを変えたら結果がどう変わったかを確認するには、 プログラムの複数のバージョンを比較できる形で持つ必要がある。 エディタで単純に保存すると上書きになってしまうので、 変更前のプログラムをコピーしておいてファイル名で区別するなどのくふうを、各自してほしい。 (バージョン管理を支援するソフトウェアを使うとよいのかもしれないが、わたしは(まだ)使っていない。)
Rのプログラム中のコメントや文字定数などに日本語文字を使ってもかまわないはずだが、 文字コードの問題にぶつかるおそれがあるので、例文中の文字はASCIIに限り、 コメントは(日本語ローマ字でもよいのだが、ここでは)英語でつけた。
R言語では各行の「#」という字からあとはコメントであり計算処理上は無視される。
例文のプログラムを部分に分けて説明するために、コメントに(1)(2)のような番号を入れた。(式の番号ではない。)
【R言語で「<-」は代入をあらわす記号である。右側の式の計算結果を左側の変数にいれる。】
【R言語で「^」は指数演算の記号である。「a^b」で abをあらわす。】
(1a) 物理定数。シュテファン・ボルツマンの定数 σB (変数名 stefanc)の数値を SI単位で与える。
(1b) 地球全体に対するパラメータ。 「太陽定数」S0 (変数名 solconst)、 「惑星アルベド」αp (変数名 albedo)。 そして、気候システムのエネルギー収入 Fin (定数扱いで値を保持するための変数名 inflowc)と 有効放射温度 Te (変数名 tempeq)を計算している。
(1c) 温室効果の強さのパラメータ g (変数名 ghfac)。 例示としては、地表面温度が現在の実際の世界平均値(288 K)に近くなるような値 (1.13)を与えておく。
(1d) 熱容量に関するパラメータ。 水の密度 ρ (変数名 denswater)、水の比熱容量 c (変数名 sphcwater)、 層の厚さ h (単位 m、変数名 depthlayer)を与えて、 地表面積あたりの熱容量 (変数名 hcaplayer)を計算している。 層の厚さの例示としては、海洋表層を想定した 500 m をとった。
差分の1ステップの時間間隔Δt (単位を年として与えるための変数名 deltatyear、 秒単位に換算した値のための変数名 deltat)を5年、 シミュレーションのステップ数(変数名 nstep)を20 (つまり継続時間を100年)とする。 (シミュレーションの考えかたのうえでは継続時間から考えて割り算でステップ数を決めるべきだが、 割り算の結果が整数にならないとやっかいなので、この手順にした。) また、初期(t = 0)での地表面温度Ts(0) (変数名 tempinit) を 280 Kとしてみる。 (これを Te と同じにするとつまらない結果になるので、 それからいくらか離れた値を入れてみる。)
時間とともに変化する変数(および、この例では一定だが、別の問題では変化させたい、流入量inflow)に ついて、それぞれ、(nstep+1)個の記憶場所を確保する。 この手続きの書きかたは、プログラム言語によって違う。
シミュレーションの本体。まず地表面温度の初期値を与えたあと、 for文によるループで、時間をΔtずつ進めながら、 nstep回、次の時間ステップでの地表面温度を計算する。
【この部分のプログラムは、R言語の「for(k in 1:length(time))」のループで、 kの値を1から順にひとつずつ length(time) [この例文では21] まで変えて それぞれ実行されることを前提としている。 もし並列に、あるいは順不同に実行される場合は 期待どおりの動作をしない。 その場合は確実に順に実行されるようにプログラムを変更する必要がある。】
ここでは、地表面温度(temp)の時系列を表示する手続きを用意した。 エネルギー流出量の時系列を表示する場合は、 ソースプログラムを編集して、「temp」を「outflow」に変更してほしい。
数値を数字で表示する手続きも用意したほうがよいのだが、 プログラム実行直後に、Rの処理系内のコマンド行として(プロンプト「>」のあとに) 「temp」や「outflow」などの変数名を単独で打てば、 その変数の値が表示される (配列の場合には列挙される)ので、 この例題の結果を確認することに限ったとりあえずの目的にはそれで充分かと思う。
そのまま実行するほかに、(1), (2)の部分を次のように変更して結果を比較してみてほしい。
【ひとつのプログラムの中で複数のパラメータセットによる シミュレーションを実行し、結果を比較できる形で表示できるようにしたいが、 そのためのプログラムの変更のしかたの説明は、まだ用意していない。追って用意したい。】
# 0-dim. energy balance model
# Simple case: greenhouse effect is represented by a single value
# R program 2015-06-09 MASUDA Kooiti
##### (1) conditions of the model
##### ##### (1a) physical constant
stefanc <- 5.67e-8
##### ##### (1b) parameters of the earth as a whole
solconst <- 1360
albedo <- 0.30
inflowc <- solconst * (1 - albedo) / 4 # inflow rate, constant (W/m2)
tempeq <- (inflowc / stefanc) ^ (1/4)
##### ##### (1c) parameter of greenhouse effect
ghfac <- 1.13
##### ##### (1d) parameters of the heat capacity of the climate system
denswater <- 1.00e3 # density of water (kg/m3)
sphcwater <- 4186 # specific heat capacity of water (J/(kg K))
depthlayer <- 500 # depth of layer (m)
hcaplayer <- denswater * sphcwater * depthlayer
# heat capacity of water layer per unit area (J/(m2 K))
##### (2) conditions of a run of simulation
tempinit <- 280 # initial level of temperature (K)
yearsec <- 86400*365.25 # unit conversion factor year/s
deltatyear <- 5 # time interval of a step (year)
deltat <- deltatyear*yearsec # time interval of a step (s)
nstep <- 20 # number of time steps
##### (3) preparing array variables
time <- seq(0, deltatyear*nstep, deltatyear) # time (year)
temp <- time * 0 # temperature (K)
inflow <- time * 0 # inflow rate (W/m2)
outflow <- time * 0 # outflow rate (W/m2)
##### (4) time-stepping simulation
temp[1] <- tempinit
for(k in 1:length(time)){
inflow[k] <- inflowc
outflow[k] <- stefanc * (temp[k]/ghfac)^4
if(k < length(time)) {
temp[k+1] <- temp[k] + deltat * (inflow[k] - outflow[k]) / hcaplayer
}
}
##### (5) time-series graph of temperature
plot(time, temp, tck=1)