NJOY2016.69

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.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.
+  - MF6 LAW=2 represents discrete two body scattering in which only angular distribution data is given (knowing that the outgoing energy of the secondary particle can be determined through kinematics when the angle is known). When calculating heating numbers based on LAW=2, ACER assumes that the yield of the secondary particle is 1, which is correct in all cases except when MT5 is used as a lumped reaction. Heating numbers in ACER for photonuclear files using LAW=2 in an MT5 entry are now correctly multiplied by the yield. A warning message is printed out whenever this situation is detected. Test 78 was added as part of this correction.
+  - Previously, ERRORR would segfault for LRF=7 resonance evaluations when MF33 was present without MF32.  A check for this situation now avoids this.
+  - Fixed an issue in GROUPR when reading some of the FENDL3.2 evaluations.
+
+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.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.
@@ -63,8 +72,7 @@ This release addresses issue [\#201](https://github.com/njoy/NJOY2016/issues/201
 ## [NJOY2016.63](https://github.com/njoy/NJOY2016/pull/193)
 This fixes a bug in ERRORR when using the `999` option. When using this option, the input and output tapes are not closed in Fortran. This causes problems in NJOY21 as the output from ERRORR isn't completely written to disk before the next module starts. This update simply closes files `nitape` and `notape` which resolves the issue.
 
-In addition, some logic statements in the GROUPR conver and ACER convr subroutines have been corrected to fix failures in the test
-suite when building in Debug mode.
+In addition, some logic statements in the GROUPR conver and ACER convr subroutines have been corrected to fix failures in the test suite when building in Debug mode.
 
 ## [NJOY2016.62](https://github.com/njoy/NJOY2016/pull/191)
 This adds a number of changes to NJOY2016 contributed by Toshihiko Kawano, Bob McFarlane, IAEA and CIEMAT. In particular, the following changes were made:

src/acepn.f90

@@ -52,14 +52,14 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
    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
-   integer::izarec,nl,iil,nexn,nexd,ki
+   integer::izarec,nl,iil,nexn,nexd,ki,ilaw2mt5
    integer::nle
    integer::imu,intmu,nmu
    real(kr)::emc2,e,enext,s,y,ynext,heat,en,ep,g,h,epl
    real(kr)::tneut,tphot,tprot,tdeut,ttrit,the3,the4,thresh
    real(kr)::ss,tt,ubar,sum,renorm,ebar,hh,u,theta,x,anorm
    real(kr)::ee,amass,avadd,avlab,avll,test,rkal,akal
-   real(kr)::eavi,avl,avcm,sign,dele,avav,zaid,gl,awp,awr,q
+   real(kr)::eavi,avl,avcm,sign,dele,avav,zaid,gl,awp,awr,q,yld
    real(kr)::av,del
    real(kr)::awprec,awpp
    integer,parameter::mmax=1000
@@ -78,6 +78,8 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
    character(66)::text
    emc2=amassn*amu*clight*clight/ev/emev
    tvn=1
+   ilaw2mt5=0
+   yld=1.
 
    nxsd=0
    jxsd=0
@@ -1095,8 +1097,8 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
                              sigfig(renorm*xss(nex+1+2*n+ii),9,0)  ! CDF(ii)
                         enddo
                         nex=nex+2+3*n       ! index for the next distribution
-                        e=xss(ie+iie)
-                        scr(llht+6+2*iie)=e
+                        e=xss(ie+iie)       ! e from xss, thus MeV
+                        scr(llht+6+2*iie)=e ! store it in scr, so energy in MeV
                         scr(llht+7+2*iie)=(awr-awpp)*(e+q)/awr
                      enddo
                      ! add in contribution to heating
@@ -1105,13 +1107,32 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
                      do ie=it,nes
                         e=xss(esz+ie-1)/emev
                         call terpa(h,e,en,idis,scr(llht),npp,nrr)
+                        !--get the yield and modify heating
+                        !--beware: e is already in MeV, scr(llht) has been
+                        !--modified to be in MeV already (see above)
+                        call terpa(yld,e*emev,en,idis,scr,npp,nrr)
+                        h=yld*h ! just in case, multiply by yield
+                        if (ilaw2mt5.eq.0.and.mt.eq.5.and.yld.ne.one) then
+                           ilaw2mt5=1
+                           write(text,'(''yield different from 1 for outgoing particle '',i5,''.'')')ip
+                           if (izarec.eq.ip) then
+                              call mess('acephn',&
+                                        'law=2 (discrete two-body scattering) used in mt=5 recoil',&
+                                        text)
+                           else
+                              call mess('acephn',&
+                                        'law=2 (discrete two-body scattering) used in mt=5',&
+                                        text)
+                           endif
+                        endif
                         ss=0
                         if (ie.ge.iaa) ss=xss(2+k+ie-iaa)
                         xss(phn+2+ie-it)=xss(phn+2+ie-it)+h*ss
                         if (xss(tot+ie-1).ne.zero)&
                           xss(thn+ie-1)=xss(thn+ie-1)&
                           +h*ss/xss(tot+ie-1)
                      enddo
+                     ilaw2mt5=0
                   else
                      call skip6(nin,0,0,scr,law)
                   endif

src/errorr.f90

@@ -3723,7 +3723,7 @@ subroutine rpxlc0(nwscr,a)
    integer::nwscr
    real(kr)::a(nwscr)
    ! internals
-   integer::nl,iloc,nrs,i,nr,ig,j,loop,igind,inow,id,ig1,ig2,ij,ii,k
+   integer::nl,iloc,nrs,i,nr,ig,j,loop,igind,inow,id,ig1,ig2,ij,ii
    integer::nb,nw,jd
    integer,parameter::maxnls=10
    integer,parameter::maxe=600000
@@ -4098,7 +4098,7 @@ subroutine rpxlc12(nwscr,a,iest,ieed)
    integer::nwscr,iest,ieed
    real(kr)::a(nwscr)
    ! internals
-   integer::nb,nw,lru1,lrf1,lb,i,ig,ig2,il,ii2,ii1,ipp,j,k
+   integer::nb,nw,lru1,lrf1,lb,i,ig,ig2,il,ii2,ii1,ipp,j
    integer::loopm,loopn,loop,ns,nrs1,nsmax,ipos,l1,l2,l3
    integer::lb2,lldum,igind,ipara,itmp,ilnum,ind,ii
    integer::imess,inow,jj,lll,nrs
@@ -4773,7 +4773,7 @@ subroutine rpxunr(a,amur,mxlru2,iest,ieed,nwscr)
    integer::mxlru2,iest,ieed,nwscr
    real(kr)::a(nwscr),amur(3,mxlru2)
    ! internals
-   integer::nb,nw,l1,l2,l3,nl,ig,i,j,k,ig2,ii,igind,njs,inow
+   integer::nb,nw,l1,l2,l3,nl,ig,i,j,ig2,ii,igind,njs,inow
    integer::loop,loopn,l0
    integer,parameter::maxb=4000
    integer,parameter::mxnpar=100
@@ -7439,7 +7439,7 @@ subroutine covout
   390 continue
 
    ! add contribution from resonance-parameter uncertainty
-   if (mfcov.ne.34.and.mfcov.ne.35) then
+   if (mfcov.eq.33.and.mf32.ne.0) then
       call rescon(ix,ixp,csig,cova,izero,ngn,nmt1)
    endif
 

src/groupr.f90

@@ -6430,7 +6430,7 @@ subroutine getyld(e,enext,idis,yld,mat,mf,mt,lfs,itape)
    integer::idis,mat,mf,mt,lfs,itape,i,nr,nc
    real(kr)::e,enext,yld,term
    ! internals
-   integer::mft,nb,nw,nk,ik,lnu,iza,lnd,izn,lfn,ip,ir,na,loc
+   integer::mft,nb,nw,nk,ik,lnu,iza,lnd,izn,lfn,ip,ir,na,loc,ll
    integer::ntmp
    real(kr),dimension(:),allocatable::tmp
    real(kr),parameter::emax=1.e10_kr
@@ -6459,12 +6459,23 @@ subroutine getyld(e,enext,idis,yld,mat,mf,mt,lfs,itape)
    if (mt.eq.455) go to 110
    lnd=0
    if (lnu.ne.1) go to 130
-   call listio(itape,0,0,tmp(loc),nb,nw)
+   ll=loc
+   call listio(itape,0,0,tmp(ll),nb,nw)
    na=nw
+   do while (nb.ne.0)
+      ll=ll+nw
+      call moreio(itape,0,0,tmp(ll),nb,nw)
+      na=na+nw
+   enddo
    enext=emax
    go to 190
   110 continue
-   call listio(itape,0,0,tmp(loc),nb,nw)
+   ll=loc
+   call listio(itape,0,0,tmp(ll),nb,nw)
+   do while (nb.ne.0)
+      ll=ll+nw
+      call moreio(itape,0,0,tmp(ll),nb,nw)
+   enddo
    lnd=nint(tmp(loc+4))
    if (lnd.gt.8) call error('getyld','illegal lnd.',' ')
    do i=1,lnd
@@ -7441,7 +7452,7 @@ subroutine getmf6(ans,ed,enext,idisc,yld,eg,ng,nl,iglo,ng2,nq,&
    integer,dimension(nssm)::iyss,izss,jjss
    integer,dimension(maxss)::jloss
    real(kr),dimension(:),allocatable::tmp
-   integer,parameter::ncmax=350
+   integer,parameter::ncmax=1000
    real(kr),parameter::emax=1.e10_kr
    real(kr),parameter::small=1.e-10_kr
    real(kr),parameter::shade=1.1999e0_kr
@@ -9576,12 +9587,12 @@ subroutine getfle(e,enext,idis,fle,nle,lcd,matd,mfd,mtd,nin)
    real(kr)::e,enext,fle(*)
    ! internals
    integer::iso,lidp,ltt,iraw,ir,nne,ne,int,nlo,nhi,ltt3,lttn
-   integer::nb,nwc,lvt,i,nlmax,il
+   integer::nb,nwc,lvt,i,nlmax,il,ll
    real(kr)::elo,ehi
    character(60)::strng
    integer,parameter::mxlg=65
    real(kr)::flo(mxlg),fhi(mxlg)
-   integer,parameter::ncmax=350
+   integer,parameter::ncmax=1000
    real(kr)::fls(ncmax)
    real(kr),parameter::emax=1.e10_kr
    real(kr),parameter::small=1.e-10_kr
@@ -9639,14 +9650,24 @@ subroutine getfle(e,enext,idis,fle,nle,lcd,matd,mfd,mtd,nin)
       ! read in raw data at first two incident energies.
       ! retrieve or compute legendre coefficients.
       iraw=1+nwc
-      if (ltt.eq.1) call listio(nin,0,0,fls(iraw),nb,nwc)
-      if (ltt.eq.2) call tab1io(nin,0,0,fls(iraw),nb,nwc)
+      ll=iraw
+      if (ltt.eq.1) call listio(nin,0,0,fls(ll),nb,nwc)
+      if (ltt.eq.2) call tab1io(nin,0,0,fls(ll),nb,nwc)
+      do while (nb.ne.0)
+         ll=ll+nwc
+         call moreio(nin,0,0,fls(ll),nb,nwc)
+      enddo
       elo=fls(iraw+1)
       nlo=nle
       if (lidp.ge.0) fls(iraw+3)=lidp
       call getco(flo,nlo,lcd,fls(iraw),lct,ltt,idis)
-      if (ltt.eq.1) call listio(nin,0,0,fls(iraw),nb,nwc)
-      if (ltt.eq.2) call tab1io(nin,0,0,fls(iraw),nb,nwc)
+      ll=iraw
+      if (ltt.eq.1) call listio(nin,0,0,fls(ll),nb,nwc)
+      if (ltt.eq.2) call tab1io(nin,0,0,fls(ll),nb,nwc)
+      do while (nb.ne.0)
+         ll=ll+nwc
+         call moreio(nin,0,0,fls(ll),nb,nwc)
+      enddo
       ehi=fls(iraw+1)
       nhi=nle
       if (lidp.ge.0) fls(iraw+3)=lidp
@@ -9691,8 +9712,13 @@ subroutine getfle(e,enext,idis,fle,nle,lcd,matd,mfd,mtd,nin)
    else if (nne.eq.ne) then
       call error('getfle','desired energy above highest given.',' ')
    endif
-   if (ltt.eq.1) call listio(nin,0,0,fls(iraw),nb,nwc)
-   if (ltt.eq.2) call tab1io(nin,0,0,fls(iraw),nb,nwc)
+   ll=iraw
+   if (ltt.eq.1) call listio(nin,0,0,fls(ll),nb,nwc)
+   if (ltt.eq.2) call tab1io(nin,0,0,fls(ll),nb,nwc)
+   do while (nb.ne.0)
+      ll=ll+nwc
+      call moreio(nin,0,0,fls(ll),nb,nwc)
+   enddo
    ehi=fls(iraw+1)
    nhi=nle
    if (lidp.ge.0) fls(iraw+3)=lidp
@@ -9788,6 +9814,9 @@ subroutine getaed(aed,e,enext,idis,nl,nq,mat,mf,mt,nin,nlg)
       call tab1io(nin,0,0,p,nb,nw)
       itt=-l1h
       law=l2h
+      do while (nb.ne.0)
+         call moreio(nin,0,0,p,nb,nw)
+      enddo
 
       !--incoherent elastic or inelastic with E-E'-mu ordering
       if (law.eq.1) then
@@ -10806,7 +10835,7 @@ subroutine conver(nin,nout,nscr)
    integer::mt0,mt0old,nb,nw,jzar,lg,n,kk,k,kp1,l,lm1,mtl,no455
    integer::nnu,lnu,mf,mt,l1,l2,n1,n2,nk,jzap,ne,ltp,nm
    integer::m1,m2,mtnow,mttst,nn
-   integer::ltt,lcd,nl
+   integer::ltt,lcd,nl,ll
    real(kr)::g,ei,eja,ejb,e1,enext,p,za2,za,ysum,yy,tsave,w,elow,etop,enxt,x
    character(60)::strng
    integer::ngam(440)
@@ -11125,7 +11154,12 @@ subroutine conver(nin,nout,nscr)
    lg=l2h
    g=1
    l2flg=1
-   call listio(nin,0,0,scr,nb,nw)
+   ll=1
+   call listio(nin,0,0,scr(ll),nb,nw)
+   do while (nb.ne.0)
+      ll=ll+nw
+      call moreio(nin,0,0,fls(ll),nb,nw)
+   enddo
 
    !--make sure mt0 is correct for this range of mt's.
    if (mth.ge.51.and.mth.le.90.and.mt0.ne.49) mt0=49
@@ -11441,10 +11475,20 @@ subroutine conver(nin,nout,nscr)
       spi=c1h
       ne=n2h
       do ie=1,ne
-         call listio(nin,0,0,scr,nb,nw)
+         ll=1
+         call listio(nin,0,0,scr(ll),nb,nw)
          ltp=l1h
+         do while (nb.ne.0)
+            ll=ll+nw
+            call moreio(nin,0,0,scr(ll),nb,nw)
+         enddo
          if (ltp.le.2) then
-            call listio(0,nout,nscr,scr,nb,nw)
+            ll=1
+            call listio(0,nout,nscr,scr(ll),nb,nw)
+            do while (nb.ne.0)
+               ll=ll+nw
+               call moreio(0,nout,nscr,scr(ll),nb,nw)
+            enddo
          else
             e1=c2h
             call getsig(e1,enext,idis,sig,1,1)
@@ -11464,7 +11508,12 @@ subroutine conver(nin,nout,nscr)
             do il=1,nl
                scr(6+il)=fl(il)
             enddo
-            call listio(0,nout,nscr,scr,nb,nw)
+            ll=1
+            call listio(0,nout,nscr,scr(ll),nb,nw)
+            do while (nb.ne.0)
+               ll=ll+nw
+               call moreio(0,nout,nscr,scr(ll),nb,nw)
+            enddo
          endif
       enddo
       izap=-izap
@@ -11586,7 +11635,12 @@ subroutine conver(nin,nout,nscr)
       ik=0
       do while (ik.lt.nk)
          ik=ik+1
-         call listio(nin,0,0,scr,nb,nw)
+         ll=1
+         call listio(nin,0,0,scr(ll),nb,nw)
+         do while (nb.ne.0)
+            ll=ll+nw
+            call moreio(nin,0,0,scr(ll),nb,nw)
+         enddo
          izan=nint(c1h)
          imf=l1h
          iis=l2h

src/purr.f90

@@ -37,7 +37,7 @@ module purm
    real(kr),dimension(:),allocatable::sb
 
    ! probability table globals
-   real(kr),dimension(:,:,:),allocatable::bval,sigb
+   real(kr),dimension(:,:,:),allocatable::bval
    real(kr),dimension(:),allocatable::tmin,tmax,tsum
 
    ! storage array for unresolved resonance parameters
@@ -643,7 +643,6 @@ subroutine purr
    deallocate(sigpl)
    if (allocated(bval)) then
       deallocate(bval)
-      deallocate(sigb)
       deallocate(tmin)
       deallocate(tmax)
       deallocate(tsum)
@@ -1856,13 +1855,11 @@ subroutine unrest(bkg,sig0,sigf,tabl,tval,er,gnr,gfr,ggr,gxr,gt,&
       elow=10
       if (allocated(bval)) then
          deallocate(bval)
-         deallocate(sigb)
          deallocate(tmin)
          deallocate(tmax)
          deallocate(tsum)
       endif
       allocate(bval(8,nsig0,ntemp))
-      allocate(sigb(5,nsig0,ntemp))
       allocate(tmin(ntemp))
       allocate(tmax(ntemp))
       allocate(tsum(ntemp))
@@ -2493,38 +2490,22 @@ subroutine unrest(bkg,sig0,sigf,tabl,tval,er,gnr,gfr,ggr,gxr,gt,&
                bval(7,i,itemp)=bval(7,i,itemp)+ttt*den*den
             endif
          enddo
-         sigb(1,i,itemp)=bval(1,i,itemp)/bval(6,i,itemp)
-         sigb(2,i,itemp)=bval(2,i,itemp)/bval(6,i,itemp)
-         sigb(3,i,itemp)=bval(3,i,itemp)/bval(6,i,itemp)
-         sigb(4,i,itemp)=bval(4,i,itemp)/bval(6,i,itemp)
-         sigb(5,i,itemp)=bval(5,i,itemp)/bval(7,i,itemp)
+         sigf(1,i,itemp)=bval(1,i,itemp)/bval(6,i,itemp)
+         sigf(2,i,itemp)=bval(2,i,itemp)/bval(6,i,itemp)
+         sigf(3,i,itemp)=bval(3,i,itemp)/bval(6,i,itemp)
+         sigf(4,i,itemp)=bval(4,i,itemp)/bval(6,i,itemp)
+         sigf(5,i,itemp)=bval(5,i,itemp)/bval(7,i,itemp)
       enddo
    enddo
    if (iprint.gt.0) then
       do itemp=1,ntemp
          do i=1,nsig0
             write(nsyso,'(3x,1p,2e10.3,5e12.4)')&
-              temp(itemp),sig0(i),(sigb(j,i,itemp),j=1,5)
+              temp(itemp),sig0(i),(sigf(j,i,itemp),j=1,5)
          enddo
       enddo
    endif
 
-   !--NOTE
-   !--sigf contains direct calculations of self shielded cross sections.
-   !--sigb contains self shielded cross sections computed from the
-   !--probability table.  copy sigb to sigf in order to make the MF152
-   !--values more consistent with the MF153 probability tables (even
-   !--if slightly less accurate).
-   !
-   ! do itemp=1,ntemp
-   !    do i=1,nsig0
-   !       do j=1,5
-   !!         sigf(i,j,itemp)=sigb(i,j,itemp)  !(i,j) wrong?!?!?
-   !          sigf(j,i,itemp)=sigb(j,i,itemp)  !(j,i) right? ... check w/Bob
-   !       enddo
-   !    enddo
-   ! enddo
-
    !--renormalize probability table and bondarenko
    !--cross sections to the computed infinitely-dilute values
    do itemp=1,ntemp