Merge pull request #264 from njoy/develop

NJOY2016.68

.github/workflows/ContinuousIntegration.yml

@@ -18,7 +18,7 @@ jobs:
     runs-on: ${{matrix.os}}
     strategy:
       matrix:
-        os: [ubuntu-18.04, macos-10.15]
+        os: [ ubuntu-20.04, macos-11 ]
         build_type: [ Debug, Release ]
         vFortran: [ gfortran-9 ]
 

ReleaseNotes.md

@@ -1,6 +1,15 @@
 # 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.68](https://github.com/njoy/NJOY2016/pull/264)
+This update fixes a number of minor issues:
+  - there was an indexing error in the calculation of inelastic thermal scattering mubar in ACER for IFENG=2 ACE files. Test 74 was added to track this issue.
+  - there was a division by zero that caused a segmentation fault in Debug mode, also for IFENG=2 ACE files.
+  - removed an erroneously added FEND record in the PENDF file when no MF12 data gets added to the PENDF file (related to issue #250).
+  - correction in ACER to fix issue #250. Previously it could happen that ACER picked the wrong temperature from a multi-temperature PENDF file when the materials in the PENDF only contained linearised MF3 data (no MF10, MF12 or MF23 data was added by RECONR).
+  - fixed an indexing problem in CCCCR for delayed neutron data and added test 76.
+  - there was a use of wrong memory addresses in CCCR not accounting for GROUPR's energy structure offset.
+
 ## [NJOY2016.67](https://github.com/njoy/NJOY2016/pull/243)
 This update resolves a number of issues encountered when processing the newly released JENDL5 nuclear data library and TENDL-2021 library for evaluations with sub-actinide fission, and provides fixes for a few minor issues.
 

docs/testDescription.md

@@ -363,3 +363,21 @@ Tests 72 was added as a consequence of issue [\#188](https://github.com/njoy/NJO
 [[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/73/input)]
 
 Tests 73 was added as a consequence of issue [\#241](https://github.com/njoy/NJOY2016/issues/241). It verifies that GROUPR automated processing of MF10 for an evaluation with subactinide fission functions properly.
+
+## Test Problem 74
+
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/74/input)]
+
+Tests 74 was added to track an issue found in ACER while calculating inelastic thermal scattering mubar for IFENG=2 ACE files. This is a test for H in ZrH and is similar to test 69 (test 69 is for Zr in ZrH and is an IFENG=1 file).
+
+## Test Problem 75
+
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/75/input)]
+
+Tests 75 was added to track an issue found in ACER when multiple temperatures were present in the PENDF file. See issue [\#250](https://github.com/njoy/NJOY2016/issues/241) for more information.
+
+## Test Problem 76
+
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/76/input)]
+
+Tests 76 was added to track a few CCCCR issues. Unfortunately, these issues were related to a binary file so we cannot verify the content of the file just yet. At the very least, we now have a test to verify if CCCCR actually runs (there was previously a crash on this input).

src/acefc.f90

@@ -594,6 +594,8 @@ subroutine first(nendf,npend,nout,nethr,itab,matd,tempd,newfor,iopp)
    enddo
    if (jethr.gt.1) call aordr(jethr,nethr,ethr)
    nethr=jethr
+   ! we're on a fend record, go back one line (this avoids the next findf
+   ! call to see a mend record if the pendf tape only contains mf3)
    call skiprz(npend,-2)
 
    !--for new format...
@@ -3943,7 +3945,12 @@ subroutine convr(nin,npend,nout,nscr,nedis,nethr,matd)
    call findf(matd,3,0,npend)
    nnth=0
   101 call contio(npend,0,0,scr,nb,nw)
-   if (mfh.eq.0) go to 102
+   if (mfh.eq.0) then
+      ! we're on a fend record, go back one line (this avoids the next findf
+      ! call to see a mend record if the pendf tape only contains mf3)
+      call skiprz(npend,-2)
+      go to 102
+   endif
    e=0
    call gety1(e,enxt,jdis,x,npend,scr)
    nnth=nnth+1

src/aceth.f90

@@ -1454,7 +1454,7 @@ subroutine tplots(nout,hk)
 !                       xss(loc+3+k)
 !              endif
               ubar=ubar+xss(loc+3+k)*p/2
-              ubar=ubar+xss(loc+3+k-nang-2)*p/2
+              ubar=ubar+xss(loc+3+k-nang-3)*p/2
               sum=sum+p
             enddo
             cdl=xss(loc+3)
@@ -1526,7 +1526,7 @@ subroutine tplots(nout,hk)
             p=xss(loc+3)-cdl
             do k=1,nang
               ubar=ubar+xss(loc+3+k)*p/2
-              ubar=ubar+xss(loc+3+k-nang-2)*p/2
+              ubar=ubar+xss(loc+3+k-nang-3)*p/2
               sum=sum+p
             enddo
             cdl=xss(loc+3)
@@ -2068,7 +2068,7 @@ subroutine tplots(nout,hk)
       do j=2,nbini-1
          ep=(xss(loc+1)+epl)/2
          if (xss(loc+1).gt.epl) then
-            if (j.eq.2.or.ep/epl.gt.skip) then
+            if (j.eq.2.or.ep.gt.skip*epl) then
                write(nout,'(1p,e14.6,''/'')') ep
                ul=-1
                write(nout,'(1p,2e14.6,''/'')') -1.,zmin

src/ccccr.f90

@@ -6,11 +6,11 @@ module ccccm
    public ccccr
 
    ! equivalenced arrays for CCCC input and output
-   integer::isiza=50000    ! container array size
-   real(k4)::a(50000)      ! reals are 4-byte
-   integer(k4)::ia(50000)  ! integers are 4-byte
-   real(k8)::ha(25000)     ! Hollerith data are 8-byte
-   character(8)::ta(25000) ! text equivalent to Hollerith
+   integer,parameter:: isiza4=100000, isiza8=isiza4/2
+   real(k4)::a(isiza4)      ! reals are 4-byte
+   integer(k4)::ia(isiza4)  ! integers are 4-byte
+   real(k8)::ha(isiza8)     ! Hollerith data are 8-byte
+   character(8)::ta(isiza8) ! text equivalent to Hollerith
    integer::mult=2         ! used for counting 8-byte entries
    equivalence(a(1),ia(1),ha(1),ta(1))
 
@@ -1061,7 +1061,7 @@ subroutine isomtx(ia,b,l21,l22,loci,nsz,irec2)
    ! for ifopt=2 and nsblk=ngps, nrec and npass are not used
 
    !--set up sizes
-   nsiza=isiza-l22+1
+   nsiza=isiza4-l22+1
    irsize=6+(ngps+1)*(maxord+1)*nsz
    if (irsize.lt.6*ngps) irsize=6*ngps
    if (irsize.gt.nsiza)&
@@ -1543,7 +1543,7 @@ subroutine pisotx(n4)
    nwi=niso
    nwds=mult*nwh+nwr+nwi
    irec1=irec1+1
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pisotx','input record too large.',' ')
    read(n4)(ha(i),i=1,nwh),(a(llr+i),i=1,nwr),(ia(lli+i),i=1,nwi)
    if (ird(2).eq.1) then
@@ -1582,7 +1582,7 @@ subroutine pisotx(n4)
       llr=npt-1
       lli=llr+nwr
       irec1=irec1+1
-      if (nwds.ge.isiza)&
+      if (nwds.ge.isiza4)&
         call error('pisotx','input record too large.',' ')
       read(n4)(a(llr+i),i=1,nwr),(ia(lli+i),i=1,nwi)
       if (ird(3).eq.1) then
@@ -1600,7 +1600,7 @@ subroutine pisotx(n4)
    nwi=11+2*nscmax+2*ngroup*nscmax
    nwds=mult*nwh+nwr+nwi
    irec1=irec1+1
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pisotx','input record too large.',' ')
    read(n4)(ha(j),j=1,nwh),(a(llr+j),j=1,nwr),(ia(lli+j),j=1,nwi)
    jread=1+nwds
@@ -1701,7 +1701,7 @@ subroutine pisotx(n4)
    if (ind.eq.1) read(hsnd,'(a6)') htype(j)
    if (ind.eq.1) j=j+1
    if (int.eq.1) read(hsnt,'(a6)') htype(j)
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pisotx','input record too large.',' ')
    read(n4)(a(jread+k-1),k=1,nwds)
    if (ird(5).eq.1) then
@@ -1732,7 +1732,7 @@ subroutine pisotx(n4)
       llr=npt-1
       lli=llr+nwr
       irec1=irec1+1
-      if (nwds.ge.isiza)&
+      if (nwds.ge.isiza4)&
          call error('pisotx','input record too large.',' ')
       read(n4)(a(llr+i),nn=1,nwr),(ia(lli+nn),nn=1,nwi)
       if (ird(6).eq.1) then
@@ -1798,7 +1798,7 @@ subroutine pisotx(n4)
    nwds=nwds*lordn
    if (nwds.eq.0) go to 180
    irec1=irec1+1
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pisotx','input record too large.',' ')
    read(n4)(a(npt+ii-1),ii=1,nwds)
    if (ird(7).ne.1) go to 180
@@ -2104,7 +2104,7 @@ subroutine brkdat
    !--iffile--length of entire self shielding file
    !--nsblk--number of s-s records to be written on scratch tape
    !--lrsize--length of one s-s record
-   nosize=isiza-l20+1
+   nosize=isiza4-l20+1
    iffile=nreact*ngpn*ntl*nzl
    nsblk=1
    if (iffile.gt.nosize) nsblk=nreact
@@ -2128,7 +2128,7 @@ subroutine brkdat
    call repoz(-nscrt2)
    jbl=ngpn
    jbh=0
-   do j=l18,isiza
+   do j=l18,isiza4
       a(j)=0
    enddo
    do j=1,10
@@ -2838,7 +2838,7 @@ subroutine pbrkxs(n4)
    nwi=4*nisosh
    nwds=mult*nwh+nwr+nwi
    irec1=irec1+1
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pbrkxs','input record too large.',' ')
    read(n4)(ha(i),i=1,nwh),(a(llr+i),i=1,nwr),(ia(lli+i),i=1,nwi)
    if (ird(2).eq.1) then
@@ -2904,7 +2904,7 @@ subroutine pbrkxs(n4)
    nwds=nreact*npro
    if (nwds.le.0) go to 100
    irec1=irec1+1
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pbrkxs','input record too large.',' ')
    read(n4)(a(jssf+j-1),j=1,nwds)
    if (ird(3).ne.1) go to 100
@@ -2979,7 +2979,7 @@ subroutine pbrkxs(n4)
    !--cross sections
    nwds=6*ngroup
    irec1=irec1+1
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pbrkxs','input record too large.',' ')
    read(n4)(a(jssf+n-1),n=1,nwds)
    if (ird(4).eq.1) then
@@ -3049,7 +3049,7 @@ subroutine cdlyxs(idlay,idlayt)
       nwi=2*nisod
       nwds=mult*nwh+nwr+nwi
       write(idlay)(ha(l1h+i),i=1,nwh),(a(l1r+i),i=1,nwr),&
-        (ia(lli+i),i=1,nwi)
+        (ia(l1i+i),i=1,nwi)
 
       !--delayed neutron precursor yield data by isotope
       l3=l1+nwds
@@ -3103,7 +3103,7 @@ subroutine dldata(nin)
    integer::nin
    ! internals
    integer::nisoo,nwds,lzero,l3,l4,l5,l6,l7,l8,l9,la,last
-   integer::i,ipr,nb,nw,ig,nl,nz,ng,ng2,loc,ngn,ngp1,loca
+   integer::i,ipr,nb,nw,ig,ig1,nl,nz,ng,ng2,loc,ngn,ngp1,loca
    integer::ifam,nd,jj,locb,j,ndg
    integer::imusd(200)
    integer,parameter::ndmax=8
@@ -3125,6 +3125,7 @@ subroutine dldata(nin)
          call contio(nin,0,0,e(1),nb,nw)
       enddo
       if (mfh.eq.5) then
+         ! number of delayed neutron groups
          ndg=nint(e(3))
       else
          nisod=0
@@ -3172,6 +3173,9 @@ subroutine dldata(nin)
    la=l9+5
    last=la+10*ngps
 
+   if (last.gt.isiza4) call error('dldata',&
+        'max4a for container a exceeded',' ')
+
    !--read through input tape for desired data
    do i=1,200
       imusd(i)=0
@@ -3244,6 +3248,10 @@ subroutine dldata(nin)
 
    !--process delayed neutron spectra record
   410 nd=nl
+   ! lowest energy group that holds delayed data
+   ! subtraction of -1 because first position in gendf mf5mt455
+   ! is decay constant of that precursor family
+   ig1=l1h-1
    ! hisnm--isotope name
    jj=(l1-1)/mult+nisod
    ta(jj)=hisnm
@@ -3263,7 +3271,7 @@ subroutine dldata(nin)
       loca=l3+ngn+ngn*(ifam-1)
       locb=6+i
       fract(i)=0
-      do j=1,ngn
+      do j=1,ig1
          a(loca-j)=e(nl*j+locb)
          fract(i)=e(nl*j+locb)+fract(i)
       enddo
@@ -3339,7 +3347,7 @@ subroutine pdlyxs(idlay)
    lli=llr+nsp
    nwds=nwh*mult+nsp+nint
    irec1=irec1+1
-   if (nwds.ge.isiza)&
+   if (nwds.ge.isiza4)&
      call error('pdlyxs','input record too large.',' ')
    read(n4)(ha(i),i=1,nwh),(a(llr+i),i=1,nsp),(ia(lli+i),i=1,nint)
    write(nsyso,'(/3x,''isotope'',5x,''name'')')
@@ -3393,7 +3401,7 @@ subroutine pdlyxs(idlay)
       llr=jwd1-1
       lli=llr+nsp
       irec1=irec1+1
-      if (nwds.ge.isiza)&
+      if (nwds.ge.isiza4)&
         call error('pdlyxs','input record too large.',' ')
       read(n4)(a(llr+j),j=1,nsp),(ia(lli+j),j=1,nint)
       jwd2=jwd1+ngroup*kk
@@ -3516,4 +3524,3 @@ subroutine stow(a,b,nwds)
    end subroutine stow
 
 end module ccccm
-

src/groupr.f90

@@ -662,7 +662,7 @@ subroutine groupr
          izam=-1
          izar=-1
       else
-         izar=(mfd-((mfd/1000000)*1000000))/10
+         izar=(mfd-((mfd/10000000)*10000000))/10
          if (lfs.lt.10) then
             izam=mod(mfd,10000000)+lfs
          else
@@ -1444,14 +1444,14 @@ subroutine namer(izad,izam,mfd,mtd,mtname)
       izaa=izam/10
       imm=isom
       if (imm.eq.0) then
-         write(azam,'(i5)') izaa
-         dummy='('//proj//','//reac//')-'//azam(1:5)//'-production.'
+         write(azam,'(i6)') izaa
+         dummy='('//proj//','//reac//')-'//azam(1:6)//'-production.'
       else if (imm.lt.10) then
-         write(azam,'(i5,''m'',i1)') izaa,imm
+         write(azam,'(i6,''m'',i1)') izaa,imm
          dummy='('//proj//','//reac//')-'//azam//'-production.'
       else
          izaa=izaa/10
-         write(azam,'(i5,''m'',i2)') izaa,imm
+         write(azam,'(i6,''m'',i2)') izaa,imm
          dummy='('//proj//','//reac//')-'//azam//'-production.'
       endif
    else

src/reconr.f90

@@ -1773,7 +1773,7 @@ subroutine lunion(nin,nout,ngrid,ngr)
    ! internals
    integer::ndim,iold,inew,i,j,mfl,nb,nw,nss,iss
    integer::l,lsave,in,ig,ir,jr,lr,ibase,nbta,inta,intn
-   integer::ipwr,inx,ngn,isave,ngneg,jg,igt,ngpos,npr
+   integer::ipwr,inx,ngn,isave,ngneg,jg,igt,ngpos,npr,nmf12
    character(40)::text
    real(kr)::stpmax,awrx,eg,qx,thrx,thrxx,et
    real(kr)::er,sr,ernext,srnext,enl,test,snl,erl,srl
@@ -1799,6 +1799,7 @@ subroutine lunion(nin,nout,ngrid,ngr)
    !--initialize.
    if (eresr.lt.elim) elim=eresr
    ndim=50
+   nmf12=0
    allocate(bufo(nbuf))
    allocate(bufn(nbuf))
    allocate(x(ndim))
@@ -1841,14 +1842,19 @@ subroutine lunion(nin,nout,ngrid,ngr)
    go to 410
   120 continue
    if (mfh.gt.0) go to 130
-   if (mfl.eq.3.or.mfl.eq.10.or.mfl.eq.12.or.mfl.eq.13.or.mfl.eq.23)&
+   ! always add a fend record for MF3, MF10, MF13 and MF23
+   if (mfl.eq.3.or.mfl.eq.10.or.mfl.eq.13.or.mfl.eq.23) then
      call afend(nout,0)
+   endif
+   ! only add a fend record for MF12 if we have written sections to nout
+   if (mfl.eq.12.and.nmf12.ne.0) then
+     call afend(nout,0)
+   endif
    mfl=0
    go to 110
-   ! process mf3, mf10, mf13, and mf23 only.
+   ! process mf3, mf10, mf12, mf13 and mf23 only.
   130 continue
-   if (mfh.eq.3.or.mfh.eq.10.or.mfh.eq.13.or.mfh.eq.23) go to 140
-   if (mfh.eq.12) go to 140
+   if (mfh.eq.3.or.mfh.eq.10.or.mfh.eq.12.or.mfh.eq.13.or.mfh.eq.23) go to 140
    call tofend(nin,0,0,scr)
    go to 110
    ! skip redundant reactions and non-cross-sections in n files
@@ -1889,6 +1895,7 @@ subroutine lunion(nin,nout,ngrid,ngr)
 
    !--process this reaction
   180 continue
+   if (mfh.eq.12) nmf12=nmf12+1
    call contio(0,nout,0,scr,nb,nw)
    if (mfh.ne.23.and.awin.ne.zero) awrx=c2h/awin
   181 continue

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.67'
-   character(8),public::vday='26May22'
+   character(8),public::vers='2016.68'
+   character(8),public::vday='30Sep22'
 end module version

tests/21/referenceTape24

@@ -66863,6 +66863,5 @@ Processed by NJOY                                                 2637 1451    6
  1.300000+8 0.000000+0 1.400000+8 0.000000+0 1.500000+8 0.000000+026371010366478
                                                                   263710  066479
                                                                   2637 0  066480
-                                                                  2637 0  066481
-                                                                     0 0  066482
+                                                                     0 0  066481
                                                                     -1 0  0    0