Fixes #383: robustness and portability fixes; WIMSR WLUP conventions for WIMSD-5 libraries

src/acepa.f90

@@ -23,7 +23,8 @@ module acepa
    integer::eszg,jinc,jcoh,jflo,lhnm,jxsd(27)
 
    ! parameter for scratch array
-   integer,parameter::nwscr=50000
+   ! modified (IAEA, D. Lopez Aldama, dlopezaldama@gmail.com)
+   integer,parameter::nwscr=250000
 
 contains
 

src/endf.f90

@@ -525,8 +525,12 @@ subroutine tpidio(nin,nout,nscr,a,nb,nw)
    integer::inin,inout,inscr,i
 
    !--input.
+   ! the tape identification is the first record: the input unit is
+   ! rewound before it is read (IAEA, D. Lopez Aldama,
+   ! dlopezaldama@gmail.com)
+   inin=iabs(nin)
+   if (inin.ne.0) rewind inin
    if (nin.lt.0) then
-      inin=iabs(nin)
       read(inin) math,mfh,mth,nb,nw,(a(i),i=1,17)
    else if (nin.gt.0) then
       read(nin,'(16a4,a2,i4,i2,i3,i5)') (hb(i),i=1,17),math,mfh,mth

src/leapr.f90

@@ -521,8 +521,10 @@ subroutine contin(temp,itemp,np,maxn)
    endif
    if (maxn.gt.maxnphon) then
       call timer(time)
+      ! no comma before the output list (non-standard) (IAEA,
+      ! D. Lopez Aldama, dlopezaldama@gmail.com)
       write(nsyse,'(/'' performing phonon expansion sum'',&
-            &37x,f8.1,''s'')'),time
+            &37x,f8.1,''s'')')time
    endif
    do n=2,maxn
       npn=np+npl-1
@@ -559,7 +561,7 @@ subroutine contin(temp,itemp,np,maxn)
    if (maxn.gt.maxnphon) then
       call timer(time)
       write(nsyse,'(/'' done with phonon expansion sum'',&
-            &38x,f8.1,''s'')'),time
+            &38x,f8.1,''s'')')time
    endif
 
    !--print out start of sct range for each beta
@@ -1139,6 +1141,8 @@ real(kr) function terps(sd,nsd,delta,be)
    terps=0
    if (be.gt.delta*nsd) return
    i=int(be/delta)
+   ! negative index: out of the table (DLA, IAEA)
+   if (i.lt.0) return
    if (i.lt.nsd-1) then
       bt=i*delta
       btp=bt+delta

src/main.f90

@@ -2,7 +2,7 @@ program njoy
 !-----------------------------------------------------------------------
 !
 !    NJOY Nuclear Data Processing System
-!    Version 2016.20
+!    Version 2016 (see vers.f90)
 !
 !-----------------------------------------------------------------------
 !

src/plotr.f90

@@ -383,9 +383,10 @@ subroutine plotr
    real(kr)::st1,st2,st3,st4,stmax,stmin,ststp,temp,dleth,dener
    real(kr)::sigig,flag,factx1,facty1,temper2,reset,enxt,zz,z2,itypx
    real(kr)::z(15)
-   integer,parameter::mmax=20000   !same in plotr and viewr
+   ! modified (IAEA, D. Lopez Aldama, dlopezaldama@gmail.com)
+   integer,parameter::mmax=50000   !same in plotr and viewr
    integer,parameter::nwamax=45000
-   integer,parameter::maxaa=200000
+   integer,parameter::maxaa=500000
    real(kr),dimension(nwamax)::a
    real(kr),dimension(maxaa)::aa
    real(kr),dimension(mmax)::x,y,b,dxm,dxp,dym,dyp

src/util.f90

@@ -193,6 +193,11 @@ subroutine openz(lun,new)
    if (nun.ge.10.and.nun.le.19) then
       ! scratch units
       age='scratch'
+      ! a scratch unit still open from a previous module is closed
+      ! before it is opened again (IAEA, D. Lopez Aldama,
+      ! dlopezaldama@gmail.com)
+      inquire(nun,opened=there)
+      if (there) close (nun)
       open(nun,form=for,status=age)
    else
       ! regular units

src/viewr.f90

@@ -339,8 +339,9 @@ subroutine viewr
    real(kr)::factx,facty
    real(kr)::xx,yy,zz
    character(80)::text
-   integer,parameter::mmax=20000   !same in plotr and viewr
-   integer,parameter::maxaa=500000
+   ! modified (IAEA, D. Lopez Aldama, dlopezaldama@gmail.com)
+   integer,parameter::mmax=50000   !same in plotr and viewr
+   integer,parameter::maxaa=5000000
    real(kr),dimension(15)::z
    real(kr),dimension(maxaa)::aa
    real(kr),dimension(mmax)::x,y,b,dxm,dxp,dym,dyp
@@ -689,7 +690,7 @@ subroutine viewr
                      idone=1
                   else
                      if (iskip.eq.0) then
-                       if (l+5000.ge.maxaa) then
+                       if (l+mmax.ge.maxaa) then
                            call mess('viewr','too much 3d data',&
                              'data truncated')
                            iskip=1
@@ -1270,16 +1271,18 @@ subroutine set3d(iplot,xyz,nxyz)
    real(kr)::xyz(nxyz)
    ! internals
    integer,parameter::length=2000
+   ! modified (DLA, IAEA)
+   integer,parameter::lenxyz=10000
    integer::i,j,nn,k,ncurv,major,minor,itop,ibot,l,n
    real(kr)::wt,xlo,xhi,ylo,yhi,zlo,zhi,xn,yn,zn,yy
    real(kr)::top,bot,xop,yop,zop
    integer::lll(length)
-   real(kr)::x(2000),y(2000),z(2000)
+   real(kr)::x(lenxyz),y(lenxyz),z(lenxyz)
    real(kr),parameter::big=1.e10_kr
    real(kr),parameter::d0=.001e0_kr
    real(kr),parameter::d3=.301e0_kr
    real(kr),parameter::d7=.699e0_kr
-   integer::kmax=1999
+   integer::kmax=lenxyz-1  !  x(:),y(:),z(:) arrays  dimension - 1
    real(kr),parameter::zero=0
    real(kr),parameter::one=1
    real(kr),parameter::ten=10

src/wimsr.f90

@@ -30,7 +30,10 @@ module wimsm
    integer,parameter::nymax=100
    real(kr)::yield(nymax)
    integer::ifisp(nymax)
-   integer::ifiss,nfiss
+   ! ifissr: fission flag of the resonance tables (WIMSD-5 library
+   ! conventions of the IAEA WIMS-D Library Update Project, WLUP)
+   ! (IAEA, D. Lopez Aldama, dlopezaldama@gmail.com)
+   integer::ifiss,ifissr,nfiss
    integer::isg
    integer::ixs
    real(kr),dimension(:),allocatable::snu
@@ -434,7 +437,7 @@ subroutine resint
    integer::i,jg,is,jtemp,jfiss,nl,nz,ntw,l,ig,kg,lim,iadd
    integer::loca,jz,loc,locn,jtem,iterm,ioff,it,index
    integer::indexl,iz
-   real(kr)::xid,siglam,sigb,siga,sig
+   real(kr)::xid,siglam,sigb,siga,sigs,sig
    real(kr),dimension(:),allocatable::sabs
    real(kr),dimension(:),allocatable::snux
    real(kr),dimension(:),allocatable::snsf
@@ -447,6 +450,7 @@ subroutine resint
    real(kr),parameter::zero=0
 
    !--allocate storage.
+   ifissr=0
    if (ires.eq.0) go to 510
    write(nsyso,'(/'' ***resonance integrals***'')')
    ntsr=ires*nsigz*nrg
@@ -519,6 +523,7 @@ subroutine resint
    if (mfh.ne.3) go to 300
    if (mth.eq.1) go to 162
    if (mth.eq.2) go to 162
+   if ((mth.gt.50).and.(mth.le.91)) go to 162
    if (mth.lt.18.or.mth.gt.150) go to 300
    if (mth.gt.21.and.mth.lt.102.and.mth.ne.38) go to 300
   162 continue
@@ -537,7 +542,8 @@ subroutine resint
    jg=kg-nfg
    lim=nsigz
    if (nsigz.gt.nz) lim=nz
-   if (mth.ne.1.and.mth.ne.2.and.mth.ne.18.and.mth.ne.102) go to 300
+   if (mth.ne.1.and.mth.ne.2.and.mth.ne.18.and.mth.ne.102 &
+        .and.((mth.le.50).or.(mth.gt.91))) go to 300
 
    !--absorption
    if (mth.eq.102) then
@@ -578,8 +584,10 @@ subroutine resint
          sabs(iadd-jz+1)=sabs(iadd-jz+1)+scr(nl*jz+loca)
       enddo
 
-   !--elastic
-   else if (mth.eq.2) then
+   !--elastic or inelastic scattering: the discrete inelastic levels
+   !--(mt51-91) are added to the resonance scattering, as in the
+   !--WLUP libraries (DLA, IAEA)
+   else if ((mth.eq.2).or.((mth.gt.50).and.(mth.le.91))) then
       iadd=nsigz+nsigz*(jtemp-1+ires*(jg-1))
       loca=l+lz+nl*(nz-1)
       do jz=1,lim
@@ -661,13 +669,17 @@ subroutine resint
             sigb=sigz(iz)+siglam
             siga=sabs(loc+iz-1)
             sig=snux(loc+iz-1)
+            sigs=elas(loc+iz-1)
             sabs(loc+iz-1)=sigb*siga/(sigb+siga)
             snux(loc+iz-1)=sigb*sig/(sigb+siga)
+            ! wims-e/wimsd-5: scattering resonance integral too (DLA, IAEA)
+            if (iverw.ne.4) elas(loc+iz-1)=sigb*sigs/(sigb+siga)
          enddo
       enddo
    enddo
 
    !--write out results
+   ifissr=jfiss
    call rsiout(xid,jfiss,fa,sabs,snux,flxr,sigz,elas)
 
    !--resint is finished.
@@ -929,6 +941,9 @@ subroutine xsecs
    call tpidio(ngendf,0,0,scr,nb,nw)
    isg=0
    if (sgref.lt.dilinf) isg=1
+   ! wims-e/wimsd-5: data taken at the reference dilution sgref
+   ! also when sgref is infinite (DLA, IAEA)
+   if (iverw.eq.5) isg=1
    jtemp=0
    i318=0
    jfisd=0
@@ -944,11 +959,16 @@ subroutine xsecs
    xxi=1+log(alf)*alf/(1-alf)
 
    !--loop over temperatures
+   ! snu, abs1, ab0, xtr and sdp were not initialized (DLA, IAEA)
+   do i=1,ngnd
+      snu(i)=0
+   enddo
   140 continue
    do i=1,ngnd
       csp1(i)=0
       spot(i)=sigp
       xi(i)=xxi
+      abs1(i)=0
       abs2(i)=0
       sf0(i)=0
       sfi(i)=0
@@ -957,6 +977,9 @@ subroutine xsecs
       chi(i)=0
       scat(i)=0
       sn2n(i)=0
+      ab0(i)=0
+      xtr(i)=0
+      sdp(i)=0
       l1e(i)=ngnd
       l2e(i)=1
       l1(i)=ngnd
@@ -2040,6 +2063,8 @@ subroutine wimout(jcc)
 
    !--write material identifier data
    awt=awr*amassn
+   ! no fission in the resonance tables: ifis=2 (DLA, IAEA)
+   if((ifis.eq.3).and.(ifissr.ne.3)) ifis=2
    write(nout,'(i6,1p,e15.8,5i6)')&
      ident,awt,iznum,ifis,ntemp,nrestb,isof
 

tests/06/referenceTape32

@@ -1247,8 +1247,10 @@ newpath
    375.27   629.20 lineto
    368.52   629.32 lineto
    365.62   629.37 lineto
+   335.84   629.39 lineto
    319.61   629.40 lineto
    321.20   629.42 lineto
+   334.52   629.43 lineto
    337.54   629.44 lineto
    358.85   629.48 lineto
    371.15   629.54 lineto