问题遇到的现象和发生背景
Fortran77代码在VScode中运行结果出现浮点数下沉，但是精度已经是双精度（real*16）
用代码块功能插入代码，请勿粘贴截图
c Appendix Fortran program for the collision model
c TDHF-motivated classical collision model, by G. Bertsch
c reconstructed by A.Bonasera, Texas A&M, August 2020.
c include imaginary time propagation, TAMU september 2020.
c include CC resonances from Tumino Nature paper. Only natural widths
integer stage,td
parameter (nmax=100)
common dndep(nmax),dndea(nmax),dndeg(nmax),dndebe(nmax)
c real*8 m,m12,l,li
c real*16 prob,sigmatot,sstar
real*16 m,m12,l,li
real*16 prob,sigmatot,sstar
dimension vrt(100),d(8)
equivalence (z1,d(1)),(z2,d(2)),(a1,d(3)),(a2,d(4)),(ecm,d(5)),
1 (ra,d(6)),(b,d(7))
c data m,hbar,rz,vsnap,d/931,197.3,0,.06,42,13,100,27,227.5, ! Original A.Bonasera prescription
c 1 30.,0.,200/
data m,hbar,rz,vsnap,d/931.,197.3,0.,.06,6.,6.,12.,12.,2,
1 30.,0.,200./
open(unit=7,file='CcWEres.txt',status='unknown')
open(unit=8,file='28Si28Si.txt',status='unknown')
open(unit=9,file='CcSstarWEres.txt',status='unknown')
write(8,*)'Ecm(MeV) sigtot(mb)'
c ecm=0.5 !modify
c do 202 j=1,12 !modify
c ecm=ecm+0.5 !modify
c print*,ecm !modify
dt=1.0 !time step calculation
xpi=3.1415
cobmax=1.
c print *,log(2.71)
c stop
c1 read *, n !see equivalence statement for parameter assignment
1 n=0
if(n.eq.0) go to 2
if(n.lt.0) stop
if(n.le.8) read *,d(n)
if(n.le.8) print*,n,d(n)
go to 1
2 r1=1.15*a1**.3333 !initialize the collision
r2=1.15*a2**.3333
rbar=r1*r2/(r1+r2)
n1=a1-z1
n2=a2-z2
ntot=n1+n2
ztot=z1+z2
xiso=(ntot-ztot)/(a1+a2)
sigma=0.9517*(1-xiso**2) !surface tension Angius2 PLB sigma=0.951699972
c sigma=1.0 !modify
c do 101 il=1,10
c print *,'a1,z1,a2,z2,ecm,b,sigma=',a1,z1,a2,z2,ecm,b,sigma
nstep=d(8)
rbass=1.0087*(r1+r2)-1.598*(1./r1+1./r2)
rpot=1.087*(r1+r2)
z2e2=z1*z2*1.44
m12=m*a1*a2/(a1+a2)
xh2=xpi*hbar**2/2./m12/ecm
rinner=r1+r2
c energy loop
deine=1.
ecm0=ecm-deine
do 312 ine=1,1
ecm0=ecm0+deine
ecm=ecm0
ecmini=ecm
c print*, 'NEW ENERGY=',ecmini
sigmatot=0.
sigmatota=0.
sigmatotn=0.
sigmatotp=0.
sigmatotg=0.
c do 313 lll=40,40,1 !original
do 313 lll=0,0 !modify
print*,'lll',lll
sss=r1+r2
ell=l*(l+1)*hbar**2/(2*2*a1*m*(r1**2+(r/2)**2)) !modify
if(r.gt.sss)then
ell=l*(l+1)*hbar**2/(2*m12*r**2 ) !modify
end if
if(ell.gt.ecm)then !modify,Centrifugal potential is less than total energy.
print*,'ell.gt.ecm'
goto 313
c
endif
xll=lll
ximg=1. !real time ximg=1.
img=0
action=0.
rz=0.
r=ra
ecml=ecm-xll*(xll+1)*hbar**2/2/m12/rinner**2 !shift the c.m. energy to take into account l
if(ecml.le.0.)go to 313
c if(ecm.le.z2e2/r)go to 313 !modify by sun
c li=b*sqrt(2*ecm*m12)/hbar
li=xll
b=xll/(sqrt(2*ecm*m12)/hbar)
epsilon=sqrt(1.+(2*ecm*b/z2e2)**2)
temp=(1.+2*ecm*b**2/z2e2/r)/epsilon
if(temp.lt.1.00001) go to 3
print*,'unphysical initial conditions'
go to 313
3 theta=acos(-temp)-acos(-1./epsilon)
x=r*cos(theta)
y=r*sin(theta)
vt=li*hbar/m12/r
vr=-sqrt(2.*(ecm-z2e2/r)/m12-vt**2)
vx=(vr*x-vt*y)/r
vy=(vr*y+vt*x)/r
td=8./3.*(r1+r2)/0.28/2.*2.0 ! time delay, see NPA85 Neck model, needs a correction between 1 and 2
do 4 i=1,td
4 vrt(i)=vr
t=0
stage=1
nt=0
k=1
ecou2=z2e2/ra
ebass2=0.
eneck=0.
xdiss=0.
etang=1./2.*m12*vt**2 !initial rotational energy
etang2=0.
cobass=1.
c5 read *,n
5 n=0
if(n.ne.0) go to 1 !carry the collision through several time steps
do 20 i=1,nstep !print after nsteps
t=t+dt
if(t.gt.30000.)go to 313
rold=sqrt(x**2+y**2)
xold=x
yold=y
x=x+vx*dt
y=y+vy*dt
r=sqrt(x**2+y**2)
deltar=r-rold
deltax=x-xold
deltay=y-yold
k=k+1
if(k.gt.td) k=1
th3=atan(vt*td*.67/r)
vr1=vrt(k)*cos(th3)+vt*sin(th3)
rx1=r*r2/(r1+r2)
if(rx1.lt.0.75*r2) rx1=0.75*r2
vr2=-sqrt((-0.28+vr1)**2+2*z2e2/a1/m*(1./(rx1+r1)-1./rx1))+.28
fc=z2e2/r**2
z=r-rpot
s1=r-r1-r2
go to (6,6,8,9,10), stage !stage of collision
6 if(z.lt.0) stage=2 ! 1: before touching
if(s1.lt.0) stage=3
s=amin1(r-rbass,20.)
c s=r-rbass
e1=exp(s/3.3)
e2=exp(s/.65)
fn=-r1*r2*(.00909*e1+.00938*e2)/(.03*e1+.0061*e2)**2/(r1+r2)
c cobass=0.9 !comment if you want resonances to modify the Bass potential
fn=1.08*fn*cobass !Bass force
fnbass=fn
ebass=-1.08*r1*r2/(.03*e1+.0061*e2)/(r1+r2) ! Bass potential
ebass2=-fn*deltar*ximg+ebass2
vr1=vr
dndt=7.42*r1*r2*(.0142*exp(-1.9*z)-.006*exp(-2.98*z))/rpot !sub barrier nucleon exchange
go to (15,7,15,1), stage
7 nt=nt+1
if(nt.gt.9) stage=3
go to 15
8 ry=amax1((r-r1-r2)/2.+.1,.1) ! 3: contact made here
stage=4
if(stage.eq.4)print*, 'nuclei approaching,e=',e ! 4: nuclei approaching
if(t.gt.3000.and.vr.le.0.)go to 30 !fusion
9 if(vr.gt.0)stage=5
if(stage.eq.5)print*,'nuclei rebounding'
c ss=(r1+r2+r+2*ry)*(r1+r2+2*ry-r)*(r+r1-r2)*(r+r2-r1)/4.
c rz=sqrt(ss)/r-ry ! Original Bertsch prescription.
if(vr.le.0.)then
rz=0. ! calculates the neck radius as section of spheres joined by a cylinder PLB84
do irz=1,1000000
rz=rz+0.0001
if(rz.gt.r1.or.rz.gt.r2)go to 121
rxx1=sqrt(r1**2-rz**2)
rxx2=sqrt(r2**2-rz**2)
hx1=r1-rxx1
hx2=r2-rxx2
rcyl=r-rxx1-rxx2
vcs1=hx1**2*(2.*r1+rxx1)/3.*xpi !volume spherical section nucleus 1
vcs2=hx2**2*(2.*r2+rxx2)/3.*xpi !volume spherical section nucleus 1
vcyl=rz**2*rcyl*xpi !volume cylinder
p=vcs1+vcs2-vcyl
if(p.lt.0.)go to 121
enddo
121 rzold=rz
c print *,'rNeck,r, approaching,vr=',rz,r,vr
endif
c ry=ry+.04/ry
go to 14
c10 rz=rz-.33*vr ! 5: nuclei rebounding. Original Bertsch prescription.
10 ssq=sqrt(4.*(r*rbar)**2-4.*r*(r1*r2*r-2.*(r1**3+r2**3)))-2.*rbar*r
rz=ssq/2./r !volume conservation PLB84
if(rz.gt.r1)rz=r1
if(rz.gt.r2)rz=r2 !neck radius smaller than initial nuclei radii
if(rz.gt.rzold)rz=rzold
c print *,'rNeck,r,rebounding,vr=',rz,r,vr
if(rz.lt.1.0) go to 23 !23:scission
c if(vr.gt.vsnap.or.vr1.gt.vsnap) go to 22
if(vr.lt.0.) go to 30 !30:fusion
if(s1.le.0.) go to 14
zx=.08*3.14*rz**2*s1/2
ra1=(-s1*zx/4.+z2*(r2+s1/2.))/(z2-zx)
ra2=(-s1*zx/4.+z1*(r1+s1/2.))/(z1-zx)
c fc=1.44*((z1-zx)*(z2-zx)/(r1+r2)**2+(z1-zx)*zx/(ra1+s1/2)**2
c 1 +(z2-zx)*zx/(ra2+s1/2)**2+zx**2/(s1**2/4.+rz**2/2.)) ! refined Coulomb force
fc=z2e2/r**2 !monopole approximation
14 fn=-6.28*sigma*rz !sigma
c if(rz.le.3.and.stage.eq.4)fn=fnbass
eneck=-fn*deltar*ximg+eneck
q1=1
vr1=vr2
dndt=0.009*rz**2*3.14 !window formula
15 continue
xpans=1.-exp(-vr1/0.28)
c if(ximg.eq.-1.)xpans=1.-cos(-vr1/0.28) !the exponent becomes as cosine in imagianry times
if(ximg.eq.-1.)xpans=1.-cos(-vr1/0.28) !modify by sun
fr=fn+fc-2.*m*0.28*xpans*dndt ! radial force, notice the exponential modification
xdiss=xdiss+2.*m*0.28*xpans*dndt*deltar
epot=epot-fn*deltar*ximg
ecou2=-fc*deltar*ximg+ecou2
fp=2*m*0.28*xpans*dndt
finertia=(r1**2/a1/(r1**2+(r*r1/(r1+r2))**2)+r2**2/a2/(r2**2+
1 (r*r2/(r1+r2))**2))*5./m/(r1+r2)
ft=-(vt-(li-1)*hbar*finertia)*m*dndt !tangential force
etang2=etang2-ft*deltay*ximg
vx=vx+(x*fr-y*ft)/m12/r*dt*ximg
vy=vy+(y*fr+x*ft)/m12/r*dt*ximg !Newton s equation
vr=(vx*x+vy*y)/r
prel=m12*vr/hbar
c if(vr.gt.0.and.img.eq.1)then
c ximg=1 ! return to real times
c rinner=r
c print*, 'real times t,2nd t.point=',t,r
c img=img+1
c endif
c if(vr.gt.0.and.img.eq.0)then
cc comment below if you do not want imaginary times
c ximg=-1 ! go to imaginary times
c print*, 'imaginary times t,1st t.point=',t,r
c img=img+1
c endif
if(vr.gt.0.and.img.eq.1)then !modify by sun
ximg=1 ! return to real times !modify by sun
rinner=r !modify by sun
print*, 'real times t,2nd t.point=',t,r
img=img+1 !modify by sun
endif
if(vr.gt.0.and.img.eq.0)then !modify by sun
ximg=-1 ! go to imaginary times !modify by sun
print*, 'imaginary times t,1st t.point=',t,r
img=img+1 !modify by sun
endif
c if(ximg.eq.-1.)then
c if(prel*deltar.ge.0.)action=action+prel*deltar !prel=m12*vr/hbar
c if(prel.eq.0.)ximg=1.
c endif
if(ximg.eq.-1.)then !modify by sun
if(prel*deltar.ge.0.)action=action+prel*deltar !modify by sun
if(prel.eq.0.)ximg=1. !modify by sun
endif
c if(r.lt.580.and.lll.eq.4)write(7,*)lll,r,prel
vt=(-vx*y+vy*x)/r
vrt(k)=vr
l=m12*vt*r/hbar
20 continue
c ecc=z2e2/r
vt=l*hbar/m12/r
c etang=1./2.*m12*vt**2
detang=etang2-etang
c impose energy conservation
105 e=ecou2+m12*(vx**2+vy**2)/2.+epot+xdiss-etang2
c if(e.gt.ecm+0.5)then
c vx=0.99*vx
c vy=0.99*vy
c print *,vx,vy
c go to 105
c endif
c if(e.lt.ecm-0.5)then
c vx=1.01*vx
c vy=1.01*vy
c print *,vx,vy
c go to 105
c endif
c
c print*,'t,r,rz,vr,vr1,l,e,ebass=',t,r,rz,vr,vr1,l,e,ebassq1
c print*,'enck,ebs2,eco2,edis=',eneck,ebass2,ecou2,xdiss
c print *,'dndt,stage,etang,etang2=',dndt,stage,etang,etang2
112 format(8f10.5)
go to 5
22 print*,'neck snaps'
c go to 24
c change l-value
go to 313
23 print*,'scission'
go to 313
24 eps=sqrt(1.+(1*hbar/z2e2)**2*2*e/m12)
thd=acos(-((l*hbar)**2/(m12*z2e2*r)+1)/eps)-acos(-1./eps)
thl=atan(y/x)
th2=atan(vy/vx)
c print*,x,y,thl,vx,vy,th2
theta=90.-(atan(-x/y)+thd)*57.3
c print*,theta,l,e
stop
c go to 1
30 print*, 'fusion'
c if(lll.eq.1.or.lll.eq.3.or.lll.eq.5.or.lll.eq.7)go to 313 !ang momentum and par cons in C+C fusion
c prob=1/(1.+exp(2.*action))
c prob=1. !modify by sun
c if(lll.eq.1.or.lll.eq.3.or.lll.eq.5.or.lll.eq.7)go to 313 !modify by sun
prob=1/(1.+exp(2.*action)) !modify by sun
c sigma=xpi*rinner**2*prob/100.
sigmatot=sigmatot+(2.*lll+1)*prob*xh2/100.
print*,'ecm,Action, prob,lll,sigma=',ecm,action,prob,lll,sigmatot
c call CCstatdecay(ecm,lll,proba,probn,probp)
c sigmatota=sigmatota+(2.*lll+1)*(prob*proba)*xh2/100
cc sigmatotp=sigmatotp+(2.*lll+1)*(prob*probp)*xh2/100
c sigmatotn=sigmatotn+(2.*lll+1)*(prob*probn)*xh2/100
c print *,'sigma(b) a,n,p=',sigmatota,sigmatotn,sigmatotp
313 continue
print*,'action',action
314 print *,'sigma(b),l_max=',sigmatot,lll
c print *,'sigma a,n,p=',sigmatota,sigmatotn,sigmatotp
c sigmareso=sigmatota+sigmatotn+sigmatotp
c print*,'sigma_resonance=',sigmareso
if(cobass.gt.cobmax)then
ecmax=ecmini
cobmax=cobass
endif
c write(8,*)ecmini,sigmatot,sigmatota,sigmatotn,sigmatotp
c cross section in mb the factor of 2 should be included for equal mass bosons
c write(8,*)ecmini,sigmatot*1000.*2 !modify by sun
c xlsstar=log(ecmini*sigmatot)+(87.12/sqrt(ecmini)+0.46*ecmini) !original
xlsstar=ecmini*sigmatot*exp(87.12/sqrt(ecmini)+0.46*ecmini) !modify by sun
print*,'ecmini,xlsstar',ecmini,xlsstar
write(9,*)ecmini,xlsstar
312 continue
c202 continue
print *,'comax,E=',cobmax,ecmax
stop
c go to 1
end
运行结果及报错内容
PS C:\Users\20499\Desktop\lunwen\tioashi> .\a.exe
lll 0
imaginary times t,1st t.point= 789.000000 25.9162922
real times t,2nd t.point= 2256.00000 6.26128387
nuclei approaching,e= 1.98791027
nuclei rebounding
scission
action 11.6816387
sigma(b),l_max= 0.00000000000000000000000000000000000 1
ecmini,xlsstar 2.00000000 0.00000000
comax,E= 1.00000000 0.00000000
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
PS C:\Users\20499\Desktop\lunwen\tioashi>
