降雨流出過程を下図のように簡単にモデル化した場合の計算方法をプログラムを作成する。
降雨流出過程の簡単モデル化(貯留型モデル)
有効降雨、雨水貯留、直接流出はすべて時間の関数であり、の時間変化は、
・・・(1)
で与えられる。また、との間には次のような関係を仮定する。
・・・(2)
ここに、とは定数であり、流域の過去の雨量・流量のデータから決定される。
時刻における初期流出()とでの有効降雨() ()が与えられる。
での直接流出量、(),(),,()を求める。
計算で用いる降雨と定数値
()()()()
降雨時系列
時間 (h) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
---|---|---|---|---|---|---|---|---|---|---|
降雨量(mm/h) | 3.5 | 5.4 | 12.0 | 7.2 | 9.7 | 5.3 | 1.2 | 0.0 | 0.0 | 0.0 |
(1)式を差分化すると
・・・(3)
次に(2)式を(3)式に代入すると、
・・・(4)
となる。式(4)で時刻における未知数は()だけであり、この方程式を解くことにより()を求める。
簡単のため、式(4)で()と置くと、
・・・(5)
の時は、1次方程式、の時はに関する2次方程式となるので、式(5)は代数的に解ける。
一般には、何らかの計算法を用いて、式(5)を近似的に満足するような解(近似解)を求める。ここでは、Newton-Raphson法1)を用いて解くことにする。
式(5)にNewton-Raphson法を適用するプログラムを考える。簡単のため、は等時間間隔ごとの各時点を表すものとし、,()での平均有効降雨
・・・(6)
が与えられるものとする。このとき、における式(5)の最終項は、
・・・(7)
となる。各時間において式(5)を解く際に必要な情報は、()と収束条件だけである。プログラム中では式(5)の右辺の第3項から第5項をまとめて変数Cとおいている。
プログラム中の変数
主要変数一覧
Newton-Raphson法による単一流域の降雨流出計算のソース
program runoff_model !プログラムの宣言 implicit none real(8) k, p, a, t0, q0, dt, eps, t, qq, p !倍精度実数型変数の宣言 real(8) c, xi, xi1, fx, fdx !倍精度実数型変数の宣言 real(8) r(50), q(50) !倍精度実数型配列の宣言 integer n, j !整数型変数の宣言 read(*,*) k, p, a !標準入力からの読み込み read(*,*) t0, q0 !標準入力からの読み込み read(*,*) dt, eps !標準入力からの読み込み read(*,*) n !標準入力からの読み込み read(*,*) (r(j), j = 1, n) !標準入力からの読み込み t = t0 !初期時刻の設定 r(0) = 0.0d0 !初期降雨強度の設定 write(*,*) 'time r q qq' ! q0 = q0 * a / 3.6d0 !初期流出高の流量へ変換 write(*,*) t, r(0), q0, qq ! !計算開始 do j = 1, n !doループによりデータ数だけ計算する t = t + dt !時刻dtだけ進める c = k * q0 ** p / dt -q0 / 2.0d0 + r(j) !定数項Cの計算 xi1 = 0.001 !初期解を0.001に設定する 10 xi = xi1 !f(x)を求めるための準備 fx = k * xi ** p / dt + xi / 2.0d0 - c !f(x)を求める if ( abs(fx) < eps ) then !収束判定のためのif文 q(j) = xi !収束したならば、配列qに格納 qq = xi * a / 3.6d0 !流量に変換 write(*,*) t, r(j), q(j), qq !結果の出力 else !収束していない場合 fdx = k * p * xi ** (p - 1.0d0) / dt + 0.5 !fdxを計算 xi1 = xi - fx / fdx !新たにxi1を計算 goto 10 !文番号10へ飛び、反復計算を繰り返す end if !if文の終了 q0 = q(j) !求められたq(j)を次のステップの初期値とする end do !doループ終了 end program runoff_model !プログラム終了文
入力データ
----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+ 7.0 0.6 250 0 0 1 0.001 10 3.5 5.4 12.0 7.2 9.7 5.3 1.2 0.0 0.0 0.0 ----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+
データ数が多くなった場合は、入力データのファイルを作成して、プログラムにて読み込むように変更することも必要である。