大気中の温室効果気体(二酸化炭素など)の濃度が大きいほど 地上気温あるいは地表面温度が高くなることを、 準定常状態どうしの比較として説明するための、 単純化したモデルがある。
これは、真鍋(1985)の図3として説明されているものである。
実際の気候システムにくらべて、次のような単純化がされている。
この単純化のもとで、温室効果気体がふえることは、 地球放射を出す代表高さが、高いほう(宇宙空間に近いほう)にずれることだと考えられる。 それは、同じ地表面温度に対して、地球放射が出ていきにくくなることなので、 ふろおけモデルの、出口がつまりぎみになることに似ている。
【ふろおけモデルの場合と同様、 コメントやグラフに表示される語句に、日本語が含まれている。】
R言語では各行の「#」という字からあとはコメントであり計算処理上は無視される。
R言語で「<-」は代入をあらわす記号である。右側の式の計算結果を左側の変数にいれる。
例文のプログラムを部分に分けて説明するために、コメントに(1)(2)のような番号を入れた。(式の番号ではない。)
Version 1とVersion 2とのちがいは、(1)で熱容量をあたえている部分、(2)のdeltatyearの値、 (6)のタイトル文字列、(6a)と(6b)の出力ファイル名の文字列だけである。
いわゆる太陽定数 (変数名solconst)、 気候システムの太陽光反射率(惑星アルベド、変数名 albedo)を与える。 気候システムの、地球の表面積あたり・単位時間あたりのエネルギー収入は、 inflowc = solconst * (1 - albedo) / 4 となる (単位は W/m2)。
地球放射を射出する代表高さを heightemit{old, new} (単位 m、例として 5000 m と 5500 m)、 対流圏の気温減率を lapserate (プログラム内での単位は K / m、6.5 K/kmから換算してあたえている)、 Stefan-Boltzmannの定数を sigmasb (5.67×10-8 W/(m2 K4) )とする。
気候システムの単位面積あたりの熱容量を capac (単位は J/(m2 K) )とする。 それをどうあたえるかは、下の別項目でのべる。
ここでついでに、モデルの定常状態の地上気温(=地表面温度) tempsteady{old, new} を計算しておく。
差分の1ステップの時間間隔Δt を、年単位 (変数名deltatyear)であたえ、 秒単位 (変数名 deltat) に換算する。1年を 365.25 日とみなした。 Δtの実際の値のあたえかたは、熱容量とあわせてのべる。 シミュレーションのステップ数(変数名 nstep)を 40 としてみた。 (シミュレーションの考えかたのうえでは継続時間から考えて割り算でステップ数を決めるべきだが、 割り算の結果が整数にならないとやっかいなので、この手順にした。) また、初期(t = 0)での地上気温 (変数名 tempinit) には、旧状態の定常状態の地上気温 tempsteadyold の値を与える。
時間とともに変化する変数(および、この例では一定だが、別の問題では変化させたい、流入量inflow)に ついて、それぞれ、(nstep+1)個の記憶場所を確保する。 この手続きの書きかたは、プログラム言語によって違う。
シミュレーションの本体。まずエネルギーのたまりの量の初期値を与えたあと、 for文によるループで、時間をΔtずつ進めながら、 nstep回、次の時間ステップでのエネルギーのたまりの量を計算する。
【この部分のプログラムは、R言語の「for(k in 1:length(time))」のループで、 kの値を1から順にひとつずつ length(time) [この例文では21] まで変えて それぞれ実行されることを前提としている。 もし並列に、あるいは順不同に実行される場合は 期待どおりの動作をしない。 その場合は確実に順に実行されるようにプログラムを変更する必要がある。】
結果の図のうち地上気温のグラフにあわせて表示するために、定常状態の地上気温を配列変数に入れておく。 【ふろおけモデルのばあいとちがって、時間発展する変数の理論値ではない。】
結果をグラフにする。ただし、画面表示ではなく、PNG画像ファイルで書き出すことにした。 気候システムの(単位面積・単位時間あたり)エネルギー収入・支出の時系列のグラフと、 地上気温の時系列のグラフとを、2つの別々の画像ファイルにする。 計算結果を記号(まるや三角)、定常状態を折れ線で表示することにする。
「太陽定数」、射出高さ、(単位面積あたり)熱容量の値を、Rのplot機能の「subtitle」として図の下に書くようにした。 現状の例文では、それと横軸のラベルとの間がせますぎる。 対策がわかったら例文を改訂するつもりだが、ひとまずこのまま使ってみてほしい。
大気の対流圏の熱容量を概算であたえた。 大気の対流圏の上端の気圧を 200 hPa、下端の気圧を 1000 hPa、重力加速度を 9.8 m/s2とし、 静水圧のつりあいを仮定して、単位面積あたりの質量をもとめる。
大気の体積は固定されていないので、圧力一定での比熱容量を使う。概算で 1000 J/(kg K)とした。
結果が約1年でほぼ定常におちつくと予想されるので、Δtは (1/40)年にした。
海洋の表層(風成循環のおよぶ層)の熱容量を概算であたえた。 概算で、海洋の表層の深さを 500 m、水の密度を 1 g/cm3、水の比熱容量を 4.2 J/(g K)とし、 SI単位に換算してあたえるようにした。
結果が数十年でほぼ定常におちつくと予想されるので、Δtは 2年にした。
Version 1, 2 では、出ていく地球放射の代表高さ heightemit を突然変えたあとは一定としたが、 Version 3 では、これが徐々に変化することを考える。 プログラム上では heightemit を配列変数にしたうえで、 (1)の部分でその初期値 heightemitinit を与え、 (2)の部分で、1年あたりの代表高さの変化量 heightemitchangerate を与えて、 (4)の k を順に変えたループの中で heightemit[k] の値を決めている。 なお tempsteady も変化するので、配列変数に値を入れていく。
作図は、エネルギーの流れの量のグラフ、温度のグラフのほか、heightemitのグラフもつくる。
自分の実習用のディレクトリ(フォルダ)を用意したうえで、 そこに、プログラム例 ondanka_m2.R をダウンロードして、保存してほしい。
そのまま実行するほかに、(1), (2)の部分を次のように変更して結果を比較し、特徴をのべてほしい。
自分の実習用のディレクトリ(フォルダ)を用意したうえで、 そこに、プログラム例 ondanka_m3.R をダウンロードして、保存してほしい。
そのまま実行するほかに、(1), (2)の部分を、たとえば次のように変更して結果を比較し、特徴をのべてほしい。
注意については、ふろおけモデルの説明を参照してほしい。
プログラム言語 Python (パイソン)で Version 2 と同じ計算をする例文を追加した。
Windowsでは、PythonのIDLEというエディタに読みこんで、 (操作メニューが英語ならば)「Run」から「Run module」で実行できるはず。
文字コードは UTF-8である。
残念ながら、結果を出力するところがまだしっかりできておらず、 配列変数の内容をそのまま数字で書き出すだけになっている。
そのほか、ふろおけモデルのページの Pythonによる例題と同様なので そちらを参照してほしい。
今回の問題の計算は、表計算ソフトウェア(Microsoft ExcelやLibreOffice Calc)のワークシートでも可能であるが、 ふろおけモデルのページに書いたような理由で、おすすめしない。