m0_62881912 2022-09-25 18:51 采纳率: 33.3%
浏览 58
已结题

Fortran77运行结果出现浮点数下沉

问题遇到的现象和发生背景

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>
  • 写回答

3条回答 默认 最新

  • _GX_ 2022-09-26 06:57
    关注

    浮点运算出现underflow,一般是除法运算时被除数很小且除数很大。

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论
查看更多回答(2条)

报告相同问题?

问题事件

  • 系统已结题 10月5日
  • 已采纳回答 9月27日
  • 创建了问题 9月25日

悬赏问题

  • ¥15 基于RTKLIB框架写的精密单点定位-AR
  • ¥15 PFENet的预训练权重
  • ¥15 程序哪有错误怎么改?
  • ¥15 交换机和交换机之间的链路带宽以及主机带宽的理解
  • ¥15 ai创想家对战模式代码
  • ¥15 集合A由3个2行4列二维数组构成,从集合A中任意取一个二维数组元素、如果该二维数组元素的对应列位置的上、下两数都是奇数,而且仅有2个列是奇数/奇数,则该数组有意义,并放入集合B中打印输出。
  • ¥15 电信IPV6 无法外网访问吗
  • ¥15 有偿求效果比较好的遥感影像匹配的c++代码
  • ¥15 博主,你好,我下载了你的智能网联汽车辅助驾驶安全信息检测系统,现在不会运行,可以教我吗,
  • ¥15 怎么在excle输入下列公式