njoy · whaeck · Oct 16, 2023 · Oct 5, 2023 · Oct 9, 2023 · Oct 10, 2023
diff --git a/ReleaseNotes.md b/ReleaseNotes.md
@@ -4,6 +4,7 @@ Given here are some release notes for NJOY2016. Each release is made through a f
 ## [NJOY2016.73](https://github.com/njoy/NJOY2016/pull/xxx)
 This update fixes the following issues:
   - Fix an issue in ACER for thermal scattering leading to energy values being out of order when plotting the coherent elastic scattering cross section (this issue only affects plots, the thermal scattering ACE files do not change).
+  - Increased allocation of an array in LEAPR to accommodate ENDF/B-VIII.1 thermal scattering evaluations and added a check to avoid an infinite loop when using a very fine beta grid. In addition, LEAPR will now warn the user about potential excessive calculation times and print out progression in the phonon expansion sum when the phonon expansion order is large.
 
 ## [NJOY2016.72](https://github.com/njoy/NJOY2016/pull/308)
 This update fixes the following issues:

diff --git a/src/leapr.f90 b/src/leapr.f90
@@ -47,6 +47,9 @@ module leapm
    real(kr),dimension(:),allocatable::dwpix,dwp1
    real(kr),dimension(:),allocatable::tempf,tempf1
 
+   ! min phonon expansion for time warning message
+   integer,parameter,public::maxnphon=250
+
 contains
 
    subroutine leapr
@@ -220,6 +223,7 @@ subroutine leapr
    integer::isym,mscr,maxb,isecs
    real(kr)::time
    character(4)::title(20)
+   character(60)::strng
   real(kr)::temp,emax
    character::text*80
    real(kr),dimension(:),allocatable::bragg
@@ -294,6 +298,12 @@ subroutine leapr
    read(nsysi,*) (alpha(i),i=1,nalpha)
    read(nsysi,*) (beta(i),i=1,nbeta)
 
+   ! warn for excessive computation time
+   if (nphon.gt.maxnphon) then
+      write(strng,'('' phonon expansion order is larger than '',i3)') maxnphon
+      call mess('leapr',strng,'calculation time may be excessive')
+   endif
+
    !--open the output unit
    call openz(nout,1)
 
@@ -373,11 +383,11 @@ subroutine leapr
          if (nd.gt.0) call discre(itemp)
 
          !--check for special hydrogen and deuterium options
-          if (ncold.gt.0) call coldh(itemp,temp)
+         if (ncold.gt.0) call coldh(itemp,temp)
 
          !--check for skold option for correlations
-          if ((nsk.eq.2) .and. (ncold.eq.0))&
-             call skold(itemp,temp,ssm,nalpha,nbeta,ntempr)
+         if ((nsk.eq.2) .and. (ncold.eq.0))&
+            call skold(itemp,temp,ssm,nalpha,nbeta,ntempr)
 
       !--continue temperature loop
       enddo
@@ -445,14 +455,15 @@ subroutine contin(temp,itemp,np,maxn)
    !--------------------------------------------------------------------
    use physics ! provides pi
    use mainio  ! provides nsyso
+   use util    ! provides timer
    ! externals
    real(kr)::temp
    integer::itemp,np,maxn
    ! internals
    integer::i,j,k,n,npn,npl,iprt,jprt
-   integer,dimension(1000)::maxt
+   integer,dimension(10000)::maxt
    character(3)::tag
-   real(kr)::al,be,bel,ex,exx,st,add,sc,alp,alw,ssct,ckk
+   real(kr)::al,be,bel,ex,exx,st,add,sc,alp,alw,ssct,ckk,time
    real(kr)::ff0,ff1,ff2,ff1l,ff2l,sum0,sum1
    real(kr),dimension(:),allocatable::p,tlast,tnow,xa
    real(kr),parameter::therm=0.0253e0_kr
@@ -498,8 +509,14 @@ subroutine contin(temp,itemp,np,maxn)
    do j=1,nbeta
       maxt(j)=nalpha+1
    enddo
-   if (iprint.eq.2)&
-     write(nsyso,'(/'' normalization check for phonon expansion'')')
+   if (iprint.eq.2) then
+      write(nsyso,'(/'' normalization check for phonon expansion'')')
+   endif
+   if (maxn.gt.maxnphon) then
+      call timer(time)
+      write(nsyse,'(/'' performing phonon expansion sum'',&
+            &37x,f8.1,''s'')'),time
+   endif
    do n=2,maxn
       npn=np+npl-1
       call convol(p,tlast,tnow,np,npl,npn,deltab,ckk)
@@ -525,7 +542,18 @@ subroutine contin(temp,itemp,np,maxn)
          tlast(i)=tnow(i)
       enddo
       npl=npn
+      if (mod(n,maxnphon).eq.0) then
+         call timer(time)
+         write(nsyse,'(2x,i5,'' of '',i5,&
+               &'' loops done for phonon expansion sum'',17x,f8.1,''s'')')&
+               &n,maxn,time
+      endif
    enddo
+   if (maxn.gt.maxnphon) then
+      call timer(time)
+      write(nsyse,'(/'' done with phonon expansion sum'',&
+            &38x,f8.1,''s'')'),time
+   endif
 
    !--print out start of sct range for each beta
    if (iprint.ne.0) then
@@ -865,8 +893,8 @@ subroutine trans(itemp)
       if (ded.lt.delta) delta=ded
       nu=1
       if (iprt.eq.1.and.iprint.eq.2) write(nsyso,&
-        '(/'' delta d='',f12.6,5x,''delta b='',f12.6,&
-        &10x,''delta='',f12.6)') ded,deb,delta
+        '(/'' delta d='',e18.5,5x,''delta b='',e18.5,&
+        &10x,''delta='',e18.5)') ded,deb,delta
 
       !--make table of s-diffusion or s-free on this interval
       call stable(ap,sd,nsd,al,delta,iprt,nu,ndmax)
@@ -1206,6 +1234,10 @@ subroutine sbfill(sb,nbt,delta,be,s,betan,nbeta,nbe,ndmax)
       else
          sb(i)=0
       endif
+      ! if delta is to small for the current value of beta, increase it
+      do while (bet.eq.(bet+delta))
+         delta=delta*10
+      end do
       bet=bet+delta
    enddo
    return

diff --git a/tests/80/CMakeLists.txt b/tests/80/CMakeLists.txt
@@ -0,0 +1,9 @@
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/input"
+               "${CMAKE_CURRENT_BINARY_DIR}/input" COPYONLY )
+
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/referenceTape24"
+               "${CMAKE_CURRENT_BINARY_DIR}/referenceTape24" COPYONLY )
+
+add_test( NAME "Test80"
+          COMMAND ${Python3_EXECUTABLE} "../execute.py"
+          WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}" )