Fortran参数传递问题

a1100200420 · 发表于 2017-5-4 22:11:13

参数传递时，个别元素变成了NaN，恳请各位解惑。下面是我的主程序

[Fortran] 纯文本查看 复制代码

01program main !顶盖驱动boltzmann 程序
02    implicit none
03    real*8 dx,dy,Lx,Ly,dt,c,cs
04    integer Nx,Ny
05    real*8 Re,niu,U0,rho0
06    real*8,allocatable:: rho(:,:)
07    real*8,allocatable:: w(:),e(:,:),un(:,:,:),un_1(:,:,:)
08    real*8 tau
09    integer Q
10    real*8,allocatable::Fn(:,:,:),Fn_1(:,:,:)
11    integer i,j,k,n,n_max
12    real*8 feq,feq_n,du,u,err_e,err
13 
14    !初始化111111111111111111111111111111111111111111111111111111111111111111111111111111111111
15    n_max=20000 !迭代次数
16    err_e=10**(-6) !终止误差
17    dx=1.0 !格子的流向长度
18    dy=1.0 !格子的垂向长度
19    Nx=256 !流向格子数量
20    Ny=256
21    Lx=dx*Nx !流向长度
22    Ly=dy*Ny
23    dt=1.0 !时间步长
24    c=dx/dt !格子速度
25    cs=c/sqrt(3.0) !格子声速
26    Re=1000.0 !雷诺数
27    U0=0.1 !顶盖速度
28    rho0=1.0 !初始密度
29    niu=U0*Lx/Re !运动粘性系数
30    tau=3.0*niu+0.5 !松弛时间
31    Q=9 !速度离散数量
32    allocate(rho(Nx,Ny)) !各节点的密度
33    allocate(w(Q)) !权系数
34    allocate(e(Q,2)) !离散速度
35    allocate(un(Nx,Ny,2)) !各节点迭代后的速度
36    allocate(un_1(Nx,Ny,2)) !各节点迭代前的速度
37    allocate(Fn(Nx,Ny,Q)) !迭代后的分布函数
38    allocate(Fn_1(Nx,Ny,Q)) !迭代前的分布函数
39    rho(:,:)=rho0
40    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/)
41    !open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
42    !write(10,*) w !正常情况下应当屏蔽*******************************************************************
43    e(:,1)=(/ 0,1,0,-1,0,1,-1,-1,1 /)
44    e(:,2)=(/ 0,0,1,0,-1,1,1,-1,-1/)
45    !open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
46    !write(10,*) (e(i,:),i=1,Q,1) !正常情况下应当屏蔽***************************************************************
47    un(:,:,:)=0 !各节点速度
48    un(:,Ny,1)=U0
49    !open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
50    !write(10,*) ((un(i,j,1),i=1,Nx,1),j=1,Ny,1) !正常情况下应当屏蔽****************************************
51    rho(:,:)=0 !各节点密度
52    do i=1,Nx,1
53        do j=1,Ny,1
54            do k=1,Q,1
55                call sub_equi(rho(i,j),w(k),e(k,:),un(i,j,:),cs,feq)
56                Fn_1(i,j,k)=feq !分布函数初始化为平衡态
57            enddo
58        enddo       
59    enddo
60end
61 
62下面是我的子程序
63subroutine sub_equi(rho,w,e,u,cs,feq)
64    !该程序用于求解Maxwell 分布函数
65    !rho 密度
66    !w 权系数
67    !e 离散速度矢量
68    !u 格子宏观速度矢量
69    !cs 格子声速
70    !feq Maxwell 分布函数
71 
72    implicit none
73    real*8 rho,w,cs,feq,eu,uu
74    real*8 e(2),u(2)
75    integer i
76 
77    open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
78    write(10,*) u !正常情况下应当屏蔽***************************************************************
79 
80    eu=e(1)*u(1)+e(2)*u(2) !向量e、u的点乘结果
81    uu=u(1)**2+u(2)**2 !向量u、u的向量点乘
82 
83    feq=rho*w*(1+eu/cs**2+(eu/cs**2)**2/2-uu/(2*cs**2))
84  end