tys-kinema_v_1_1_0

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

本マクロは、各条件シートに条件を入力するだけで、流域内の斜面と河道の流出解析計算が行えます。

tys-kinema_v_1_1_0モデルの特徴は

1. 支川の合流も可能（流域内の河道網を表現して、一括で流出解析ができる）。
2. 斜面から降雨による流入量を計算する。
3. 表層（中間層の貯留効果を計算が可能である。

雨量データ数や河道数などの計算条件を入力します。

• パラメータ
  • 総雨量データ数
  • 河道数
  • 継続時間区分数
  • 雨量の表面流出率
  • 雨量観測所数
  • 斜面表層の空隙率
  • 表層厚(m)

• 継続時間区分に応じた時間緒元
  • 番号
  • 継続時間(s)
  • 計算刻み時間(s)

• 出力設定（流量のみ）
  • 出力箇所数のデータ数
  • 出力ピッチ(s) : 出力時間間隔を秒で指定
• 出力箇所（10カ所まで）
  • 番号
  • 河道or斜面
  • 河道or斜面番号
  • 断面or斜面番号
  • 名称
• 流出解析ボタンをクリックすると計算を開始する。
• 計算開始後、計算河道流量、計算河道水深、計算斜面流量、モニター地点流量、モニター地点流量時系列のシートが自動的に作成されます。　

雨量データを入力するシートです。最大１０箇所の雨量データ観測所場所を登録する事ができます。

• 継続時間区分番号
  • 「計算条件」シートの継続時間区分に応じた時間緒元の番号を入力します。

シートの黄色の網掛けセル部分にデータを入力します。水色のの網掛けセル部分は自動的に計算されます。

• 河道番号
• 斜面番号
• 流域名称等
• 流入断面番号(下流側）
• 流入断面番号(上流側）
• 流域面積
• 流路長
• 高低差
• 下端幅
• 断面数
• 等価粗度
• 飽和透水係数
• 観測所番号
• 斜面標高上
• 斜面標高下

河道網データを作成します。

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

Option Explicit

Public ttt As Long
Public fit As Long
Public fidt  As Long

Public windowflag As Boolean

Public Sub kinemamain()
    'キネマによる流出解析
    Dim ir As Integer, it As Integer
    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 sqcount As Long
   
    '計算データの読み込み
    Call datainput(flag)
    If flag = False Then End

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

    irr = 0
    dtl = 0#    '経過時間(s)
    qcount = 0
    hcount = 0
    sqcount = 0
       
    '初期値の出力
    Call SheetWriteMatrixQ(irr, dtl, qcount)
    Call SheetWriteMatrixH(irr, dtl, hcount)
    Call SheetWriteMatrixSQ(irr, dtl, sqcount)
    Call sheetwrite1(irr, dtl)
    
    For ir = 1 To nndata   '雨量データ数

        ii = khy1(ir)
        dtt = dt(ii)
        idt = CInt(delt(ii) / dtt)
        
        ttt = ir
        
        For it = 1 To idt      '計算刻み時間

            dtl = dtl + dtt
            
            'フォームの計算時刻表示
            fit = it
            fidt = idt          
            
            '流出解析
            For kr = km To 1 Step -1     '河道数
                Call slope(ir, it, kr) '斜面流出解析
                Call river(ir, it, kr) '河道流出解析
            Next kr
            Call reinitial      '変数の更新
            
            '指定時間ピッチで出力
            If (CLng(dtl) Mod dtout) = 0 And (dtl - CLng(dtl) = 0#) Then
                irr = irr + 1
                Call sheetwrite1(irr, dtl)
                Call SheetWriteMatrixQ(irr, dtl, qcount)
                Call SheetWriteMatrixH(irr, dtl, hcount)
                Call SheetWriteMatrixSQ(irr, dtl, sqcount)
            End If
        Next it
    Next ir
    
    '最終結果の出力
    Call モニター地点ハイドログラフ作成(irr, dtl)
    
 
    Close

End Sub


Private Sub initial(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)
            qrssig(k, l) = 0#
        Next l
        
        For j = 1 To ns(k)
            ar1(k, j) = 0#
            qqsig(k, j) = 0#
            qqq0(k, j) = 0#
            qqq1(k, j) = 0#
        Next j
    Next k
    
    flag = True
    Exit Sub
dataerror:
    MsgBox "初期化ルーチンでエラーが発生しました。" & vbNewLine _
                    & "処理を中断します。", vbCritical
    flag = False
End Sub

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


Private Sub slope(ir As Integer, it As Integer, k As Integer)
    '斜面の計算
    Dim ddz(ksd) As Double, sns(ksd) As Double, wbu(ksd) As Double, wbd(ksd) As Double
    Dim j As Integer, l As Integer
    Dim iii As Integer
    Dim dx As Double, xn As Double, dk As Double
    Dim dl1 As Double, dl2 As Double, qrain As Double
    Dim re1 As Double, re2 As Double, dhs1 As Double
    Dim aaa As Double, timerate As Double, hh2 As Double
        
    For j = 1 To ks(k)
        '降雨がなく､前回ステップで斜面断面全部の水深が0の時､計算をスキップ
        If rain(ir, khy2(k, j)) <= 0# Then
            For l = 1 To ids(k, j)
                If hs1(k, j, l) > 0# Then GoTo Skip1
                If hs2(k, j, l) > 0# Then GoTo Skip1
            Next l
            For l = 1 To ids(k, j)
                qrs(k, j) = 0#
                qs1(k, j) = 0#
                qs2(k, j) = 0#
            Next l
            GoTo skip2
        End If
Skip1:
        '-------------- setting data -----------------------
        dx = xls(k, j) / CDbl(ids(k, j))    '斜面セル長さ
        xn = xns(k, j)                      '等価粗度
        dk = xns2(k, j)                     '飽和透水係数
        '----- setting mesh data
        For l = 1 To ids(k, j)
            dl1 = dx * CDbl(l - 1)
            dl2 = dx * CDbl(l)
            '単位セルの下端・上端幅
            '台形形状考慮-------------------------------------------
            wbd(l) = bds(k, j) + (bus(k, j) - bds(k, j)) * dl1 / xls(k, j)
            wbu(l) = bds(k, j) + (bus(k, j) - bds(k, j)) * dl2 / xls(k, j)
            '-------------------------------------------------------
            ddz(l) = (zus(k, j) - zds(k, j)) / CDbl(ids(k, j))      '斜面セル高さ
            sns(l) = ddz(l) / Sqr(dx * dx + ddz(l) * ddz(l))        '斜面セル正弦(sin)
        Next l
        
        '-----calculation of discharge
        qrain = rain(ir, khy2(k, j)) / CDbl(idt) / 1000#        '降雨強度(mm/h-->m/s)

        For l = ids(k, j) To 1 Step -1
            ' -----前回の計算値より変数定義
            hs(l) = hs1(k, j, l)
            qs(l) = hs(l) ^ (5# / 3#) * Sqr(sns(l)) / xn
            qs0(l) = qs(l) * wbd(l)

            hsm(l) = hs2(k, j, l)
            qsm(l) = hsm(l) * sns(l) * dk
            qsm0(l) = qsm(l) * wbd(l)
            aaa = (wbd(l) + wbu(l)) * dx / 2#   '単位斜面セル面積

            re1 = qrain
            re2 = 0#
            
            '中間流出成分
            If dmm > 0 Then
                If l <> ids(k, j) Then
                    hs2(k, j, l) = hsm(l) - ((qsm0(l) - qsm0(l + 1)) / aaa * dtt - re2) / ramda
                Else
                    hs2(k, j, l) = hsm(l) - (qsm0(l) / aaa * dtt - re2) / ramda
                End If
                If hs2(k, j, l) < 0# Then hs2(k, j, l) = 0#
                '表面・中間の水深の調節
                If hs2(k, j, l) > dmm Then
                    '中間分を表面へ(浸み出し)
                    dhs1 = (hs2(k, j, l) - dmm) * ramda
                    hs2(k, j, l) = dmm
                Else
                    dhs1 = 0#
                End If
            Else
                hs2(k, j, l) = 0#
                dhs1 = 0#
            End If

            '表面流出成分
            If l <> ids(k, j) Then
                hs1(k, j, l) = hs(l) - (qs0(l) - qs0(l + 1)) / aaa * dtt + re1 + dhs1
            Else
                hs1(k, j, l) = hs(l) - qs0(l) / aaa * dtt + re1 + dhs1
            End If
            If hs1(k, j, l) < 0# Then hs1(k, j, l) = 0#
        Next l
        
        qs1(k, j) = qs0(1)
        qs2(k, j) = qsm0(1)
        qrs(k, j) = qs0(1) + qsm0(1)
        qrssig(k, j) = qrssig(k, j) + qrs(k, j) * dtt
        
skip2:
    Next j

End Sub

Private Sub river(ir As Integer, it As Integer, k As Integer)
    '河道の計算
    Dim j As Integer, kr2 As Integer
    Dim volv As Double, dx As Double
    Dim volw0 As Double, volw1 As Double
    Dim voldsl As Double, voldsf As Double
    Dim dsl As Double, dsf As Double
    Dim timerate As Double
    Dim vdosyal As Double, vdosyaf As Double
    Dim qiil As Double, qiif As Double
    Dim voltq As Double, deltz As Double
    Dim aa As Double, bb As Double, cc As Double, ee As Double
    
    '流量と流積の関係を求める - ---計算の開始時のみ
    If ir = 1 And it = 1 Then Call secof(k)

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

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

        If sita < 0# Then sita = 0.00001    '逆勾配になっても計算が止まらないように
        

        '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        
        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) - 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))
        '--------------------------------
    
    Next j

    For j = 1 To ns(k)
        If ar1(k, j) < 0# Then
            dep(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
        End If
    Next j
    
End Sub
'
Private Function calq(k As Integer, j As Integer, aa As Double)
    '運動方程式より流量qを計算する
    Dim rr As Double, ccl As Double, ss As Double, dd As Double, ee As Double
    Dim ccc As Double, pp As Double
    
    rr = ck1(k, j) * (aa ^ (cz(k, j))) '径深
    ss = sita
    dd = dep(k, 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# 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 kds As Integer
    
    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
    
    '斜面からの流入量
    For l = 1 To ks(k)
        If ksj1(k, l) <= j And ksj2(k, l) >= j Then
            If j = ns(k) Then
                qpls(k, j) = qpls(k, j) + qrs(k, l)
            ElseIf ksj1(k, l) = ksj2(k, l) Then
                qpls(k, j) = qpls(k, j) + qrs(k, l)
            Else
                qpls(k, j) = qpls(k, j) + qrs(k, l) * rl(k, j + 1) / slen(k, l)
            End If
        End If
    Next l
    
End Sub