Fortran高斯烟羽模型代码do循环仅执行第一步卡住问题排查

• 赋值顺序错误
  你在代码中先执行了dy = dx、dz = dx，之后才给dx赋值为2 * kx * t。此时dx未初始化，值为内存随机数，极大概率出现无效浮点数（比如极大值、NaN），导致后续计算卡死。
• 整数除法错误
  代码中(1/2)是整数运算，Fortran中整数相除结果会截断为整数，最终得到的结果是0而非预期的0.5，会导致所有f2计算结果恒为exp(0) = 1，属于逻辑错误。
• 输出格式冗余导致写入速度极慢
  你使用了F50.48的输出格式，要求每个浮点数输出48位小数，但单精度REAL类型的有效位数仅为6~7位，冗余的格式会让单次写入的内容大幅增加。你的三层循环总执行次数为101 * 101 * 51 = 520251次，过慢的写入速度会让程序看起来处于卡住状态。
• 逻辑缺失（非卡死直接原因，但不符合高斯烟羽计算要求）
  你的三层循环中x变量完全没有参与浓度计算，标准高斯烟羽模型的扩散参数σy、σz都和下风距离x相关，当前代码所有x位置的浓度输出完全一致，不符合需求。

PROGRAM plume
IMPLICIT NONE
REAL,PARAMETER :: pi=3.141592, q=1000, u=5, t=10, &
                  kx=1, ky=1, kz=1
REAL :: dx, dy, dz, f2, C
INTEGER :: x, y, z
! 声明函数返回值类型
REAL :: f1

! 先给dx赋值，再给dy、dz赋值
dx = (2 * kx * t)
dy = dx
dz = dx

OPEN(10, file='./results/plume.txt')
! 修改输出格式为合理范围
WRITE(10,'(3(A4),2(x,A15))') "x","y","z", "f2", "C"

PRINT*, "f1=", f1(q,pi,u,dx)

DO x = 0, 1000, 10
 ! 可在这里补充σy、σz和x的关联逻辑，比如 sigma_y = ky*sqrt(real(x)) 之类的，根据你的模型调整
 DO y = -500,500,10 
  DO z = 0,500,10 
   ! 把1/2改成1.0/2.0，用浮点数运算
   f2 = exp( -(1.0/2.0) * ( ( (real(y)**2.) + (real(z)**2.) ) / dx ) ) 
   C = f1(q,pi,u,dx) * f2
   ! 修改输出格式为合理的浮点格式
   WRITE(10,'(3(I4),2(x,E15.7))') x, y, z, f2, C
  END DO
 END DO
END DO
     

END PROGRAM


FUNCTION f1(q,pi,u,dx)
IMPLICIT NONE
REAL :: q, pi, u
REAL :: f1, dx
 f1 = q / (2 * pi * u * dx)
END FUNCTION f1

内容的提问来源于stack exchange，提问作者won