====== tys-sfm_v10 ======

tys-sfm_v10.lzhをダウンロードして解凍すると、「<fc #FF0000>**貯留関数マクロ.xls**</fc>」ができます。
この「<fc #FF0000>**貯留関数マクロ.xls**</fc>」の使用方法を説明します。
ダウンロードについては、[[http://www.disaster-software.net|水理水文ソフトウェア]]のHPを参照してください。


===== マクロの説明 =====

本マクロは、各条件シートに条件を入力するだけで、<fs medium><fc #FF0000>**１流域についての貯留関数による流域流出解析計算**</fc></fs>が行えます。

  - 実績流量より、貯留関数係数の回帰を行う。
  - 任意の貯留関数係数を用いて流出解析を行う。
の２種類の計算を行うことができます。

==== 計算条件とデータの入力 ====

=== 計算条件 ===

「<fc #008080>計算条件</fc>」シートで貯留関数係数の回帰を行うかどうかをチェックしてください。


{{ :流域流出計算:image01.jpg?650|計算条件シート}}
;#;
計算条件シート
;#;


<note important>
再度計算実行した場合は、計算結果シート（複数）の値がすべて上書きされます。
</note>
=== 雨量・流量データ ===

{{ :流域流出計算:image03.jpg?650|雨量・流量データシート}}
;#;
雨量・流量データシート
;#;
\\

  - **実績流量より、貯留関数係数の回帰を行う場合**　「<fc #008080>計算条件</fc>」シートの流域面積・基底流量・流路長さ・遅滞時間計算ピッチを入力してください。　「<fc #008080>雨量・流量データ</fc>」シートの流域平均雨量・実績流量・計算対象範囲指定のデータを入力してください。__計算対象範囲指定は入力した実績流量のうち、係数決定に用いるデータの範囲を指定します。開始点を「1」、終了点を「2」と入力してください。__
  - **任意の貯留関数係数を用いて流出解析を行う場合** 「<fc #008080>計算条件</fc>」シートの流域面積・基底流量・貯留関数係数（Ｋ・Ｐ値、遅滞時間、流出係数）を入力してください。「<fc #008080>雨量・流量データ</fc>」シートの流域平均雨量を入力してください。実績流量は、計算には使用しませんが、雨量・流量結果図で計算結果の評価を行うために入力する方が望ましいでしょう。
==== 計算結果 ====
  - 実績流量より、貯留関数係数の回帰を行う場合
    - 「<fc #008080>計算結果</fc>」シートに遅滞時間の設定値ごとの貯留量と流出量の関係の計算結果を書き込みます。この中で、貯留量と流出量のループの積算値の最小値を遅滞時間として採用します。尚、遅滞時間の繰り返し計算は、０時間～最大値（流路長を木村の式に適用した場合の遅滞時間＋１）を遅滞時間の計算ピッチ分です。ここで求めた遅滞時間より実績流量を用いた最小自乗法による回帰を行い、Ｋ・Ｐを算出します。
    - 「<fc #008080>雨量・流量データ</fc>」シートの計算流量の列に上記で求めた貯留関数係数を用いて計算した流量を書き込みます。計算結果は「<fc #008080>雨量・流量結果図</fc>」に反映されます。
    - 「<fc #008080>貯留量と流出量の関係</fc>」シートに決定した貯留関数係数を用いた貯留量と流出量の関係を出力します。その結果は「<fc #008080>貯留量と流出量の関係図</fc>」シートに反映されます。
  -任意の貯留関数係数を用いて流出解析を行う場合
    - 「<fc #008080>雨量・流量データ</fc>」シートの計算流量の列に入力した貯留関数係数を用いて計算した流量を書き込みます。計算結果は「<fc #008080>雨量・流量結果図</fc>」シートに反映されます。
    - 「<fc #008080>貯留量と流出量の関係</fc>」シートに貯留量と流出量の関係を出力します。その結果は「<fc #008080>貯留量と流出量の関係図</fc>」シートに反映されます。
\\

{{ :流域流出計算:image02.jpg?650|計算結果}}
;#;
計算結果
;#;

{{ :流域流出計算:image04.jpg?650|雨量・流量結果図}}
;#;
雨量・流量結果図
;#;

{{ :流域流出計算:image05.jpg?650|貯留量と流出量の関係図}}
;#;
貯留量と流出量の関係図
;#;


==== 注意事項 ====
　本マクロは、通常の貯留関数法についての計算用であり、総合貯留関数法には適用していません。通常の貯留関数法と総合貯留関数法では、Ｋ・Ｐ値はそれぞれ全く異なりますので注意してください。
　ここで、通常の貯留関数法とは、貯留関数を貯留量(m3)と流出量(m3/s)の関係で表現したもの、総合貯留関数法とは貯留関数を貯留高(mm)と流出高(mm/hr)の関係で表現したものです。総合貯留関数法は貯留関数法を単位流域面積で表現したものと言い換えることができます。
　本や論文によっては、__総合貯留関数法を単に「貯留関数法」__と表記されている場合がありますので、基礎式による確認が必要です。

===== VBAソースファイル =====

　本マクロのソースを公開しますので、参考にしてください。

<sxh vb tys-sfm_v10.bas>
Option Explicit

Public m As Integer     'ハイドロデータ数
Public m1 As Integer    'ハイドロ区間１
Public m2 As Long       'ハイドロ区間２
Public l As Double      '流域の流路延長(km)
Public a As Double      '流域面積(km2)
Public qi As Double     '基底流量(m3/s)
'
Public Const mnum As Integer = 1000  '配列数
Public r(0 To mnum) As Double       '雨量(mm/h)
Public q(0 To mnum) As Double       '実測流量(m3/s)

Public fr(0 To mnum) As Double      '雨量を流量に換算した値
Public fi(0 To mnum) As Double      '雨量を流量に換算した値の積算値

Public ql(0 To mnum) As Double      '遅滞時間を考慮した流量(基底流量を差し引いたもの)
Public qs(0 To mnum) As Double      '遅滞時間を考慮した流量の積分値(基底流量を差し引いたもの)
Public s(0 To mnum) As Double
Public sa(0 To 10000) As Double     'ループの積算値(時計回りか反時計回りかを判定)
'
Public rt As Double         '対象範囲の累積雨量(mm)
Public f As Double          '流出係数
Public tlk As Double        '遅滞時間(設定値)
Public qlm As Double        '流量の最大値(基底流量を差し引いたもの)
Public tpl As Double        '流量が最大の時間(遅滞時間分差し引いたもの)
Public qtm As Double        '流量の積算(遅滞時間分差し引いたもの)
Public qsm As Double        'Qsの積算
Public fim As Double        'Fiの最大値
Public sam As Double        'ループの最大
Public sas As Double        'ループの最小
Public tld As Double        '遅滞時間（決定分）
Public tl0 As Integer       '遅滞時間（切り上げ）
Public tl1 As Integer       '遅滞時間（整数部）
Public tl2 As Double        '遅滞時間小数部
Public t_Step As Double     '遅滞時間の計算ステップ
Public num_Cal As Integer    '遅滞時間計算回数

'
Public e As Double          '最小自乗法の対象となるデータの個数
Public sm As Double         '貯留高さの最大値
Public sar As Double        '上昇ループの積算値
Public sad As Double        '下降ループの積算値
'
Public k As Double          'Ｋ値
Public p As Double          'Ｐ値

Public Sub Main()
'貯留関数法１流域のみ、係数を回帰で決定
    Dim ws1 As Worksheet
    Dim ws2 As Worksheet
    Dim ws3 As Worksheet
    Dim ws4 As Worksheet
    Dim i As Integer
    Dim tl As Integer
    Dim t As Integer
    Set ws1 = Sheets("計算条件")
    Set ws2 = Sheets("雨量・流量データ")
    Set ws3 = Sheets("計算結果")
    Set ws4 = Sheets("貯留量と流出量の関係")
    'データ読み込み
    DataRead1 ws1
    DataRead2 ws2
    
    If Sheet1.CheckBox1.Value = False Then GoTo CalSkip '係数回帰をするかどうかのチェック
    
    '遅滞時間の初期定義（木村の式）
    If l < 11.9 Then
        tl = 0
    Else
        tl = Int(0.5 + 0.047 * l - 0.56)   '切り上げ
    End If
    '適当な遅滞時間の判別
    '**********************
        num_Cal = CInt(CDbl(tl + 1) / t_Step)
    '**********************
    For t = 0 To num_Cal
        tlk = t_Step * CDbl(t)  '遅滞時間設定値
        tl1 = Int(tlk)          '遅滞時間整数部
        tl2 = tlk - tl1         '遅滞時間少数部
        If tl2 = 0 Then
            tl0 = tl1
        Else
            tl0 = tl1 + 1
        End If
    
        CalQl
        
        'DataWrite0 T   'デバック
        
        CalQtm
        f = qtm / ((a / 3.6) * rt)
        CalFr
        CalSa t
        '
        DataWrite1 ws3, t   '計算経過出力
        '
    Next t
    'ループの最大・最小判定
    sam = 0
    For t = 0 To num_Cal
        If Abs(sa(t)) >= sam Then sam = Abs(sa(t))
    Next t
    sas = sam
    For t = 0 To num_Cal
        If Abs(sa(t)) <= sas Then
            sas = Abs(sa(t))
            tld = t_Step * CDbl(t)  '遅滞時間設定値
        End If
    Next t
    tl1 = Int(tld)          '遅滞時間整数部
    tl2 = tld - tl1         '遅滞時間少数部
    If tl2 = 0 Then
        tl0 = tl1
    Else
        tl0 = tl1 + 1
    End If
    
    CalQl
    
    CalQtm
    
    f = qtm / ((a / 3.6) * rt)  '流出係数
    
    fim = 0
    For t = 1 To m
        fr(t) = a * f * r(t) / 3.6
        fim = fim + fr(t)
        fi(t) = fim
    Next t
    
    qsm = 0#
    For t = 2 To m - tl0
        qsm = qsm + (ql(t - 1) + ql(t)) / 2#
        qs(t) = qsm
        If fi(t) - qs(t) >= 0 Then s(t) = fi(t) - qs(t)
    Next t
    
    sm = 0
    For t = 2 To m - tl0 - 1
        If s(t) >= sm Then sm = s(t)
    Next t
    
    Cal_KP  '最小自乗法によるK・Pの計算
    
    DataWrite2 ws1      '貯留関数式書き込み
    DataWrite3 ws4      '貯留量と流出高書き込み
    Graph_Refresh       'グラフ修正
    
CalSkip:    '流量計算
    
    DataRead3 ws1
        
    Cal_Qt ws2          '求めた係数で流量を計算
        
    Graph_Refresh2       'グラフ修正
    
End Sub

Private Sub DataRead1(ws As Worksheet)
'データ読み込み１（流域諸元）
    ws.Select
    a = ws.Cells(3, "D")        '流域面積(km2)
    qi = ws.Cells(4, "D")       '基底流量(m3/s)
    l = ws.Cells(5, "D")        '流路長(km)
    t_Step = ws.Cells(6, "D")   '遅滞時間計算ピッチ(hr)
    If t_Step = 0 Then t_Step = 1
End Sub

Private Sub DataRead2(ws As Worksheet)
'データ読み込み１（ハイドロ・ハイエトデータ）
    Dim row_num As Integer
    Dim i As Integer
    '①ハイドロ・ハイエト読み込み
    row_num = 3
    i = 0
    Do Until ws.Cells(row_num, "C") = ""    'データがあるまでループ
        i = i + 1
        r(i) = ws.Cells(row_num, "C")
        q(i) = ws.Cells(row_num, "D")
        row_num = row_num + 1
    Loop
    m = i       'ハイドロ・ハイエトデータ数定義
    '計算対象範囲指定
    For i = 1 To m
        If ws.Cells(i + 2, "E") = 1 Then m1 = i         '開始点
        If ws.Cells(i + 2, "E") = 2 Then m2 = i         '終了点
    Next i
    '期間の累積雨量
    rt = 0
    For i = m1 To m2
        rt = rt + r(i)
    Next i
End Sub

Private Sub CalFr()
    Dim t As Integer
    fim = 0
    For t = m1 To m2
        fr(t) = a * f * r(t) / 3.6
        fim = fim + fr(t)
        fi(t) = fim
    Next t
End Sub

Private Sub CalQl()
'遅滞時間分戻した流量の定義（基底流量含まず）
    Dim t As Integer
    qlm = 0
    For t = 1 To m2
        If t - tl1 >= 0 Then
            ql(t - tl1) = q(t) + (q(t + 1) - q(t)) * tl2 - qi  '遅滞時間を考慮した流量(基底流量を差し引いたもの)
            If ql(t - tl1) >= qlm Then
                qlm = ql(t - tl1)    '最大流量
                tpl = t - tlk        '最大流量の発生時間
            End If
        End If
    Next t
End Sub

Private Sub CalQtm()
'遅滞時間を考慮した流量の積算（計算区間内）
    Dim t As Integer
    qtm = 0
    For t = m1 To m2 - tl0
        qtm = qtm + ql(t)
    Next t
End Sub

Private Sub CalSa(tl As Integer)
'貯留高さと流出高をプロットした場合のループ形状
    Dim t As Integer
    qsm = 0: sar = 0: sad = 0
    
    For t = m1 To m2 - tl0 - 1
        qsm = qsm + (ql(t - 1) + ql(t)) / 2#
        qs(t) = qsm
        If fi(t) - qs(t) < 0 Then
            s(t) = 0#
        Else
            s(t) = fi(t) - qs(t)
        End If
    Next t
    '上昇ループの積算値計算
    For t = m1 To tpl - 1
        sar = sar + (s(t) + s(t + 1)) * Abs(ql(t + 1) - ql(t)) / 2#
    Next t
    '下降ループの積算値計算
    For t = tpl To m2 - tl0 - 1
        sad = sad - (s(t) + s(t + 1)) * Abs(ql(t) - ql(t + 1)) / 2#
    Next t
    'ループ全体の積算値計算
    sa(tl) = sar + sad
End Sub

Private Sub Cal_KP()
'最小自乗法による貯留関数係数の計算
    Dim M3 As Integer
    Dim t As Integer
    Dim Xi As Double
    Dim Xx As Double
    Dim Yi As Double
    Dim XY As Double
    Dim P0 As Double
    Dim K0 As Double
    Dim K1 As Double
    'M3 = Int(0.8 * M)
    M3 = m2 - tl0
    e = 0
    'For T = 2 To M3
    For t = m1 To M3
        If s(t) <> 0 Then e = e + 1
    Next t
    Xi = 0: Xx = 0: Yi = 0: XY = 0
    'For T = 2 To M3
    For t = m1 To M3
        Xi = Xi + Log(ql(t))
        Xx = Xx + Log(ql(t)) * Log(ql(t))
        Yi = Yi + Log(s(t))
        XY = XY + Log(ql(t)) * Log(s(t))
    Next t
    P0 = (Xi * Yi - e * XY) / (Xi * Xi - e * Xx)
    K0 = (Yi - P0 * Xi) / e
    K1 = Exp(K0)
    k = K1
    p = P0
    k = Int(0.5 + K1 * 100#) / 100#     '桁合わせ
    p = Int(0.5 + P0 * 1000#) / 1000#   '桁合わせ
End Sub
Private Sub DataWrite1(ws As Worksheet, t As Integer)
'結果書き込み１（遅滞時間計算経過）
    Dim Sa_State As String
    If t = 0 Then
        ws.Cells(2, "A") = "遅滞時間(hr)"
        ws.Cells(2, "B") = "最大流量(m3/s)"
        ws.Cells(2, "C") = "最大流量時間(hr)"
        ws.Cells(2, "D") = "流量積算値(m3/s)"
        ws.Cells(2, "E") = "流出係数"
        ws.Cells(2, "F") = "上昇ループの積算値"
        ws.Cells(2, "G") = "下降ループの積算値"
        ws.Cells(2, "H") = "ループ全体の積算値"
        ws.Cells(2, "I") = "ループの状態"
    End If
    '
    Sa_State = ""
    If sa(t) > 0 Then Sa_State = "時計回り"
    If sa(t) < 0 Then Sa_State = "反時計回り"
    If sa(t) = 0 Then Sa_State = "ぴったりOK"
    '
    ws.Cells(t + 3, "A") = tlk
    ws.Cells(t + 3, "B") = qlm
    ws.Cells(t + 3, "C") = tpl
    ws.Cells(t + 3, "D") = qtm
    ws.Cells(t + 3, "E") = f
    ws.Cells(t + 3, "F") = sar
    ws.Cells(t + 3, "G") = sad
    ws.Cells(t + 3, "H") = sa(t)
    ws.Cells(t + 3, "I") = Sa_State
    
End Sub

Private Sub DataWrite2(ws As Worksheet)
'結果書き込み２（貯留関数式出力）
    ws.Select
    ws.Cells(10, "B") = k
    ws.Cells(11, "B") = p
    ws.Cells(12, "B") = tld
    ws.Cells(13, "B") = f
End Sub

Private Sub DataWrite3(ws As Worksheet)
'貯留量と流出高書き込み
    Dim t As Integer
    ws.Select
    ws.Cells(3, "D") = tld
    For t = 1 To m - tl0
        ws.Cells(t + 2, "A") = ql(t)
        ws.Cells(t + 2, "B") = s(t)
    Next t
End Sub

Private Sub DataRead3(ws As Worksheet)
'データ読み込み3（貯留関数式元）
    ws.Select
    k = ws.Cells(10, "B")
    p = ws.Cells(11, "B")
    tld = ws.Cells(12, "B")
    f = ws.Cells(13, "B")
End Sub

Private Sub Cal_Qt(ws As Worksheet)
'求めた貯留関数を用いて流量を再度計算
    Dim t As Integer
    Dim tm As Integer
    Dim mr As Integer   '雨量・流量データ数
    Dim qql(0 To mnum) As Double
    Dim qqt(0 To mnum) As Double
    Dim j As Integer
    Dim c As Double
    Dim d As Double
    Dim x As Double
    Dim x2 As Double
    Dim y As Double
    mr = m
    tl1 = Int(tld)          '遅滞時間整数部
    tl2 = tld - tl1         '遅滞時間少数部
    'If Tl2 <> 0 Then Tl1 = Tl1 + 1
     'ニュートン法による流量計算
    qql(1) = 0
    For t = 1 To mr
        c = 2# * f * a * r(t) / 3.6 - qql(t) + 2# * k * qql(t) ^ p
        x = 1
        j = j + 1
Step1:  y = x + (2# * k * x ^ p) - c
        d = 1# + 2# * p * k * x ^ (p - 1#)
        x2 = x - y / d
        If Abs(x2 - x) <= 0.0001 Then GoTo Step2
        x = x2
        GoTo Step1
Step2:  qql(t + 1) = x
        If t <= tl1 Then qqt(t) = qi
        qqt(t + tl1) = qi + qql(t + 1) + (qql(t + 1) - qql(t)) * tl2
        If j <= 100 Then GoTo Step3
        If qqt(t) < 1.1 * qqt(99) Then GoTo Skip1
Step3:
    Next t
Skip1:
    For t = 1 To mr
        ws.Cells(t + 2, "G") = qqt(t)
    Next t
End Sub

Private Sub Graph_Refresh()
'グラフの修正
    Dim g_range As String
    Sheets("貯留量と流出量の関係図").Select
    g_range = "A3:B" & CStr(m - tl0 + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheets("貯留量と流出量の関係").Range(g_range), PlotBy:=xlColumns

End Sub

Private Sub Graph_Refresh2()
'グラフの修正
    Dim g_range As String
    Sheets("雨量・流量結果図").Select
    g_range = "B2:D" & CStr(m + 2) & ",F2:G" & CStr(m + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheets("雨量・流量データ").Range(g_range), PlotBy:=xlColumns
End Sub

Private Sub DataWrite0(t As Integer)
'結果書き込み（デバック）
    Dim ws As Worksheet
    Dim i As Integer
    Set ws = Sheets("Ql")
    ws.Cells(2, t + 2) = t
    ws.Cells(3, t + 2) = tld
    ws.Cells(4, t + 2) = tl1
    ws.Cells(5, t + 2) = tl2
    For i = 1 To m2 - tl0
        'If I - Tl1 >= 0 Then
        '    ws.Cells(I + 5, T + 2) = Ql(I - Tl1)
        'End If
        ws.Cells(i + 5, t + 2) = ql(i)
    Next i
End Sub
</sxh>


====== tys-sfm_v_2_0_1 ======

tys-sfm_v_2_0_1.lzhをダウンロードして解凍すると、「<fc #FF0000>**tys-sfm_v_2_0_1.xls**</fc>」ができます。
この「<fc #FF0000>**tys-sfm_v_2_0_1.xls**</fc>」の使用方法を説明します。
ダウンロードについては、[[http://www.disaster-software.net|水理水文ソフトウェア]]のHPを参照してください。

===== マクロの説明 =====

本マクロは、各条件シートに条件を入力するだけで、<fs medium><fc #FF0000>**１流域についての総合貯留関数法と貯留関数法の２種類の流域流出解析計算**</fc></fs>が行えます。

  - 実績流量より、貯留関数係数の回帰を行う。
  - 任意の貯留関数係数を用いて流出解析を行う。
の２種類の計算を行うことができます。


==== 計算条件とデータの入力 ====

=== 計算条件 ===

「<fc #008080>計算条件</fc>」シートで貯留関数係数の回帰を行うかどうかをチェックしてください。


{{ :流域流出計算:image06.jpg?650|計算条件シート}}
;#;
計算条件シート
;#;

=== 雨量・流量データ ===

{{ :流域流出計算:image08.jpg?650|雨量・流量データシート}}
;#;
雨量・流量データシート
;#;
\\

  - **実績流量より、貯留関数係数の回帰を行う場合**　「<fc #008080>計算条件</fc>」シートの流域面積・基底流量・流路長さ・遅滞時間計算ピッチを入力してください。　「<fc #008080>雨量・流量データ</fc>」シートの流域平均雨量・実績流量・計算対象範囲指定のデータを入力してください。__計算対象範囲指定は入力した実績流量のうち、係数決定に用いるデータの範囲を指定します。開始点を「1」、終了点を「2」と入力してください。__
  - **任意の貯留関数係数を用いて流出解析を行う場合** 「<fc #008080>計算条件</fc>」シートの流域面積・基底流量・貯留関数係数（Ｋ・Ｐ値、遅滞時間、流出係数）を入力してください。「<fc #008080>雨量・流量データ</fc>」シートの流域平均雨量を入力してください。実績流量は、計算には使用しませんが、雨量・流量結果図で計算結果の評価を行うために入力する方が望ましいでしょう。
==== 計算結果 ====
  - 実績流量より、貯留関数係数の回帰を行う場合
    - 「<fc #008080>計算結果</fc>」シートに遅滞時間の設定値ごとの貯留量と流出量の関係の計算結果を書き込みます。この中で、貯留量と流出量のループの積算値の最小値を遅滞時間として採用します。尚、遅滞時間の繰り返し計算は、０時間～最大値（流路長を木村の式に適用した場合の遅滞時間＋１）を遅滞時間の計算ピッチ分です。ここで求めた遅滞時間より実績流量を用いた最小自乗法による回帰を行い、Ｋ・Ｐを算出します。
    - 「<fc #008080>雨量・流量データ</fc>」シートの計算流量の列に上記で求めた貯留関数係数を用いて計算した流量を書き込みます。計算結果は「<fc #008080>雨量・流量結果図</fc>」に反映されます。
    - 「<fc #008080>貯留量と流出量の関係</fc>」シートに決定した貯留関数係数を用いた貯留量と流出量の関係を出力します。その結果は「<fc #008080>貯留量と流出量の関係図</fc>」シートに反映されます。
  -任意の貯留関数係数を用いて流出解析を行う場合
    - 「<fc #008080>雨量・流量データ</fc>」シートの計算流量の列に入力した貯留関数係数を用いて計算した流量を書き込みます。計算結果は「<fc #008080>雨量・流量結果図</fc>」シートに反映されます。
    - 「<fc #008080>貯留量と流出量の関係</fc>」シートに貯留量と流出量の関係を出力します。その結果は「<fc #008080>貯留量と流出量の関係図</fc>」シートに反映されます。
\\

{{ :流域流出計算:image07.jpg?650|計算結果}}
;#;
計算結果
;#;

{{ :流域流出計算:image09.jpg?650|雨量・流量結果図}}
;#;
雨量・流量結果図
;#;

{{ :流域流出計算:image10.jpg?650|貯留量と流出量の関係図}}
;#;
貯留量と流出量の関係
;#;

{{ :流域流出計算:image11.jpg?650|貯留量と流出量の関係図}}
;#;
貯留量と流出量の関係図
;#;


==== 注意事項 ====

　本プログラムは、総合貯留関数法と通常の貯留関数法の両方が行えますが、総合貯留関数法と通常の貯留関数法では、Ｋ・Ｐ値はそれぞれ全く異なりますので注意してください。
　ここで、総合貯留関数法とは貯留関数を貯留高(mm)と流出高(mm/hr)の関係で表現したもの、通常の貯留関数法とは、貯留関数を貯留量(m3)と流出量(m3/s)の関係で表現したものです。総合貯留関数法は貯留関数法を単位流域面積で表現したものと言い換えることができます。
　本や論文によっては、総合貯留関数法を単に「貯留関数法」と表記されている場合があり、特殊な手法ではなく、かなり一般的に利用されていますので、基礎式による確認が必要です。
===== VBAソースファイル =====

　本マクロのソースを公開しますので、参考にしてください。

<sxh vb tys-sfm_v_2_0_1.bas>
Option Explicit

Public m As Integer     'ハイドロデータ数
Public m1 As Integer    'ハイドロ区間１
Public m2 As Long       'ハイドロ区間２
Public l As Double      '流域の流路延長(km)
Public a As Double      '流域面積(km2)
Public qi As Double     '基底流量(m3/s)
'
Public Const mnum As Integer = 1000 '配列数
Public r(0 To mnum) As Double       '雨量(mm/h)
Public q(0 To mnum) As Double       '実測流量(m3/s)
Public fr(0 To mnum) As Double      '雨量を流量に換算した値
Public fi(0 To mnum) As Double      '雨量を流量に換算した値の積算値
Public ql(0 To mnum) As Double      '遅滞時間を考慮した流量(基底流量を差し引いたもの)
Public qs(0 To mnum) As Double      '遅滞時間を考慮した流量の積分値(基底流量を差し引いたもの)
Public s(0 To mnum) As Double
Public sa(0 To 10000) As Double     'ループの積算値(時計回りか反時計回りかを判定)
'
Public rt As Double         '対象範囲の累積雨量(mm)
Public f As Double          '流出係数
Public tlk As Double        '遅滞時間(設定値)
Public qlm As Double        '流量の最大値(基底流量を差し引いたもの)
Public tpl As Double        '流量が最大の時間(遅滞時間分差し引いたもの)
Public qtm As Double        '流量の積算(遅滞時間分差し引いたもの)
Public qsm As Double        'Qsの積算
Public fim As Double        'Fiの最大値
Public sam As Double        'ループの最大
Public sas As Double        'ループの最小
Public tld As Double        '遅滞時間（決定分）
Public tl0 As Integer       '遅滞時間（切り上げ）
Public tl1 As Integer       '遅滞時間（整数部）
Public tl2 As Double        '遅滞時間小数部
Public t_step As Double     '遅滞時間の計算ステップ
Public num_cal As Integer    '遅滞時間計算回数
'
Public e As Double          '最小自乗法の対象となるデータの個数
Public sm As Double         '貯留高さの最大値
Public sar As Double        '上昇ループの積算値
Public sad As Double        '下降ループの積算値
Public k As Double          'Ｋ値
Public p As Double          'Ｐ値
'
Public cal_flag As Boolean  'Trueであれば総合貯留関数法

Public Sub Main()
'貯留関数法１流域のみ、係数を回帰で決定
    Dim ws1 As Worksheet
    Dim ws2 As Worksheet
    Dim ws3 As Worksheet
    Dim i As Integer
    Dim tl As Integer
    Dim t As Integer
    
    'シートのセット
    Set ws1 = Sheets("計算条件")
    Set ws2 = Sheets("雨量・流量データ")
    Set ws3 = Sheets("計算結果")
    
    'データ読み込み
    data_read1 ws1
    data_read2 ws2
    
    If Sheet1.CheckBox1.Value = True Then
        '係数回帰をするかどうかのチェック
        If Sheet1.OptionButton1.Value = True Then
            '総合貯留関数法
            cal_flag = True
        Else
            '貯留関数法
            cal_flag = False
        End If
    
        '遅滞時間の初期定義（木村の式）
        If l < 11.9 Then
            tl = 0
        Else
            tl = Int(0.5 + 0.047 * l - 0.56)   '切り上げ
        End If
    
        '適当な遅滞時間の判別
        '**********************
        num_cal = CInt(CDbl(tl + 1) / t_step)
        '**********************
        For t = 0 To num_cal
            tlk = t_step * CDbl(t)  '遅滞時間設定値
            tl1 = Int(tlk)          '遅滞時間整数部
            tl2 = tlk - tl1         '遅滞時間少数部
            If tl2 = 0 Then
                tl0 = tl1
            Else
                tl0 = tl1 + 1
            End If
    
            cal_ql
            cal_qtm
        
            'data_write0 t   'デバック
            
            If cal_flag = True Then
                '総合貯留関数法
                f = qtm / rt
            Else
                '貯留関数法
                f = qtm / ((a / 3.6) * rt)
            End If
        
            cal_fr
        
            cal_sa t
        
            data_write1 ws3, t   '計算経過出力
        '
        Next t
        
        'ループの最大・最小判定
        sam = 0
        For t = 0 To num_cal
            If Abs(sa(t)) >= sam Then sam = Abs(sa(t))
        Next t
        sas = sam
        For t = 0 To num_cal
            If Abs(sa(t)) <= sas Then
                sas = Abs(sa(t))
                tld = t_step * CDbl(t)  '遅滞時間設定値
            End If
        Next t
        tl1 = Int(tld)          '遅滞時間整数部
        tl2 = tld - tl1         '遅滞時間少数部
        If tl2 = 0 Then
            tl0 = tl1
        Else
            tl0 = tl1 + 1
        End If
    
        cal_ql
    
        cal_qtm
    
        If cal_flag = True Then
            '総合貯留関数法
            f = qtm / rt                '流出係数
        Else
            '貯留関数法
            f = qtm / ((a / 3.6) * rt)  '流出係数
        End If
    
        fim = 0
        For t = 1 To m
            If cal_flag = True Then
                '総合貯留関数法
                fr(t) = f * r(t)
            Else
                '貯留関数法
                fr(t) = a * f * r(t) / 3.6
            End If
            fim = fim + fr(t)
            fi(t) = fim
        Next t
    
        qsm = 0#
        For t = 2 To m - tl0
            qsm = qsm + (ql(t - 1) + ql(t)) / 2#
            qs(t) = qsm
            If fi(t) - qs(t) >= 0 Then s(t) = fi(t) - qs(t)
        Next t
    
        sm = 0
        For t = 2 To m - tl0 - 1
            If s(t) >= sm Then sm = s(t)
        Next t
    
        cal_kp  '最小自乗法によるK・Pの計算
    
        data_write2 ws1      '貯留関数式書き込み
        data_write3          '貯留高と流出高書き込み
        graph_refresh        'グラフ修正
    End If
    
    '流量計算
    
    data_read3 ws1
        
    cal_qt ws2          '求めた係数で流量を計算
        
    graph_refresh2       'グラフ修正
    
End Sub

Private Sub data_read1(ws As Worksheet)
'データ読み込み１（流域諸元）
    ws.Select
    a = ws.Cells(3, "D")        '流域面積(km2)
    qi = ws.Cells(4, "D")       '基底流量(m3/s)
    l = ws.Cells(5, "D")        '流路長(km)
    t_step = ws.Cells(6, "D")   '遅滞時間計算ピッチ(hr)
    If t_step = 0 Then t_step = 1
End Sub

Private Sub data_read2(ws As Worksheet)
'データ読み込み１（雨量・流量データ）
    Dim row_num As Integer
    Dim i As Integer
    
    '雨量・流量読み込み
    row_num = 3
    i = 0
    Do Until ws.Cells(row_num, "C") = ""    'データがあるまでループ
        i = i + 1
        r(i) = ws.Cells(row_num, "C")
        q(i) = ws.Cells(row_num, "D")
        row_num = row_num + 1
    Loop
    m = i       '雨量・流量データ数定義
    '計算対象範囲指定
    For i = 1 To m
        If ws.Cells(i + 2, "E") = 1 Then m1 = i         '開始点
        If ws.Cells(i + 2, "E") = 2 Then m2 = i         '終了点
    Next i
    '期間の累積雨量
    rt = 0
    For i = m1 To m2
        rt = rt + r(i)
    Next i
End Sub

Private Sub cal_fr()
'frの決定
    Dim t As Integer
    fim = 0
    For t = m1 To m2
        If cal_flag = True Then
            '総合貯留関数法
            fr(t) = f * r(t)
        Else
            '貯留関数法
            fr(t) = a * f * r(t) / 3.6
        End If
        '
        fim = fim + fr(t)
        fi(t) = fim
    Next t
End Sub

Private Sub cal_ql()
'遅滞時間分戻した流出高の定義（基底流量含まず）
    Dim t As Integer
    qlm = 0
    For t = 1 To m2
        If t - tl1 >= 0 Then
            ql(t - tl1) = q(t) + (q(t + 1) - q(t)) * tl2 - qi  '遅滞時間を考慮した流量(基底流量を差し引いたもの)
            
            '総合貯留関数法
            If cal_flag = True Then
                ql(t - tl1) = ql(t - tl1) / a * 3.6 '流量を流出高へ変換
            End If
            
            If ql(t - tl1) >= qlm Then
                qlm = ql(t - tl1)    '最大流量
                tpl = t - tlk        '最大流量の発生時間
            End If
        End If
    Next t
End Sub

Private Sub cal_qtm()
'遅滞時間を考慮した流量の積算（計算区間内）
    Dim t As Integer
    qtm = 0
    For t = m1 To m2 - tl0
        qtm = qtm + ql(t)
    Next t
End Sub

Private Sub cal_sa(tl As Integer)
'貯留高さと流出高をプロットした場合のループ形状
    Dim t As Integer
    qsm = 0: sar = 0: sad = 0
    
    For t = m1 To m2 - tl0 - 1
        qsm = qsm + (ql(t - 1) + ql(t)) / 2#
        qs(t) = qsm
        If fi(t) - qs(t) < 0 Then
            s(t) = 0#
        Else
            s(t) = fi(t) - qs(t)
        End If
    Next t
    '上昇ループの積算値計算
    For t = m1 To tpl - 1
        sar = sar + (s(t) + s(t + 1)) * Abs(ql(t + 1) - ql(t)) / 2#
    Next t
    '下降ループの積算値計算
    For t = tpl To m2 - tl0 - 1
        sad = sad - (s(t) + s(t + 1)) * Abs(ql(t) - ql(t + 1)) / 2#
    Next t
    'ループ全体の積算値計算
    sa(tl) = sar + sad
End Sub

Private Sub cal_kp()
'最小自乗法による貯留関数係数の計算
    Dim m3 As Integer
    Dim t As Integer
    Dim xi As Double
    Dim xx As Double
    Dim yi As Double
    Dim xy As Double
    Dim p0 As Double
    Dim k0 As Double
    Dim k1 As Double
    
    m3 = m2 - tl0
    e = 0
    For t = m1 To m3
        If s(t) <> 0 Then e = e + 1
    Next t
    xi = 0: xx = 0: yi = 0: xy = 0
    For t = m1 To m3
        xi = xi + Log(ql(t))
        xx = xx + Log(ql(t)) * Log(ql(t))
        yi = yi + Log(s(t))
        xy = xy + Log(ql(t)) * Log(s(t))
    Next t
    p0 = (xi * yi - e * xy) / (xi * xi - e * xx)
    k0 = (yi - p0 * xi) / e
    k1 = Exp(k0)
    k = k1
    p = p0
    k = Int(0.5 + k1 * 100#) / 100#     '桁合わせ
    p = Int(0.5 + p0 * 1000#) / 1000#   '桁合わせ
End Sub

Private Sub data_write1(ws As Worksheet, t As Integer)
'結果書き込み１（遅滞時間計算経過）
    Dim sa_state As String
    If t = 0 Then
        ws.Select
        Cells.Select
        Selection.ClearContents 'セルの値をクリア
        ws.Cells(2, "A") = "遅滞時間(hr)"
        ws.Cells(2, "B") = "最大流量(m3/s)"
        ws.Cells(2, "C") = "最大流量時間(hr)"
        ws.Cells(2, "D") = "流量積算値(m3/s)"
        ws.Cells(2, "E") = "流出係数"
        ws.Cells(2, "F") = "上昇ループの積算値"
        ws.Cells(2, "G") = "下降ループの積算値"
        ws.Cells(2, "H") = "ループ全体の積算値"
        ws.Cells(2, "I") = "ループの状態"
    End If
    '
    sa_state = ""
    If sa(t) > 0 Then sa_state = "時計回り"
    If sa(t) < 0 Then sa_state = "反時計回り"
    If sa(t) = 0 Then sa_state = "ぴったりOK"
    '
    ws.Cells(t + 3, "A") = tlk
    ws.Cells(t + 3, "B") = qlm
    ws.Cells(t + 3, "C") = tpl
    ws.Cells(t + 3, "D") = qtm
    ws.Cells(t + 3, "E") = f
    ws.Cells(t + 3, "F") = sar
    ws.Cells(t + 3, "G") = sad
    ws.Cells(t + 3, "H") = sa(t)
    ws.Cells(t + 3, "I") = sa_state
End Sub

Private Sub data_write2(ws As Worksheet)
'結果書き込み２（貯留関数式出力）
    ws.Select
    ws.Cells(10, "B") = k
    ws.Cells(11, "B") = p
    ws.Cells(12, "B") = tld
    ws.Cells(13, "B") = f
End Sub

Private Sub data_write3()
'貯留高と流出高書き込み
    Dim t As Integer
    Sheet8.Select
    If cal_flag = True Then
        '総合貯留関数法
        Sheet8.NAME = "貯留高と流出高の関係"
        Sheet8.Cells(1, "A") = "流出高と貯留高の関係"
        Sheet8.Cells(2, "A") = "流出高(mm/hr)"
        Sheet8.Cells(2, "B") = "貯留高(mm)"
    Else
        '貯留関数法
        Sheet8.NAME = "貯留量と流出量の関係"
        Sheet8.Cells(1, "A") = "流出量と貯留量の関係"
        Sheet8.Cells(2, "A") = "流出高(m3/s)"
        Sheet8.Cells(2, "B") = "貯留高(m3)"
    End If
    Sheet8.Cells(3, "D") = tld
    For t = 1 To m - tl0
        Sheet8.Cells(t + 2, "A") = ql(t)
        Sheet8.Cells(t + 2, "B") = s(t)
    Next t
End Sub

Private Sub data_read3(ws As Worksheet)
'データ読み込み3（貯留関数式元）
    ws.Select
    k = ws.Cells(10, "B")
    p = ws.Cells(11, "B")
    tld = ws.Cells(12, "B")
    f = ws.Cells(13, "B")
End Sub

Private Sub cal_qt(ws As Worksheet)
'求めた貯留関数を用いて流量を再度計算
    Dim t As Integer
    Dim tm As Integer
    Dim mr As Integer   '雨量・流量数
    Dim qql(0 To mnum) As Double
    Dim qqt(0 To mnum) As Double
    Dim qi2 As Double
    Dim j As Integer
    Dim c As Double
    Dim d As Double
    Dim x As Double
    Dim x2 As Double
    Dim y As Double
    
    mr = m
    tl1 = Int(tld)          '遅滞時間整数部
    tl2 = tld - tl1         '遅滞時間少数部
    
    If cal_flag = True Then
        '総合貯留関数法
        qi2 = qi / a * 3.6
    Else
        '貯留関数法
        qi2 = qi
    End If
    
     'ニュートン法による流量計算
    qql(1) = 0
    For t = 1 To mr
        If cal_flag = True Then
            '総合貯留関数法
            c = 2# * f * r(t) - qql(t) + 2# * (k * qql(t) ^ p)
        Else
            '貯留関数法
            c = 2# * f * a * r(t) / 3.6 - qql(t) + 2# * k * qql(t) ^ p
        End If
        '
        x = 0.0001
        j = j + 1
Step1:  y = x + (2# * k * x ^ p) - c
        d = 1# + 2# * p * k * x ^ (p - 1#)
        x2 = x - y / d
        If Abs(x2 - x) <= 0.0001 Then GoTo Step2
        x = x2
        GoTo Step1
Step2:  qql(t + 1) = x
        If t <= tl1 Then qqt(t) = qi2
        qqt(t + tl1) = qi2 + qql(t + 1) + (qql(t + 1) - qql(t)) * tl2
        If j <= 100 Then GoTo Step3
        If qqt(t) < 1.1 * qqt(99) Then GoTo Skip1
Step3:
    Next t
Skip1:
    For t = 1 To mr
        '総合貯留関数法
        If cal_flag = True Then
            qqt(t) = qqt(t) * a / 3.6
        End If
        ws.Cells(t + 2, "G") = qqt(t)
    Next t
End Sub

Private Sub graph_refresh()
'グラフの修正
    Dim g_range As String
    グラフ1.Select
    If cal_flag = True Then
        '総合貯留関数法
        グラフ1.NAME = "貯留高と流出高の関係図"
        ActiveChart.ChartTitle.Select
        Selection.Characters.Text = "貯留高と流出高の関係図"
        ActiveChart.Axes(xlCategory).AxisTitle.Select
        Selection.Characters.Text = "流出高(mm/hr)"
        ActiveChart.Axes(xlValue).AxisTitle.Select
        Selection.Characters.Text = "貯留高(mm)"
    Else
        '貯留関数法
        グラフ1.NAME = "貯留量と流出量の関係図"
        ActiveChart.ChartTitle.Select
        Selection.Characters.Text = "貯留量と流出量の関係図"
        ActiveChart.Axes(xlCategory).AxisTitle.Select
        Selection.Characters.Text = "流出量(m3/s)"
        ActiveChart.Axes(xlValue).AxisTitle.Select
        Selection.Characters.Text = "貯留量(m3)"
    End If
    g_range = "A3:B" & CStr(m - tl0 + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheet8.Range(g_range), PlotBy:=xlColumns

End Sub

Private Sub graph_refresh2()
'グラフの修正
    Dim g_range As String
    Sheets("雨量・流量図").Select
    g_range = "B2:D" & CStr(m + 2) & ",F2:G" & CStr(m + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheets("雨量・流量データ").Range(g_range), PlotBy:=xlColumns
End Sub

Private Sub data_write0(t As Integer)
'結果書き込み（デバック）
    Dim ws As Worksheet
    Dim i As Integer
    Set ws = Sheets("デバッグQl")
    ws.Cells(2, t + 2) = t
    ws.Cells(3, t + 2) = tld
    ws.Cells(4, t + 2) = tl1
    ws.Cells(5, t + 2) = tl2
    For i = 1 To m2 - tl0
        'If i - tl1 >= 0 Then
        '    ws.Cells(i + 5, t + 2) = ql(i - tl1)
        'End If
        ws.Cells(i + 5, t + 2) = ql(i)
    Next i
End Sub

</sxh>


====== tys-sfb_v1_0_0 ======

tys-sfb_v1_0_0.lzhをダウンロードして解凍すると、「<fc #FF0000>**tys-sfb_v1_0_0.xls**</fc>」ができます。
この「<fc #FF0000>**tys-sfb_v1_0_0.xls**</fc>」の使用方法を説明します。
ダウンロードについては、[[http://www.disaster-software.net|水理水文ソフトウェア]]のHPを参照してください。


===== マクロの説明 =====

本マクロは、各条件シートに条件を入力するだけで、<fs medium><fc #FF0000>**複数流域の総合貯留関数による流出解析計算**</fc></fs>が行えます。

  - 複数流域を一括して流出解析を行う。
==== 計算条件とデータの入力 ====

=== 計算条件 ===

「<fc #008080>計算条件</fc>」シートで貯留関数係数の回帰を行うかどうかをチェックしてください。


{{ :流域流出計算:image31.jpg?650|計算条件シート}}
;#;
計算条件シート
;#;

=== 河道諸元 ===

{{ :流域流出計算:image32.jpg?650|河道諸元シート}}
;#;
河道諸元シート
;#;

=== 流域諸元 ===

{{ :流域流出計算:image33.jpg?650|流域諸元シート}}
;#;
流域諸元シート
;#;

=== 雨量データ ===

{{ :流域流出計算:image34.jpg?650|雨量データシート}}
;#;
雨量データシート
;#;
\\

==== 計算結果 ====

{{ :流域流出計算:image35.jpg?650|流域計算結果}}
;#;
流域計算結果
;#;

{{ :流域流出計算:image36.jpg?650|河道計算結果}}
;#;
河道計算結果図
;#;

{{ :流域流出計算:image37.jpg?650|作図データ}}
;#;
作図データ
;#;

{{ :流域流出計算:image38.jpg?650|計算結果図}}
;#;
計算結果図
;#;


==== 注意事項 ====

　本プログラムは、総合貯留関数法と通常の貯留関数法の両方が行えますが、総合貯留関数法と通常の貯留関数法では、Ｋ・Ｐ値はそれぞれ全く異なりますので注意してください。
　ここで、総合貯留関数法とは貯留関数を貯留高(mm)と流出高(mm/hr)の関係で表現したもの、通常の貯留関数法とは、貯留関数を貯留量(m3)と流出量(m3/s)の関係で表現したものです。総合貯留関数法は貯留関数法を単位流域面積で表現したものと言い換えることができます。
　本や論文によっては、総合貯留関数法を単に「貯留関数法」と表記されている場合があり、特殊な手法ではなく、かなり一般的に利用されていますので、基礎式による確認が必要です。
===== VBAソースファイル =====

　本マクロのソースを公開しますので、参考にしてください。

<sxh vb tys-sfb_v_1_0_0.bas>
Option Explicit
Option Base 1

'貯留関数による洪水流量計算
Sub main()
    Dim ws1 As Worksheet, ws2 As Worksheet, ws3 As Worksheet, ws4 As Worksheet
    Dim ws5 As Worksheet, ws6 As Worksheet, ws7 As Worksheet
    Set ws1 = Sheets("計算条件")
    Set ws2 = Sheets("流域諸元")
    Set ws3 = Sheets("河道諸元")
    Set ws4 = Sheets("雨量データ")
    Set ws5 = Sheets("流域計算結果")
    Set ws6 = Sheets("河道計算結果")
    Set ws7 = Sheets("作図データ")
    
    Application.StatusBar = "データ読み込み中"
    data_read ws1, ws2, ws3, ws4
    inital_set ws5, ws6
    
    For m = 1 To m1     '流域数のループ
        Application.StatusBar = "流域計算中" & Str(m) & "/" & m1
        cal_basin m, ws5
    Next m
    For n = 1 To n1     '河道数のループ
        Application.StatusBar = "河道計算中" & Str(n) & "/" & n1
        cal_river n, ws6, ws7
    Next n
    
    Graph_Refresh2
    
    Application.StatusBar = False
End Sub

Sub inital_set(ws5 As Worksheet, ws6 As Worksheet)
    Dim i As Integer, j As Integer, l As Integer
    For i = 1 To 50
        qb(i) = 0#
        qi(i) = 0#
        For n = 1 To 500
            q8sum(i, n) = 0#
        Next n
    Next i
    
    '   降雨量の計算   ---------------------------------------------------
    dn = d3 / d2
    For i = 1 To nt
        For l = 1 To dn
            i2 = (i - 1) * dn + l
            For j = 1 To nhyeto
                r1(j, i2) = r(j, i) / d3 * d2
            Next j
        Next l
    Next i
    i1 = Int(nt * dn)
    '結果書き込み用のシートのクリア
    ws5.Select
    Cells.Select
    Selection.ClearContents 'セルの値をクリア
    ws6.Select
    Cells.Select
    Selection.ClearContents 'セルの値をクリア
End Sub

Sub cal_basin(m As Integer, ws5 As Worksheet)
    Dim i As Integer, l As Integer
        '流域の計算
        '初期値の消去
        For i = 1 To 500
            r2(i) = 0#
            r3(i) = 0#
            qr(i) = 0#
            qs(i) = 0#
            q1(i) = 0#
            q2(i) = 0#
            q3(i) = 0#
            q4(i) = 0#
            q5(i) = 0#
            q6(i) = 0#
            q7(i) = 0#
            q8(i) = 0#
        Next i
        
        '流域降雨     , 流出域降雨の計算
        r3(1) = r1(norain(m), 1)
        r2(1) = r1(norain(m), 1)
        
        For i = 2 To i1
            r3(i) = r3(i - 1) + r1(norain(m), i)
            r2(i) = r1(norain(m), i)
        Next i
        
        '浸透域有効降雨の計算
        rs(1) = 0#
        For i = 2 To i1
            rs(i) = 0#
            If r3(i) > sr(m) Then GoTo S1
        Next i
        
S1:     For i = i To i1
            rs(i) = 0#
            If r3(i) > sr(m) Then rs(i) = r3(i) - r3(i - 1)
            If r3(i - 1) < sr(m) Then rs(i) = r3(i) - sr(m)
        Next i
        
        '流出域　流出高さの計算
        i = 0
S2:     i = i + 1
        d5 = d1
S3:     i2 = Int(d2 / d5) + 1
        qt = 0#
        If i > 1 Then qt = qr(i - 1)
        rt = r2(i) / d2
        tq(1) = qt
        t1(i2) = rt
        t1(1) = 0#
        If i > 1 Then t1(1) = r2(i) / d2
        dt = d5
        qk = qt
        t1(1) = t1(i2)
        l = 1
S4:     l = l + 1
        t1(l) = (t1(i2) - t1(1)) * CDbl(l - 1) / CDbl(i2 - 1) + t1(1)
        qt = tq(l - 1)
        qk = qt
        t1(1) = t1(i2)
        rt = t1(l)
        k_p_cal 0#, kr(m), pr(m)
        If tq(1) = 0# Then GoTo S5
        If dt > 2# * kr(m) * pr(m) / tq(1) ^ (1# - pr(m)) Then GoTo S6
        If dt > 2# * kr(m) / tq(1) ^ (1# - pr(m)) Then GoTo S6
S5:     If qk = 0# Then GoTo S7
        If dt < 2# * kr(m) * pr(m) / qk ^ (1# - pr(m)) Then GoTo S7
        If dt < 2# * kr(m) / qk ^ (1# - pr(m)) Then GoTo S7
S6:     d5 = d5 / 2#
        GoTo S3
S7:     tq(l) = qk
        If l < i2 Then GoTo S4
        qr(i) = tq(i2)
        If i < i1 Then GoTo S2
        
        '浸透域流出ﾃﾞｰﾀの計算
        i = 0
S8:     i = i + 1
        d5 = d1
S9:     i2 = Int(d2 / d5) + 1
        qt = 0#
        If i > 1 Then qt = qs(i - 1)
        rt = rs(i) / d2
        dt = d5
        tq(1) = qt
        t1(i2) = rt
        t1(1) = 0#
        If i > 1 Then t1(1) = rs(i - 1) / d2
        qk = qt
        t1(1) = t1(i2)
        l = 1
S10:    l = l + 1
        t1(l) = (t1(i2) - t1(1)) * CDbl(l - 1) / CDbl(i2 - 1) + t1(1)
        qt = tq(l - 1)
        qk = qt
        rt = t1(l)
        t1(1) = t1(i2)
        k_p_cal 0#, kr(m), pr(m)
        If tq(1) = 0 Then GoTo S11
        If dt > 2# * kr(m) * pr(m) / tq(1) ^ (1# - pr(m)) Then GoTo S12
        If dt > 2# * kr(m) / tq(1) ^ (1# - pr(m)) Then GoTo S12
S11:    If qk = 0# Then GoTo S13
        If dt < 2# * kr(m) * pr(m) / qk ^ (1# - pr(m)) Then GoTo S13
        If dt < 2# * kr(m) / qk ^ (1# - pr(m)) Then GoTo S13
S12:    d5 = d5 / 2#
        GoTo S9
S13:    tq(l) = qk
        If l < i2 Then GoTo S10
        qs(i) = tq(i2)
        If i < i1 Then GoTo S8
        
        '遅滞時間降雨量
        tl = tr(m) / d2
        For i = 1 To i1
            z1 = Int(tl)
            z2 = tl - z1
            q1(i + z1) = qr(i)
            q2(i + z1) = qs(i)
            q3(i + z1) = q1(i + z1)
            q4(i + z1) = q2(i + z1)
            If (i + z1 - 1) > 0 Then q3(i + z1) = q1(i + z1 - 1) + (q1(i + z1) - q1(i + z1 - 1)) * (1# - z2)
            If (i + z1 - 1) > 0 Then q4(i + z1) = q2(i + z1 - 1) + (q2(i + z1) - q2(i + z1 - 1)) * (1# - z2)
        Next i
        
        '流域　流出高ﾃﾞｰﾀ
        For i = 1 To i1
            q5(i) = q3(i) * f(m) + q4(i) * (1# - f(m)) + 3.6 * bq(m) / aa(m)
            q6(i) = q5(i) * aa(m) / 3.6
        Next i
        
        '流入する河道への流量の受け渡し
        For i = 1 To i1
            q8sum(er(m), i) = q8sum(er(m), i) + q6(i)
        Next i
        
        qb(er(m)) = qb(er(m)) + bq(m)
        
        print_basin m             '計算結果出力（ファイル）
        
        print_basin2 m, ws5       '計算結果出力（ワークシート）
End Sub

Sub cal_river(n As Integer, ws6 As Worksheet, ws7 As Worksheet)
    Dim i As Integer, l As Integer
        '河道の計算
        For i = 1 To i1
            q7(i) = q8sum(n, i)
        Next i
        
         '河道の上流に連結する河道の基底流量を足し合わせてその河道の基底流量を定義
        For i = 1 To n1
            If ek(i) = n Then qi(n) = qi(n) + qi(i)
        Next i
        qi(n) = qi(n) + qb(n)
        q7(1) = qi(n)
        q8(1) = qi(n)
        q2(1) = qi(n)
        For i = 2 To i1
            q7(i) = q7(i)
            q2(i) = q7(i - 1) + (q7(i) - q7(i - 1)) * (d2 - tk(n)) / d2
        Next i
        
        '繰り返し時間数のループ
        For i = 2 To i1
            d5 = d1
S1:         i2 = Int(d2 / d5) + 1
            tq(1) = q8(1)
            If i > 1 Then tq(1) = q8(i - 1)
            t1(i2) = q2(i)
            t1(1) = qi(n)
            If i > 1 Then t1(1) = q2(i - 1)
            dt = d5
            l = 1
            l = l + 1
            For l = 2 To i2
                t1(l) = (t1(i2) - t1(1)) * CDbl(l - 1) / CDbl(i2 - 1) + t1(1)
                qt = tq(l - 1)
                rt = (t1(l) + t1(l - 1)) / 2#
                qk = qt
                k_p_cal qi(n), kk(n), pk(n)
                If tq(1) = 0# Then GoTo S2
                If dt > 2# * kk(n) * pk(n) / tq(1) ^ (1# - pk(n)) Then GoTo S3
                If dt > 2# * kk(n) / tq(1) ^ (1# - pk(n)) Then GoTo S3
S2:             If qk = 0# Then GoTo S4
                If dt < 2# * kk(n) * pk(n) / qk ^ (1# - pk(n)) Then GoTo S4
                If dt < 2# * kk(n) / qk ^ (1# - pk(n)) Then GoTo S4
S3:             d5 = d5 / 2#
                GoTo S1
S4:             tq(l) = qk
            Next l
            q8(i) = tq(i2)
            If q8(i) < qi(n) Then q8(i) = qi(n)
        Next i
        For i = 1 To i1
            q8sum(ek(n), i) = q8sum(ek(n), i) + q8(i)
        Next i
        
        '計算結果出力（ファイル）
        print_river n
        
        '計算結果出力（ワークシート）
        print_river2 n, ws6
        
        If n = numgraph Then print_plot ws7  '作図データの書き込み

End Sub

'*************************  K-Pの計算　********************************************
Sub k_p_cal(qi1 As Double, k As Double, p As Double)
    Dim rq1 As Double, qb1 As Double
    Dim i As Integer
    
    qp1 = 0#
    qq1 = 0#
    qm1 = 0#
    qn1 = 0#
    qk1 = qk
    qb1 = 0#
    rq1 = 0#
    
    For i = 1 To 750
        qe1 = 2# * rt - qt - 2# * (k * (qk1 ^ p) - k * (qt ^ p)) / dt
        qg1 = qk1 - qe1
        If Abs(qg1) < 0.0001 Then GoTo S1
        rq1 = qk1
        If qg1 > 0# Then qp1 = qg1
        If qg1 > 0# Then qq1 = rq1
        If qg1 < 0# Then qm1 = qg1
        If qg1 < 0# Then qn1 = rq1
        qa = qk1
        qk1 = qe1
        If i = 1 Then GoTo S2
        qk1 = qq1 - qp1 * (qq1 - qn1) / (qp1 - qm1)
        qb1 = qk1
        If qa = qb1 Then GoTo S1
S2:     If qk1 < 0# Then qk1 = 0#
    Next i
S1: qk = qk1
    If qk < qi1 Then qk = qi1
End Sub

'*********************** データの読み込み *******************************
Sub data_read(ws1 As Worksheet, ws2 As Worksheet, ws3 As Worksheet, ws4 As Worksheet)
    Dim i As Integer, j As Integer
    '計算条件の読み込み
        d1 = ws1.Cells(2, 3)
        d2 = ws1.Cells(3, 3)
        d3 = ws1.Cells(4, 3)
        nt = ws1.Cells(5, 3)
        m1 = ws1.Cells(6, 3)
        n1 = ws1.Cells(7, 3)
        nhyeto = ws1.Cells(8, 3)
    
    'グラフ描画河道の読み込み
        numgraph = ws1.Cells(9, 3)
        
    '流域諸元の読み込み
        For i = 1 To m1
            kr(i) = ws2.Cells(i + 2, 2)
            pr(i) = ws2.Cells(i + 2, 3)
            tr(i) = ws2.Cells(i + 2, 4)
            f(i) = ws2.Cells(i + 2, 5)
            sr(i) = ws2.Cells(i + 2, 6)
            aa(i) = ws2.Cells(i + 2, 7)
            bq(i) = ws2.Cells(i + 2, 8)
            er(i) = ws2.Cells(i + 2, 9)
            norain(i) = ws2.Cells(i + 2, 10)
        Next i
    '河道諸元の読み込み
        For i = 1 To n1
            kk(i) = ws3.Cells(i + 2, 2)
            pk(i) = ws3.Cells(i + 2, 3)
            tk(i) = ws3.Cells(i + 2, 4)
            ek(i) = ws3.Cells(i + 2, 5)
        Next i
    '雨量データの読み込み
        For i = 1 To nt
            For j = 1 To nhyeto
                r(j, i) = ws4.Cells(i + 2, j + 1)
            Next j
        Next i

End Sub

'******************** 計算結果のファイルへの出力 ******************************
Sub print_basin(m As Integer)
    Dim NAME As String, OUFL1 As String
    Dim i As Integer
    Dim D2z As String, R1z As String, R3z As String, RSz As String, QRz As String, QSz As String
    Dim Q5z As String, Q6z As String
    NAME = ActiveWorkbook.Path
    OUFL1 = "RUN-R.OUT"
    '   流域計算出力　　　　------------------------------------------------
    If m = 1 Then Open NAME + "\" + OUFL1 For Output As #2
        Print #2, "流域 ("; m; ")"
        Print #2, "流域の諸元    "
        Print #2, "K 値               = "; Format(Format(kr(m), "##0.00"), "@@@@@@")
        Print #2, "P 値               = "; Format(Format(pr(m), "##0.00"), "@@@@@@")
        Print #2, "遅滞時間(hr)    TR = "; Format(Format(tr(m), "##0.00"), "@@@@@@")
        Print #2, "流出率          F  = "; Format(Format(f(m), "##0.00"), "@@@@@@")
        Print #2, "飽和雨量(mm/hr) RS = "; Format(Format(sr(m), "##0.00"), "@@@@@@")
        Print #2, "流域面積(km2)   A  = "; Format(Format(aa(m), "##0.00"), "@@@@@@")
        Print #2, "基底流量(m3/s)  BQ = "; Format(Format(bq(m), "##0.00"), "@@@@@@")
        Print #2, "連立する河道 NO    = "; Format(er(m), "@@@@@")
        Print #2, "     時間   降雨量  累加雨量  浸透 R   流出 Q   浸透 Q   総合 Q      流出量 "
        Print #2, "     (時)    (㎜)     (㎜)     (㎜)     (㎜)    (㎜)     (㎜)         (m3/s)"
        For i = 1 To i1
            D2z = Format(Format(i * d2, "#####0.00"), "@@@@@@@@@")
            R1z = Format(Format(r1(norain(m), i), "#####0.00"), "@@@@@@@@@")
            R3z = Format(Format(r3(i), "#####0.00"), "@@@@@@@@@")
            RSz = Format(Format(rs(i), "#####0.00"), "@@@@@@@@@")
            QRz = Format(Format(qr(i), "#####0.00"), "@@@@@@@@@")
            QSz = Format(Format(qs(i), "#####0.00"), "@@@@@@@@@")
            Q5z = Format(Format(q5(i), "#####0.00"), "@@@@@@@@@")
            Q6z = Format(Format(q6(i), "######0.00"), "@@@@@@@@@@")
            Print #2, D2z; R1z; R3z; RSz; QRz; QSz; Q5z; Q6z
        Next i
     If m = m1 Then Close #2
End Sub

'******************** 計算結果のファイルへの出力 ******************************
Sub print_river(n As Integer)
    Dim NAME As String, OUFL2 As String
    Dim i As Integer
    Dim D2z As String, Q7z As String, Q2z As String, Q8z As String
    NAME = ActiveWorkbook.Path
    OUFL2 = "RUN-K.OUT"
    If n = 1 Then Open NAME + "\" + OUFL2 For Output As #3
    '   河道計算結果の出力 -----------------------------------------------
        Print #3, "河道 ("; n; ")"
        Print #3, "河道の諸元 "
        Print #3, "K 値           = "; Format(Format(kk(n), "##0.00"), "@@@@@@")
        Print #3, "P 値           = "; Format(Format(pk(n), "##0.00"), "@@@@@@")
        Print #3, "遅滞時間(hr)TK = "; Format(Format(tk(n), "##0.00"), "@@@@@@")
        Print #3, "     時間     流入量    流入量    流出量   "
        Print #3, "                     (遅滞時間)"
        Print #3, "     (時)     (m3/s)   (m3/s)    (m3/s)"
        For i = 1 To i1
            D2z = Format(Format(i * d2, "#####0.00"), "@@@@@@@@@")
            Q7z = Format(Format(q7(i), "######0.00"), "@@@@@@@@@@")
            Q2z = Format(Format(q2(i), "######0.00"), "@@@@@@@@@@")
            Q8z = Format(Format(q8(i), "######0.00"), "@@@@@@@@@@")
            Print #3, D2z; Q7z; Q2z; Q8z
        Next i
    If n = n1 Then Close #3
End Sub

'******************** 計算結果のワークシートへの出力 ******************************
Sub print_basin2(m As Integer, ws5 As Worksheet)
    Dim i As Integer, j As Integer
    ws5.Select
    '   流域計算出力　　　　------------------------------------------------
    If m = 1 Then
        ws5.Cells(2, 1) = "時間(hr)"
        For i = 1 To i1
            ws5.Cells(i + 2, 1) = i * d2
        Next i
        For j = 1 To nhyeto
            ws5.Cells(2, (j - 1) * 2 + 2) = "降雨量(" & j & ") (mm)"
            ws5.Cells(2, (j - 1) * 2 + 3) = "累積降雨量(" & j & ") (mm)"
            For i = 1 To i1
                ws5.Cells(i + 2, (j - 1) * 2 + 2) = r1(norain(m), i)
                ws5.Cells(i + 2, (j - 1) * 2 + 3) = r3(i)
            Next i
        Next j
    End If
    ws5.Cells(1, nhyeto * 2 + 2 + (m - 1)) = "流域 (" + CStr(m) + ")"
    ws5.Cells(2, nhyeto * 2 + 2 + (m - 1)) = "流域流出量(m3/s)"
    For i = 1 To i1
        ws5.Cells(i + 2, nhyeto * 2 + 2 + (m - 1)) = q6(i)
    Next i
    If m = m1 Then
        For i = 1 To (nhyeto * 2 + 1)
            Columns(i).AutoFit                        '列幅を整える
            If i = 1 Then
                Columns(i).NumberFormatLocal = "0.0_ "              '桁数を揃える
            Else
                Columns(i).NumberFormatLocal = "0.00_ "   '桁数を揃える
            End If
        Next i
        For i = (nhyeto * 2 + 1) To m1 + (nhyeto * 2 + 1)
            Columns(i).NumberFormatLocal = "0.000_ "        '桁数を揃える
            Columns(i).AutoFit                              '列幅を整える
        Next i
    End If
End Sub

'******************** 計算結果のワークシートへの出力 ******************************
Sub print_river2(n As Integer, ws6 As Worksheet)
    Dim i As Integer
    ws6.Select
    If n = 1 Then
        ws6.Cells(2, 1) = "時間(hr)"
        For i = 1 To i1
            ws6.Cells(i + 2, 1) = i * d2
        Next i
    End If
    ws6.Cells(1, (n - 1) * 2 + 2) = "河道 (" + CStr(n) + ")"
    ws6.Cells(2, (n - 1) * 2 + 2) = "河道流入量(m3/s)"
    ws6.Cells(2, (n - 1) * 2 + 3) = "河道流出量(m3/s)"
    For i = 1 To i1
        ws6.Cells(i + 2, (n - 1) * 2 + 2) = q7(i)
        ws6.Cells(i + 2, (n - 1) * 2 + 3) = q8(i)
    Next i
    If n = n1 Then
        Columns(1).NumberFormatLocal = "0.0_ "          '桁数を揃える
        Columns(1).AutoFit                              '列幅を整える
        For i = 2 To n1 * 2 + 1
            Columns(i).NumberFormatLocal = "0.000_ "    '桁数を揃える
            Columns(i).AutoFit                          '列幅を整える
        Next i
    End If
End Sub


'******************** 計算結果のワークシートへの出力（グラフ用） *************************
Sub print_plot(ws7 As Worksheet)
    Dim i As Integer
    
    ws7.Cells(1, 3) = "NO." & numgraph
    ws7.Select
    For i = 1 To i1
        ws7.Cells(i + 1, 4) = q8(i)
    Next i
End Sub


'******************** 計算結果図の修正 *************************
Private Sub Graph_Refresh2()
'グラフの修正
    Dim g_range As String
    Sheets("計算結果図").Select
    
    
    ActiveChart.SeriesCollection(2).Select
    Selection.Delete
    ActiveChart.SeriesCollection(1).Select
    Selection.Delete
    
    g_range = "B1:C" & CStr(nt + 1) & ",D1:D" & CStr(nt + 1)
'
    Sheets("作図データ").Select
    Range(g_range).Select
    Selection.Copy
    Sheets("計算結果図").Select
    ActiveChart.SeriesCollection.Paste Rowcol:=xlColumns, SeriesLabels:=True, _
        CategoryLabels:=True, Replace:=True
'    ActiveChart.SeriesCollection.Paste Rowcol:=xlColumns, SeriesLabels:=True, _
        CategoryLabels:=True, Replace:=False, NewSeries:=True
        
    With ActiveChart
        .HasTitle = True
        .ChartTitle.Characters.Text = "計算結果 " & "河道" & "NO." & numgraph
    End With

End Sub

</sxh>