*/
*/ OECD/NEA compilation of NJOY updates
*/ compatible with the official patch distribution up304
*/
*/ 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
*/
*/ *ident upnea050 is included in up289
*/
*ident upnea051
*/ errorr Tsuneo Nakagawa, Go Chiba, July 2009
*/        (original idents provided by the authors: nkgw, cbg)
*/        More robust treatment to identify resonance parameters
*/        with given uncertainties (the list can be reduced from
*/        the end, but no skipping is allowed).
*i errorj.6734
      data mpid /1,3,4,5,6,0/
*d up279.38,50
*i up279.56
      mm=0   
*d up279.61,67
*d up287.9
         do m=1,5
           mm=mm+1
           ii=l3+mpid(m)-1
           cov(mm,mm)=a(ii)
         enddo
*d up279.74,76
      l3=lbg+6 
      l2=lbg+6 
      do n2=1,nrb
*i up279.100
c
      mpar=nnn/nrb
      if(mpar.ne.5)then
        mm=0
        do nr=1,nrb
          mbase=5*(nr-1)
          do m=1,mpar
            mm=mm+1
            mmm=mbase+m
            cov(mm,mm)=cov(mmm,mmm)
          enddo
        enddo
      endif
c
      if(nm.ne.0)then
      nw=ndigit
*i up279.122
      endif
*/ 
*/ Modification for high-speed calculation option of errorr
*/
*d up290.49
            e2=1.0+(eskip1-1.0)*5.
*/
*/ Bug-fix for scattering radius uncertainty treatment
*/
*d upnea048.27,28
           if(apl.eq.0) then
             apl=ap 
           else 
             apl=apl*1.01                                              
           endif
*i upnea048.99
         itmp=1
*d upnea048.105
           inow=inow+6
           do jj = 1 , nrs
             b(inow+nrs*6+(jj-1)*3)=pneorg(itmp)               
             b(inow+nrs*6+(jj-1)*3+1)=pneorg(itmp+1)           
             itmp=itmp+2                                       
           enddo
           inow=inow+6*nrs+3*nrs
*/
*ident upnea052
*/ covr  R. Perel, July 2009
*/       Patches up273 and up283 have been erroneously applied
*/       inside an update "*if sw" statement. Restore the
*/       integrity of the *if statements for long-word machines.
*i covr.947
*else
      xsize=5.00e0
      ysize=3.38e0
*endif
*d covr.948,951
*i covr.1418
*endif
*d covr.1419
*/
*/ the following correction has been performed as part of up273
*/ but it was applied there to sw-mode only
*/ here: apply it if not sw-mode
*d up192.7
      data tlev/0.001e0,0.1e0,0.2e0,0.3e0,0.6e0,1.0e0/
*/
*ident upnea053
*/ errorr A. Trkov, July 2009
*/        Fix trivial error in printout labels (the "if" statement
*/        was wrong, the printout was fixed incorrectly at some
*/        previous step).
*d errorj.1867,1871
            if (mats(ixp).ne.0) then
               if(irelco.eq.0)write(nsyso,40) mt,mats(ixp),mts(ixp),time
               if(irelco.eq.1)write(nsyso,45) mt,mats(ixp),mts(ixp),time
*/
*ident upnea054
*/ errorr D.L. Aldama, July 2009
*/        Fix weighting function option
*d errorj.660
      if (iwt.eq.1.or.iwt.eq.4.or.iwt.eq.5.or.
     &     (iwt.ge.7.and.iwt.lt.11))
*d errorj.1492
      if (iwt.eq.1.or.iwt.eq.4.or.iwt.eq.5.or.
     &     (iwt.ge.7.and.iwt.lt.11))
*/
*ident upnea055
*/ purr  A. Trkov, August 2009
*/       intunr is already defined if all parameters are 
*/       energy-independent.
*/       move definition of intunr inside the if statement
*/       for other cases.
*/       The fix is for consistency and has no influence on the
*/       results
*d up293.44
*b purr.768
            intunr=5
*/
*ident upnea056
*/ heatr A. Trkov, October 2009
*/       external declaration missing
*d heatr.5005
      external tofend,tosend,tomend,closz ,error                             
*/
*ident upnea057
*/ errorr A. Trkov, October 2009
*/        Sign-flag on AJ was incorrectly interpreted.
*/        Absolute value is taken.
*i errorj.6537
                  ajres=abs(ajres)
*/
*ident upnea058
*/ acer  A. Trkov, October 2009
*/       Redo update to initialise the jxs filed
*/       (proposed as upiaea20 on 10aug2006 but got lost)
*/       and add initialisation of the xss field because it
*/       is used recursively to accumulate cross sections.
*/       Failure to initialise these fields may sometimes
*/       result in corrupted ace files when more than one
*/       are produced in the same NJOY run (e.g. generating
*/       fast and thermal ACE files in the same run or 
*/       making several ACE files at different temperatures).
*i acer.227
      common/xsst/xss(4000000),n3
      common/jxst/jxs(32)
*i acer.265
      n3=0
      do j=1,max3
        xss(j)=0
      end do
      do j=1,32
        jxs(j)=0
      end do
*/
*ident upnea059
*/ errorr A. Trkov, October 2009
*/        Replace relevant coding with J-sign dependent equivalent
*/        in csrmat of the reconr module
*d errorj.8350,8462
c
c           ***loop over possible channel spins
            kchanl=0
            idone=0
            do while (kchanl.lt.2.and.idone.eq.0)
               kchanl=kchanl+1
               inow=inowb
               kpstv=0
               kngtv=0
c              initialize matrix
               do j=1,3
                  do i=1,3
                     s(j,i)=0
                     r(j,i)=0
                  enddo
               enddo
c
c              ***loop over resonances
               inow=inow+6
               in=inow+nrs*6
               do i=1,nrs
                  aj=abs(a(inow+1))
c                 select only resonances with current j value
                  if (abs(aj-ajc).le.quar) then
                     if (a(inow+1).lt.zero) kngtv=kngtv+1
                     if (a(inow+1).gt.zero) kpstv=kpstv+1
                     iskip=0
                     if (kchanl.eq.1.and.a(inow+1).lt.zero) iskip=1
                     if (kchanl.eq.2.and.a(inow+1).gt.zero) iskip=1
                     if (iskip.eq.0) then
c                       retrieve parameters
                        er=a(inow)
                        gn=a(inow+2)
                        gg=a(inow+3)
                        gfa=a(inow+4)
                        gfb=a(inow+5)
                        per=a(in+1)
c                       gc=a(in+2)
                        a1=sqrt(gn*pe/per)
                        a2=0
                        if (gfa.ne.zero) a2=sqrt(abs(gfa))
                        if (gfa.lt.zero) a2=-a2
                        a3=0
                        if (gfb.ne.zero) a3=sqrt(abs(gfb))
                        if (gfb.lt.zero) a3=-a3
c                       compute energy factors
                        diff=er-e
                        den=diff*diff+quar*gg*gg
                        de2=haf*diff/den
                        gg4=quar*gg/den
c                       calculate r-function, or
c                       calculate upper triangular matrix terms
                        r(1,1)=r(1,1)+gg4*a1*a1
                        s(1,1)=s(1,1)-de2*a1*a1
                        if (gfa.ne.zero.or.gfb.ne.zero) then
                           r(1,2)=r(1,2)+gg4*a1*a2
                           s(1,2)=s(1,2)-de2*a1*a2
                           r(1,3)=r(1,3)+gg4*a1*a3
                           s(1,3)=s(1,3)-de2*a1*a3
                           r(2,2)=r(2,2)+gg4*a2*a2
                           s(2,2)=s(2,2)-de2*a2*a2
                           r(3,3)=r(3,3)+gg4*a3*a3
                           s(3,3)=s(3,3)-de2*a3*a3
                           r(2,3)=r(2,3)+gg4*a2*a3
                           s(2,3)=s(2,3)-de2*a2*a3
                           gf=1
                        endif
                     endif
                  endif
                  inow=inow+ncyc
                  in=in+3
               enddo
c
c              ***take care of extra channel spin as defined
c              ***by the sign of aj:
c              *** kkkkkk = 0 => do not add anything in here
c              *** kkkkkk = 1 => add resonance contribution but
c              ***  not extra hard-sphere
c              *** kkkkkk = 2 => add resonance plus hard-sphere
c              ***  phase shift contribution
               kkkkkk = 0
               if (kchanl.eq.1) then
                  if (kpstv.gt.0) then
                     if (kngtv.eq.0) then
                        if (jj.gt.jjl.and.jj.lt.numj) then
                           kkkkkk=2
                        else
                           kkkkkk=1
                        endif
                     else if (kngtv.gt.0) then
                        kkkkkk=1
                     endif
                  else if (kpstv.eq.0) then
                     if (kngtv.eq.0) then
                        if (jj.gt.jjl.and.jj.lt.numj) then
                           kkkkkk=2
                        else
                           kkkkkk=1
                        endif
                     else if (kngtv.gt.0) then
                        kkkkkk=0
                     endif
                  endif
               else if (kchanl.eq.2) then
                  if (kpstv.gt.0) then
                     if (kngtv.eq.0) then
                     else if (kngtv.gt.0) then
                        kkkkkk=1
                     endif
                  else if (kpstv.eq.0) then
                     if (kngtv.eq.0) then
                     else if (kngtv.gt.0) then
                        if (jj.gt.jjl.and.jj.lt.numj) then
                           kkkkkk=2
                        else
                           kkkkkk=1
                        endif
                     endif
                  endif
               endif
               if (kkkkkk.ne.0) then
c
c                 ***r-matrix path -- make symmetric matrix
                  if (gf.ne.zero) then
                     r(1,1)=uno+r(1,1)
                     r(2,2)=uno+r(2,2)
                     r(3,3)=uno+r(3,3)
                     r(2,1)=r(1,2)
                     s(2,1)=s(1,2)
                     r(3,1)=r(1,3)
                     s(3,1)=s(1,3)
                     r(3,2)=r(2,3)
                     s(3,2)=s(2,3)
c                    invert the complex matrix
                     call frobns(r,s,ri,si)
c                    fission term for r-matrix path
                     t1=ri(1,2)
                     t2=si(1,2)
                     t3=ri(1,3)
                     t4=si(1,3)
                     termf=four*gj*(t1*t1+t2*t2+t3*t3+t4*t4)
                     u11r=p1*(two*ri(1,1)-uno)+two*p2*si(1,1)
                     u11i=p2*(uno-two*ri(1,1))+two*p1*si(1,1)
                     termt=two*gj*(uno-u11r)
                     termn=gj*((uno-u11r)**2+u11i**2)
c
c                 ***r-function path
                  else
                     dd=r(1,1)
                     rr=uno+dd
                     ss=s(1,1)
                     amag=rr**2+ss**2
                     rri=rr/amag
                     ssi=-ss/amag
                     uur=p1*(two*rri-uno)+two*p2*ssi
                     uui=p2*(uno-two*rri)+two*p1*ssi
                     if (abs(dd).lt.small.and.
     &                 abs(phid).lt.small) then
                        xx=2*dd
                        xx=xx+2*(dd*dd+ss*ss+phid*phid+p2*ss)
                        xx=xx-2*phid*phid*(dd*dd+ss*ss)
                        xx=xx/amag
                        termt=two*gj*xx
                        termn=gj*(xx**2+uui**2)
                     else
                        termt=two*gj*(uno-uur)
                        termn=gj*((uno-uur)**2+uui**2)
                     endif
                     termf=0
                  endif
c
c                 ***cross sections contributions
                  if (kkkkkk.eq.2) then
                     termn=termn+two*gj*(1-p1)
                     termt=termt+two*gj*(1-p1)
                  endif
                  termg=termt-termf-termn
                  sigp(2)=sigp(2)+termn
                  sigp(4)=sigp(4)+termg
                  sigp(3)=sigp(3)+termf
                  sigp(1)=sigp(1)+termt
               endif
            enddo
*/
*ident upnea060
*/ errorr A. Trkov, October 2009
*/        A correction from D. Rochman was missing in upnea048
*d upnea048.61
         if(lrf.eq.3)then
            erap=a(l-5)
         else
            erap=a(l+1)
         endif
*/
*ident upnea061
*/ acer    A. Trkov, January 2010
*/         Correct errors in upnea027 (noticed by Skip Kahler)
*d acer.119,120
c     *               7   read fast ace files to print or edit         *
c     *               8   read thermal ace files to print or edit      *
c     *               9   read dosimetry ace files to print or edit    *
*d upnea027.108
      else if (iopt.eq.6.or.iopt.gt.9) then
*d upnea027.114
      else if (itype.eq.3) then
*/
*ident upnea062
*/ errorr  A. Trkov, Jan-2010
*/      1. Add MF40 processing capability
*/         WARNING:
*/         The coding was tested on the IRDF-2002 dosimetry library.
*/         If mfcov=40, the GENDF file must contain procesed mf=10 data.
*/         The procedure was tested for a single LFS state in MF10
*/         and a single covariance set for that state in MF40. 
*/         There may be more final states (e.g. in an activation library)
*/         but there are no suitable examples for testing.
*/      2. There is a change of convention in GENDF files generated by ERRORR
*/         Originally the C1 and C3 records were zero. For consistency
*/         with GENDF files produced by groupr, C1 was set to ZA and C2
*/         to 10*ZAP+LFS. This allows picking proper group data for MF40
*/         covariance processing.
*/      3. Correction of a trivial typing error correction in up307
*/         (statement beyond column 72)
*/
*/         NOTE: All testing was done with group data on ngout present.
*/               Further modifications might be needed if group data are
*/               to be generated internally.
*/
*d errorj.122
c     *    mfcov   endf covariance file (31, 33, 34, 35 or 40) to be   *
*i up272.28
     & ,mzap(80),lfs
*i up272.37   
     & ,mzap(80),lfs
*i up272.50   
     & ,mzap(80),lfs
*i up272.69   
     & ,mzap(80),lfs
*i up272.73   
     & ,mzap(80),lfs
*i up272.102  
     & ,mzap(80),lfs
*i up272.108  
     & ,mzap(80),lfs
*i up272.110  
     & ,mzap(80),lfs
      common/temper/tempin
*i up272.141  
     & ,mzap(80),lfs
*i up272.178  
     & ,mzap(80),lfs
*i up272.181  
     & ,mzap(80),lfs
*i up272.185  
     & ,mzap(80),lfs
*i up272.187  
     & ,mzap(80),lfs
*i up272.189  
     & ,mzap(80),lfs
*i up272.191  
     & ,mzap(80),lfs
*i up272.193  
     & ,mzap(80),lfs
*i up272.216  
     & ,mzap(80),lfs
*i up272.222  
     & ,mzap(80),lfs
*i up272.230  
     & ,mzap(80),lfs
*i up272.236  
     & ,mzap(80),lfs
*d errorj.424,425
         if(mfcov.ne.31.and.mfcov.ne.33.and.
     &      mfcov.ne.34.and.mfcov.ne.35.and.
     &      mfcov.ne.40) then
*i errorj.620
c
c     *** check if relevant covariance data are available
      nd=nw/6
      icov=0
      do i=1,nd
        mfi=nint(a(idict+2+(i-1)*6))
        if(mfcov.eq.30 .and. (mfi.ge.30 .and. mfi.le.33)) icov=icov+1
        if(mfcov.eq.31 .and.  mfi.eq.31)                  icov=icov+1
        if(mfcov.eq.32 .and.  mfi.eq.32)                  icov=icov+1
        if(mfcov.eq.33 .and. (mfi.eq.32 .or.  mfi.eq.33)) icov=icov+1
        if(mfcov.eq.34 .and.  mfi.eq.34)                  icov=icov+1
        if(mfcov.eq.35 .and.  mfi.eq.35)                  icov=icov+1
        if(mfcov.eq.40 .and.  mfi.eq.40)                  icov=icov+1
      end do
      if(icov.eq.0) then
         write(strng,'(''no data on file for mfcov='',i3)') mfcov
         call mess('errorr',strng,'processing terminated')
c        -- skip remaining errorr input (if any)
         if (ign.eq.1.or.ign.eq.19) then
            read(nsysi,*) ng
            ngp=ng+1
            read(nsysi,*) (dmy,i=1,ngp)
         end if
         go to 330
      end if
*i errorj.685
  330 continue
*i errorj.816  
      if (iverf.eq.4) then
         nl=l2h
      else
         nl=n2h
      end if
      izap=0
*d errorj.821,824
      elseif (mfcov.eq.40) then
         za=c1h
         awr=c2h
         nl=n1h
         call contio(nendf,0,0,a(iscr),nb,nw)
         izap=10*l1h+l2h
*i errorj.1480
      else if (mfcov.eq.40) then
         call sigc(ngn,a(isum),a(icflx),a(iscr),a(iun),a(iflx),a(isig))
*i errorj.1566 
         izap=mzap(ix)
*d errorj.1576
            a(iscr+1)=izap
*d errorj.1590
         a(iscr+1)=izap
*d errorj.1873
               if (mfcov.eq.31.or.mfcov.eq.33.or.mfcov.eq.35.or.
     &             mfcov.eq.40) then
*i errorj.3793 
         nfs =1
*i errorj.3795 
         nfs =1
         if(mfcov.eq.40) then
           nfs=n1h
           nsub=1
           call contio(nendf,0,0,a(iscr),nb,nw)
           lfs=l2h
           izap=10*l1h+l2h
         end if
*i errorj.3815
      mzap(nmt)=izap
*i errorj.3831
      do 310 ilfs=1,nfs
*i errorj.3978
  310 continue
*i errorj.1875
                  if (mfcov.eq.40) write(nsyso,38)lfs
c...              if (mfcov.eq.40) write(nsyso,38)lfs/10,lfs-10*(lfs/10)
*i errorj.2126 
   38 format(' final metastable state lfs',i3/)
c..38 format(' final metastable state zap,lfs',i6,i3/)
*/
*/ search gout tape by product identifier
*d errorj.2985
      subroutine rdgout(ngout,matd,mfd,mti,izap,b,sig)
*i errorj.3076 
      jzap=0
      if(izap.ne.0) jzap=c2h+0.01
*d errorj.3078
      if (mf.eq.mfd.and.mt.eq.mtd.and.jzap.eq.izap) go to 230
*d errorj.5320
      subroutine rdsig(mat,mt,izap,b,sig)
*i errorj.5329
      izero=0
*d errorj.5342
            call rdgout(ngout,matrd,mfri,mtri,izero,b,sig)
*d errorj.5345
         call rdgout(ngout,matrd,mfrd,mtrd,izap,b,sig)
*i errorj.734
      izero=0
*d errorj.752
      call rdgout(ngout,matd,mfd,mtd,izero,a(ib),a(iegt))
*d errorj.784
      call rdgout(ngout,matd,mfd,mtd,izero,a(ib),a(iflx))
*d errorj.786
      call rdsig(matd,izero,izero,a(ib),a(iscr))
*d errorj.874
         call rdsig(mat,mt,izap,a(ib),a(isig))
*d errorj.986
            if (mt1.lt.851) call rdsig(mat1,mt1,izero,a(ib),a(isig1))
*i errorj.3142 
      izero=0
*d errorj.3220 
            call rdsig(mats(ix),mtd,izero,b,a(isg))
*d errorj.3226 
  250    call rdsig(mats(ix),mts(ix),mzap(ix),b,sig)
*d errorj.3243,3244
         b(1)=za
         b(2)=mzap(ix)
*i errorj.5217
      izero=0
*d errorj.5246
         call rdsig(matd,mtd,izero,a(ib),a(iscr2))
*i errorj.7604
      izero=0
*d errorj.7617
  200 call rdsig(matd,mt1,izero,b,sig)
*d errorj.7628
      call rdsig(matd,mt2,izero,b,sig)
*i errorj.7782 
      izero=0
*d errorj.7841
      call rdsig(mat,mt,izero,b,sig)
*/
*/ add tmeperature to the output tape
*i errorj.3134
      common/temper/tempin
*d errorj.3163
         b(1)=tempin
*i errorj.3520 
      tempin=c1h
*/
*/ prevent printing "undefined" cross sections and covariances
*d errorj.1921
*i errorj.1922,1924
                     if(denom.gt.zero) then
                        denom=max(denom,eps)
                        a(iscr+ibase+ip-1)=a(iscr+ibase+ip-1)/denom*
     &                         (egn(ig+1)-egn(ig))*(egn(igp+1)-egn(igp))
                     else
                        a(iscr+ibase+ip-1)=0
                     end if
*d errorj.1926
                     if(denom.gt.zero) then
                        denom=max(denom,eps)
                        a(iscr+ibase+ip-1)=a(iscr+ibase+ip-1)/denom
                     else
                        a(iscr+ibase+ip-1)=0
                     end if
*d errorj.3235
            ff=cflx(ig)
            if(ff.gt.0) then
              csig(ig,ix)=csig(ig,ix)/ff
            else
              csig(ig,ix)=0
            end if
*/
*ident upnea063
*/ reconr  A. Trkov, November 2009
*/         When no MF3 data present (e.g. dosimetry library with MF10)
*/         make sure that the group structure (MF3 MT1)
*/         is written with zero cross sections and that N1, N2
*/         parameters from the CONT record are transferred consistently.
*i reconr.4560 
      ns1=n1h
      nr1=n2h
*i reconr.4561
      mf1=mfh
      if(mfh.eq.10) mfh=3
*i reconr.4623
      mfh=mf1
*i reconr.4635
      ns1=n1h
      nr1=n2h
      mf1=mfh
*i reconr.4617 
      if(mfh.ne.mf1) call afend(nout,0)
*d reconr.4655,4656
      a(iscr+4)=ns1
      a(iscr+5)=nr1
      mfh=mf1
*/
*ident upnea064
*/ acer    A. Trkov, February 2010
*/         Particle emission spectra in the acer plots are corrupted
*/         for MF6 Law 1 LANG 1 when more than one interpolation range
*/         is specified for the distributions on incident particle
*/         energies. In some cases the MCNP calculations are also wrong.
*/         A fix is done when writing the data to a temporary file.
*/         The multiple ranges are suppressed. The first interpolation
*/         law is prescribed over the entire incident energy range
*/         and a message is printed.
*/         WARNING:
*/         This is a temporary patch before a proper solution is found.
*/         The true error probably occurs somewhere near acer.6651
*/         or later, and/or possibly in MCNP.
*/         Implications:
*/         JENDL-HE files truncated to 150 MeV (as proposed for FENDL-3)
*/         led to strongly discrepant results in a benchmark model
*/         calculation representative of the ITER device. The discrepancy
*/         is greatly reduced with the use of this patch. The assumption
*/         in the patch has no influence on the evaluated data because
*/         the interpolation law changes above 20 MeV where the yield
*/         drops to zero due to the reaction representation in MT5.
*d acer.2373,2374
                  call tab2io(nin,0,0,b,nb,nw)
                  nr=n1h
                  ne=nint(b(6))
                  if(nr.gt.1) then
                    write(string,'(a,i3)')
     &                   'multiple interp. ranges for mf6, mt',mt
                    call mess('topfil',string
     &                    ,'first law applied everywhere')
                    nr=1
                    b(5)=nr
                    b(7)=ne
                  end if
                  call tab2io(0,nout,0,b,nb,nw)
*/