Merge remote-tracking branch 'origin/develop' into fix/leapr_array

ReleaseNotes.md

@@ -1,6 +1,12 @@
 # Release Notes—NJOY2016
 Given here are some release notes for NJOY2016. Each release is made through a formal [Pull Request](https://github.com/njoy/NJOY2016/pulls) made on GitHub. There are links in this document that point to each of those Pull Requests, where you can see in great details the changes that were made. Often the Pull Requests are made in response to an [issue](https://github.com/njoy/NJOY2016/issues). In such cases, links to those issues are also given.
 
+## [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.
+  - Fixing a few thing related to intel compiler warnings and errors.
+
 ## [NJOY2016.72](https://github.com/njoy/NJOY2016/pull/308)
 This update fixes the following issues:
   - Fixed an issue in GROUPR related to an error coming up in production matrix calculations. Depending on when a user asks for a production matrix associated to a reaction, it is possible that the reference frame of the previous reaction is used instead (caused by erronously defining an already declared global variable as local with a "save" attribute). In some circumstances, this causes NJOY2016 to error out (with a message related to unsupported reference frames). No test results had to be updated due to this change.

src/aceth.f90

@@ -1209,7 +1209,7 @@ subroutine tplots(nout,hk)
       endif
       write(nout,'(''0/'')')
       if (idpnc.eq.4.or.idpnc.eq.5) then
-         e=xss(itce+1)-xss(itce+1)/1000
+         e=xss(itce+1)
          xs=xss(itce+nee+1)/e
          xs=xs/100
          write(nout,'(1p,2e14.6,''/'')') e,xs
@@ -1220,7 +1220,7 @@ subroutine tplots(nout,hk)
          if (idpnc.eq.4.or.idpnc.eq.5) xs=xs/e
          write(nout,'(1p,2e14.6,''/'')') e,xs
          if ((idpnc.eq.4.or.idpnc.eq.5).and.i.lt.nee) then
-            e=xss(itce+i+1)-xss(itce+i+1)/1000
+            e=xss(itce+i+1)
             xs=xss(itce+nee+i)/e
             write(nout,'(1p,2e14.6,''/'')') e,xs
          endif
@@ -1280,7 +1280,7 @@ subroutine tplots(nout,hk)
          endif
       enddo
       if (idpnc.eq.4.or.idpnc.eq.5) then
-         e=xss(itce+1)-xss(itce+1)/1000
+         e=xss(itce+1)
          j=0
          do while (j.lt.nie)
             j=j+1
@@ -1333,7 +1333,7 @@ subroutine tplots(nout,hk)
          tot=xs+xn+xielas
          write(nout,'(1p,2e14.6,''/'')') e,tot
          if (i.lt.nee.and.(idpnc.eq.4.or.idpnc.eq.5)) then
-            e=xss(itce+i+1)-xss(itce+i+1)/1000
+            e=xss(itce+i+1)
             xs=xss(itce+nee+i)/e
             j=0
             do while (j.lt.nie)
@@ -2039,7 +2039,7 @@ subroutine tplots(nout,hk)
                if (k.eq.nang.and.un-u.gt.5*(u-ul)) un=u+3*(u-ul)
                p=1
                p=p/nang
-               p=p/(un-ul)
+               if (un.ne.ul) p=p/(un-ul)
                ul=un
             enddo
          endif
@@ -2083,7 +2083,7 @@ subroutine tplots(nout,hk)
                   if (k.eq.nang.and.un-u.gt.5*(u-ul)) un=u+3*(u-ul)
                   p=1
                   p=p/nang
-                  p=p/(un-ul)
+                  if (un.ne.ul) p=p/(un-ul)
                   if (k.eq.1) write(nout,'(1p,2e14.6,''/'')') ul,zmin
                   write(nout,'(1p,2e14.6,''/'')') u,p
                   if (k.eq.nang) write(nout,'(1p,2e14.6,''/'')') un,zmin

src/groupr.f90

@@ -7444,7 +7444,7 @@ subroutine getmf6(ans,ed,enext,idisc,yld,eg,ng,nl,iglo,ng2,nq,&
    real(kr)::zad,elo,ehi,apsx,enow,eihi,ep,epnext,en
    real(kr)::pspmax,yldd,el,eh,e0,g0,e1,e2,test,pe,disc102
    real(kr)::val,fx,ex,cx,cxx,rn,dx
-   integer(kr)::nx,ncyc,n,ix
+   integer::nx,ncyc,n,ix
    integer,parameter::mxlg=65
    real(kr)::term(mxlg),terml(mxlg)
    integer,parameter::maxss=500

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,allocatable,dimension(:)::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
@@ -499,8 +510,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)
@@ -526,7 +543,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
@@ -866,8 +894,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)
@@ -1207,6 +1235,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

src/thermr.f90

@@ -2494,7 +2494,7 @@ real(kr) function sig(e,ep,u,tev,nalpha,alpha,nbeta,beta,sab)
    ! internals
    integer::nb1,na1,i,ib,ia
    real(kr)::rtev,bb,a,sigc,b,c,bbb,s,s1,s2,s3,arg
-   real(kr)::tfff,tfff2,rat
+   real(kr)::tfff,tfff2
    real(kr),parameter::sigmin=1.e-10_kr
    real(kr),parameter::sabflg=-225.e0_kr
    real(kr),parameter::amin=1.e-6_kr

src/vers.f90

@@ -3,6 +3,6 @@ module version
    ! These values are updated during the NJOY revision-control process.
    implicit none
    private
-   character(8),public::vers='2016.72'
-   character(8),public::vday='29Sep23'
+   character(8),public::vers='2016.73'
+   character(8),public::vday='xxDec23'
 end module version