【[2019-02-01] データファイルを用意することについて補足しました。】

気候システム論

0次元エネルギー収支モデルのRによる計算プログラム

時間微分の差分近似

時間微分を含むモデルが少し複雑になると、 解析解(数式としての解)が得られなくなる。 その場合に有効なのは、数値計算によって近似解を求めることである。 有限の時間区間(時間差Δt)の間のSの変化 ΔS = S(t+Δt) - S(t) を考える。 dS/dtは ΔS/Δt の Δt→0 での極限だが、これを、極限をとらないΔS/Δtで近似する。 たとえば

dS/dt = F

の近似は次のようになる。

ΔS/Δt = F

これを変形して、S(t)がわかっているときにS(t+Δt)を求める式が得られる。

S(t+Δt) = S(t) + F Δt

この式を使おうとすると、時間区間のうちのどの時刻でのFの値を使うかが問題になる。 ここでは、区間の初めの時刻 t での値で代表させることにする。

【微分を有限区間の差分で近似する方法にはいろいろなものがある。 気象学で応用される方法の概略は、別の教材ページ 「数値計算の教材: 時間発展問題の差分計算」で紹介した。 区間の中央で代表させる方法もあり、うまく使えればそのほうが精度がよいのだが、 数値があばれること(「計算不安定」という)をふせぐのにくふうがいるので、 ここでは、区間の初めで代表させる(専門用語で「前方差分」という)方法に限って紹介する。】

Δtが大きすぎると、差分による計算結果は微分方程式の解の近似として不適切になる。 Δtを小さくすると、計算量がふえて、計算機に対する負荷が重くなるほか、誤差の累積が問題になる場合もある。 適切なΔtを選ぶには、試行錯誤を必要とすることもあるが、 まず重要なのは、起こる現象の時間規模よりもやや細かくとることである。 たとえば、exp(-a t)に比例するような指数関数型の変化が起こるならば、Δtは 1/a よりも小さく(たとえばその5分の1くらいに)する。


R言語によるプログラム例

わたしが教材として示すプログラム例は、 百行くらいまでの分量のソースプログラム(人が読み書きする形でのプログラム)を 必要に応じて部分的に書きかえて実行することができる人向けである。 プログラミングをしない人に使ってもらうために配慮したプログラムではない。 また、計算のプロフェッショナルの観点ですぐれたプログラムでもない。

わたしはこれまで、数値計算のプログラムをおもにFortran言語で書いてきた。 「数値計算の教材: 時間発展問題の差分計算」の例文は、Fortran, Pascal, C言語によるものを用意した。

ここでは、R言語によるものを整備したいと考えている。 Rの処理系は、 オープンソースソフトウェアとして提供され、 MS Windows, Mac, Linuxをはじめとする多くのOSの計算機上で動くので、 Rに慣れることは多くの進路の人にとって損にならないと思う。 【農工大にパソコン教室があってOSがWindowsのマシンがはいっていたときには そこに常備されていた。】 Rは統計的データ解析のための言語であり、 ここで行なうのは数値シミュレーションの基礎となる理論的計算なので、 Rの典型的な使いかたではないのだが、なんとか使うことができる。

なお、Rを含むほとんどのプログラム言語は、数値を扱うことはできるが、 単位のついた物理量を扱う機能をもっていない。 物理量の次元と単位についてはプログラミングの前の段階で人が注意して扱わなければならない。

教材用のプログラムを作るにあたって、次のような方針をとることにする。

プログラム例を、下の「ソースプログラム」という見出しのあとに用意した。 その部分だけの内容をテキストファイルとして、たとえば ebex.R という名まえで保存したとする。 R処理系を起動し、 Windows用の処理系の日本語メニューで言えば「ファイル」の中の「Rコードのソースを読みこみ」で そのファイルを指定すれば、実行されるだろう。 しかしむしろ、「ディレクトリの変更」をしたあとで、「R console」のコマンドとして(「>」のプロンプトから) 「source "ebex.R"」のような形で実行するのがよいだろう。 プログラムの編集は、Rの外のテキストエディタ (Windowsの「メモ帳」を含む)ですることもできるが、 メニューの「ファイル」の中の「スクリプトを開く」からRのエディタで編集して「保存」するのが便利だろう。

プログラムのどこを変えたら結果がどう変わったかを確認するには、 プログラムの複数のバージョンを比較できる形で持つ必要がある。 エディタで単純に保存すると上書きになってしまうので、 変更前のプログラムをコピーしておいてファイル名で区別するなどのくふうを、各自してほしい。 (バージョン管理を支援するソフトウェアを使うとよいのかもしれないが、わたしは(まだ)使っていない。)

Rのプログラム中のコメントや文字定数などに日本語文字を使ってもかまわないはずだが、 文字コードの問題にぶつかるおそれがあるので、例文中の文字はASCIIに限り、 コメントは(日本語ローマ字でもよいのだが、ここでは)英語でつけた。

R言語では各行の「#」という字からあとはコメントであり計算処理上は無視される。

例文のプログラムを部分に分けて説明するために、コメントに(1)(2)のような番号を入れた。(式の番号ではない。)

(1) モデルのパラメータ

【R言語で「<-」は代入をあらわす記号である。右側の式の計算結果を左側の変数にいれる。】

(1a) 熱容量に関するパラメータ。 水の密度 ρ (変数名 denswater)、水の比熱容量 c (変数名 sphcwater)、 層の厚さ h (単位 m、変数名 depthlayer)を与えて、 地表面積あたりの熱容量 (変数名 hcaplayer)を計算している。 層の厚さの例示としては、海洋表層を想定した 500 m をとった。

(1b) 気候システムとその「外」との熱交換に関するパラメータ。 熱交換は気候システムの温度と「外の温度」との差に比例するとしたので、 その比例定数(問題文では「k」)を変数 coef とした。単位は (W/m2)/K。

(1c) 温室効果に関するパラメータ。 二酸化炭素倍増に対する定常応答を変数名 tdouble (単位は Kまたは℃)とし、 数値は 3.0 Kとした。 産業革命前(気候システムが定常状態にあったと仮定する)の二酸化炭素濃度を 変数名 co2preind (単位は ppm)とし、数値は 280 ppm とした。

(2) シミュレーション実行のパラメータ

差分の1ステップの時間間隔Δt (単位を年として与えるための変数名 deltatyear、 秒単位に換算した値のための変数名 deltat)を10年、 シミュレーションのステップ数(変数名 nstep)は、 観測された二酸化炭素濃度にもとづく計算が1850年から2000年までの16ステップ、 RCPシナリオの二酸化炭素濃度にもとづく計算が2010年から2100年までの10ステップである。 (「ステップ数」は時間区間の数をさすのと区間の切れめの時点の数をさすのとで 意味がすこしちがうが、ここでは時点のほうをかぞえている。) (なお、シミュレーションの考えかたのうえでは、 まずシミュレーションの対象とする期間とΔtを考えて割り算でステップ数を決めるべきだが、 割り算の結果が整数にならないとやっかいなので、整数であたえることにした。)

初期(t = 0)は、西暦1850年とし、yearinitにはそれを西暦の年であたえてある。

初期(t = 0)での温度 (変数名 tempinit) は、ひとまず 0 があたえてある。 今回の問題では、初期の温度がなんであっても、 それを基準とした温度差は同じなので、 実際の温度のほうがわかりやすければ、 その初期値を K あるいは ℃ 単位で入れても、ほかの部分を変える必要はない。 (他の問題、とくに黒体放射の式がかかわる問題や、 水の相変化がかかわる問題では、温度の絶対値が必要になる。)

(3) データの読みこみ

R言語の data frame という機能を使い、read.table 関数で読みこむことにした。 データフレーム変数 concobs に観測された二酸化炭素濃度、 concrcp にRCPシナリオの二酸化炭素濃度を読みこんでおく。

【[2019-02-01補足] プログラム実行まえに、 データファイル ar5_co2conc_obs.txt と ar5_co2conc_rcp.txt を、 [期末レポート問題]のページに示したリンク先から ダウンロードして、ebex.R と同じディレクトリに置いてあることを想定している。】

(4) 配列変数を確保する

時間とともに変化する変数について、それぞれ、(nstep+1)個の記憶場所を確保する。 この手続きの書きかたは、プログラム言語によって違う。

(5) 差分法による計算

シミュレーションの本体。まず温度の初期値 temp[1] を与えたあと、 for文によるループで、時間をΔtずつ進めながら、 nstep回、次の時間ステップでの温度 temp[k] を計算する。

【この部分のプログラムは、R言語の「for(k in 1:nstep)」のループで、 kの値を1から順にひとつずつ nstep まで変えて それぞれ実行されることを前提としている。 もし並列に、あるいは順不同に実行される場合は 期待どおりの動作をしない。 その場合は確実に順に実行されるようにプログラムを変更する必要がある。】

k が nstepobs (=16)以下ならば、観測にもとづく二酸化炭素濃度、 そのさきは、RCP 4.5シナリオにもとづく二酸化炭素濃度を使うようにしてある。 「RCP4_5」と書いてあるところを変えれば、別のRCPシナリオを使うことができる。

二酸化炭素濃度から、その濃度への定常応答を計算し、 このモデルの「外の温度」(変数名 tempx)とみなす。

単位時間あたりのエネルギー供給量 F は、coef * (tempx[k] - temp[k]) である。

温度に、deltat * F をかけたものをたし、つぎのステップの温度をもとめる。 (ただし、nstepに達したらこの計算はしない。)

(6) 結果を数値で出力する

ステップ番号、時刻(年)、温度、「外の温度」を横にならべて書きだすようにしてみた。

出力さきファイル名は仮に "ebex_RCP4_5.txt" としてあるが、適宜、変えてほしい。

ファイルへの書き出しは「追加」(append)にしてある。 出力さきファイル名を変えずにくりかえし実行すると、 前回の結果のあとにつけたされるので、注意してほしい。

書式指定は最小限なので、文字数はそろわない。 数値のあいだは空白がはいるはず。

(7) 結果を図示する

時刻(年)を横軸、温度を縦軸にとったグラフをかくようにしてある。 ただし、軸の目盛りや説明ラベルなどの指定は、まだできていない。


ソースプログラム

# ebex.R
# 0-dim. energy balance model
# Simple case: energy supply is like external "heat bath"
# R program 2019-01-23 MASUDA Kooiti

##### (1) conditions of the model
##### ##### (1a) parameters of the heat capacity of the climate system
denswater  <-  1.00e3 # density of water (kg/m3)
sphcwater  <-  4.20e3 # 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)

##### ##### (1b) parameter of heat exchange
coef       <-  1.5    # ((W/m2)/K)

##### ##### (1c) parameters of greenhouse effect
tempdouble <-  3.0    # steady-state response to CO2 doubling (K)
concpreind <-  280    # concentration of CO2 before industrialization (ppm)

##### (2) conditions of a run of simulation
tempinit   <-     0         # initial level of temperature (K)
                            # (relative values are used in this example)
yearsec    <-  86400*365.25 # unit conversion factor year/s
yearinit   <-  1850
deltatyear <-    10         # time interval of a step (year)
deltat     <-  deltatyear*yearsec # time interval of a step (s)
nstepobs   <-    16   # number of time steps with "obs" ((2000-1850)/10 +1)
nsteprcp   <-    10   # number of time steps with RCP   ((2100-2010)/10 +1)
nstep      <-  nstepobs + nsteprcp

##### (3) read data
concobs <- read.table("ar5_co2conc_obs.txt", skip=5, header=TRUE)
concrcp <- read.table("ar5_co2conc_rcp.txt", skip=5, header=TRUE)

##### (4) preparing array variables
time     <- seq(yearinit, yearinit+deltatyear*(nstep-1), deltatyear) # time (year)
temp     <- time * 0  # temperature (K)
tempx    <- time * 0  # "external" temperature (K)

##### (5) time-stepping simulation
temp[1] <- tempinit
for(k in 1:nstep){
  if(k <= nstepobs) {
    concnow  <- concobs[concobs$year == time[k], ]$co2conc
  } else {
    concnow  <- concrcp[concrcp$year == time[k], ]$RCP4_5
  }
  tempx[k] <- tempinit + tempdouble*log(concnow/concpreind)/log(2)
  if(k < nstep) {
    temp[k+1] <- temp[k] + deltat * coef * (tempx[k] - temp[k]) / hcaplayer
  }
}

##### (6) print
for(k in 1:nstep){
  resultline <- sprintf("%d %d %f %f\n", k, time[k], temp[k], tempx[k])
  cat(resultline, file="ebex_RCP4_5.txt", append=TRUE)
}

##### (7) time-series graph of temperature
plot(time, temp, tck=1)

2019-01-23, 2019-02-01
増田 耕一 (MASUDA Kooiti)
[2018年度期末レポート問題へ] [授業日程と内容のページ] [2018年度授業の目次ページ]