*cpl all
*set sw
*ident up1
*/ leapr -- 04apr00
*/ this fixes errors that destroy the calculation of s(alpha,beta)
*/ for cold para-hydrogen and for both cold deuteriums.  the first
*/ was introduced in 97.53 while making the constants symbolic,
*/ and the second was introduced in the cleanup process for njoy99.
*d leapr.1902
      if (law.gt.3) de=ded
*d leapr.2121
         sum2=(sjbes(jp,y)*cn(j,jp,jp))**2
*ident up2
*/ groupr -- 04apr00
*/ there is an error in the indexing for the xmas 172-group structure
*/ that throws all the group bounds off by one.  this structure is
*/ used in europe for thermal reactor calculations.
*d groupr.1835
            eg(ig)=eg18(174-ig)
*/ fix a problem introduced with the activation patch of 97.102.  the
*/ nk parameter is only used when doing mf values for activation
*/ products, and it is not appropriate for fission nubar values.
*d groupr.3824
      nk=0
      if (mf.gt.99) nk=nint(a(iyld+4))
*i groupr.3853
      if (nk.eq.0) go to 180
*/ don't strip off the upscatter groups for the neutron spectra
*/ coming from delayed neutron emission (mt=455).
*d groupr.8913
                  if (mfd.ne.15.and.mtd.ne.455) then
*/ add a missing save statement and fix an unset variable in anased.
*/ these problems affect delayed neutron spectra.
*d groupr.8983
      save new,theta,xc,rc,bot,ca,loct
*i groupr.9050
         np=nint(a(loct+6))
*ident up3
*/ acer -- 04apr00
*/ the acer consistency checks include an option to readjust
*/ eprime values that are greater than e, when appropriate.
*/ there are some problems with the logic, especially for
*/ cases that use histogram interpolation for the distribution.
*i acer.18017
                           ishift=j-nn-1
*d acer.18019
                           xss(j+loci)=sigfig(epmax,7,ishift)
*i acer.18022
                              xss(j-1+nn+loci)=p
                              xss(j+nn+loci)=p
*i acer.18026
                              xss(j+nn+loci)=p
*d acer.18028
*i acer.18131
                              xss(j-1+nn+loci)=p
                              xss(j+nn+loci)=p
*i acer.18135
                              xss(j+nn+loci)=p
*d acer.18137
*/ when using the old format (mcnp4b and earlier), some angle-energy
*/ distributions from file 6 are converted into the law 67 format,
*/ because these earlier versions of mcnp couldn't use all the 
*/ file 6 representations.  when converting from the cm to the lab, 
*/ the methods used in subroutine fix6 are a little crude.  they get
*/ confused when cm energies are so small that lab cosines of -1
*/ are not reached.  this patch tries to fix that in a rough way,
*/ but evaluations that use the cm frame in file 6 will work best
*/ if most of the cm energies are greater than e/(awr+1)**2 for
*/ each incident energy e.  for mcnp4c and later, the code can
*/ sample directly from tabulated cm representations, and the
*/ approximations of the fix6 routine are avoided.  this patch
*/ is needed for one evaluation from JEFF-3.
*d acer.3238
*i acer.3250
      data namax/1000/
*d acer.3353
                  if (lct.ne.1) then
*d acer.3372
                     if (ep.gt.zero) then
                        csn=clb*sqrt(elb/ep)-sqrt(ein/ep)/aw1
                     endif
*i acer.3398,3420
                  if (j.le.l2+8.or.elb.gt.a(j-2)) then
                     a(j)=elb
                     a(j+1)=fmu*drv
                     j=j+2
                  endif
                  if (j.ge.namax-1) call error('fix6',
     &              'storage in a exceeded',' ')
*i acer.3421
               nnep=(j-(l2+8))/2
               if (nnep.eq.1) then
                  a(l2+10)=2*a(l2+8)
                  a(l2+11)=0
                  nnep=2
               endif
               a(l2+5)=nnep
               a(l2+6)=nnep
               j2=l2
               call tab1io(0,nout,ndebug,a(j2),nb,nw)
               do while (nb.ne.0)
                  j2=j2+nw
                  call moreio(0,nout,ndebug,a(j2),nb,nw)
               enddo
*/ increase the available storage to handle the very large
*/ mf6/mt16 tabulation in JEFF-3 Be-9.
*d acer.226
      common/astore/a(80000)
*d acer.235
      data namax/80000/, nidmax/27/
*d acer.460
      common/astore/a(80000)
*d acer.2130
      data nwmaxn/65000/
*d acer.4672
      common/astore/a(80000)
*d acer.5604
      common/astore/a(80000)
*d acer.5765
      common/astore/a(80000)
*d acer.5954
      common/astore/a(80000)
*d acer.6326
      common/astore/a(80000)
*d acer.7385
      common/astore/a(80000)
*d acer.8058
      common/astore/a(80000)
*d acer.8068
      data namax/40000/
*d acer.9762
      common/astore/a(80000)
*d acer.10677
      common/astore/a(80000)
*d acer.13068
      common/astore/a(80000)
*d acer.13464
      common/astore/a(80000)
*d acer.14300
      common/astore/a(80000)
*d acer.14665
      common/astore/a(80000)
*d acer.15215
      common/astore/a(80000)
*d acer.21814
      common/astore/a(80000)
*/ increase the space available for discontinuities in convr
*/ to allow for JENDL-3.2 si-nat
*d acer.254
      nned=50
*ident up4
*/ reconr -- 05apr00
*/ be sure to count subsections of file 12 before allocating
*/ storage for the elements of the new directory.  otherwise,
*/ some materials with many sections of file 12 will overflow.
*/ this is a longstanding problem that we never noticed before.
*i reconr.419
            nxn=nxn+1
*ident up5
*/ purr -- 7may00
*/ fix a problem introduced while installing the heating part
*/ of the probability tables.  it shows up when doing elements
*/ that have unresolved data.  also, increase the number of
*/ resonance sections allowed to handle the very large cd-nat
*/ evaluation from JENDL.
*d purr.1076
                     e=abs(eunr(ie))
*d purr.1106,1108
      common/sigcon/e,t,cth(50),csz(50),cc2p(50),cs2p(50), 
     &  cgn(50),cgg(50),cgf(50),cgx(50),cgt(50),dbar(50),
     &  spot,dbarin,sigi(4),ndfn(50),ndff(50),ndfx(50),nseq0
*d purr.1139
*d purr.1187,1189
*d purr.1247
                  if (nseq0.gt.50) call error('unresx',
*d purr.1501,1503
      common/sigcon/e,t,cth(50),csz(50),cc2p(50),cs2p(50),
     &  cgn(50),cgg(50),cgf(50),cgx(50),cgt(50),dbar(50),
     &  spot,dbarin,sigi(4),ndfn(50),ndff(50),ndfx(50),nseqz
*d purr.1621,1623
      common/sigcon/e,t,cth(50),csz(50),cc2p(50),cs2p(50),
     &  cgn(50),cgg(50),cgf(50),cgx(50),cgt(50),dbar(50),
     &  spot,dbarin,sigi(4),ndfn(50),ndff(50),ndfx(50),nseqz
*d purr.1739,1741
      common/sigcon/e,t,cth(50),csz(50),cc2p(50),cs2p(50),
     &  cgn(50),cgg(50),cgf(50),cgx(50),cgt(50),dbar(50),
     &  spot,dbarin,sigi(4),ndfn(50),ndff(50),ndfx(50),nseq0
*ident up6
*/ acer -- 30may00
*/ fix a typo in up3 (reported by bunde, anl)
*d up3.93
      data namax/80000/
*/ acer -- 30may00
*/   fix problems with converting cm distributions to law=7
*/   and problems reading law=7 into the ace file.  these problems
*/   show up when running newfor=0 with njoy2000, especially on
*/   some materials from jef-2.2.
*d acer.3342
*d acer.3347,3348
*d acer.3364
c                       if(imu.lt.nmu.and.amu(imu+1).le.cmn) drv=0
*i acer.3377
c                    include jacobian for cm-to-lab transformation
                     if (ep.ne.zero) drv=drv*sqrt(elb/ep)
*d up3.46,up3.50
                  if (j.le.l2+8) then
                     a(j)=elb
                     a(j+1)=fmu*drv
                     j=j+2
                  else if (elb.gt.a(j-2)) then
                     a(j)=elb
                     a(j+1)=fmu*drv
                     j=j+2
                  endif     
*d acer.3399,3420
                  if (iep.eq.nep) idone=1
*d acer.6330
      external listio,terpa,terp1,bachaa,mess,fndar1,fndar2,skip6a
*d acer.6380
               call skip6a(nin,0,0,a(jscr),law)
*d acer.6412
               call skip6a(nin,0,0,a(jscr),law)
*d acer.6497
      call skip6a(nin,0,0,a(jscr),law)
*i acer.6935
c
      subroutine skip6a(nin,nout,nscr,a,law)
c     ******************************************************************
c     special version of skip6 for special version of File 6 used
c     in acer.  law=7 has a tab1 containing the angular distribution
c     instead of the normal tab2 for each incident energy.
c     skip the next subsection in the current section (mt).
c     ******************************************************************
*if sw
      implicit real*8 (a-h,o-z)
*endif
      common/cont/c1h,c2h,l1h,l2h,n1h,n2h,math,mfh,mth,nsh,nsp,nsc
      dimension a(*)
c
      if (law.eq.6) then
         call contio(nin,nout,nscr,a(1),nb,nw)
      else if (law.eq.1.or.law.eq.2.or.law.eq.5) then
         call tab2io(nin,nout,nscr,a(1),nb,nw)
         ne=n2h
         do ie=1,ne
            call listio(nin,nout,nscr,a(1),nb,nw)
            do while (nb.ne.0)
               call moreio(nin,nout,nscr,a(1),nb,nw)
            enddo
         enddo
      else if (law.eq.7) then
         call tab2io(nin,nout,nscr,a(1),nb,nw)
         ne=n2h
         do ie=1,ne
            call tab1io(nin,nout,nscr,a(1),nb,nw)
            nmu=n2h
            do imu=1,nmu
               call tab1io(nin,nout,nscr,a(1),nb,nw)
               do while (nb.ne.0)
                  call moreio(nin,nout,nscr,a(1),nb,nw)
               enddo
            enddo
         enddo
      endif
      return
      end
*ident up7
*/ viewr -- 30may00
*/ increase the allowed size of 3d plots.
*/ pushed by pb-nat from jef-2.2.
*d viewr.680
                          if (l+5000.ge.maxaa) then
*d viewr.1295
      dimension x(2000),y(2000),z(2000)
*d viewr.1304
      kmax=1999
*ident up8
*/ acer -- 3jun00
*/ subroutine ptleg2 does not need the dynamic array xat.
*/ this problem was first noted by Waclaw Gudowski for ENDF/B-VI 
*/ si-nat.  it shows up as "id xat not defined".
*d acer.5628
*d acer.5632
                        call ptleg2(a(iscr))
*d acer.5646
                     call ptleg2(a(iscr))
*d acer.6838,6839
                     call ptleg2(a(jscr))
*d acer.6937
      subroutine ptleg2(a)
*d acer.6950
*d acer.8470,8471
                        call ptleg2(a(lld))
*d acer.8720
*d acer.8731
                              call ptleg2(a(lld)) 
*d acer.16034
*d acer.16037
                              call ptleg2(a(lld))
*ident up9
*/ acer -- 07jun00
*/  add the capability for processing anisotropic charged particle
*/  emission using tabulated legendre coefficients into the 
*/  mcnp4c law61 format.  this is needed for jeff-3 cr-52.
*/  allow for multiple interpolation ranges in file 6.  this
*/  also occurs for jeff-3 cr-52.  currently, the neutron
*/  energy-angle distribution only allows for combinations of
*/  histogram and linear linear interpolation, but the 
*/  charged-particle sections allow for general combinations of 
*/  all allowed interpolation laws.
*d acer.6391
               jnt=nint(a(jscr+5+2*m))
*i acer.6456
                  if (idis.eq.1.and.xx.lt..9999*xn) xn=sigfig(xn,7,-1)
*d acer.8278,8285
                        next=next+2
                        nrint=nint(a(iscr+4))
                        if (nrint.eq.1.and.nint(a(iscr+7)).eq.2) then
                           xss(next)=0
                        else
                           xss(next)=nrint
                           do i=1,nrint
                              xss(next+i)=nint(a(iscr+4+2*i))
                              xss(next+nrint+i)=nint(a(iscr+5+2*i))
                           enddo
                           next=next+2*nrint
                        endif
                        next=next+1
                        ne=nint(a(iscr+5))
                        xss(next)=ne
                        do i=1,ne
                           xss(next+i)=
     &                       sigfig(a(iscr+4+2*nrint+2*i)/emev,7,0)
                           xss(next+i+ne)=
     &                       sigfig(a(iscr+5+2*nrint+2*i),7,0)
                        enddo
                        next=next+1+2*ne
*d acer.9031,9032
*i acer.9033
                        lang=nint(a(ll+2))
                        lawnow=0
                        if (law.eq.1.and.lang.eq.1) lawnow=61
                        if (law.eq.1.and.lang.eq.2) lawnow=44
                        if (law.eq.2) lawnow=33
                        if (lawnow.eq.0) call error('acelcp',
     &                    'unsupported law and lang',' ')
                        xss(last+1)=lawnow
*i acer.9090
                              nexcd=next+4*ng+2
*d acer.9121
c                                kalbach distribution
                                 if (lang.eq.2) then
*i acer.9126
c                                legendre distribution
                                 else if (lang.eq.1) then
                                    ep=xss(next+1+ig)
                                    a(ll)=0
                                    a(ll+1)=ep
                                    a(ll+2)=0
                                    a(ll+3)=0
                                    a(ll+4)=na
                                    a(ll+5)=0
                                    do ia=1,na
                                       lll=lld+7+ncyc*(ig-1)
                                       a(ll+5+ia)=0
                                       if (a(lll).ne.zero) then
                                          a(ll+5+ia)=a(lll+ia)/a(lll)
                                       endif
                                    enddo
                                    call ptleg2(a(ll))
                                    xss(next+1+3*ng+ig)=nexcd-dlwh+1
                                    intmu=2
                                    xss(nexcd)=intmu
                                    nmu=nint(a(ll+5))
                                    xss(nexcd+1)=nmu
                                    do imu=1,nmu
                                       xss(nexcd+1+imu)=sigfig(
     &                                   a(ll+6+2*imu),7,0)
                                       xss(nexcd+1+nmu+imu)=sigfig(
     &                                   a(ll+7+2*imu),7,0)
                                       if (imu.eq.1) then
                                          xss(nexcd+1+2*nmu+imu)=0
                                       else if (imu.eq.nmu) then
                                          xss(nexcd+1+2*nmu+imu)=1
                                       else
                                          del=a(ll+6+2*imu)
     &                                      -a(ll+4+2*imu)
                                          av=(a(ll+7+2*imu)
     &                                      +a(ll+5+2*imu))/2
                                          xss(nexcd+1+2*nmu+imu)
     &                                      =xss(nexcd+1+2*nmu+imu-1)
     &                                      +del*av
                                          xss(nexcd+1+2*nmu+imu)=sigfig
     &                                      (xss(nexcd+1+2*nmu+imu),7,0)
                                       endif
                                    enddo
                                    nexcd=nexcd+2+3*nmu
*d acer.9160
                              if (lang.eq.1) then
                                 next=nexcd
                              else
                                 next=next+2+(2*na+3)*ng 
                              endif
*d acer.11158,11163
            l2=sigh+l1+1
            nrint=nint(xss(l2))
            write(nsyso,'(4x,''nr ='',i4)') nrint
            if (nrint.ne.0) then
               write(nsyso,'(4x,''nbt(i=1,nr) = '',20i5)')
     &           (nint(xss(l2+ii)),ii=1,nrint)
               l2=l2+nrint
               write(nsyso,'(4x,''int(i=1,nr) = '',20i5)')
     &           (nint(xss(l2+ii)),ii=1,nrint)
               l2=l2+nrint
            endif
            l2=l2+1
            ne=nint(xss(l2))
            write(nsyso,'(4x,''ne ='',i4)') ne
*d acer.11169
     &         xss(l2+ii),xss(l2+ne+ii)
*d acer.11255,11257
                  write(nsyso,'(4x,''nr ='',i4)') nrint
                  if (nrint.ne.0) then
                     write(nsyso,'(4x,''nbt(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                     write(nsyso,'(4x,''int(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                  endif
                  l3=l3+1
                  ne=nint(xss(l3))
                  write(nsyso,'(4x,''ne ='',i4)') ne
*d acer.11259,11260
                     e2=xss(l3+ie)
                     loci=nint(xss(l3+ne+ie))+dlwh-1
*d acer.11286,11288
                  write(nsyso,'(4x,''nr ='',i4)') nrint
                  if (nrint.ne.0) then
                     write(nsyso,'(4x,''nbt(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                     write(nsyso,'(4x,''int(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                  endif
                  l3=l3+1
                  ne=nint(xss(l3))
                  write(nsyso,'(4x,''ne ='',i4)') ne
*d acer.11290,11291
                     e2=xss(l3+ie)
                     loci=nint(xss(l3+ne+ie))+dlwh-1
*i acer.11311
c
c              ***law=61
               else if (law.eq.61) then
                  nrint=nint(xss(l3))
                  write(nsyso,'(4x,''nr ='',i4)') nrint
                  if (nrint.ne.0) then
                     write(nsyso,'(4x,''nbt(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                     write(nsyso,'(4x,''int(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                  endif
                  l3=l3+1
                  ne=nint(xss(l3))
                  write(nsyso,'(4x,''ne ='',i4)') ne
                  do ie=1,ne
                     e2=xss(ie+l3)
                     loci=nint(xss(ie+ne+l3)+dlwh-1)
                     intt=mod(nint(xss(loci)),10)
                     nd=nint(xss(loci)/10)
                     nn=nint(xss(loci+1))
                     loci=loci+1
                     write(nsyso,'(/6x,'' incident energy = '',
     &                 1p,e14.6,''   intt ='',i2,''    nd ='',i4,
     &                 ''    np ='',i3)') e2,intt,nd,nn
                     do ip=1,nn
                        locj=nint(xss(ip+3*nn+loci)+dlwh-1)
                        intmu=nint(xss(locj))
                        nmu=nint(xss(locj+1))
                        write(nsyso,'(/
     &                    6x,'' secondary energy = '',1p,e14.6/
     &                    6x,''              pdf = '',e14.6/
     &                    6x,''              cdf = '',e14.6/
     &                    6x,''            intmu = '',i8/
     &                    6x,''              nmu = '',i8/
     &               ''         cosine           pdf           cdf'',
     &               ''        cosine           pdf           cdf''/
     &               ''   ------------  ------------  ------------'',
     &               ''  ------------  ------------  ------------'')')
     &                    xss(ip+loci),xss(ip+nn+loci),
     &                    xss(ip+2*nn+loci),intmu,nmu
                        do imu=1,nmu,2
                           if (imu.eq.nmu) then
                              write(nsyso,'(1x,1p,3e14.6)')
     &                          xss(locj+1+imu),xss(locj+1+nmu+imu),
     &                          xss(locj+1+2*nmu+imu)
                           else
                              write(nsyso,'(1x,1p,6e14.6)')
     &                          xss(locj+1+imu),xss(locj+1+nmu+imu),
     &                          xss(locj+1+2*nmu+imu),xss(locj+1+imu+1),
     &                          xss(locj+1+nmu+imu+1),
     &                          xss(locj+1+2*nmu+imu+1)
                           endif
                        enddo
                     enddo
                  enddo
*d acer.11333,11335
                  write(nsyso,'(4x,''nr ='',i4)') nrint
                  if (nrint.ne.0) then
                     write(nsyso,'(4x,''nbt(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                     write(nsyso,'(4x,''int(i=1,nr) = '',20i5)')
     &                 (nint(xss(l3+ii)),ii=1,nrint)
                     l3=l3+nrint
                  endif
                  l3=l3+1
                  ne=nint(xss(l3))
                  write(nsyso,'(4x,''ne ='',i4)') ne
*d acer.11337,11339
                     e2=xss(l3+ie)
                     loci=nint(xss(l3+ne+ie))+dlwh-1
                     intmu=nint(xss(loci))
*i acer.12617
                  if (nout.ne.1) nr=nint(xss(l))
                  if (nout.eq.1) nr=iss(l)
*i acer.12619
                  if (nr.gt.0) then
                     n=2*nr
                     do jj=1,n
                        call typen(l,nout,1)
                        l=l+1
                     enddo
                  endif
*i acer.12791
                     else if (law.eq.61) then
                        if (nout.ne.1) nr=nint(xss(l))
                        if (nout.eq.1) nr=iss(l)
                        call typen(l,nout,1) 
                        l=l+1
                        if (nr.ne.0) then
                           n=2*nr
                           do j=1,n
                              call typen(l,nout,1)
                              l=l+1
                           enddo
                        endif
                        if (nout.ne.1) ne=nint(xss(l))
                        if (nout.eq.1) ne=iss(l)
                        call typen(l,nout,1)
                        l=l+1
                        do j=1,ne
                           call typen(l,nout,2)
                           l=l+1
                        enddo
                        do j=1,ne
                           call typen(l,nout,1)
                           l=l+1
                        enddo
                        do j=1,ne
                           call typen(l,nout,1)
                           l=l+1
                           if (nout.ne.1) np=nint(xss(l))
                           if (nout.eq.1) np=iss(l)
                           call typen(l,nout,1)
                           l=l+1
                           n=3*np
                           do k=1,n
                              call typen(l,nout,2)
                              l=l+1
                           enddo
                           do k=1,np
                              call typen(l,nout,1)
                              l=l+1
                           enddo
                           do k=1,np
                              call typen(l,nout,1)
                              l=l+1
                              if (nout.ne.1) nmu=nint(xss(l))
                              if (nout.eq.1) nmu=iss(l)
                              call typen(l,nout,1)
                              l=l+1
                              nw=3*nmu
                              do kk=1,nw
                                 call typen(l,nout,2)
                                 l=l+1
                              enddo
                           enddo
                        enddo                  
*d acer.18174
                        locj=nint(xss(j+3*nn+loci)+dlw-1)
*d acer.18179
                           cc=xss(locj+1+2*nmu+k)
*d acer.18353
                  j=nint(xss(l3))
                  if (j.ne.0) then
                     l3=l3+2*j
                  endif
                  l3=l3+1
                  ne=nint(xss(l3))
*d acer.18355,18356
                     e=xss(l3+ie)
                     loci=nint(xss(l3+ne+ie))+dlwh-1
*d acer.18384
                  j=nint(xss(l3))
                  if (j.ne.0) then
                     l3=l3+2*j
                  endif
                  l3=l3+1
                  ne=nint(xss(l3))
*d acer.18386,18387
                     e=xss(l3+ie)   
                     loci=nint(xss(l3+ne+ie))+dlwh-1
*d acer.18424
                  j=nint(xss(l3))
                  if (j.ne.0) then
                     l3=l3+2*j
                  endif
                  l3=l3+1
                  ne=nint(xss(l3))
*d acer.18426,18427
                     e=xss(l3+ie)   
                     loci=nint(xss(l3+ne+ie))+dlwh-1
*d acer.18449
                        locj=nint(xss(j+3*nn+loci)+dlwh-1)
*d acer.18454
                          cc=xss(locj+1+2*nmu+k)
*d acer.18459
     &                          '' at'',1p,e14.6,'' ->'',e13.6,e14.6)')
*ident up10
*/ leapr -- 13jun00
*/  fix two incorrect constants in leapr.  one affects cases with
*/  diffusive effects, and it has been incorrect since njoy97.0 
*/  (oct 97).  the other affects cold hydrogen calculations, and it 
*/  has been incorrect since njoy97.53 (dec98).
*d leapr.1186
      data c0/.125d0/
*d leapr.1864
      data amassh/3.3465d-24/
*ident up11
*/ acer -- 26jun00
*/  fix an error in determinining which reactions have to be
*/  converted into law=7 format when using newfor=0.  because of 
*/  overzealous code cleanup, acer is trying to convert sections 
*/  with the kalbach representation in addition to sections with
*/  tabulated angular distributions.
*d acer.2324,2330
                  do while(nb.ne.0)
                     call moreio(nin,0,0,a(iscr),nb,nw)
                  enddo
                  if (lf.eq.6) then
                     call contio(nin,0,0,a(iscr),nb,nw)
                  else if (lf.eq.1.or.lf.eq.2.or.lf.eq.5) then
                     call tab2io(nin,0,0,a(iscr),nb,nw)  
                     lang=l1h
                     if (dzap.le.test.and.lf.eq.1.and.lang.ne.2) new6=1
                     ne=n2h
                     do ie=1,ne
                        call listio(nin,0,0,a(iscr),nb,nw)
                        do while (nb.ne.0)
                           call moreio(nin,0,0,a(iscr),nb,nw)
                        enddo
                     enddo
                  else if (lf.eq.7) then
                     call tab2io(nin,0,0,a(iscr),nb,nw)  
                     ne=n2h
                     do ie=1,ne
                        call tab1io(nin,0,0,a(iscr),nb,nw)
                        nmu=n2h
                        do imu=1,nmu
                           call tab1io(nin,0,0,a(iscr),nb,nw)
                           do while (nb.ne.0)
                              call moreio(nin,0,0,a(iscr),nb,nw)
                           enddo
                        enddo
                     enddo
                  endif
*/ acer -- 26jun00
*/  during the cleanup of the topfil routine, the logic to process
*/  sections of file 6 using law=2 (two-body distributions) into
*/  equally probable bins for newfor=0 was omitted.  this shows up
*/  for evaluations that use mf6/mt51, etc., to represent inelastic
*/  levels.
*d acer.2379
                     if (lf.eq.1) then
*d acer.2387,2388
c                    law=2 for newfor=1 - copy the subsection
                     else if (lf.eq.2.and.newfor.eq.1) then
*i acer.2394
c                    law=2 for newfor=0 - convert to probability bins
                     else if (lf.eq.2.and.newfor.eq.0) then
                        call listio(nin,0,0,a(iscr),nb,nw)
                        now=iscr+nw
                        do while (nb.ne.0)
                           call moreio(nin,0,0,a(now),nb,nw)
                           now=now+nw
                        enddo
                        now=now-1
                        lang=nint(a(iscr+2))
                        if (lang.eq.0) then
c                          legendre coefficients
                           call ptleg(nout,a)
                        else
c                          tabulated angular distribution
                           do i=iscr,now
                              a(now+2-i+iscr)=a(now-1+iscr)
                           enddo
                           np=nint(a(iscr+7))
                           a(iscr)=a(iscr+2)
                           a(iscr+1)=a(iscr+3)
                           a(iscr+2)=0
                           a(iscr+3)=0
                           a(iscr+4)=1
                           a(iscr+5)=np
                           a(iscr+6)=np
                           a(iscr+7)=lang-10
                           call pttab(ltt,a(iscr),nout)
                        endif
c                    law=5 for - copy the subsection
                     else if (lf.eq.5) then
                        call listio(nin,nout,0,a(iscr),nb,nw)
                        now=iscr+nw
                        do while (nb.ne.0)
                           call moreio(nin,nout,0,a(now),nb,nw)
                           now=now+nw
                        enddo
*/ acer -- 26jun00
*/  fix an error in the reaction naming.  it affects mt=44 (n,n2p) and
*/  mt=45 (n,npa).  this problem was introduced in may of 1995.
*/  examples of cases that use these reactions are later releases of 
*/  endf/b-vi al-27.
*d acer.11415,11416
     &  '(n,x)     ', '(n,2np)   ', '(n,3np)   ', '(n,x)     ',
     &  '(n,n2p)   ', '(n,npa)   ', '(n,2/2*1) ', '(n,2/2*2) ',
*/ acer -- 26jun00
*/  add missing external statement.  reported by bokyun seo (kaeri)
*i acer.3250
      external error
*/ acer -- 26jun00
*/  add missing line for the sequential (n,2n) reactions for be-9.
*/  this line was accidentally removed in njoy 94.19 (jan96).  the
*/  error was continued through njoy 97 and 99.  discovered by
*/  bob little (lanl).
*i acer.5102
                  if (mth.ge.46.and.mth.le.49) s=sigfig(s/2,7,0)
*/ acer -- 26jun00
*/  fix anisotropic photon production (e.g., endf c,n,o)
*d acer.7559
         if (lff.le.1) then
*ident up12
*/ njoy -- 27jul00
*/  fix two typographical errors in the 64-bit version of the
*/  slatec math library.  reported by piet de leege (delft).
*d njoy.4617
         if (a.ge.(-0.5).or.aeps.ne.0.0) then
*d njoy.4935
         gamma=0.9375+csevl(2.*y-1.,gamcs,ngamcs)
*ident up13
*/ reconr -- 12jul00
*/  if a reaction uses histogram interpolation, reconr tries to
*/  change it to linear interpolation by moving each point down by 
*/  one in the seventh place and and adding a point higher by one in 
*/  the seventh place.  if there is already a point in the evaluator's 
*/  grid higher by one in the seventh place, the algorithm gets 
*/  confused.  this currently occurs for carbon from release 6 of 
*/  endf/b-vi.  the symptom is an infinite loop while processing 
*/  mf=12,mt=51.  we found this at los alamos, and skipped over the 
*/  problem by temporarily patching the evaluation.  more recently, 
*/  it was re-reported by waclaw gudowski, and now we are making a 
*/  real fix for the problem.
*i reconr.1830
      if (er.lt.(1+small)*enl) go to 255
*ident up14
*/ acer -- 20jul00
*/  acer fails if you run it on a pendf tape that only has the
*/  single reaction mt=2 (elastic).  this can happen for he-4 if you
*/  don't run heatr, thermr, or gaspr first.  found by gudowski.
*d acer.5121
      mt=2
*/ acer -- 20jul00
*/  acer fails for mf=6, law=2, lang>0 (angular distribution with
*/  tabulated cosines).  the only known example is mt=51 for pb-208 from 
*/  release 6 of endf/b-vi.  found by waclaw gudowski.
*d up11.67
                              a(now+2-i+iscr)=a(now-i+iscr)
*/ acer -- 21jul00
*/  an error was included in up9, which added a capability to handle
*/  anisotropic charged-particle emission represented using legendre
*/  polynomials.  the update disabled the case of isotropic
*/  charged-particle emission, which occurs in a number of important
*/  materials from release 6 of endf/b-vi.  the symptom is a serious
*/  clobbering of the ace file, such that it cannot even be read into
*/  mcnp or even printed using acer.  also reported by gudowski.
*d up9.56
                                 else if (lang.eq.1.and.na.gt.0) then
*d up9.100
                              if (lang.eq.1.and.na.gt.0) then
*ident up15
*/ heatr -- 31jul00
*/  incorrect initial value found by m.mattes (u.stuttgart).
*d heatr.2586
         ir=1
*/  increase the allowed number of legendre terms in h6ddx
*/  to handle the new jeff-3t fe-56 evaluation.
*d heatr.3284
      dimension cnow(*),p(15)
*i heatr.3292
      data nlmax/15/
*i heatr.3315
      if (nl.gt.nlmax) call error('h6ddx',
     &  'too many legendre terms',' ')
*/ watch for ill-defined vertical segments in distributions.  these
*/ have been seen in zr90 from cendl3 and fe56 from jeff3.  actually,
*/ the evaluations should be fixed to avoid such features, because
*/ we don't really know what y value to select in the vertical
*/ segment.  we choose to just move the second energy of the double
*/ point up a little.  we only print the diagnostic once to keep the
*/ output cleaner, but there could be more than one vertical segment.
*i heatr.3286
      external mess
*i heatr.3287
      save illdef
*i heatr.3318
         illdef=0
*d heatr.3352,3353
               x1=cnow(lnow-ncnow)
               x2=cnow(lnow)
               if (x1.eq.x2.and.lep.gt.1) then
                  x2=sigfig(x2,6,1)
                  if (illdef.eq.0) then
                     call mess('h6ddx',
     &                 'vertical segment(s) in distribution',
     &                 'y(x) is ill defined')
                     illdef=1
                  endif
               endif
               y1=cnow(lll-ncnow)
               y2=cnow(lll)
               call terp1(x1,y1,x2,y2,ep,tt,lep)
*d heatr.3364,3367
            x1=cnow(lnow-ncnow)
            x2=cnow(lnow)
            if (x1.eq.x2.and.lep.gt.1) then
               x2=sigfig(x2,6,1)
               if (illdef.eq.0) then
                  call mess('h6ddx',
     &              'vertical segment(s) in distribution',
     &              'y(x) is ill defined')
                  illdef=1
               endif
            endif
            y1=cnow(lnow-ncnow+1)
            y2=cnow(lnow+1)
            call terp1(x1,y1,x2,y2,ep,s,lep)
            y1=cnow(lnow-ncnow+2)
            y2=cnow(lnow+2)
            call terp1(x1,y1,x2,y2,ep,r,lep)
*d heatr.3380
         x1=cnow(ii)
         x2=cnow(jj)
         if (x1.eq.x2.and.lep.gt.1) then
            x2=sigfig(x2,6,1)
            if (illdef.eq.0) then
               call mess('h6ddx',
     &           'vertical segment(s) in distribution',
     &           'y(x) is ill defined')
               illdef=1
            endif
         endif
         y1=cnow(ii+1)
         y2=cnow(jj+1)
         call terp1(x1,y1,x2,y2,ep,s,lep)
*d heatr.3395
         call terp1(x1,tii,x2,tjj,ep,t,lep)
*ident up16
*/ groupr -- 31jul00
*/ watch for ill-defined vertical segments in distributions.  these
*/ have been seen in zr90 from cendl3 and fe56 from jeff3.  actually,
*/ the evaluations should be fixed to avoid such features, because
*/ we don't really know what y value to select in the vertical
*/ segment.  we choose to just move the second energy of the double
*/ point up a little.  we only print the diagnostic once to keep the
*/ output cleaner, but there could be more than one vertical segment.
*i groupr.5588
      save illdef
*i groupr.5591
      external mess
*i groupr.5633
         illdef=0
*d groupr.5675,5676
                  x1=cnow(lnow-ncnow)
                  x2=cnow(lnow)
                  if (x1.eq.x2.and.lep.gt.1) then
                     x2=sigfig(x2,6,1)
                     if (illdef.eq.0) then
                        call mess('h6ddx',
     &                    'vertical segment(s) in distribution',
     &                    'y(x) is ill defined')
                        illdef=1
                     endif
                  endif
                  y1=cnow(lll-ncnow)
                  y2=cnow(lll)
                  call terp1(x1,y1,x2,y2,ep,tt,lep)
*d groupr.5688,5691
               x1=cnow(lnow-ncnow)
               x2=cnow(lnow)
               if (x1.eq.x2.and.lep.gt.1) then
                  x2=sigfig(x2,6,1)
                  if (illdef.eq.0) then
                     call mess('h6ddx',
     &                 'vertical segment(s) in distribution',
     &                 'y(x) is ill defined')
                     illdef=1
                  endif
               endif
               y1=cnow(lnow-ncnow+1)
               y2=cnow(lnow+1)
               call terp1(x1,y1,x2,y2,ep,s,lep)
               y1=cnow(lnow-ncnow+2)
               y2=cnow(lnow+2)
               call terp1(x1,y1,x2,y2,ep,r,lep)
*d groupr.5722
            x1=cnow(ii)
            x2=cnow(jj)
            if (x1.eq.x2.and.lep.gt.1) then
               x2=sigfig(x2,6,1)
               if (illdef.eq.0) then
                  call mess('h6ddx',
     &              'vertical segment(s) in distribution',
     &              'y(x) is ill defined')
                  illdef=1
               endif
            endif
            y1=cnow(ii+1)
            y2=cnow(jj+1)
            call terp1(x1,y1,x2,y2,ep,s,lep)
*d groupr.5737
            call terp1(x1,tii,x2,tjj,ep,t,lep)
*ident up17
*/ acer -- 31jul00
*/ fix an incorrect index in the law=61 section for the primary
*/ particle.  the effect of this error is to give an incorrect
*/ angular distribution for the energy points with scattering
*/ probability zero (which should be isotropic).  this change
*/ is strictly cosmetic and shouldn't affect any results.
*d acer.6833
                        a(jscr+5+ia)=0
*/ acer -- 31jul00
*/ we want to use the compact law=4 for isotropic charged particle
*/ distributions, and the more general law=61 for anisotropic cp
*/ distributions.  unfortunately, we can't tell which is which
*/ without reading past the first few energies for some evaluations.
*/ the need for this patch was first noticed by jeff3 fe56.
*i acer.8982
c
c                 ***first check the subsection to see whether
c                 ***the distribution is isotropic or not.
                  isocp=1
                  call findf(matd,mf,mt,nin)
                  call contio(nin,0,0,a(iscr),nb,nw)
                  nk=n1h
                  ik=0
                  idone=0
                  do while (ik.lt.nk.and.idone.eq.0)
                     ik=ik+1
                     ll=iscr
                     lly=ll
                     call tab1io(nin,0,0,a(ll),nb,nw)
                     izap=nint(c1h)
                     awp=c2h
                     law=l2h
                     ll=ll+nw
                     do while (nb.ne.0)
                        call moreio(nin,0,0,a(ll),nb,nw)
                        ll=ll+nw
                     enddo
c
c                    ***if not the desired particle, skip the subsection
                     if (izap.ne.ip) then
                        call skip6(nin,0,0,a(iscr),law)
c
c                    ***we only need to check law 1 subsections
                     else if (law.eq.1) theN
                        call tab2io(nin,0,0,a(ll),nb,nw)
                        lang=nint(a(ll+2))
                        lep=nint(a(ll+3))
                        ne=nint(a(ll+5))
                        do ie=1,ne
                           ll=lld
                           call listio(nin,0,0,a(ll),nb,nw)
                           ll=ll+nw
                           do while (nb.ne.0)
                              call moreio(nin,0,0,a(ll),nb,nw)
                              ll=ll+nw
                           enddo
                           na=nint(a(lld+3))
                           if (na.gt.0) isocp=0
                        enddo
                     endif
                  enddo
c
c                 ***go back and process the subsection
*d up9.43
                        if (law.eq.1.and.lang.eq.1.and.isocp.eq.0)
     &                    lawnow=61
                        if (law.eq.1.and.lang.eq.1.and.isocp.eq.1)
     &                    lawnow=4
*i up9.48
                        if (law.eq.1.and.isocp.eq.1) xss(landh+jp-1)=0 
*d acer.9079,9082
*d acer.9085
                              if (lawnow.eq.4) then
*d up14.23
                                 else if (lawnow.eq.61) then
*d up14.25
                              if (lawnow.eq.61) then
*ident up18
*/ acer -- 03aug00
*/ some jef, eff, and jeff evaluations contain a redundant reaction
*/ mt=10 that gives the continuum neutron production.  it is 
*/ necessary to exclude this reaction from the reconstructed total
*/ cross section and to omit the associated energy-angle distribution.
*/ otherwise, the continuum neutron production will be counted twice.
*i acer.1941
               if (mt.eq.10) then
                  idone=0
                  call mess('unionx','redundant mt=10 found',
     &               'cross section and distribution excluded')
               endif
*i acer.2193
     &     (mf.eq.6.and.mt.eq.10).or.
*/ acer -- 03aug00
*/ this change can fix an infinite loop during acer plotting
*d acer.22571
      if (ep.lt.zero) then
*/ acer -- 03aug00
*/ for file 6 sections with only one subsection, the mt number is
*/ set to zero for the messages about energy-dependent yields.
*/ this is a trivial cosmetic patch and doesn't affect results.
*i acer.6414
         if (ikk.eq.nk) idone=1
*/ acer -- 03aug00
*/ this change is needed to handle nubar for jendl-3.2 u-235.  it
*/ was originally made at los alamos in april, but somehow didn't
*/ make it to the official update file.
*d acer.1090
         if (int.gt.2) nonlin=1
*ident up19
*/ heatr -- 03aug00
*/ as noted above, some jef, eff, and jeff evaluations use the
*/ redundant mt=10.  this value needs to be excluded from the
*/ heating and damage calculations.
*d heatr.639
                                 else if (mt.ne.10) then
*d heatr.690
                     else if (mt.ne.10) then
*i heatr.855
      if (mt.eq.10) go to 110
*ident up20
*/ acer -- 16aug00
*/ there is an error in processing angular distributions using ltt=3
*/ when newfor=0 (mcnp4b compatibility).  the extra tosend causes the
*/ code to skip over the first reaction after the elastic mf=4.
*/ this leads to a bad tyr=0 value in the ace file, which causes
*/ mcnp to issue a confusing error message about "sabcol," even
*/ when s(alpha,beta) is not being used.  this problem occurs only
*/ when processing the 150-mev evalutions from endf/b-vi.6.  it
*/ is probably best to use release 5 with mcnp4b anyway.  the
*/ release 5 and 6 data are identical below 20 mev.
*d acer.2292
*/ acer -- 16aug00
*/ the code is finding the wrong value for the lct parameter (lab
*/ or cm frame) when processing file 4 angular distributions if
*/ the section is fully isotropic.  this can result in an incorrect 
*/ value for the ace tyr parameter, which can result in an apparent 
*/ error from sabcol, even with no s(alpha,beta) data in the problem.
*/ this problem was introduced while the njoy97 coding was being
*/ converted to block structuring.
*d acer.5342
*i acer.5359
               lct=nint(a(iscr+3))
*d acer.5361
*/ acer -- 16aug00
*/ the consistency check for incorrect reference frame should take 
*/ place for isotropic distributions also.  sometimes, this check is 
*/ not a real error.  users should check the evaluation to see if the
*/ reference frame is really as intended by the evaluator.  as fixed,
*/ this check would have found the two problems above!
*d acer.17785
         if (na.ge.0) then
*/ acer-- 16aug00
*/ we are not currently handling law=5 for energy distributions.
*/ this occurs for u-233 fission from jef-2.2.  the evaluation can
*/ be patched by converting the lf=5 part of the distribution to
*/ lf=1, which is sampled much better by mcnp using cummulative
*/ probability distributions anyway.
*i acer.6098
            call error('acelf5','sorry. acer cannot handle lf=5.',
     &        'you will have to patch the evaluation to use lf=1.')
*/ acer -- 16aug00
*/ there are some additional places where skip6 should be skip6a.
*/ see up6 above.  this shows up for endf be-9 with newfor=1.
*d acer.8066
      external error,findex,skip6a,contio,listio,tab1io,moreio,tab2io
*d acer.8289
                     call skip6a(nin,0,0,a(iscr),law)
*d acer.8880
                        call skip6a(nin,0,0,a(iscr),law)
*d up17.42
                        call skip6a(nin,0,0,a(iscr),law)
*d acer.9004
                        call skip6a(nin,0,0,a(iscr),law)
*d acer.9027
                        call skip6a(nin,0,0,a(iscr),law)
*d acer.9725
                  call skip6a(nin,0,0,a(iscr),law)
*/ acer -- 16aug00
*/ there is an error in the law=7 part of up11.  this shows up
*/ when processing endf be-9 using newfor=0
*d up11.29
                        call tab2io(nin,0,0,a(iscr),nb,nw)
*/ acer -- 16aug00
*/ there is an error in the new skip6a routine introduced by up6
*/ that shows up when processing sections with law=7 with newfor=0.
*d up6.69
            nmu=nint(a(4))
*/ acer -- 16aug00
*/ missing initialization in ptlegc (this could affect incident 
*/ charged particles on some systems).
*i acer.2217
      dco=0
*ident up21
*/ groupr -- 28sep00
*/ the self-shielded cross sections are not being printed out
*/ correctly for the reactions, but total is ok.  the gendf
*/ file is ok, so libraries made with njoy99 are ok.
*d groupr.3613
                     call a10(ans(il,i,2),field(i))
*ident up22
*/ reconr -- 28sep00
*/ add capability to handle the new extension to the reich-moore
*/ resonance format that uses the sign of aj to designate which
*/ channel spin to use for a particular resonance.  based on
*/ coding provided by nancy larson, ornl.
*d reconr.2828,2942
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              ***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                 ***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                 ***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 up23
*/ gaminr -- 28sep00
*/  allow for up to 400 groups (added by request)
*d gaminr.78
      common/groupg/igg,ngg,egg(400)
*d gaminr.87
      dimension a(250000)  
*d gaminr.91
      dimension ng2s(400),ig2s(400) 
*d gaminr.455
      common/groupg/igg,ngg,egg(400)
*d gaminr.602
      data ngmax/400/
*d gaminr.521
      common/groupg/igg,ngg,egg(400)
*d gaminr.1138
      common/groupg/igg,ngg,egg(400)
*ident up24
*/ dtfr -- 28sep00
*/  allow for up to 400 groups (added by request)
*d dtfr.105,107
      common/dgrpn/egn(400),ngn
      common/dgrpg/egg(400),ngg
      common/dstore/a(20000),sig(200000) 
*d dtfr.110,111
      dimension spect(400)
      dimension fcap(400),ffis(400) 
*d dtfr.114
      data nwamax/20000/, nwsmax/200000/ 
*d dtfr.928
      common/dgrpn/egn(400),ngn
*d dtfr.932
      common/dstore/x(3500),y(3500),z(1000),a(212000)
*d dtfr.1262
      common/dgrpn/egn(400),ngn
*d dtfr.1409,1410
      common/dgrpn/egn(400),ngn
      common/dgrpg/egg(400),ngg
*d viewr.1294
      dimension lll(400)
*ident up25
*/ groupr -- 11oct00
*/  fix the section that reduces the number of sig figs in getdis.
*/  it was only acting on the in-group probabilities.  this helps
*/  to make the results for elastic and discrete inelastic matrices
*/  the same on different machines.  the basic idea is that these
*/  numbers are obtained by subtraction of numbers on the order of
*/  unity, so any results less than about 1e-7 are just random
*/  numbers and can be removed.
*d groupr.6637,6642
      ndig=7
      fact=ten**ndig
      do il=1,nl
         do ii=1,ng
            iii=nint(fact*ff(il,ii)+ten**(ndig-11))
            ff(il,ii)=iii/fact
         enddo
      enddo
*/ groupr -- 12oct00
*/  change the size of common groupg to agree with the changes
*/  made in gaminr above.
*d groupr.229
      common/groupg/igg,ngg,egg(400)
*d groupr.773
      common/groupg/igg,ngg,egg(400)
*d groupr.1919
      common/groupg/igg,ngg,egg(400)
*d groupr.3075
      common/groupg/igg,ngg,egg(400)
*d groupr.4275
      common/groupg/igg,ngg,egg(400)
*d groupr.7780
      common/groupg/igg,ngg,egg(400)
*ident up26
*/ acer -- 12oct00
*/  the current coding sometimes gets the threshold for charged
*/  particle production off by one point.
*i acer.8075
      data delt/1.d-10/
*i acer.8082
      data delt/1.e-10/
*d acer.8166
         do while (xss(esz+it-1).lt.thresh*(1-delt))
*ident up27
*/ dtfr -- 27oct00
*/  fix problem with finding right material and temperature
*/  on the pendf tape.  the goto loop was not translated correctly!        
*d dtfr.220
      idone=0
      do while (idone.eq.0)
*i dtfr.239
         else
            idone=1
*/ dtfr -- 27oct00
*/  fix error made in up24
*d up24.21
      common/dgrpg/egg(400),ngp
*ident up28
*/ acer -- 05nov00
*/  the pointer into the a array is not being correctly incremented
*/  for the "call moreio" line.  this only affects the new JEFF
*/  evaluation for beryllium, which has exceptionally detailed
*/  angulur tabulations.  found by fischer (fzk).
*i acer.2439
                              l=l+nw
*/  fix an indexing error in adjusting the normalization and
*/  precision for the pdf of angular distributions for law67 charged
*/  particle production that causes the pdf to be the same as the
*/  cdf.  this problem shows up for beryllium (n,2n) alpha production 
*/  in endf/b-vi, for example.  identified by konno (jaeri).
*d acer.8864
     &                          sigfig(renorm*xss(next+1+nx+ix),7,0)
*ident up29
*/ ccccr -- 05nov00
*/  the pointer in the e array for moreio is wrong.  the result of
*/  this is that larger group structures cannot be handled correctly
*/  for delayed neutrons.  found by broeders (karlsruhe).
*d ccccr.3140
         call moreio(nin,0,0,e(loc),nb,nw)
*ident up30
*/ heatr -- 05nov00
*/  the insert of the data value for nlmax was incorrectly done into
*/  the "sw" conditional block instead of after the conditional
*/  block was complete.  thus, it was only available to 32-bit
*/  versions of the code.  this was discovered by deleege (delft)
*/  when running in 64-bit mode on a vax/alpha.
*d up15.11
*i heatr.3298
      data nlmax/15/
*ident up31
*/ groupr -- 05nov00
*/  fix two problems with the ltt3 option for 150 mev evaluations.
*/  the incorrect index for the c array leads to findex problems
*/  caused by clobbering the index for the dynamic storage system.
*/  you also have to make sure that the "over" option that allows
*/  getfle to extrapolate to energies slightly higher than the
*/  upper limit of the table doesn't act at the break between
*/  the two energy ranges with ltt3.  this problem was reported 
*/  by wienke (sck-cen).
*d groupr.6838
      if (nne.eq.ne.and.e.lt.over*ehi) then
         if (ltt3.eq.3.and.lttn.eq.1) go to 210
         go to 300
      endif
*d groupr.6850,6851
         call tab2io(nin,0,0,c(ifls),nb,nwc)
         ne=nint(c(ifls+5))
*ident up32
*/ reconr -- 05nov00
*/  some fission products from the jendl-3.2 library include
*/  an unresolved resonance range with no corresponding resolved
*/  range.  trkov (iaea) proposed the following fix.
*i reconr.672
      if (eresr.lt.eresl) eresr=eresl
*ident up33
*/ gaminr -- 18jan01
*/  the photoatomic group cross sections are not printed out
*/  correctly for a p-order greater than 5.
*d gaminr.1075
     &              write(nsyso,'(13x,1p,6a11)') (field(i),i=7,nl)
*ident up34
*/ groupr -- 29jan01
*/  need more storage in groupr to handle mt=91 for am243 from
*/  endf/b-vi release 5, which goes to 30 mev.  the symptom was
*/  "storage exceeded" from cm2lab.
*d groupr.248
      dimension a(150000)
*d groupr.273
      iamax=150000
*ident up35
*/ groupr -- 08feb01
*/  when we increased the common block for photon group structures
*/  to allow as many as 400 groups (see up25), we forgot to update
*/  the parameter ngmax.  this causes a "too many groups" error if
*/  you run with more than 150 gamma groups.
*d groupr.2000
      data ngmax/400/
*ident up36
*/ acer -- 08feb01
*/  in up17, we checked for isotropic distributions in order to use
*/  a more compact presentation.  the logic misses one special case,
*/  namely, pb208 from endf/b-vi release 6.
*d up17.41
                     if (izap.ne.ip.or.law.ne.1) then
*d up17.45
                     else
*ident up37
*/ reconr -- 09feb01
*/  all through njoy, we have been using 1e10 ev as our idea of
*/  an infinite energy.  progress happens, and red cullen at llnl
*/  is putting out an endf version of the evaluated photon data 
*/  library (epdl97), which contains data to 100 gev.  the following
*/  change prevents reconr from going into an infinite loop in the
*/  emerge routine with 100 gev data. 
*d reconr.4126
      data finity/.99d12/
*d reconr.4130
      data finity/.99e12/
*ident up38
*/ njoy -- 09feb01
*/  keep on increasing infinity for the 100 gev data.  the routines
*/  gety1, gety2, and terpa return an "infinite" energy at the end
*/  of the table, and we now increase that to 1e12 ev.  this doesn't
*/  seem to cause any problems in njoy modules (such as groupr) that 
*/  still check for return values of 1e10 or more; all the standard
*/  test problems still work fine.
*d njoy.2204
      data xbig/1.d12/
*d njoy.2208
      data xbig/1.e12/
*d njoy.2418
      data xbig/1.d12/
*d njoy.2422
      data xbig/1.e12/
*d njoy.2532
      data xbig/1.d12/
*d njoy.2536
      data xbig/1.e12/
*ident up39
*/ gaminr -- 09feb01
*/  keep on increasing infinity for the 100 gev data.
*d gaminr.106
      data emax/1.d12/
*d gaminr.110
      data emax/1.e12/
*d gaminr.779
      data emax/1.d12/
*d gaminr.782
      data emax/1.e12/
*d gaminr.1164
      data emax/1.d12/
*d gaminr.1186
      data emax/1.e12/
*ident up40
*/ acer -- 23mar01
*/  due to a bad if clause, the contribution to heating from charged
*/  particles is not being included for mf=6, law 3 or 4.  this was
*/  noticed in the run for endf/b-vi be-9 by lanl/x-5.  the errors in
*/  this particular case are quite small because of the small cross
*/  sections for charged-particle emission.  this error will only
*/  effect mcnpx calculations for coupled neutron-proton transport.
*d acer.9220,9231
c                       add in contribution to heating
                        naa=nint(xss(hpd+1))
                        do ie=it,nes
                           e=xss(esz+ie-1)
                           ss=0
                           if (ie.ge.iaa) ss=xss(2+k+ie-iaa)
                           tt=xss(next+1)*(e-xss(next))*ss
                           xss(hpd+2+naa+ie-it)=xss(hpd+2+naa+ie-it)
     &                       +tt
                        enddo
*ident up41
*/ acer -- 27mar01
*/  the value "nr = 0", implying linear interpolation over all points,
*/  is not printed on the acer output listing for two cases, as reported
*/  by lanl/x-5.  these errors do not affect mcnp results, but the
*/  repair makes the printout for photon yields and energy distributions
*/  match those for other types of data.
*i acer.10808
               write(nsyso,'(12x,''nr ='',i4)') m
*d acer.10810
*i acer.10998
               write(nsyso,'(12x,''nr ='',i4)') m
*d acer.11000,11002
*ident up42
*/ purr -- 27mar01
*/  remove the timers that are given as each ladder is processed
*/  in order to reduce the number of diffs that show up when
*/  succesive runs are checked for qa purposes using the same
*/  sequence of random numbers.  for lanl/x-5.
*d purr.1746
      external fsort,ladr2,fsrch
*d purr.1798
     &  ''capture'')') e,spot,dbart,sigx
*d purr.2145,2147
      if (iprint.gt.0) write(nsyso,'(i6,1p,4e12.4)')
     &  iladr,totf,elsf,fisf,capf
*ident up43
*/ heatr -- 27mar01
*/  the roundup applied to the first energy grid point should 
*/  be smaller now that we are routinely working with 7-digit 
*/  energies.  the effect if this in current files is that the 
*/  first energy in any of the heating and damage reactions is 
*/  a little larger than the normal 1e-5.  this shows up as a 
*/  zero heating or damage value for the first point in the mcnp 
*/  ace files, which is strange looking, but of little significant 
*/  impact on real calculations.  reported by lanl/x-5.
*d heatr.425
      data rup/1.0000001d0/
*ident up44
*/ acer -- 29mar01
*/  lanl/x-5 has requested that the main container array be increased
*/  in size to allow bigger ace files to be generated. it is also
*/  necessary to increase the i7 length field on the xsdir cards to i8
*/  to accomodate the larger ace files.
*d acer.257
      max3=1500000
*d acer.4662
      common/xsst/xss(1500000),n3
*d acer.5601
      common/xsst/xss(1500000),n3
*d acer.5762
      common/xsst/xss(1500000),n3
*d acer.5951
      common/xsst/xss(1500000),n3
*d acer.6322
      common/xsst/xss(1500000),n3
*d acer.7383
      common/xsst/xss(1500000),n3
*d acer.8055
      common/xsst/xss(1500000),n3
*d acer.9754
      common/xsst/xss(1500000),n3
*d acer.10202
      common/xsst/xss(1500000),n3
*d acer.10675
      common/xsst/xss(1500000),n3
*d acer.11068
      common/xsst/xss(1500000),n3
*d acer.11588
      common/xsst/xss(1500000),n3
*d acer.11649
     &        '(a10,f12.6,'' filename route'',i2,i4,i8,2i6,1p,e10.3,
*d acer.11653
     &        '(a10,f12.6,'' filename route'',i2,i4,i8,2i6,1p,e10.3)')
*d acer.11659
     &        '(a13,f12.6,'' file route'',i2,i4,i8,2i6,1p,e10.3,
*d acer.11663
     &        '(a13,f12.6,'' file route'',i2,i4,i8,2i6,1p,e10.3)')
*d acer.11689
      common/xsst/xss(1500000),n3
*d acer.12854
      common/xsst/xss(1500000),n3
*d acer.13452
      common/xsst/xss(1500000),n3
*d acer.13591
      common/xsst/xss(1500000),n3
*d acer.13771
      common/xsst/xss(1500000),n3
*d acer.14170
      common/xsst/xss(1500000),n3
*d acer.14274
     &     '(a10,f12.6,'' filename route'',i2,2h 1,i8,2i6,1p,e10.3)')
*d acer.14278
     &     '(a13,f12.6,'' filename route'',i2,2h 1,i8,2i6,1p,e10.3)')
*d acer.14305
      common/xsst/xss(1500000),n3
*d acer.14462
      common/xsst/xss(1500000),n3
*d acer.14548
      common/xsst/xss(1500000),n3
*d acer.14640
     &     '(a10,f12.6,'' filename route'',i2,2h 1,i8,2i6,1p,e10.3)')
*d acer.14644
     &     '(a13,f12.6,'' filename route'',i2,2h 1,i8,2i6,1p,e10.3)')
*d acer.14674
      common/xsst/xss(1500000),n3
*d acer.15012
      common/xsst/xss(1500000),n3
*d acer.15107
      common/xsst/xss(1500000),n3
*d acer.15187
     &     '(a10,f12.6,'' filename route'',i2,'' 1'',i8,2i6,1p,e10.3)')
*d acer.15191
     &     '(a13,f12.6,'' filename route'',i2,'' 1'',i8,2i6,1p,e10.3)')
*d acer.15216
      common/xsst/xss(1500000),n3
*d acer.16604
      common/xsst/xss(1500000),n3
*d acer.17057
      common/xsst/xss(1500000),n3
*d acer.17436
     &     '(a10,f12.6,'' filename route'',i2,'' 1'',i8,2i6,1p,e10.3)')
*d acer.17440
     &     '(a13,f12.6,'' filename route'',i2,'' 1'',i8,2i6,1p,e10.3)')
*d acer.17459
      common/xsst/xss(1500000),n3
*d acer.17727
      common/xsst/xss(1500000),n3
*d acer.18534
      common/xsst/xss(1500000),n3
*d acer.19545
      common/xsst/xss(1500000),n3
*d acer.19817
      common/xsst/xss(1500000),n3
*d acer.19934
      common/xsst/xss(1500000),n3
*d acer.20164
      common/xsst/xss(1500000),n3
*d acer.20610
      common/xsst/xss(1500000),n3
*d acer.21222
      common/xsst/xss(1500000),n3
*d acer.21815
      common/xsst/xss(1500000),n3
*ident up45
*/ acer -- 08apr01
*/  as discovered by jean christophe sublet, sun forte6 f95 is 
*/  finiky about opening a scratch file that is already open,
*/  although all other compilers used for njoy thus far were more
*/  accepting.  we just have to be careful to close a unit used
*/  as a scratch file before reusing the unit for another purpose.
*i acer.2082
      call closz(nscr)
*ident up46
*/ gaspr -- 09apr01
*/  close another scratch unit.
*i gaspr.838
      call closz(nscr1)
*ident up47
*/ acer -- 09apr01
*/  the length published for thermal data files is too long by one
*/  for cases including incoherent elastic scattering.  for endf,
*/  this is poly, h(zrh), and cold solid methane.  discovered by
*/  roberto orsi (enea-bologna).
*d acer.13517
*ident up48
*/ acer -- 09apr01
*/  the landh parameter should be zero (not -1) for isotropic 
*/  subsections of mf=6 described using law=3.  this occurs for 
*/  the reactions (n,p0) through (n,a0) in be-9 from endf/b-vi.  
*/  Noted by bob little (lanl/x-5).
*i acer.9209
                        if (law.eq.3) xss(landh+jp-1)=0
*ident up49
*/ acer -- 09apr01
*/  the representation for ace law3/33 should use -q instead of 
*/  abs(q) in order to handle two-body reactions for isomeric
*/  targets.  this change in the ace specifications was recommended
*/  by bob little (lanl/x-5) after a query by waclaw gudowski.  it
*/  only affects a small number of evaluations.
*d acer.5465
               xss(next+9)=sigfig(x*(-q),7,0)
*d acer.8964
                  xss(next)=sigfig((1+amass)*(-q)/amass,7,0)
*d acer.9167
                           xss(next)=sigfig((1+amass)*(-q)/amass,7,0)
*d acer.9216
                        xss(next)=sigfig((1+amass)*(-q)/amass,7,0)
*d acer.16168
               xss(nex)=sigfig(-q,7,0)
*d acer.16526
                           xss(nex)=sigfig(-q,7,0)
*ident up50
*/ acer -- 09apr01
*/  most of the jendl photonuclear evaluations currently available
*/  from http://iaeand.iaea.or.at/photonuclear/ crash with an i/o
*/  error because they use a non-conforming format where mf=6,
*/  mt=201-27 are used to represent particle production.  we are
*/  providing a clearer error message for the user's convenience.
*/  these evaluations cannot be used in njoy or mcnpx in their
*/  current form.
*i acer.15310
            if (mfd.eq.6.and.mtd.ge.201.and.mtd.le.207)
     &        call error('acephn','mf=6/mt=201-207 not supported.',
     &        'does not conform to endf format.')
*ident up51
*/ acer -- 12apr01
*/  add a capability to handle a two-body recoil subsection of mf=6
*/  for photonuclear files.  this may be useful for representing the
*/  photodisintegration of the deuteron with full distributions for
*/  both neutron and proton.  we tested the patch using a modified
*/  version of the g+2H evaluation from JENDL.
*d acer.16011
c
c                    ***special steps for two-body recoil
c                    ***back up to the corresponding law=2 distr.
                     izarec=0
                     awprec=0
                     if (izap.eq.ip.and.law.eq.4) then
                        izarec=izap
                        awprec=awp
                        mf=6
                        call findf(matd,mf,mt,nin)
                        call contio(nin,0,0,a(iscr),nb,nw)
                        call tab1io(nin,0,0,a(iscr),nb,nw)
                        izap=nint(c1h)
                        awp=c2h
                        law=l2h
                        jscr=iscr+nw
                        do while (nb.ne.0)
                           call moreio(nin,0,0,a(jscr),nb,nw)
                           jscr=jscr+nw
                        enddo
                     endif
c     
c                    ***law2 angular distribution
c                    ***also used for law 4 two-body recoils
                     if ((izap.eq.ip.and.law.eq.2).or.
     &                 (izarec.eq.ip.and.law.eq.2)) then
                        lld=jscr
*i acer.16030
                        if (izarec.eq.0) then
                           awpp=awp
                        else
                           awpp=awprec
                        endif
*i acer.16036
                              if (izarec.ne.0) then
                                 nl=nint(a(lld+5))
                                 do iil=1,nl
                                    if (mod(iil,2).eq.1) then
                                       a(lld+5+iil)=-a(lld+5+iil)
                                    endif
                                 enddo
                              endif
*d acer.16104
                           a(llht+7+2*iie)=(awr-awpp)*(e+q)/awr
*d acer.16354
                        if (law.ne.1.and.law.ne.2.and.law.ne.4) then
*d acer.16356
     &                       'law=2, or law=4 currently')
*i acer.16540
                        else if (law.eq.4) then
                           xss(last+1)=33
                           xss(nex)=0
                           xss(nex+1)=2
                           nnr=nint(a(iscr+4))
                           nnp=nint(a(iscr+5))
                           xss(nex+2)=sigfig(a(iscr+6+2*nnr)/emev,7,0)
                           xss(nex+3)= 
     &                        sigfig(a(iscr+4+2*nnr+2*nnp)/emev,7,0)
                           xss(nex+4)=1
                           xss(nex+5)=1
                           nex=nex+2+2*2
                           xss(last+2)=nex-dlwp+1
                           xss(nex)=sigfig(-q,7,0)
                           xss(nex+1)=sigfig(awr/(1+awr),7,0)
                           nex=nex+2
*ident up52
*/ acer -- 13apr01
*/  add a capability to handle tabulated sections of File 5 (lf=1)
*/  for photonuclear files.  Such sections are used in the Russian
*/  evaluations for three isotopes of plutonium included in the
*/  iaea photonuclear compilation.  this also fixes a bug in the
*/  storage of fission nubar.  the first point for energy distributions
*/  often has a nonrealistic sharp triangle given for the spectrum.
*/  this can cause problems with the vertical scale for plots
*/  because the emission probabilities get very large for small
*/  ranges of secondary energy.  therefore, we ignore the first
*/  incident energy in determining the vertical scale for the plot.
*d acer.15858,15859
                        xss(nex+3+j)=sigfig(fnubar(5+2*nr+2*j)/emev,7,0)
                        xss(nex+3+ne+j)=sigfig(fnubar(6+2*nr+2*j),7,0)
*d acer.16247
                  if (lf.eq.1) then
                     call tab2io(nin,0,0,a(iscr),nb,nw)
                     m=nint(a(iscr+4))
                     n=nint(a(iscr+5))
                     jnt=nint(a(iscr+7))
                     jnt=mod(jnt,10)
                     if (jnt.gt.2) jnt=2
                     if (m.ne.1.or.jnt.ne.2) then
                        xss(nex)=m
                        do j=1,m
                           xss(j+nex)=a(2*j+4+iscr)
                           jnt=nint(a(2*j+5+iscr))
                           jnt=mod(jnt,10)
                           if (jnt.gt.2) jnt=2
                           xss(j+m+nex)=jnt
                        enddo
                        nex=nex+1+2*m
                     else
                        xss(nex)=0
                        nex=nex+1
                     endif
                     xss(nex)=n
                     nexn=nex+n
                     nexd=nexn+n+1
                     ne=n
                     do j=1,ne
                        call tab1io(nin,0,0,a(iscr),nb,nw)
                        jscr=iscr
                        do while (nb.ne.0)
                           jscr=jscr+nw
                           call moreio(nin,0,0,a(jscr),nb,nw)
                        enddo
                        e=c2h
                        xss(nex+j)=sigfig(e/emev,6,0)
                        xss(nexn+j)=nexd-dlwp+1
                        m=n1h
                        n=n2h
                        jnt=nint(a(iscr+5+2*m))
                        xss(nexd)=jnt
                        xss(nexd+1)=n
                        nexd=nexd+1
                        xss(nexd+1+2*n)=0
                        do ki=1,n
                           ep=a(iscr+4+2*m+2*ki)
                           ll=iscr+4+2*m+2*ki
                           xss(ki+nexd)=sigfig(a(ll)/emev,7,0)
                           xss(ki+n+nexd)=sigfig(a(ll+1)*emev,7,0)
                           if (xss(ki+n+nexd).lt.rmin) xss(ki+n+nexd)=0
                           if (ki.gt.1.and.jnt.eq.1) xss(ki+2*n+nexd)=
     &                       xss(ki+2*n-1+nexd)+a(ll-1)*(a(ll)-a(ll-2))
                           if (ki.gt.1.and.jnt.eq.2) xss(ki+2*n+nexd)=
     &                       xss(ki+2*n-1+nexd)+((a(ll-1)
     &                       +a(ll+1))/2)*(a(ll)-a(ll-2))
                        enddo
c                       renormalize
                        renorm=1
                        if (xss(3*n+nexd).ne.zero) 
     &                    renorm=1/xss(3*n+nexd)
                        do ki=1,n
                           xss(ki+n+nexd)=
     &                       sigfig(xss(ki+n+nexd)*renorm,7,0)
                           xss(ki+2*n+nexd)=
     &                       sigfig(xss(ki+2*n+nexd)*renorm,9,0)
                        enddo
                        nexd=nexd+3*n+1
                     enddo
                     nex=nexd
                  else if (lf.eq.7.or.lf.eq.9) then
*d acer.22424
               do ie=2,ne
*ident up53
*/ groupr -- 11jun01
*/  if the file6 distribution is fully isotropic (law=3), the getfle 
*/  routine doesn't realize that when doing a discrete recoil (law=4).  
*/  we create a special flag of law=-4 to pass the fact of isotropy 
*/  into getfle.  this problem only affects runs that compute a 
*/  transfer matrix for the recoil particle when the first particle
*/  emitted is given as totally isotropic (for example, mt=701 for 
*/  endf be-9).  the error message is "desired energy above highest 
*/  given."  found by dieter leichtle (fzk).
*i groupr.4859
         lf=nint(c(l+3))
*i groupr.4860
         if (lf.eq.3) law=-4
*i groupr.4869
      if (law.eq.-4) go to 194
*i groupr.6786
         if (law.eq.-4) iso=1
*ident up54
*/ reconr -- 12jun01
*/  allow for the series of mt numbers 875-891 that can be used
*/  to represent different levels of the (n,2n) reaction in the
*/  same way that 600-649 are used to represent different levels
*/  of the (n,p) reaction.  the code expects that mf=3/mt=16
*/  contains the sum of mt=875 through 891 in the same way that
*/  mt=103 contains the sum of 600-649.  this representation is 
*/  used for be-9 for eff-3.1 and jeff-3.0.
*d reconr.1696
         if (mth.ge.900) go to 150
*ident up55
*/ heatr -- 12jun01
*/  if mt=875-891 appears in the file, mt=16 is redundant.  this
*/  is analogous to the way mt=107 is redundant if mt=800-850 
*/  is present.
*d heatr.412
      common/heat4/mt103,mt104,mt105,mt106,mt107,mt16
*i heatr.440
      mt16=0
*i heatr.499
                        if (mtd.ge.875.and.mtd.lt.891) mt16=1
*d heatr.783
      common/heat4/mt103,mt104,mt105,mt106,mt107,mt16
*i heatr.865
      if (mt.eq.16.and.mt16.gt.0) go to 110
*/  the integration over secondary energy for law 7 in getsix
*/  must allow for histogram interpolation as used in be-9
*/  from eff-3.1.  the effect of this is to get especially
*/  bad particle energies for the discrete neutron in mt=876.
*i heatr.3008
         iint=nint(c(l+7))
*d heatr.3020
            if (i.gt.1) then
               if (iint.eq.1) then
                  h=h+(xx-xl)*el
               else
                  h=h+(xx-xl)*(en+el)/2
               endif
            endif
*d heatr.3022
            if (i.gt.1) then
               if (iint.eq.1) then
                  d=d+(xx-xl)*fl
               else
                  d=d+(xx-xl)*(fn+fl)/2
               endif
            endif
*ident up56
*/ acer 12jun01
*/  the angle-energy law in file 6 is always causing trouble.
*/  it is especially difficult when more than one subsection
*/  is used to describe the emission for a particle, because
*/  an overall angular distribution for the reaction must be
*/  contructed.  with the new formats, it is easy to eliminate
*/  law=7 sections by converting them to law=1 with tabulated
*/  angular distributions.
*d acer.2111,2112
c     mf4 and 5.  mf6 is also copied, unless law=7 is found,
c     in which case the law=7 data are converted to law=1.
c     all other conversions of the distributions will be
c     done in acelod.
*i acer.2132
      zero=0
      one=1
*i acer.2365
               if (newfor.eq.1.and.lf.eq.7) a(iscr+3)=1
*i acer.2371
               else if (lf.eq.7.and.newfor.eq.1) then
c                 law=7 for newfor=1 -- convert the law7
c                 data into law1 format.
                  call tab2io(nin,0,0,b,nb,nw)
                  ne=nint(b(6))
                  do ie=1,ne
c                    read in the data
                     call tab2io(nin,0,0,a(iscr),nb,nw)
                     ei=a(iscr+1)
                     intmu=nint(a(iscr+7))
                     nmu=n2h
                     loc=iscr+nmu
                     do imu=1,nmu
                        a(iscr+imu-1)=loc
                        call tab1io(nin,0,0,a(loc),nb,nw)
                        intep=nint(a(loc+7))
                        loc=loc+nw
                        do while (nb.ne.0)
                           call moreio(nin,0,0,a(loc),nb,nw)
                           loc=loc+nw
                        enddo
                     enddo
c                    fix up the tab2 for law1
                     if (ie.eq.1) then
                        b(3)=10+intmu
                        b(4)=intep
                        call tab2io(0,nout,0,b,nb,nw)
                        ncs(nxc)=ncs(nxc)+2 
                     endif
c                    construct a union grid for eprime
                     igrd=loc
                     ngrd=0
                     do imu=1,nmu
                        loc=nint(a(iscr+imu-1))
                        m=nint(a(loc+4))
                        n=nint(a(loc+5))
                        do iep=1,n
                           ngrd=ngrd+1
                           a(igrd+ngrd-1)=a(loc+4+2*m+2*iep)
                        enddo
                     enddo
                     call ordr(a(igrd),ngrd)
c                    interpolate for angular distributions
c                    on the union eprime grid to construct
c                    the law1 distribution.
                     ians=igrd+ngrd
                     a(ians)=0
                     a(ians+1)=ei
                     a(ians+2)=0
                     a(ians+3)=2*nmu
                     a(ians+4)=ngrd*(2+2*nmu)
                     a(ians+5)=ngrd
                     ll=ians+6
                     do iep=1,ngrd
                        ep=a(igrd+iep-1)
                        a(ll)=ep
                        ss=0
                        do imu=1,nmu
                           loc=nint(a(iscr+imu-1))
                           ipp=2
                           irr=1
                           call terpa(ff,ep,epn,idis,a(loc),
     &                       ipp,irr)
                           a(ll+2*imu)=a(loc+1)
                           a(ll+1+2*imu)=ff
                           if (imu.gt.1) then
                              dmu=a(ll+2*imu)-a(ll+2*imu-2)
                              if (intmu.eq.1) then
                                 ss=ss+dmu*a(ll+1+2*imu-2)
                              else
                                 ss=ss+dmu*
     &                             (a(ll+1+2*imu)+a(ll+1+2*imu-2))/2
                              endif
                           endif
                        enddo
                        a(ll+1)=ss
                        do imu=1,nmu
                           if (ss.ne.zero) then
                              a(ll+1+2*imu)=a(ll+1+2*imu)/ss
                           else
                              a(ll+1+2*imu)=one/2
                           endif
                        enddo
                        ll=ll+2+2*nmu
                     enddo
                     call listio(0,nout,0,a(ians),nb,nw)
                     ll=ians+nw
                     do while (nb.ne.0)
                        call moreio(0,nout,0,a(ll),nb,nw)
                        ll=ll+nw
                     enddo
                     nw=ngrd*(2+2*nmu)
                     nw=(nw+5)/6
                     ncs(nxc)=ncs(nxc)+1+nw 
                  enddo
*i acer.2457
c
      subroutine ordr(x,n)
c     ******************************************************************
c     sort the n elements of x into ascending order
c     removing any duplicate elements
c     ******************************************************************
*if sw
      implicit real*8 (a-h,o-z)
*endif
      dimension x(*)
*if sw
      data small/1.d-10/
*else
      data small/1.e-10/
*endif
c
      if (n.le.2) return
c     sort
      i=0
  110 i=i+1
      j=i
  120 j=j+1
      if (x(j).lt.x(i)) then
         tsave=x(j)
         x(j)=x(i)
         x(i)=tsave
      endif
      if (j.lt.n) go to 120
      if (i.lt.n-1) go to 110
c     remove duplicates
      m=n
      i=1
      do while (i.lt.m)
         i=i+1
         if (abs(x(i)-x(i-1)).le.small*x(i)) then
            m=m-1
            do k=i,m
               x(k)=x(k+1)
            enddo
            i=i-1
         endif
      enddo
      n=m
      return
      end
*d acer.6792
                  else if (lang.gt.2.and.newfor.eq.0) then 
*i acer.6861
c
c                 ***convert tabulated distribution to law 61
                  else if (lang.gt.2.and.newfor.eq.1) then 
                     xss(ki+3*n+nexd)=nexcd-dlw+1
                     ll=iscr+6+ncyc*(ki-1)
                     intmu=lang-10
                     xss(nexcd)=intmu
                     nmu=na/2
                     xss(nexcd+1)=nmu
                     do imu=1,nmu
                        xss(nexcd+1+imu)=a(ll+2*imu)
                        xss(nexcd+1+nmu+imu)=a(ll+2*imu+1)
                        if (imu.eq.1) then
                           sum=0
                           xss(nexcd+1+2*nmu+imu)=0
                        else
                           del=a(ll+2*imu)-a(ll+2*imu-2)
                           if (intmu.eq.1) then
                              sum=sum+del*a(ll+1+2*imu-2)
                              xss(nexcd+1+2*nmu+imu)=sum
                           else
                              av=(a(ll+1+2*imu)+a(ll+1+2*imu-2))/2
                              sum=sum+del*av
                              xss(nexcd+1+2*nmu+imu)=sum
                           endif
                        endif
                     enddo
                     do imu=1,nmu
                        xss(nexcd+1+imu)=
     &                    sigfig(xss(nexcd+1+imu)/sum,7,0)
                        xss(nexcd+1+nmu+imu)=
     &                    sigfig(xss(nexcd+1+nmu+imu)/sum,7,0)
                        xss(nexcd+1+2*nmu+imu)=
     &                    sigfig(xss(nexcd+1+2*nmu+imu)/sum,7,0)
                     enddo
                     nexcd=nexcd+2+3*nmu
*/  add support for law1 with tabulated angular distributions
*/  to the charged-particle section.
*i up9.45
                        if (law.eq.1.and.lang.ge.11) lawnow=61
*d up9.55
c                                legendre or tabulated distribution
*d up9.58,up9.97
                                    if (lang.eq.1) then
                                       a(ll)=0
                                       a(ll+1)=ep
                                       a(ll+2)=0
                                       a(ll+3)=0
                                       a(ll+4)=na
                                       a(ll+5)=0
                                       do ia=1,na
                                          lll=lld+7+ncyc*(ig-1)
                                          a(ll+5+ia)=0
                                          if (a(lll).ne.zero) then
                                             a(ll+5+ia)=a(lll+ia)/a(lll)
                                          endif
                                       enddo
                                       call ptleg2(a(ll))
                                       intmu=2
                                       nmu=nint(a(ll+5))
                                       llx=ll+6
                                    else
                                       intmu=lang-10
                                       nmu=na/2
                                       llx=lld+6
                                    endif
                                    xss(next+1+3*ng+ig)=nexcd-dlwh+1
                                    xss(nexcd)=intmu
                                    xss(nexcd+1)=nmu
                                    do imu=1,nmu
                                       xss(nexcd+1+imu)=
     &                                   a(llx+2*imu)
                                       xss(nexcd+1+nmu+imu)=
     &                                   a(llx+1+2*imu)
                                       if (imu.eq.1) then
                                          sum=0
                                          xss(nexcd+1+2*nmu+imu)=0
                                       else
                                          del=a(llx+2*imu)
     &                                      -a(llx+2*imu-2)
                                          if (intmu.eq.1) then
                                             sum=sum
     &                                         +del*a(llx+1+2*imu-2)
                                             xss(nexcd+1+2*nmu+imu)=sum
                                          else
                                             av=(a(llx+1+2*imu)
     &                                         +a(llx+1+2*imu-2))/2
                                             sum=sum+del*av
                                             xss(nexcd+1+2*nmu+imu)=sum
                                          endif
                                       endif
                                    enddo
                                    do imu=1,nmu
                                       xss(nexcd+1+imu)=sigfig(
     &                                   xss(nexcd+1+imu)/sum,7,0)
                                       xss(nexcd+1+nmu+imu)=sigfig(
     &                                   xss(nexcd+1+nmu+imu)/sum,7,0)
                                       xss(nexcd+1+2*nmu+imu)=sigfig(
     &                                   xss(nexcd+1+2*nmu+imu)/sum,7,0)
                                    enddo
*ident up57
*/ acer -- 12jun01
*/  changes to acer needed to support the eff-3.1/jeff-3.0
*/  representation for be-9.  we have to allow for the series
*/  of mt numbers 875-891.  if present, mt=16 is redundant and
*/  must appear after the distributions in the reaction list.
*/  the section for mt=876 has two subsections for neutron
*/  emission.  this problem is handled by the previous update.
*i acer.551
      common/ace9/mt16
*i acer.642
      mt16=0
*i acer.663
               if (mtd.ge.875.and.mtd.le.891) mt16=1
*d acer.1939
     &           (iverf.ge.6.and.mt.gt.900)) then
*i acer.4671
      common/ace9/mt16
*i acer.4760
                     if (mt16.gt.0.and.mt.eq.16) nr=nr-1 
*d acer.4759
                  if ((mt.ge.5.and.mt.le.91).or. 
     &                (mt.ge.875.and.mt.le.899)) then
*i acer.5021
            if (mt.gt.91.and.mt.le.849) iskip=1
            if (mt16.gt.0.and.mt.eq.16) iskip=1
*d acer.5118,5120
      call findf(matd,3,2,nin)
*d acer.5127
            if (mt.gt.91.and.mt.le.849) iskip=0
            if (mt16.gt.0.and.mt.eq.16) iskip=0
*d acer.9481
                              renorm=1
                              if (xss(next+3*npep).ne.zero)
     &                          renorm=1/xss(next+3*npep)
*/  add the new reaction names to the mtname routine
*d acer.11389
      character*10 hndf(457)
*d acer.11393,11394
      character*10 hndf9(50)
      character*10 hndf10(7)
      character*10 hndf11(1)
*d acer.11403,11404
      equivalence (hndf9(1),hndf(400))
      equivalence (hndf10(1),hndf(450))
      equivalence (hndf11(1),hndf(457))
*i acer.11514
      data hndf9/'(n,x)     ',
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     ',
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     ',
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     ',
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     ',
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     ',
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     ', 
     &  '(n,2n*0)  ',
     &  '(n,2n*1)  ', '(n,2n*2)  ', '(n,2n*3)  ', '(n,2n*4)  ',
     &  '(n,2n*5)  ', '(n,2n*6)  ', '(n,2n*7)  ', '(n,2n*8)  ',
     &  '(n,2n*9)  ', '(n,2n*10) ', '(n,2n*11) ', '(n,2n*12) ',
     &  '(n,2n*13) ', '(n,2n*14) ', '(n,2n*15) ', '(n,2n*c)  ',
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     ', 
     &  '(n,x)     ', '(n,x)     ', '(n,x)     ', '(n,x)     '/ 
*d acer.11520,11522
      data hndf10/'(n,xn)    ','(n,xgma)  ','(n,xp)    ',
     &  '(n,xd)    ','(n,xt)    ','(n,xhe3)  ','(n,xa)    '/
      data hndf11(1)/'damage    '/
*d acer.11526,11528
         if (i.ge.201.and.i.le.207) i=i+249
         if (i.ge.600) i=i-450
         if (mt.eq.444) i=457
*d acer.11534
            name=hndf(mt+250)
*d acer.11536
            name=hndf(457)
*/  need more storage for eff-3.1 be-9
*d up44.8
      max3=3000000
*d up44.10
      common/xsst/xss(3000000),n3
*d up44.12
      common/xsst/xss(3000000),n3
*d up44.14
      common/xsst/xss(3000000),n3
*d up44.16
      common/xsst/xss(3000000),n3
*d up44.18
      common/xsst/xss(3000000),n3
*d up44.20
      common/xsst/xss(3000000),n3
*d up44.22
      common/xsst/xss(3000000),n3
*d up44.24
      common/xsst/xss(3000000),n3
*d up44.26
      common/xsst/xss(3000000),n3
*d up44.28
      common/xsst/xss(3000000),n3
*d up44.30
      common/xsst/xss(3000000),n3
*d up44.32
      common/xsst/xss(3000000),n3
*d up44.42
      common/xsst/xss(3000000),n3
*d up44.44
      common/xsst/xss(3000000),n3
*d up44.46
      common/xsst/xss(3000000),n3
*d up44.48
      common/xsst/xss(3000000),n3
*d up44.50
      common/xsst/xss(3000000),n3
*d up44.52
      common/xsst/xss(3000000),n3
*d up44.58
      common/xsst/xss(3000000),n3
*d up44.60
      common/xsst/xss(3000000),n3
*d up44.62
      common/xsst/xss(3000000),n3
*d up44.68
      common/xsst/xss(3000000),n3
*d up44.70
      common/xsst/xss(3000000),n3
*d up44.72
      common/xsst/xss(3000000),n3
*d up44.78
      common/xsst/xss(3000000),n3
*d up44.80
      common/xsst/xss(3000000),n3
*d up44.82
      common/xsst/xss(3000000),n3
*d up44.88
      common/xsst/xss(3000000),n3
*d up44.90
      common/xsst/xss(3000000),n3
*d up44.92
      common/xsst/xss(3000000),n3
*d up44.94
      common/xsst/xss(3000000),n3
*d up44.96
      common/xsst/xss(3000000),n3
*d up44.98
      common/xsst/xss(3000000),n3
*d up44.100
      common/xsst/xss(3000000),n3
*d up44.102
      common/xsst/xss(3000000),n3
*d up44.104
      common/xsst/xss(3000000),n3
*d up44.106
      common/xsst/xss(3000000),n3
*ident up58
*/ broadr -- 10jul01
*/  increase the storage area in broadr to reduce paging and make
*/  comparisons between njoy99 and njoy2001 easier.  there will
*/  normally be a small difference in the grids produced by broadr
*/  each time paging takes place, and this makes it hard to compare
*/  files using diff.
*d broadr.113
      dimension a(95000)
*d broadr.137
      namax=95000
*ident up59
*/ acer -- 20jul01
*/  in charged-particle emission, the first point for energy 
*/  distributions often has a nonrealistic sharp triangle given for 
*/  the spectrum.  this can cause problems with the vertical scale 
*/  incident energy in determining the vertical scale for the plot.
*d acer.21529
               do ie=2,ne
*ident up60
*/ reconr -- 24sep01
*/  occasionally, reactions are given with a nonzero cross section
*/  at threshold, even though this violates endf procedures.  reconr
*/  had some logic for handling this that was being overwritten by
*/  another change.  we fix it here by inserting an extra energy
*/  point just above the threshold and zeroing the cross section at
*/  the threshold.  a diagnostic message is provided.  one example 
*/  of a place were this occurs is gd158 from endf/b-vi.  reported 
*/  by frankle (lanl).
*i reconr.1588
      character*40 text
*d reconr.1716,1719
         write(text,'(''xsec nonzero at threshold for mt='',i3)') mth
         call mess('lunion',text,'adusted using jump in xsec')
*d reconr.1767
      er=sigfig(er,7,0)
*d reconr.1783
      enl=sigfig(er,7,0)
*d reconr.4204
*/ reconr -- 24sep01
*/  reconr contains some logic that tries to avoid doing work on
*/  very small charged-particle cross sections by defining a
*/  "pseudo-threshold" when the cross section rises to more than
*/  1e-15 barns.  however, this scheme isn't carried out completely,
*/  and it only results in the omission of the threshold energy for
*/  reactions that have less than this cross section just above
*/  the threshold.  this effect shows up for the (n,n't) reaction
*/  mt=33 for cd-110 from endf/b-vi.4.  at the request of bob
*/  little (lanl), we are changing the constant "ssmall" that
*/  triggers this effect to a smaller number.  in the long term,
*/  we should reconsider this logic.
*d reconr.1601
      ssmall=1.d-30
*d reconr.1614
      ssmall=1.e-30
*ident up61
*/ acer -- 25sep01
*/  kisako kazuaki (sumimoto) has observed that the common variable
*/  coeff in eval is not set.  actually, eval is not really used in
*/  tabize anymore.  it is just leftover as an intialization for a(iy).
*/  the only other place it is used is in islin2, and islin2 is not 
*/  called anymore!  This update removes these leftover remnants.
*d acer.1019
*d acer.1022
      external loada,finda,error,sigfig
*d acer.1128
         a(iy)=0
*d acer.1216,1263
*/  kazuaki also noticed that the photoatomic heating value was
*/  being stored in ev instead of mev and that the atomic number
*/  aw0 was not being set.
*i acer.14730
      aw0=c2h
*d acer.14839
         xss(lhnm-1+i)=heat/emev
*ident up62
*/ purr -- 28sep01
*/  lanl group x-5 has noted that the conditional heating cross 
*/  section in the mcnp probability tables is not quite what they
*/  expected.  we change the calculation here to get results that
*/  are consistently given as eV/reaction for lssf=0 and fluctuation
*/  factors for lssf=1. 
*i purr.79
      zero=0
*d purr.454,458
            if (sigu(2,1,1).ne.zero) h=h*a(n1+j+2*nbin)/sigu(2,1,1)
*d purr.461,465
            if (sigu(3,1,1).ne.zero) h=h*a(n1+j+3*nbin)/sigu(3,1,1)
*d purr.468,472
            if (sigu(4,1,1).ne.zero) h=h*a(n1+j+4*nbin)/sigu(4,1,1)
*d purr.477
               if (a(n1+j+nbin).ne.zero) a(l)=a(l)/a(n1+j+nbin)
*ident up63
*/ acer -- 15oct01
*/  add a capability to include delayed neutron data in the ace
*/  file as allowed for mcnp4c.
*d up57.10
      common/ace9/mt16,mt455
*i up57.12
      mt455=0
*i up57.14
               if (mfd.eq.5.and.mtd.eq.455) mt455=1
*d up57.18
      common/ace9/mt16,mt455
*d acer.2194
     &     (mf.eq.5.and.mt.gt.900)) then
*d acer.4645
      integer dndat,dnd,ptype,ploct
*d acer.4658,4661
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*i acer.4674
      dimension dntc(6)
*i acer.4680
      data shake/1.d8/
*i acer.4684
      data shake/1.e8/
*i acer.4743
      nnud=0
*d acer.4754
         if (mf.eq.1.and.mt.eq.455) kfis=2
*i acer.4869
               if (mta.eq.455) then
                  call listio(nin,0,0,a(iscr),nb,nw)
                  nnf=n1h
                  do i=1,nnf
                     dntc(i)=a(iscr+5+i)
                  enddo
               endif
*i acer.4883
               if (mta.eq.455) call reserv('nud',nw,inud,a)
*i acer.4885
               if (mta.eq.455) in=inud
*i acer.4904
               if (mta.eq.455) nnud=nw
*d acer.4907,4908
            if (mta.eq.452.and.kfis.eq.2) then
               mta=455
            else if (mta.eq.455) then
               mta=456
            else
               idone=1
            endif
*d acer.5474
c     ***after energy distributions
*i acer.5515
c
c     ***store delayed neutron data
      ndnf=0
      if (nnud.gt.0.and.mt455.eq.0) write(nsyso,
     &  '(/''  a delayed nubar section was found, but''/
     &     ''   no delayed neutron spectra were found:''/
     &     ''   delayed neutron data supressed'')')
      if (nnud.gt.0.and.mt455.eq.1) then
         write(nsyso,'(/''  adding delayed neutron data'')') 
         nud=next
         l=next-1
c
c        ***fission delayed nubar data
         do i=1,nnud
            l=l+1
            xss(l)=a(i-1+inud)
         enddo
         next=l+1

c        ***locate the delayed neutron data in file 5
         call findf(matd,5,455,nin)
         call contio(nin,0,0,a(iscr),nb,nw)
         ndnf=n1h
c
c        ***dndat block
         dndat=next
c
c        ***read through the section to load the dndat block
         lff=dndat
         do i=1,ndnf
            call tab1io(nin,0,0,a(iscr),nb,nw)
            law=l2h
            n=n2h
c           dndat entry
            xss(lff)=dntc(i)/shake
            xss(lff+1)=0
            xss(lff+2)=n
            do j=1,n
               xss(lff+2+j)=sigfig(a(iscr+6+2*j)/emev,7,0)
               xss(lff+2+n+j)=sigfig(a(iscr+6+2*j+1),7,0)
            enddo
            lff=lff+3+2*n
            if (law.eq.1) then
c              law=1
               call tab2io(nin,0,0,a(iscr),nb,nw)
               ne=n2h
               do ie=1,ne
                  call tab1io(nin,0,0,a(iscr),b,nw)
                  do while (nb.ne.0)
                     call moreio(nin,0,0,a(iscr),nb,nw)
                  enddo
               enddo
            else if (law.eq.5) then
c              law=5
               call tab1io(nin,0,0,a(iscr),nb,nw)
               call tab1io(nin,0,0,a(iscr),nb,nw)
               do while (nb.ne.0)
                  call moreio(nin,0,0,a(iscr),nb,nw)
               enddo
            endif
         enddo
         next=lff
c
c        ***ldnd block
         ldnd=next
         next=ldnd+ndnf
c
c        ***dnd block
         dnd=next
c
c        ***go back to the start of the sections
         call repoz(nin,-2)
         call findf(matd,5,455,nin)
         call contio(nin,0,0,a(iscr),nb,nw)
c
c        ***store the data
         do i=1,ndnf
            call tab1io(nin,0,0,a(iscr),nb,nw)
            law=l2h
            n=n2h
c           ldnd entry
            xss(ldnd-1+i)=next-dnd+1
c           dnd data
c           there is only one law per family
            xss(next)=0
            xss(next+1)=4
            xss(next+2)=10
            xss(next+3)=0
            xss(next+4)=2
            xss(next+5)=sigfig(xxmin/emev,7,0)
            xss(next+6)=sigfig(xxmax/emev,7,0)
            xss(next+7)=1
            xss(next+8)=1
            if (law.eq.1) then
c              law=1
               call tab2io(nin,0,0,a(iscr),nb,nw)
               ne=n2h
               xss(next+9)=0
               xss(next+10)=ne
               next=next+11
               lxx=next
               next=next+2*ne
               do ie=1,ne
                  call tab1io(nin,0,0,a(iscr),b,nw)
                  nn=n1h
                  mm=n2h
                  iint=nint(a(iscr+7)) 
                  xss(lxx+ie-1)=c2h
                  xss(lxx+ne+ie-1)=next-dnd+1
                  xss(next)=iint
                  xss(next+1)=mm
                  loc=iscr+nw
                  do while (nb.ne.0)
                     call moreio(nin,0,0,a(loc),nb,nw)
                     loc=loc+nw
                  enddo
                  l=next+1
                  sumup=0
                  do j=1,mm
                     xss(l+j)=sigfig(a(iscr+4+2*nn+2*j)/emev,7,0)
                     xss(l+j+2+3*mm)=
     &                 sigfig(a(iscr+4+2*nn+2*j)/emev,7,0)
                     xss(l+mm+j)=sigfig(a(iscr+4+2*nn+2*j+1)*emev,7,0)
                     xss(l+mm+j+2+3*mm)=
     &                 sigfig(a(iscr+4+2*nn+2*j+1)*emev,7,0)
                     xss(l+2*mm+j)=sigfig(sumup,9,0)
                     xss(l+2*mm+j+2+3*mm)=sigfig(sumup,9,0)
                     ll=iscr+4+2*nn+2*j
                     if (j.lt.mm.and.iint.eq.1) then
                        sumup=sumup+(a(ll+2)-a(ll))*a(ll+1)
                     else if (j.lt.mm.and.iint.eq.2) then
                        sumup=sumup+(a(ll+2)-a(ll))*(a(ll+3)+a(ll+1))/2
                     endif
                  enddo
                  next=next+2+3*mm
               enddo
            else if (law.eq.5) then
c              law=5
               call tab1io(nin,0,0,a(iscr),nb,nw)
               xxmin=a(iscr+8)
               xxmax=a(iscr+10)
               call tab1io(nin,0,0,a(iscr),nb,nw)
               loc=iscr+nw
               do while (nb.ne.0)
                  call moreio(nin,0,0,a(loc),nb,nw)
                  loc=loc+nw
               enddo
               nn=n1h
               mm=n2h
c              there is no incident energy dependence, we represent
c              this by two energies with duplicated distributions
               xss(next+9)=0
               xss(next+10)=2
               xss(next+11)=sigfig(xxmin/emev,7,0)
               xss(next+12)=sigfig(xxmax/emev,7,0)
               xss(next+13)=next+15-dnd+1
               xss(next+14)=next+15+2+3*mm-dnd+1
               iint=nint(a(iscr+7))
               xss(next+15)=iint
               xss(next+15+2+3*mm)=iint
               xss(next+16)=mm
               xss(next+16+2+3*mm)=mm
               l=next+16
               sumup=0
               do j=1,mm
                  xss(l+j)=sigfig(a(iscr+4+2*nn+2*j)/emev,7,0)
                  xss(l+j+2+3*mm)=sigfig(a(iscr+4+2*nn+2*j)/emev,7,0)
                  xss(l+mm+j)=sigfig(a(iscr+4+2*nn+2*j+1)*emev,7,0)
                  xss(l+mm+j+2+3*mm)=
     &              sigfig(a(iscr+4+2*nn+2*j+1)*emev,7,0)
                  xss(l+2*mm+j)=sigfig(sumup,9,0)
                  xss(l+2*mm+j+2+3*mm)=sigfig(sumup,9,0)
                  ll=iscr+4+2*nn+2*j
                  if (j.lt.mm.and.iint.eq.1) then
                     sumup=sumup+(a(ll+2)-a(ll))*a(ll+1)
                  else if (j.lt.mm.and.iint.eq.2) then
                     sumup=sumup+(a(ll+2)-a(ll))*(a(ll+3)+a(ll+1))/2
                  endif
               enddo
               next=next+15+2*(2+3*mm)
            endif
         enddo
      endif
*d acer.7378,7381
      integer dndat,dnd,ptype,ploct
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*d acer.8049,8053
      integer dndat,dnd,ptype,ploct,hpd,tyrh,sigh,andh,dlwh,yh
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*d acer.9746
      integer dndat,dnd,ptype,ploct
*d acer.9750,9753
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*d acer.9784,9797
     &  6x,''*                     *'',9x,''ndnf'',i10/
     &  6x,''*        njoy         *'',10x,''esz'',i10/
     &  6x,''*                     *'',11x,''nu'',i10/
     &  6x,''***********************'',10x,''mtr'',i10/
     &  39x,''lqr'',i10/39x,''tyr'',i10/38x,''lsig'',i10/
     &  39x,''sig'',i10/38x,''land'',i10//39x,''and'',i10/
     &  38x,''ldlw'',i10/39x,''dlw'',i10/39x,''gpd'',i10/
     &  38x,''mtrp'',i10/37x,''lsigp'',i10/38x,''sigp'',i10/
     &  37x,''landp'',i10/38x,''andp'',i10/37x,''ldlwp'',i10/
     &  38x,''dlwp'',i10/40x,''yp'',i10/39x,''fis'',i10/
     &  39x,''end'',i10/37x,''iurpt'',i10/39x,''nud'',i10/
     &  37x,''dndat'',i10/38x,''ldnd'',i10/39x,''dnd'',i10/
     &  37x,''ptype'',i10/38x,''ntro'',i10/
     &  37x,''ploct'',i10///6x,''hk---'',a70)')
     &  hz,aw0,tz,hd,hm,nxs(1),(nxs(i),i=3,8),
     &  (jxs(i),i=1,27),(jxs(i),i=30,32),hk
*i acer.10147
c
c     ***print delayed neutron data
      if (nud.gt.0) then
c
c        ***delayed nubar
         write(nsyso,'(''1''/'' delayed nubar data''/
     &                       '' ------------------''/)')
         l=nud
         j=nint(xss(l))
         write(nsyso,'(12x,''lnu =  '',i3,25x,''tabular nu'')') j
         l=l+1
         j=nint(xss(l))
         write(nsyso,'(12x,''nr ='',i4)') j
         if (j.ne.0) then
            write(nsyso,'(12x,''nbt(i=1,nr) = '',20i5)')
     &        (nint(xss(l+i)),i=1,j)
            l=l+j
            write(nsyso,'(12x,''int(i=1,nr) = '',20i5)')
     &        (nint(xss(l+i)),i=1,j)
            l=l+j
         endif
         l=l+1
         j=nint(xss(l))
         write(nsyso,'(12x,''ne ='',i4)') j
         write(nsyso,'(12x,''e(i=1,ne) =   '',1p,6e14.6/
     &     (12x,7e14.6))')
     &     (xss(l+i),i=1,j)
         l=l+j
         write(nsyso,'(12x,''nu(i=1,ne) =  '',1p,6e14.6/
     &     (12x,7e14.6))')
     &     (xss(l+i),i=1,j)
         l=l+j
c
c        ***precursor information
         write(nsyso,'(/'' precursor information''/
     &                  '' ---------------------'')')
         l=dndat
         do i=1,ndnf
            write(nsyso,'(/6x,''decay constant'',i3,'' of'',i3,
     &        '' (per shake) ='',1p,e13.5)') i,ndnf,xss(l)
            l=l+2
            j=nint(xss(l))
            write(nsyso,'(/6x,''delayed fraction'')')
            write(nsyso,'(12x,''ne ='',i4)') j
            write(nsyso,'(12x,
     &        ''e(i=1,ne) =   '',1p,6e14.6/(12x,7e14.6))')
     &        (xss(l+ii),ii=1,j)
            l=l+j
            write(nsyso,'(12x,
     &        ''p(i=1,ne) =   '',1p,6e14.6/(12x,7e14.6))')
     &        (xss(l+ii),ii=1,j)
            l=l+j+1
         enddo
c
c        ***precursor energy distributions
         write(nsyso,'(/
     &     '' delayed neutron energy distributions by precursor''/
     &     '' -------------------------------------------------'')')
         l=0
         k=3
         do i=1,ndnf
            nlaw=1
            loct=nint(xss(i-1+ldnd)+dnd-1)
            law=nint(xss(loct+1))
            if (law.eq.4) then
               l=l+1
               if (l.gt.1) write(nsyso,'(/)')
               if (l.gt.1) k=1
               write(nsyso,'(//
     &           ''   energy distribution for delayed neutrons from '',
     &           ''precursor '',i3,'' of'',i3)') i,ndnf
               write(nsyso,'(/
     &           ''    law ='',i2,i5,''st of'',i2,'' laws''/)')
     &           law,nlaw,nlaw
               k=k+3
               m=nint(xss(loct+3))
               loct=loct+3
               write(nsyso,'(8x,''probability of law'')')
               write(nsyso,'(12x,''nr ='',i4)') m
               if (m.ne.0) then
                  write(nsyso,'(12x,''nbt(i=1,nr) = '',20i5)')
     &              (nint(xss(j+loct)),j=1,m)
                  write(nsyso,'(12x,''int(i=1,nr) = '',20i5)')
     &              (nint(xss(j+m+loct)),j=1,m)
                  k=k+4
                  loct=loct+2*m
               endif
               loct=loct+1
               n=nint(xss(loct))
               write(nsyso,'(12x,''ne ='',i4)') n
               write(nsyso,'(12x,''e(i=1,ne) =   '',1p,6e14.6
     &           /(12x,7e14.6))') (xss(j+loct),j=1,n)
               write(nsyso,'(12x,''p(i=1,ne) =   '',1p,6e14.6
     &           /(12x,7e14.6))') (xss(j+n+loct),j=1,n)
               k=k+3
               loct=loct+1+2*n
               write(nsyso,'(/)')
               write(nsyso,'(8x,''data for law'')')
               k=k+2
               m=nint(xss(loct))
               write(nsyso,'(12x,''nr ='',i4)') m
               if (m.ne.0) then
                  k=k+4
                  write(nsyso,'(12x,''nbt(i=1,nr) = '',20i5)')
     &              (nint(xss(j+loct)),j=1,m)
                  write(nsyso,'(12x,''int(i=1,nr) = '',20i5)')
     &              (nint(xss(j+m+loct)),j=1,m)
                  loct=loct+2*m
               endif
               loct=loct+1
               ne=nint(xss(loct))
               write(nsyso,'(12x,''ne ='',i4)') ne
               if (m.eq.0) k=k+1
               do ie=1,ne
                  eg=xss(ie+loct)
                  loci=nint(xss(ie+ne+loct))+dnd-1
                  intt=nint(xss(loci))
                  n=nint(xss(loci+1))
                  loci=loci+1
                  if (ie.ne.1.and.k+6+n.ge.57) then
                     write(nsyso,'(/)')
                     k=1
                  endif
                  write(nsyso,'(/6x,'' incident energy = '',1p,e14.6,
     &              ''   intt ='',i2,''    np ='',i3//
     &              1x,
     &              2(''        energy           pdf           cdf'')/
     &              1x,
     &              2(''  ------------  ------------  ------------'')/
     &              (1x,1p,6e14.6))')
     &              eg,intt,n,(xss(j+loci),xss(j+n+loci),
     &              xss(j+2*n+loci),j=1,n)
                  k=k+n+6
               enddo
            endif      
         enddo
      endif
*i acer.10197
      integer dndat,dnd,ptype,ploct
*d acer.10198,10201
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*d acer.10671,10674
      integer dndat,dnd,ptype,ploct
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*d acer.11063
      integer dndat,dnd,ptype,ploct,hpd,tyrh,sigh,andh,dlwh,yh
*d acer.11064,11067
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*i acer.11579
      integer dndat,dnd,ptype,ploct
*d acer.11583,11587
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*d acer.11685,11688
      integer dndat,dnd,ptype,ploct
      common/nxst/len2,izaid,nes,ntr,nrx,ntrp,ntype,ndnf,nxsd(8)
      common/jxst/esz,nu,mtr,lqr,tyr,lsig,sig,land,and,ldlw,dlw,
     &  gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     &  iurpt,nud,dndat,ldnd,dnd,jxsd(2),ptype,ntro,ploct
*i acer.12303
c
c     ***delayed neutron block
      if (ndnf.ne.0) then
c        delayed nubar
         l=nud
         if (nout.ne.1) lnu=nint(xss(l))
         if (nout.eq.1) lnu=iss(l)
         call typen(l,nout,1)
         l=l+1
         if (nout.ne.1) nr=nint(xss(l))
         if (nout.eq.1) nr=iss(l)
         call typen(l,nout,1)
         l=l+1
         if (nr.ne.0) then
            n=2*nr
            do j=1,n
               call typen(l,nout,1)
               l=l+1
            enddo
         endif
         if (nout.ne.1) ne=nint(xss(l))
         if (nout.eq.1) ne=iss(l)
         call typen(l,nout,1)
         l=l+1
         n=2*ne
         do j=1,n
            call typen(l,nout,2)
            l=l+1
         enddo
c        precursor data
         l=dndat
         do i=1,ndnf
            call typen(l,nout,2)
            l=l+1
            if (nout.ne.1) nr=nint(xss(l))
            if (nout.eq.1) nr=iss(l)
            call typen(l,nout,1)
            l=l+1
            if (nr.ne.0) then
               n=2*nr
               do j=1,n
                  call typen(l,nout,1)
                  l=l+1
               enddo
            endif
            if (nout.ne.1) ne=nint(xss(l))
            if (nout.eq.1) ne=iss(l)
            call typen(l,nout,1)
            l=l+1
            n=2*ne
            do j=1,n
               call typen(l,nout,2)
               l=l+1
            enddo
         enddo
c        precursor energy distribution locators
         do i=1,ndnf
            call typen(l,nout,1)
            l=l+1
         enddo
c        precursor energy distributions
         do i=1,ndnf
            call typen(l,nout,1)
            l=l+1
            call typen(l,nout,1)
            l=l+1
            call typen(l,nout,1)
            l=l+1
            if (nout.ne.1) nr=nint(xss(l))
            if (nout.eq.1) nr=iss(l)
            call typen(l,nout,1)
            l=l+1
            if (nr.ne.0) then
               n=2*nr
               do j=1,n
                  call typen(l,nout,1)
                  l=l+1
               enddo
            endif
            if (nout.ne.1) ne=nint(xss(l))
            if (nout.eq.1) ne=iss(l)
            call typen(l,nout,1)
            l=l+1
            n=2*ne
            do j=1,n
               call typen(l,nout,2)
               l=l+1
            enddo
c           law=4 data
            if (nout.ne.1) nr=nint(xss(l))
            if (nout.eq.1) nr=iss(l)
            call typen(l,nout,1)
            l=l+1
            if (nr.gt.0) then
               n=2*nr
               do j=1,n
                  call typen(l,nout,1)
                  l=l+1
               enddo
            endif
            if (nout.ne.1) ne=nint(xss(l))
            if (nout.eq.1) ne=iss(l)
            call typen(l,nout,1)
            l=l+1
            do j=1,ne
               call typen(l,nout,2)
               l=l+1
            enddo
            do j=1,ne
               call typen(l,nout,1)
               l=l+1
            enddo
            do j=1,ne
               call typen(l,nout,1) 
               l=l+1
               if (nout.ne.1) np=nint(xss(l))
               if (nout.eq.1) np=iss(l)
               call typen(l,nout,1)
               l=l+1
               n=3*np
               do k=1,n
                  call typen(l,nout,2)
                  l=l+1
               enddo
            enddo
         enddo
      endif
*ident up64
*/ acer -- 07dec01
*/  there was a mistake introduced with up57 that only affects
*/  64-bit compiles (when "*set sw" is not used).  In addition, there
*/  were three mistakes related to the delayed neutron patch of up63.
*/  it's bad practice to use different names in common blocks in
*/  different places!
*d up56.16,17
*i acer.2135
      zero=0
      one=1
*d up63.255
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntyph,ndnf,nxsd(8)
*d up63.430
      common/nxst/len2,izaid,nes,ntr,nr,ntrp,ntyph,ndnf,nxsd(8)
*i up63.440
      common/ace7/awi,izai,mcnpx,newfor
*ident up65
*/ groupr -- 10dec01
*/  fix groupr to handle radionuclide production using the new
*/  endf file 8 and data from eaf2001.
*d groupr.501,502
      if (mfd.gt.36.and.mfd.lt.10000000) go to 381
*d groupr.516
      if (mfd.gt.10000000) izam=mod(mfd,10000000)
*d groupr.525
      if (mfd.le.10000000) go to 405
*d groupr.621,622
            if (mfd.lt.1000000) then
               write(nsyso,'('' for mf'',i3,'' and mt'',i3,1x,15a4)')
     &           mfd,mtd,(mtname(i),i=1,15)
            else
               mfdn=mfd/10000000
               jzam=mod(mfd,10000000)
               if (mfdn.eq.1) then
                  write(nsyso,
     &              '('' for mf3 mt'',i3,'' zam'',i8,1x,15a4)')
     &              mtd,jzam,(mtname(i),i=1,15)
               else if (mfdn.eq.2) then
                  write(nsyso,
     &              '('' for mf3*mf6 mt'',i3,'' zam'',i8,1x,15a4)')
     &              mtd,jzam,(mtname(i),i=1,15)
               else if (mfdn.eq.3) then
                  write(nsyso,
     &              '('' for mf3*mf9 mt'',i3,'' zam'',i8,1x,15a4)')
     &              mtd,jzam,(mtname(i),i=1,15)
               else if (mfdn.eq.4) then
                  write(nsyso,
     &              '('' for mf10 mt'',i3,'' zam'',i8,1x,15a4)')
     &              mtd,jzam,(mtname(i),i=1,15)
               endif
            endif
*i groupr.633
            if (mfd.gt.10000000) mfh=3
*i groupr.657
            if (mfd.gt.10000000) mfh=3
*d groupr.910
      mfd=mf10i(ir)
*d groupr.1163
      else if (mfd.gt.10000000) then
*d groupr.3123
      if (mfd.ne.3.and.mfd.ne.8.and.mfd.ne.18.and.mfd.lt.10000000) then
*d groupr.3857
      if (mft.eq.9.or.mft.eq.10) lfn=nint(a(iyld+3))
      if (mft.eq.6) lfn=nint(a(iyld+2))
*d groupr.3862,3871
      call skip6(itape,0,0,a(loc),law)
*d groupr.3859
      if (mft.gt.6.and.izn.eq.0.and.lfs.eq.lfn) go to 180
*d groupr.3975,3982
         if (mf.eq.10) then
            nfs=n1h
            jfs=-1
            do i=1,nfs
               call tab1io(nsig,0,0,a(isig),nb,nw)
               if (l2h.eq.lfs) jfs=i
               do while (nb.ne.0) 
                  call moreio(nsig,0,0,a(isig),nb,nw)
               enddo
            enddo
            if (jfs.lt.0) call error('getsig',
     &        'desired lfs not found',' ')
            nskip=jfs-1
            call skiprz(nsig,-1)
            call findf(matd,mf,mt,nsig)
            call contio(nsig,0,0,a(isig),nb,nw)
            if (nskip.gt.0) then
               do i=1,nskip
                  call tab1io(nsig,0,0,a(isig),nb,nw)
                  do while (nb.ne.0)
                     call moreio(nsig,0,0,a(isig),nb,nw)
                  enddo
               enddo
            endif
*d groupr.4345
      if (mfd.eq.12.or.(mfd.gt.20000000.and.mfd.lt.40000000))
*ident up66
*/ groupr -- 12feb02
*/  if you mix automatic reactions with manual reactions where
*/  the mtname string is not given, the mtname on the manual
*/  reaction will have whatever string was left from the previous
*/  case.  we fix that here.
*d groupr.507
*/ groupr -- 13feb02
*/  there is an error in the calculation of the kalbach a factor
*/  for the photonuclear case.  it is necessary to convert e to
*/  mev for this formula.  the symptom is results for a that are
*/  so large that sinh(a) overflows with a floating point error.
*d groupr.5980
         bb=bb*sqrt((tomev*e)/(2*emc2))*fact
*ident up67
*/ matxsr -- 12feb02
*/  add photonuclear capability
*i matxsr.30
c     *   ngen8     photonuclear data from groupr (default=0)          *
*d matxsr.210
cd                     gscat     gamma scattering (atomic)             -
cd                     gg        gamma scattering (photonuclear)       -
*d matxsr.392
     &  nscrt5,nscrt6,nscrt7,ngen3,ngen4,ngen5,ngen6,ngen7,ngen8
*d matxsr.422
      ngen8=0
      read(nsysi,*)
     &  ngen1,ngen2,nmatx,ngen3,ngen4,ngen5,ngen6,ngen7,ngen8
*d matxsr.433,434
     &  '' incident alpha unit .................. '',i10/
     &  '' photonuclear unit .................... '',i10)')
     &  ngen3,ngen4,ngen5,ngen6,ngen7,ngen8
*d matxsr.441
      nscrt8=17
*i matxsr.448
      call openz(ngen8,0)
*i matxsr.475
      call closz(ngen8)
*d matxsr.495
     &  nscrt5,nscrt6,nscrt7,ngen3,ngen4,ngen5,ngen6,ngen7,ngen8
*d matxsr.892
     &  nscrt5,nscrt6,nscrt7,ngen3,ngen4,ngen5,ngen6,ngen7,ngen8
*i matxsr.909
      character*8 hgg
*d matxsr.917
      data hgsct/'gscat '/, hgg/'gg    '/, hnthr/'ntherm'/
*i matxsr.945
      if (nin.eq.0) nin=ngen8
*d matxsr.1012
      if (hprt(ip1).eq.hgm) nin=ngen8
      if (hprt(ip1).eq.hgm.and.htyp.eq.hgsct) nin=ngen2
*d matxsr.1027,1028
      if (hprt(ip2).eq.hgm.and.htyp.ne.hgg) mfv=13
      if (hprt(ip1).eq.hgm.and.htyp.eq.hgsct) mfv=23
*i matxsr.1035
      if (hprt(ip2).eq.hgm.and.htyp.eq.hgg) mfm=16
*d matxsr.1478
     &  nscrt5,nscrt6,nscrt7,ngen3,ngen4,ngen5,ngen6,ngen7,ngen8
*d matxsr.1809
     &  nscrt5,nscrt6,nscrt7,ngen3,ngen4,ngen5,ngen6,ngen7,ngen8
*ident up68
*/ reconr -- 18apr02
*/  the name "pi" is being used for both the imaginary part of
*/  p and the pi constant.  the former is changed to "pim" as
*/  used in the version of this subroutine in purr.  the symptom
*/  is a floating point error when sqrt(pi) is calculated with
*/  pi changed to a negative value.
*d reconr.4868
      pim=aimw
*d reconr.4873
      if (abs(aimw-pim).ge.eps) go to 380
*d reconr.4908
      pim=aimw
*d reconr.4912
      if (abs(aimw-pim).ge.eps) go to 470
*/  allow for more digits in the temperature printout to
*/  handle the usage of reconr being made for eaf-2001.  the
*/  temperature field is not being used for resonance 
*/  reconstruction but only passed to the output file for
*/  later use.
*d reconr.348
     &  '' reconstruction temperature ........... '',f10.2,''k''/
*ident up69
*/ acer -- 07dec01
*/  add a capability to generate fluorescence data for mcnp using
*/  the existing cashwell-everett format with new numbers obtained
*/  from the endf versions of eadl and epdl.  the data produced with
*/  this method should give reasonable results for transport and
*/  heating for energies above 1 kev.  the new evaluations allow
*/  for lower incident photon energies and for more detail in photon
*/  and electron distributions from photoabsorption, and future
*/  versions of njoy and mcnp can take advantage of this.
*d acer.66
c     * data.  the input photoatomic data is mounted on nendf.         *
c     * fluorescence data can be generated from atomic relaxation      *
c     * data in endf format mounted on npend.                          * 
*i acer.186
c     *             photoatomic data on nendf                          *
c     *             atomic relaxation data on npend                    *
*d acer.428
         call acepho(nendf,npend)
*d acer.14651
      subroutine acepho(nin,nlax)
*d acer.14690
      data emax/1.01d11/
*d acer.14702
      data emax/1.01e11/
*i acer.14780
c
c     ***set number of fluorescence lines
      iz=matd/100
      if (iz.lt.12) then
         nflo=0
      else if (iz.lt.20) then
         nflo=2
      else if (iz.lt.31) then
         nflo=4
      else if (iz.lt.37) then
         nflo=5
      else 
         nflo=6
      endif
*d acer.14786
      lhnm=jflo+4*nflo
*d acer.14843
c     ***for fluorescence photons
*d acer.14845,14846
      if (nlax.eq.0) then
         call mess('acepho','no atomic relaxation data',
     &     'fluorescence data not processed')
      else if (nflo.gt.0) then
         call alax(nin,nlax,xss(jflo),a(iscr))
      endif
*i acer.14997
c
      subroutine alax(nin,nlax,fluor,a)
c     ******************************************************************
c     generate fluorescence data in the cashwell-everett format.
c     ******************************************************************
      implicit real*8 (a-h,o-z)
      common/util/npage,iverf
      common/cont/c1h,c2h,l1h,l2h,n1h,n2h,math,mfh,mth,nsh,nsp,nsc
      common/mainio/nsysi,nsyso,nsyse,ntty
      common/ace1/tempd,err,matd,nbina,nbinp,negn,iprint,iopt,ndigit
      dimension fluor(*)
      dimension a(*)
      dimension loc(50)
      dimension enl(3),rhol(3),wtl(3)
*if sw
      data dn/.9999d0/
      data up/1.0001d0/
      data emev/1.d6/
*else
      data dn/.9999e0/
      data up/1.0001e0/
      data emev/1.e6/
*endif
c
c     ***charge for desired material
      iz=matd/100
c
c     ***read in the atomic relaxation file for the desired material
      call openz(nlax,0)
      call tpidio(nlax,0,0,a,nw,nb)
  110 call contio(nlax,0,0,a,nw,nb)
      if (math.gt.0) go to 120
      call error('alax','mat not found',' ')
  120 if (math.eq.matd) go to 130
      call tomend(nlax,0,0,a)
      go to 110
  130 call tofend(nlax,0,0,a)
      call contio(nlax,0,0,a,nw,nb)
      nss=n1h
      ll=1
      do iss=1,nss
         loc(iss)=ll
         call listio(nlax,0,0,a(ll),nb,nw)
         ntr=n2h
         ll=ll+nw
         do while (nb.ne.0)
            call moreio(nlax,0,0,a(ll),nb,nw)
            ll=ll+nw
         enddo
      enddo
c
c     ***read in the photoionization cross section for the material
      kk=ll
      call openz(nin,0)
      call tpidio(nin,0,0,a(ll),nw,b)
  210 call contio(nin,0,0,a(ll),nw,b)
      if (math.gt.0) go to 220
      call error('spect','mat not found',' ')
  220 if (math.eq.matd) go to 230
      call tomend(nin,0,0,a(ll))
      go to 210
  230 call tofend(nin,0,0,a(ll))
      call findf(matd,23,522,nin)
      call contio(nin,0,0,a(ll),nw,nb)
      e=0
      call gety1(e,en,idis,sig,nin,a(ll))
c
c     ***for z>30, get the l1, l2, and l3 edges and jumps
      if (iz.gt.30) then
         do i=1,3
            jj=loc(5-i)
            enl(4-i)=a(jj+6)
            e=dn*enl(4-i)
            call gety1(e,en,idis,slo,nin,a(ll))
            e=up*enl(4-i)
            call gety1(e,en,idis,shi,nin,a(ll))
            rhol(4-i)=slo/shi
         enddo
      endif
c
c     ***get the energy and jump of the k edge
      ek=a(7)
      e=dn*ek
      call gety1(e,en,idis,slo,nin,a(ll))
      e=up*ek
      call gety1(e,en,idis,shi,nin,a(ll))
      rhok=slo/shi
c
c     ***case of 11<z<20
      if (iz.gt.11.and.iz.lt.20) then
c
c        ***average all the transitions to the k shell
         n=nint(a(6))
         sum1=0
         sum2=0
         do i=1,n
            jj=8+6*i
            if (nint(a(jj)).eq.0) then
               sum1=sum1+a(jj+2)
               sum2=sum2+a(jj+1)*a(jj+2)
            endif
         enddo
         sum2=sum2/sum1
c
c        ***store the results
         fluor(1)=ek/emev
         fluor(2)=ek/emev
         fluor(3)=rhok
         fluor(4)=1
         fluor(5)=0
         fluor(6)=(1-rhok)*sum1
         fluor(7)=0
         fluor(8)=sum2/emev
c
c     ***case of z>19 and z<31
      else if (iz.gt.19.and.iz.lt.31) then
c
c        ***extract l2, l3, and total for higher shells
         n=nint(a(6))
         sum1=0
         sum2=0
         do i=1,n
            jj=8+6*i
            if (nint(a(jj)).eq.0.and.nint(a(jj-1)).eq.3) then
               el2=a(jj+1)
               pl2=a(jj+2)
            else if (nint(a(jj)).eq.0.and.nint(a(jj-1)).eq.4) then
               el3=a(jj+1)
               pl3=a(jj+2)
            else if (nint(a(jj)).eq.0.and.nint(a(jj-1)).gt.4) then
               sum1=sum1+a(jj+2)
               sum2=sum2+a(jj+1)*a(jj+2)
            endif
         enddo
         sum2=sum2/sum1
c
c        ***store the results
         tot=(pl2+pl3+sum1)/(1-rhok)
         y=0
         phi=rhok
         fluor(1)=ek/emev
         fluor(5)=phi
         fluor(9)=y
         fluor(13)=0
         phi=phi+pl3/tot
         y=y+(1-rhok)*pl3
         fluor(2)=ek/emev
         fluor(6)=phi
         fluor(10)=y
         fluor(14)=el3/emev
         phi=phi+pl2/tot
         y=y+(1-rhok)*pl2
         fluor(3)=ek/emev
         fluor(7)=phi
         fluor(11)=y
         fluor(15)=el2/emev
         phi=1
         y=y+(1-rhok)*sum1
         fluor(4)=ek/emev
         fluor(8)=phi
         fluor(12)=y
         fluor(16)=sum2/emev
c
c     ***all other z values
      else
         rholt=rhol(1)*rhol(2)*rhol(3)
         elav=(enl(1)+enl(2)+enl(3))/3
         wtl(1)=1/rhol(1)
         wtl(2)=wtl(1)/rhol(2)
         wtl(3)=wtl(2)/rhol(3)
         denom=wtl(3)-1
         wtl(3)=(wtl(3)-wtl(2))/denom
         wtl(2)=(wtl(2)-wtl(1))/denom
         wtl(1)=(wtl(1)-1)/denom
c
c        ***compute the average yield and energy for l fluorescence
         sum1=0
         sum2=0
         do iss=2,4
            jj=loc(iss)
            n=nint(a(jj+5))
            wt=wtl(iss-1)
            do i=1,n
               if (nint(a(jj+7+6*i)).eq.0) then
                  sum1=sum1+a(jj+9+6*i)*wt
                  sum2=sum2+a(jj+8+6*i)*a(jj+9+6*i)*wt
               endif
            enddo
         enddo
         sum2=sum2/sum1
         ylt=sum1
         flt=sum2
         if (flt.gt.enl(1)) then
            write(nsyso,'('' L edge problem'')')
            write(nsyso,'(1x,3f10.4)') flt,enl(1),elav
         endif
c
c        ***extract kalpha1, kalpha2, kbeta1, and kbeta2
         n=nint(a(6))
         sum11=0
         sum12=0
         sum21=0
         sum22=0
         do i=1,n
            jj=8+6*i
            if (nint(a(jj)).eq.0) then
               mm=nint(a(jj-1))
               if (mm.eq.3) then
                  el2=a(jj+1)
                  pl2=a(jj+2)
               else if (mm.eq.4) then
                  el3=a(jj+1)
                  pl3=a(jj+2)
               else if (mm.ge.5.and.mm.le.9) then
                  sum11=sum11+a(jj+2)
                  sum12=sum12+a(jj+1)*a(jj+2)
               else if (mm.ge.10.and.mm.le.16) then
                  sum21=sum21+a(jj+2)
                  sum22=sum22+a(jj+1)*a(jj+2)
               endif
            endif
         enddo
         if (iz.ge.37) then
            sum22=sum22/sum21
         else
            sum11=sum11+sum21
            sum12=sum12+sum22
            sum21=0
         endif
         sum12=sum12/sum11
c
c        ***store the results
         n=5
         if (iz.gt.36) n=6
         fluor(1)=elav/emev
         fluor(1+n)=rholt
         fluor(1+2*n)=0
         fluor(1+3*n)=0
         y=(1-rholt)*ylt
         fluor(2)=elav/emev
         fluor(2+n)=1
         fluor(2+2*n)=y
         fluor(2+