我正在 Fortran 中运行一个程序，该程序进行大量计算。基本上，它运行几个计算周期的次数与实验中的粒子数一样多。
我已经能够成功测试并获得 500 000、1 000 000 和 1 500 000 个粒子的输出数据，但是一旦我使用更多输入粒子（例如，2 000 000 或 5 000 000）运行程序，它就会崩溃接近结尾，在它正在读取的文件中显示“文件结尾”错误。读取提到的文件是我正在执行的计算的必要条件，并且程序在每个周期中都会执行此操作。
我关心的是程序如何能够准确地进行所有计算 1 500 000 次，而不是 2 000 000 次？在执行大量计算时，fortran 代码是否存在使程序不稳定的限制？
该程序适用于低于一定数量的任意数量的粒子这一事实表明一切都应该井井有条。文件在读取操作结束时倒带，依此类推。如何确定问题？
谢谢你。
更新：好的，我附上代码，以便人们可以查看它并希望有人能帮助我。这是第二个计算周期（第一个计算周期即使在计算大量粒子时也没有显示错误）。
NSTAT = 0
NSTAP = int(NSTATN/10)
NDSTAT = NSTAP ! Number STAT for visualisation
10 continue
NSTAT=NSTAT+1
x = 50.*sqrt(RNDM(IY))
y = 0.
z = 1.E-20
ct = 1. ! initial state photon
st = sqrt(1.-ct*ct)
cf = 1.
sf = 0.
ves = 1.
E = YNT2001()
Call FGAMMA() ! -> calculation of penetration photons
!** zahler
IF(NSTAT.GE.NDSTAT) GO TO 15
GO TO 11
15 PRINT*, NSTAT
IAWA = NSTAT*100/NSTATN
WRITE (TEXTW,103) IAWA
PRINT*, TEXTW ! print STAT
NDSTAT=NDSTAT+NSTAP
11 IF(NSTAT.GE.NSTATN) GO TO 20
GO TO 10
!** zahler
20 call RESULT ! Results of Calculation
这部分在这里选择随机启动能量并处理计算以及在输出窗口（zahler）中显示进度。让我们进一步追踪计算。法格玛：
SUBROUTINE FGAMMA()
COMMON /SOS/ E,X,Y,Z,CT,ST,CF,SF,W
COMMON /RMX/ Rmax,Rkmax
COMMON/OTC/ GOTC,WOTC
COMMON /ktt/ ktot,kcom,kpho,kpar
COMMON /RES/ VS,VSI,TD
COMMON /wpp/ Wp
COMMON /w51/ w511,IP
VS = 0.
VSI = 0.
ktot = 1;kcom = 1;kpho = 1;kpar = 1
CALL SIG_phots(E,Stots,SSKs,SSFs,SSPs)
if(ct > 0.) then
dl = (Rmax - z)/ct
else
dl = - z/ct
endif
sss = dl*Stots
WFZ = 1- exp(-sss) ! Wes wzaimodeystwiy
RN = RNDM(IY)
IF(RN <= 0.) RN=1.E-4
DLPR = - LOG(1.-RN*WFZ)/Stots ! probeg of photon
1 Z=Z+DLPR*CT
RR=DLPR*ST
X=X+RR*CF ! new coordianats of photon
Y=Y+RR*SF
100 continue
!*************** para ****************************
If(E > 1.022) then
Wp= SSPs/STOTs*W*WFZ ! Probabilaty of paar creation
IP = 1
GP = GOTC
GOTC = 0.06
CALL BLAN() !!!!!!!!!!!!!! WKLAD of effekt
GOTC = GP
endif
!************** compton **************************
W= SSKs/STOTs*W*WFZ ! Probabilaty of Compton creation
IP = 0
CALL BLPS() !!!!!!!!!!!!!!WKLAD Local Estimation
EP = E
CALL COMKL3(Ep,E)
IF(E.LT.GOTC) RETURN
!NEXT COLLISIONS
ktot = 1;kcom = 1;kpho = 1;kpar = 1
CALL SIG_phots(E,Stots,SSKs,SSFs,SSPs)
if(ct > 0.) then
dl = (Rmax - z)/ct
else
dl = - z/ct
endif
sss = dl*Stots
WFZ = 1- exp(-sss) ! Wes wzaimodeystwiy
RN = RNDM(IY)
IF(RN <= 0.) RN=1.E-4
DLPR = - LOG(1.-RN*WFZ)/Stots
Z=Z+DLPR*CT
RR=DLPR*ST
X=X+RR*CF
Y=Y+RR*SF
GO TO 100
END
正如我们所看到的，这里有更多的子程序。我也在提供它们。
SUBROUTINE BLAN()
COMMON /SOS/ E,X,Y,Z,CT,ST,CF,SF,VES
COMMON /DET/ XDET,YDET,ZDET
COMMON /RES/ VS,VSI,TD
COMMON /OTC/ GOTC,WOTC
COMMON /A511/ S511s
COMMON /RMX/ Rmax,Rkmax
COMMON /ktt/ ktot,kcom,kpho,kpar
COMMON /wpp/ Wp
COMMON /w51/ w511,IP
TW=E
XW=X
YW=Y
ZW=Z
COSTW=CT !Memory 1 photon
SINTW=ST
COSFW=CF
SINFW=SF
VESW=VES
ves = wp
XD=XDET-X
YD=YDET-Y
ZD=ZDET-Z
DW = SQRT(XD*XD+YD*YD+ZD*ZD)
DWW = 1./DW
AL2=(XDET-X)*DWW
AM2=(YDET-Y)*DWW
AN2=(ZDET-Z)*DWW
if(ZDET > 0.) then
DWZ = (Rmax-Z)/AN2
else
DWZ = -Z/AN2
endif
!+++++++++++++++++++++++++++++++++++++++++++ RANGE IN MATTER
CT = AN2
if(CT > 1. ) CT = 1.
if(CT < -1.) CT = -1.
ST = SQRT(1.-CT*CT)
IF(ST.LE.0.) THEN
CF = 1.
SF = 0.
ELSE
if(AN2 > 1. ) AN2 = 1.
if(AN2 < -1.) AN2 = -1.
CF = AL2/SQRT(1.-AN2*AN2)
if(CF < -1.) CF =-1.
if(CF > 1. ) CF = 1.
WSFK = SQRT(1.-CF*CF)
SF = AM2/SQRT(1.-AN2*AN2)
IF(SF > 0.)SF = WSFK
IF(SF < 0.)SF = -WSFK
ENDIF
! FIRST
RQ=DWZ*S511s
w511 = EXP(-RQ)*2.
!+++++++++++++++++++++++++++++++++++++++++++ RANGE IN MATTER
VS = EXP(-RQ)*0.15915*Wp
VSI=VS*0.511
TD=0.511
CALL ANALYS
!!+++++++++++++++++++++++++++++++++++++++++++ first 0.511
CT=1.-2.*RNDM(IY)
ST=sqrt(1.-CT*CT)
Cf=cfp(Sf)
CSW=CT
STW=ST
CFW=CF
SFW=SF
E=0.511
ves = wp
IP = 0.
CALL GAM_511()
!!+++++++++++++++++++++++++++++++++++++++++++ second 0.511
E=0.511
CT=-CSW
ST=-STW
CF=-CFW
SF=-SFW
X=XW
Y=YW
Z=ZW
ves = wp
IP = 0.
CALL GAM_511()
E=TW
X=XW
Y=YW
Z=ZW
CT=COSTW !Memory 1 photon
ST=SINTW
CF=COSFW
SF=SINFW
VES=VESW
kcom = 1
kpar = 1
ktot = 1
RETURN
END
!************************************************************************
SUBROUTINE BLPS()
COMMON /RES/ VS,VSI,TD
COMMON /SOS/ T,X,Y,Z,COST,SINT,COSF,SINF,VES
COMMON /DET/ XDET,YDET,ZDET
COMMON /OTC/ GOTC,WOTC
COMMON /ktt/ ktot,kcom,kpho,kpar
COMMON /UG1/ SINP,COSP,SFK,CFK
COMMON /w51/ w511,IP
COMMON /wpp/ Wp
COMMON /RMX/ Rmax,Rkmax
parameter (pi=3.14159,re=2.8179E-13,rek12=0.5*re*re,pire2=pi*re*re)
parameter (amc2=0.511,aimc2=1./amc2)
COSTW=COST
SINTW=SINT
COSFW=COSF
SINFW=SINF
VESW=VES
AL1=SINT*COSF
AM1=SINT*SINF
AN1=COST
DW=SQRT((XDET-X)**2+(YDET-Y)**2+(ZDET-Z)**2)
DWW = 1./DW
AL2=(XDET-X)*DWW
AM2=(YDET-Y)*DWW
AN2=(ZDET-Z)*DWW
COSK=AL1*AL2+AM1*AM2+AN1*AN2
AL=T/0.511
TD=T/(1.+AL*(1.-COSK))
IF(TD.lt.GOTC) GO TO 12
W=TD/T
e0 = T*aimc2
twoe01 = 2.*e0+1.
Wl =alog(twoe01)
compt1 = 2.*pire2*( (1.+e0)/(e0*e0)*( 2.*(1.+e0)/twoe01-&
Wl/e0 )+Wl/(2.*e0)-(1.+3.*e0)/(twoe01*twoe01))
dsdo1 = rek12*W*W*(W+1./W-1.+COSK*COSK) !Dsigma/Domega/sig_tot -> Compton
ot_sec = dsdo1/compt1
!+++++++++++++++++++++++++++++++++++++++++++ RANGE IN MATTER
COST = AN2
if(COST > 1. ) COST = 1.
if(COST < -1.) COST = -1.
SINT = SQRT(1.-COST*COST)
IF(SINT.LE.0.) THEN
COSF = 1.
SINF = 0.
ELSE
if(AN2 > 1. ) AN2 = 1.
if(AN2 < -1.) AN2 = -1.
COSF = AL2/SQRT(1.-AN2*AN2)
if(COSF < -1.) COSF =-1.
if(COSF > 1. ) COSF = 1.
WSFK = SQRT(1.-COSF*COSF)
SFK = AM2/SQRT(1.-AN2*AN2)
IF(SFK > 0.)SINF = WSFK
IF(SFK < 0.)SINF = -WSFK
ENDIF
ktot = 1;kcom = 0;kpho =0; kpar = 0
! FIRST
CALL SIG_phots(TD,SZs,TSK,TSF,TPP)
if(ZDET > 0.) then
DWZ = (Rmax-Z)/AN2
else
DWZ = -Z/AN2
endif
RQ=DWZ*SZs
!+++++++++++++++++++++++++++++++++++++++++++ RANGE IN MATTER
IF(RQ.GE.35.) then
VS = 0.
VSI = 0.
GO TO 12
endif
IP=0
if(IP==0) then
VS=VES*ot_sec*EXP(-RQ)
else
VS=wp*ot_sec*EXP(-RQ)
endif
IF(VS.LT.WOTC) then
VS = 0.
VSI = 0.
GO TO 12
endif
VSI=VS*TD
CALL ANALYS
12 COST=COSTW
SINT=SINTW
COSF=COSFW
SINF=SINFW
VES=VESW
kcom = 1
kpar = 1
ktot = 1
END
!*******************************************************
SUBROUTINE GAM_511()
COMMON /SOS/ E,X,Y,Z,CT,ST,CF,SF,W
COMMON /RMX/ Rmax,Rkmax
COMMON/OTC/ GOTC,WOTC
COMMON /ktt/ ktot,kcom,kpho,kpar
COMMON /w51/ w511,IP
ktot = 1;kcom = 1;kpho = 1;kpar = 0
CALL SIG_phots(E,Stots,SSKs,SSF,SSPs)
RN = RNDM(IY)
IF(RN <= 0.) RN=1.E-4
RNl = - ALOG(RN)
DLPR = RNl/Stots
Z=Z+DLPR*CT
IF(Z < 0.) RETURN
IF(Z > Rmax) RETURN
RR=DLPR*ST
X=X+RR*CF
Y=Y+RR*SF
100 continue
RN=RNDM(IY)
STET=1./STOTs
IF(RN.LT.SSKs*STET) GO TO 54 !COMPTON
!*************** photo ***************************
RETURN
!*************************************************
!************** compton **************************
54 continue
IP = 1.
CALL BLPS() !!!!!!!!!!!!!!WKLAD
EP =E
CALL COMKL3(Ep,E)
IF(E.LT.GOTC) RETURN
!NEXT COLLISIONS
! ktot = 1;kcom = 1;kpho = 1;kpar = 0
CALL SIG_phots(E,Stots,SSKs,SSFs,SSPs)
RN = RNDM(IY)
IF(RN.LE.0.) RN=1.E-4
RNl = - ALOG(RN)/Stots
Z=Z+DLPR*CT
if(z < 0.) return
IF(Z > Rmax) RETURN
RR=DLPR*ST
X=X+RR*CF
Y=Y+RR*SF
GO TO 100
END
Subroutine SIG_phots(E_ph,Stot,Scom,Spho,Spai)
!******************************************************
! ktot -> total if = 1
! kcom -> Compton if = 1
! kpho -> Photo if = 1
! kpar -> Pair if = 1
!******************************************************
COMMON /SBTgs/etabls(1000),SFLs(1000),SKKs(1000),SPPs(1000),SGTs(1000)
COMMON /SB/ dem,Emin,deb
COMMON /ktt/ ktot,kcom,kpho,kpar
IF(E_ph.LT.1.) THEN
i = int((E_ph-Emin)*dem)+1
!print*,E_ph,Emin,dem
!read(*,*)
DX=E_ph-etabls(I)
AA=DX/(etabls(I+1)-etabls(I))
Spai = 0.
IF(ktot.eq.1) Stot = SGTs(I)+(SGTs(I+1)-SGTs(I))*AA
IF(kcom.eq.1) Scom = SKKs(I)+(SKKs(I+1)-SKKs(I))*AA
IF(kpho.eq.1) Spho = SFLs(I)+(SFLs(I+1)-SFLs(I))*AA
return
ENDIF
i = int((E_ph-1.)*deb)+200
DX=E_ph-etabls(I)
AA=DX/(etabls(I+1)-etabls(I))
IF(ktot.eq.1) Stot = SGTs(I)+(SGTs(I+1)-SGTs(I))*AA
IF(kcom.eq.1) Scom = SKKs(I)+(SKKs(I+1)-SKKs(I))*AA
IF(kpho.eq.1) Spho = SFLs(I)+(SFLs(I+1)-SFLs(I))*AA
Spai = 0.
IF(E_ph.gT.1.2) THEN
IF(kpar.eq.1) Spai = SPPs(I)+(SPPs(I+1)-SPPs(I))*AA
ENDIF
return
END
SUBROUTINE COMKL3(T,TR)
COMMON /SOS/ E,X,Y,Z,CT,ST,CF,SF,WW
COMMON /UG1/ SINTP,COSTP,SINFP,COSFP
TAU = T*1.9569
TAUS = TAU
ab = 1./(1.+2.*TAU)
A = TAU*ab
AA = A+TAU
FA = 1.+2.*TAU+ab
FAK=FA+2.
GT=TAU-A
6 CSIS=A+RNDM(IY)*GT
R1=RNDM(IY)
ETA=RNDM(IY)*FAK
FIL=FA-2.*(CSIS-A)
IF(ETA.GT.FIL) THEN
CSI=AA-CSIS
ETA=FAK-ETA
ELSE
CSI=CSIS
ENDIF
AB=1./TAU-1./CSI
FKA=TAU/CSI+CSI/TAU+AB*(2.+AB)
IF(ETA.GT.FKA) THEN
GO TO 6
ELSE
TR=CSI*0.511
TAU=CSI
ENDIF
COSTP = 1.-1./TAU+1./TAUS
S = 1.-COSTP*COSTP
IF(S.lt.0.) S = 0.
COSFP = CFP(SINFP)
SINTP=SQRT(S)
CALL UGOLKY(CT,ST,CF,SF,A1,A2,A3,A4,1)
END
SUBROUTINE COMKL(T)
COMMON/UG1/ STDOP,CTDOP,SFDOP,CFDOP
COMMON /SOS/ E,X,Y,Z,CT,ST,CF,SF,WW
TAU=T*1.9569E0
TS=T
TAUS=TAU
1 CONTINUE
ALFR=TAU*(1.+2.D0*RNDM(iy)*TAU)/(1.+2.*TAU)
AN=1.0/TAU
AR=1.0/ALFR
FF=AN*ALFR+TAU*AR+(AN-AR)*(2.0+AN-AR)
PRW=RNDM(iy)*(1.0+2.0*TAU+1.0/(1.0+2.0*TAU))
IF(FF.GT.PRW) THEN
T=ALFR*0.511
TAU = ALFR
ELSE
GO TO 1
ENDIF
COSTGP=1.0-(TAUS-TAU)/(TAU*TAUS)
S=1.0-COSTGP*COSTGP
if(s.lt.0.0)S=0.0
CALL azim(CFDOP,SFDOP)
CTDOP= COSTGP
STDOP= SQRT(S)
L=1
CALL UGOLKY(CT,ST,CF,SF,A1,A2,A3,A4,L)
END
SUBROUTINE UGOLKY(CTE,STE,CFE,SFE,COST,SINT,COSF,SINF,L)
COMMON/UG1/ ST,CT,SF,CF
A1=ST*CF
COST=CTE*CT-STE*A1
IF(COST.GT.1.)COST=1.
if(SQRT(1.-COST*COST)<0.) then
print*,COST,COST*COST
read(*,*)
endif
SINT=SQRT(1.-COST*COST)
IF(abs(SINT).GT.1.E-5) GO TO 2
1 SINF=0.
COSF=1.
GO TO 3
2 A2=STE*CT+CTE*A1
A1=ST*SF
COSF=(A2*CFE-A1*SFE)/SINT
SINF=(A2*SFE+A1*CFE)/SINT
3 GO TO (4,5),L
4 STE=SINT
CTE=COST
SFE=SINF
CFE=COSF
5 RETURN
END
那么……确定大量粒子计数的误差就足够了吗？