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

気候システム論

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

[0次元エネルギー収支モデルのRによる計算プログラム]と 基本的に同じ計算をする、Python (パイソン)言語によるプログラムをつくってみました。

ソースプログラムは ebex.py です。

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

Rの例文の(7)にあたるグラフ出力はありません。 また、Rの例文の(6)ではファイル出力にしましたが、 いまの段階でのPythonの例文の結果は、標準出力に書き出されます。 Pythonを「IDLE」から起動した場合は、結果は「Shell」ウィンドウに出ます。 Shellウィンドウから、エディタ(「メモ帳」でもよい)に cut and paste して、 テキストファイルとして保存するのがよいでしょう。 (わたしがファイル出力の方法を理解したら例文をなおそうと思いますが、 いつまでにできると約束できません。)

RCPシナリオの選択は、 ixrcp という変数の値によることにしました。 ixrcpとシナリオとの対応は、 1:RCP2.6, 2:RCP4.5, 3:RCP6.0, 4:RCP8.5 です。

言語仕様のちがいで、Rでは配列変数の要素番号が 1 からはじまりますが、 Pythonではリスト変数の要素番号が 0 からはじまります。 そのため、Python の例文での変数 k の値は、R の例文での k よりも 1 だけ小さい値をとります。


ソースプログラム

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

##### (0) use Python library (function "log" is in the "math" library)
from math import log

##### (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
##### ##### initialize array variables
years = []
conc  = []
##### ##### read observed CO2 concentration
klineskip = 6
kline = 0 # initial
with open('ar5_co2conc_obs.txt') as f:
  for f_line in f:
    kline = kline + 1
    if (kline > klineskip):
      sss = f_line.strip().split()
      years.append(int(sss[0]))
      conc.append(float(sss[1]))
nkobs = kline - klineskip
##### ##### read CO2 concentration of RCP scenarios
ixrcp = 1  # 1:RCP2.6, 2:RCP4.5, 3:RCP6.0, 4:RCP8.5
klineskip = 6
kline = 0 # initial
with open('ar5_co2conc_rcp.txt') as f:
  for f_line in f:
    kline = kline + 1
    if (kline > klineskip):
      sss = f_line.strip().split()
      years.append(int(sss[0]))
      conc.append(float(sss[ixrcp]))
nkrcp = kline - klineskip
nkall = nkobs + nkrcp
##### ##### output for checking
## for k in range(nkall):
##   print (k, years[k], conc[k])

##### (4) preparing array variables
time   = [] # time (year)
temp   = [] # temperature (K)
tempx  = [] # "external" temperature (K)

##### (5) time-stepping simulation
# temp[0] = ...
temp.append(tempinit)
for k in range(nstep): # k = 0, ..., k = nstep-1
  year = yearinit + deltatyear * k
# time[k] = ...
  time.append(year)
  concnow = conc[ years.index(year) ]
# tempx[k] = ...
  tempx.append(tempinit + tempdouble*log(concnow/concpreind)/log(2))
# temp[k+1] = ...
  if(k < nstep-1):
    temp.append(temp[k] + deltat * coef * (tempx[k] - temp[k]) / hcaplayer)

##### (6) print
for k in range(nstep):
   print(k, time[k], temp[k], tempx[k])

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