调试泥沙动力学的旧Fortran代码

Question

我正在看一些Fortran代码从旧的扫描纸。扫描质量不是很好，所以我可能复制错了。我试着用一个在线Fortran编译器来运行这个，但是它爆炸了。不熟悉Fortran，我想知道是否有人可以指出语法不合理的地方?代码来自一篇关于泥沙动力学的论文:

Komar, P.D.和Miller, m.c.， 1975。波浪作用下泥沙运动阈值与单向水流作用下泥沙运动阈值的比较及对阈值实际评价的讨论。沉积研究，45(1).

PROGRAM TSHOLD 
REAL LI, LO
G = 981.0 
PIE = 3.1416 
RHOW  = 1.00 
READ (6O,1) DIAM, RHOS
1 FORMAT (2X, F6.3,2X, F5.3) 
IF(DIAM .LT. 0.05) GO TO 5 
A = 0.463 * PIE 
B = 0.25 
GO TO 7 
5 A = 0.21 
B = 0.50 
7 PWR = 1.0 / (2.0 - B) 
FAC = (A * (RHOS - RHOW) * G/(RHOW * PIE**B))**PWR
FAC1 = FAC * DIAM**((1.0 - B) * PWR) 
T = 1.0 
15 J = 1.20 
LD = 156.13 * (T**2) 
UM = FAC1 * T**(B*PWR) 
WRITE(61,9) DIAM, T, UM
9 FORMAT(1H0, 10X, 17HGRAIN DIAMETER = ,F6.3,1X,2HCM //
1 11X, 14HWAVE PERIOD = ,F5.2, 1X, 3HSEC // 
2 11X, 22HORBITAL VELOCITY, UM = ,F6.2, 1X, 6HCM/SECl //
3 20X, 6HHEIGHT, 5X, 5HDEPTH, 8X, 3HH/L, 6X, 7HH/DEPTH //
4 22X, 2HCM, 8X, 2HCM /)
C INCREMENT WAVE HEIGHT, CALCULATE DEPTH 
H = 10.0 
DO 12 K = 1.60
SING = PIE * H / (UM * T) 
X = SING 
IF(X.LT.1.0) GO TO 30 
30 ASINH = X - 0.16666*X**3.0 + 0.07500* X ** 5.0 - 0.04464 * X ** 7.0
1 + 0.03038 * X ** 9.0 - 0.02237 * X ** 11.0 
32 LI = LD * (SINH(ASINH)/COSH(ASINH)) 
OPTH = ASINH * LI / 6.2832 
C CHECK WAVE STABILITY 
RATIO = H / DPTH 
IF(RATIO.GE.0.78) GO TO 11 
STEEP = H / LI 
TEST = 0.142 * (SINH(ASINH)/COSH(ASINH)) 
IF(STEEP.GE.TEST) GO TO 11 
WRITE(61,10) H, OPTH, STEEP, RATIO 
I0 FORMAT(IH0, 20X, F5.1, 4X, E9.3, 4X, F5.3, 4X, F4.2) 
11 H = H + 10.0 
12 CONTINUE 
T = T + 1.0 
15 CONTINUE 
END

Answer 1

问题更可能是旧Fortran需要固定的格式代码格式，其中语句前的空格数量非常重要。

以下是一些通用规则

• 正常语句从第7列及以上开始
• 行不能超过72列
• 列6上的任何字符表示该行是上一行的延续。我在上面的代码中看到9 FORMAT(..
后面的行
• 放在列1-5之间的数字表示标签，它可以是GO TO语句、DO语句或格式规范的目标。
• 第一列上的字符C，有时第一列上的任何字符表示该行是注释行。

参见https://people.cs.vt.edu/~asandu/Courses/MTU/CS2911/fortran_notes/node4.html获取更多信息。

根据上述规则，以下是如何输入代码，并使用正确的空格。我通过转换器运行F77代码，使其同时与F90和F77兼容。下面的代码可能现在使用在线编译器编译。

PROGRAM TSHOLD 
REAL LI, LO 
G = 981.0 
PIE = 3.1416 
RHOW  = 1.00 
READ (60,1) DIAM, RHOS 
1 FORMAT (2X, F6.3,2X, F5.3) 
IF(DIAM .LT. 0.05) GO TO 5 
A = 0.463 * PIE 
B = 0.25 
GO TO 7 
5 A = 0.21 
B = 0.50 
7 PWR = 1.0 / (2.0 - B) 
FAC = (A * (RHOS - RHOW) * G/(RHOW * PIE**B))**PWR 
FAC1 = FAC * DIAM**((1.0 - B) * PWR) 
T = 1.0 
DO 15 J=1,20 
LD = 156.13 * (T**2) 
UM = FAC1 * T**(B*PWR) 
WRITE(61,9) DIAM, T, UM 
9 FORMAT(1H0, 10X, 17HGRAIN DIAMETER = ,F6.3,1X,2HCM //             &
& 11X, 14HWAVE PERIOD = ,F5.2, 1X, 3HSEC //                        &
& 11X, 22HORBITAL VELOCITY, UM = ,F6.2, 1X, 6HCM/SECl //           &
& 20X, 6HHEIGHT, 5X, 5HDEPTH, 8X, 3HH/L, 6X, 7HH/DEPTH //          &
& 22X, 2HCM, 8X, 2HCM /)                                           
! INCREMENT WAVE HEIGHT, CALCULATE DEPTH                                
H = 10.0 
DO 12 K = 1,60 
SING = PIE * H / (UM * T) 
X = SING 
IF(X.LT.1.0) GO TO 30 
30 ASINH = X - 0.16666*X**3.0 + 0.07500* X ** 5.0 - 0.04464 * X ** 7.&
& + 0.03038 * X ** 9.0 - 0.02237 * X ** 11.0                       
32 LI = LD * (SINH(ASINH)/COSH(ASINH)) 
OPTH = ASINH * LI / 6.2832 
! CHECK WAVE STABILITY                                                  
RATIO = H / DPTH 
IF(RATIO.GE.0.78) GO TO 11 
STEEP = H / LI 
TEST = 0.142 * (SINH(ASINH)/COSH(ASINH)) 
IF(STEEP.GE.TEST) GO TO 11 
WRITE(61,10) H, OPTH, STEEP, RATIO 
10 FORMAT(G14.4, 20X, F5.1, 4X, E9.3, 4X, F5.3, 4X, F4.2) 
11 H = H + 10.0 
12 CONTINUE 
T = T + 1.0 
15 CONTINUE 
END                                           

我发现了几个转录错误，用点代替逗号，用字母O代替零，以及缺少DO语句。