[推荐]促淤围垦的龙口水力计算BAS程序

促淤围垦的龙口水力计算BAS程序
n1-涨潮落潮数据个数，n2-水位库容数据个数，n3,n4糙率等，不改，m流量系数
5 : CLEAR
6 : PRINT "*********************************************"
8 : INPUT "工程名称"; a$
10 : PRINT a$; "水库截流水力计算"
20 : n1 = 13: n2 = 8: n3 = 21: n4 = 21: m = .32
30 : DIM o(n4), y(n4)
40 : INPUT "立堵/平堵(1/2)="; a1
45 : IF a1 = 1 THEN PRINT "截流方式:立堵"
48 : IF a1 = 2 THEN PRINT "截流方式:平堵"
50 : INPUT "堰顶高程yd="; yd
51 : INPUT "打印时间dt="; dt
52 : INPUT "打印时间段内计算次数n5="; n5
53 : INPUT "堰前水宽b0="; b0
54 : INPUT "涨潮历时s1="; s1
55 : INPUT "落潮历时s2="; s2
60 : INPUT "计算开始潮位h0="; h0
61 : INPUT "计算开始潮时t0="; t0
62 : PRINT "堰顶高程yd="; yd
63 : PRINT "堰前水宽b0="; b0
65 : PRINT "涨潮历时s1="; s1; TAB(9); "落潮历时s2="; s2: s1 = s1 - t0
68 : INPUT "堰宽b="; b
69 : PRINT "堰宽b="; b
70 : PRINT "计算开始潮位h0="; h0; TAB(9); "计算开始潮时t0="; t0
71 : LET a2 = 0: IF b < 0 THEN END
72 : PRINT "打印时间内计算次数n5="; n5
74 : PRINT "*+-+-+-+-+-+-+-+-+-+-+-+-+*"
75 : T1 = t0: H1 = h0: x1 = h0
80 : FOR i1 = 1 TO INT((s1 + s2) / dt) * n5
90 : T1 = T1 + dt / n5
100 : LET ss = 1
110 : LET q = T1
112 : GOSUB 3900
115 : LET h2 = y
120 : IF i1 > 1 GOTO 150
130 : LET ss = 2: LET q = H1: GOSUB 3900
140 : LET v1 = y
150 : K4 = (x1 * (1 - a2) + h2 * a2 - yd) / (x1 * a2 + h2 * (1 - a2) - yd)
160 : IF K4 < .8 THEN GOTO 210
180 : LET ss = 3
190 : LET q = (x1 * (1 - a2) + h2 * a2 - yd) / (x1 * a2 + h2 * (1 - a2) - yd): GOSUB 3900
200 : K1 = y
210 : K2 = (h2 * (1 - a2) + x1 * a2 - yd) / b
220 : IF K2 > 1 THEN LET K2 = 1
230 : K3 = 1 - .2 * K2 * .7
240 : IF K4 < .8 THEN LET K1 = 1
250 : QQ = K1 * K3 * m * b * (2 * 9.81) ^ .5 * (h2 * (1 - a2) + x1 * a2 - yd) ^ 1.5
260 : HK = (1.05 * QQ ^ 2 / b ^ 2 / 9.81) ^ (1 / 3)
270 : IF K4 < .8 THEN LET V = QQ / (b * HK)
280 : IF K4 < .8 GOTO 300
290 : V = QQ / (b * (x1 * (1 - a2) + h2 * a2 - yd))
300 : IF V < vm GOTO 308
305 : vm = V: tm = T1: xm = x1: qm = QQ: hm = h2
308 : IF T1 / dt <> INT(T1 / dt) AND i1 > 1 GOTO 342
310 : PRINT "V="; INT(V * 1000) / 1000; TAB(9); "T="; INT(T1 * 1000) / 1000
320 : PRINT "XS="; INT(x1 * 1000) / 1000; TAB(9); "SS="; INT(h2 * 1000) / 1000
330 : PRINT "Q="; INT(QQ * 1000) / 1000
342 : IF i1 <> INT((s1 / dt) * n5 + 1) GOTO 350
344 : PRINT "涨潮时段最大流速VM="; vm
345 : PRINT "最大流速对应潮时TM="; tm
346 : PRINT "涨潮最大流速出现时的库内水位XM="; xm
347 : PRINT "最大流速对应流量QM="; qm
348 : PRINT "最大流速对应潮位HM="; hm
349 : vm = 0: tm = 0: xm = 0: qm = 0: hm = 0
350 : IF a2 = 1 GOTO 370
352 : IF i1 > INT(s1 / dt) * n5 AND x1 > h2 THEN a2 = 1
370 : v1 = v1 + QQ * dt / n5 * (-1) ^ a2 * 60
380 : LET q = v1
390 : LET ss = 4: GOSUB 3900
394 : x1 = y
396 : IF (x1 * (1 - a2) + h2 * a2 - x1 * a2 - h2 * (1 - a2)) <= 0 GOTO 400
397 : LET ss = 2: LET q = h2: GOSUB 3900
398 : LET v1 = y: LET x1 = h2
400 : NEXT i1
420 : PRINT "落潮时段最大流速vm="; vm
430 : PRINT "最大流速对应潮时tm="; tm
440 : PRINT "涨潮最大流速出现时的库内水位xm="; xm
450 : PRINT "最大流速对应流量qm="; qm
460 : PRINT "最大流速对应潮位hm="; hm
470 : IF a1 = 1 THEN b = b - 20
480 : IF a1 = 1 GOTO 70
500 : END
3900 : IF ss = 1 THEN RESTORE 5000
3910 : IF ss = 1 THEN ii = n1
3920 : IF ss = 2 THEN ii = n2
3930 : IF ss = 2 THEN RESTORE 6000
3940 : IF ss = 3 THEN RESTORE 7000
3950 : IF ss = 3 THEN ii = n3
3960 : IF ss = 4 THEN ii = n2
3970 : IF ss = 4 THEN RESTORE 6000
4000 : FOR i = 0 TO ii
4002 : IF ss <> 4 GOTO 4010
4005 : READ y(i), o(i): GOTO 4020
4010 : READ o(i), y(i)
4020 : NEXT i
4030 : FOR i = 2 TO ii - 1
4040 : IF i = ii - 1 GOTO 4080
4050 : IF q < o(i) GOTO 4070
4060 : NEXT