Microsoft ExcelのVBAのプログラム
Option Explicit Private Const nix As Integer = 31 Private Const ndx As Integer = 8 Private Const nout As Integer = 20 Private h(nix) As Double Private z(nix) As Double Private hs(nix) As Double Private b(nix) As Double Private x(nix) As Double Private us(nix) As Double Private sib(nix, ndx) As Double Private sib0(nix, ndx) As Double Private a(nix) As Double Private u(nix) As Double Private si(nix) As Double Private tsm(nix) As Double Private dd(ndx) As Double Private d(ndx) As Double Private p0(nix, ndx) As Double Private p(nix, ndx) As Double Private w0(ndx) As Double Private uci(ndx) As Double Private tsci0(ndx) As Double Private tsi(nix, ndx) As Double Private wc(nix) As Double Private d10(nix) As Double Private d50(nix) As Double Private d90(nix) As Double Private dm(nix) As Double Private dzdt(nix) As Double Private ct(nix) As Double Private tscm(nix) As Double Private qb(nix) As Double Private qbi(nix, ndx) As Double Private qs(nix) As Double Private qsi(nix, ndx) As Double Private c(nix, ndx) As Double Private cb(nix, ndx) As Double Private dm0(nix) As Double Private z0(nix) As Double Private sn(nix) As Double Private dx(nix) As Double Public Sub pro1d() 'パラメータ Const g As Double = 9.8 Const s As Double = 1.65 Const ram As Double = 0.4 Const delta As Double = 0.1 Const bs0 As Double = 0.143 Const snu As Double = 0.000001 Const rs As Double = 1.436 Const sk As Double = 0.008 Const als As Double = 0.14 Dim i As Integer, j As Integer, k As Integer Dim ni As Integer ', dx As Double Dim ax As Double Dim sechour As Double Dim tuk As Double Dim etime As Double Dim dtmax As Double Dim al23 As Double Dim nd As Integer Dim ttime As Double Dim qq As Double Dim h00 As Double Dim d10_0 As Double, d50_0 As Double, d90_0 As Double Dim tt As Double Dim ucm As Double Dim aa As Double, bb As Double, cc As Double Dim fh As Double, fdh As Double, dh As Double Dim xi As Double Dim tsci As Double Dim usci As Double, xxi As Double Dim bs As Double Dim ome As Double Dim bet As Double Dim alfx As Double Dim amax As Double Dim dqb As Double Dim ah As Double Dim bh1 As Double, bh2 As Double, bh3 As Double Dim bh As Double, abh As Double Dim tmax As Double Dim dt As Double Dim dz As Double Dim dqbi As Double, dpi As Double Dim ncount As Long, mcount As Long, mcol As Integer Dim ws_out As Worksheet Dim ws_zd As Worksheet Dim ws_h As Worksheet Dim ws_z As Worksheet Dim ws_dm As Worksheet, ws_d10 As Worksheet, ws_d50 As Worksheet, ws_d90 As Worksheet Set ws_out = Sheets("計算結果") Set ws_zd = Sheets("河床変動量") Set ws_h = Sheets("水位") Set ws_z = Sheets("河床標高") Set ws_dm = Sheets("平均粒径") Set ws_d10 = Sheets("d10粒径") Set ws_d50 = Sheets("d50粒径") Set ws_d90 = Sheets("d90粒径") 'Call param '粒径分布は対数正規分布(8区分) For i = 1 To ndx dd(1) = 0.1 dd(2) = 0.2 dd(3) = 0.5 dd(4) = 0.7 dd(5) = 1# dd(6) = 2# dd(7) = 4# dd(8) = 10# Next i '初期値 Dim ws1 As Worksheet Set ws1 = Sheets("初期値") qq = ws1.Cells(3, "A").Text '500# h00 = ws1.Cells(3, "B").Text '3 'sn = ws1.Cells(3, "C").Text '0.025 tuk = ws1.Cells(3, "D").Text '5# etime = ws1.Cells(3, "E").Text '30# d10_0 = ws1.Cells(7, "A").Text '0.2 d50_0 = ws1.Cells(7, "B").Text '0.5 d90_0 = ws1.Cells(7, "C").Text '1# 'Call datainput '河道断面数と初期河床標高の読み込み Dim ws0 As Worksheet Set ws0 = Sheets("河道データ") ni = ws0.Cells(1, "B").Text '31 'dx = ws0.Cells(1, "C").Text '200 For i = 1 To ni 'read(17,*) b(i),z(i) b(i) = ws0.Cells(i + 2, "E").Text z(i) = ws0.Cells(i + 2, "D").Text dx(i) = ws0.Cells(i + 2, "C").Text sn(i) = ws0.Cells(i + 2, "F").Text z0(i) = z(i) d10(i) = d10_0 d50(i) = d50_0 d90(i) = d90_0 Next i '計算時間 ax = 1# / delta sechour = 3600# tuk = sechour * tuk etime = etime * sechour dtmax = tuk al23 = Log10(23#) nd = ndx - 1 For k = 2 To nd d(k) = (dd(k) + dd(k - 1)) * 0.5 / 1000# Next k For i = 1 To ni x(i) = (i - 1) * dx(i) Next i ' ' ------ initial grain size distribution ------ ' 'Call initgrain For i = 1 To ni p0(i, 1) = 0# p0(i, nd) = 1# For k = 2 To nd - 1 p0(i, k) = ppx(dd(k), d10(i), d50(i), d90(i)) Next k Next i For i = 1 To ni For k = 2 To nd sib(i, k) = p0(i, k) - p0(i, k - 1) sib0(i, k) = sib(i, k) Next k For k = 1 To nd p(i, k) = p0(i, k) Next k Next i ' ' ------ fall velocity ----- ' 'Call fallvel For k = 2 To nd If d(k) > 0.001 Then w0(k) = 32.8 * Sqr(d(k) * 100#) * 0.01 Else w0(k) = Sqr(2# / 3# * s * g * d(k) + _ (6# * snu / d(k)) ^ 2) - 6# * snu / d(k) End If Next k ' ' ------ critical shear velocity ------ ' 'Call critcalshear For k = 2 To nd Call usc(d(k), uci(k), s, snu, g) tsci0(k) = uci(k) ^ 2 / (s * g * d(k)) Next k ' '------ initial mean diameter ------ ' 'Call meandia For i = 1 To ni dm0(i) = 0# For k = 1 To nd dm0(i) = dm0(i) + d(k) * sib(i, k) Next k Next i ttime = 0# tt = tuk ' *********************************** ' 計算開始 ' *********************************** ncount = 0 '出力用のカウンタ mcol = 3 '出力用のカウンタ g100: Debug.Print "tt = "; tt ' '------ dm ------ ' 'Call meandia For i = 1 To ni dm(i) = 0# For k = 1 To nd dm(i) = dm(i) + d(k) * sib(i, k) Next k Next i ' ' ------ tau*cm ------ ' 'Call taucm For i = 1 To ni Call usc(dm(i), ucm, s, snu, g) tscm(i) = ucm ^ 2 / (s * g * dm(i)) Next i ' '------ h ------ ' 'Call height h(1) = h00 hs(1) = h(1) - z(1) For i = 2 To ni cc = -z(i - 1) + z(i) - hs(i - 1) - _ qq ^ 2 / (2# * g * b(i - 1) ^ 2 * hs(i - 1) ^ 2) _ - dx(i - 1) * sn(i - 1) ^ 2 * qq ^ 2 / _ - (2# * b(i - 1) ^ 2 * hs(i - 1) ^ (10# / 3#)) hs(i) = hs(i - 1) g120: aa = qq ^ 2 / (2# * g * b(i) ^ 2) bb = -dx(i) * sn(i) ^ 2 * qq ^ 2 / (2# * b(i) ^ 2) fh = hs(i) + aa / hs(i) ^ 2 + bb / hs(i) ^ (10# / 3#) + cc fdh = 1# - 2# * aa / hs(i) ^ 3 - 10# / 3# * bb / hs(i) ^ (13# / 3#) dh = -fh / fdh 'Debug.Print cc, hs(i), dh If Abs(dh) < 0.0001 Then GoTo g135: hs(i) = hs(i) + dh GoTo g120: g135: h(i) = z(i) + hs(i) Next i ' ' ------ a,u,si ------ ' 'Call ausi For i = 1 To ni a(i) = b(i) * hs(i) u(i) = qq / a(i) si(i) = (sn(i) * qq / b(i)) ^ 2 / hs(i) ^ 3.3333 us(i) = Sqr(g * hs(i) * si(i)) tsm(i) = us(i) ^ 2 / (s * g * dm(i)) Next i ' ' ------ qb (bed load) and qs (suspended load) ------ ' 'Call bedload For i = 1 To ni qb(i) = 0# qs(i) = 0# For k = 2 To nd tsi(i, k) = us(i) ^ 2 / (s * g * d(k)) xi = (al23 / Log10(21# * d(k) / dm(i) + 2#)) ^ 2 tsci = xi * tscm(i) If tsi(i, k) <= tsci Then qbi(i, k) = 0# qsi(i, k) = 0# Else usci = Sqr(tsci * s * g * d(k)) qbi(i, k) = sib(i, k) * 17# * Sqr(s * g * d(k) ^ 3) * _ tsi(i, k) ^ 1.5 * (1# - tsci / tsi(i, k)) * (1# - usci / us(i)) If qbi(i, k) < 1E-20 Then qbi(i, k) = 0# qb(i) = qb(i) + qbi(i, k) xxi = tsci / tsci0(k) bs = bs0 * xxi Call omega(bs, tsi(i, k), ome) qsi(i, k) = sib(i, k) * sk * (als * ome * s * g * d(k) / (rs * us(i)) - w0(k)) If qsi(i, k) < 1E-20 Then qsi(i, k) = 0# qs(i) = qs(i) + qsi(i, k) End If Next k Next i ' ' ------ cal. of c ------ ' 'Call calc For k = 2 To nd bet = 15# * w0(k) / us(ni) alfx = alf(bet) c(ni, k) = qsi(ni, k) / (w0(k) * alfx) cb(ni, k) = c(ni, k) * alfx For i = ni - 1 To 1 Step -1 bet = 15# * w0(k) / us(i) alfx = alf(bet) c(i, k) = (qq * c(i + 1, k) + dx(i) * b(i) * qsi(i, k)) / _ (qq + dx(i) * b(i) * w0(k) * alfx) If c(i, k) <= 0.0000000001 Then c(i, k) = 0# End If cb(i, k) = c(i, k) * alfx Next i Next k ' ' ------ wc, ct------ ' 'Call wcct For i = 1 To ni ct(i) = 0# wc(i) = 0# For k = 2 To nd ct(i) = ct(i) + c(i, k) wc(i) = wc(i) + w0(k) * cb(i, k) Next k Next i ' ' ------ cal. of dt ------ ' 'Call caldt amax = 0# For i = 1 To ni - 1 dqb = (qb(i + 1) * b(i + 1) - qb(i) * b(i)) / (b(i) * dx(i)) dzdt(i) = (dqb + wc(i) - qs(i)) / (1# - ram) For k = 2 To nd If (sib(i, k) = 0# Or sib(i + 1, k) = 0#) Then GoTo g130: ah = ax / (1# - ram) * qbi(i, k) / sib(i, k) bh1 = (qbi(i, k) / sib(i, k) - qbi(i + 1, k) / sib(i + 1, k)) / dx(i) bh2 = qbi(i, k) / (b(i) * sib(i, k) * dx(i)) * (b(i) - b(i + 1)) bh3 = qsi(i, k) / sib(i, k) bh = ax / (1# - ram) * (bh1 + bh2 + bh3) If dzdt(i) > 0# Then bh = bh + ax * dzdt(i) abh = ah / dx(i) + bh * 0.5 amax = max(amax, abh) g130: Next k Next i If amax > 0# Then tmax = 1# / amax dt = min(tmax, dtmax) Else dt = dtmax End If If tt + dt > tuk Then dt = tuk - tt ' '-------- output ------ ' 'Call dataoutput mcount = 0 '出力用のカウンタ If tt >= tuk Then Debug.Print "tt=", tt, "dt=", dt, "tuk=", tuk tt = 0# ' ' ------ p ------ ' For i = 1 To ni For k = 2 To nd - 1 p(i, k) = p(i, k - 1) + sib(i, k) Next k Next i ' ------ d10,d50,d90 ------ ' For i = 1 To ni For k = 2 To nd If (p(i, k) >= 0.1 And p(i, k - 1) <= 0.1) Then d10(i) = 10# ^ ( _ Log10(dd(k - 1)) + _ (0.1 - p(i, k - 1)) / (p(i, k) - p(i, k - 1)) * _ (Log10(dd(k)) - Log10(dd(k - 1)))) End If If (p(i, k) >= 0.5 And p(i, k - 1) <= 0.5) Then d50(i) = 10# ^ ( _ Log10(dd(k - 1)) + _ (0.5 - p(i, k - 1)) / (p(i, k) - p(i, k - 1)) * _ (Log10(dd(k)) - Log10(dd(k - 1)))) End If If (p(i, k) >= 0.9 And p(i, k - 1) <= 0.9) Then d90(i) = 10# ^ ( _ Log10(dd(k - 1)) + _ (0.9 - p(i, k - 1)) / (p(i, k) - p(i, k - 1)) * _ (Log10(dd(k)) - Log10(dd(k - 1)))) End If Next k Next i Debug.Print "time=", ttime / sechour '計算結果のシート出力 ncount = ncount + 1 ws_out.Cells(ncount, "A") = "time:" ws_out.Cells(ncount, "B") = Format(ttime / sechour, "######.###0") '計算結果の出力 ncount = ncount + 1 ws_out.Cells(ncount, "C") = "断面i" ws_out.Cells(ncount, "D") = "断面距離(m)" ws_out.Cells(ncount, "E") = "河床標高(m)" ws_out.Cells(ncount, "F") = "河川幅(m)" ws_out.Cells(ncount, "G") = "水位(m)" ws_out.Cells(ncount, "H") = "平均粒径(mm)" ws_out.Cells(ncount, "I") = "d10粒径(mm)" ws_out.Cells(ncount, "J") = "d50粒径(mm)" ws_out.Cells(ncount, "K") = "d90粒径(mm)" ws_out.Cells(ncount, "L") = "初期河床からの変動量(m)" '河床変動量のみ出力 mcount = mcount + 1 ws_zd.Cells(mcount, "A") = "断面i" ws_zd.Cells(mcount, "B") = "断面距離(m)" ws_zd.Cells(mcount, "C") = "河床標高(m)" '河床標高のみ出力 ws_z.Cells(mcount, "A") = "断面i" ws_z.Cells(mcount, "B") = "断面距離(m)" ws_z.Cells(mcount, "C") = "河床標高(m)" '水位のみ出力 ws_h.Cells(mcount, "A") = "断面i" ws_h.Cells(mcount, "B") = "断面距離(m)" ws_h.Cells(mcount, "C") = "河床標高(m)" '平均粒径のみ出力 ws_dm.Cells(mcount, "A") = "断面i" ws_dm.Cells(mcount, "B") = "断面距離(m)" 'd10粒径のみ出力 ws_d10.Cells(mcount, "A") = "断面i" ws_d10.Cells(mcount, "B") = "断面距離(m)" 'd50粒径のみ出力 ws_d50.Cells(mcount, "A") = "断面i" ws_d50.Cells(mcount, "B") = "断面距離(m)" 'd90粒径のみ出力 ws_d90.Cells(mcount, "A") = "断面i" ws_d90.Cells(mcount, "B") = "断面距離(m)" mcol = mcol + 1 For j = 1 To ni ncount = ncount + 1 mcount = mcount + 1 ws_out.Cells(ncount, "C") = j ws_out.Cells(ncount, "D") = Format(x(j), "####.##0") ws_out.Cells(ncount, "E") = Format(z(j), "####.##0") ws_out.Cells(ncount, "F") = Format(b(j), "####.##0") ws_out.Cells(ncount, "G") = Format(h(j), "####.##0") ws_out.Cells(ncount, "H") = Format(dm(j) * 1000#, "###.##0") ws_out.Cells(ncount, "I") = Format(d10(j), "####.##0") ws_out.Cells(ncount, "J") = Format(d50(j), "####.##0") ws_out.Cells(ncount, "K") = Format(d90(j), "####.##0") ws_out.Cells(ncount, "L") = Format(z(j) - z0(j), "####.##0") '河床変動量のみ出力 ws_zd.Cells(mcount, "A") = j ws_zd.Cells(mcount, "B") = Format(x(j), "####.##0") ws_zd.Cells(mcount, "C") = Format(z0(j), "####.##0") ws_zd.Cells(1, mcol) = "time:" & Format(ttime / sechour, "######.##0") ws_zd.Cells(mcount, mcol) = Format(z(j) - z0(j), "####.##0") '河床標高のみ出力 ws_z.Cells(mcount, "A") = j ws_z.Cells(mcount, "B") = Format(x(j), "####.##0") ws_z.Cells(mcount, "C") = Format(z0(j), "####.##0") ws_z.Cells(1, mcol) = "time:" & Format(ttime / sechour, "######.###0") ws_z.Cells(mcount, mcol) = Format(z(j), "####.##0") '水位のみ出力 ws_h.Cells(mcount, "A") = j ws_h.Cells(mcount, "B") = Format(x(j), "####.##0") ws_h.Cells(mcount, "C") = Format(z0(j), "####.##0") ws_h.Cells(1, mcol) = "time:" & Format(ttime / sechour, "######.###0") ws_h.Cells(mcount, mcol) = Format(h(j), "####.##0") '平均粒径のみ出力 ws_dm.Cells(mcount, "A") = j ws_dm.Cells(mcount, "B") = Format(x(j), "####.##0") ws_dm.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0") ws_dm.Cells(mcount, mcol - 1) = Format(dm(j) * 1000#, "####.##0") 'd10粒径のみ出力 ws_d10.Cells(mcount, "A") = j ws_d10.Cells(mcount, "B") = Format(x(j), "####.##0") ws_d10.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0") ws_d10.Cells(mcount, mcol - 1) = Format(d10(j), "####.##0") 'd50粒径のみ出力 ws_d50.Cells(mcount, "A") = j ws_d50.Cells(mcount, "B") = Format(x(j), "####.##0") ws_d50.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0") ws_d50.Cells(mcount, mcol - 1) = Format(d50(j), "####.##0") 'd90粒径のみ出力 ws_d90.Cells(mcount, "A") = j ws_d90.Cells(mcount, "B") = Format(x(j), "####.##0") ws_d90.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0") ws_d90.Cells(mcount, mcol - 1) = Format(d90(j), "####.##0") Next j End If ' ' ------ cal. of z, sib ------ ' 'Call calzsib For i = 1 To ni - 1 dz = dt * dzdt(i) z(i) = z(i) + dz For k = 2 To nd If dz > 0# Then xi = sib(i, k) Else xi = sib0(i, k) End If dqbi = (qbi(i + 1, k) * b(i + 1) - qbi(i, k) * b(i)) _ / (b(i) * dx(i)) dpi = ax * (dt / (1# - ram) * _ (dqbi + w0(k) * cb(i, k) - qsi(i, k)) - xi * dz) sib(i, k) = sib(i, k) + dpi If sib(i, k) < 0.0000000001 Then sib(i, k) = 0# Next k Next i ttime = ttime + dt tt = tt + dt If ttime > etime Then GoTo g999: GoTo g100: g999: Close ' *********************************** ' 計算終了 ' *********************************** End Sub ' ' *********************************** ' Private Sub usc(d As Double, uc As Double, s As Double, snu As Double, g As Double) Dim rst As Double rst = Sqr(s * g) * d ^ (3# / 2#) / snu If rst <= 2.14 Then uc = 0.14 * s * g * d If rst > 2.14 Then uc = (0.1235 * s * g) ^ (25# / 32#) *_ snu ^ (7# / 16#) * d ^ (11# / 32#) If rst > 54.2 Then uc = 0.034 * s * g * d If rst > 162.7 Then uc = (0.01505 * s * g) ^ (25# / 22#) * _ snu ^ (-3# / 11#) * d ^ (31# / 22#) If rst > 671# Then uc = 0.05 * s * g * d uc = Sqr(uc) End Sub ' ' *********************************** ' Private Sub omega(bs As Double, tsi As Double, ome As Double) Dim pai As Double Dim a1 As Double, a2 As Double, a3 As Double Dim ad As Double Dim t As Double, zx As Double, px As Double Dim er1 As Double, er As Double Dim xxx As Double a1 = 0.4361836 a2 = -0.1201676 a3 = 0.937298 pai = 3.141592 If tsi <= 0.0000001 Then ome = 0# Exit Sub End If ad = bs / tsi - 2# If ad >= 0# Then xxx = ad * Sqr(2#) End If If ad < 0# Then xxx = -ad * Sqr(2#) End If t = 1# / (1# + 0.33627 * xxx) zx = 1# / Sqr(2# * pai) * Exp(-xxx ^ 2# / 2#) px = 1# - zx * (a1 * t + a2 * t ^ 2 + a3 * t ^ 3) er1 = 2# - 2# * px If ad >= 0# Then er = er1 / 2# End If If ad < 0# Then er = (2# - er1) / 2# End If If (bs = 0# Or er = 0#) Then ome = 0# Else ome = tsi / bs / (2# * Sqr(pai)) * Exp(-ad ^ 2) / er + tsi * 2# / bs - 1# End If 'if(ome.lt.1e-10) ome=0. End Sub ' ' *********************************** ' Private Function alf(bet As Double) 'if bet > 20. then ' alf = bet 'else alf = bet / (1# - Exp(-bet)) 'end if End Function ' ' *********************************** ' Private Function ppx(dd As Double, d10 As Double, d50 As Double, d90 As Double) Dim s1 As Double, s2 As Double Dim x As Double s1 = 1.282 / (Log10(d90) - Log10(d50)) s2 = -1.282 / (Log10(d10) - Log10(d50)) If dd >= d50 Then x = (Log10(dd) - Log10(d50)) * s1 Else x = (Log10(dd) - Log10(d50)) * s2 End If ppx = pc(x) End Function ' ' *********************************** ' Private Function tfunc(x As Double) If x >= 0# Then tfunc = 1# / (1# + 0.33267 * x) Else tfunc = 1# / (1# - 0.33267 * x) End If End Function ' ' *********************************** ' Private Function pc(x As Double) Dim pai As Double Dim a1 As Double, a2 As Double, a3 As Double Dim tx As Double pai = 3.141592 a1 = 0.4361836 a2 = -0.1201676 a3 = 0.937298 tx = tfunc(x) If x >= 0# Then pc = 1# - 1# / Sqr(2# * pai) * _ Exp(-x ^ 2 / 2#) * (a1 * tx + a2 * tx ^ 2 + a3 * tx ^ 3) Else pc = 1# / Sqr(2# * pai) * _ Exp(-x ^ 2 / 2#) * (a1 * tx + a2 * tx ^ 2 + a3 * tx ^ 3) End If End Function ' ' *********************************** ' Private Function max(a As Double, b As Double) max = a If max < b Then max = b End Function ' ' *********************************** ' Private Function min(a As Double, b As Double) min = a If min > b Then min = b End Function Private Function Log10(ByVal Val As Double) Log10 = Log(Val) / Log(10) End Function