[気候システム全体の(0次元)エネルギー収支]のページで、 0次元(緯度・経度・高さのいずれの方向の不均一も省略する)、 しかも定常(時間による変化がないとする)の、 気候システムのエネルギー収支を考え、 大気の効果がないばあいの気候システムの温度が「有効放射温度」であるとのべた。
ここではまず、0次元のままで、定常は仮定せず、 時間とともに気候システムのもつエネルギーが変化し、 それとともに温度が変化する問題について、計算プログラムを、 R言語と、表計算(ExcelまたはLibreOffice)のワークシートとで表現してみる。
便宜上、単位面積あたりのエネルギーのたまりと、 単位面積・単位時間あたりのエネルギーの流れについて考える。 地球(あるいは他の惑星)の形を球とする近似で、球の表面積あたりで考えるのだ (断面積あたりではない)。
エネルギー収支の式は、つぎのように書ける。 【分数の記号は横線であらわしたほうがよいのだが、 HTMLのつごうで便宜上斜線であらわしている。 ノートをとる際には、横線で書きなおすことをおすすめする。】
d(C T)/dt = Fin - Fout ... (1)
ここで、t は時間 (SI単位は秒)、 C は、ここで考える気候システムの単位面積あたりの熱容量 (SI単位は J/(m2 K) )、Tは温度 (SI単位は K)である。 気候システムのもつエネルギーは、単純に C かける T とかけると仮定しよう。 はじめから定常状態を考えるならば、Cの値はいくつでも答えは変わらないのだが、 時間変化のあるばあいは、Cの値によって結果がかわってくる。 Cは、モデルの実行ごとにちがう値をとりうるが、実行中には変化しない数量だと仮定しよう (このような性格の数量を、モデルのパラメータ(parameter)という。)
式(1)は微分方程式であり、Fin、Foutが 時間tの関数または温度Tの関数としてあつかいやすい形をしていれば積分でとけるかもしれないが、 多くの応用の場面では、数値計算で近似解をもとめることになる。 なにかの量 X の時間による微分 dX/dt は、有限の時間区間 Δ[デルタ]t での X の変化をΔtでわったもの 「(X(t+Δt) - X(t)) / Δt」の、Δtをゼロにちかづけた極限と考えることができる。 近似計算では、極限をとらずに、時間ステップ Δt での変化を考えることが多い。 このような考えかたを「差分法」という。 (1)を差分近似すると、つぎのようになる。 ただし、Cは変化しない量だとしたので、差分の外にくくりだした。 【かけざんの記号は省略するのがふつうだが、ここでは意識的に「×」とかいておく。】
C × {T(t＋Δt) - T(t)} / Δt = Fin - Fout ... (2)
ここで、右辺がいつの時点の数量であるかが、まだ指定されていない。 本来は、時刻t から時刻 t+Δt までの時間区間での数量だが、 近似として、時間区間のはじめの時点である時刻 t での数量とみなすことにする。 (このような近似を、差分法の分類のうえで「前方差分」ということがある。 [時間発展問題の数値計算 -- 差分法]のページで説明した。) この仮定のもとで、式(2)を、T(t)がわかっているとき、T(t+Δt)をもとめる式に変形する。
T(t+Δt) = T(t) + (1/C) × Δt × { Fin(t) - Fout(t) } ... (3)
ここで、時刻 t = 0 での T の値 T(0) 【これを「Tの初期値」という】がわかっていれば、つぎのように時間ステップを進めて、 原理的にはいくらでも長い時間の変化をシミュレートすることができる。 (ただし、差分近似による誤差が蓄積して、微分方程式の解から離れていく。 また、差分のやりかたによっては計算結果があばれることがある。)
T(Δt) = T(0) + (1/C) × Δt × { Fin(0) - Fout(0) }
T(2Δt) = T(Δt) + (1/C) × Δt × { Fin(Δt) - Fout(Δt) }
T(3Δt) = T(2Δt) + (1/C) × Δt × { Fin(2Δt) - Fout(2Δt) } ... (4)
Fin、Foutのきめかたは、ひとまず、 [気候システム全体の(0次元)エネルギー収支]のページの 例題と同じにしてみよう。
Fin = (1/4) × {1 - αp} × S ... (5)
ここで S は、太陽と地球を結ぶ直線に垂直な平面の単位面積あたりに、単位時間あたりに到達する 太陽放射のエネルギー (SI単位は W/m2)、いわゆる「太陽定数」である。 この太陽放射を、地球(大気をふくむ)は吸収するか反射するかのいずれかに分配する。 α[アルファ]pは、このうち反射するほうの割合で、「惑星アルベド」とよばれることがある。 4分の1は、球の断面積と表面積との比である。 Sとαpは、変量ではあるが、 モデルのひとつの実行のあいだは変化しない定数とみなすことにしよう。 (Sとαpを「モデルのパラメータ」としたことになる。)
Fout = σB T4 ... (6)
ひとまず、気候システムのエネルギー支出は、温度 T での黒体放射だと考えると、 この式になる。 σ[シグマ]Bは、シュテファン・ボルツマンの定数という物理定数で、 有効数字3けたで、5.67×10-8 W m-2 K-4 である。
これで、物理定数とモデルパラメータと時間間隔Δtの数値と、 初期値T(0)があたえられれば、(4)のように時間ステップを追った計算ができるはずだ。 R言語と、表計算ワークシートで、それぞれやってみることにしよう。
モデルパラメータのうち、面積あたりの熱容量 C は、 考える対象によってもちがってくる。 ここでは、まず、数十年から百年の気候変化を想定することにする。 そのばあいは、海洋の表層 (風成循環がおよぶ層) の熱容量が重要であることがわかっている。 それを単純化して、深さ h = 500 m の水の層を考えることにしよう。 液体の水の比熱容量 (単位質量の水を単位温度だけ上げるのに必要なエネルギー)は cwater = 4.2 × 103 J / (kg K)、 液体の水の密度 (単位体積あたりの質量)は ρ[ロー]water = 1.0 × 103 kg / m3、 いずれも定数とみなすことにしよう。 単位面積あたりの熱容量は、C = cwater ρwater h である。
計算プログラムを組むまでに、もうひとつ考えておく必要があるのは、 計算機上の記憶領域(プログラム言語での変数、表計算でのワークシート)の使いかたである。 3次元の気候モデルでは、1つの時点での気候システムの状態をあらわすだけでも多数の数値をつかうので、 それをさらに多数の時間ステップのぶんだけ記憶したまま計算をつづけることはできない。 計算機の主記憶装置上には、 現在計算中の時間ステップの値と、そのひとつ前の時間ステップの値だけを持ち、 過去の値があとで必要になるならば、計算できしだいファイルに書き出しておくのがふつうである。
0次元モデルは、あまり記憶領域をくわないので、 全部の時間ステップでの気候システムの状態の量を持ちつづけることができる。 例文はこの形でつくっている。 できれば(別の機会に)、計算できしだい書き出す方式も紹介したい。
結果をグラフにするプログラム、結果を数値としてファイルに書きだすプログラム、 そのファイルを読んでグラフにするプログラムの例をふくめた。
こちらは、ワークシート中の計算式を入れたセルで数値が得られるところまで。 グラフ化については、各自、表計算ソフトウェアのグラフ機能 (そのうちではおそらく「散布図」機能をつかうのがよい)、 または、そこから書き出して他の道具をつかってほしい。 その方法の教材は(わたしは)用意しない。
例題プログラムをそのまま動かすことができたら、 つぎには、温度の初期値を 50℃ から、ほかのなんとおりかの値に変えて、 結果がどう変わるか、見てほしい。 さらにできれば、層の深さを 500 m からほかの値に変えてみてほしい。 ただし、その場合は、時間刻みを 1年からほかの値に変える必要があるかもしれない。