Merge pull request #191 from njoy/develop

NJOY2016.62
njoy · Mar 16, 2021 · d157060 · d157060
2 parents ce8ff6a + 9f08877
commit d157060
Show file tree

Hide file tree

Showing 11 changed files with 205,543 additions and 79 deletions.
diff --git a/.github/workflows/ContinuousIntegration.yml b/.github/workflows/ContinuousIntegration.yml
@@ -4,9 +4,11 @@ name: Continuous Integration
 # This workflow is triggered on pushes and pull requests to the repository.
 on:
   push:
-    branches: '*'
+    branches: '**'
   pull_request:
-    branches: 'master'
+    branches:
+      - 'master'
+      - 'develop'
 
 jobs:
   build:

diff --git a/ReleaseNotes.md b/ReleaseNotes.md
@@ -1,6 +1,14 @@
 # Release Notes&mdash;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.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:
+  - fixed an issue in PURR when nunx is not set to the default value (when nunx is set to anything else but 0 or the number of energy points in the unresolved energy region, purr will now properly calculate the probability tables for the given number of tables spread out over the entire unresolved resonance region)
+  - fixed an issue in RECONR where the wrong channel radius was used for the calculation of the phase shift
+  - array sizes were increased
+  - additional allocated arrays are now deallocated at the end of the subroutine in which they were created
+  - added consistency checking for probability tables in ACER for continuous energy libraries
+
 ## [NJOY2016.61](https://github.com/njoy/NJOY2016/pull/140)
 This updates the physical constants used in NJOY2016 to the CODATA2018 values recommended in the [ENDF-6 Formats Manual](https://www.nndc.bnl.gov/csewg/docs/endf-manual.pdf). The manual was updated through a format proposal at CSEWG and merged in the official repository in commit [28aca81c](https://git.nndc.bnl.gov/endf/format/endf6man/-/commit/28aca81c23d96ed783efde7f3759477972040efb).
 

diff --git a/docs/testDescription.md b/docs/testDescription.md
@@ -288,3 +288,9 @@ This test was added to following issue [\#163](https://github.com/njoy/NJOY2016/
 ## Test Problem 62
 
 This test was added to following issue [\#173](https://github.com/njoy/NJOY2016/issues/173). For a new d+He3 evaluation, the ACE file produced by NJOY2016 still has a few NaN values appearing in it. This problem was due to an array index overflow in acecpe.
+
+## Test Problem 63
+
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/63/input)]
+
+Tests 63 was added as a consequence of issue [\#178](https://github.com/njoy/NJOY2016/issues/178). It verifies that setting nunx in PURR to anything other than the default value does not break downstream processing (in this case up to ACER and VIEWR). The input file is equivalent to the input file for test 35 (with the exception of nunx which is set to 2).
diff --git a/src/acefc.f90 b/src/acefc.f90
@@ -840,6 +840,7 @@ subroutine tabize(mti,ln,nin,nout,itape,err,b)
          ncs(nxc)=ncs(nxc)+1+nint((scr(5)+5)/6)
       endif
       call tosend(nin,nout,0,scr)
+      deallocate(scr)
       return
    endif
 
@@ -900,6 +901,7 @@ subroutine tabize(mti,ln,nin,nout,itape,err,b)
          deallocate(p1456nu)
       endif
       call contio(nin,nout,0,scr,nb,nw)
+      deallocate(scr)
       return
    endif
 
@@ -1028,6 +1030,11 @@ subroutine tabize(mti,ln,nin,nout,itape,err,b)
    endif
    call contio(nin,nout,0,scr,nb,nw)
    call closz(-itape)
+   deallocate(y)
+   deallocate(x)
+   deallocate(buf)
+   deallocate(scr2)
+   deallocate(scr)
 
    !--finished
    return
@@ -1939,7 +1946,7 @@ subroutine topfil(nin,nout,matd,newfor)
    real(kr)::b(50)
    real(kr),dimension(:),allocatable::tab1
    real(kr),dimension(:),allocatable::scr
-   integer,parameter::nwmaxn=65000
+   integer,parameter::nwmaxn=1000000
    real(kr),parameter::big=1.e9_kr
    real(kr),parameter::zero=0
    real(kr),parameter::one=1
@@ -1953,7 +1960,7 @@ subroutine topfil(nin,nout,matd,newfor)
    jpp=0
    jpnu=0
    nkk=0
-   nt1w=4200
+   nt1w=20000
    allocate(tab1(nt1w))
    nsix=19
    nin0=0
@@ -2517,6 +2524,7 @@ subroutine topfil(nin,nout,matd,newfor)
 
    !--topfil is finished.
    deallocate(scr)
+   deallocate(tab1)
    return
    end subroutine topfil
 
@@ -3825,6 +3833,10 @@ subroutine gamsum(npend,nout,nf12c,matd,iopp)
    call closz(-iold)
    call closz(-inew)
    call repoz(nf12c)
+   deallocate(scr2)
+   deallocate(scr)
+   deallocate(bufn)
+   deallocate(bufo)
    return
    end subroutine gamsum
 
@@ -4457,6 +4469,15 @@ subroutine convr(nin,npend,nout,nscr,nedis,nethr,matd)
      &'' energy discontinuities found in gamma files''/&
      &(10x,1p,e13.5))')&
      (disc(iedis),iedis=1,nedis)
+   deallocate(tot)
+   deallocate(yold)
+   deallocate(scr)
+   deallocate(rr)
+   deallocate(aa)
+   deallocate(yy)
+   deallocate(es)
+   deallocate(eg)
+   deallocate(ee)
    return
    end subroutine convr
 
@@ -4836,6 +4857,10 @@ subroutine gamout(ngend,nendf,nout,nf12c,matd)
    !--the gamma angle and energy distributions are ready
    call amend(nout,0)
    call atend(nout,0)
+   if (allocated(ee)) deallocate(ee)
+   if (allocated(eb)) deallocate(eb)
+   deallocate(scr)
+   deallocate(sig)
    deallocate(egn)
    deallocate(egg)
    return
@@ -4878,7 +4903,7 @@ subroutine acelod(nin,suff,matd,tempd,newfor,mcnpx,ismooth)
    !--initialize
    inow=0
    nnex=0
-   nwscr=10000
+   nwscr=1000000
    allocate(scr(nwscr))
    do i=1,8
       nxsd(i)=0
@@ -6233,6 +6258,11 @@ subroutine acelod(nin,suff,matd,tempd,newfor,mcnpx,ismooth)
 
    !--load particle production information
    call acelcp(next,matd,nin,za,awr)
+   if (allocated(urd)) deallocate(urd)
+   if (allocated(nudn)) deallocate(nudn)
+   if (allocated(nup)) deallocate(nup)
+   if (allocated(nut)) deallocate(nut)
+   deallocate(scr)
 
    return
    end subroutine acelod
@@ -7011,7 +7041,7 @@ subroutine acelf6(next,i,matd,mt,q,iza,izai,nin,newfor,ismooth)
    integer::loc(5)
    character(60)::strng
    real(kr),dimension(:),allocatable::scr
-   integer,parameter::nwscr=18000
+   integer,parameter::nwscr=1000000
    real(kr),parameter::emev=1.e6_kr
    real(kr),parameter::small=1.e-30_kr
    real(kr),parameter::etop=1.e10_kr
@@ -8075,7 +8105,7 @@ subroutine acelpp(next,matd,ngmt,nin)
    real(kr),dimension(:),allocatable::dise
    real(kr),dimension(:),allocatable::tdise
    character(66)::strng
-   integer,parameter::nwscr=5000
+   integer,parameter::nwscr=1000000
    integer,parameter::ndise=5000
    real(kr),dimension(:),allocatable::phot
    real(kr),parameter::emev=1.e6_kr
@@ -9007,7 +9037,7 @@ subroutine acelcp(next,matd,nin,za,awr)
    real(kr)::eav,suml,chkl,ebar,ad,amuu,amulst,chklst,heat,g
    real(kr)::epl,gl,gammsq,amun,eavlst
    real(kr),dimension(:),allocatable::scr
-   integer,parameter::nwscr=800000
+   integer,parameter::nwscr=10000000
    real(kr),parameter::awr1=.99862e0_kr
    real(kr),parameter::awr2=1.99626e0_kr
    real(kr),parameter::awr3=2.98960e0_kr
@@ -9075,6 +9105,7 @@ subroutine acelcp(next,matd,nin,za,awr)
       ptype=0
       ntro=0
       ploct=0
+      deallocate(scr)
       return
    endif
 
@@ -10704,7 +10735,10 @@ subroutine acelcp(next,matd,nin,za,awr)
    !--continue loop over production types
    enddo
    len2=next-1
-   if (izai.le.1) return
+   if (izai.le.1) then
+     deallocate(scr)
+     return
+   endif
 
    !--for incident charged particles,
    !--go back through file 6 and get heating from recoils
@@ -12338,6 +12372,7 @@ subroutine aceppp(izai,nbina)
         &'' ---------------------------------------------------''//&
         &(6x,12i6))') (nint(xss(i+yp)),i=1,nyp)
    endif
+   if (allocated(indx)) deallocate(indx)
    return
    end subroutine aceppp
 
@@ -14117,6 +14152,8 @@ subroutine consis
    integer::ipt,mtrh,nmtr
    integer::hpd,lsigh,sigh,ldlwh,dlwh,nj
    integer::locl(20)
+   integer::nure,nurb,lurt,luri,lura,lurf,ib
+   integer::nerrtb,nerrxt,nerrxe,nerrxf,nerrxg,nerrxh
    real(kr)::thresh,elast,e,aprime,q,epmax,clast,ep,c,p,r
    real(kr)::cclast,co,cc,gsum,ss,y,gg,sum,frac,x,xlast
    character(10)::name
@@ -14637,6 +14674,71 @@ subroutine consis
       enddo
    endif
 
+   !--check unresolved-range probability tables
+   if (iurpt.ne.0) then
+      write(nsyso,'(/'' check probability tables''/&
+      &'' ------------------------'')')
+      nure=nint(xss(iurpt))
+      nurb=nint(xss(iurpt+1))
+      lurt=nint(xss(iurpt+2))
+      luri=nint(xss(iurpt+3))
+      lura=nint(xss(iurpt+4))
+      lurf=nint(xss(iurpt+5))
+      write(nsyso,'(/''     number of energies: '',i6/&
+                    &''     number of bins:     '',i6/&
+                    &''     interpolation law:  '',i6/&
+                    &''     inelastic reaction: '',i6/&
+                    &''     absorption reaction:'',i6)')&
+                    & nure,nurb,lurt,luri,lura
+      if (lurf.eq.0) write(nsyso,'(&
+                    &''     tables are cross sections (lssf=0)'')')
+      if (lurf.eq.1) write(nsyso,'(&
+                    &''     tables are factors (lssf=1)'')')
+      do ie=1,nure
+         write(nsyso,'(/'' energy='',1p,e14.6)') xss(iurpt+5+ie)
+         nerrtb=0
+         nerrxt=0
+         nerrxe=0
+         nerrxf=0
+         nerrxg=0
+         nerrxh=0
+         ll=iurpt+5+nure+(ie-1)*6*nurb
+         do ib=1,nurb
+           if (xss(ll+ib).lt.0.0d0) then
+             write(nsyso,'('' consis: negative probability value='',&
+             & 1pe14.6,'' for bin '',i4)')xss(ll+ib),ib
+             nerrtb=nerrtb+1
+           endif
+           if (xss(ll+nurb+ib).lt.0.0d0) then
+             write(nsyso,'('' consis: negative total cross section/factor='',&
+             & 1pe14.6,'' for bin '',i4)')xss(ll+nurb+ib),ib
+             nerrxt=nerrxt+1
+           endif
+           if (xss(ll+2*nurb+ib).lt.0.0d0) then
+             write(nsyso,'('' consis: negative elastic cross section/factor='',&
+             & 1pe14.6,'' for bin '',i4)')xss(ll+2*nurb+ib),ib
+             nerrxe=nerrxe+1
+           endif
+           if (xss(ll+3*nurb+ib).lt.0.0d0) then
+             write(nsyso,'('' consis: negative fission cross section/factor='',&
+             & 1pe14.6,'' for bin '',i4)')xss(ll+3*nurb+ib),ib
+             nerrxf=nerrxf+1
+           endif
+           if (xss(ll+4*nurb+ib).lt.0.0d0) then
+             write(nsyso,'('' consis: negative capture cross section/factor='',&
+             & 1pe14.6,'' for bin '',i4)')xss(ll+4*nurb+ib),ib
+             nerrxg=nerrxg+1
+           endif
+           if (xss(ll+5*nurb+ib).lt.0.0d0) then
+             write(nsyso,'('' consis: negative heating cross section/factor='',&
+             & 1pe14.6,'' for bin '',i4)')xss(ll+5*nurb+ib),ib
+             nerrxh=nerrxh+1
+           endif
+         enddo
+         nerr=nerr+nerrtb+nerrxt+nerrxe+nerrxf+nerrxg+nerrxh
+      enddo
+   endif
+
    !--check delayed neutron data
    if (ndnf.gt.0) then
       write(nsyso,'(/'' check delayed neutron fractions'')')