*/
*/ OECD/NEA compilation of NJOY updates
*/ compatible with the official patch distribution up296
*/
*/ The patches include those carried-over from the set compatible
*/ with the official up259 but not adopted in up296.
*/
*ident upnea001
*/ njoy - C. Broeders, 6-Oct-2006
*/        Open scratch files in local disc area
*/        Ref.: Comments by D.E. Cullen in PrePro-2004 (Scratcha.f) on p.19 of:
*/        "http://www-nds.iaea.org/ndspub/endf/prepro/DOCUMENT/PDF/Overview.pdf"
*/        Status "scratch" does not work with Lahey compiler on Linux,
*/        therefore 'age" is set "unknown"
*d njoy.470
         age='unknown'
*d njoy.471
         write(fn,'(a,i2.2)') 'temp',nun
         open(nun,file=fn,form=for,status=age)  
*/
*ident upnea004
*/ groupr - C. Broeders, 6-Oct-2006
*/          extend IWT=5 spectrum definition up to 200MeV
*/ -------------------------------------------
*d groupr.2145
      dimension w1(92),w2(92),w3(10)
*d groupr.2150
      data w1/0d0,0d0,0d0,0d0,1.d0,93.d0,93.d0,5.d0,
*i groupr.2175
      data w3/3.d7,1.0318d-10,5.d7,6.1908d-11,1.d8,3.0954d-11,
     &1.5d8,2.0636d-11,2.d8,1.5477d-11/
*d groupr.2358
      iw=194
*i groupr.2365
      do i=1,10
         a(i+183+iwght)=w3(i)
      enddo
*/
*/ *ident upnea019 is included in up265
*/
*/ *ident upnea020 is included in up235
*/
*/ *ident upnea021 is included in up261
*/
*/ *ident upnea022 is included in up260
*/
*/ *ident upnea023 is included in up272
*/
*/ *ident upnea024 is included in up273
*/
*ident upnea025
*/ reconr A. Trkov 09 Dec 2007
*/    Allow processing of MT261, usen in the dosimetry library IRDF
*/    to store reference neutron fields (spectra).
*/    The patch has no impact on cross section processing.
*d reconr.1689      
      if ((mth.ge.251.and.mth.le.300).and.mth.ne.261) go to 150
*/
*ident upnea026
*/ groupr A. Trkov 09 Dec 2007
*/    Allow processing of MT261, usen in the dosimetry library IRDF
*/    to store reference neutron fields (spectra).
*/    The patch has no impact on cross section processing.
*i groupr.3967
         if (mtd.eq.261) mt=261
*/
*ident upnea027
*/ acer A. Trkov 14 Jan 2008
*/      Implement enhancements that allow processing of dosimetry data
*/      including metastable nuclide excitation functions in MF10.
*/      Internally the ace MT numbers were extended according to the
*/      following convention:
*/        MT* = MT + 1000*(10+LFS)
*/      where LFS is the final state of the nuclide.
*i acer.11525
c        strip the leading digits from MT in dosimetry reactions
         if (i.gt.999) i=i-1000*(i/1000)
*d acer.14319
c     reserve all available space for scratch
      nwscr=-1
*d acer.14357
      do while (mfh.ne.3 .and. mfh.ne.10)
         call tofend(nin,0,0,a(iscr))
         call contio(nin,0,0,a(iscr),nb,nw)
         if (mfh.gt.10 .or. math.le.0)
     &     call error('acedos','no x-sect. data for desired mat.',' ')
      end do
*d acer.14367
      if(mfh.gt.3) go to 110
*d acer.14371
*/ guard against array overflow
*i acer.14380
               if(jscr.gt.nwscr)
     &           call error('acedos','aray storage exceeded'
     &                     ,' Execution terminated')
*i acer.14383
                 if(jscr.gt.nwscr)
     &             call error('acedos','aray storage exceeded'
     &                       ,' Execution terminated')
*i acer.14409
            call contio(nin,0,0,a(iscr),nb,nw)
*i acer.14411
c     ***locate first reaction in file 10
      call contio(nin,0,0,a(iscr),nb,nw)
      if (mfh.gt.10 .or. math.le.0) goto 120
      do while (mfh.ne.10)
         call tofend(nin,0,0,a(iscr))
         call contio(nin,0,0,a(iscr),nb,nw)
         if (mfh.gt.10 .or. math.le.0) goto 120
      end do
  110 continue
c     ***loop over reactions on nin
      do while (mfh.ne.0)
         ns=max(1,n1h)
         xss(lsig-1+j)=l
         if (mfh.ne.0) then
            if (mth.ne.1) then
               do is=1,ns
                  jscr=iscr
                  call tab1io(nin,0,0,a(jscr),nb,nw)
                  lfs=l2h
                  xss(mtr-1+j)=mth+1000*(10+lfs)
                  nr=nint(a(iscr+4))
                  ne=nint(a(iscr+5))
                  intr=nint(a(iscr+7))
                  jscr=jscr+nw
                  if(jscr.gt.nwscr)
     &              call error('acedos','aray storage exceeded'
     &                        ,' Execution terminated')
                  do while (nb.ne.0)
                     call moreio(nin,0,0,a(jscr),nb,nw)
                     jscr=jscr+nw
                    if(jscr.gt.nwscr)
     &              call error('acedos','aray storage exceeded'
     &                        ,' Execution terminated')
                  enddo
                  if (nr.ne.1.or.intr.ne.2) then
                     xss(l)=nr
                     l=l+1
                     do i=1,nr
                        xss(l+2*i-2)=a(iscr+4+2*i)
                        xss(l+2*i-1)=a(iscr+5+2*i)
                     enddo
                     l=l+2*nr
                  else
                     xss(l)=0
                     l=l+1
                  endif
                  xss(l)=ne
                  k=iscr+6+2*nr
                  l=l+1
                  do i=1,ne
                     xss(l)=a(k)/emev
                     xss(l+ne)=a(k+1)
                     l=l+1
                     k=k+2
                  enddo
                  l=l+ne
                  j=j+1
               enddo
            endif
            call tosend(nin,0,0,a(iscr))
            call contio(nin,0,0,a(iscr),nb,nw)
         endif
      enddo
  120 continue
*/ allow more digits for dosimetry mt's
*d acer.14496
     &        '' reaction mt = '',i6,3x,a10/'' interpolation: '',12i6)')
*d acer.14501
     &        '' reaction mt = '',i6,3x,a10/'' linear interpolation'')')
*/ allow processing of dosimetry files
*d acer.407  
      else if (iopt.lt.7.or.iopt.gt.9) then
*d acer.435  
      else if (iopt.ge.7.and.iopt.le.9) then
*i acer.17639
c
c     ***read type 3 ace format file
      else if (itype.eq.2) then
         if (mcnpx.eq.0) then
           read(nin) hz(1:10),aw0,tz,hd,hko,hm,(izo(i),awo(i),i=1,16),
     &       (nxs(i),i=1,16),(jxs(i),i=1,32)
         else
           read(nin) hz(1:13),aw0,tz,hd,hko,hm,(izo(i),awo(i),i=1,16),
     &       (nxs(i),i=1,16),(jxs(i),i=1,32)
         endif
         len2=nxs(1)
         n=(len2+ner-1)/ner
         l=0
         do i=1,n
            max=len2-l
            if (max.gt.ner) max=ner
            read (nin) (xss(l+j),j=1,max)
            l=l+max
         enddo
         call closz(-nin)
c        flag incident particle "undefined"
         izai=-1
         awi =-1
*d acer.17692
      else if (ht.eq.'d'.or.ht.eq.'y') then
*/
*ident upnea028
*/ acer A. Trkov 22 Jan 2008
*/      Increase the size of the work array from 180000 to 510000
*/      to process U-238 and Gd-nat from IRDF-2002.
*d up123.5
      common/astore/a(510000)
*d up123.7
      data namax/510000/, nidmax/27/
*d up123.9
      common/astore/a(510000)
*d up123.11
      common/astore/a(510000)
*d up123.13
      common/astore/a(510000)
*d up123.15
      common/astore/a(510000)
*d up123.17
      common/astore/a(510000)
*d up123.19
      common/astore/a(510000)
*d up123.21
      common/astore/a(510000)
*d up123.23
      common/astore/a(510000)
*d up123.25
      data namax/510000/
*d up123.27
      common/astore/a(510000)
*d up123.29
      common/astore/a(510000)
*d up123.31
      common/astore/a(510000)
*d up123.33
      common/astore/a(510000)
*d up123.35
      common/astore/a(510000)
*d up123.37   
      common/astore/a(510000)
*d up123.39
      common/astore/a(510000)
*d up123.41
      common/astore/a(510000)
*/
*ident upnea029
*/ groupr  A. Trkov, 6-Feb-2008
*/         Define minimum cross section for printout as 1.e-9
*/         for consistency with other parts of the code.
*i groupr.4286
      data smin/1.d-9/
*i groupr.4289
      data smin/1.e-9/
*d groupr.4432 
      if (a(ised-1+ik+nk*(ig-1)).gt.smin) go to 280
*/
*/ *ident upnea030 included in up272
*/
*/ *ident upnea031 included in up289
*/
*ident upnea032
*/ acer     A. Trkov, 30 May 2008
*/          Plots of discrete inelastic make sense below 5 MeV
*d acer.19586
c                       limit discrete inelastic plots to 5 MeV
                        if (x.gt.xmax) xmax=min(x,5.)
*/
*/ *ident upnea033 included in up278
*/
*/ *ident upnea034 included in up277
*/
*/ *ident upnea035 included in up279
*/
*/ *ident upnea036 included in up279
*/
*/ *ident upnea037 included in up277
*/
*ident upnea038
*/ covr    A. Trkov, 20 October 2008
*          Fix diagnostic message
*d covr.1793,1805
      zero=0
      if (abs(crit).gt.zero) write(strng,                                                      covr.1794 
     &      '(''i'',i3,'' j'',i3,'' xa(i,j)='',1p,e12.4,'' xa(j,i)='',
     &     1p,e12.4)') i,im1,xa(im1ind,i),xa(iind,im1)
      if (abs(crit).gt.test .and. abs(xa(im1ind,i)).gt.1e-20) then                                          covr.1793 
         call error('press','matrix not symmetric',strng)
      endif
      test=1
      test=test/1000000
      if (abs(crit).lt.test) go to 250
      call mess('press','matrix not symmetric',strng)
*/
*/ *ident upnea039 included in up272
*/
*ident upnea040
*/
*/ GROUPR - R. Perry, 15 October 2008
*/          (Original update name: upwe6)
*/          Problem with thermal scattering matricies in GROUPR            
*/          There were small discrepancies between the sums thermal scatter                   
*/          matrices and the group cross-sections. This update renormalises                  
*/          the interpolated spectra (produced during integration)                           
*/          to sum to unity.                                                                 
*/          (guarding against array overflow added by A. Trkov)
*i groupr.6915
      integer nsi
      parameter (mxnsi=15000)
*i groupr.6922
      dimension sint(2,mxnsi)
*i groupr.7030
      nsi = 0
*i groupr.7146
c
c store interpolated function. 
c
      nsi = nsi+1
      if(nsi.gt.mxnsi)
     &   call error('getaed','mxnsi array limit exceeded',' ')
      sint(1,nsi) = ei
      sint(2,nsi) = fi(1)
*i groupr.7155
c
c sum interpolated spectrum and adjust aed. 
c
      sisum = 0.0
      do il = 1, nsi - 1
        sisum = sisum + (sint(2,il)+sint(2,il+1))*
     &                                  (sint(1,il+1)-sint(1,il))/2
      enddo
c
c don't re-normalise if scatter source above top energy group. 
c
      if(i .lt. ngn .or. (k1 .eq. nlo .and. k2 .eq. nhi) ) then 
        do il = 1, i
          aed(1,il) = aed(1,il)/sisum
        enddo
      endif
*/
*/ *ident upnea041 included in up282
*/
*ident upnea042
*/ groupr  A. Trkov, 6-Feb-2009
*/         Truncation of the minimum cross section causes problems
*/         in the case of the fission spectrum when processing
*/         covariance matrices. Update upnea029 is refined to
*/         decrease the tolerance by a factor 1000, and abolish
*/         it if processing mt=18 spectrum
*d upnea029.10
      if (a(ised-1+ik+nk*(ig-1)).gt.smin/1000.or.mtd.eq.18) go to 280
*/
*/ *ident upnea043 included in up282
*/
*/ *ident upnea044 included in up273
*/
*/ *ident upnea045 included in up282
*/
*/ *ident upnea046 withdrawn
*/
*/ *ident upnea047 included in up287
*/
*ident upnea048
*/ errorr  Go Chiba, 3-Apr-2009
*/         (implemented by D. Rochman and A. Hogenbirk)
*/         (ident name changed from roch01 to upnea048)
*/         - include treatment of scattering radius uncertainty
*i errorj.6360
      dimension pneorg(10000)
*i errorj.6491
c
c --- For scattering radius uncertainty treatment
c
c +++ reference (no perturbation)
         call resprx_cal_pendf(ii,99,0.d0,a,sigr,-1.d0,b,maxb)
c +++ perturbated system
c     (perturbation for scattering radius of MF=2)
c     (perturbation for penetration factor)
         inow=1
         ap=b(inow+7)
         dap=ap*0.01
         ap=ap*1.01
         b(inow+7)=ap
         nls=nint(b(inow+10))
         inow=inow+12
         itmp=1
         do lll = 1 , nls
           apl=b(inow+1)
           if(apl.eq.0) apl=ap
           apl=apl*1.01
           b(inow+1)=apl
           ll=nint(b(inow+2))
           nrs=nint(b(inow+5))
           inow=inow+6
           do jj = 1 , nrs
             rho=cwaven*arat*sqrt(abs(b(inow+6*(jj-1))))*apl
             call facts(ll,rho,ser,per)
             pneorg(itmp)=b(inow+nrs*6+(jj-1)*3)
             pneorg(itmp+1)=b(inow+nrs*6+(jj-1)*3+1)
             itmp=itmp+2
             b(inow+nrs*6+(jj-1)*3)=ser
             b(inow+nrs*6+(jj-1)*3+1)=per
           enddo
           inow=inow+6*nrs+3*nrs
         enddo
c
         call resprx_cal_pendf(ii,99,0.d0,a,sigp,-1.d0,b,maxb)
c     (Sensitivity calculation)
         dap=1./dap
         do ii1=1,4
         do ii2=1,ii
           tmp=((sigp(ii2,ii1)-sigr(ii2,ii1)))*dap
           sigp(ii2,ii1)=tmp
         enddo
         enddo
         call Resprx_grping(ngn,egn,sigp,ii,gsig(1,1),a)
         do ii1=1,ngn
           tmp=a(icflx-1+ii1)*abn
           do j = 1 , 4
             gsig(j,ii1)=gsig(j,ii1)*tmp
           enddo
         enddo
c     (To get absolute standard deviation for scattering radius)
         erap=a(l+1)
         erap=erap*erap
c     (Error propagation)
         igind=0
         do ig = 1 , ngn
           do ig2 = ig , ngn
             igind=igind+1
             itmp1=icff+igind-1
             itmp2=icgg+igind-1
             itmp3=icee+igind-1
             itmp4=ictt+igind-1
             a(itmp1)=a(itmp1)+erap*gsig(3,ig)*gsig(3,ig2)
             a(itmp2)=a(itmp2)+erap*gsig(4,ig)*gsig(4,ig2)
             a(itmp3)=a(itmp3)+erap*gsig(2,ig)*gsig(2,ig2)
             a(itmp4)=a(itmp4)+erap*gsig(1,ig)*gsig(1,ig2)
           enddo
         enddo
c
         igind=0
         do ig = 1 , ngn
           do ig2 = 1 , ngn
             igind=igind+1
             itmp1=icef+igind-1
             itmp2=iceg+igind-1
             itmp3=icfg+igind-1
             a(itmp1)=a(itmp1)+erap*gsig(2,ig)*gsig(3,ig2)
             a(itmp2)=a(itmp2)+erap*gsig(2,ig)*gsig(4,ig2)
             a(itmp3)=a(itmp3)+erap*gsig(3,ig)*gsig(4,ig2)
           enddo
         enddo
c
c     (Go back to the reference situation)
         inow=1
         ap=b(inow+7)
         ap=ap*0.9901
         b(inow+7)=ap
         nls=nint(b(inow+10))
         inow=inow+12
         do lll = 1 , nls
           apl=b(inow+1)
           apl=apl*0.9901
           b(inow+1)=apl
           nrs=nint(b(inow+5))
           inow=inow+6+(nrs*9)
         enddo
c --- End of scattering radius uncertainty treatment
*/
*d errorj.8345
            if((abs(ajc-valspi).gt.0.01.or.l.ne.npnls).and.
     c         npnls.ne.99) then
*/
*/ *ident upnea049 is included in up290
*/