tys-kf1d_v_1_0_1

tys-kf1d_v_1_0_1.lzhをダウンロードして解凍すると、「tys-kf1d_v_1_0_1.xls」ができます。 この「tys-kf1d_v_1_0_1.xls」の使用方法を説明します。 ダウンロードについては、水理水文ソフトウェアのHPを参照してください。

本マクロは、各条件シートに条件を入力するだけで、Kinematic Wave法による1次元不定流計算が行えます。

tys-kf1d_v_1_0_1モデルの特徴は

1. 支川の合流も可能（流域内の河道網を表現して、一括で流出解析ができる）。
2. 任意の河道断面に斜面から降雨による流入量を設定できる。

「計算条件」シートに計算パラメータを入力します。

• 基本パラメータ
  • 総ハイエト数
  • 河道数
  • 継続時間区分数
  • ハイドロ数
  • 流量判定値
• 出力設定
  • 出力箇所数
  • 出力ピッチ(s)
• 出力箇所（50カ所まで）
  • 河道番号
  • 断面番号
  • 名称
• 継続時間区分に応じた時間緒元
  • 継続時間(s)
  • 計算刻み時間(s)

　計算結果は、以下の８つのシートが作成され、結果が書き出されます。　グラフについては、使用者各自に作成して下さい。

• 計算フルード数シート
• 計算流速シート
• 計算水深シート
• 計算流量シート
• モニター地点水深シート
• モニター地点水深図シート
• モニター地点流量シート
• モニター地点流量図シート

Kinematic Wave法による流出解析の計算部分のソースコードを公開します。

Option Explicit

Public Sub kinemaMain()
    'キネマによる流出解析
    Dim ir As Integer, it As Long
    Dim ii As Integer, kr As Integer, irr As Long
    Dim dtl As Double
    Dim flag As Boolean
    Dim qcount As Long
    Dim hcount As Long
    Dim vcount As Long
    Dim fcount As Long
      
    '計算データの読み込み
    Call DataInput(flag)
    If flag = False Then End

    '初期設定
    Call setInitial(flag)
    If flag = False Then End

    irr = 0
    dtl = 0#    '経過時間(s)
    qcount = 0
    hcount = 0
    vcount = 0
    fcount = 0
    
    For ir = 1 To nndata

        ii = khy1(ir)
        dtt = dt(ii)
        idt = CLng(delt(ii) / dtt)
        For it = 1 To idt

            dtl = dtl + dtt
            '洪水解析
            For kr = km To 1 Step -1
                'Call setSlope(ir, it, kr) '斜面からの流入
                Call calRiver(ir, it, kr) '河道洪水解析
            Next kr
            Call resetData      '変数の更新
            
        Next it
    Next ir
    
End Sub

Private Sub setInitial(flag As Boolean)
    '計算開始時の初期化
    Dim i As Integer, j As Integer, k As Integer, l As Integer, m As Integer
    Dim kk As Integer
    
    On Error GoTo DataError
    
    '各河道断面の基底流量の設定
    For k = km To 1 Step -1
        bqr2(k, ns(k)) = bqr(k)
        For j = ns(k) - 1 To 1 Step -1
          bqr2(k, j) = bqr2(k, j + 1)
          If kgs(k, j) > 0 Then
                For kk = 1 To km
                    If kg(kk) = k And kg1(kk) = j Then bqr2(k, j) = bqr2(k, j) + bqr2(kk, 1)
                Next kk
          End If
        Next j
    Next k

    '変数の初期化
    For k = 1 To km
        For l = 1 To ks(k)
            qrs(k, l) = 0#
        Next l
        
        For j = 1 To ns(k)
            ar1(k, j) = 0#
            qqq0(k, j) = 0#
            qqq1(k, j) = 0#
            qmax(k, j) = 0#
            dep(k, j) = 0#
            vel(k, j) = 0#
            frd(k, j) = 0#
        Next j
    Next k
    
    flag = True
    Exit Sub
    
DataError:
    MsgBox "初期化ルーチンでエラーが発生しました。" & vbNewLine _
                    & "処理を中断します。", vbCritical
    flag = False
End Sub

Private Sub resetData()
    '変数の入れ替え
    Dim j As Integer, k As Integer
    
    For k = 1 To km
        For j = 1 To ns(k)
            qqq0(k, j) = qqq1(k, j)
            If qmax(k, j) < qqq1(k, j) Then qmax(k, j) = qqq1(k, j)
        Next j
    Next k
End Sub

Private Sub calRiver(ir As Integer, it As Long, k As Integer)
    '河道の計算
    Dim j As Integer, kr2 As Integer
    Dim dx As Double
    Dim aa As Double, bb As Double, cc As Double, ee As Double
    
    '流量と流積の関係を求める - ---計算の開始時のみ
    If ir = 1 And it = 1 Then Call secoef(k)

    For j = ns(k) To 1 Step -1

        If j = 1 Then
            sita(j) = Atn((z(k, j + 1) - z(k, j)) / rl(k, j + 1))
        Else
            sita(j) = Atn((z(k, j) - z(k, j - 1)) / rl(k, j))
        End If

        '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        'エネルギー勾配で計算する方法も検討が必要
        '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++       
        If sita(j) < 0# Then sita(j) = 0.00000001    '逆勾配になっても計算が止まらないように
        
        ar(j) = ar1(k, j)
        If j = 1 Then
            dx = rl(k, j + 1)
        ElseIf j = ns(k) Then
            dx = rl(k, j)
        Else
            dx = rl(k, j + 1)
        End If

        '支川横流入流量の計算
        Call qPlus(k, j)

        ' ------新しい時間の流積
        If j <> ns(k) Then
            ar1(k, j) = ar(j) - (qqq0(k, j) - qqq0(k, j + 1) - bqr2(k, j) - qpls(k, j)) / dx * dtt
        Else
            ar1(k, j) = ar(j) - (qqq0(k, j) - bqr2(k, j) - qpls(k, j)) / dx * dtt
        End If
        If ar1(k, j) < 0# Then
            MsgBox "河道番号" & CStr(k) & "、断面番号" & CStr(j) & "の流積が負になりました。強制終了します。", vbCritical
            'Stop
            End
        End If
        
        '--------------------------------
        qqq1(k, j) = calQ(k, j, ar1(k, j))
        '--------------------------------
        'b(k, j) = b0(k, j)
        
        If qqq1(k, j) < qeps Then qqq1(k, j) = bqr(k)     'qeps以下の流量の場合は基底流量とする。
        
    
        If ar1(k, j) <= 0# Then
            dep(k, j) = 0#
            vel(k, j) = 0#
            frd(k, j) = 0#
        Else
            '-----台形形状の流積から水深を逆算
            aa = b1(k, j) + b2(k, j)
            bb = 2# * b(k, j)
            cc = -2# * ar1(k, j)
            ee = bb ^ 2# - 4# * aa * cc
            dep(k, j) = (-bb + Sqr(ee)) / 2# / aa
            If qqq1(k, j) > 0 Then
                'vel(k, j) = qqq1(k, j) / dep(k, j)
                vel(k, j) = 1# / rn(k, j) * dep(k, j) ^ (2# / 3#) * sita(j) ^ 0.5
                frd(k, j) = vel(k, j) / Sqr(g * dep(k, j))
            End If
        End If
    Next j
End Sub
'
Private Function calQ(k As Integer, j As Integer, aa As Double)
    '運動方程式より流量Qを計算する
    Dim rr As Double, ss As Double, ee As Double
    Dim ccc As Double, pp As Double
    
    rr = ck1(k, j) * (aa ^ (cz(k, j))) '径深
    ss = sita(j)
    ee = Sin(ss)


    'キネマの式
    pp = 1# / 2#
    ccc = 1# / rn(k, j) * rr ^ (1# / 6#)
    calQ = aa * ccc * Sqr(ee) * rr ^ (pp)
    'If CalQ<0.0 then CalQ=0.0
    If calQ < 0# Then
        MsgBox "流量が負になりました。強制終了します。" & vbNewLine & _
                " 河道番号：" & CStr(k) & "、断面番号" & CStr(j), vbCritical
        Stop
        End
    End If
End Function

Private Sub qPlus(k As Integer, j As Integer)
    '特定河道断面への支川からの横流入量を計算
    Dim l As Integer, kr As Integer
    Dim DosyaA As Double, DosyaRate As Double
    
    qpls(k, j) = 0#
    '流入支川の量
    If kgs(k, j) > 0 Then
        For kr = 1 To km
            If kg(kr) = k And kg1(kr) = j Then
                qpls(k, j) = qpls(k, j) + qqq0(kr, 1)
            End If
        Next kr
    End If
End Sub