Merge pull request #321 from njoy/develop

NJOY2016.73

ReleaseNotes.md

@@ -1,6 +1,14 @@
 # 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/321)
+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.
+  - Added logic to MODER to read background R-matrix element information from LRF=7 resonance parameter data.
+  - Updated RECONR to use background R-matrix element information from LRF=7 and added test 81 using ENDF/B-VIII.1 Sr88.
+  - 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
@@ -413,7 +423,13 @@ subroutine leapr
    isym=0
    if (ncold.ne.0) isym=1
    if (isabt.eq.1) isym=isym+2
-   mscr=4000
+
+   ! Based on endout, to write the actual TSL data, the max number of entries
+   ! needed in scr is either 8+2*nalpha, or 8+2*nedge. However, we have no way
+   ! of knowing how many comment lines were added to the leaper input. The
+   ! previous hard coded limit of 4000 is also used as a possible max as this
+   ! has apparently been sufficient to hold all comments in the past.
+   mscr = max(8 + 2*nalpha, 8 + 2*nedge, 4000)
    allocate(scr(mscr))
    call endout(ntempr,bragg,nedge,maxb,scr,mscr,isym,ilog)
 
@@ -445,14 +461,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,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
@@ -468,6 +485,7 @@ subroutine contin(temp,itemp,np,maxn)
    allocate(tlast(nphon*np1))
    allocate(tnow(nphon*np1))
    allocate(xa(nalpha))
+   allocate(maxt(nbeta))
 
    !--calculate various parameters for this temperature
    call start(itemp,p,np,deltab,tev)
@@ -498,8 +516,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 +549,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 +900,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 +1241,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/moder.f90

@@ -693,18 +693,25 @@ subroutine file2a(nin,nout,nscr,a)
                      call moreio(nin,nout,nscr,a,nb,nw)
                   enddo
                   if (kbk.gt.0) then
-                     call listio(nin,nout,nscr,a,nb,nw)
-                     lbk=n1h
-                     if (lbk.eq.1) then
-                        call tab1io(nin,nout,nscr,a,nb,nw)
-                        do while (nb.ne.0)
-                           call moreio(nin,nout,nscr,a,nb,nw)
-                        enddo
-                        call tab1io(nin,nout,nscr,a,nb,nw)
-                        do while (nb.ne.0)
-                           call moreio(nin,nout,nscr,a,nb,nw)
-                        enddo
-                     endif
+                     do l=1,kbk
+                        call contio(nin,nout,nscr,a,nb,nw)
+                        lbk=l2h
+                        if (lbk.eq.1) then
+                           call tab1io(nin,nout,nscr,a,nb,nw)
+                           do while (nb.ne.0)
+                              call moreio(nin,nout,nscr,a,nb,nw)
+                           enddo
+                           call tab1io(nin,nout,nscr,a,nb,nw)
+                           do while (nb.ne.0)
+                              call moreio(nin,nout,nscr,a,nb,nw)
+                           enddo
+                        else if (lbk.eq.2.or.lbk.eq.3) then
+                           call listio(nin,nout,nscr,a,nb,nw)
+                           do while (nb.ne.0)
+                              call moreio(nin,nout,nscr,a,nb,nw)
+                           enddo
+                        endif
+                     enddo
                   endif
                   if (kps.eq.1)then
                      call listio(nin,nout,nscr,a,nb,nw)