fortran程序运行无结果，且界面闪退
自己编写的交错网格有限差分的fortran程序，没有语言错误但是运行无结果且界面闪退，请求各位帮忙解决，急，谢谢！
贴代码，别截图
大概率是代码没问题，运行结束后自动关闭，看看这个解释
FAQ之 Console程序一闪而过
http://fcode.cn/guide-32-1.html
PS：贴代码没有fortran选项，可以当c代码贴出来
vs操作可以看视频https://www.bilibili.com/video/BV1oh411o7AT?p=2
program main
real dx,dz !网格间距
integer pml,modx,modz,n !吸收层厚度，模型网格大小（包含PML层），差分阶数
integer mx,mz !计算模型网格大小（包含PML层）
integer sx,sz,ns,ds,rd,off,dj !震源坐标,ns炮数，ds炮间隔,rd接收道数，最小偏移距,道距
integer tt,nt,out_t !采样时长,采样数
real mf,dt !子波主频,间隔（ms）
integer,allocatable::nsx(:),nsz(:)
real,allocatable::v(:,:),st(:,:)
real,allocatable::jl(:,:,:),jlw(:,:),out_w(:,:)
real start, finish
character*40 outfile,input_filename,outfile2
integer::i,j,k
!输入参数
input_filename='cmd.txt'
call input_cmd(input_filename,dx,dz,pml,n,modz,modx,tt,dt,ns,ds,rd,off,dj,sz,sx,mf,out_t)
!差分网格大小
mx=n*2+modx
mz=n*2+modz
!采样点数
nt=int(tt/dt)
allocate(nsx(ns),nsz(ns))
do i=1,ns
nsx(i)=sx+(i-1)*ds
nsz(i)=sz+n
end do
allocate(v(mz,mx))
allocate(jl(ns,nt,rd),jlw(nt,rd),out_w(mz,mx)) !jl(炮数，采样数，接收道数)
allocate(st(mz,mx))
!输入速度模型并标记地层
v=0
call vr(v,st,mz,mx,n)
do i=1,ns
call calculate_jcwg(i,pml,mx,mz,n,nsx(i),nsz(i),rd,off,dj,dt,nt,mf,dx,dz,v,st,j1w,out_w,out_t)
forall(j=1:nt,k=1:rd)
jl(i,j,k)=jlw(j,k)
end forall
end do
outfile='jlw.sgy'
outfile2='w.sgy'
!sgy格式输出数据
call out_sgy(nt,rd,ns,jl,outfile,dt)
call out_sgy(mz,mx,1,out_w,outfile2,dz)
end program main
subroutine input_cmd(input_filename,dx,dz,pml,n,modz,modx,tt,dt,ns,ds,rd,off,dj,sz,sx,mf,out_t)
character*40 input_filename
real dx,dz
integer pml,modx,modz,n
integer sx,sz,ns,ds,rd,off,dj
integer tt,nt,out_t
real mf,dt
input_filename = 'cmd.txt'
open(111,file=input_filename)
read(111,*)dx,dz
read(111,*)pml
read(111,*)n
read(111,*)modz,modx
read(111,*)tt,dt
read(111,*)ns,ds
read(111,*)rd,off,dj
read(111,*)sz,sx
read(111,*)mf
read(111,*)out_t
out_t=out_t/0.2
close(111)
end subroutine
subroutine calculate_jcwg(ins,pml,mx,mz,n,sx,sz,rd,off,dj,dt,nt,mf,dx,dz,v,st,j1w,out_w,out_t)
real dx,dz !节点长度
integer pml,n !吸收层厚度，空间阶数n*2
integer mx,mz !立方单元数
integer sx,sz !震源节点号
integer nt !总步数
integer rd,off,dj,out_t
real mf,dt !子波主频，采样间隔
real wavelet(nt) !地震子波值
real cn(n)
real vmax
integer ins,m
real v(mz,mx),st(mz,mx),dtdx
real U(mz,mx),W(mz,mx)
real p(mz,mx)
real decz(mz),decx(mx)
real pzw,pxu
real ppzw(mz,mx),ppxu(mz,mx)
real ppxp(mz,mx),ppzp(mz,mx)
real pxp,pzp
character*20 outfile,fft
character*4 filenumb
real d1,d2,a
real jlw(nt,rd),out_w(mz,mx) !切片输出
U=0.
W=0.
p=0.
wavelet=0.
ppzw=0
ppxu=0
ppxp=0
ppzp=0
call calcn(n,cn) !计算空间差分系数
call calwave(nt,wavelet,mf,dt) !计算地震子波
open(11,file='wavelet.dat')
do i=1,nt
write(11,*)wavelet(i)
end do
close(11)
call calpmldec(decz,decx,mx,mz,dx,n,pml,v) !计算pml衰减因子
if(maxval(v)*dt/1000/dx>sqrt(2.0)/2)then
write(*,*)'模型不满足稳定性条件'
pause
end if
a=sqrt(-1.0*log(0.01)/((10.0*dx)**2.0)); !空间衰减因子，其值越大则子波沿空间衰减越迅速
dtdx=dt/1000/dx
!开始计算
do it=1,nt
forall(i=n+1:mz-n,j=n+1:mx-n)
p(i,j) = p(i,j) + wavelet(it) * exp(-a**2 * (((i-sz)*dz)**2 + ((j-sx)*dx)**2))
end forall
do m=1,mx
p(n+1,m) = 0 ! 自由表面应力赋0
end do
do m=1,n
do k=1,mx
p(m,k) = -p(2*(n+1)-m,k) ! 应力关于自由表面上下反对称
end do
end do
do j=n+1,mx-n
do i=n+1,mz-n
if(st(i,j)==1)then !进行声波波动方程
!普通交错网格
pxp=0
do m=1,n
pxp=pxp+cn(m)*(p(i,j+m)-p(i,j-m+1)) !对x方向压力进行中心差分
end do
pzp=0
do m=1,n
pzp=pzp+cn(m)*(p(i+m,j)-p(i-m+1,j)) !对z方向压力进行中心差分
end do
ppxp(i,j)=ppxp(i,j)+(pxp-ppxp(i,j))*decx(j)*dt/1000
ppzp(i,j)=ppzp(i,j)+(pzp-ppzp(i,j))*decz(i)*dt/1000
u(i,j)=u(i,j)+(pxp-ppxp(i,j))*dtdx
w(i,j)=w(i,j)+(pzp-ppzp(i,j))*dtdx
!普通交错网格
!u
pxu=0
do m=1,n
pxu=pxu+cn(m)*(u(i,j+m-1)-u(i,j-m)) !对x方向速度进行中心差分
end do
!w
pzw=0
do m=1,n
pzw=pzw+cn(m)*(w(i+m-1,j)-w(i-m,j)) !对z方向速度进行中心差分
end do
ppxu(i,j)=ppxu(i,j)+(pxu-ppxu(i,j))*decx(j)*dt/1000
ppzw(i,j)=ppzw(i,j)+(pzw-ppzw(i,j))*decz(i)*dt/1000
p(i,j)=p(i,j)+(v(i,j)**2*(pxu-ppxu(i,j))+v(i,j)**2*(pzw-ppzw(i,j)))*dtdx !根据速度反演压力
else !进行弹性波波动方程
!普通交错网格
!p
pxp=0
do m=1,n
pxp=pxp+cn(m)*(p(i,j+m-1)-p(i,j-m)) !对x方向应力进行中心差分
end do
pzp=0
do m=1,n
pzp=pzp+cn(m)*(p(i+m,j)-p(i-m+1,j)) !对z方向应力进行中心差分
end do
ppxp(i,j)=ppxp(i,j)+(pxp-ppxp(i,j))*decx(j)*dt/1000
ppzp(i,j)=ppzp(i,j)+(pzp-ppzp(i,j))*decz(i)*dt/1000
u(i,j)=u(i,j)+(pxp-ppxp(i,j))*dtdx
w(i,j)=w(i,j)+(pzp-ppzp(i,j))*dtdx
!求取速度分量结束
!普通交错网格
!u
pxu=0
do m=1,n
pxu=pxu+cn(m)*(u(i,j+m)-u(i,j-m+1))
end do
!w
pzw=0
do m=1,n
pzw=pzw+cn(m)*(w(i+m-1,j)-w(i-m,j))
end do
ppzw(i,j)=ppzw(i,j)+(pzw-ppzw(i,j))*decz(i)*dt/1000
ppxu(i,j)=ppxu(i,j)+(pxu-ppxu(i,j))*decx(j)*dt/1000
p(i,j)=p(i,j)+(v(i,j)*v(i,j)*(pxu-ppxu(i,j))+v(i,j)*v(i,j)*(pzw-ppzw(i,j)))*dtdx
!求取应力分量结束
end if
end do
end do
do i=1,rd
!jlw(it,i)=w(sz,sx+(i-1)*dj-off)
jlw(it,i)=w(25,sx+(i-1)*dj-off) !从某一平面x方向输出
!起伏海底输出
!m=(i-1)*6+35
!if(m<335) then
! jlw(it,i)=w((1265-m)/6,m)
!else
! m=(i-1)*6+35
! jlw(it,i)=w((m+595)/6,m)
!end if
end do
if(it==out_t) then
out_w=w
!!输出某一偏移距的波形图
!open(333,file='250道.txt')
!do i=1,mz
!write(333,*) w(i,250)
!end do
!close(333)
else
end if
end do
open(555,file='20米深-50道.txt')
do i=1,nt
write(555,*) jlw(i,50) !输出某一道检波器单道记录
end do
close(555)
end subroutine
subroutine calcn(n,cn)
integer n,l
real cn(n)
real A(n,n),B(n)
do jc=1,n
do ic=1,n
A(ic,jc)=(2*jc-1)**(2*ic-1)
end do
end do
B=0.
B(1)=1.0
call AGAUS(A,B,n,cn,l)
if(l==0) then
write(*,*)'矩阵奇异，计算cn错误'
stop
end if
end subroutine
subroutine calwave(wl,wavelet,mf,dt) !计算地震子波
integer wl
real mf,pi,dt
real wavelet(wl)
real delay
pi=3.141592653
delay=1.5*sqrt(6.0)/(pi*mf)
do iw=1,wl
wavelet(iw)=(1-2.0*pi**2*mf**2*(1.0*(iw*dt/1000-delay))**2)* &
exp(-pi**2*mf**2*(1.0*(iw*dt/1000-delay))**2)
end do
end subroutine
subroutine out_sgy(m,n,ns,a,filename,dt)
integer m,n,ns
real dt
real a(ns,m,n)
character*20 filename
integer*2 seg_y_head(120)
integer*4 seg_y_head_4(60)
equivalence(seg_y_head,seg_y_head_4)
integer rec_long,i,j
seg_y_head(58)=m !采样点数
seg_y_head(59)=int(dt*1000) !采样间隔
rec_long=240+m*4 !计算道长
open(unit=13,file=filename,access='direct',status='replace',form='binary',recl=rec_long)
do is=1,ns
do i=1,n,1 !读记录到矩阵形式
write(13,rec=i+(is-1)*n)seg_y_head,a(is,:,i)
end do
end do
close(13)
end subroutine
subroutine calpmldec(decz,decx,mx,mz,dx,n,pml,v) !计算pml衰减系数
integer mx,mz,n,pml,x
real decz(mz),decx(mx),dx
real v(mz,mx)
real r,pi,c
pi=3.141592653
r=0.001
decz=0.
decx=0.
!do iz=n+1,pml+n+1
! x=n+pml+1-iz
! decz(iz)=v*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
!end do
do iz=mz-pml-n,mz-n
x=iz-(mz-pml-n)
c=v(iz,n+1)!选择z方向的速度作为衰减系数的计算
decz(iz)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jx=n+1,pml+n+1
x=n+pml+1-jx
c=v(n+1,jx)!选择x方向的速度作为衰减系数的计算
decx(jx)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jx=mx-pml-n,mx-n
x=jx-(mx-pml-n)
c=v(n+1,jx)!选择x方向的速度作为衰减系数的计算
decx(jx)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
end subroutine
SUBROUTINE AGAUS(A,B,N,X,L) !全选主元高斯消去法
real A(N,N),X(N),B(N),JS(N)
!DIMENSION A(N,N),X(N),B(N),JS(N)
!DOUBLE PRECISION A,B,X,T
L=1
DO 50 K=1,N-1
D=0.0
DO 210 I=K,N
DO 210 J=K,N
IF (ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
JS(K)=J
IS=I
END IF
210 CONTINUE
IF (D+1.0.EQ.1.0) THEN
L=0
ELSE
IF (JS(K).NE.K) THEN
DO 220 I=1,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
220 CONTINUE
END IF
IF (IS.NE.K) THEN
DO 230 J=K,N
T=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=T
230 CONTINUE
T=B(K)
B(K)=B(IS)
B(IS)=T
END IF
END IF
IF (L.EQ.0) THEN
WRITE(*,100)
RETURN
END IF
DO 10 J=K+1,N
A(K,J)=A(K,J)/A(K,K)
10 CONTINUE
B(K)=B(K)/A(K,K)
DO 30 I=K+1,N
DO 20 J=K+1,N
A(I,J)=A(I,J)-A(I,K)*A(K,J)
20 CONTINUE
B(I)=B(I)-A(I,K)*B(K)
30 CONTINUE
50 CONTINUE
IF (ABS(A(N,N))+1.0.EQ.1.0) THEN
L=0
WRITE(*,100)
RETURN
END IF
X(N)=B(N)/A(N,N)
DO 70 I=N-1,1,-1
T=0.0
DO 60 J=I+1,N
T=T+A(I,J)*X(J)
60 CONTINUE
X(I)=B(I)-T
70 CONTINUE
100 FORMAT(1X,' FAIL ')
JS(N)=N
DO 150 K=N,1,-1
IF (JS(K).NE.K) THEN
T=X(K)
X(K)=X(JS(K))
X(JS(K))=T
END IF
150 CONTINUE
RETURN
END SUBROUTINE
subroutine vr(v,st,mz,mx,n)
integer mz,mx,i,n,j
real v(mz,mx),st(mz,mx)
! 设定地层模型
do j=1,mx
do i=1,mz
if (i <= 200) then
v(i,j) = 1500
st(i,j) = 1
elseif (i > 200 .and. i <= 300) then
if (j >= 500) then
v(i,j) = 1500
st(i,j) = 1
elseif(j<=400) then
v(i,j)=2000
st(i,j)=2
elseif(j > (400+100.*(i-200)/100.) .and. j <= 500) then
v(i,j)=1500
st(i,j)=1
else
v(i,j)=2000
st(i,j)=2
end if
elseif(i > 300.and.i <= 400) then
v(i,j)=2000
st(i,j)=2
else
v(i,j)=2500
st(i,j)=2
end if
end do
end do
return
end subroutine
```贴出来的代码
```c