[数值问题] Fortran参数传递问题

program main !顶盖驱动boltzmann 程序
    implicit none
    real*8 dx,dy,Lx,Ly,dt,c,cs
    integer Nx,Ny
    real*8 Re,niu,U0,rho0
    real*8,allocatable:: rho(:,:)
    real*8,allocatable:: w(:),e(:,:),un(:,:,:),un_1(:,:,:)
    real*8 tau
    integer Q
    real*8,allocatable::Fn(:,:,:),Fn_1(:,:,:)
    integer i,j,k,n,n_max
    real*8 feq,feq_n,du,u,err_e,err

    !初始化111111111111111111111111111111111111111111111111111111111111111111111111111111111111
    n_max=20000 !迭代次数
    err_e=10**(-6) !终止误差
    dx=1.0 !格子的流向长度
    dy=1.0 !格子的垂向长度
    Nx=256 !流向格子数量
    Ny=256
    Lx=dx*Nx !流向长度
    Ly=dy*Ny
    dt=1.0 !时间步长
    c=dx/dt !格子速度
    cs=c/sqrt(3.0) !格子声速
    Re=1000.0 !雷诺数
    U0=0.1 !顶盖速度
    rho0=1.0 !初始密度
    niu=U0*Lx/Re !运动粘性系数
    tau=3.0*niu+0.5 !松弛时间
    Q=9 !速度离散数量
    allocate(rho(Nx,Ny)) !各节点的密度
    allocate(w(Q)) !权系数
    allocate(e(Q,2)) !离散速度
    allocate(un(Nx,Ny,2)) !各节点迭代后的速度
    allocate(un_1(Nx,Ny,2)) !各节点迭代前的速度
    allocate(Fn(Nx,Ny,Q)) !迭代后的分布函数
    allocate(Fn_1(Nx,Ny,Q)) !迭代前的分布函数
    rho(:,:)=rho0
    w(:)=(/4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36/)
    !open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
    !write(10,*) w !正常情况下应当屏蔽*******************************************************************
    e(:,1)=(/ 0,1,0,-1,0,1,-1,-1,1 /)
    e(:,2)=(/ 0,0,1,0,-1,1,1,-1,-1/)
    !open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
    !write(10,*) (e(i,:),i=1,Q,1) !正常情况下应当屏蔽***************************************************************
    un(:,:,:)=0 !各节点速度
    un(:,Ny,1)=U0 
    !open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
    !write(10,*) ((un(i,j,1),i=1,Nx,1),j=1,Ny,1) !正常情况下应当屏蔽****************************************
    rho(:,:)=0 !各节点密度
    do i=1,Nx,1
        do j=1,Ny,1
            do k=1,Q,1
                call sub_equi(rho(i,j),w(k),e(k,:),un(i,j,:),cs,feq)
                Fn_1(i,j,k)=feq !分布函数初始化为平衡态
            enddo
        enddo        
    enddo
end

下面是我的子程序
subroutine sub_equi(rho,w,e,u,cs,feq)
    !该程序用于求解Maxwell 分布函数
    !rho 密度
    !w 权系数
    !e 离散速度矢量
    !u 格子宏观速度矢量
    !cs 格子声速
    !feq Maxwell 分布函数

    implicit none
    real*8 rho,w,cs,feq,eu,uu
    real*8 e(2),u(2)
    integer i

    open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
    write(10,*) u !正常情况下应当屏蔽***************************************************************

    eu=e(1)*u(1)+e(2)*u(2) !向量e、u的点乘结果
    uu=u(1)**2+u(2)**2 !向量u、u的向量点乘

    feq=rho*w*(1+eu/cs**2+(eu/cs**2)**2/2-uu/(2*cs**2))
  end