[求助] 编译没错运行输出数据文件却是空的？

module global
implicit none
real EK,BX,BY,PTR,AK
common EK,BX,BY,PTR,AK
dimension AK(0:65,0:65)
end module

program main
use global
implicit none
real X,Y,H,P,W,DX
dimension X(65),Y(65),H(65,65),P(65,65),W(65,65)
real*8:: RX=0.3,RY=0.5,KA,KB,AA,BB,D,C,PH,RAD,W1,H0,EA,EB !,HMIN
integer :: N=65,MM,I,J
real :: E1=2.21E11,X0=-2.0,XE=2.0,W0=200.0,PAI=3.14159265,PAI1=0.2026423,Y0=-2.0
open(4,file="OUT1.DAT",status="unknown")
open(8,file="FILM1.DAT",status="unknown")
open(10,file="PRESS1.DAT",status="unknown")
EK=RX/RY
AA=0.5*(1./RX+1./RY)
BB=0.5*abs(1./RX-1./RY)
call HERTZELLIPTIC(RX,RY,KA,KB)
EA=KA*(1.5*W0/AA/E1)**(1.0/3.0)
EB=KB*(1.5*W0/AA/E1)**(1.0/3.0)
PH=1.5*W0/(EA*EB*PAI)
write(4,*)N,X0,XE,W0,PH,E1,RX,RY
MM=N-1
BX=EB
BY=EA
if(RX.GT.RY)then
  BX=EA
  BY=EB
end if
PTR=3.0*W0*(RX/BY)*(RX/BX)**2/(PAI**2)/(E1*RX**2)
call SUBAK(MM)
DX=(XE-X0)/(N-1.0)
do I=1,N
  X(I)=X0+(I-1)*DX
  Y(I)=Y0+(I-1)*DX
end do
do I=1 ,N
  D=1.0-X(I)*X(I)
  do J=1,N
    C=D-(BX/BY)**2*Y(J)*Y(J)
    if(C.LE.0.0)then
      P(I,J)=0.0
    end if
    if(C.GT.0.0)then
      P(I,J)=sqrt(C)
    end if
  end do
end do
call VI(N,DX,P,W)
do I=1,N
  do J=1,N
    RAD=X(I)*X(I)+EK*Y(J)*Y(J)
    W1=0.5*RAD
    H0=W1+W(I,J)
    !if(H0.LT.HMIN)then
      !HMIN=H0
    !end if
    H(I,J)=H0
  end do
end do
call OUTPUT(N,X,Y,H,P)
stop
end

subroutine VI(N,DX,P,V)
use global
implicit none
real DX,P,V
real*8 H0
integer N,I,J,K,L,IK,JL
dimension P(N,N),V(N,N)
do I=1,N
  do J=1,N
    H0=0.0
    do K=1,N
      IK=iabs(I-K)
      do L=1,N
        JL=iabs(J-L)
        H0=H0+AK(IK,JL)*P(K,L)
      end do
    end do
    V(I,J)=H0*DX*PTR
  end do
end do
return
end

subroutine SUBAK(MM)
use global
implicit none
real S,X,Y,XP,XM,YP,YM,A1,A2,A3,A4
integer I,J,MM
S(X,Y)=X+sqrt(X**2+Y**2)
do I=0,MM
  XP=I+0.5
  XM=I-0.5
  do J=0,I
    YP=J+0.5
    YM=J-0.5
    A1=S(YP,XP)/S(YM,XP)
    A2=S(XM,YM)/S(XP,YM)
    A3=S(YM,XM)/S(YP,XM)
    A4=S(XP,YP)/S(XM,YP)
    AK(I,J)=XP*alog(A1)+YM*alog(A2)+XM*alog(A3)+YP*alog(A4)
    AK(J,I)=AK(I,J)
  end do
end do
return
end

subroutine OUTPUT(N,X,Y,H,P)
implicit none
real X,Y,H,P,A
integer N,I,J
dimension X(N),Y(N),H(N,N),P(N,N)
A=0.0
write(8,110)A,(Y(I),I=1,N)
do I=1,N
  write(8,110)X(I),(H(I,J),J=1,N)
end do
write(10,110)A,(Y(I),I=1,N)
do I=1,N
  write(10,110)X(I),(P(I,J),J=1,N)
end do
110 format(66(E12.6,1X))
return
end

subroutine HERTZELLIPTIC(RX,RY,KA,KB)
implicit none
real*8,external :: EE,KE
real*8 RX,RY,BPA,BMA,CTH,THT,E1,KA,KB
real*8,parameter :: PAI=3.1415926
BPA=0.5*(1.0/RX+1.0/RY)
BMA=0.5*abs(1.0/RX-1.0/RY)
CTH=BMA/BPA
THT=ACOS(CTH)*180.0/PAI
call CACUE(CTH,E1)
KA=(2.0*EE(E1)/(PAI*(1-E1**2)))**(1/3.0)
KB=KA*(1.0-E1**2)**(1/2.0)
return
end

subroutine CACUE(CTH,E1)
implicit none
real*8,external :: EE,KE
real*8,external :: FAB
integer FLG,I
real*8 PAI,CTH,E1,E11,E12,DX,A,B,A1,A2,A3,A4,A5,T1,T2,T3,T4,T5,ER0
data PAI,DX,FLG,I,T1,T5,ER0/3.1415926,0.0001,1,1,1.E-30,1.,1.E-12/
if(CTH.LT.1.E-6)then
  write(*,*)"!NOTE:COS(THETA) IS TOO SAMLL TO CALCULATE,E1 IS SETTED TO 0."
  E1=0.0
end if
if(CTH.GT.0.9999999999)then
   E1=1.0
end if
!在e1和e5之间找e3，使f(e1)*f(e3)<0，并且f(e3)*f(e5)<0
A1=FAB(T1,CTH)
A5=FAB(T5,CTH)
do while(FLG.EQ.1)
  T3=T1+I*DX
  A3=FAB(T3,CTH)
  I=I+1
  if((A1*A3.LT.0.).AND.(A3*A5.LT.0.))then
    FLG=0
  end if
end do
!二分法找e1和e3之间的e2，并将e4赋值给e11
do while((T3-T1).GT.ER0)
  T2=(T1+T3)/2.0
  A2=FAB(T2,CTH)
  if(A2.GT.0.)T1=T2  
  if(A2.LT.0.)T3=T2
  if(A2.EQ.0.)then
    E11=T2
    exit
  end if
end do
E11=T2
!二分法找e3和e5之间的e4，并将e4赋值给e12
do while((T5-T3).GT.ER0)
  T4=(T3+T5)/2.0
  A4=FAB(T4,CTH)
  if(A4.GT.0.) T5=T4
  if(A4.LT.0.) T3=T4
  if(A4.EQ.0.)then
    E12=T2
    exit
  end if
end do
E12=T4

E1=E11
if(E11.LT.E12) E1=E12
return
end

!定义方程FAB
real*8 function FAB(E1,CTH)
implicit none
real*8 E1,CTH,T1,T2
real*8,external :: EE,KE
T1=EE(E1)
T2=KE(E1)
FAB=2*(1.0-E1**2)*(T1-T2)+(1.0-CTH)*E1**2*T1
return
end

!定义方程EE
real*8 function KE(E1)
implicit none
integer N,I,FLG
real*8 E1,H,T,T1,T2,S1,S2,P,Q
real*8,parameter :: PAI=3.1415926
!如果e趋近于1，K(e)趋向于无穷大，赋值为1.E10
if(E1.EQ.1)then
  KE=1.E10
end if
!如果e趋近于0，K(e)等于AI/2.0
if(E1.LT.1.E-20)then
  KE=PAI/2.0
end if
!当在0与1之间时，使用辛普森积分法求解积分
N=1
H=PAI/2.0
Q=sqrt(1.-E1*E1*sin(H)*sin(H))
if(Q.LT.1.E-35)then
  Q=1.E35
  Q=1./Q
end if
T1=0.5*H*(1+Q)
S1=T1
FLG=1
do while(FLG.EQ.1)
  P=0
  do I=0,N-1
    T=(I+0.5)*H
    Q=sqrt(1.-E1*E1*sin(T)*sin(T))
    if(Q.LT.1.E-35)then
      Q=1.E35
      Q=1./Q
    end if   
    P=P+Q
  end do
  T2=(T1+H*P)/2.
  S2=(4.*T2-T1)/3.
  if(abs(S2-S1).LT.abs(S2)*1.E-7)then
    FLG=0
  end if
  T1=T2
  S1=S2
  N=N+N
  H=0.5*H
end do
KE=S2
return
end

!定义方程EE
real*8 function EE(E1)
implicit none
integer N,I,FLG
real*8 E1,H,T,T1,T2,S1,S2,P,Q
real*8,parameter :: PAI=3.1415926
N=1
H=PAI/2.0
if(E1.EQ.1)then
  EE=1.0
end if
if(E1.LT.1.E-20)then
  EE=PAI/2.0
end if
Q=sqrt(1.-E1*E1*sin(H)*sin(H))
T1=.5*H*(1+Q)
S1=T1
FLG=1
do while(FLG.EQ.1)
  P=0
  do I=0,N-1
    T=(I+0.5)*H
    Q=sqrt(1.-E1*E1*sin(T)*sin(T))      
    P=P+Q
  end do
  T2=(T1+H*P)/2.
  S2=(4.*T2-T1)/3.
  if(abs(S2-S1).LT.abs(S2)*1.E-7)then
    FLG=0
  end if
  T1=T2
  S1=S2
  N=N+N
  H=0.5*H
end do
EE=S2
return
end