====== 簡単な流出計算 ======

===== 合理式法 =====

===== 単位図法 =====

===== 単一流域の流出計算 =====

降雨流出過程を下図のように簡単にモデル化した場合の計算方法をプログラムを作成する。

{{ :fortran:4水理計算:1流出計算:model_4_1_1.jpg?nolink& |}}
;#;
**降雨流出過程の簡単モデル化（貯留型モデル）**
;#;

 \\

有効降雨<m>r</m>、雨水貯留<m>S</m>、直接流出<m>q</m>はすべて時間<m>t</m>の関数であり、<m>S</m>の時間変化は、

<m>{dS}/{dt}~=~r~-~q</m>　・・・(1)

で与えられる。また、<m>S</m>と<m>q</m>の間には次のような関係を仮定する。

<m>S=K ・ q^p</m>　・・・(2)

ここに、<m>K</m>と<m>p</m>は定数であり、流域の過去の雨量・流量のデータから決定される。

時刻<m>t</m>における初期流出<m>q_0=q</m>(<m>t_0</m>)と<m>t_0<=t<t_f</m>での有効降雨<m>r</m>(<m>t</m>) (<m>~>=0</m>)が与えられる。

<m>t_0<t_1<t_2< cdots cdots <t_n <= t_f</m>での直接流出量、<m>q</m>(<m>t_1</m>),<m>q</m>(<m>t_2</m>),<m>cdots cdots</m>,<m>q</m>(<m>t_n</m>)を求める。

計算で用いる降雨と定数値

<m>K=7.0,~ p=0.6,~ A=250.0</m>(<m>km^2</m>)<m>,~ t_0=0.0</m>(<m>h</m>)<m>,~ q_0=0.0</m>(<m>mm/h</m>)<m>,~ Delta t=1.0</m>(<m>h</m>)<m>,~ e=0.001,~ n=10</m>

降雨時系列
^時間 (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)式を差分化すると

<m>{s(t+Delta{t})-s(t)}/{Delta{t}}~=~{r(t+Delta{t})+r(t)}/{2}~-~{q(t+Delta{t})+q(t)}/{2}</m>　・・・(3)

次に(2)式を(3)式に代入すると、

<m>K/{Delta{t}}q(t+Delta{t})^p~+~1/2q(t+Delta{t})~-~K/{Delta{t}}q(t)^p~+~1/2q(t)~-~{r(t+Delta{t})+r(t)}/{2}~=~0</m>　・・・(4)

となる。式(4)で時刻<m>t</m>における未知数は<m>q</m>(<m>t+Delta{t}</m>)だけであり、この方程式を解くことにより<m>q</m>(<m>t+Delta{t}</m>)を求める。

簡単のため、式(4)で<m>x=q</m>(<m>t+Delta{t}</m>)と置くと、

<m>K/{Delta{t}}x^p~+~1/2x~-~K/{Delta{t}}q(t)^p~+~1/2q(t)~-~{r(t+Delta{t})+r(t)}/{2}~=~0</m>　・・・(5)

<m>p=1.0</m>の時は、1次方程式、<m>p=0.5</m>の時は<m>sqrt{x}</m>に関する2次方程式となるので、式(5)は代数的に解ける。

一般には、何らかの計算法を用いて、式(5)を近似的に満足するような解（近似解）を求める。ここでは、Newton-Raphson法((Newton-Raphson法については、数値解析等の参考書を参照してください。ここでは詳細は省きます。))を用いて解くことにする。

 \\

式(5)に**Newton-Raphson法**を適用するプログラムを考える。簡単のため、<m>t_0,t_1,cdots cdots t_n</m>は等時間間隔<m>Delta{t}</m>ごとの各時点を表すものとし、<m>t_j<=t<t_{j+1}</m>,(<m>j=0,1, cdots cdots cdots</m>)での平均有効降雨

<m>r_{j+1}=1/{Delta t}int{t_j}{t_{j+1}}{r(t)}{dt}</m>　・・・(6)

が与えられるものとする。このとき、<m>t_j<=t<t_{j+1}</m>における式(5)の最終項は、

<m>{r(t+Delta{t})+r(t)}/{2}~=~r_{j+1}</m>　・・・(7)

となる。各時間<m>t_j</m>において式(5)を解く際に必要な情報は、<m>K, p, Delta{t},q</m>(<m>t_j</m>)<m>,r_{j+1}</m>と収束条件<m>varepsilon</m>だけである。プログラム中では式(5)の右辺の第3項から第5項をまとめて変数**C**とおいている。

<note warning>
<fc #FF0000>**（単位に関する注意）**</fc> \\ 
降雨r<sub>j</sub>は降雨強度(mm/h)で与えられることが多い。したがって、流出qは流出高(mm/h)で求めておく。流出流量Q(m<sup>3</sup>/s)に変換する場合は、流域面積A(km<sup>2</sup>)を用いて、Q=q・A/3.6とする。
</note>

==== プログラム ====

プログラム中の変数

<box 80% |主要変数一覧>
  * k   :定数　式(2)
  * p   :定数　式(2)
  * a   :流域面積 (km<sup>2</sup>)
  * t0  :初期時刻t<sub>0</sub> (h)
  * q0  :初期流出高q<sub>0</sub> (mm/h)
  * dt  :単位時間間隔⊿t (h)
  * eps :収束判定規準e
  * n   :時間ステップ数n
  * r   :第jステップの降雨強度r<sub>j</sub> (mm/h)
  * q   :時刻t<sub>j</sub>の流出高q<sub>j</sub> (mm/h)
</box>

 \\ 
 \\ 
 \\ 
Newton-Raphson法による単一流域の降雨流出計算のソース            \\ 

<code fortran runoff_model.f90>
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                          !プログラム終了文

</code>

入力データ
<code>
----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+
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
----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+
</code>


データ数が多くなった場合は、入力データのファイルを作成して、プログラムにて読み込むように変更することも必要である。