Merge pull request #308 from njoy/develop

NJOY2016.72

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.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.
+  - Fixed issues in acer to properly print already existing dosimetry and photoatomic ace files when running a stand-alone acer iopt=7 job.
+  - The meaning of legord and ifissp in the ERRORR input file has been repurposed when mfcov=34. The values now represent the L,L1 values of the MF34 sub-subsection to be calculated (instead of the first one). By default, the L=1,L1=1 sub-subsection will be calculated which in almost all cases will correspond to the first sub-subsection in the MF34 data (as a result, the default behaviour of NJOY2016 will not change).
+  - Increased allocation of arrays to accommodate ENDF/B-VIII.1 evaluations.
+  - Fixed a typo in the name for MT195 in ACER.
+
 ## [NJOY2016.71](https://github.com/njoy/NJOY2016/pull/301)
 This update adds the new MF7 MT451 (thermal scattering general information) ENDF format to MODER so that this module will be able to interpret the new MF7 section. No other capability in NJOY2016 currently uses the information in this section.
 

src/acecm.f90

@@ -67,7 +67,7 @@ subroutine mtname(mt,name,izai)
      '(n,2nhe3) ','(n,3nhe3) ','(n,4nhe3) ','(n,3n2p)  ','(n,3n2a)  ',&
      '(n,3npa)  ','(n,dt)    ','(n,npd)   ','(n,npt)   ','(n,ndt)   ',&
      '(n,nphe3) ','(n,ndhe3) ','(n,nthe3) ','(n,nta)   ','(n,2n2p)  ',&
-     '(n,phe3)  ','(n,dhe3)  ','(n,he3a)  ','(n,4n2p)  ','(n,2n2a)  ',&
+     '(n,phe3)  ','(n,dhe3)  ','(n,he3a)  ','(n,4n2p)  ','(n,4n2a)  ',&
      '(n,4npa)  ','(n,3p)    ','(n,n3p)   ','(n,3n2pa) ','(n,5n2p)  ',& !200
      '(n,p*0)   ',&
      '(n,p*1)   ','(n,p*2)   ','(n,p*3)   ','(n,p*4)   ','(n,p*5)   ',&
@@ -146,7 +146,15 @@ subroutine mtname(mt,name,izai)
      '(n,xhe3)  ','(n,xa)    '/)
    character(10)::hndf10(1)='damage    '
 
-   if (iverf.ge.6) then
+   !-- default name value is blank unless reset below.
+   name=''
+
+   !--when iverf=-1 the endf version used to create this file is
+   !  unknown.  we assume it comes from version 6 but if the file
+   !  originated from version 5 or earlier the name assigned here
+   !  might be incorrect.
+
+   if (iverf.ge.6.or.iverf.eq.-1) then
       if (mt.ge.201.and.mt.le.207) then
          name=hndf9(mt-200)
       else if (mt.eq.444) then

src/acedo.f90

@@ -334,7 +334,7 @@ subroutine dosfix(itype,nin,nout,ndir,iprint,nplot,mcnpx,&
       read(nin,'(8i9)')&
         len2,za,nxs3,ntr,nxsd(1:12),&
         lone,jxs2,mtr,jxs4,jxs5,lsig,sigd,jxsd(1:14),end,jxsd2(1:10)
-      n=(lone+3)/4
+      n=(len2+3)/4
       l=0
       do i=1,n
          read (nin,'(4e20.0)') (xss(l+j),j=1,4)

src/acepa.f90

@@ -804,6 +804,7 @@ subroutine phoprt(hk)
    !-------------------------------------------------------------------
    ! Print ACE photon interaction data from memory.
    !-------------------------------------------------------------------
+   use endf   ! provides iverf
    use mainio ! provides nsyso
    ! externals
    character(70)::hk
@@ -840,6 +841,11 @@ subroutine phoprt(hk)
    iabs=icoh+nes
    ipair=iabs+nes
    ihtng=lhnm-1
+
+   !--if this is a stand-alone acer iopt=7 job then iverf retains
+   !  its default -1 value and the values written as natural logs
+   !  need to be converted back to real numbers prior to printing.
+
    do i=1,nes
       if (mod(i,57).eq.1) write(nsyso,'(''1''/&
         &''     i'',8x,''energy'',4x,''incoherent'',&
@@ -849,18 +855,23 @@ subroutine phoprt(hk)
         &4x,''----------'',4x,''----------'')')
       col(1)=blank
       x=xss(ieg+i)
+      if (iverf.eq.-1.and.x.ne.zero) x=exp(x)
       if (x.ne.zero) write(col(1),'(1p,e14.4)') x
       col(2)=blank
       x=xss(iinc+i)
+      if (iverf.eq.-1.and.x.ne.zero) x=exp(x)
       if (x.ne.zero) write(col(2),'(1p,e14.4)') x
       col(3)=blank
       x=xss(icoh+i)
+      if (iverf.eq.-1.and.x.ne.zero) x=exp(x)
       if (x.ne.zero) write(col(3),'(1p,e14.4)') x
       col(4)=blank
       x=xss(iabs+i)
+      if (iverf.eq.-1.and.x.ne.zero) x=exp(x)
       if (x.ne.zero) write(col(4),'(1p,e14.4)') x
       col(5)=blank
       x=xss(ipair+i)
+      if (iverf.eq.-1.and.x.ne.zero) x=exp(x)
       if (x.ne.zero) write(col(5),'(1p,e14.4)') x
       col(6)=blank
       x=xss(ihtng+i)
@@ -906,6 +917,7 @@ subroutine phoout(itype,nout,ndir,mcnpx,hk,izn,awn)
    !-------------------------------------------------------------------
    ! Write photo-atomic ACE data to output and directory files.
    !-------------------------------------------------------------------
+   use endf  ! provides iverf
    use util  ! provides openz,closz,error
    use acecm ! provides write routines
    ! externals
@@ -941,11 +953,12 @@ subroutine phoout(itype,nout,ndir,mcnpx,hk,izn,awn)
         len2,z,nes,nflo,nxsd(1:12),&
         eszg,jinc,jcoh,jflo,lhnm,jxsd(1:27)
 
-      !--eszg block
+      !--eszg block.
+      !  convert to natural log, if not already done, prior to writing
       l=eszg
       n=5*nes
       do i=1,n
-         if (xss(l).ne.0.) xss(l)=log(xss(l))
+         if (xss(l).ne.0..and.iverf.ne.-1) xss(l)=log(xss(l))
          call typen(l,nout,2)
          l=l+1
       enddo

src/endf.f90

@@ -17,7 +17,7 @@ module endf
 
    !--Public variables
    integer,public::npage=306
-   integer,public::iverf
+   integer,public::iverf=-1
    real(kr),public::c1h,c2h
    integer,public::l1h,l2h,n1h,n2h,math,mfh,mth,nsh,nsp,nsc
    real(kr),public::thr6

src/errorr.f90

@@ -229,14 +229,14 @@ subroutine errorr
    !            (mf=32) (default=1)
    !            0 = area sensitivity method
    !            1 = 1% sensitivity method
-   !    legord  legendre order for calculating covariances (default=1)
+   !    l       Legendre order for reaction MT (default=1)
    !            (if mfcov is not 34, legord is ignored)
-   !    ifissp  subsection of the fission spectrum covariance
-   !            matrix to process (default=-1 which means process
+   ! l1/ifissp  if mfcov is 34: Legendre order for reaction MT1 (default=1)
+   !            if mfcov is 35: matrix to process (default=-1 which means process
    !            the subsection that includes efmean).  The value
    !            for ifissp that appears in njoy's standard output
    !            will equal the subsection containing efmean.
-   !            (if mfcov is not 35, ifissp is ignored)
+   !            (if mfcov is not 34 or 35, l1/ifissp is ignored)
    !    efmean  incident neutron energy (eV).  Process the covar-
    !            iance matrix subsection whose energy interval in-
    !            cludes efmean.  if ifissp=-1 and efmean is not
@@ -515,16 +515,21 @@ subroutine errorr
       mfcov=33
       irespr=1
       legord=1
-      ifissp=-1
+      ifissp=-1000
       isru=0
       dap=0
       read(nsysi,*) iread,mfcov,irespr,legord,ifissp,efmean,dap
+      if (mfcov.eq.35.and.ifissp.eq.-1000) ifissp=-1
+      if (mfcov.eq.34.and.ifissp.eq.-1000) ifissp=1
 
       !--only allow legord=1 at this time
-      if (legord.ne.1) then
-         write(strng,'(''reset legord from '',i2,'' to 1'')')legord
-         call mess('errorr',strng,'')
-         legord=1
+      if (mfcov.eq.34.and.legord.lt.0) then
+         write(strng,'(''Legendre order l for mf34 cannot be less than zero'')')
+         call error('errorr',strng,' ')
+      endif
+      if (mfcov.eq.34.and.ifissp.lt.0) then
+         write(strng,'(''Legendre order l1 for mf34 cannot be less than zero'')')
+         call error('errorr',strng,' ')
       endif
 
       !--input check for legal and/or consistent options
@@ -681,8 +686,12 @@ subroutine errorr
         '('' irespr processing option for mf=33 ... '',i10)') irespr
       if (isru.ne.0) write(nsyso,&
         '('' user dap for scat radius unc override '',f11.4)') dap
-      if (mfcov.eq.34) write(nsyso,&
-        '('' legendre order for mf=34 ............. '',i10)') legord
+      if (mfcov.eq.34) then
+        write(nsyso,&
+        '('' legendre order l for mf=34 ........... '',i10)') legord
+        write(nsyso,&
+        '('' legendre order l1 for mf=34 .......... '',i10)') ifissp
+      endif
       if (mfcov.eq.35) then
          write(nsyso,&
            '('' igflag ............................... '',i10/&
@@ -1771,7 +1780,7 @@ subroutine covcal
    integer::jh,loci,lt,lb,locip4,locip6,nk1,k,k2,locnec,nl1
    integer::ifloc,nlt,nk,nlt1,locl,lend,loclp4,loclp6
    integer::l2,m,m2,jgend,ih,ibase,ij,kmtb,isrrr,namx
-   integer::mat2,mt2,nlg1,nlg2,ld,ld1,izap,izero
+   integer::mat2,mt2,nlg1,nlg2,ld,ld1,izap,izero,ifoundmf34
    real(kr)::eg,xcv,time,de,flux,dne
    character(70)::strng,strng2
    integer,parameter::locm=100
@@ -1784,6 +1793,7 @@ subroutine covcal
    real(kr),parameter::zero=0
 
    !--initialize
+   ifoundmf34=0
    izero=0
    nscr2=13
    if (ngout.lt.0) nscr2=-nscr2
@@ -2051,7 +2061,8 @@ subroutine covcal
   320 continue
    if (iok.eq.0) go to 600
    if (mfcov.eq.34) then
-      if (ld.gt.legord.or.ld1.gt.legord) go to 650
+      if (.not.(ld.eq.legord.and.ld1.eq.ifissp)) go to 650
+      ifoundmf34=1
    endif
 
    !--retrieve sigma for mt1, either from ngout or sig.
@@ -2363,6 +2374,12 @@ subroutine covcal
   660 continue
    call asend(0,nscr1)
 
+   !--check if we found the MF34 l,l1 sub-subsection
+   if (mfcov.eq.34.and.ifoundmf34.eq.0) then
+      write(strng,'(''no sub-subsection for l='',i2,'' and l1='',i2)') legord,ifissp
+      call error('covcal',strng,'in mf=34.')
+   endif
+
    !--close loop over sections of mfcov
    go to 140
   700 continue
@@ -5653,7 +5670,11 @@ subroutine grpav4(mprint)
          z(6)=ig
          z(7)=ans(1,1)
          do i=1,legord
-            z(i+7)=ans(i,2)
+           if (ans(i,2).lt.1e-30) then
+             z(i+7)=0.
+           else
+             z(i+7)=ans(i,2)
+           endif
          enddo
          nwds=legord+7
          call listio(0,nscr4,0,z,nb,nwds)
@@ -6049,7 +6070,6 @@ subroutine rdlgnd(nscr4,matd,mtd,npl,b,alp)
    return
    end subroutine rdlgnd
 
-
    subroutine fssigc(ncg,ncm,nun,csig,cflx,b,egt,flux,sig)
    !--------------------------------------------------------------------
    ! Calculate the coarse group fission spectrum chi.

src/groupr.f90

@@ -6440,7 +6440,7 @@ subroutine getyld(e,enext,idis,yld,mat,mf,mt,lfs,itape)
 
    !--initialize
    if (e.gt.zero) go to 200
-   ntmp=10000
+   ntmp=20000
    allocate(tmp(ntmp))
    mft=mf
    if (mft.ge.40000000) mft=10
@@ -7436,7 +7436,7 @@ subroutine getmf6(ans,ed,enext,idisc,yld,eg,ng,nl,iglo,ng2,nq,&
    real(kr)::ed,enext,yld
    real(kr)::ans(nlg,*),eg(*)
    ! internals
-   integer::mfn,nb,nw,lct,lct3,ik,nne,ne,int,nss
+   integer::mfn,nb,nw,lct3,ik,nne,ne,int,nss
    integer::ie,ilo,jlo,jhi,ii,nn,nnn,langn,lepn,idis,jzap
    integer::nk,jzad,lang,lep,i,npsx,irr,npp,nmu,l1
    integer::j,iss,ip,ir,jgmax,jj,jg,ndlo,nplo,nclo,nphi,nchi
@@ -7466,7 +7466,7 @@ subroutine getmf6(ans,ed,enext,idisc,yld,eg,ng,nl,iglo,ng2,nq,&
    save nne,ne,int
    save jlo,elo,jhi,ehi,terml
    save pspmax,langn,lepn,disc102
-   save idis,iyss,izss,jjss,jloss,nss,jzap,lct3,lct
+   save idis,iyss,izss,jjss,jloss,nss,jzap,lct3
 
    !--initialize
    if (ed.gt.zero) go to 200
@@ -11385,15 +11385,15 @@ subroutine conver(nin,nout,nscr)
    call contio(nin,nout,nscr,scr,nb,nw)
    if (mfh.eq.0) go to 110
    if (mth.ne.452) go to 595
-   nnu=8000
+   nnu=20000
    allocate(nu(nnu))
    l=1
    lnu=l2h
    if (lnu.eq.1) call listio(nin,nout,nscr,nu(l),nb,nw)
    if (lnu.eq.2) call tab1io(nin,nout,nscr,nu(l),nb,nw)
    do while (nb.ne.0)
       if (l+nw.gt.nnu) call error('conver',&
-        'storage for fission nu exceeded',' ')
+        'storage for fission nu exceeded (increase the size of the nu array)',' ')
       l=l+nw
       call moreio(nin,nout,nscr,nu(l),nb,nw)
    enddo

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.71'
-   character(8),public::vday='14Jul23'
+   character(8),public::vers='2016.72'
+   character(8),public::vday='29SepXX'
 end module version

tests/15/input

@@ -47,7 +47,7 @@ errorr
  -21 0 91 27 0 0 /
  9237 3 6 1 1 /
  1 300./
- 0 34 1 1 -1 /
+ 0 34 1 1 1 /
 --
 -- make mf31 plot file.
 covr