[求助] 各位大佬，最近刚学lbm，利用Fortran进行编程，但是代码输...

parameter (n=1000,m=40)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
real ssum, usum, vsum
integer i
open(2,file='uvfield')
open(3,file='uvely')
open(4,file='vvelx')
open(8,file='timeu')
uo=0.20
sumvelo=0.0
rhoo=5.00 ! 虽然不知道为什么在后边密度的计算没有加上rhoo
dx=1.0
dy=dx
dt=1.0
alpha=0.02
Re=uo*m/alpha
print *, "Re=", Re
omega=1.0/(3.*alpha+0.5)
mstep=40000      ! 时间步长的确定需要考虑
cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
! 初始条件的确定  对于初始条件的设定，是流体没有进入空间之前的设定
do j=21,m
do i=0,40
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
    end do
do j=0,m
    do i=41,n
        rho(i,j)=rhoo
        u(i,j)=0.0
        v(i,j)=0.0
    end do
end do
! 边界条件设定   这个是流体没有进入前，也就是流体状态的设定
do j=21,m
    u(0,j)=uo
    v(0,j)=0.0
    end do
! main loop
1 do kk=1,mstep
call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
call streaming(f,n,m)
! ——————————–
call sfbound(f,n,m,uo)
call rhouv(f,rho,u,v,cx,cy,n,m)
!print *, u(0,m/2),v(0,m/2),rho(0,m/2),u(n,m/2),v(n,m/2),rho(n,m/2)
write(8,*) kk,u(n/2,m/2),v(n/2,m/2)   
    end do
! end of the main loop
call result(u,v,rho,uo,n,m)
stop
end  
subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
do j=21,m    ! j=20是代表障碍物，21代表了流体
    do i=0,40
        t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
DO k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
END DO
END DO
END DO
do j=0,m
    do i=41,n
        t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
DO k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
END DO
END DO
END DO
return
    end 
subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
! 没有障碍物的流情况 这一部分存在着命名的重叠，先写后改 这一部分的节点是有问题的
! 存在障碍物的部分流
do j=22,m
    do i=40,1,-1
        f(1,i,j)=f(1,i-1,j)
    end do
    do i=0,39
        f(3,i,j)=f(3,i+1,j)
    end do
end do
do j=40,22,-1
    do i=0,40
        f(2,i,j)=f(2,i,j-1)
    end do
    do i=40,1,-1
        f(5,i,j)=f(5,i-1,j-1)
    end do
    do i=0,39
        f(6,i,j)=f(6,i+1,j-1)
    end do
end do
do j=21,39
    do i=0,40
        f(4,i,j)=f(4,i,j+1)
    end do
    do i=0,39
        f(7,i,j)=f(7,i+1,j+1)
    end do
    do i=40,1,-1
        f(8,i,j)=f(8,i-1,j+1)
    end do
end do
DO j=0,m
DO i=n,42,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=41,n-1 !LEFT TO RIGHT
f(3,i,j)=f(3,i+1,j)
END DO
END DO
DO j=m,1,-1 !TOP TO BOTTOM
DO i=41,n
f(2,i,j)=f(2,i,j-1)
END DO
DO i=n,42,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO i=41,n-1
f(6,i,j)=f(6,i+1,j-1)
END DO
END DO
DO j=0,m-1 !BOTTOM TO TOP
DO i=41,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=41,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,42,-1
f(8,i,j)=f(8,i-1,j+1)
END DO
END DO
return
    end 
subroutine sfbound(f,n,m,uo)
real f(0:8,0:n,0:m)
real uo,rhow
! 有障碍物和没有障碍物的底边条件    这一部分没验证输出数据
do i=0,n
    f(2,i,0)=f(4,i,0)
    f(5,i,0)=f(7,i,0)
    f(6,i,0)=f(8,i,0)
    end do
! 存在障碍物北边界条件
do i=0,40
    f(2,i,20)=f(4,i,20)
    f(5,i,20)=f(7,i,20)
    f(6,i,20)=f(8,i,20)
end do
! 存在障碍物右边界条件
do j=0,20
    f(1,40,j)=f(3,40,j)
    f(5,40,j)=f(7,40,j)
    f(8,40,j)=f(6,40,j)
end do
! 左边没有障碍物的已知速度边界条件
do j=21,m
    rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))/(1.-uo)
    f(1,0,j)=f(3,0,j)+2.*rhow*uo/3.
    f(5,0,j)=f(7,0,j)+rhow*uo/6.
    f(8,0,j)=f(6,0,j)+rhow*uo/6.
end do
! 上方边界
do i=0,n
    f(4,i,m)=f(2,i,m)
    f(7,i,m)=f(5,i,m)
    f(8,i,m)=f(6,i,m)
end do
! 没有障碍物的右边界
do j=1,m
    f(1,n,j)=2.*f(1,n-1,j)-f(1,n-2,j)
    f(5,n,j)=2.*f(5,n-1,j)-f(5,n-2,j)
    f(8,n,j)=2.*f(8,n-1,j)-f(8,n-2,j)
end do
return
end 
subroutine rhouv(f,rho,u,v,cx,cy,n,m)
real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m),cx(0:8),cy(0:8)
integer :: i, j, k
real :: ssum, usum, vsum
! 有障碍物的一部分
do j=21,m
    do i=0,40
        ssum=0.0
            do k=0,8
            ssum=ssum+f(k,i,j)
            end do
            rho(i,j)=ssum    
            usum=0.0 
            vsum=0.0
            do k=0,8
            usum=usum+f(k,i,j)*cx(k)
            vsum=vsum+f(k,i,j)*cy(k)  
            !print*,usum! 为什么？？？？？？
            end do
            u(i,j)=usum/rho(i,j)
            v(i,j)=vsum/rho(i,j)
            !print*,u(i,j),v(i,j)
        end do
end do
! 这一部分错了
do j=0,m
    do i=41,n
        ssum=0.0
        usum=0.0 
        vsum=0.0
        do k=0,8
            ssum=ssum+f(k,i,j)
        end do
        rho(i,j)=ssum  
        do k=0,8
            usum=usum+f(k,i,j)*cx(k)
            vsum=vsum+f(k,i,j)*cy(k)
        end do
        u(i,j)=usum/rho(i,j)
        v(i,j)=vsum/rho(i,j)
    end do
end do
do j=1,m
    v(n,j)=0.0
end do
do j=0,20
    do i=0,40
        u(i,j)=0.0
        v(i,j)=0.0
    end do
end do
return
    end 
subroutine result(u,v,rho,uo,n,m)  
! 流函数这一部分我太会
   real u(0:n,0:m),v(0:n,0:m),rho(0:n,0:m),strf(0:n,0:m)
    integer :: i, j
    real :: rhoav, rhom
    open(5, file='streamf')
    strf(0,0)=0.
    do i=1,n
        rhoav=0.5*(rho(i-1,0)+rho(i,0))
        strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
    end do
    do i=0,40
        do j=21,m
            rhom=0.5*(rho(i,j)+rho(i,j-1))
            strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
        end do
    end do
    do i=41,n
        do j=1,m
            rhom=0.5*(rho(i,j)+rho(i,j-1))
            strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
        end do
    end do
    ! 写入数据到输出文件
    write(2,*)"VARIABLES =X, Y, U, V, S"
    write(2,*)"ZONE ","I=",n+1,"J=",m+1,",","F=BLOCK"
    do j=0,m
        do i=0,n
            write(2,*)i,j,u(i,j),v(i,j),strf(i,j)
        end do
    end do
    do j=0,m
        write(3,*)j/float(m),u(n/2,j)/uo,u(n/4,j)/uo,u(3*n/4,j)/uo
    end do
    do i=0,n
        write(4,*) i/float(n),v(i,m/2)/uo
    end do
end subroutine result

do j=21,m
    do i=0,40
        ssum=0.0
            do k=0,8
            ssum=ssum+f(k,i,j)
            end do
            rho(i,j)=ssum 
            usum=0.0 
            vsum=0.0
            do k=0,8
            usum=usum+f(k,i,j)*cx(k)
            vsum=vsum+f(k,i,j)*cy(k)  
            !print*,usum! 为什么？？？？？？
            end do
            !print*,usum !这个可以正常输出数据
            u(i,j)=usum/rho(i,j)
            v(i,j)=vsum/rho(i,j)
            !print*,u(i,j),v(i,j)
        end do
end do
! 这一部分错了
do j=0,m
    do i=41,n
        ssum=0.0
        do k=0,8
            ssum=ssum+f(k,i,j)
        end do   
        rho(i,j)=ssum  
        usum=0.0 
        vsum=0.0
        do k=0,8
            usum=usum+f(k,i,j)*cx(k)
            vsum=vsum+f(k,i,j)*cy(k)    
            !print*,usum
        end do
        !print*,usum
        u(i,j)=usum/rho(i,j)
        v(i,j)=vsum/rho(i,j) 
    end do
end do