表計算ソフトウェアのワークシートファイルは、複数のワークシートからなる「ブック」であることが多くなったが、 ここでは、ひとつのワークシート内で計算する。 ただし、モデルの中で時間とともに変化する量と、 変化しない量 (定数およびモデルパラメータ) がある。 ワークシートを、 左のほう(A-D列)は変化しない量、 右のほう(E列から)は変化する量を置く場所として、区別してつかうことにする。
ワークシートの縦の列には左からA, B, C, ...、 横の行には上から 1, 2, 3, ... のように番号がついていて、 セル(升目)は、たとえば「A1」のような番号でよぶことができる。 セルの中には、計算式を書くことができる。 計算式は「=」で始まり、定数やセル番号をふくむ演算式がつづく。 演算式中では、かけざんは「*」、わりざんは「/」であらわす。 計算式の中でセル番号を参照する際には、「絶対参照」と「相対参照」がある。 絶対参照は、行、列それぞれの番号のまえに「$」をつける。相対参照はつけない。 このちがいの効果は、計算式をコピーしたばあいにあらわれる。 コピーすると、式の中の相対参照さきは、新しい位置に相対的につけかわるが、 絶対参照さきは変わらない。 わたしのプログラムでは、時間とともに変化しない量は絶対参照し、変化する量は相対参照する。
【セルにラベル(名まえ)をつけて、それによって参照することもできるようだが、 わたしがまだその方法を理解していないので、セル番号による参照という従来の方法をとっている。 ワークシートには、[Rプログラムによる計算]でつかう「変数名」も入れてあるが、 いまのところワークシートではそれを使っていない。 この変数名をラベルとし、それを計算式の中でセル番号のかわりにつかうことができれば、プログラムの見通しがだいぶよくなるだろう。】
数量は原則として SI 単位であらわすが、「年」や「℃」などの慣用の単位をつかうこともある。 その場合は、秒やKなどのSI単位とのあいだの換算をおこなう。 ワークシートに書きこんだ単位は、人が見るためのもので、演算では使っていない。 単位のうちで、べき乗をあらわす上つき添え字を、上つきでない通常の字で代用している。
例文のワークシートは、MS Windows上の LibreOffice Calc バージョン 7.0.0 で動作確認した。 Microsoft Excelではまだ動作確認していないが、たぶん動くと思う。
この例題で、数値を変えることを想定しているモデルパラメータは、「層の深さ」である。 これを変えると、「時間刻み(年)」も変える必要があるかもしれない。 層の深さを小さくすると面積あたりの熱容量も小さくなり、 温度の変化が速くなるので、時間刻みを短くする必要があるだろう。 シミュレーション対象となる時間は、時間刻みとステップ数の積である。 (ただし、現在のワークシートでは、ワークシート上の「ステップ数」の値は使っておらず、 実際に計算するステップ数は「変わる量の部分」の行を何行コピーするかできめている。) 「初期時刻(年)」は、この例題では 0 のまま変えないのがよいだろう。 (別の例で、たとえば西暦2020年からのシミュレーションをすることを考えて、変えられるようにしておいた。)
また、モデルパラメータの「太陽定数」と「惑星アルベド」も、 例題としては一定値を与えたが、変更しうる量である。
| A | B | C | D | |
|---|---|---|---|---|
| 1 | シミュレーションの例題1: 0次元エネルギー収支気候モデル、時間変化あり、大気の効果なし (増田 2020-08-09版) | |||
| 2 | 定数とモデルパラメータ | |||
| 3 | 意味 | 変数名 | 数値 | 単位 |
| 4 | 初期時刻 | tinit | =$C$5*$C$8 | 秒 |
| 5 | 初期時刻(年) | tinityear | 0 | 年 |
| 6 | 時間刻み | deltat | =$C$7*$C$8 | 秒 |
| 7 | 時間刻み(年) | deltatyear | 1 | 年 |
| 8 | 1年の秒数 | yearsecond | =60*60*24*365.25 | 秒 |
| 9 | ステップ数 | nstep | 200 | 無次元整数 |
| 10 | ステファン・ボルツマン定数 | sigmab | 5.67e-8 | W m-2 K-4 |
| 11 | 水の凝固温度 | tzero | 273.15 | K |
| 12 | 「太陽定数」 | solconst | 1.36e3 | W m-2 |
| 13 | 惑星アルベド | alphap | 0.3 | 無次元量 |
| 14 | ||||
| 15 | 面積あたり熱容量 | capac | =$C$16*$C$17*$C$18 | J m-2 K-1 |
| 16 | 層の深さ | depth | 500 | m |
| 17 | 密度 | dens | 1.0e3 | kg m-3 |
| 18 | 比熱容量 | sphcapac | 4.2e3 | J kg-1 K-1 |
「↓」は、その上のセルにある計算式をコピーしてはりつけることをさす。
モデルの初期値を、時間ステップ番号0 での「温度(℃)」の値として、 セル J6 にあたえることにした。 例として、50℃という値をあたえている。この値を変えたシミュレーションをためしてほしい。 なお、初期値を、K [ケルビン] 単位の温度としてあたえたいならば、 セル I6 に数値をあたえ、セル J6 には セル J7 にある計算式をコピーすればよい。
数値書式を、G, H列は指数表現(小数点下6けた)、I..L列は小数表現(小数点下2けた)とするのがよい。
結果を、時間軸を「年」として表示したければ、 「時刻(年)」の列をふやして、6行めのセルに =$C$4+F6*$C$7 を入れ、 その計算式を下にコピーしておけばよいだろう。 ただし、ワークシートの行や列の数を変えるときには、 計算式の意味が変わらないように、計算式から参照されるセル番号を点検する必要がある。
| E | F | G | H | I | J | K | L | |
|---|---|---|---|---|---|---|---|---|
| 2 | 時間発展する変数 | |||||||
| 3 | 意味 | ステップ番号 | 時刻 | エネルギーたまり | 温度 | 温度(℃) | エネルギー収入 | エネルギー支出 |
| 4 | 変数名 | kt | tim | estock | temp | tdegc | fin | fout |
| 5 | 単位 | 無次元整数 | 秒 | J m-2 | K | ℃ | W m-2 | W m-2 |
| 6 | 0 | =$C$4+F6*$C$6 | =I6*$C$15 | =J6+$C$11 | 50.00 | =$C$12*(1-$C$13)/4 | =$C$10*(I6^4) | |
| 7 | =F6+1 | ↓ | =H6+$C$6*(K6-L6) | =H7/$C$15 | =I7-$C$11 | ↓ | ↓ | |
| 8 | ↓ | ↓ | ↓ | ↓ | ↓ | ↓ | ↓ | |
| ... | ||||||||
| 206 | ↓ | ↓ | ↓ | ↓ | ↓ | ↓ | ↓ | |
この例題で、数値を変えることを想定しているモデルパラメータは、 「層の深さ」depth ($C$16) の値である。 これを変えると、「時間刻み(年)」 deltatyear ($C$7) の値も変える必要があるかもしれない。 層の深さを小さくすると面積あたりの熱容量も小さくなり、 温度の変化が速くなるので、時間刻みを短くする必要があるだろう。 シミュレーション対象となる時間は、時間刻みとステップ数 nstep ($C$9) の積である (ただしこのワークシートでは実際には nstep の値をつかわず、計算式を何行コピーするかできめている)。 「初期時刻(年)」tinityear ($C$5) は、この例題では 0 のまま変えないのがよいだろう。 (別の例で、たとえば西暦2020年からのシミュレーションをすることを考えて、変えられるようにしておいた。)
また、モデルパラメータの「太陽定数」solconst ($C$12) と「惑星アルベド」alphap ($C$13) も、 例題としては一定値を与えたが、変更しうる量である。
時間ステップごとの計算の内容は、エネルギー保存の法則にもとづいて、 単位面積あたりのエネルギーのたまりの量の変化をもとめているのだが、 便宜上、「エネルギーのたまりの量」の変数をもちださずに、 温度の変化をもとめる形に書いた。 このモデルでは、単位面積あたりのエネルギーのたまりの量は、 温度 (I列) に、単位面積あたりの熱容量 capac ($C$15) をかけたものであり、 capac は時間とともに変わらないと仮定したから、このような単純化ができた。 (エネルギーのたまりの量がほしければ、I列の値 に capac をかければよい。) 熱容量が時間とともにかわりうるようなモデルでは、 このような単純化ができず、 温度とエネルギーのたまりの量のそれぞれの変数をもつ必要があるだろう。
まず、モデルパラメータを変更せずに、温度の初期値 ($J$6 に入れる数値)を なんとおりか変えて、結果を比較してみよう。
つぎに、層の深さ depth ($C$16) の数値を、なんとおりか変えて、 初期値のほうは変えないで計算してみよう。 ただし、層の深さを小さくすると、時間刻み (「年」単位) の値 deltatyear ($C$7) も 小さくする必要があるかもしれない。