====== 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**」の使用方法を説明します。
ダウンロードについては、[[http://www.disaster-software.net|水理水文ソフトウェア]]のHPを参照してください。
===== マクロの説明 =====
本マクロは、各条件シートに条件を入力するだけで、**流域内の斜面と河道の流出解析計算**が行えます。
tys-kinema_v_1_1_0モデルの特徴は
- 支川の合流も可能(流域内の河道網を表現して、一括で流出解析ができる)。
- 斜面から降雨による流入量を計算する。
- 表層(中間層の貯留効果を計算が可能である。
==== 計算条件とデータの入力 ====
=== 計算条件 ===
雨量データ数や河道数などの計算条件を入力します。
{{ :流域流出計算:image12.jpg?650|計算条件シート}}
;#;
計算条件シート
;#;
* パラメータ
* 総雨量データ数
* 河道数
* 継続時間区分数
* 雨量の表面流出率
* 雨量観測所数
* 斜面表層の空隙率
* 表層厚(m)
* 継続時間区分に応じた時間緒元
* 番号
* 継続時間(s)
* 計算刻み時間(s)
* 出力設定(流量のみ)
* 出力箇所数のデータ数
* 出力ピッチ(s) : 出力時間間隔を秒で指定
* 出力箇所(10カ所まで)
* 番号
* 河道or斜面
* 河道or斜面番号
* 断面or斜面番号
* 名称
* 流出解析ボタンをクリックすると計算を開始する。
* 計算開始後、**計算河道流量**、**計算河道水深**、**計算斜面流量**、**モニター地点流量**、**モニター地点流量時系列**のシートが自動的に作成されます。
再度解析を行った場合は、計算結果シートの値が上書きされます。
=== 雨量データ ===
雨量データを入力するシートです。最大10箇所の雨量データ観測所場所を登録する事ができます。
{{ :流域流出計算:image13.jpg?650|雨量シート}}
;#;
計算条件シート
;#;
* 継続時間区分番号
* 「計算条件」シートの継続時間区分に応じた時間緒元の番号を入力します。
=== 河道データ ===
{{ :流域流出計算:image15.jpg?650|河道シート}}
;#;
流域諸元シート
;#;
=== 斜面データ ===
{{ :流域流出計算:image16.jpg?650|斜面シート}}
;#;
斜面シート
;#;
シートの黄色の網掛けセル部分にデータを入力します。水色のの網掛けセル部分は自動的に計算されます。
* 河道番号
* 斜面番号
* 流域名称等
* 流入断面番号(下流側)
* 流入断面番号(上流側)
* 流域面積
* 流路長
* 高低差
* 下端幅
* 断面数
* 等価粗度
* 飽和透水係数
* 観測所番号
* 斜面標高上
* 斜面標高下
=== 流域諸元データ ===
河道網データを作成します。
{{ :流域流出計算:image14.jpg?650|流域諸元シート}}
;#;
流域諸元シート
;#;
==== 計算結果 ====
{{ :流域流出計算:image17.jpg?650|計算河道流量シート}}
;#;
計算河道流量シート(本川)
;#;
{{ :流域流出計算:image18.jpg?650|計算河道水深シート}}
;#;
計算河道水深シート(本川)
;#;
{{ :流域流出計算:image19.jpg?650|計算斜面流量シート}}
;#;
計算斜面流量シート
;#;
{{ :流域流出計算:image20.jpg?650|モニター地点流量シート}}
;#;
モニター地点流量シート
;#;
{{ :流域流出計算:image21.jpg?650|モニター地点流量時系列図シート}}
;#;
モニター地点流量時系列図
;#;
===== VBAソースファイル =====
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