[子程序] call subroutine sub() 的问题

program main
 implicit none
 call subroutine sub()
 stop
 end program main

 subroutine sub()
 implicit none
 real::PI=3.1415926
 real::PM=1
 real::M=1
 real::PL=1
 real::PT=1
 real,parameter::X0=20
 real,parameter::Y0=20
 real,parameter::T0=20
 real::a(X0,Y0,T0)
 real::b(X0,Y0,T0)
 real::c(X0,Y0,T0)
 real::d(X0,Y0,T0)
 real::ans(X0,Y0,T0)
 real X
 real Y
 real T 
 real::x1,x2,dx,y1,y2,dy,t1,t2,dt
 integer::i,j,k,Nx,Ny,Nt
 x1=0.5
 x2=10
 y1=0.5
 y2=10
 t1=0.5
 t2=10
 dx=0.5
 dy=0.5
 dt=0.5
 Nx=(x2-x1)/dx
 Ny=(y2-y1)/dy
 Nt=(t2-t1)/dt
 open(unit=10,file="data1.txt")
 do i=1,Nx+1
   X=x1+(i-1)*dx
     do j=1,Ny+1
           y=y1+(j-1)*dy
       do k=1,Nt+1 
       T=t1+(k-1)*dt
       a(i,j,k)=PM*M/(4*(PI**1)*(PL**0.5)*(PT**0.5)*(T**1))
       b(i,j,k)=-1/(2*PL*T)
       c(i,j,k)=((X-T)/(2*PL*T))**2
       d(i,j,k)=((X-T)**2)/(4*PL*T)+(y**2/(4*PT*T))
       ans(i,j,k)=(PM*M/(4*(PI**1)*(PL**0.5)*(PT**0.5)*(T**2)))*((-1/(2*PL*T))+((((X-T)/(2*PL*T))**2)**2))*(exp(-(((X-T)**2)/(4*PL*T)+(y**2/(4*PT*T)))))
       if(ans(i,j,k)>-0.01.and.ans(i,j,k)<0.01) then
           write(10,"('(',f5.2,',',f5.2,','f5.2')=',(2X,F9.6),(2X,F9.6),(2X,F9.6),(2X,F9.6),(2X,F9.6))") X,Y,T,a(i,j,k),b(i,j,k),c(i,j,k),d(i,j,k),ans(i,j,k)
       end if
          end do
   end do
 end do
 return
 end