Merge pull request #295 from njoy/develop

NJOY2016.70

.github/workflows/ContinuousIntegration.yml

@@ -27,7 +27,7 @@ jobs:
       run: |
         which ${{matrix.vFortran}}
         ${{matrix.vFortran}} --version
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - name: mkdir bin
       run: mkdir bin
     - name: cmake

ReleaseNotes.md

@@ -1,6 +1,18 @@
 # 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.70](https://github.com/njoy/NJOY2016/pull/295)
+This update fixes a number of minor issues:
+  - Fixed an issue in HEATR when reading evaluations with large multiplicity tables in MF6.
+  - Fixed an issue in HEATR when calculating the average outgoing energy from a distribution that uses multiple interpolation ranges in TAB1 records (test 79 was added to detect this issue in the future). Mainly nuclides using MF5 instead of MF6 are impacted by this change (e.g. Sn119 and Sn122 from ENDF/B-VIII.0).
+  - Fixed an issue in HEATR where the photon recoil needed to be multiplied by the photon multiplicity to obtain the photon recoil per interaction.
+  - Fixed a crash in THERMR when asking for S(a,b) processing (iinc=2) while no ENDF tape is given (nendf=0).
+  - Multiple ERRORR calls can now be made in the same input file without crashing. This is of interest to users that wish to process MF34 and MF35 (where ERRORR needs to be called for each sub-subsection and incident energy group). The issue was related to arrays being allocated but not unallocated in the previous ERRORR run in NJOY's Sammy routines (evaluations using MF2 LRF=7 had this issue).
+  - Fixed an issue in ACER where the number of photons given in the ENDF file was larger than the hardcoded limit. The new limit is now adaptive.
+  - Fixed an issue in ACER where NaN values were produced in the postscript file for the checking plots.
+
+A few compiler warnings have been resolved as well (unused variables). For source files that were corrected in this way, the remaining warnings relate to equality comparisons for real values, unused dummy arguments in subroutines and potential 0 indices into arrays (in all cases, if statements prevented this from happening).
+
 ## [NJOY2016.69](https://github.com/njoy/NJOY2016/pull/281)
 This update fixes a number of minor issues:
   - PURR now writes Bondarenko data obtained from the probability tables to MF2 MT152 instead of the Bondarenko data obtained from the direct sampled cross sections (for very low dilutions, the Bondarenko data obtained using these two methods does not align, with the direct sampled data leading to extremely low P1 values). When comparing with the Bondarenko data at low dilutions obtained with UNRESR, the Bondarenko data obtained from the probability table directly seems to be the best.

src/acefc.f90

@@ -3894,7 +3894,7 @@ subroutine convr(nin,npend,nout,nscr,nedis,nethr,matd)
 !  imax=49
    imax=50
    imaxsq=imax*imax
-   lmax=100
+   lmax=imax*(imax-1)/2
    nned=50
    allocate(disc(nned))
    allocate(ee(imax))
@@ -18724,7 +18724,7 @@ subroutine aploxp(nout,iwcol,hk)
    integer::ipt,mtrh,nmtr
    integer hpd,lsigh,sigh,landh,andh,ldlwh,dlwh,yh
    real(kr)::test,e,xs,xstep,ystep,xtag,ytag,thin,xlast
-   real(kr)::ep,pd,ylast,break,cc,pp,rat,stepm,elast
+   real(kr)::ep,pd,epl,pdl,ylast,break,cc,pp,rat,stepm,elast
    real(kr)::xmin,xmax,ymin,ymax,zmin,zmax,zmax1,zmax2,heat
    character(10)::name
    character(1)::qu=''''
@@ -19141,13 +19141,15 @@ subroutine aploxp(nout,iwcol,hk)
                   intt=nint(xss(loci))
                   nn=nint(xss(loci+1))
                   loci=loci+1
+                  epl = xss(loci+1)
+                  pdl = xss(loci+nn+1)
                   do j=1,nn
                      ep=xss(loci+j)
                      pd=xss(loci+nn+j)
-                     if (pd.ge.zmin.and.pd.le.zmax) then
-                        if (ep.lt.xmin) xmin=ep
-                        if (ep.gt.xmax) xmax=ep
-                     endif
+                     if(pd  .ge. zmin .and. epl .lt. xmin) xmin = epl
+                     if(pdl .ge. zmin .and. ep  .gt. xmax) xmax = ep
+                     epl = ep
+                     pdl = pd
                   enddo
                endif
             enddo

src/acepn.f90

@@ -48,7 +48,7 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
    integer::nneut,nphot,nprot,ndeut,ntrit,nhe3,nhe4
    integer::k,ia,iaa,nk,ik,lly,izai,izap,law,jscr,nrr,npp
    integer::ll,lll,lep,ne,llh,lld,ie,np,ip,mtt,lct,ii
-   integer::icapt,jj,itype,it,jp,nr,il,llht,iie,lang,lleg,ileg
+   integer::icapt,jj,itype,it,jp,nr,il,llht,iie,lang
    integer::iint,nn,kk,m,intt,last,lf,jnt,ja,jb,ipp,irr
    integer::lee,lle,nd,na,ncyc,ng,ig,nnr,nnp,mf,mt
    integer::ipt,ntrp,pxs,phn,mtrp,tyrp,lsigp,sigp,landp,andp,ldlwp,dlwp
@@ -74,7 +74,6 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
    real(kr),parameter::eps=1.e-10_kr
    real(kr),parameter::zero=0
    real(kr),parameter::one=1
-   integer,parameter::ni=64
    character(66)::text
    emc2=amassn*amu*clight*clight/ev/emev
    tvn=1
@@ -1861,7 +1860,6 @@ subroutine phnfix(itype,nin,nout,ndir,iprint,nplot,mcnpx,suff,&
    real(kr),parameter::zero=0
 
    integer,parameter::ner=1
-   integer,parameter::nbw=1
 
    !--read type 1 ace format file
    call openz(nin,0)
@@ -2455,8 +2453,8 @@ subroutine phnout(itype,nout,ndir,mcnpx,hk,izn,awn)
    real(kr)::awn(16)
    character(70)::hk
    ! internals
-   integer::l,n,ne,ip,mftype,nr,li,ir,nn,ll,k,np,nw,nmu,nrr
-   integer::ii,lnw,law,kk,nern,lrec,j,i
+   integer::l,n,ne,ip,mftype,nr,nn,ll,k,np,nmu,nrr
+   integer::lnw,law,nern,lrec,j,i
    integer::ipt, ntrp, pxs, phn, mtrp, tyrp, lsigp, sigp, landp, andp, ldlwp, dlwp ! IXS
    integer::rlocator  ! locator index for reaction data
    integer::plocator  ! locator index for the particle IXS array

src/errorr.f90

@@ -796,6 +796,8 @@ subroutine errorr
             if (mmtres(i).eq.106) mmtres(i)=750
             if (mmtres(i).eq.107) mmtres(i)=800
          enddo
+      else
+         call desammy
       endif
       call findf(matd,2,0,nendf)
       deallocate(a)
@@ -1048,6 +1050,7 @@ subroutine errorr
    !--errorr is finished.
    call atend(nout,0)
   330 continue
+   if (nmtres.ne.0) call desammy
    if (allocated(flx))  deallocate(flx)
    if (allocated(sig))  deallocate(sig)
    if (allocated(cov))  deallocate(cov)
@@ -3225,7 +3228,6 @@ subroutine resprx(nwscr,a)
       enddo
    enddo
    call closz(nscr6)
-   if (nmtres.gt.0) call desammy
 
    return
    end subroutine resprx

src/heatr.f90

@@ -644,7 +644,13 @@ subroutine hinit(iold,nend4,nend6,nscr)
          ifc4=0
          ifc5=0
          ifc6=0
-         call listio(nendf,0,0,scr,nb,nw)
+         l=1
+         call listio(nendf,0,0,scr(l),nb,nw)
+         l=l+nw
+         do while (nb.ne.0)
+            call moreio(nendf,0,0,scr(l),nb,nw)
+            l=l+nw
+         enddo
          nply=nint(scr(4))
          if (lfc.eq.0) then
             allocate(cpoly(0:nply))
@@ -693,7 +699,7 @@ subroutine hinit(iold,nend4,nend6,nscr)
             qdel=0
             do i=1,nfc
                l=1
-               call tab1io(nendf,0,0,scr,nb,nw)
+               call tab1io(nendf,0,0,scr(l),nb,nw)
                do while (nb.ne.0)
                   l=l+nw
                   call moreio(nendf,0,0,scr(l),nb,nw)
@@ -810,9 +816,16 @@ subroutine hinit(iold,nend4,nend6,nscr)
             awrr=c2h+1
             ielem=mod(nint(c1h),1000)
             do ik=1,nk
-               call tab1io(nendf,nend6,0,scr,nb,nw)
+               l=1
+               call tab1io(nendf,nend6,0,scr(l),nb,nw)
+               l=l+nw
+               do while (nb.ne.0)
+                  call moreio(nendf,nend6,0,scr(l),nb,nw)
+                  l=l+nw
+               enddo
                nr=n1h
                zap=c1h
+               law=l2h
                if (nint(zap).eq.0) then
                   i6p=i6p+1
                   mt6yp(i6p)=mth
@@ -853,10 +866,6 @@ subroutine hinit(iold,nend4,nend6,nscr)
                   mt6no(ii6)=nk
                endif
                if (zap.eq.zero) mgam=10+mod(mgam,10)
-               law=l2h
-               do while (nb.ne.0)
-                  call moreio(nendf,nend6,0,scr,nb,nw)
-               enddo
                if (law.eq.6) then
                   call contio(nendf,nend6,0,scr,nb,nw)
                else if (law.eq.1.or.law.eq.2.or.law.eq.5) then
@@ -1362,7 +1371,7 @@ subroutine nheat(iold,inew,nscr,nend4,nend6,local)
    elst=elst-elst/10000000
   193 continue
    if (e.lt.thresh) go to 290
-   call gety1(e,enext,idis,y,nin,b)
+   call gety1(e,enext,idis,y,nin,b) ! y is the reaction cross section
    ! check for energy-dependent q
    if (iimt.gt.0) then
       if (qa(iimt).ge.qtest) then
@@ -2311,7 +2320,13 @@ subroutine hgtyld(e,enext,idis,yld,mat,mf,mt,itape,a,na)
       call contio(itape,0,0,a,nb,nw)
       lnu=nint(a(4))
       if (mth.eq.455) then
-         call listio(itape,0,0,a,nb,nw)
+         loc=1
+         call listio(itape,0,0,a(loc),nb,nw)
+         loc=loc+nw
+         do while (nb.ne.0)
+            call moreio(itape,0,0,a(loc),nb,nw)
+            loc=loc+nw
+         enddo
          lnd=nint(a(5))
          if (lnd.ne.6.and.lnd.ne.8)&
            call error('hgtyld','illegal lnd, must be 6 or 8',' ')
@@ -2347,8 +2362,15 @@ subroutine hgtyld(e,enext,idis,yld,mat,mf,mt,itape,a,na)
             nr=nint(a(5))
             enext=a(7+2*nr)
          else
-            call listio(itape,0,0,a,nb,nw)
+            loc=1
+            call listio(itape,0,0,a(loc),nb,nw)
+            loc=loc+nw
             na=nw
+            do while (nb.ne.0)
+               call moreio(itape,0,0,a(loc),nb,nw)
+               loc=loc+nw
+               na=na+nw
+            enddo
             enext=emax
          endif
       endif
@@ -2641,8 +2663,8 @@ subroutine tabbar(f,a,law)
       yh=a(ibase+ncyc*(i-1)+2)
       if (law.gt.0.and.i.gt.nbt) then
          ir=ir+1
-         nbt=nint(a(ibase+2*ir-1))
-         inn=nint(a(ibase+2*ir))
+         nbt=nint(a(6+2*ir-1))
+         inn=nint(a(6+2*ir))
       endif
       if (xl.ne.xh) then
 
@@ -2778,6 +2800,11 @@ subroutine sixbar(e,ebar,yld,dame,nin,c,ncmax,nscr,b,nbmax,&
    do while (ik.lt.irec-1)
       ik=ik+1
       call tab1io(nin,0,0,c(l),nb,nw)
+      l=l+nw
+      do while (nb.ne.0)
+         call moreio(nin,0,0,c(l),nb,nw)
+         l=l+nw
+      enddo
       law=l2h
       call skip6(nin,0,0,c(l),law)
    enddo
@@ -2790,10 +2817,14 @@ subroutine sixbar(e,ebar,yld,dame,nin,c,ncmax,nscr,b,nbmax,&
   110 continue
    l=1
    call tab1io(nin,0,0,c(l),nb,nw)
+   l=l+nw
+   do while (nb.ne.0)
+      call moreio(nin,0,0,c(l),nb,nw)
+      l=l+nw
+   enddo
    zap=c1h
    awp=c2h
    law=l2h
-   l=l+nw
 
    !--error for awp=0. for non photons
    if (zap.ne.0.and.awp.eq.0) then
@@ -2802,10 +2833,6 @@ subroutine sixbar(e,ebar,yld,dame,nin,c,ncmax,nscr,b,nbmax,&
      call error('sixbar',strng,' ')
    endif
 
-   do while (nb.ne.0)
-      call moreio(nin,0,0,c(l),nb,nw)
-      l=l+nw
-   enddo
    iflag=0
    disc102=0
    if (zap.eq.zero) then
@@ -2917,7 +2944,7 @@ subroutine sixbar(e,ebar,yld,dame,nin,c,ncmax,nscr,b,nbmax,&
       call getsix(elo,flo,dlo,c(iraw),law,lang,lep,irec)
    else
       ztt=int(zat/1000)
-      call tabsq6(flo,dlo,c(iraw),law,ztt,awrt)
+      call tabsq6(flo,dlo,c(iraw),law,ztt,awrt,elo,c(1))
    endif
    if ((mth.ge.18.and.mth.le.21).or.mth.eq.38) then
       matd=math
@@ -3001,7 +3028,7 @@ subroutine sixbar(e,ebar,yld,dame,nin,c,ncmax,nscr,b,nbmax,&
       call getsix(ehi,fhi,dhi,c(iraw),law,lang,lep,irec)
    else
       ztt=int(zat/1000)
-      call tabsq6(fhi,dhi,c(iraw),law,ztt,awrt)
+      call tabsq6(fhi,dhi,c(iraw),law,ztt,awrt,ehi,c(1))
    endif
    go to 305
 
@@ -3228,7 +3255,7 @@ subroutine getsix(e,ebar,dame,c,law,lang,lep,irec)
    fl=0
    do i=1,nep
       l=7+ncyc*(i-1)
-      xx=c(l)
+      xx=abs(c(l))
       yy=c(l+1)
       x2=xx
       if (irec.eq.0) then
@@ -4105,18 +4132,18 @@ real(kr) function h6psp(ep,epnext,epmax,w,e,c)
    return
    end function h6psp
 
-   subroutine tabsq6(g,h,a,law,z,awr)
+   subroutine tabsq6(g,h,a,law,z,awr,e,yld)
    !-------------------------------------------------------------------
    ! Compute average of photon recoil energy from capture and
    ! corresponding damage energy for File 6 capture photons
    !-------------------------------------------------------------------
    use endf ! provides terp1
    ! externals
    integer::law
-   real(kr)::g,h,a(*),z,awr
+   real(kr)::g,h,a(*),z,awr,e,yld(*)
    ! internals
-   integer::nd,np,ncyc,ibase,inn,nc,i,j
-   real(kr)::ein,rein,x,y,xr,awc,xh,yh,xl,yl,dx,s
+   integer::nd,np,ncyc,ibase,inn,nc,i,j,ip,ir,idis
+   real(kr)::ein,rein,x,y,xr,awc,xh,yh,xl,yl,dx,s,yield,enext
    integer,parameter::nq=4
    real(kr),dimension(nq),parameter::qp=(/&
      -.86114e0_kr,-.33998e0_kr,.33998e0_kr,.86114e0_kr/)
@@ -4135,6 +4162,11 @@ subroutine tabsq6(g,h,a,law,z,awr)
    ibase=6
    inn=2
 
+   !--interpolate the photon yield
+   ip=2
+   ir=1
+   call terpa(yield,e,enext,idis,yld(1),ip,ir)
+
    !--accumulate contributions from discrete levels
    if (nd.ne.0) then
       do i=1,nd
@@ -4174,8 +4206,8 @@ subroutine tabsq6(g,h,a,law,z,awr)
    endif
 
    !--finished
-   g=g/s
-   h=h/s
+   g=yield*g/s
+   h=yield*h/s
    return
    end subroutine tabsq6
 
@@ -4721,7 +4753,13 @@ subroutine hconvr(nin,nout,nscr)
    lg=l2h
    g=1
    l2flg=1
-   call listio(nin,0,0,scr,nb,nw)
+   l=1
+   call listio(nin,0,0,scr(l),nb,nw)
+   l=l+nw
+   do while (nb.ne.0)
+      call moreio(nin,0,0,scr(l),nb,nw)
+      l=l+nw
+   enddo
 
    !--set base value for mt0
    if (mth.ge.51.and.mth.le.91.and.mt0.ne.49) mt0=49
@@ -5610,8 +5648,8 @@ subroutine tabsqr(g,h,a,law,z,awr)
       yh=a(ibase+2*i)
       if (i.gt.nbt) then
          ir=ir+1
-         nbt=nint(a(ibase+2*ir-1))
-         inn=nint(a(ibase+2*ir))
+         nbt=nint(a(6+2*ir-1))
+         inn=nint(a(6+2*ir))
       endif
       if (xl.ne.xh) then
          dx=xh-xl

src/reconr.f90

@@ -289,6 +289,8 @@ subroutine reconr
                     if (mmtres(i).eq.106) mmtres(i)=750
                     if (mmtres(i).eq.107) mmtres(i)=800
                  enddo
+              else
+                 call desammy
               endif
               call repoz(nin)
               call findf(mata,2,0,nin)