*cpl njoy,reconr,broadr,unresr,heatr,thermr,groupr,errorr
*cpl covr,moder,ccccr,matxsr,acer,viewr,purr,gaspr,wimsr
*cpl plotr,leapr
*set sw
*ident up1
*/ reconr -- 04jul98 
*/   fix an error made in enhancing the equality tests in 
*/   reconr when building njoy97.  this error causes the code 
*/   to miss discontinuities, such as those at the ends of 
*/   resonance ranges, where it normally tries to put in points 
*/   just below and just above the discontinuity.  problem 
*/   noticed by little (lanl).
*d reconr.1456
      if (abs(srnext-sr).lt.small*sr) go to 260
*/ reconr -- 04jul98
*/   change the wide constant to prevent infinite loops
*/   (this showed up for u-232).
*d reconr.2808
      data wide/3.01d0/
*d reconr.2814
      data wide/3.01/
*d reconr.3018
      data wide/3.01d0/
*d reconr.3024
      data wide/3.01/
*/ reconr -- 02jul98 
*/   fix up some mistakes in the use of the cmplx intrinsic 
*/   function, fix some variables left as real when they should 
*/   be double, and make sure that the complex variables go to 
*/   double precision properly when requested.
*d reconr.2582,2583
      dimension qre(4),gc(4),sigc(4)
*if sw
      real*8 k,kfac
      complex*16 one,eye,rfcnn,slsjnn,phi0,phi,r0,exphi
*else
      real k,kfac
      complex one,eye,rfcnn,slsjnn,phi0,phi,r0,exphi
*endif
*d reconr.2588
      data zero,uno/0.d0,1.d0/
*d reconr.2592
      data zero,uno/0.,1./
*d reconr.2595,2596
*if sw
      eye=dcmplx(zero,uno)
      one=dcmplx(uno,zero)
*else
      eye=cmplx(zero,uno)
      one=cmplx(uno,zero)
*endif
*d reconr.2668
*if sw
      phi0=dcmplx(phir,zero)
*else
      phi0=cmplx(phir,zero)
*endif
*d reconr.2673
*if sw
      rfcnn=dcmplx(zero,zero)
*else
      rfcnn=cmplx(zero,zero)
*endif
*d reconr.2729
*if sw
      r0=dcmplx(zero,zero)
*else
      r0=cmplx(zero,zero)
*endif
*d reconr.2741
*if sw
         r0=dcmplx(rr0,ri0)
*else
         r0=cmplx(rr0,ri0)
*endif
*d reconr.2758
*if sw
         phi=dcmplx(phr0,phri)
*else
         phi=cmplx(phr0,phri)
*endif
*d reconr.2762
      exphi=exp(-2*eye*phi)
*d reconr.2765,2766
*if sw
      rpart=dble(slsjnn)
*else
      rpart=real(slsjnn)
*endif
      sabs=abs(slsjnn)
      sigp(1)=sigp(1)+2*pifac*gj*rpart
*/ reconr -- 02jul98
*/   fix up the use of float.  for double precision version, 
*/   float will put out a real*4, which is not what we want
*/   when the user has asked for "*set sw".
*d reconr.2040
      spot=4*(2*ll+1)*pifac*sinsq
*d reconr.2201
      fl=ll
*d reconr.2348
      fl=ll
*d reconr.3123
      spot=4*pi*(2*ll+1)*(sin(ps)/k)**2
*/ reconr -- 02jul98
*/   fix some places where mixed integer and real variables 
*/   could cause problems.
*d reconr.1309
      l=iscr+6+2*nint(a(iscr+4))
*d reconr.2686
      lc(ii)=nint(a(inow+ii+4))
*d reconr.3100
      jnow=nint(a(inow+4))+inow+6
*ident up2
*/ broadr -- 4jul98
*/   patches an error in the calculation of the thermal
*/   quantities involving fission, and then further
*/   improves the energy integrals so that the values 
*/   computed in broadr match those that would be 
*/   computed by groupr.  we also print out the
*/   1/v-equivalent version of k1, that is, the k1
*/   integral times 2 over sqrt(pi).
*i broadr.118
      data pi/3.14159265d0/
*i broadr.120
      data pi/3.14159265/
*i broadr.453
      picon=2/sqrt(pi)
*d broadr.473
 2240 step=1.01
      enow=tt(1)
      if (enow.gt.step*elast.and.elast.gt.0.) then
         enow=step*elast
         sf=slf+(tt(llf)-slf)*(enow-elast)/(tt(1)-elast)
         sc=slc+(tt(llc)-slc)*(enow-elast)/(tt(1)-elast)
      else
         sf=tt(llf)
         sc=tt(llc)
      endif
*d broadr.477
            ss=slf+(therm-elast)*(sf-slf)/(enow-elast)
*d broadr.491
            ss=slc+(therm-elast)*(sc-slc)/(enow-elast)
*d broadr.499
            ss=slf+(tev-elast)*(sf-slf)/(enow-elast)
*d broadr.514
            ss=slc+(tev-elast)*(sc-slc)/(enow-elast)
*i broadr.533
         fnow=enow*exp(-enow/tev)
*d broadr.535
*d broadr.537
               cint=cint+(fnow*sc+flast*slc)*(enow-elast)/2
*d broadr.539
*d broadr.542
*d broadr.544,551
               fint=fint+(fnow*sf+flast*slf)*(enow-elast)/2
               if (sf.gt.0..and.slf.gt.0.)
     1           alint=alint+(fnow*sc/sf+flast*slc/slf)
     2           *(enow-elast)/2
               etint=etint+(fnow*fnu*sf/(sf+sc)
     1           +flast*fnul*slf/(slf+slc))*(enow-elast)/2
               v1int=v1int+(fnow*(fnu*sf-sf-sc)
     1           +flast*(fnul*slf-slf-slc))*(enow-elast)/2

*d broadr.553
*i broadr.554
         flast=fnow
*d broadr.559,560
            ss=slf+(eone-elast)*(sf-slf)/(enow-elast)
            ssf=(sf/enow+ss/elast)*(enow-eone)/2
*d broadr.563,564
            ss=slc+(eone-elast)*(sc-slc)/(enow-elast)
            ssc=(sc/enow+ss/elast)*(enow-eone)
*d broadr.568,571
         if (llf.gt.0) ssf=ssf+(sf/enow+slf/elast)*(enow-elast)/2
         if (llc.gt.0) ssc=ssc+(sc/enow+slc/elast)*(enow-elast)/2
*d broadr.576,577
      if (llf.gt.0) slf=sf
      if (llc.gt.0) slc=sc
      if (enow.lt.tt(1)) go to 2240
*d broadr.586
     2 picon*cint/ctev
*d broadr.598
     2    picon*fint/ftev
*i broadr.604
      write(nsyso,'(
     1    ''                  equivalant k1:'',1p,e12.4)') 
     2     picon*v1int
*ident up3
*/ heatr -- 02jul98
*/   fix some variables inadvertantly typed as real.  
*/   they should become double if "*set sw" was used.
*d heatr.3263
*if sw
      real*8 nc,nb
*else
      real nc,nb
*endif
*/ heatr -- 02jul98
*/   fix some places where mixed integer and real variables
*/   could cause problems.
*d heatr.1746,1747
      nr=nint(c(lnow+4))
      np=nint(c(lnow+5))
*d heatr.2749
      l=7+2*nint(c(5))
*/ heatr -- 02jul98
*/   change the name of the 939 mev constant.  it should
*/   an implicit floating type.
*d heatr.654
      data emc2/939.512d6/
*d heatr.659
      data emc2/939.512e6/
*d heatr.732
      tm=emc2*(awr+1)
*d heatr.1245
      data emc2/939.512d6/
*d heatr.1253
      data emc2/939.512e6/
*d heatr.1265
         tm=emc2*(awr+1)
*d heatr.3444
      data emc2/939.512d6/
*d heatr.3448
      data emc2/939.512e6/
*d heatr.3453
      ein=2*emc2*(awr+1)
*d heatr.4108
      data emc2/939.512d6/
*d heatr.4111
      data emc2/939.512e6/
*d heatr.4120
      tm=emc2*(awr+1)
*d heatr.4546
      data emc2/939.512d6/
*d heatr.4550
      data emc2/939.512e6/
*d heatr.4555
      ein=2*emc2*(awr+1)
*d heatr.4605
      data emc2/939.512d6/
*d heatr.4607
      data emc2/939.512e6/
*d heatr.4611
      tm=emc2*(awr+1)
*ident up4
*/ groupr -- 19jun98
*/   fix a bad write format and a bad message.
*d groupr.1917
     1  '' gamma group structure......read in'')')
*d groupr.3733
     1  call error('stounr','storage exceeded.',' ')
*/ groupr -- 02jul98
*/   fix problems with complex arithmetic when double 
*/   precisioning is used.
*d groupr.6846
*if sw
      complex*16 carg1,carg2,cs1,cs2
*else
      complex carg1,carg2,cs1,cs2
*endif
*i groupr.6887
      data zero,uno/0.d0,1.d0/
*i groupr.6925
      data zero,uno/0.,1./
*d groupr.7024
*if sw
      cs1=dcmplx(c(8+np),c(9+np))/2
*else
      cs1=cmplx(c(8+np),c(9+np))/2
*endif
*d groupr.7026
*if sw
      cs1=cs1+(2*it+1)*p(it+1)*dcmplx(c(8+np+2*it),c(9+np+2*it))/2
*else
      cs1=cs1+(2*it+1)*p(it+1)*cmplx(c(8+np+2*it),c(9+np+2*it))/2
*endif
*d groupr.7028,7029
*if sw
      carg1=dcmplx(zero,uno)*eta*log((1-x)/2)
      sigi=(-2*eta/(1-x))*dble(cs1*exp(carg1))  
*else
      carg1=cmplx(zero,uno)*eta*log((1-x)/2)
      sigi=(-2*eta/(1-x))*real(cs1*exp(carg1))  
*endif
*d groupr.7037
*if sw
      cs1=dcmplx(c(8+nt),c(9+nt))/2
*else
      cs1=cmplx(c(8+nt),c(9+nt))/2
*endif
*d groupr.7041,7042
*if sw
      cs1=cs1+(2*it+1)*p(it+1)*dcmplx(c(8+nt+2*it),c(9+nt+2*it))/2
      cs2=cs2+sgn*(2*it+1)*p(it+1)*dcmplx(c(8+nt+2*it),c(9+nt+2*it))/2
*else
      cs1=cs1+(2*it+1)*p(it+1)*cmplx(c(8+nt+2*it),c(9+nt+2*it))/2
      cs2=cs2+sgn*(2*it+1)*p(it+1)*cmplx(c(8+nt+2*it),c(9+nt+2*it))/2
*endif
*d groupr.7045,7048
*if sw
      carg1=dcmplx(zero,uno)*eta*log((1-x)/2)
      carg2=dcmplx(zero,uno)*eta*log((1+x)/2)
      sigi=(-2*eta/(1-x*x))*dble(cs1*(1+x)*exp(carg1)
     1   +cs2*(1-x)*exp(carg2))
*else
      carg1=cmplx(zero,uno)*eta*log((1-x)/2)
      carg2=cmplx(zero,uno)*eta*log((1+x)/2)
      sigi=(-2*eta/(1-x*x))*real(cs1*(1+x)*exp(carg1)
     1   +cs2*(1-x)*exp(carg2))
*endif
*/ groupr -- 02jul98
*/   fix some variables inadvertantly typed as real.  
*/   they should become double if desired.
*d groupr.5413
*if sw
      real*8 nc,nb
*else
      real nc,nb
*endif
*/ groupr -- 02jul98
*/   fix some places where mixed integer and real variables 
*/   could cause problems.
*d groupr.4551
      npsx=nint(c(l+5))
*d groupr.4645
      law=nint(c(l1+3))
*d groupr.6367
      nbrag=nint(p(6))
*d groupr.6753,6754
      il=iflg+nint(a(iloca+i-1))
      im=iflg+nint(a(iloca+nipj-1))
*d groupr.7158,7160
      lp=nint(a(l+2))
      nr=nint(a(l+4))
      np=nint(a(l+5))
*d groupr.7913
      nktot=nint(c(ic+4))
*/ groupr -- 02jul98 
*/   change the name of the mev constant.
*d groupr.1803
      data emev/1.d6/
*d groupr.1833
      data emev/1.e6/
*d groupr.1854
      egg(ig)=eg2(ig)*emev
*d groupr.1861
      egg(ig)=eg3(ig)*emev
*d groupr.1868
      egg(ig)=eg4(ig)*emev
*d groupr.1875
      egg(ig)=eg5(ig)*emev
*d groupr.1882
      egg(ig)=eg6(ig)*emev
*d groupr.1898
      egg(ig)=eg8(ig)*emev
*d groupr.1903
      egg(ig)=eg8(ig+1)*emev
*ident up5
*/ acer -- 04jul98 
*/   fix an uninitialized variable that always used to come
*/   through as zero on the type1 xsdir "stub" cards.
*d acer.7902
      lrec=0
*/ acer -- 22may98
*/   improve the printout, especially for secondary photons.
*/   there are no differences in calculated values caused
*/   by this change.  requested by frankle, lanl/x-ci.
*d acer.7061
     1  ''   intt ='',i2,''    nd ='',i4,''    np ='',i3//
*d acer.7065
     5  e2,intt,nd,nn,(xss(j+loci),xss(j+nn+loci),
*d acer.7186
     1  ''   intt ='',i2,''    nd ='',i4,''    np ='',i3//
*d acer.7192
     7  e2,intt,nd,nn,(xss(j+loci),xss(j+nn+loci),
*i acer.7586
      mtrn=nint(xss(i-1+mtrp))
      call mtname(mti,title(1))
      nlaw=1
*d acer.7593
      write(nsyso,'(//
     1  '' energy distribution for secondary photons from '',
     2  ''reaction '',a10,'' (mt='',i6,'') with '',i1,'' law(s)'')')
     3  title(1),mtrn,nlaw
      write(nsyso,'(/
     1  ''    law ='',i2,i5,''st of'',i2,'' laws''/)')
     2  law,nlaw,nlaw
*i acer.7631
      if (iprint.eq.1.and.m.ne.0) write(nsyso,'(12x,''ne ='',i4)') ne
*d acer.7647
     1  ''   intt ='',i2,''    nd ='',i3,''    np ='',i3//
*d acer.7651
     5  eg,intt,nd,n,(xss(j+loci),xss(j+n+loci),xss(j+2*n+loci),j=1,n)
*/ acer -- 22may98
*/   fix reaction names on plots for inelastic cross sections.  
*/   we're having trouble with quotes insides quotes, so we will 
*/   leave the classical ace stars, e.g., (n,n*2).
*/   fix one bad y-axis label where the quotes were omitted.
*/   fix the line type for the resonance fission plots.
*d acer.10844,10846
*d acer.10881,10883
*d acer.11019,11021
*d acer.11563,11565
*d acer.11882,11884
*d acer.12008,12010
*d acer.12103,12105
*d acer.12210,12212
*d acer.12266,12268
*d acer.12359,12361
*d acer.12728,12730
*d acer.11626
      write(nout,'(a,''<c>ross section (barns)'',a,''/'')') qu,qu 
*d acer.11366
*i acer.11367
      write(nout,'(''0 0 1/'')')
*/ acer -- 19jun98 
*/   be sure to compute enough legendre polynomials for 
*/   the required limits.
*d acer.2452
      call legndr(u,pl,nmax+1)
*/ acer -- 19jun98 
*/   fix a bad message format.
*d acer.3282
     1  '' ---message from summer---for distribution as per '',
*/ acer -- 02jul98 
*/   fix some variables inadvertantly typed as real.  they 
*/   should become double if desired.
*d acer.8902
*if sw
      real*8 nc,nb
*else
      real nc,nb
*endif
*/ acer -- 02jul98
*/   change the name of the mev constant.
*d acer.4577
      data emev/1.d6/
*d acer.4581
      data emev/1.e6/
*d acer.4680
      a(i+j)=sigfig(a(2*i+2*m+4+iscr)/emev,7,0)
*d acer.4805
      xss(lqr+ir)=sigfig(a(iscr+1)/emev,7,0)
*d acer.4828
  230 xss(ih+j)=sigfig(s/emev/xss(it+j),7,0)
*d acer.4833
      if (mth.eq.444) s=s/emev
*d acer.4950
      xss(ie+j)=sigfig(a(iscr+1)/emev,7,0)
*d acer.5010,5011
      xss(next+5)=sigfig(xss(esz+ja)/emev,7,0)
      xss(next+6)=sigfig(xss(esz+jb)/emev,7,0)
*d acer.5053
      xss(next-1+ie)=sigfig(a(jj)/emev,7,0)
*d acer.5106
      a(iegn-1+i)=sigfig(a(iscr+1)/emev,7,0)
*d acer.5108
      xss(next)=sigfig(a(iscr+5+j)/emev,7,0)
*d acer.5126
      xss(esz+i)=sigfig(xss(esz+i)/emev,9,0)
*d acer.5157
      data emev/1.d6/
*d acer.5159
      data emev/1.e6/
*d acer.5201
      xss(j+next)=sigfig(a(2*j+l)/emev,7,0)
*d acer.5252
      xss(next+j)=sigfig(e/emev,7,0)
*d acer.5271,5272
      xss(ki+nexd)=sigfig(a(iscr+6+2*ki)/emev,7,0)
      xss(ki+n+nexd)=sigfig(a(iscr+7+2*ki)*emev,7,0)
*d acer.5310,5311
      xss(j+next)=sigfig(a(2*j+2*m+4+iscr)/emev,7,0)
      xss(j+n+next)=sigfig(a(2*j+2*m+5+iscr)/emev,7,0)
*d acer.5321
      xss(j+next)=sigfig(a(2*j+2*m+4+iscr)/emev,7,0)
*d acer.5346,5347
      xss(j+next)=sigfig(a(2*j+2*m+4+iscr)/emev,7,0)
      xss(j+n+next)=sigfig(a(2*j+2*m+5+iscr)/emev,7,0)
*d acer.5350
      xss(next)=u/emeV
*d acer.5364,5365
      xss(j+next)=sigfig(a(2*j+4+iscr)/emev,7,0)
      xss(j+m+next)=sigfig(a(2*j+5+iscr)/emev,7,0)
*d acer.5373,5374
      xss(j+next)=sigfig(a(2*j+2*m+4+iscr)/emev,7,0)
      xss(j+n+next)=sigfig(a(2*j+2*m+5+iscr)/emev,7,0)
*d acer.5386,5387
      xss(j+next)=sigfig(a(2*j+2*m+4+iscr)/emev,7,0)
      xss(j+n+next)=sigfig(a(2*j+2*m+5+iscr)/emev,7,0)
*d acer.5395,5396
      xss(j+next)=sigfig(a(2*j+2*m+4+iscr)/emev,7,0)
      xss(j+n+next)=sigfig(a(2*j+2*m+5+iscr)*emev,7,0)
*d acer.5399
      xss(next)=sigfig(u/emev,7,0)
*d acer.5428
      xs(2)=emev
*d acer.5459
      xss(next+ie)=sigfig(tme/emev,7,0)
*d acer.5466,5467
      xss(ki+nexd)=sigfig(xxs(ki)/emev,7,0)
      xss(ki+js+nexd)=sigfig(yys(ki)*renorm*emev,7,0)
*d acer.5505
      data emev/1.d6/
*d acer.5507
      data emev/1.e6/
*d acer.5631
      xss(next+1+j)=sigfig(eyl/emev,7,0)
*d acer.5683
      xss(next+k)=sigfig(a(iscr+4+2*m+2*k)/emev,7,0)
*d acer.5694
      eyl=xss(igyl+1+j)*emev
*d acer.5700
     1  xss(lgyl+1),xss(lgyl+1+ngyl),eyl/emev,gyl,2)
*d acer.5703
     1  xss(lgyl),xss(lgyl+ngyl),eyl/emev,gyl,2)
*d acer.5709
      xss(next+j)=sigfig(eyl/emev,7,0)
*d acer.5763,5764
  450 xss(next+j)=sigfig(c2h/emev,7,0)
      ee=c2h/emev
*d acer.5780
      e=ee*emev
*d acer.5785
     3     6x,''patching...'')') mt,e/emev,ep/emev
*d acer.5792
     3     6x,''leaving it as is...'')') mt,e/emev,ep/emev
*d acer.5794
      xss(ki+nexd)=sigfig(a(iscr+6+ncyc*(ki-1))/emev,7,0)
*d acer.5796,5797
      if (ki.gt.nd) xss(ki+n+nexd)
     1  =sigfig(a(iscr+7+ncyc*(ki-1))*emev,7,0)
      bzro=xss(nexd+ki+n)/emev
*d acer.5898
      ee=c2h/emev
*d acer.5922,5923
      xss(nexd+ki)=sigfig(a(iscr+8+2*(ki-1))/emev,7,0)
      xss(nexd+npep+ki)=sigfig(a(iscr+9+2*(ki-1))*emev,7,0)
*d acer.6023
      data emev/1.d6/
*d acer.6025
      data emev/1.e6/
*d acer.6098
      xss(nex+i)=sigfig(a(jscr+2*i)/emev,7,0)
*d acer.6214
      xss(nex+ie)=sigfig(c2h/emev,7,0)
*d acer.6250,6251
      xss(nex+5)=sigfig(ef/emev,7,0)
      xss(nex+6)=sigfig(el/emev,7,0)
*d acer.6255
      xss(nex+10)=sigfig(eg/emev,7,0)
*d acer.6455
      xss(nex+10)=sigfig(eg/emev,7,0)
*d acer.6292
      xss(j+nex)=sigfig(a(2*j+2*m+4+iscr)/emev,7,0)
*d acer.6324
  550 xss(nex-1+ie)=sigfig(c2h/emev,7,0)
*d acer.6334,6335
      xss(k+nexd)=sigfig(a(iscr+6+2*k)/emev,7,0)
      xss(k+n+nexd)=a(iscr+7+2*k)*emev
*d acer.6375,6376
      xss(nex+5)=sigfig(a(iscr+6+2*m)/emev,7,0)
      xss(nex+6)=sigfig(a(iscr+4+2*m+2*n)/emev,7,0)
*d acer.6408
  670 xss(nex-1+ie)=sigfig(c2h/emev,7,0)
*d acer.6416
      xss(ki+nexd)=sigfig(a(iscr+6+ncyc*(ki-1))/emev,7,0)
*d acer.6418
      if (ki.gt.nd) xss(ki+n+nexd)=a(iscr+7+ncyc*(ki-1))*emev
*d acer.9353
      data emev/1.d6/
*d acer.9356
      data emev/1.e6/
*d acer.9411,9412
      xss(i+itce)=a(2*(i-1)+iscr)/emev
      xss(i+nee+itce)=a(1+2*(i-1)+iscr)/emev/nmix
*d acer.9425
      xss(itce+i)=a(indx)/emev
*d acer.9443
      xss(itie+i)=a(iscr)/emev
*d acer.9448
      xss(indx+k)=a(iscr+1+k)/emev
*d acer.9759
      data emev/1.d6/
*d acer.9761
      data emev/1.e6/
*d acer.9841
      xss(l)=a(k)/emev
*d acer.10095
      data emev/1.d6/
*d acer.10105
      data emev/1.e6/
*d acer.10175
      xss(eszg+i)=xss(eszg+i)/emev
*/ acer -- 02jul98
*/   fix some places where mixed integer and real variables 
*/   could cause problems.
*d acer.1204
      ln=nint(b(4))
*d acer.1219
      int=nint(a(iscr+5+2*i))
*d acer.2134
      iso=nint(a(iscr+2))
*d acer.13187,13188
      intmu=nint(a7(loci))
      nmu=nint(a7(loci+1))
*d acer.13202,13203
  210 intmu=nint(a7(loci))
      nmu=nint(a7(loci+1))
*/ acer -- 02jul98
*/   fix up the use of float.  for double precision version, 
*/   float will put out a real*4, which is not what we want.
*/   these floats were left over from the old days.
*d acer.2597
      if (n.gt.1) afac=afac*(2*n-1)/n
*d acer.2604,2605
      aa=2*n+1
      aa=aa/2
      a(ixat+n-1)=cf*afac*aa/(n+jfi)
*d acer.2610,2613
      b=(n-ist)*(n-ifi)
      c=ivar*(2*n-ifi)
      d=n+jfi
      cf=cf*asign*b/c
*d acer.2754
      delmu1=(amumid+1)/mfmid
*d acer.2760
      delmu2=(1-amumid)/mfmid
*d acer.2783
      dprob=sum/mcoars
*d acer.2794
      ap=a(ipro+k-1)-dprob*k
*d acer.2803
  330 b=dprob*k
*d acer.3010
      dprob=sum/mcoars
*ident up6
*/ purr -- 26may98
*/   add the capability of producing conditional probabilities 
*/   for heating in the probability tables.  these are best 
*/   if partial heating cross sections are available for
*/   elastic (MT=302), fission (318), and capture (402). The
*/   coding will leave the heating unselfshielded if only the
*/   total heating (mt=301) is present.
*i purr.20
c     *  a conditional probability for heating is added to the table.  *
c     *  if partial heating cross sections for elastic (302), fission  *
c     *  (318), and capture (402) are available from heatr, full       *
c     *  fluctuations will be provided for the total heating.          *
c     *  otherwise, the same value will be provided for each bin.      *
*i purr.59
      external rdheat
*i purr.184
c
c     ***read in the total and partial heating cross sections
      nw=ntemp*nunr*4
      call reserv('heat',nw,iheat,a)
      call rdheat(a,a(iheat),a(ieunr),temp,ntemp,nunr,ihave,matd)
      if (ihave.eq.0) call mess('purr',
     1  'no heating found on pendf',
     2  'ur heating set to zero')
      if (ihave.eq.1) call mess('purr',
     1  'no partial heating xsecs found on pendf',
     2  'ur heating will not selfshield')
*d purr.259
      nc153=(1+6*nbin)*nunr
*d purr.366
      if (12+(1+6*nbin)*nunr.gt.maxscr) call error('purr',
*d purr.381
      a(n+10)=(1+6*nbin)*nunr
*i purr.402
      do 306 j=1,nbin
      n=n+1
      a(n)=0
  306 continue
*i purr.415
      k=iheat+4*nunr*(it-1)+4*(ie-1)-1
      do 319 j=1,nbin
      l=n1+j+5*nbin
      if (ihave.eq.1) then
         if (lssf.eq.1) then
            a(l)=1
         else
            a(l)=a(k+1)/sigu(1,1,1)
         endif
      else if (ihave.eq.2) then
         a(l)=a(l)+(a(k+1)-a(k+2)-a(k+3)-a(k+4))
         h=a(k+2)
         if (lssf.eq.1) then
            h=h*a(n1+j+2*nbin) 
         else if (sigu(2,1,1).ne.0) then
            h=h*a(n1+j+2*nbin)/sigu(2,1,1)
         endif
         a(l)=a(l)+h
         h=a(k+3)
         if (lssf.eq.1) then
            h=h*a(n1+j+3*nbin) 
         else if (sigu(3,1,1).ne.0) then
            h=h*a(n1+j+3*nbin)/sigu(3,1,1)
         endif
         a(l)=a(l)+h
         h=a(k+4)
         if (lssf.eq.1) then
            h=h*a(n1+j+4*nbin) 
         else if (sigu(4,1,1).ne.0) then
            h=h*a(n1+j+4*nbin)/sigu(4,1,1)
         endif
         a(l)=a(l)+h
         if (lssf.eq.1) then
            a(l)=a(l)/a(k+1)
         else
            a(l)=a(l)/sigu(1,1,1)
         endif
      endif
  319 continue
*d purr.419
      nwds=6*nbin*nunr
*i purr.897
c
      subroutine rdheat(a,heat,eunr,temp,ntemp,nunr,ihave,matd)
c     ******************************************************************
c     read the total heating (mt=301) and partial heating cross 
c     sections for elastic (302), fission (318), and capture (402) 
c     from the pendf tape on the unresolved energy grid.
c     ******************************************************************
*if sw
      implicit real*8 (a-h,o-z)
*endif
      common/mainio/nsysi,nsyso,nsyse,ntty
      common/units/nendf,nin,nout
      common/cont/c1h,c2h,l1h,l2h,n1h,n2h,math,mfh,mth,nsh,nsp,nsc
      dimension a(1)
      dimension heat(4,nunr,ntemp)
      dimension eunr(nunr),temp(ntemp)
c
c     ***allocate storage
      call findex('scr',iscr,a)
c
c     ***loop over the temperatures
      call repoz(nin)
      call tpidio(nin,0,0,a(iscr),nb,nw)
      itemp=0
  110 itemp=itemp+1
      call findf(matd,3,1,nin)
      ihave=0
  120 call contio(nin,0,0,a(iscr),nb,nw)
      if (mth.eq.0) go to 150
      ix=0
      if (mth.eq.301) ix=1
      if (mth.eq.302) ix=2
      if (mth.eq.318) ix=3
      if (mth.eq.402) ix=4
      if (ix.gt.0) go to 130
      call tosend(nin,0,0,a(iscr))
      go to 120
  130 e=0.
      call gety1(e,enext,idis,h,nin,a(iscr))
      if (ix.eq.1) ihave=1
      if (ix.gt.1) ihave=2
      do 140 ie=1,nunr
      e=eunr(ie)
      call gety1(e,enext,idis,h,nin,a(iscr))
      heat(ix,ie,itemp)=h
  140 continue
      call tosend(nin,0,0,a(iscr))
      go to 120
  150 continue
      call tomend(nin,0,0,a(iscr))
      if (itemp.lt.ntemp) go to 110
      call repoz(nin)
      return
      end
*ident up7
*/ acer -- 26may98
*/   adjust acer to accept the new probability tables that 
*/   contain a conditional probability for heating.
*d acer.5039
      nurb=(nurb/nure-1)/6
*d acer.5052
      jj=iurd+6+(ie-1)*(1+6*nurb)
*d acer.5062
      xss(ll+5*nurb+ib)=a(jj+5*nurb+ib)/emev
*ident up8
*/ acer -- 20jun98
*/   fix an error introduced by up5.  this problem damages the
*/   list of photon production yield mt values.
*d up5.196
      nex=nex+1
*/ acer -- 6aug98 
*/   fix some problems with the ace files for photon angular
*/   distributions and with the corresponding output listing.
*/   noted by bob little (lanl).
*d acer.6216
      xss(nex+ne+ie)=lc-andp+2
*d acer.7552
      if (k.gt.0) k=k+andp-1
*d acer.7559,7560
      write(nsyso,'(6x,8(4x,a6,a4))') (ek,blank,ii=1,ic)
      write(nsyso,'(5x,1p,8e14.5)') (xss(ii+na),ii=iaa,ib)
*d acer.7568
      write(kk(m),'(1p,e14.5)') xss(ii)
*d acer.7572
      write(nsyso,'(1x,i4,8a14)') j,(kk(ii),ii=1,nkk)
*ident up9
*/ gaspr -- 14aug98
*/   fix an error in the construction of the new file1
*/   dictionary on the pendf tape after the insertion of the 
*/   gas production reactions.  this problem doesn't affect 
*/   the numbers in file3.
*d gaspr.628
      nm=6*nx-j+k+5
      do 400 i=1,nm
*ident up10
*/ reconr -- 15aug98
*/   modify the reconstruction process in both lunion and resxs
*/   to use a tighter tolerance than err in the thermal range.
*/   this means that the .0253 ev cross section will be computed
*/   better, and the thermal parameters generated in broadr can
*/   be computed more accurately.  requested by lubitz (kapl).
*i reconr.1211
      trange=.5d0
*i reconr.1218
      trange=.5
*d reconr.1510
      errn=err
      if (a(ix+i-2).lt.trange) errn=err/5
      test=errn*ym+ssmall
*i reconr.1608
      data trange/.5d0/
*i reconr.1612
      data trange/.5/
*d reconr.1736,1741
      errn=err
      if (a(ix+i-2).lt.trange) errn=err/5
      errm=errmax
      if (a(ix+i-2).lt.trange) errm=errmax/5
      if (dm1.le.errn*sig(2).and.
     1    dm2.le.errn*sig(3).and.
     2    dm3.le.errn*sig(4)) go to 145
      if (dm1.gt.errm*sig(2).or.
     1    dm2.gt.errm*sig(3).or.
     2    dm3.gt.errm*sig(4)) go to 175
*/ reconr -- 20aug98
*/   fix problem with reading in user's energy grid.  this is
*/   a rarely used option.
*d reconr.323
      read(nsysi,*) (a(ienode+i-1),i=1,ngrid)
*ident up11
*/ broadr -- 15aug98
*/   modify the reconstruction tests to use a tighter
*/   tolerance in the thermal range to be consistent with
*/   the changes to reconr in up10.  change the upper limit
*/   for the thermal maxwellian integration to a more
*/   reasonable value.  change the display of the thermal
*/   parameters so that it only appears if the temperature
*/   is close to 293.6 K (0.0253 ev).
*d broadr.28,31
c     * if the temperature is close to 293.6 K (.0253 ev), broadr      *
c     * computes and displays thermal cross sections, maxwellian       *
c     * integrals (one-group thermal cross sections), g-factors,       *
c     * integral ratios (eta, alpha), the k1 integral and the          *
c     * corresponding 1/v-equivalent, and resonance integrals.         *
*d broadr.177
     8  '' starting material temperature ........ '',f10.1,''k''/
*d broadr.179,181
     a  '' max. energy .......................... '',1p,e10.3/
     b  '' errmax for thinning .................. '',1p,e10.3/
     c  '' errint for thinning .................. '',1p,e10.3)')
*d broadr.185,186
     1  '' final temperatures ................... '',1p,e10.3/
     2  (40x,e10.3))') (temp2(i),i=1,ntemp2)
*d broadr.388
         write(nsyse,'(1p,e11.3,'' deg'',54x,0p,f8.1,''s'')')
     1     tempk,time
*d broadr.446
c     ***write out thermal quantities if t is close to 293.6 kelvin.
*d broadr.449
      if (abs(tev-therm).gt.therm/100) go to 2250
*d broadr.479
     1      ''           thermal fission xsec:'',1p,e12.4)') ss
*d broadr.488
     1      ''          thermal fission nubar:'',1p,e12.4)') fnu
*d broadr.493
     1      ''           thermal capture xsec:'',1p,e12.4)') ss
*d broadr.501,502
*d broadr.510,511
*d broadr.516,517
*d broadr.532
      break=1.
*d up2.74
     1    ''                  equivalent k1:'',1p,e12.4)') 
*i broadr.1002
      data trange/.5d0/
*i broadr.1009
      data trange/.5/
*i broadr.1087
      errt=errthn
      if (es(is-1).lt.trange) errt=errt/5
      errm=errmax
      if (es(is-1).lt.trange) errm=errm/5
*d broadr.1099,1100
      if (dy.le.errt*abs(sn(i))+errmin) go to 145
      if (dy.gt.errm*abs(sn(i))+errmin) go to 170
*/ broadr -- 18aug98
*/   modify the doppler-broadening method to extrapolate
*/   below the lowest energy point assuming a 1/v cross section
*/   rather than a constant.  this is consistent with the newer
*/   versions of the sigma1 method from red cullen.  otherwise,
*/   the 1/v cross sections at the lowest energies have a
*/   tendency to sag slightly from the 1/v shape.
*d broadr.1245
c     continue cross sections as 1/v to x=0.
*d broadr.1248
      aa=y
*d broadr.1251
      sbt(i)=sbt(i)-s(i,klow)*e(klow)*(oy**2*h(2)+oy*h(1))
*d broadr.1290
c     extend cross section to x=0.0 as 1/v.
*d broadr.1298
      sbt(i)=sbt(i)-s(i,klow)*e(klow)*(oy**2*h(2)+oy*h(1))
*i broadr.1452
      if (a.eq.alast) go to 120
*ident up12
*/ unresr -- 18aug98
*/   fix error made in rdunf3 when eliminating an arithmetic if.
*d unresr.745,746
      if (mtx(ix).eq.mth) go to 105
      if (mtx(ix).gt.mth) go to 100
*/ unresr -- 18aug98
*/   fix the temperature format to allow for 293.6 degrees.
*d unresr.107,108
     1  '' temperatures ......................... '',1p,e10.3/
     2  (40x,e10.3))') (temp(i),i=1,ntemp)
*d unresr.110,111
     1  '' sigma zero values .................... '',1p,e10.3/
     2  (40x,e10.3))') (sigz(i),i=1,nsigz)
*d unresr.120
     1  '' mat = '',i4,3x,'' temp = '',1p,e10.3,37x,0p,f8.1,''s'')')
*/ unresr -- 24aug98
*/   fix the wide parameter to match reconr.
*d unresr.370
      data wide/3.01d0/
*d unresr.376
      data wide/3.01/
*ident up13
*/ purr -- 18aug98
*/   increase the storage space to make room for endf/b-vi pu-239.
*/   patch purr for elemental evaluations, and allow space for
*/   more spin sequences to handle the case of endf/b-vi mo-nat.
*d purr.51
      common/pustore/a(70000)
*d purr.65
      data namax/70000/, nidmax/24/, ipr/1/
*d purr.910,912
      common/sigcon/e,t,cth(40),csz(40),cc2p(40),cs2p(40),
     1   cgn(40),cgg(40),cgf(40),cgx(40),cgt(40),dbar(40),
     2   spot,dbarin,sigi(4),ndfn(40),ndff(40),ndfx(40),nseq0
*d purr.1050
      if (nseq0.gt.40) call error('unresx','too many sequences',' ')
      csz(nseq0)=abn*ab*gj
*d purr.1070
      temp=abn*pi*ab*gj*gnx/(2*dx)
*d purr.1278,1280
      common/sigcon/e,t,cth(40),csz(40),cc2p(40),cs2p(40),
     1   cgn(40),cgg(40),cgf(40),cgx(40),cgt(40),dbar(40),
     2   spot,dbarin,sigi(4),ndfn(40),ndff(40),ndfx(40),nseqz
*d purr.1381,1383
      common/sigcon/e,t,cth(40),csz(40),cc2p(40),cs2p(40),
     1   cgn(40),cgg(40),cgf(40),cgx(40),cgt(40),dbar(40),
     2   spot,dbarin,sigi(4),ndfn(40),ndff(40),ndfx(40),nseqz
*d purr.1487,1489
      common/sigcon/e,t,cth(40),csz(40),cc2p(40),cs2p(40),
     1   cgn(40),cgg(40),cgf(40),cgx(40),cgt(40),dbar(40),
     2   spot,dbarin,sigi(4),ndfn(40),ndff(40),ndfx(40),nseq0
*/ purr -- 18aug98
*/   change the temperature format to allow for 293.6.
*d purr.114
     1   '' temperatures ......................... '',1p,e10.3)')
*d purr.116
      if (ntemp.gt.1) write(nsyso,'(40x,1p,e10.3)') (temp(i),i=2,ntemp)
*d purr.120
      if (nsigz.gt.1) write(nsyso,'(40x,1p,e10.3)') (sigz(i),i=2,nsigz)
*d purr.2020
      write(nsyso,'(3x,1p,2e10.3,5e12.4)')
*d purr.2071
      write(nsyso,'(1x,a,1p,e10.3,10e11.3/(14x,10e11.3))')
*d purr.2110
      write(nsyso,'(3x,1p,2e10.3,5e12.4)')
*/ purr -- 24aug98
*/   fix the wide parameter to match reconr.
*d purr.483
      data wide/3.01d0/
*d purr.489
      data wide/3.01/
*ident up14
*/ acer -- 24aug98
*/   don't divide the heating by 1e6 ev when it is really
*/   a factor.  this results in unresolved heating values too 
*/   small by a factor of 1e6 for two endf/b-vi materials.  
*/   found by little (los alamos).
*d up7.10
      if (lssf.eq.0) then
         xss(ll+5*nurb+ib)=a(jj+5*nurb+ib)/emev
      else
         xss(ll+5*nurb+ib)=a(jj+5*nurb+ib)
      endif
*/ acer -- 24aug98
*/   fix some problems in the length of strings in common
*d acer.11938
      character hz*13,hd*10,hk*70,hm*10
*d acer.12141
      character hz*13,hd*10,hk*70,hm*10
*d acer.12244
      character hz*13,hd*10,hk*70,hm*10
*d acer.12334
      character hz*13,hd*10,hk*70,hm*10
*d acer.12706
      character hz*13,hd*10,hk*70,hm*10
*/ acer -- 22sep98
*/   fix the initialization of e for gety2 for subsections after
*/   the first.  this was a byproduct of taking constants out of 
*/   subroutine calls.  it causes problems for materials with
*/   anisotropic photons, such as n and o from endf/b-vi.
*i acer.3595
      e=0
*/ acer -- 22sep98
*/   fix the logic for identifying which photon an angular
*/   distribution corresponds to.  the old scheme only used
*/   energy, and if two different reactions produced photons
*/   with the same energy, the angular distribution would
*/   always be indexed to the first photon in the list with
*/   that energy.
*d acer.6197
      eps=1.e-6
      mmm=nint(xss(mtrp+j-1))
      mmm=mmm/10000
      if (mmm.ne.mtd) go to 340
      if (abs(eg-a(iphot+5*(j-1))).gt.eps*eg) go to 340
*d acer.6207
      eps=1.e-6
      mmm=nint(xss(mtrp+j-1))
      mmm=mmm/10000
      if (mmm.ne.mtd) go to 370
      if (abs(eg-a(iphot+5*(j-1))).gt.eps*eg) go to 370
*/ acer -- 22sep98
*/   allow for more angles in the contour plots to handle the case
*/   of endf/b-vi oxygen.
*d acer.12247
      dimension loce(1200)
*d acer.12271
      if (ne.gt.1200) call error('aplof4',
*ident up15
*/ covr -- 23sep98
*/   fix some problems with the covariance plots that arose when
*/   we changed to using quotes to delimit character strings.
*/   the reaction names containing "prime" do not print out 
*/   right.  we remove the prime, counting on the level subscript
*/   do tell the user that it is an inelastic reaction.  also,
*/   change the blank labels to ".".  this is a flag to viewr
*/   to use a blank label, but to go ahead with the numbering
*/   of the axis.  these problems are visible in the file plot05
*/   from the original njoy97 test set.  noted by muir (iaea).
*d covr.966
      write(nplot,'(a,''.'',a,''/'')') qu,qu
*d covr.968
      write(nplot,'(a,''.'',a,''/'')') qu,qu
*d covr.1022
      write(nplot,'(a,''.'',a,''/'')') qu,qu
*d covr.1024
      write(nplot,'(a,''.'',a,''/'')') qu,qu
*d covr.1557
      data nmeb1/'(n,n'/, nmee1/'(delayed'/, nmef1/'(prompt'/
*d covr.1573
      inamel=4
*ident up16
*/ viewr -- 23sep98
*/   when using quotes as character delimiters, it is hard to
*/   tell the difference between a blank label and a "default"
*/   label.  we change the input so that a label of '.' results
*/   in a blank label on the plot with numbering on the axis
*/   ticks.  a label of ' ' results in a "default" label, that
*/   is, no text label, and no axis numbering.  these features
*/   are used for the covariance plots from covr (see up15).
*i viewr.2438
      if (label.eq.'.') return
*ident up17
*/ thermr -- 23sep98
*/   modify output format to allow for 4-digit temperatures.
*d thermr.2564
     1  '' wrote thermal data for temp='',1p,e10.3,30x,0p,f8.1,''s'')')
*d thermr.2567
     1  '' wrote thermal data for temp='',1p,e10.3,30x,0p,f8.1,''s'')')
*ident up18
*/ acer -- 24sep98
*/   move some of the legend blocks so that they don't interfere
*/   with the curves as much.  this only works for curves with
*/   very predictable locations.
*d acer.11134
      xtag=5*xmin
      ytag=0.6*log10(ymin)+0.4*log10(ymax)
      ytag=10.**ytag
      write(nout,'(''4 0 2 1'',2e12.4,''/'')') xtag,ytag
*d acer.11588,11610
  270 xmin=0.
      xmax=xss(esz-1+nes)
      ymin=0.
      ymax=xss(esz+nes-1+nes)
      if (gpd.gt.0) then
         if (xss(gpd-1+nes).gt.ymax) ymax=xss(gpd-1+nes)
      endif
      ymax=2*ymax
*d acer.11622
      xtag=0.45*xmin+0.55*xmax
      ytag=0.1*ymin+0.9*ymax
      write(nout,'(''1 0 2 1'',2e12.4,''/'')') xtag,ytag
*d acer.11987
      xtag=0.3*xmin+0.7*xmax
      ytag=0.1*ymin+0.9*ymax
      write(nout,'(''1 0 2 1'',2e12.4,''/'')') xtag,ytag
*/ acer -- 24sep98
*/   fix the expanded resonance plots to be more clever about
*/   the energy ranges plotted.  check the number of points in
*/   each decade.  if is is more than some cutoff, make a plot
*/   showing that decade.
*d acer.11216,11230
      maxii=120
      e1=xss(esz)
      ii1=1
  169 ii2=ii1
      e2=10*e1
 1169 ii2=ii2+1
      if (ii2.ge.nes) go to 2169
      if (xss(esz-1+ii2).ge.e2) go to 2169
      go to 1169
 2169 if (ii2-ii1.gt.maxii) go to 3169 
      if (ii2.ge.nes) go to 1174
      e1=e2
      ii1=ii2
      go to 169
 3169 e2=xss(esz-1+ii2)
*i acer.11279
         ii1=ii2+1
*i acer.11281
 1174 continue
*d acer.11284,11297
      maxii=100
      e1=xss(esz)
      ii1=1
  168 ii2=ii1
      e2=10*e1
 1168 ii2=ii2+1
      if (ii2.ge.nes) go to 2168
      if (xss(esz-1+ii2).ge.e2) go to 2168
      go to 1168
 2168 if (ii2-ii1.gt.maxii) go to 3168 
      if (ii2.ge.nes) go to 183
      e1=e2
      ii1=ii2
      go to 168
 3168 e2=xss(esz-1+ii2)
*i acer.11385
         ii1=ii2+1
*/ acer -- 24sep98
*/   separate out the gas-production reactions onto their own page
*i acer.12060
      if (mt.lt.203.and.nint(xss(mtr+i)).ge.203) nlev=5
*i acer.12120
      if (mt.lt.203.and.nint(xss(mtr+i)).ge.203) nlev=5
*d acer.12125
      go to 210
*/ acer -- 24sep98
*/   fix another variable that is not defined except for type2
*/   format.  see up5 for the other one.
*i up5.6
      nern=0
*/ acer -- 24sep98 
*/   fix a numerical problem that shows up for endf/b-vi h-2.
*d acer.5976,5977
      if (xx.gt.test) then
         xx=1
         pn=0
      else
         pn=sqrt(xx)*(1-xx)**(3*npsx/2-4)
      endif
*/ acer -- 24sep98
*/   fix the temperature print format to allow for 293.6k.
*d acer.252
     2  '' temperature .......................... '',1p,e10.3/
*d acer.309
     2  '' temperature .......................... '',1p,e10.3/
*ident up19
*/ thermr -- 24sep98
*/   fix a numerical problem, and allow for larger data sets
*/   (e.g., h in zrh from endf/b-vi rel. 3).
*d thermr.94
      common/tstore/a(80000)
*d thermr.127
      namax=80000
*d thermr.130
      nwscr=40000
*d thermr.2057
      data shade/.99999999d0/
*d thermr.2063
      data shade/.99999999/
*ident up20
*/ njoy -- 25sep98
*/   fix some typing errors in njoy
*d njoy.497
      if (hdate(1:1).eq.' ') hdate(1:1)='0'
*d njoy.510
      if (htime(1:1).eq.' ') htime(1:1)='0'
*ident up21
*/ errorr -- 25sep98
*/   permit user energies above 20 mev in errorr when iread=0
*/   and correct a related minor code bug
*d errorr.2986,2996
      etop=a(iun+nunion)
      etop=sigfig(etop,ndig,0)
      if (a(iela+ng).ge.etop) go to 155
      a(iela+ng)=etop
      write(strng,'(''since iread='',i2)') iread
      call mess('colaps',
     1 'raising top of ngout structure to top of union structure',' ')
  155 continue
      if (a(iela).gt.a(iun)) call error('colaps',
     1  'ngout group structure does not span union grid.',' ')
*d errorr.3147,3152
*i errorr.3164
      if (ign.ne.1.and.ign.ne.19.and.egn(ngnp1).lt.a(ie+neni-1)) then
         call mess('uniong',
     1   'raising top built-in multigroup boundary to top of file 33',
     2   ' ')
         egn(ngnp1)=a(ie+neni-1)
      endif
*d errorr.3181,3182
c
c     ***finished with endf energies
      in=i
      go to 170
c
c     ***treat endf energies below first group boundary
*i errorr.3200
      do 125 i=2,k
      if (a(iun+i-1).le.a(iun+i-2)) then
         write(strng,'(1p,e12.4,'' le '',1p,e12.4)')
     1      a(iun+i-1),a(iun+i-2)
         call error('uniong','union energies out of order',strng)
      endif
  125 continue
*d errorr.3709,3717
      call uniong(nendf,a)
*ident up22
*/ covr -- 25sep98
*/   fix two additional places where a space needs to be a dot
*/   in order to get axis values labeled.  See up15 for others.
*d covr.1258
      write(nplot,'(a,''.'',a,''/'')') qu,qu
*d covr.1328
      write(nplot,'(a,''.'',a,''/'')') qu,qu
*/ covr -- 25sep98
*/   restore a deleted line and fix a typing error in covr
*i covr.1750
      ndig=nvf-6
*d covr.1838
      test=1.e-10*abs(xctest)
*ident up23
*/ moder -- 21jun98 
*/   add new ltt=3 format.  this format modification was introduced 
*/   to help with adding high-energy data to endf evaluations.  
*/   it allows for the use of legendre polynomials for elastic 
*/   scattering at lower incident energies, with a change to
*/   tabulated distributions at higher energies.
*d moder.595
      if (lvt.eq.0.and.l1h.eq.1) go to 170
*d moder.618
      go to 170
*i moder.632
  170 if (ltt.lt.3) go to 200
      call tab2io(nin,nout,nscr,a,nb,nw)
      ne=n2h
      do 190 ie=1,ne
      call tab1io(nin,nout,nscr,a,nb,nw)
      if (nb.eq.0) go to 190
  180 call moreio(nin,nout,nscr,a,nb,nw)
      if (nb.gt.0) go to 180
  190 continue
*ident up24
*/ heatr -- 21jun98
*/   add new ltt=3 format to heatr. 
*d heatr.3517
      save elo,ehi,nlo,nhi,flo,fhi,ltt,ltt3,lttn
*i heatr.3532
      ltt3=ltt
      if (ltt.eq.3) then
         ltt=1
         lttn=1
      endif
*d heatr.3592,3593
      if (nne.eq.ne.and.ltt3.eq.3.and.lttn.eq.1) then
         call tab2io(nin,0,0,c(1),nb,nwc)
         ne=nint(c(6))
         nne=0
         ir=1
         ltt=2
         lttn=2
      else if (nne.eq.ne) then
         call error('hgtfle','desired energy above highest given.',' ')
      endif
*/ heatr -- 21jun98
*/   make more room for lo=2 gammas and provide a protective
*/   error message.
*i heatr.3827
      lmax=100
*d heatr.3866,3868
      call reserv('eg',lmax,ieg,a)
      call reserv('es',lmax,ies,a)
      call reserv('y',lmax,iy,a)
*i heatr.3930
      if (l.gt.lmax) call error('hconvr',
     1  'too many lo=2 gammas',' ')
*ident up25
*/ groupr -- 21jun98
*/   add new ltt=3 format to groupr. 
*i groupr.6142
      save ltt3,lttn
*i groupr.6165
      ltt3=ltt
      if (ltt.eq.3) then
         ltt=1
         lttn=1
      endif     
*d groupr.6236,6237
      if (nne.eq.ne.and.ltt3.eq.3.and.lttn.eq.1) then
         call tab2io(nin,0,0,c(1),nb,nwc)
         ne=nint(c(6))
         nne=0
         ir=1
         ltt=2
         lttn=2
      else if (nne.eq.ne) then
         call error('getfle','desired energy above highest given.',' ')
      endif
*/ groupr -- 21jun98
*/   make more room for lo=2 gammas and provide a protective 
*/   error message.
*i groupr.7243
      lmax=100
*d groupr.7420,7422
      call reserv('eg',lmax,ieg,a)
      call reserv('es',lmax,ies,a)
      call reserv('y',lmax,iy,a)
*i groupr.7483
      if (l.gt.lmax) call error('hconvr',
     1  'too many lo=2 gammas',' ')
*ident up26
*/ njoy -- 25sep98
*/   fix the typing errors in njoy again.
*d up20.5
      if (hdate(i:i).eq.' ') hdate(i:i)='0'
*d up20.7
      if (htime(i:i).eq.' ') htime(i:i)='0'
*ident up27
*/ acer -- 11oct98
*/   add the word ptable to the xsdir cards when unresolved
*/   probability tables are present
*d acer.7873
     1 gpd,mtrp,lsigp,sigp,landp,andp,ldlwp,dlwp,yp,fis,end,
     2 iurpt,jxsd(9)
*d acer.8053,8055
      if (iurpt.gt.0) then
         write(ndir,
     1    '(a10,f12.6,'' filename route'',i2,i4,i7,2i6,1p,e10.3,
     2    '' ptable'')') hz,aw0,itype,irec1,len2,lrec,nern,tz
      else
          write(ndir,
     1     '(a10,f12.6,'' filename route'',i2,i4,i7,2i6,1p,e10.3)')
     2     hz,aw0,itype,irec1,len2,lrec,nern,tz
      endif
*/   add more storage space to take care of the large probability
*/   tables in endf/b-vi pu-239.
*d acer.184
      common/astore/a(35000)
*d acer.193
      data namax/35000/, nidmax/27/
*d acer.402
      common/astore/a(35000)
*d acer.4570
      common/astore/a(35000)
*d acer.5151
      common/astore/a(35000)
*d acer.5499
      common/astore/a(35000)
*d acer.6019
      common/astore/a(35000)
*d acer.6625
      common/astore/a(35000)
*d acer.7346
      common/astore/a(35000)
*d acer.8999
      common/astore/a(35000)
*d acer.9347
      common/astore/a(35000)
*d acer.9747
      common/astore/a(35000)
*d acer.10072
      common/astore/a(35000)
*ident up28
*/ purr -- 15oct98
*/   fix some tests that are supposed to protect the sqrt
*/   function against small negative arguments.
*d purr.1977
      if (argt.lt.0.) argt=0
*d purr.1980
      if (arge.lt.0.) arge=0
*d purr.1983
      if (argf.lt.0.) argf=0
*d purr.1986
      if (argc.lt.0.) argc=0
*ident up29
*/ covr -- 5nov98
*/   correct typos in 2 reaction names and add 7 new names.
*/   contributed by muir (iaea).
*d covr.1542,1544
      character*8 lnamel(2),hira(32),blank,nmea1,nmeb1,lnamer,
     1  nmed1,nmee1,nmef1,nmeh1
      dimension ira1(32),ira2(32)
*d covr.1547
     2  207,780,781,452,455,456,25,24,32,33,41,112,115,116/
*d covr.1549
     1  5,4,7,6,3,3,2,5,3,1,6,4,2,2,2,6,4,3,5,5,5,5,5,5,6,6,3,3,4,5,3,3/
*d covr.1552,1554
     2  'n]a<)','np)',')','cont.)',']g<)','p)','d)','t)',
     3  '<h>e3)',']a<)','2p)','<h>e)',']a<0)',']a<1)',' ]n<)',' ]n<)',
     4  ' ]n<)','3n]a<)','2n]a<)','nd)','nt)','2np)','p]a<)','pd)',
     5  'pt)'/
*d covr.1556
     1  nmed1/'(total'/
*d covr.1563
*d covr.1567
*d covr.1585,1587
*/ covr -- 5nov98
*/   correct error in color plot option (muir, iaea).
*d covr.255
      if (nplot.ne.0) write(nplot,'(''1 2 .22'',i3,''/'')') icolor
*/ covr -- 12nov98
*/   correct diagnostic for empty plots and fix shading
*/   for color plots, b-w plots, and ndiv greater than 1.
*d covr.1223,1230
      call patlev(jpat,ilevel)
      if (jpat.ne.0) ishade=ishade+1
*d covr.1260,1318
c     ***positive part
      tlow=0.
      do 420 ilevel=1,nlev
      call patlev(jpat,ilevel)
c     ***switch orientation for positive part of legend on b-w plots
      if(icolor.eq.0.and.jpat.ne.0) jpat=jpat+10
      write(nplot,'(''0 0 0 0 0'',i3,''/'')') jpat
      write(nplot,'(''0'')')
      thi=a(ixlev-1+ilevel)
      write(nplot,'(''0. '',f6.3,''/'')') tlow
      write(nplot,'(''1. '',f6.3,''/'')') tlow
      write(nplot,'(''1. '',f6.3,''/'')') thi
      write(nplot,'(''0. '',f6.3,''/'')') thi
      write(nplot,'(''0. '',f6.3,''/'')') tlow
      write(nplot,'(''/'')')
      if (ilevel.lt.nlev) write(nplot,'(''2/'')')
      if (ilevel.lt.nlev) write(nplot,'(''/'')')
      tlow=thi
  420 continue
*d covr.1330,1388
c     ***negative part
      tlow=0.
      do 440 iloop=1,nlev
      ilevel=-iloop
      call patlev(jpat,ilevel)
      write(nplot,'(''0 0 0 0 0'',i3,''/'')') jpat
      write(nplot,'(''0'')')
      thi=-a(ixlev-1+iloop)
      write(nplot,'(''0. '',f6.3,''/'')') tlow
      write(nplot,'(''1. '',f6.3,''/'')') tlow
      write(nplot,'(''1. '',f6.3,''/'')') thi
      write(nplot,'(''0. '',f6.3,''/'')') thi
      write(nplot,'(''0. '',f6.3,''/'')') tlow
      write(nplot,'(''/'')')
      if (iloop.lt.nlev) write(nplot,'(''2/'')')
      if (iloop.lt.nlev) write(nplot,'(''/'')')
      tlow=thi
  440 continue
*i covr.1432
      subroutine patlev(jpat,ilevel)
c     ******************************************************************
c     convert the correlation level to a color or b-w shading pattern
c     following the conventions of viewr
c     ******************************************************************
*if sw
      implicit real*8 (a-h,o-z)
*endif
      common/cov0/icolor
      common/cov2/xcycle,nlev,ndiv,ixmin,ixmax,
     1  izero,irelco,ismall,ishade
c
      jpat=0
      if (icolor.eq.0) then
         if (ilevel.gt.1) jpat=20-(nlev-ilevel)*2/ndiv
         if (ilevel.lt.-1) jpat=40-(nlev+ilevel)*2/ndiv
      else
         if (ilevel.gt.1) jpat=50-(nlev-ilevel)*2/ndiv
         if (ilevel.lt.-1) jpat=60-(nlev+ilevel)*2/ndiv
      endif
      return
      end
*ident up30
*/ njoy -- 5nov98
*/   add a comment-card capability for njoy input decks.  any
*/   lines starting with "-- " will be ignored (from muir, iaea).
*i njoy.142
c
c     ***comment card
c     ***treat comment card as a null module
      else if (module.eq.'--') then
         go to 100
*/ njoy -- 2nov98
*/   increase size of field in storage-error diagnostic (muir, iaea)
*d njoy.2517
     1   '(''need '',i7,'' more words for id '',a4,''.'')') need,id
*ident up31
*/ reconr -- 5nov98
*/   fix typographical error (noted by p.deleege, tu/delft)
*d up1.73,77
*if sw
         phi=dcmplx(phro,phri)
*else
         phi=cmplx(phro,phri)
*endif
*ident up32
*/ unresr -- 5nov98
*/   fix misplaced update from up12.  Noticed by
*/   pelloni (psi) and deleege (tu/delft).
*d up12.22
     4  63000.,80000.,100000./
*d unresr.377
      data wide/3.01/
*ident up33
*/ purr -- 5nov98
*/   fix misplaced update from up13.  Noticed by
*/   pelloni (psi) and deleege (tu/delft).
*d up13.50
     4  63000.,80000.,100000./
*d purr.490
      data wide/3.01/
*ident up34
*/ reconr -- 5nov98 
*    allow for more unresolved-range grid energies to handle
*/   the case of fe-56 from eff-3.0 (fendl-2.0).  it corrects
*/   a fatal error when the ur range exceeds 100 kev.  this 
*/   patch was provided by trkov (ijs/slovenia).
*d reconr.491  
      dimension egridu(58)
*d reconr.494  
      data ngridu/58/
*d reconr.501  
     5  63000.d0,80000.d0,100000.d0,125000.d0,150000.d0,200000.d0,
     6  320000.d0,400000.d0,500000.d0,630000.d0,800000.d0,1.d6,
     7  1.25d6,1.5d6,2.0d6,3.2d6,4.0d6,5.0d6,6.3d6,8.0d6/
*d reconr.512  
     5  63000.,80000.,100000.,125000.,150000.,200000.,
     6  320000.,400000.,500000.,630000.,800000.,1.e6,
     7  1.25e6,1.5e6,2.0e6,3.2e6,4.0e6,5.0e6,6.3e6,8.0e6/
*ident up35
*/ unresr -- 5nov98
*/   corresponding changes to unresr (see up35)
*d unresr.358  
      dimension egridu(58)
*d unresr.362  
      data ngridu/58/
*d unresr.369  
     5  63000.d0,80000.d0,100000.d0,125000.d0,150000.d0,200000.d0,
     6  320000.d0,400000.d0,500000.d0,630000.d0,800000.d0,1.d6,
     7  1.25d6,1.5d6,2.0d6,3.2d6,4.0d6,5.0d6,6.3d6,8.0d6/
*d up32.6  
     5  63000.,80000.,100000.,125000.,150000.,200000.,
     6  320000.,400000.,500000.,630000.,800000.,1.e6,
     7  1.25e6,1.5e6,2.0e6,3.2e6,4.0e6,5.0d6,6.3e6,8.0e6/
*ident up36
*/ purr -- 5nov98
*/   corresponding changes to purr (see up35)
*d purr.472  
      dimension egridu(58)
*d purr.475  
      data ngridu/58/
*d purr.482  
     5  63000.d0,80000.d0,100000.d0,125000.d0,150000.d0,200000.d0,
     6  320000.d0,400000.d0,500000.d0,630000.d0,800000.d0,1.d6,
     7  1.25d6,1.5d6,2.0d6,3.2d6,4.0d6,5.0d6,6.3d6,8.0d6/
*d up33.6  
     5  63000.,80000.,100000.,125000.,150000.,200000.,
     6  320000.,400000.,500000.,630000.,800000.,1.e6,
     7  1.25e6,1.5e6,2.0e6,3.2e6,4.0e6,5.0d6,6.3e6,8.0e6/
*ident up37
*/ groupr -- 5nov98
*/   the addition of the logic to automatically determine the
*/   version of endf being used damaged the writing of the id
*/   record on the gendf file.  in addition, the count for the
*/   title field was undefined from subroutine ruinb.  this
*/   patch was modified from one provided by trkov (ijs/sovenia). 
*d groupr.292,298
*d groupr.316,317
      math=1                                                             
      mfh=0
      mth=0
      nsh=0
      itend=0
      nwds=17
      text=' '
      read(text,'(16a4,a2)') (a(iscr+i-1),i=1,nwds)
      if (ngout1.eq.0) call tpidio(0,ngout2,0,a(iscr),nb,nwds)
*d groupr.749  
      common/gout1/ngout1,ngout2,nt
*/ groupr -- 5nov98
*/   fix a problem with the input of a weighting function.
*/   it was introduced by the changeover from infree to read*.
*/   discovered by trkov (ijs/slovenia).
*d groupr.2140,2142
*d groupr.2144,2148
      nr=nint(a(iwght+4))
      np=nint(a(iwght+5))
      iw=6+2*nr+2*np
*ident up38
*/ acer -- 5nov98 
*/   previous versions of acer converted cm legendre polynomials
*/   in file6/law1 into kalbach distributions, but this does not
*/   always provide a good representation of the evaluator's
*/   intent.  this patch converts these data into a laboratory
*/   file6/law7 form instead.  it is needed for fe-56 from
*/   eff-3.1.  the patch was contributed by trkov (ijs/slovenia).
*/   mcnp-4 will provide direct sampling from the original cm
*/   distribution, and this patch will no longer be effective.
*d acer.2063 
c     set flag to convert file-6 Law-7 data representation to Law-7
      if (lf.eq.1) new6=1
*d acer.2093 
c     ***convert mf6,law1 lab or cm distributions to law7
*d acer.3299,3300
c     convert a section of file 6 using legendre or tabulated
c     angular distributions into law 7 format with 33 angles.
*d acer.3317 
      ncos=33
*i acer.3321 
      awr=c2h
      lct=l2h
c     force laboratory coordinate system
      a(4)=1
*i acer.3352 
      ein=a(l1+1)
*d acer.3380 
*d acer.3398 
      csn=acos
      elb=ep
      drv=1
      if (lct.eq.1) go to 296
      if (ep.le.0) go to 296
      aw1=awr+1
      clb=csn
c     minimum lab cosine (= zero particle energy in C.M.)
      cmn=-1
      qq=1-aw1*aw1*ep/ein
      if (qq.gt.0) cmn=sqrt(qq)
      if (clb.lt.cmn) then
         clb=cmn
c        zero distribution when far outside valid cosine range
         if(imu.lt.nmu.and.amu(imu+1).le.cmn) drv=0
      endif
c     outgoing particle energy in the laboratory system
      qq=ep-ein*(1-clb*clb)/(aw1*aw1)
      if (qq.lt.0) qq=0
      elb=clb*sqrt(ein)/aw1+sqrt(qq)
      elb=elb*elb
c     calculate corresponding cosine in the cm system
      csn=clb*sqrt(elb/ep)-sqrt(ein/ep)/aw1
c     check the limits
      qq=-1
      if (csn.lt.qq) csn=qq
      qq=1
      if (csn.gt.qq) csn=qq
  296 if (lang.gt.10) go to 301
c     calculate the probability from Legendre polynomials
      call legndr(csn,p,nl)
*i acer.3403 
c     calculate the probability from pointwise representation
*d acer.3412,3413
  305 a(j)=elb
      a(j+1)=fmu*drv
*ident up39
*/ acer -- 5nov98 
*/   in up14, the length of the hz string in common/mis1 was
*/   set to 13 in preparation for the new mcnpx format, but 
*/   it also needs to be set in other routines.  this was
*/   being done in the mcnp4/mcnpx jumbo acer patch (to be
*/   released later), but these changes are being moved here
*/   to keep the compilers happy while waiting for the
*/   jumbo patch.  recommended by trkov (ijs/slovenia).
*d acer.187  
      character hz*13,hd*10,hk*70,hm*10
*d acer.4553 
      character hz*13,hd*10,hk*70,hm*10
*d acer.6609 
      character hz*13,hd*10,hk*70,hm*10
*d acer.7869 
      character hz*13,hd*10,hk*70,hm*10
*d acer.9332 
      character hz*13,hd*10,hk*70,hm*10
*d acer.9467 
      character hz*13,hd*10,hk*70,hm*10
*d acer.9637 
      character hz*13,hd*10,hk*70,hm*10
*d acer.9742 
      character hz*13,hd*10,hk*70,hm*10
*d acer.9891 
      character hz*13,hd*10,hk*70,hm*10
*d acer.9971 
      character hz*13,hd*10,hk*70,hm*10
*d acer.10066
      character hz*13,hd*10,hk*70,hm*10
*d acer.10384
      character hz*13,hd*10,hk*70,hm*10
*d acer.10476
      character hz*13,hd*10,hk*70,hm*10
*d acer.10558
      character hz*13,hd*10,hk*70,hm*10
*d acer.10761
      character hz*13,hd*10,hk*70,hm*10
*d acer.11081
      character hz*13,hd*10,hk*70,hm*10
*ident up40
*/ wimsr -- 5nov98
*/   the following wimsr updates were provided by a.trkov,
*/   (ijs/slovenia).  define nfid from rdfid, and make
*/   sure that some flags are consistent.
*i wimsr.168  
      nfid=rdfid+0.001
*i wimsr.181
      if (ifprod.gt.0) then
         ifprod=1
         if (ires.gt.0) ifprod=2
      endif
*d wimsr.1861,1862
      jfis=ifis
      if (ifprod.eq.1) jfis=-1
      if (ifprod.eq.2) jfis=-2
*/ wimsr -- 5nov98
*/   extend absorption and multiple neutron emission reactions.
*d wimsr.929  
      if (mth.eq.16.or.mth.eq.24) go to 236
      if (mth.eq.17.or.mth.eq.25) go to 237
*i wimsr.982  
c
c     ***n3n, n3n-alpha cross section
  237 loc=in2n+jg-1
      loca=l+lz+nl*nz
      a(loc)=a(loc)+2.*a(loca)
      go to 360
*/ wimsr -- 5nov98
*/   improve input instructions for burnup data.
*/   add delimiter records when burnup data are present.
*d wimsr.79   
c     *   the following cards 5 and 6 are for iburn gt 0 only          *
*i wimsr.81
c     *            for burnable materials ntis=2                       *
c     *            for fissile materials ntis>2 when fission product   *
c     *            yields are given.
*d wimsr.84
c     * card  6a                                                       *
c     *    identa  ident of capture product isotope                    *
c     *    yield   yield of product identa from capture                *
c     *                                                                *
c     * card  6b                                                       *
c     *    identa  ident of decay product isotope (zero if stable)     *
c     *    lambda  decay constant (s-1)                                *
c     *                                                                *
c     * card  6c  (repeated ntis-2 times, if necessary)
*i wimsr.1875 
      write(nout,'(2i15)') 999999999,3
*i wimsr.1879 
      write(nout,'(2i15)') 999999999,4
*d wimsr.1875 
      if (jcc.ge.8) ifisp(4)=jfis
*ident up41
*/ acer -- 05nov98
*/   fix a problem with the dimensions of two arrays that was
*/   discovered by broeders (fzk).
*d acer.5153
      dimension xs(20),ys(20),xxs(200),yys(200)
*d acer.5155
      data ismax/20/, jsmax/200/
*d acer.5448
      if (js.gt.jsmax) js=jsmax
*ident up42
*/ plotr -- 17nov98
*/   fix some problems noted in plotr by m.mattes (u.stuttgart).
*/   the values of xtag and ytag are not being read in.  factx and
*/   and facty are being applied to the x and y values in plotr,
*/   so they need to be reset to unity before going on to viewr.
*/   this is necessary to prevent the factors from being applied
*/   to the coordinates twice.
*d plotr.452
      read(nsysi,*) itype,jtype,igrid,ileg,xtag,ytag
*d plotr.1218,1219
      factx1=1
      facty1=1
      write(nplt,'(i4,i8,7f7.2,''/ 3d plot'')') iplot,iwcol,
     1   factx1,facty1,xll,yll,ww,wh,wr
*d plotr.1357,1358
      factx=1
      facty=1
      write(nplt,'(i4,i8,7f7.2,''/ 3d plot'')') iplot,iwcol,
     1   factx1,facty1,xll,yll,ww,wh,wr
*d plotr.1693,1694
      factx1=1
      facty1=1
      write(nplt,'(i4,i8,7f7.2,''/ 3d plot'')') iplot,iwcol,
     1   factx1,facty1,xll,yll,ww,wh,wr
*d plotr.1762,1763
  610 factx1=1
      facty1=1
      write(nplt,'(i4,i8,7f7.2,''/ 2d plot'')') iplot,iwcol,
     1   factx1,facty1,xll,yll,ww,wh,wr
*ident up43
*/ unresr -- 20nov98
*/   the changes to the unresolved energy grid in up34 were
*/   not quite correct for the sw machines.  trkov.
*d unresr.365,368
     1  63.d0,80.d0,100.d0,125.d0,150.d0,200.d0,250.d0,320.d0,400.d0,
     2  500.d0,630.d0,800.d0,1000.d0,1250.d0,1500.d0,2000.d0,2500.d0,
     3  3200.d0,4000.d0,5000.d0,6300.d0,8000.d0,10000.d0,12500.d0,
     4  15000.d0,20000.d0,25000.d0,32000.d0,40000.d0,50000.d0,
*ident up44
*/ purr -- 20nov98
*/   the changes to the unresolved energy grid in up35 were
*/   not quite correct for the sw machines.  trkov.
*d purr.478,481
     1  63.d0,80.d0,100.d0,125.d0,150.d0,200.d0,250.d0,320.d0,400.d0,
     2  500.d0,630.d0,800.d0,1000.d0,1250.d0,1500.d0,2000.d0,2500.d0,
     3  3200.d0,4000.d0,5000.d0,6300.d0,8000.d0,10000.d0,12500.d0,
     4  15000.d0,20000.d0,25000.d0,32000.d0,40000.d0,50000.d0,
*ident up45
*/ covr -- 20nov98
*/   one change was left out of muir's covr patches in up29.
*/   it is needed to get correct labels for plots of the
*/   high threshold reactions.
*d covr.1594
      do 300 j=1,32
*ident up46
*/ njoy -- 4dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision, and that any conversion
*/   of real numbers to integers is explicitly flagged for
*/   truncation or for nearest integer. 
*d njoy.1062
      if (nin.eq.0.and.a(5).le.npage) nw=6+nint(a(5))
*d njoy.1213
      zero=0.d0
*d njoy.1238
      f=0
*d njoy.1246
      n=int(log10(abs(xx)))
*d njoy.1281
  130 n=int(log10(abs(x)))
*d njoy.1287
      f=f/ten
*d njoy.1594
      z=0
*d njoy.1627
      z=0
*d njoy.1661
      z=0
*d njoy.1696
      z=0
*d njoy.1831
      x=j
      x=x/nl
*d njoy.1834
      inow=nint(xnow)
*i njoy.1940
      zero=0
*d njoy.1963
  500 if (y1.eq.zero) go to 100
*d njoy.1968
  600 if (y1.eq.zero) go to 100
*d njoy.1971
      b=b/(1/t-1/sqrt(x2-thr6))
*d njoy.1993
      f=0
*i njoy.2062
*if sw
      data shade/1.00001d0/
      data xbig/1.d10/
      data zero/0.d0/
*else
      data shade/1.00001e0/
      data xbig/1.e10/
      data zero/0.e0/
*endif
*d njoy.2114,2116
  150 if (x.lt.shade*a(jp)) go to 160
      y=0
      xnext=xbig
*d njoy.2119,2120
      xnext=xbig
      if (y.gt.zero) xnext=shade*a(jp)
*d njoy.2124
  170 y=0
*i njoy.2139
*if sw
      data break/0.1d0/
      data zero/0.d0/
      data one/1.d0/
*else
      data break/0.1e0/
      data zero/0.e0/
      data one/1.e0/
*endif
*d njoy.2142
      gral=0
*d njoy.2156
      gral=(x2-x1)*(a+b*(x2+x1)/2)
*d njoy.2160
  300 if ((xl.le.zero).or.(xh.le.zero)) go to 200
*d njoy.2163,2165
      if (abs(z).gt.break) go to 310
      gral=(x2-x1)*(yl+b*log(x1/xl))+b*x1*z*z*(1+
     1   z*(-one/3+z*(one/6-z/10)))/2
*d njoy.2171
  400 if ((yl.lt.zero).or.(yh.lt.zero)) go to 200
*d njoy.2175,2176
      if (abs(z).gt.break) go to 410
      gral=exp(a+b*x1)*(x2-x1)*(1+z*(one/2+z/6))
*d njoy.2178
  410 gral=exp(a+b*x1)*(exp(z)-1)/b
*d njoy.2182,2183
  500 if ((xl.le.zero).or.(xh.le.zero))go to 400
      if ((yl.le.zero).or.(yh.le.zero)) go to 300
*d njoy.2185,2187
      z=(b+1)*log(x2/x1)
      if (abs(z).gt.break) go to 510
      gral=yl*x1*((x1/xl)**b)*log(x2/x1)*(1+z*(one/2+z/6))
*d njoy.2189
  510 gral=yl*x1*((x1/xl)**b)*(((x2/x1)**(b+1))-1)/(b+1)
*i njoy.2206
*if sw
      data down/.99999d0/
      data xbig/1.d10/
      data zero/0.d0/
*else
      data down/.99999e0/
      data xbig/1.e10/
      data zero/0.e0/
*endif
*d njoy.2210
      if (x.gt.zero) go to 110
*d njoy.2227,2228
      if (ylast.ne.zero) xlast=down*xlast
      if (ylast.ne.zero) return
*d njoy.2230
      if (a(loc+3).ne.zero) return
*d njoy.2283
  200 y1=0
*d njoy.2289,2290
      xlast=xbig
      xnext=xbig
*i njoy.2309
*if sw
      data down/.99999d0/
      data zero/0.d0/
      data xbig/1.d10/
*else
      data down/.99999e0/
      data zero/0.e0/
      data xbig/1.e10/
*endif
*d njoy.2313
      if (x.gt.zero) go to 110
*d njoy.2328
      if (ylast.ne.zero) xlast=down*xlast
*d njoy.2368
  200 y2=0
*d njoy.2374,2375,2290
      xlast=xbig
      xnext=xbig
*d njoy.2408
      a(i)=0
*i njoy.2436
      zero=0
*d njoy.2455
      if (a(3*i+4).ne.zero) nidtot=nidtot+1
*d njoy.2486
      if (a(3*jd+4).eq.zero) go to 140
*i njoy.2655
      data zero/0.d0/
*i njoy.2667
      data zero/0.e0/
*d njoy.2669,2670
      xx=0
      if (x.eq.zero) go to 110
*d njoy.2673
      if (aa.lt.zero) ipwr=ipwr-1
*d njoy.2672
      ipwr=int(aa)
*ident up47
*/ reconr -- 4dec98
*/   make sure this step factor is not a multiple of 2.
*/   without this patch, the if statement can be sensitive
*/   to differences in the 10th or 11th place, thus giving 
*/   different results with different compilers.
*i reconr.1606
      data estp/4.1d0/
*i reconr.1610
      data estp/4.1e0/
*d reconr.1753
  145 est=estp*(a(ix+1-1)-res(1))
      if (in.gt.3.and.dx.gt.est) go to 175
*/ reconr -- 4dec98
*/   at low energies in nonfissionable materials, the r-function
*/   in csrmat is obtained from a degenerate r matrix that is very
*/   close to unity.  this leads to difficult numerics for low
*/   incident energies.  this patch rearranges the terms to allow
*/   unit terms to cancel out, thus giving better results for the
*/   total and capture cross sections for materials in the
*/   fe,ni,cr range.  the old approach could show differences
*/   between different machines starting in the 6th significant
*/   figure of the capture cross section.  these differences are
*/   only significant when comparing results from different systems.
*i reconr.2288
      data small/3.d-4/
*i reconr.2294
      data small/3.e-4/
*d reconr.2368
*i reconr.2418
      r(1,1)=uno+r(1,1)
      r(2,2)=uno+r(2,2)
      r(3,3)=uno+r(3,3)
*i reconr.2432
      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)
*d reconr.2435,2445
  160 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
  170 if (jj.gt.jjl.and.jj.lt.numj) then
*/ reconr -- 4dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*i reconr.102
*if sw
      data elow,ehigh,elarge/1.d-5,20.d6,9.d9/
*else
      data elow,ehigh,elarge/1.e-5,20.e6,9.e9/
*endif
*d reconr.193,194
      zain=1
      awin=1
*d reconr.198
      zain=int(nsub/10)
*d reconr.221,225
      eresl=elow
      eresr=ehigh
      eresh=ehigh
      eresu=elarge
      eresm=elarge
*i reconr.288
      zero=0
*d reconr.293,295
      tempr=0
      errmax=-1
      errint=-1
*d reconr.297
      if (errmax.le.zero) errmax=20*err
*d reconr.299
      if (errint.le.zero) errint=err/10000
*d reconr.488
      common/runits/nin,nout,nscrr,nscrs
*d reconr.3754
      common/runits/nin,nout,nscrr,nscrs
*i reconr.506
      data elow/1.d-5/
      data elarge/9.d9/
      data small/1.d-6/
      data wide/3.d0/
*i reconr.517
      data elow/1.e-5/
      data elarge/9.e9/
      data small/1.e-6/
      data wide/3.e0/
*d reconr.530,536
      eresl=elarge
      eresh=0
      eresr=0
      eresu=elarge
      eresm=elarge
      zero=0
*d reconr.569
      if (abs(el-elow).gt.small) go to 107
*d reconr.666
      if (enut.gt.ener+ener/100) go to 148
*d reconr.713
      if (ascatl.eq.zero) ascatl=ascat
*d reconr.759
      a(jj+2)=0
*d reconr.776
      if (enut.gt.ener+ener/100) go to 229
*d reconr.846
      if (enut.gt.ener+ener/100) go to 339
*d reconr.871
      hw=0
*d reconr.929
      el=0
*d reconr.1009
         se=0
*i reconr.1057
*if sw
      data big/1.d10/
*else
      data big/1.e10/
*endif
*i reconr.1064
      zero=0
*d reconr.1080
      a(isunr+12)=big
*d reconr.1083
      e=0
*d reconr.1092
      a(l+k)=0
*d reconr.1126
      e=0
*d reconr.1131
      if (e.lt.zero) go to 160
*i reconr.1205
      elow=1.d-5
      elim=.99d6
      up=1.001d0
      dn=0.999d0
*i reconr.1212
      elow=1.e-5
      elim=.99e6
      up=1.001e0
      dn=0.999e0
*i reconr.1233
      zero=0
*d reconr.1277
      if (mfh.eq.13) a(iscr+4)=1
*d reconr.1306
      if (thr6.lt.0) thr6=0
*d reconr.1310
      if (a(l+1).eq.zero) go to 183
*d reconr.1314
      a(l+1)=0
*d reconr.1372
  220 test=elow
*d reconr.1376
      if (sr.eq.zero) go to 225
*d reconr.1406
  245 et=0
*d reconr.1475,1476
      if (a(ix+i-1).gt.elim) go to 320
*d reconr.1479
      if (a(ix+i-2).gt.up*xm.and.a(ix+i-1).lt.dn*xm) go to 315
*d reconr.1482
      if (a(ix+i-2).gt.up*xm.and.a(ix+i-1).lt.dn*xm) go to 315
*d reconr.1485
      if (a(ix+i-2).gt.up*xm.and.a(ix+i-1).lt.dn*xm) go to 315
*d reconr.1488
      if (a(ix+i-2).gt.up*xm.and.a(ix+i-1).lt.dn*xm) go to 315 
*d reconr.1497,1498
      if (a(ix+i-1).gt.elim) go to 325
*d reconr.1552
      egl=0
*d reconr.1556
      if (eg.lt.zero) ngneg=ngneg+1
*d reconr.1631,1636
      crej=0
      cmax=0
      cint=0
      eint=0
      erej=0
      emax=0
      zero=0
*d reconr.1649
      e=0
*d reconr.1700,1702
      cb=(c1+c2)/2
      if (dc.lt.cb/2.and.i.ge.3) go to 120
      crej=crej+dc*dx/(2*xm)
*d reconr.1706
      crej=crej+abs(cc*dx3/(6*dx2))/xm
*d reconr.1710,1712
      cb=(c1+c2)/2
      if (dc.lt.cb/2.and.i.ge.3) go to 130
      erej=erej+dc*dx/(2*xm)
*d reconr.1716
      erej=erej+abs(cc*dx3/(6*dx2))/xm
*d reconr.1727
      fr1=1-fr2
*d reconr.1730,1731
      dm2=0
      if (sig(3).eq.zero) go to 140
*d reconr.1748,1749
      cmax=cmax+dm3*dx/(2*xm)
      emax=emax+dm1*dx/(2*xm)
*d reconr.1770
      cint=cint+(c1+c2)*dx/(2*xm)
*d reconr.1773
      eint=eint+(c1+c2)*dx/(2*xm)
*d reconr.1778,1782
      if (cint.eq.zero.or.eint.eq.zero) go to 160
      cmax=100*cmax/cint
      emax=100*emax/eint
      crej=100*crej/cint
      erej=100*erej/eint
*d reconr.1798,1803
      cint=0
      eint=0
      cmax=0
      emax=0
      crej=0
      erej=0
*d reconr.1879
      sig(i)=0
*i reconr.1880
      zero=0
*d reconr.1917
      if (sigp(j).lt.zero) sigp(j)=0
*d reconr.1951
      sigp(i)=0
*i reconr.1992,1993
      zero=0
      if (tempr.gt.zero) call findex('tr',itr,aa)
      if (tempr.gt.zero) call findex('ti',iti,aa)
*d reconr.1997
      sigp(i)=0
*d reconr.2030
      pec=0
*d reconr.2041
      if (tempr.eq.zero) go to 35
*d reconr.2157
      zero=0
      if (tempr.gt.zero)
*d reconr.2163
      sigp(i)=0
*d reconr.2192
      pec=0
*d reconr.2200
      sum=0
*d reconr.2204
      nj=nint(ajmax-ajmin+1)
*d reconr.2214
      sigj(i,ii)=0
*d reconr.2222
      j=nint(a(inow+1)-ajmin+1)
*d reconr.2298
      zero=0
      if (tempr.gt.zero)
*d reconr.2304
      sigp(i)=0
*d reconr.2324,2326
      gfa=0
      gfb=0
      gf=0
*d reconr.2337,2338
      if (apl.ne.zero) rhoc=k*apl
      if (apl.ne.zero.and.naps.eq.1) rho=k*apl
*d reconr.2351
      numj=nint(ajmax-ajmin+1)
*d reconr.2365,2366
      s(j,i)=0
      r(j,i)=0
*d reconr.2376,2377
      if (abs(aj-ajc).gt.quar) go to 140
*d reconr.2387,2392
      a2=0
      if (apl.ne.zero) rhoc=k*apl
      if (apl.ne.zero.and.naps.eq.1) rho=k*apl
      a3=0
      if (gfb.ne.zero) a3=sqrt(abs(gfb))
      if (gfb.lt.zero) a3=-a3
*d reconr.2402
      if (gfa.eq.zero.and.gfb.eq.zero) go to 140
*d reconr.2413
      gf=1
*d reconr.2417
      if (gf.eq.zero) go to 160
*d reconr.2558
      c(i,j)=0
*d reconr.2599
      if (tempr.gt.zero) call error('cshybr',
*d reconr.2604
      sigp(i)=0
*d reconr.2624
      sigc(n)=0
*d reconr.2685
      gc(ii)=0
*d reconr.2819
      zero=0
      if (e.gt.zero) go to 125
*d reconr.2830
      elast=0
*d reconr.2832
      elst(isect)=0
*d reconr.2858
      sigp(i)=0
*d reconr.2891,2892
      gfx=0
      gxx=0
*d reconr.2915
      mu=int(amun)
*i reconr.3019
      data small/1.d-8/
*i reconr.3025
      data small/1.e-8/
*d reconr.3029
      zero=0
      if (e.gt.zero) go to 125
*d reconr.3042
      elast=0
*d reconr.3044
      elst(isect)=0
*d reconr.3068
      sigp(i)=0
*d reconr.3115,3117
      if (gxx.lt.small) gxx=0
      if (gfx.lt.small) gfx=0
*d reconr.3213,3214
      zero=0
      if (tempr.gt.zero) call findex('tr',itr,aa)
      if (tempr.gt.zero) call findex('ti',iti,aa)
*d reconr.3218
      sigp(i)=0
*d reconr.3229
      if (tempr.gt.zero) delta=1/sqrt(4*bk*tempr*e/awri)
*d reconr.3240,3242
      bakt=0
      bakf=0
      bakc=0
*d reconr.3277
      if (tempr.gt.zero) go to 155
*i reconr.3319
*if sw
      data test/1.d-6/
*else
      data test/1.e-6/
*endif
*d reconr.3327
*d reconr.3328
         if ((phi/rho).lt.test) phi=0
*d reconr.3331
*d reconr.3332
         if ((phi/rho).lt.test) phi=0
*d reconr.3338
*d reconr.3339
         if ((phi/rho).lt.test) phi=0
*d reconr.3411
      zero=0
      s=0
*d reconr.3412,3417
      if (galpha.le.zero) go to 1000
      if (gamma.le.zero) go to 1000
      if (gbeta.lt.zero) go to 1000
      if (gbeta.gt.zero) go to 350
      if (df.lt.zero) go to 1000
      if (df.gt.zero) go to 250
*d reconr.3442,3443
  350 if (df.lt.zero) go to 1000
      if (df.gt.zero) go to 450
*i reconr.3519
      data finity/.99d10/
      data small/1.d-8/
*i reconr.3521
      data finity/.99e10/
      data small/1.e-8/
*d reconr.3535,3536
      zero=0
*d reconr.3540
      tot(i+1)=0
*d reconr.3560,3562
      inn=in
      if (ig.eq.ngp.and.en.eq.eg) inn=-in
      call loada(inn,tot,ntot,iold,a(ibufo),nbuf)
*d reconr.3591
      e=0
*d reconr.3594
      if (sn.ne.zero) thresh=sigfig(thresh,7,-1)
*d reconr.3620
      sn=0
*d reconr.3623
      if (thresh.gt.1..and.abs(thresh-eg).lt.test) sn=0
*d reconr.3626
      if (eg.ge.eresr.and.eg.lt.eresh.and.itype.gt.0) sn=0
*d reconr.3648,3650
      inn=in
      if (ig.eq.ngo) inn=-in
      call loada(inn,tot,2,ngrid,a(ibufg),nbufg)
*d reconr.3695
      a(incs+nxc-1)=a(incs+nxc-1)+2+int((n2h+2)/3)
*i reconr.3770
*if sw
      data emin/1.d-5/
*else
      data emin/1.e-5/
*endif
*d reconr.3802,3803
      a(iscr)=0
      a(iscr+1)=0
*d reconr.3812,3816
      a(iscr+1)=0
      a(iscr+2)=0
      a(iscr+3)=0
      a(iscr+4)=int(10*zain)
      a(iscr+5)=0
*d reconr.3854,3855
      a(j+1)=0
      a(j+2)=0
*d reconr.3859
      a(j+5)=3+int((np+2)/3)
*d reconr.3900
      a(iscr+1)=0
*d reconr.3902
      a(iscr)=emin
*d reconr.3937
      a(iscr+2)=0
*d reconr.3951
      a(iscr+1)=0
*d reconr.4055,4056
      zero=0
      aki=1
      if (ax.lt.zero) aki=-1
*d reconr.4066,4067
   10 ii=int(x*10)
      jj=int(y*10)
*i reconr.4124
*if sw
      data tenth/.1d0/
*else
      data tenth/.1e0/
*endif
*d reconr.4126,4129
      x0=-tenth
      y0=-tenth
      dx=tenth
      dy=tenth
*d reconr.4179,4180
      zero=0
      rew=0
      aimw=0
*d reconr.4183,4185
      if (abrez+aimz.ne.zero) go to 20
      rew=1
      aimw=0
*d reconr.4192,4195
      if (abrez+brk1*aimz-brk2.gt.zero) go to 350
      if (abrez+brk3*aimz-brk4.gt.zero) go to 350
      if (r2+brk5*ai2-brk6.lt.zero) go to 340
      if (r2+brk7*ai2-brk8.ge.zero) go to 340
      if (aimz-brk9.ge.zero) go to 350
*d reconr.4200
      if (aim1.ge.zero) go to 370
*d reconr.4209,4213
      h=0
      b=0
      a=0
      tempm=0
      temel=0
*d reconr.4239,4246
      c=0
      b=0
      ajsig=0
      d=0
      jsig=0
      g=0
      h=0
      el=0
*d reconr.4287
      if ((abs(tempm)+abs(temel)-up).lt.zero) go to 500
*d reconr.4297
  500 if ((abs(tempm)+abs(temel)-dn).gt.zero) go to 520
*ident up48
*/ broadr -- 11dec98
*/   make sure this step factor is not a multiple of 2.
*/   without this patch, the if statement can be sensitive
*/   to differences in the 10th or 11th place, thus giving 
*/   different results with different compilers.
*i broadr.997
      data estp/4.1d0/
*i broadr.1004
      data estp/4.1e0/
*d broadr.1105
      est=estp*(es(is)-tt(1))
      if (j.gt.3.and.dx.gt.est) go to 170
*/ broadr -- 11dec98
*/   tighten tolerance for printing thermal quantities.
*/   fix problem in capture resonance integrals.
*d up11.31
      if (abs(tev-therm).gt.therm/1000) go to 2250
*d up2.60
            ssc=(sc/enow+ss/elast)*(enow-eone)/2
*/ broadr -- 11dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*i broadr.119
      fact=.99999d0
      onemev=1.d6
      tbreak=0.5d0
      big=9.d9
      step=1.01d0
*i broadr.122
      fact=.99999e0
      onemev=1.e6
      tbreak=0.5e0
      big=9.e9
      step=1.01e0
*d broadr.130
*i broadr.131
      zero=0
*d broadr.158
      temp1=0
*d broadr.162.164
      thnmx=onemev
      errmax=0
      errint=0
*d broadr.166,167
      if (errmax.eq.zero) errmax=20*errthn
      if (errint.eq.zero) errint=errthn/10000
*d broadr.244,245
      emin=1
      if (thnmax.lt.zero) emin=-thnmax
*d broadr.276
      test=1+temp1/1000
*d broadr.290
      enow=0
*d broadr.303
      if (en.gt.zero) sig=sig-sun
*d broadr.306,307
      if (enext.ge.big) kt=-kt
*d broadr.311
      if (enext.lt.big) go to 150
*d broadr.320
      enow=0
*d broadr.340
      if (en.gt.zero) sig=sig-sun
*d broadr.352
      if (thnmax.lt.zero) thnmax=-thnmax
*d broadr.384
      break=9999
*d broadr.401
      if (mti(k).eq.0.and.tt(1+k).gt.zero) mti(k)=j
*d broadr.428
      if (mt.eq.1.or.mt.eq.3.or.mt.eq.19) a(iscr+i+3)=3+int((n2out+2)/3)
*d broadr.435
      a(iscr+i+3)=3+int((n2out-mt4i+3)/3)
*d broadr.439
      if (mt.eq.mtr(j)) a(iscr+i+3)=3+int((n2out-jjj+3)/3)
*d broadr.447
      eone=tbreak
*d broadr.460,469
      fint=0
      cint=0
      alint=0
      etint=0
      v1int=0
      ssf=0
      slf=0
      ssc=0
      slc=0
      elast=0
      flast=0
      fnul=0
*d up2.17,19
 2240 enow=tt(1)
      if (enow.gt.step*elast.and.elast.gt.zero) then
*d broadr.521
      fnu=0
*d up11.42
      break=1
*d up2.44
               if (sf.gt.zero.and.slf.gt.zero)
*d broadr.584
      if (ctev.ne.zero) write(nsyso,'(
*d broadr.596
         if (ftev.ne.zero) write(nsyso,'(
*d broadr.651
  275 ee=0
*d broadr.681
      if (en.gt.zero) a(k)=a(k)+sun
*d broadr.684
      a(k)=0
*d broadr.694
      if (en.gt.zero) a(k)=a(k)+sun
*i broadr.763
      zero=0
*d broadr.811,812
      if (i.eq.1.or.e.lt.zero) eumax=abs(e)
      if (e.lt.zero) a(l)=-e
*i broadr.834
*if sw
      data up/1.00001d0/
*else
      data up/1.00001e0/
*endif
*d broadr.837
      sun=0
*d broadr.855
      enext=eumax*up
*i broadr.883
      zero=0
*d broadr.911
  120 if (tempef.eq.zero) go to 130
*i broadr.1012
      zero=0
*d broadr.1028,1029
      if (abs(dn).lt.abs(s(i,k))/1000) dn=0
      dl(i)=1
*d broadr.1031
      if (dn.lt.zero) dl(i)=-1
*d broadr.1054
      if (abs(dn).lt.abs(s(i,k))/1000) dn=0
*d broadr.1055,1056
      if (dn.ge.zero) then
         dn=1
*d broadr.1058
         dn=-1
*i broadr.997
      one=1.d0
*i broadr.1004
      one=1.e0
*d broadr.1092
      test=one-one/100
*d broadr.1097
      si=f*ss(i,is-1)+(1-f)*ss(i,is)
*d broadr.1101
      if (dy*dx/2.gt.errint*em) go to 170
*d broadr.1122
      if (abs(dn).lt.abs(s(i,k))/1000) dn=0
*d broadr.1124
         dl(i)=1
*d broadr.1126
         dl(i)=-1
*d broadr.1203
      aa=0
*d broadr.1208
      sbt(i)=0
*d broadr.1246,1247
      x=0
      xx=0
*d broadr.1295
      xx=0
*d broadr.1323
      if (sbt(i).lt.sigmin) sbt(i)=0
*d broadr.1353
*if sw
      data stpmax/1.24d0/
*else
      data stpmax/1.24e0/
*endif
*i broadr.1354
      zero=0
*d broadr.1372,1373
      if (denom.eq.zero) go to 150
      denom=1/denom
*i broadr.1508
      data aerr,rerr/1.d-30,1.d-8/
*i broadr.1512
      data aerr,rerr/1.e-30,1.e-8/
*d broadr.1519
      sign=1
*d broadr.1537,1538
      cm(1)=1
      xk=1
*d broadr.1562
      beta=0
*d broadr.1570
  160 xk=1
*d broadr.1578
      test=aerr+rerr*abs(s)
*ident up49
*/ unresr -- 11dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*/   the common with the unit numbers was changed to be the
*/   same in every subroutine to avoid confusion.
*d unresr.187,190
         bkgz(1)=0
         bkgz(2)=0
         bkgz(3)=0
         bkgz(4)=0
*d unresr.224
      if (nx.gt.0.and.new.gt.0) a(iscr+5)=a(iscr+5)+1
*d unresr.232
      if (nx.gt.0.and.new.gt.0) a(iscr+5)=a(iscr+5)+1
*d unresr.351
      common/uunits/nendf,nin,nout,nscr
*d unresr.353
      common/unen/nunr,matd,nsigz,nmtx,mtx(5),npnts(5)
*i unresr.379
*if sw
      data onemev/1.d6/
      data small/1.d-12/
      data step/1.01d0/
*else
      data onemev/1.e6/
      data small/1.e-12/
      data step/1.01e0/
*endif
*d unresr.388,389
      zero=0
      elr=0
      ehr=onemev
*d unresr.391
      a(ieunr)=onemev
*d unresr.400
      call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.407
  115 call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.417
      call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.420
      call listio(nendf,0,0,a(iscr),nb,nw)
*d unresr.422
      call moreio(nendf,0,0,a(iscr),nb,nw)
*d unresr.426
      call listio(nendf,0,0,a(iscr),nb,nw)
*d unresr.428
      call moreio(nendf,0,0,a(iscr),nb,nw)
*d unresr.430
  126 call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.433
      call listio(nendf,0,0,a(iscr),nb,nw)
*d unresr.435
      call moreio(nendf,0,0,a(iscr),nb,nw)
*d unresr.469
      call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.476
      call listio(nendf,0,0,a(iscr),nb,nw)
*d unresr.505
  190 call listio(nendf,0,0,a(iscr),nb,nw)
*d unresr.527
      call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.537
      call listio(nendf,0,0,a(iscr),nb,nw)
*d unresr.558
      call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.567
      call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.575
      call listio(nendf,0,0,a(iscr),nb,nw)
*d unresr.578
      call moreio(nendf,0,0,a(loc),nb,nw)
*d unresr.486
      a(inow+5)=0
*d unresr.596
      if (a(jnow+k).le.zero) a(jnow+k)=small
*d unresr.629,630
      if (enext.ge.onemev) go to 380
*d unresr.635
      if (enut.gt.step*et) go to 363
*d unresr.648
      en=0
*d unresr.653
      if (et.lt.zero) iovlp=1
*d unresr.706
      common/uunits/nendf,nin,nout,nscr
*i unresr.710
*if sw
      data up,dn/1.00001d0,0.99999d0/
*else
      data up,dn/1.00001e0,0.99999e0/
*endif
*d unresr.717
  100 call contio(nendf,0,0,a(iscr),nb,nw)
*d unresr.724
      a(ie+ibase)=0
*d unresr.727
  120 call tosend(nendf,0,0,a(iscr))
*d unresr.730,731
      e=0
      call gety1(e,enext,idis,bkg,nendf,a(iscr))
*d unresr.736,738
      if (ie.eq.1) e=up*e
      if (ie.eq.nunr) e=dn*e
      call gety1(e,enext,idis,a(ibase+ie),nendf,a(iscr))
*d unresr.740
      call tosend(nendf,0,0,a(iscr))
*i unresr.796
      data small/1.d-8/
*i unresr.819
      data small/1.e-8/
*d unresr.829,830
      if (temp(1).lt.one) temp(1)=one
*d unresr.835,836
      spot=0
      sint=0
*d unresr.900,901
      amux=0
      gxx=0
*d unresr.903
      amuf=0
*d unresr.909
      gfx=0
*d unresr.950
      tk(ks,is0,itp)=0
*d unresr.952
      t(kx,is0,itp)=0
*d unresr.957
*d unresr.959
      if (gfx.le.small) mu=0
*d unresr.962
      if (gxx.lt.small) lu=0
*d unresr.1022,1026
      yy(1)=0
      yy(2)=0
      yy(3)=0
      yj=0
      yk=0
*d unresr.1028,1029
      xk=0
      xj=0
*i unresr.1118
*if sw
      data small/1.d-8/
*else
      data small/1.e-8/
*endif
*d unresr.1131,1133
      if (gxx.lt.small) gxx=0
      if (gfx.lt.small) gfx=0
*i unresr.1156
      zero=0
*d unresr.1158,1159
      aki=1
      if (ax.lt.zero) aki=-1
*d unresr.1169,1170
   10 ii=int(x*10)
      jj=int(y*10)
*d unresr.1196
      a2=2*x*y
*i unresr.1241
      zero=0
*d unresr.1244
      if (ep.gt.zero) go to 20
*d unresr.1249,1252
      y=0
      z=0
      yk=0
      zk=0
*d unresr.1254
      rki=0
*d unresr.1268
      x2=1/(1+b1/rew)
*i unresr.1297
*if sw
      data tenth/0.1d0/
*else
      data tenth/0.1e0/
*endif
*d unresr.1299,1302
      x0=-tenth
      y0=-tenth
      dx=tenth
      dy=tenth
*d unresr.1352,1353
      zero=0
      rew=0
      aimw=0
*d unresr.1356,1358
      if (abrez+aimz.ne.zero) go to 20
      rew=1
      aimw=0
*d unresr.1363,1367
      if (abrez+brk1*aimz-brk2.gt.zero) go to 350
      if (abrez+brk3*aimz-brk4.gt.zero) go to 350
      if (r2+brk5*ai2-brk6.lt.zero) go to 340
      if (r2+brk7*ai2-brk8.ge.zero) go to 340
      if (aimz-brk9.ge.zero) go to 350
*d unresr.1373
      if (aim1.ge.zero) go to 370
*d unresr.1381,1385
      h=0
      b=0
      a=0
      tempm=0
      temel=0
*d unresr.1411,1418
      c=0
      b=0
      ajsig=0
      d=0
      jsig=0
      g=0
      h=0
      el=0
*d unresr.1459
      if ((abs(tempm)+abs(temel)-up).lt.zero) go to 500
*d unresr.1469
  500 if ((abs(tempm)+abs(temel)-dn).gt.zero) go to 520
*ident up50
*/ heatr -- 11 dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*d heatr.735,739
      qs=0
      q0=0
      icon=0
      ebar=0
      dame=0
*d heatr.787
      q0=0
*d heatr.832
      yld=1
*d heatr.919
      elst=0
*d heatr.921,922
      if (icon.eq.0) yld=0
      yld0=0
*d heatr.930
      if (new6.eq.1) c(npkk)=0
*d heatr.985
      ebal6=0
*d heatr.993,996
      if (izap.eq.1) h=0
      if (izap.eq.0.and.local.eq.0) h=0
      if (izap.le.1) dame=0
      ebal6=0
*i heatr.2528
      zero=0
      one=1
*d heatr.2596
      if (y(i-1,1).ne.zero) f1=y(i-1,2)/y(i-1,1)
*d heatr.2651
      if (yy.ne.zero) b=b/yy
*d heatr.2708
      f=one/2
*d heatr.2823
      zero=0
      if (ep.gt.zero) go to 200
*d heatr.3035
      zero=0
      if (ep.gt.zero) go to 120
*d heatr.3044
      if (cnow(lnow).le.zero) lnow=lnow+ncnow
*d heatr.3090
      if (t.lt.zero) t=0
*d heatr.3122
      if (t.lt.zero) t=0
*d heatr.3180
      zero=0
      if (i.gt.zero) go to 120
*d heatr.3208
      if (t.lt.zero) t=0
*d heatr.3234
      if (t.lt.zero) t=0
*d heatr.3470
      awc=awr+1
      h=h+df(xr,z,awc,z,awr)*y
*d heatr.3526
      zero=0
      if (e.gt.zero) go to 200
*d heatr.3612
  320 fle(i)=0
*d heatr.3620
  410 fle(1)=1
*d heatr.3622
      fle(il)=0
*d heatr.3630
  440 fle(1)=1
*d heatr.3633
      fle(il)=0
*d heatr.3754
  220 fl(il)=0
*d heatr.3762
      fl(il)=0
*d heatr.3802
      if (abs(fl(il)).gt.zero) nlz=il
*d heatr.3871
      a(ie)=0
*d heatr.3875,3877
      a(ir+m)=0
      if (i.eq.j) a(ir+m)=1
      a(iaa+m)=0
*d heatr.3885
      g=1
*d heatr.3905,3906
      a(iaa+(k-1)*imax+j-1)=0
      a(ir+(k-1)*imax+j-1)=0
*d heatr.3909
      g=1
*d heatr.3921
      ysum=0
*d heatr.3937
      a(ie)=0
*d heatr.3973
      a(iscr+1)=0
*d heatr.3989
      a(iscr+2)=0
*d heatr.4123
      e=0
*d heatr.4176
      cerr=0
*d heatr.4179
      e=0
*ident up51
*/ thermr -- 15dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*/   in addition, any conversion of real numbers to integers is 
*/   explicitly flagged for truncation or for nearest integer. 
*i thermr.114
      data up/1.00001d0/
      data zero/0.d0/
*i thermr.122
      data up/1.00001e0/
      data zero/0.e0/
*d thermr.163,164
      eftemp(i)=0
      eftmp2(i)=0
*d thermr.185,186
      sz2=0
      az2=0
*d thermr.192
      if (az2.ne.zero) sb2=sz2*((az2+1)/az2)**2
*d thermr.198
      if (sz2.eq.zero) go to 110
*d thermr.238
      if (sz2.eq.zero) go to 115
*d thermr.328
      e=0
*d thermr.336
      if (e.gt.emax*(1+small)) ex(2)=0
*d thermr.342
      if (abs(e-emax).lt.small*emax) enext=up*emax
*d thermr.379
      if (sz2.gt.zero) teff2=eftmp2(itemp)
*d thermr.446
  130 if (abs(tnow-temp).lt.temp/1000+5) return
*i thermr.554
      zero=0
*d thermr.557
      if (eftemp(i).ne.zero) go to 110
*d thermr.559
      jmat=nint(tabl(1,j))
*d thermr.561
      test=5
*d thermr.565
      if (eftemp(i).eq.zero) eftemp(i)=temp(i)
*i thermr.601
      zero=0
*d thermr.604
      if (eftmp2(i).ne.zero) go to 110
*d thermr.606
      jmat=nint(tabl(1,j))
*d thermr.608
      test=5
*d thermr.612
      if (eftmp2(i).eq.zero) eftmp2(i)=temp(i)
*d thermr.627
      common/tunits/nendf,nin,nout,nscr,nscr2
*i thermr.638
      data tolmin/1.d-6/
      data eps/3.d-5/
*i thermr.641
      data tolmin/1.e-6/
      data eps/3.e-5/
*d thermr.654
      ex(nex+i)=0
*d thermr.659
*d thermr.668
      e=0
*d thermr.712
      if (a(istk+nx*(i-2))-a(istk+nx*(i-1)).lt.eps*xm) go to 160
*d thermr.740
      if (a(istk+nx*(i-1)).gt.emax*(1+small)) ej(nex)=0
*d thermr.743
      if (a(istk+nx*(i-1)).gt.emax*(1+small)) ej(nex+il)=0
*i thermr.844
      data eps/0.05d0/
      data zero/0.d0/
*i thermr.863
      data eps/0.05e0/
      data zero/0.e0/
*d thermr.868
      if (e.gt.zero) go to 210
*d thermr.873
*d thermr.911
      i1m=int(a*sqrt(phi))
*d thermr.916
      i2m=int(half*(l1+sqrt(3*(a*a*phi-l1*l1))))
*d thermr.922
      if (x.gt.zero) i3m=int(c*sqrt(x))
*d thermr.926
      w1=2
*d thermr.928
      w2=2
*d thermr.945
      if (tsq.lt.b(ifl+2*i-2).or.tsq.ge.(1+eps)*b(ifl+2*i-2))
*d thermr.1003
  200 k=int(b(ifl+5))
*d thermr.1005
      blast=0
*d thermr.1022
      s(il)=0
*d thermr.1029
      if (e.gt.emax) f=0
*i thermr.1118
      data up,dn/1.1d0,.9d0/
*i thermr.1125
      data up,dn/1.1e0,.9e0/
*d thermr.1158
      if (abs(temp-tt1).gt.temp/10) call error('iel',
*d thermr.1164
      if (temp.lt.dn*tt1.or.temp.gt.up*ttn) call error('iel',
*d thermr.1213
      a(iscr)=0
*d thermr.1237
      if (iex.eq.ne) ej(nj)=0
*d thermr.1336
      sum=0
*i thermr.1401
      data zero/0.d0/
      data tenth/0.1d0/
      data up,dn/1.1d0,0.9d0/
      data uumin/0.00001d0/
*i thermr.1418
      data zero/0.e0/
      data tenth/0.1e0/
      data up,dn/1.1e0,0.9e0/
      data uumin/0.00001e0/
*d thermr.1473,1474
      sb2=0
      az2=0
*d thermr.1476,1477
      if (a(iscr+12).eq.zero) az2=a(iscr+14)
      if (az2.ne.zero) sb2=a(iscr+13)*((az2+1)/az2)**2
*d thermr.1480
      if (a(iscr+18).eq.zero)
*d thermr.1509
      if (abs(t-temp).le.(t/50)) go to 130
*d thermr.1528
      if (t.eq.t1) l=l+2*nint(a(loc+4))+1
*d thermr.1532
         if (a(l).gt.zero) itrunc=1
*d thermr.1565
      test=5
*d thermr.1587
      if (abs(tempt-tt1).gt.tempt/10) call error('calcem',
*d thermr.1596
      if (tempt.lt.dn*tt1.or.tempt.gt.up*ttn) call error('calcem',
*d thermr.1601
  190 if (sb2.eq.zero) go to 300
*d thermr.1607
      if (abs(tempt-tt1).gt.tempt/10) call error('calcem',
*d thermr.1616
      if (tempt.lt.dn*tt1.or.tempt.gt.up*ttn) call error('calcem',
*d thermr.1630,1632
      a(ibeta)=0
      a(ibeta+1)=tenth
      a(ibeta+2)=25
*d thermr.1666
      cliq=0
*d thermr.1678,1680
      a(ixs+ie-1)=0
      ubar(ie)=0
      ep=0
*d thermr.1704
  315 if (ep.gt.zero) go to 316
*d thermr.1726,1727
      uu=0
      uum=0
*d thermr.1737
      test=2*tol*abs(uu)+uumin
*d thermr.1745
         uu=0
*d thermr.1768
      if (ylast.ne.zero) jnz=j
*d thermr.1769
      ulast=0
*d thermr.1779
      y(il,i)=0
*d thermr.1794
      uu=0
*d thermr.1807
      if (y(1,1).ne.zero) jnz=j
*d thermr.1822
      a(iscr)=0
*d thermr.1840
      ex(ie)=0
*d thermr.1854
      if (ie.eq.ne) xs=0
*d thermr.1864
      if (ie.eq.ne) xs=0
*i thermr.1894
      zero=0
*d thermr.1918
      if (lasym.eq.1.and.bb.lt.zero) bbb=-b
*d thermr.1950
      sig=0
*d thermr.1953
      if (sig.lt.sigmin) sig=0
*d thermr.1959
      s=0
*d thermr.1963
      if (sb2.eq.zero) go to 190
*d thermr.1966
      s2=0
*d thermr.1969
  190 if (sig.lt.sigmin) sig=0
*d thermr.1974
      s=0
*d thermr.1978
      if (sig.lt.sigmin) sig=0
*d thermr.1983
      sig=0
*i thermr.2052
      data zero/0.d0/
*i thermr.2058
      data zero/0.e0/
*d thermr.2074
      if (ep.ne.zero) seep=1/sqrt(e*ep)
*d thermr.2078
      sum=0
*d thermr.2083,2085
      if (ep.eq.zero) x(2)=0
      if (ep.ne.zero) x(2)=half*(e+ep-(s1bb-1)*az*tev)*seep
      if (abs(x(2)).gt.1-eps) x(2)=0
*d thermr.2117
      s(il)=0
*d thermr.2127,2128
      sum=0
      gral=0
*d thermr.2130
      s(il)=0
*d thermr.2139,2141
      if (ep.eq.zero) x(2)=0
      if (ep.ne.zero) x(2)=(e+ep-(s1bb-1)*az*tev)*seep*half
      if (abs(x(2)).gt.1-eps) x(2)=0
*d thermr.2189
      if (disc.ge.zero) go to 178
*d thermr.2193,2194
  178 if (f.gt.zero) xn=xl-(yl*rf)+sqrt(disc)
      if (f.lt.zero) xn=xl-(yl*rf)-sqrt(disc)
*d thermr.2220,2221
      sum=0
      gral=0
*i thermr.2253
      data etop/20.d6/
      data zero/0.d0/
*i thermr.2255
      data etop/20.e6/
      data zero/0.e0/
*d thermr.2261
      ex(i)=0
*d thermr.2318
      a(i-1+idicn)=0
*d thermr.2426
      a(iscr+1)=0
*d thermr.2440,2441
      ex(1)=etop
      ex(ix)=0
*d thermr.2515
      if (a(ixs+ie-1).eq.zero) 
     1   call error('tpend','cross section=0.',' ')
*ident up52
*/ groupr -- 14dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*/   in addition, any conversion of real numbers to integers is 
*/   explicitly flagged for truncation or for nearest integer. 
*i groupr.271
      zero=0
*d groupr.366
      eps=tempin/10000
*d groupr.368
      if (diff.gt.1+eps) go to 220
*d groupr.421
      a(is+1)=0
*d groupr.445
      if (ngg.eq.0) a(nw)=0
*d groupr.544
      ee=0
*d groupr.548
      test=-1
*d groupr.557
      if (econst.eq.zero) go to 510
*d groupr.580
      a(it)=0
*d groupr.614
      if (econst.gt.zero) write(nsyso,'(/
*d groupr.647
      a(iscr+1)=0
*i groupr.760
*if sw
      data sigzmx/1.d10/
      data zero/0.d0/
*else
      data sigzmx/1.e10/
      data zero/0.e0/
*endif
*d groupr.798
      title(1)=0
*d groupr.805
      if (temp(1).eq.zero) write(nsyso,'(
*d groupr.807
      if (temp(1).ne.zero) write(nsyso,'(
*d groupr.814
      sigz(1)=sigzmx
*d groupr.1390,1391
      data u187a,u187b,u187c,e187d,e187e/-15.5d0,-14.125d0,-5.875d0,
     1  2.6058d4,6.868d0/
*d groupr.1499,1500
      data u187a,u187b,u187c,e187d,e187e/-15.5e0,-14.125e0,-5.875e0,
     1  2.6058e4,6.868e0/
*d groupr.1601
      eg(60)=e187e
*d groupr.1839
      egg(1)=0
*i groupr.2046
      data small/1.d-10/
      data zero/0.d0/
      data onep5/1.5d0/
*i groupr.2106
      data small/1.e-10/
      data zero/0.e0/
      data onep5/1.5e0/
*d groupr.2112
      sigpot=0
*d groupr.2115,2119
      alpha2=0
      sam=0
      beta=0
      alpha3=0
      gamma=0
*d groupr.2127,2128
      if (alpha3.ne.zero.and.alpha2.eq.zero) alpha2=small
      if (alpha2.gt.zero) write(nsyso,'(
*d groupr.2163
      ac=1/(exp(-ec/tc)*ec**onep5)
*i groupr.2261
      data exmin/-89.d0/
      data zero/0.d0/
*i groupr.2276
      data exmin/-89.e0/
      data zero/0.e0/
*d groupr.2308
      if (wtf.eq.zero) go to 120
*d groupr.2351
      cc=1
*d groupr.2363,2364
      if (pow.gt.exmin) wtf=wtf+wt6h*exp(pow)
*i groupr.2457
*if sw
      data up,dn/1.001d0,0.999d0/
      data tenth/0.1d0/
      data sigzmx/1.d10/
      data etop/1.d10/
      data zero/0.d0/
*else
      data up,dn/1.001e0,0.999d0/
      data tenth/0.1e0/
      data sigzmx/1.e10/
      data etop/1.e10/
      data zero/0.e0/
*endif
*d groupr.2471
      e=0
*d groupr.2473
      if (enext.gt.up*egn(1))
*d groupr.2493,2494
      if (alpha2.ne.zero) nalph=2
      if (alpha3.ne.zero) nalph=3
*d groupr.2499,2500
      if (elo.lt.tenth) elo=tenth
*d groupr.2505
      e=0
*d groupr.2517
      e=0
*d groupr.2572
      f4=1-f2*log(ejp/ej)
*d groupr.2637
      e=0
*d groupr.2644
      e=0
*d groupr.2666
      sigzj=sigzmx
*d groupr.2675
      e=0
*d groupr.2709,2710
      if (en.gt.dn*etop) enext=etop
*d groupr.2719
      xtot=0
*d groupr.2721
      if (en.lt.zero) go to 340
*d groupr.2726
      fac=0
*d groupr.2737,2738
      if (enext.gt.dn*etop) nett=-nett
*d groupr.2781
      zero=0
      if (e.gt.zero) go to 200
*d groupr.2783,2788
      el=0
      wl=0
      xl=0
      en=0
      wn=0
      xn=0
*d groupr.2792
      if (flag.lt.zero) go to 110
*d groupr.2801
      if (flag.ge.zero) go to 130
*i groupr.2860
*if sw
      data ebig/1.d8/
*else
      data ebig/1.e8/
*endif
*d groupr.2897
      econst=0
*d groupr.2901,2903
      if (mtd.eq.18.or.mtd.eq.19) econst=ebig
      if (mtd.eq.452.or.mtd.eq.455.or.mtd.eq.456) econst=ebig
      if (mtd.eq.102) econst=ebig
*i groupr.2960
      data zero/0.d0/
      data smin/1.d-9/
      data ttest/0.1d0/
*i groupr.2976
      data zero/0.e0/
      data smin/1.e-9/
      data ttest/0.1e0/
*d groupr.2990
      eqw=0
*d groupr.3027
  114 if (sig(1,1).ne.zero) go to 125
*d groupr.3030,3031
      ans(il,iz,1)=1
      ans(il,iz,2)=0
*d groupr.3103,3106
      if (elo.lt.ttest) go to 220
      if (abs(sig(1,1)).gt.smin) go to 220
*d groupr.3109
      ans(il,iz,2)=0
*i groupr.3152
      data big/1.d10/
      data zero/0.d0/
*i groupr.3155
      data big/1.e10/
      data zero/0.e0/
*d groupr.3228,3230
      result(il)=0
      if (ans(il,1,ig2).eq.zero) go to 230
      if (ans(il,1,1).eq.zero) ans(il,1,1)=big
*d groupr.3232,3233
      if (result(il).ne.zero) then
         jjj=int(log10(result(1)))
*d groupr.3262
      if (ans(1,1,1).eq.zero) ans(1,1,1)=big
*d groupr.3264,3265
      if (abs(ans(1,1,j)).lt.smin/1000) ans(1,1,j)=0
      if (ans(1,1,j).eq.zero) go to 260
*d groupr.3292,3293
      result(iz)=0
      if (ans(il,iz,ig2).eq.zero) go to 330
*d groupr.3296
      if (ans(il,iz,locd).eq.zero) ans(il,iz,locd)=big
*d groupr.3300
         jjj=int(log10(result(1)))
*d groupr.3304
         iii=nint(fact*result(iz))+nint(ten**(ndig-11))
*d groupr.3416
      zero=0
      if (e.gt.zero) go to 105
*d groupr.3488
      zero=0
      if (e.gt.zero) go to 200
*d groupr.3584
      zero=0
      if (e.gt.zero) go to 300
*d groupr.3629
      thresh=0
*d groupr.3651
      if (en.lt.zero) return
*d groupr.3670
      sig(iz,il)=1
*d groupr.3750
      if (abs(sigz(i)-a(iunr+5+j)).lt.sigz(i)/100) go to 170
*i groupr.3781
      zero=0
*d groupr.3820
      if (a(locl).lt.zero) iovl=1
*d groupr.3827
      if (lssf.gt.0.and.sinf.ne.zero) sig(is)=sig(is)*s/sinf
*d groupr.3835
      if (a(loc).lt.0..and.a(loc1).lt.zero) iovl=1
*d groupr.3846
      if (lssf.gt.0.and.sinf.ne.zero) sig(is)=sig(is)*s/sinf
*i groupr.3860
*if sw
      data rlo,rhi/0.99d0,10.1d0/
*else
      data rlo,rhi/0.99e0,10.1e0/
*endif
*d groupr.3862
      sig=0
*d groupr.3866,3868
      if (sigz.lt.sigs(ns).and.sm.gt.rhi*sigs(ns)) go to 110
      if (sigz.lt.sigs(2).and.sm.gt.rhi*sigs(2)) go to 110
      if (sigz.ge.sigs(2).and.sm.lt.rlo*sigs(2)) go to 110
*d groupr.3871
      if (abs(sm-st).lt.sm/1000) go to 120
*d groupr.3873
      if (st.gt.rhi*sigs(ns)) go to 120
*d groupr.3877
      if (st.gt.rhi*sigs(2)) go to 120
*d groupr.3880
  140 if (st.lt.rlo*sigs(2)) go to 120
*i groupr.3914
      zero=0
*d groupr.3968
      yld=1
*d groupr.3975
      if (e.eq.zero) return
*d groupr.3988,3989
      if (e.eq.zero) ifirst=0
      if (e.gt.zero) ifirst=ifirst+1
*d groupr.4012,4013
      if (e.eq.zero) call reserv('sed',nwds,ised,a)
      if (e.ne.zero) call findex('sed',ised,a)
*d groupr.4018,4019
      if (e.gt.zero) go to 205
      if (econst.eq.zero) return
*d groupr.4026
      if (jconst.le.2) econst=0
*d groupr.4065
      if (ff(ik,jg).ne.zero) ighi=ig
*d groupr.4080
      if (e.gt.zero) go to 325
*d groupr.4088
      if (econ.ne.zero) go to 305
*d groupr.4099
  315 if (econst.eq.zero) go to 322
*d groupr.4106
      if (jconst.le.2) econst=0
*d groupr.4124
      if (a(ieyl-1+nyl).ne.zero) go to 330
*d groupr.4131
      if (e.gt.zero) go to 330
*d groupr.4141
      ff(il,ig)=0
*d groupr.4144,4145
      if (a(ik-1+igyl).eq.zero) go to 350
      if (a(ik-1+ieyl).eq.zero) go to 340
*d groupr.4149
      if (a(ieyl+ik-1).eq.zero) tempo=tempo*a(ised+ig-1)
*d groupr.4157
      yld=0
*d groupr.4186
      yld=0
*d groupr.4206
      yld=1
*d groupr.4208
      if (e.eq.zero) return
*d groupr.4213,4214
      if (abs(ff(il,ig)).le.small) go to 720
*i groupr.4397
      data up,dn/1.00001d0,0.99999d0/
      data eps/0.02d0/
      data zero/0.d0/
*i groupr.4401
      data up,dn/1.00001e0,0.99999e0/
      data eps/0.02e0/
      data zero/0.e0/
*d groupr.4405
      if (ed.gt.zero) go to 200
*d groupr.4418
      if (mfd.eq.18) zad=0
*d groupr.4552,4553
      c(l)=0
      c(l+1)=0
*d groupr.4670
  199 enext=-1
*d groupr.4677
      ans(l,j)=0
*d groupr.4681
      yld=0
*d groupr.4696
      if (pe.eq.zero) go to 450
*d groupr.4730
  300 ep=0
*d groupr.4763,4764
      if (el.lt.zero) el=ed*awr/(awr+1)-el
      if (eh.lt.zero) eh=ed*awr/(awr+1)-eh
*d groupr.4769
      if (j.eq.1) e1=0
*d groupr.4796,4798
  700 if (ed.le.zero) return
      if (yld.eq.zero) return
      test=0
*d groupr.4804
      if (q.gt.zero) go to 710
*d groupr.4807
      ans(l,i)=0
*d groupr.4810,4811
      if (test.eq.zero) return
*d groupr.4813
     1  .and.(ed.lt.up*elo.or.ed.gt.dn*ehi))
*d groupr.4851
      c(jnow)=0
*d groupr.4856,4857
      sum=0
      x(2)=0
*d groupr.4966
      zero=0
      if (ep.gt.zero) go to 200
*d groupr.4973
      elmax=0
*d groupr.4975
      epmax=0
*d groupr.5001
      term(l)=0
*d groupr.5007
      term(l)=0
*d groupr.5048
      if (un.lt.umin+tiny/10) un=umin
*d groupr.5141
      if (abs(term(l)).lt.test) term(l)=0
*i groupr.5175
      data tomev/1.d-6/
      data half/0.5d0/
      data bach1,bach2/0.0236d0,0.0039d0/
      data bach3,bach4/92.d0,90.d0/
      data zero/0.d0/
*i groupr.5181
      data tomev/1.e-6/
      data half/0.5e0/
      data bach1,bach2/0.0236e0,0.0039e0/
      data bach3,bach4/92.e0,90.e0/
      data zero/0.e0/
*d groupr.5186
      if (ep.gt.zero) go to 120
*d groupr.5194
      if (cnow(lnow).le.zero) lnow=lnow+ncnow
*d groupr.5200,5201
      efirst=0
      f6ddx=0
*d groupr.5223
      t=0
*d groupr.5229
      t=0
*d groupr.5240
      if (t.lt.zero) t=0
*d groupr.5245,5246
      s=0
      r=0
*d groupr.5257
      t=half
*d groupr.5259,5262
      aa=bach1+bach2*l*(l+1)
      ss=l*(l+1)
      bb=bach3-bach4/sqrt(ss)
      sa=1/(1+exp(aa*(bb-ep*tomev)))
      tt=(2*l+1)*sa*p(l+1)/2
*d groupr.5272
      if (t.lt.zero) t=0
*d groupr.5282
      tii=0
*d groupr.5289
      tjj=0
*i groupr.5326
*if sw
      data tomev/1.d-6/
      data half/0.5d0/
      data bach1,bach2/0.0236d0,0.0039d0/
      data bach3,bach4/92.d0,90.d0/
      data zero/0.d0/
*else
      data tomev/1.e-6/
      data half/0.5e0/
      data bach1,bach2/0.0236e0,0.0039e0/
      data bach3,bach4/92.e0,90.e0/
      data zero/0.e0/
*endif
*d groupr.5330
      if (i.gt.0) go to 120
*d groupr.5342
      f6dis=0
*d groupr.5352
      t=0
*d groupr.5358
      if (t.lt.zero) t=0
*d groupr.5369
      t=half
*d groupr.5371,5374
      aa=bach1+bach2*l*(l+1)
      ss=l*(l+1)
      bb=bach3-bach4/sqrt(ss)
      sa=1/(1+exp(aa*(bb-ep*tomev)))
      tt=(2*l+1)*sa*p(l+1)/2
*d groupr.5384
      if (t.lt.zero) t=0
*d groupr.5391
      t=0
*d groupr.5452
      za=int(iza/1000)
*d groupr.5454
      zc=za+int(iza1/1000)
*d groupr.5456
      zb=zc-int(iza2/1000)
*d groupr.5487,5488
      fa=1
      if (iza1.eq.2004) fa=0
*i groupr.5528
      zero=0
*d groupr.5532
      c(jnow)=0
*d groupr.5541
      epnext=0
*d groupr.5566
      term(l)=0
*d groupr.5582
      sum=0
*d groupr.5587
      if (term(1).ne.zero) fact=sum/term(1)
*d groupr.5599
      test=tol*abs(tmid)+small/100
      if (abs(tmid-c(lnow-nl-1+i)).gt.test) go to 550
*i groupr.5646
      data emin/1.d-5/
*i groupr.5655
      data emin/1.e-5/
*d groupr.5659
      zero=0
      if (ep.gt.zero) go to 200
*d groupr.5675
      if (clo(llo).le.zero) llo=llo+nclo
*d groupr.5682
      if (chi(lhi).le.zero) lhi=lhi+nchi
*d groupr.5696
      term(l)=0
*d groupr.5699
  150 epnext=emin
*d groupr.5735
      term1(l)=0
*d groupr.5741
  255 term2(l)=0
*d groupr.5751
  265 term1(l)=0
*d groupr.5757
  266 term2(l)=0
*d groupr.5774
      term1(l)=0
*d groupr.5778
      t0=0
*d groupr.5799
      term2(l)=0
*d groupr.5803
      t0=0
*i groupr.5915
      zero=0
*d groupr.5946
      ast=0
*d groupr.5961
      ff(il,ng)=0
*d groupr.5980
      nqp=int(npo*2*b)
*d groupr.5987
      flt(il)=0
*d groupr.6009
      prob=0
*d groupr.6016
      i2s=nint(2*spi)
*d groupr.6019,6020
      zt=int(izat/1000)
      zi=int(izap/1000)
*d groupr.6023
      sigc=0
*d groupr.6050
  464 ff(1,1)=1
*d groupr.6069
  600 if (e.gt.zero) go to 610
*d groupr.6080
      if (thresh.gt.zero) go to 630
*d groupr.6072,6074
      ecl=0
      ech=0
      ecn=0
*d groupr.6077
      ecl=emax
*d groupr.6090
      if (disc.lt.zero) go to 620
*d groupr.6096
      ech=emax
*d groupr.6099
      if (thresh.gt.zero) go to 650
*d groupr.6102
      test=test/10
*d groupr.6111
      if (disc.lt.zero) go to 640
*d groupr.6156
      zero=0
      if (e.gt.zero) go to 200
*d groupr.6257
  320 fle(i)=0
*d groupr.6324
      zero=0
      if (e.gt.zero) go to 200
*d groupr.6358
      ehi=0
*d groupr.6373
      clast=0
*d groupr.6382
      if ((cnow-clast).le.(eps*clast+eps)) go to 160
*d groupr.6439
      aed(il,ig)=0
*d groupr.6445
      fl(il)=0
*d groupr.6447
      eg=0
*d groupr.6451,6453
      egp=0
      egp1=0
      egp2=0
*d groupr.6479
      if (ek2.gt.up*eg2) go to 385
*d groupr.6515
      aed(il,i)=aed(il,i)+(fi(il)+fl(il))*(ei-eg)/2
*d groupr.6531
      aed(il,ig)=0
*d groupr.6546
      aed(il,i)=0
*d groupr.6613
  630 if (aed(1,ig).eq.zero) go to 650
*d groupr.6618
  650 aed(1,ig)=1
*d groupr.6641
      fl(il)=0
*i groupr.6684
      data up/1.1d0/
*i groupr.6687
      data up/1.1e0/
*d groupr.6690
      zero=0
      if (ed.gt.zero) go to 200
*d groupr.6779
  205 gfl(1,j)=1
*d groupr.6782
      gfl(i,j)=0
*d groupr.6803
      if (im.eq.np.and.ed.gt.up*ehi) call error('getgfl',
*d groupr.6949
  220 fl(il)=0
*d groupr.6961
  240 fl(il)=0
*d groupr.6968
      fl(il)=0
*d groupr.7004
  356 y=0
*d groupr.7013,7014
      zt=int(izat/1000)
      zi=int(izap/1000)
*d groupr.7100
      zero=0
      if (ed.gt.zero) go to 200
*d groupr.7162
      e=0
*d groupr.7168
      if (econst.eq.zero) go to 190
*d groupr.7191
      if (a(l).eq.zero) go to 230
*i groupr.7239
*if sw
      data etop/20.d6/
      data zero/0.d0/
*else
      data etop/20.e6/
      data zero/0.e0/
*endif
*d groupr.7287,7288
  312 za2=-1
      if (mth.eq.2.or.(mth.ge.50.and.mth.lt.91)) za2=1
*d groupr.7425
      a(ie)=0
*d groupr.7427,7428
      a(i+ir)=0
      a(i+iaa)=0
*d groupr.7431
      a(imax*(i-1)+i-1+ir)=1
*d groupr.7458,7459
      a(iaa+(k-1)*imax+j-1)=0
      a(ir+(k-1)*imax+j-1)=0
*d groupr.7474
      ysum=0
*d groupr.7482
      if (y.eq.zero) go to 245
*d groupr.7490
      a(ie)=0
*d groupr.7525,7526
      a(iscr)=0
      a(iscr+1)=0
*d groupr.7534
      a(iscr+10)=etop
*d groupr.7542
      a(iscr+2)=0
*d groupr.7550
      a(iscr+10)=etop
*d groupr.7672
      e1=0
*d groupr.7688
      a(iscr+5+2*im)=sig(1,1)*(1-w)*a(iscr+5+2*im)
*i groupr.7897
      data ebig/1.d8/
      data tenth/0.1d0/
      data zero/0.d0/
*i groupr.7903
      data ebig/1.e8/
      data tenth/0.1e0/
      data zero/0.e0/
*d groupr.7908
      if (ed.gt.zero) go to 200
*d groupr.8002
      if (ig.eq.1) e1=0
*d groupr.8004
      if (ig.eq.ng) e2=ebig
*d groupr.8043
      econ=0
*d groupr.8055
      sed(ik,ig)=0
*d groupr.8071
      if (pe.eq.zero) go to 290
*d groupr.8082
      if (ig.eq.1) e1=0
*d groupr.8084
      if (ig.eq.ng) e2=ebig
*d groupr.8121
      sigup=0
*d groupr.8129
      sed(ikt,ig)=0
*d groupr.8131
  280 test=tenth
*d groupr.8155
      sed(1,ig)=0
*i groupr.8191
      data etest/1.1d-5/
*i groupr.8199
      data etest/1.1e-5/
*d groupr.8206,8207
      if (ep1.lt.etest) new=1
*d groupr.8212
      g=0
*d groupr.8219
      zero=0
      if (epl.gt.de*(1+small).or.de.le.zero) return
*d groupr.8273
      if (top.lt.small) top=0
*d groupr.8286
      if (top.lt.small) top=0
*d groupr.8325
  200 if (top.lt.small) top=0
*d groupr.8468,8469
      cm(1)=1
      xk=1
*i groupr.8416
      data epsrel,epsabs/1.d-9,1.d-30/
      data zero/0.d0/
*i groupr.8423
      data epsrel,epsabs/1.e-9,1.e-30/
      data zero/0.e0/
*d groupr.8432
      if (aa.lt.zero) fa(1)=-fa(1)
*d groupr.8434
      if (bb.lt.zero) fb(1)=-fb(1)
*d groupr.8450
  200 sgn=1
*d groupr.8493
      beta=0
*d groupr.8509
      test=epsabs+epsrel*abs(s)
*d groupr.8554
      zero=0
      if (ep.gt.zero) go to 200
*d groupr.8563
      if (eimax.le.zero) eimax=100
*d groupr.8572
      f6psp=0
*d groupr.8578
  200 s=0
*ident up53
*/ leapr -- 2dec98
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*/   in addition, any conversion of real numbers to integers is 
*/   explicitly flagged for truncation or for nearest integer. 
*/   an error in mixed-mode conversion was fixed that affects
*/   the statistical factor for cold hydrogen or deuterium.
*d leapr.229,231
      b7=0
      aws=0
      sps=0
*d leapr.318
      call trans(itemp,a(issm),nalpha,nbeta,ntempr,a(ibeta),
*d leapr.373
      emax=5
*i leapr.420
      data tiny/1.d-30/
*i leapr.423
      data tiny/1.e-30/
*d leapr.431
      sc=1
*d leapr.447,448
      if (add.lt.tiny) add=0
*d leapr.473,474
      if (add.lt.tiny) add=0
*d leapr.478
      if (add.gt.ssm(k,j,itemp)/1000.and.j.lt.maxt(k)) maxt(k)=j
*i leapr.511
      bel=0
      ff1l=0
      ff2l=0
      sum0=0
      sum1=0
*d leapr.542.543
      sum0=0
      sum1=0
*d leapr.594
      tau=1
      tau=tau/2
*d leapr.658.662
      v=1
      an=1-2*mod(n,2)
      be=0
      fs=0
      w=1
*d leapr.685
      i=int(be/delta)
*d leapr.692
  120 terpt=0
*i leapr.707
*if sw
      data tiny/1.d-30/
*else
      data tiny/1.e-30/
*endif
*d leapr.709
      ckk=0
*d leapr.711
      tnext(k)=0
*d leapr.715
      f1=0
*d leapr.719
      f2=0
*d leapr.731,732
      if (tnext(k).lt.tiny) tnext(k)=0
*i leapr.767
      data tiny/1.d-30/
*i leapr.770
      data tiny/1.e-30/
*d leapr.777
      sc=1
*d leapr.820
      s=0
*d leapr.836,837
      if (s.lt.tiny) s=0
*i leapr.851
      sum0=0
      sum1=0
      ff1l=0
      ff2l=0
      bel=0
*d leapr.867,868
      sum0=0
      sum1=0
*d leapr.929
      be=0
*d leapr.935
         c5=0
*d leapr.957
      be=0
*d leapr.979,980
      check0=0
      check1=0
*i leapr.1006
*if sw
      data slim/-225.d0/
*else
      data slim/-225.e0/
*endif
*d leapr.1008
      i=int(be/delta)
*d leapr.1014
         st=slim
*d leapr.1019
         stp=slim
*d leapr.1024,1026
      terps=0
      if (stt.gt.slim) terps=exp(stt) 
*d leapr.1028
  120 terps=0
*i leapr.1042
*if sw
      data shade/1.00001d0/
      data slim/-225.d0/
*else
      data shade/1.00001e0/
      data slim/-225.e0/
*endif
*d leapr.1053
      if (j.eq.nbeta.and.b.lt.shade*beta(j)) go to 130
*d leapr.1065
         st=slim
*d leapr.1070
         stm=slim
*d leapr.1077,1079
      sb(i)=0
      if (arg.gt.slim) sb(i)=exp(arg)
*d leapr.1081
  135 sb(i)=0
*d leapr.1127
      test=1
*i leapr.1170
      data small/1.d-5/
      data vsmall/1.d-10/
      data tiny/1.d-20/
*i leapr.1172
      data small/1.e-5/
      data vsmall/1.e-10/
      data tiny/1.e-20/
*d leapr.1180
      sc=1
*d leapr.1184
      dwt=0
*d leapr.1189
      tsave=0
*d leapr.1228
      sexpb(j)=0
*d leapr.1232,1233
      ben(1)=0
      wtn(1)=1
*d leapr.1248,1249
      if (besn.gt.0..and.wtsn.lt.small) go to 160
*d leapr.1263,1264
      if (wtsn.lt.small) go to 170
*d leapr.1279,1280
      if (wtsn.lt.small) go to 185
*d leapr.1317,1318
      if (wts(i).lt.small) go to 207
*d leapr.1326,1327
      sumn=0
      sumr=0
*d leapr.1344,1345
      if (add.lt.tiny) go to 225
*d leapr.1354,1355
      if (dwf.lt.vsmall) go to 265
*d leapr.1359
      db=1000
*d leapr.1370,1371
      if (add.lt.tiny) go to 260
*i leapr.1395
      sum0=0
      sum1=0
      ff1l=0
      ff2l=0
      bel=0
*d leapr.1411,1412
      sum0=0
      sum1=0
*d leapr.1414
      if (twt.eq.0.) sum0=sum0/(1-exp(-al*dwpix(itemp)))
*i leapr.1457
      data big/1.d10/
      data tiny/1.d-30/
*i leapr.1469
      data big/1.e10/
      data tiny/1.e-30/
*d leapr.1495,1496
      bn(imax)=0
      bn(imax-1)=1
*d leapr.1500
      if (bn(i).lt.big) go to 180
*d leapr.1502
      bn(j)=bn(j)/big
*d leapr.1516
      bplus(i)=0
*d leapr.1518,1520
      if (bplus(i).lt.tiny) bplus(i)=0
      bminus(i)=0
*d leapr.1523,1524
      if (bminus(i).lt.tiny) bminus(i)=0
*d leapr.1530
      bplus(i)=0
*d leapr.1533,1535
      if (bplus(i).lt.tiny) bplus(i)=0
      bminus(i)=0
*d leapr.1538,1539
      if (bminus(i).lt.tiny) bminus(i)=0
*i leapr.1553
*if sw
      data small/1.e-9/
*else
      data small/1.e-9/
*endif
*d leapr.1560,1562
      if (beta(1).gt.small) go to 110
      bex(nbeta)=0
*i leapr.1589
*if sw
      data small/1.e-9/
*else
      data small/1.e-9/
*endif
*d leapr.1596,1597
      if (beta(1).gt.small) go to 110
*i leapr.1621
      data slim/-225.d0/
*i leapr.1623
      data slim/-225.e0/
*d leapr.1629
      sv=0
*d leapr.1654
         ss1=slim
*d leapr.1659
         ss3=slim
*d leapr.1664,1666
      sv=0
      if (ex.gt.slim) sv=exp(ex)
*d leapr.1688
*i leapr.1700
      data cct/0.01381d0/
      data deh/0.0147d0/
      data ded/0.0074d0/
      data amassh/3.465d0/
      data amassd/6.69d0/
      data disth/1.5d0/
      data distd/0.71d0/
      data sampch/0.356d0/
      data sampcd/0.668d0/
      data sampih/2.526d0/
      data sampid/0.403d0/
      data small/1.d-6/
*i leapr.1707
      data cct/0.01381e0/
      data deh/0.0147e0/
      data ded/0.0074e0/
      data amassh/3.465e0/
      data amassd/6.69e0/
      data disth/1.5e0/
      data distd/0.71e0/
      data sampch/0.356e0/
      data sampcd/0.668e0/
      data sampih/2.526e0/
      data sampid/0.403e0/
      data small/1.e-6/
*d leapr.1721,1722
      t=cct*abs(temp)
      sc=1
*d leapr.1725,1726
      de=deh
      if (law.eq.3) de=ded
*d leapr.1728,1735
      amassm=amassh
      if (law.gt.3) amassm=amassd
      bp=hbar/2*sqrt(2/disth/ev/pmass)
      if (law.gt.3) bp=hbar/2*sqrt(2/distd/ev/dmass)
      sampc=sampch
      if (law.gt.3) sampc=sampcd
      sampi=sampih
      if (law.gt.3) sampi=sampid
*d leapr.1745
      waven=sqrt(amassm*t*al)/hbar
*d leapr.1750
      if (nokap.eq.1) sk=1
*d leapr.1757,1764
  120 if (law.eq.2) swe=sampi**2/3
      if (law.eq.2) swo=sk*sampc**2+2*sampi**2/3
      if (law.eq.3) swe=sk*sampc**2
      if (law.eq.3) swo=sampi**2
      if (law.eq.4) swe=sk*sampc**2+5*sampi**2/8
      if (law.eq.4) swo=3*sampi**2/8
      if (law.eq.5) swe=3*sampi**2/4
      if (law.eq.5) swo=sk*sampc**2+sampi**2/4
*d leapr.1789,1790
      sn=0
      total=0
*d leapr.1804
      snlg=0
*d leapr.1807,1810
      betap=(-j*(j+1)+jp*(jp+1))*x/2
      tmp=(2*jp+1)*pj*swe*4*sumh(j,jp,y)
      if (jj.eq.1.and.tmp.ge.small) then
*d leapr.1826
      snlk=0
*d leapr.1829
      betap=(-j*(j+1)+jp*(jp+1))*x/2
*d leapr.1831,1832
      if (jj.eq.1.and.tmp.ge.small) then
*d leapr.1867
      up=0
*i leapr.1874
      sum0=0
      bel=0
      ff1l=0
      ff2l=0
*d leapr.1888
      sum0=0
*d leapr.1905
*if sw
      data half/0.5d0/
*else
      data half/0.5e0/
*endif
c
      yy=half*j*(j+1)
*d leapr.1907
      b=0
*d leapr.1911
      yy=half*k*(k+1)
*d leapr.1934
  120 sum1=0
*d leapr.1971,1972
      s=0
      fact=1
*d leapr.1974
      zi=i
*d leapr.1979,1980
      s=0
      fact=1
*d leapr.1982
      zi=i
*d leapr.1987,1988
      s=0
      fact=1
*d leapr.1990
      zi=i
*d leapr.1995,1996
      s=0
      fact=1
*d leapr.1998
      zi=i
*d leapr.2003,2004
      s=0
      b1=1
*d leapr.2006
      zi=i
*d leapr.2010,2011
      s=0
      b2=1
*d leapr.2013
      zi=i
*d leapr.2017,2018
      s=0
      b3=1
*d leapr.2020
      zi=i
*d leapr.2024,2025
      s=0
      b4=1
*d leapr.2027
      zi=i
*d leapr.2031,2033
      rat=2*nn+1
      rat=rat/(jj+ll+nn+1)
      iwign=(jj+ll-nn)/2
      wign=(-1)**iwign
      wign=wign*sqrt(rat)*b1/a1*a2/b2*a3/b3*a4/b4
*d leapr.2035
  190 wign=0
*i leapr.2048
*if sw
      data huge/1.d25/
      data small/2.d-38/
      data break1/3.d4/
      data break2/7.d-4/
      data break3/0.2d0/
      data one/1.d0/
      data ten/1.d1/
      data hund/1.d2/
*else
      data huge/1.e25/
      data small/2.e-38/
      data break1/3.e4/
      data break2/7.e-4/
      data break3/0.2e0/
      data one/1.e0/
      data ten/1.e1/
      data hund/1.e2/
*endif
*d leapr.2051,2052
      if (x.gt.break1) go to 240
*d leapr.2054,2056
      if (x.gt.break2) go to 120
      w=1
*d leapr.2060,2061
      t1=3
      t2=1
*d leapr.2064
      t1=t1+2
*d leapr.2069,2070
  120 if (x.ge.break3) go to 130
*d leapr.2072
      w=1-y*(1-y/20)/6
*d leapr.2079,2080
  160 if (x.lt.hund) go to 170
      l=int(x/50+18) 
*d leapr.2082,2083
  170 if (x.lt.ten) go to 180
      l=int(x/10+10) 
*d leapr.2085,2086
  180 if (x.le.one) go to 190
      l=int(x/2+5) 
*d leapr.2093,2095
      z=1/x
      t3=0
      t2=small
*d leapr.2098
      t1=(2*k+3)*z*t2-t3
*d leapr.2101,2105
  210 if (abs(t1).lt.huge) go to 220
      t1=t1/huge
      t2=t2/huge
      sj=sj/huge
*d leapr.2118
  260 bessel=0
*d leapr.2133
      i=int(be/delta)
*d leapr.2140
  120 terpk=1
*i leapr.2184
      data toler/1.d-6/
      data eps/.05d0/
*i leapr.2194
      data toler/1.e-6/
      data eps/.05/
*d leapr.2204
      twopis=(2*pi)**2
*d leapr.2207,2209
      tsqx=econ/20
*d leapr.2233
      wint=0
*d leapr.2241
      i1m=int(a*sqrt(phi))
*d leapr.2246
      i2m=int((l1+sqrt(3*(a*a*phi-l1*l1)))/2)
*d leapr.2252
      if (x.gt.0.) i3m=int(c*sqrt(x))
*d leapr.2260
      if (l1.eq.0.and.l2.eq.0) w2=1
      if (l1.eq.0.and.l2.eq.0) w2=w2/2
*d leapr.2275
      if (tsq.lt.b(ifl+2*i-2).or.tsq.ge.(1+eps)*b(ifl+2*i-2))
*d leapr.2296
      if (tsq.lt.b(ifl+2*i-2).or.tsq.ge.(1+eps)*b(ifl+2*i-2))
*d leapr.2328
      bel=-1
*d leapr.2333,2334
      if (be-bel.ge.toler) go to 504
*i leapr.2360
      data c1,c2,c3/7.54d0,4.24d0,11.31d0/
*i leapr.2362
      data c1,c2,c3/7.54,4.24,11.31/
*d leapr.2379,2380
  200 formf=(1+cos(2*pi*(2*l1+4*l2+3*l3)/6))*(c1+c2+c3*cos(3*pi*l3/4))
*i leapr.2409
*if sw
      data small/1.d-9/
      data tiny/-999.d0/
      data smin/2.d-38/
      data tol/0.9d-7/
      data up/1.01d0/
      data therm/.0253d0/
*else
      data small/1.e-9/
      data tiny/-999.e0/
      data smin/2.e-38/
      data tol/0.9e-7/
      data up/1.01e0/
      data therm/.0253e0/
*endif
*d leapr.2425,2426
*d leapr.2433,2434
      sb=spr*((1+awr)/awr)**2
      if (aws.ne.0.) sbs=sps*((1+aws)/aws)**2
*d leapr.2493,2497
      scr(1)=0
      scr(2)=0
      scr(3)=0
      scr(4)=0
      scr(5)=0
*d leapr.2500,2503
      scr(1)=1
      scr(2)=0
      scr(3)=0
      scr(4)=0
*d leapr.2507,2508
      scr(1)=0
      scr(2)=0
*d leapr.2532,2533
  160 scr(1)=0
      scr(2)=0
*d leapr.2540,2541
      scr(ii+1)=0
      scr(ii+2)=0
*d leapr.2551,2552
  165 scr(ii+1)=0
      scr(ii+2)=0
*d leapr.2578
      scr(2)=0
*d leapr.2611,2612
      sum=0
      suml=0
*d leapr.2618
      if (sum-suml.le.tol*sum) go to 190
*d leapr.2627
      scr(2)=0
*d leapr.2637
      sum=0
*d leapr.2639
      if (nss.gt.0.and.b7.eq.0.) w=(dbw(i)+dbw1(i))/2
*d leapr.2658
      scr(2)=0
*d leapr.2666
      sum=0
*d leapr.2668
      if (nss.gt.0.and.b7.eq.0.) w=(dbw(i)+dbw1(i))/2
*d leapr.2697,2698
      scr(1)=0
      scr(2)=0
*d leapr.2708
      scr(11)=0
*d leapr.2714,2715
      scr(16)=0
      scr(17)=0
*d leapr.2718,2719
      scr(1)=0
      scr(2)=0
*d leapr.2732
      sc=1
*d leapr.2751,2752
            if (scr(8+2*j).ge.small) then
*d leapr.2758
            scr(8+2*j)=tiny
*d leapr.2768,2769
               if (scr(8+2*j).ge.small) then
*d leapr.2775
               scr(8+2*j)=tiny
*d leapr.2784,2785
               if (scr(8+2*j).ge.small) then
*d leapr.2791
               scr(8+2*j)=tiny
*d leapr.2801,2802
           if (scr(8+2*j).ge.small) then
*d leapr.2808
            scr(8+2*j)=tiny
*d leapr.2818,2819
               if (scr(8+2*j).ge.small) then
*d leapr.2825
               scr(8+2*j)=tiny
*d leapr.2834,2835
               if (scr(8+2*j).ge.small) then
*d leapr.2841
               scr(8+2*j)=tiny
*d leapr.2849
      if (ilog.eq.0.and.scr(8+2*j).lt.smin) scr(8+2*j)=0
*d leapr.2866,2867
            if (scr(6+j).ge.small) then
*d leapr.2873
            scr(6+j)=0
*d leapr.2883,2884
               if (scr(6+j).ge.small) then
*d leapr.2890
               scr(6+j)=tiny
*d leapr.2899,2900
               if (scr(6+j).ge.small) then
*d leapr.2906
               scr(6+j)=tiny
*d leapr.2916,2917
            if (scr(6+j).ge.small) then
*d leapr.2923
            scr(6+j)=tiny
*d leapr.2933,2934
               if (scr(6+j).ge.small) then
*d leapr.2940
               scr(6+j)=tiny
*d leapr.2949,2950
               if (scr(6+j).ge.small) then
*d leapr.2956
               scr(6+j)=tiny
*d leapr.2964
      if (ilog.eq.0.and.scr(6+j).lt.smin) scr(6+j)=0
*d leapr.2970,2971
      scr(1)=0
      scr(2)=0
*d leapr.2984
  290 scr(2*i+7)=sigfig(up*scr(2*i+5),7,0)
*d leapr.2988,2989
  300 scr(1)=0
      scr(2)=0
*d leapr.3002
  305 scr(2*i+7)=sigfig(up*scr(2*i+5),7,0)
*ident up54
*/ gaspr -- 2dec98
*/   increase the maximum energy grid size.
*d gaspr.28
      dimension egas(50000),sgas(5,50000)
*d gaspr.33
*if sw
      data emax/1.d10/
*else
      data emax/1.e10/
*endif
      maxg=50000
*/   fix various issues raised by ftnchek.  these changes make
*/   sure that there are no 4-byte constants in statements that
*/   might lead to a loss in precision.  we have even eliminated
*/   occurences of "0." for completeness, even though this value
*/   probably doesn't lead to problems with known compilers.
*/   in addition, any conversion of real numbers to integers is 
*/   explicitly flagged for truncation or for nearest integer. 
*/   also, add more fortran block structures, including do blocks.
*d gaspr.74
      zain=1
*d gaspr.78
      zain=int(nsub/10)
*d gaspr.81,84
      do while (nb.ne.0)
         call moreio(nendf,0,0,a(1),nb,nw)
      enddo
      if (iverf.ne.4) nx=n2h
*d gaspr.113,117
      do while (nb.ne.0)
         call moreio(nendf,0,0,a(ll),nb,nw)
         ll=ll+nw
      enddo
      if (izap.eq.1001) then
*d gaspr.119,121
         do i=1,nw
            six(lsix+i-1)=a(i)
         enddo
*d gaspr.125.127
         do i=1,nw
            six(lsix+i-1)=a(i)
         enddo
*d gaspr.131,133
         do i=1,nw
            six(lsix+i-1)=a(i)
         enddo
*d gaspr.137,139
         do i=1,nw
            six(lsix+i-1)=a(i)
         enddo
*d gaspr.143,145
         do i=1,nw
            six(lsix+i-1)=a(i)
         enddo
*d gaspr.157,162
      do ie=1,ne
         call listio(nendf,0,0,a,nb,nw)
         do while (nb.ne.0)
            call moreio(nendf,0,0,a,nb,nw)
         enddo
      enddo
*d gaspr.166.175
      do ie=1,ne
         call tab2io(nendf,0,0,a,nb,nw)
         nmu=n2h
         do imu=1,nmu
            call tab1io(nendf,0,0,a,nb,nw)
            do while (nb.ne.0)
               call moreio(nendf,0,0,a,nb,nw)
            enddo
         enddo
      enddo
*d gaspr.205
      thrg=emax
*d gaspr.210
      en=0
*d gaspr.302
      en=0
*d gaspr.312
      if (enext.lt.emax) go to 260
*d gaspr.321
      sgas(j,i)=0
*d gaspr.338
      en=0
*d gaspr.342,512
      y203=0
      y204=0
      y205=0
      y206=0
      y207=0
      if (mth.eq.5.and.mf6mt5.eq.1) then
         if (l203.gt.0) y203=111
         if (l204.gt.0) y204=111
         if (l205.gt.0) y205=111
         if (l206.gt.0) y206=111
         if (l207.gt.0) y207=111
      else if (mth.eq.11) then
         izr=izr-1004
         y204=1
      else if (mth.eq.16) then
         izr=izr-2
      else if (mth.eq.17) then
         izr=izr-3
      else if (mth.eq.22) then
         izr=izr-2005
         y207=1
      else if (mth.eq.23) then
         izr=izr-6013
         y207=3
      else if (mth.eq.24) then
         izr=izr-2006
         y207=1
      else if (mth.eq.25) then
         izr=izr-2007
         y207=1
      else if (mth.eq.28) then
         izr=izr-1002
         y203=1
      else if (mth.eq.29) then
         izr=izr-4009
         y207=2
      else if (mth.eq.30) then
         izr=izr-4010
         y207=2
      else if (mth.eq.32) then
         izr=izr-1003
         y204=1
      else if (mth.eq.33) then
         izr=izr-1004
         y205=1
      else if (mth.eq.34) then
         izr=izr-2004
         y206=1
      else if (mth.eq.35) then
         izr=izr-5011
         y204=1
         y207=2
      else if (mth.eq.36) then
         izr=izr-5012
         y205=1
         y207=2
      else if (mth.eq.37) then
         izr=izr-4
      else if (mth.eq.41) then
         izr=izr-1003
         y203=1
      else if (mth.eq.42) then
         izr=izr-1004
         y203=1
      else if (mth.eq.44) then
         izr=izr-2003
         y203=2
      else if (mth.eq.45) then
         izr=izr-3006
         y203=1
         y207=1
      else if (mth.ge.51.and.mth.le.91) then
         izr=izr-1
         if (lr.eq.22) then
            izr=izr-2004
            y207=1
         else if (lr.eq.23) then
            izr=izr-6012
            y207=3
         else if (lr.eq.24) then
            izr=izr-2005
            y207=1
         else if (lr.eq.25) then
            izr=izr-2006
            y207=1
         else if (lr.eq.28) then
            izr=izr-1001
            y203=1
         else if (lr.eq.29) then
            izr=izr-4008
            y207=2
         else if (lr.eq.30) then
            izr=izr-4009
            y207=2
         else if (lr.eq.32) then
            izr=izr-1002
            y204=1
         else if (lr.eq.33) then
            izr=izr-1003
            y205=1
         else if (lr.eq.34) then
            izr=izr-2003
            y206=1
         else if (lr.eq.35) then
            izr=izr-5010
            y204=1
            y207=2
         else if (lr.eq.36) then
            izr=izr-5011
            y205=1
            y207=2
         else if (lr.eq.39) then
            izr=izr
         else if (lr.eq.40) then
            izr=izr
         endif
      else if (mth.eq.103) then
         izr=izr-1001
         y203=1
      else if (mth.eq.104) then
         izr=izr-1002
         y204=1
      else if (mth.eq.105) then
         izr=izr-1003
         y205=1
      else if (mth.eq.106) then
         izr=izr-2003
         y206=1
      else if (mth.eq.107) then
         izr=izr-2004
         y207=1
      else if (mth.eq.108) then
         izr=izr-4008
         y207=2
      else if (mth.eq.109) then
         izr=izr-6012
         y207=3
      else if (mth.eq.111) then
         izr=izr-2002
         y203=2
      else if (mth.eq.112) then
         izr=izr-3005
         y203=1
         y207=1
      else if (mth.eq.113) then
         izr=izr-5011
         y205=1
         y207=2
      else if (mth.eq.114) then
         izr=izr-5010
         y204=1
         y207=2
      else if (mth.eq.115) then
         izr=izr-2003
         y203=1
         y204=1
      else if (mth.eq.116) then
         izr=izr-2004
         y203=1
         y205=1
      else if (mth.eq.117) then
         izr=izr-3007
         y204=1
         y207=1
      endif
      if (izr.eq.1001) y203=y203+1
      if (izr.eq.1002) y204=y204+1
      if (izr.eq.1003) y205=y205+1
      if (izr.eq.2003) y206=y206+1
      if (izr.eq.2004) y207=y207+1
      if (izr.eq.4008) y207=y207+2
*d gaspr.596,609
      if (iverf.ge.5) then
         call contio(nscr1,0,0,a(k),nb,nw)
         k=k+nw
      endif
      if (iverf.ge.6) then
         call contio(nscr1,0,0,a(k),nb,nw)
         k=k+nw
      endif
      call hdatio(nscr1,0,0,a(k),nb,nw)
      k=k+nw
      if (iverf.ne.4) nx=n2h
      do while (nb.ne.0)
         call moreio(nscr1,0,0,a(k),nb,nw)
         k=k+nw
      enddo
      nw=nx
*d gaspr.631,671
  401 if (n203.ne.0) then
         a(j+1)=0
         a(j+2)=0
         a(j+3)=3
         a(j+4)=203
         a(j+5)=int((n203+2)/3)
         a(j+6)=1
         j=j+6
      endif
      if (n204.ne.0) then
         a(j+1)=0
         a(j+2)=0
         a(j+3)=3
         a(j+4)=204
         a(j+5)=int((n204+2)/3)
         a(j+6)=1
         j=j+6
      endif
      if (n205.ne.0) then
         a(j+1)=0
         a(j+2)=0
         a(j+3)=3
         a(j+4)=205
         a(j+5)=int((n205+2)/3)
         a(j+6)=1
         j=j+6
      endif
      if (n206.eq.0) then
         a(j+1)=0
         a(j+2)=0
         a(j+3)=3
         a(j+4)=206
         a(j+5)=int((n206+2)/3)
         a(j+6)=1
         j=j+6
      endif
      if (n207.ne.0) then
         a(j+1)=0
         a(j+2)=0
         a(j+3)=3
         a(j+4)=207
         a(j+5)=int((n207+2)/3)
         a(j+6)=1
         j=j+6
      endif
      nw=6
*d gaspr.676,688
      if (iverf.ge.5) then
         call contio(0,noutp,0,a(k),nb,nw)
         k=k+nw
      endif
      if (iverf.ge.6) then
         call contio(0,noutp,0,a(k),nb,nw)
         k=k+nw
      endif
      if (iverf.ne.4) a(k+5)=nx+nsec-nold
      call hdatio(0,noutp,0,a(k),nb,nw)
      k=k+nw
      do while (nb.ne.0)
         call moreio(0,noutp,0,a(k),nb,nw)
         k=k+nw
      enddo
      nw=nx+nsec-nold
*d gaspr.721,724
      a(3)=0
      a(4)=0
      a(5)=0
      a(6)=0
*d gaspr.728,731
      a(1)=0
      a(2)=0
      a(3)=0
      a(4)=0
*ident up55
*/ reconr -- 13jan99
*/   fix a typographical error in up47
*d up47.12
  145 est=estp*(a(ix+i-1)-res(1))
*/ reconr -- 13jan99
*/   make trange incommensurable with normal grid energies
*/   to help in getting the same grid on different machines
*d up10.9
      trange=.4999d0
*d up10.11
      trange=.4999
*d up10.17
      data trange/.4999d0/
*d up10.19
      data trange/.4999/
*/ reconr -- 14jan99
*/   fix errors made in up47 while eliminating real*4
*/   constants from the code.  from mattes (ike-stuttgart).
*d up47.294,295
      if (gfa.ne.zero) a2=sqrt(abs(gfa))
      if (gfa.lt.zero) a2=-a2
*/ reconr -- 15jan99
*/   fix an incorrect comment line, and replace one
*/   last real*4 constant in an if (trkov, ijs slovenia).
*d up47.408
      one=1
      if (thresh.gt.one.and.abs(thresh-eg).lt.test) sn=0
*d reconr.1494
c     for energies less than elim, add points by dividing
*d reconr.4191
*ident up56
*/ broadr -- 13jan99
*/   make trange incommensurable with normal grid energies
*/   to help in getting the same grid on different machines
*d up11.46
      data trange/.4999d0/
*d up11.48
      data trange/.4999/
*ident up57
*/ acer -- 13jan99
*/   fix an error in up38.  this update introduced a conversion
*/   to law=7 format for file 6 cm or lab tabulated or legendre 
*/   data (needed for jef iron).  unfortunately, the conversion
*/   is running for endf kalbach data also!
*d up38.12,13
      lang=l1h
c     set flag to convert file-6 law=1 data representation 
c     to law-7 angle-energy format except for kalbach data.
      if (lf.eq.1.and.lang.ne.2) new6=1
*ident up58
*/ heatr -- 14jan99
*/   add initialization for zero.
*/   found by both deleege (dlft) and muir (iaea).
*i heatr.3737
      zero=0
*ident up59
*/ groupr -- 14jan99
*/   fix an error introduced by up52.  the 64-bit versions of the
*/   constants bgam2 and tgam2 were clobbered by some of the
*/   other revised constants.  this does not effect "sw" versions.
*/   noted by deleege (delft).
*i groupr.1498
      data bgam2,tgam2/27.63e0,-0.53063e0/
*ident up60
*/ wimsr -- 15jan99
*/   add additional space for fission products and fix the
*/   consistency between the container array size and dimensions.
*/   requested by trkov (ijs, slovenia).
*d wimsr.101  
      common/wstore/a(400000)
*d wimsr.110  
      common/burnup/yield(100),ifisp(100)
*i wimsr.118
      nymax=100
*i wimsr.202
      if (ntis.gt.nymax) call error('wimsr',
     1   'too many time-dependent isotopes',' ')
*d wimsr.1832
      common/wstore/a(400000)
*d wimsr.1842 
      common/burnup/yield(100),ifisp(100)
*ident up61
*/ njoy -- 15jan99
*/   make a change recommended by trkov (ijs slovenia) that
*/   prevents an infinite loop for some pc compilers.
*d up46.59   
      if (y.gt.zero) xnext=shade*shade*a(jp)
*ident up62
*/ matxsr -- 15jan99
*/   remove some obsolete coding.  this shouldn't affect
*/   answers, but it stops a few annoying compiler messages.
*d matxsr.1681
*d matxsr.1684
      if (abs(a(nxsd)).lt.eps) go to 510
*d matxsr.1739,1757
*d matxsr.1851
*d matxsr.1859
      call shufl(nscrr,irefr,ia(ijgll),ia(ijgll+noutg),
*d matxsr.2007
      subroutine shufl(nscr,iref,jg1lo,jg1hi,nout3,irec3,irite)
*ident up63
*/ njoy -- 27jan99
*/   fix a typing error when adding double constants.
*/   reported by mattes, muir, panini.
*d up46.148
*i njoy.2665
      data zero/0.d0/
*ident up64
*/ broadr -- 27jan99
*/   fix a typing error when adding double constants.
*/   only affects 64-bit installations.  from mattes.
*d broadr.1200
      data sigmin/1.e-15/
*ident up65
*/ unresr -- 27jan99
*/   fix a typing error when adding double constants.
*/   only affects 64-bit installations.  from mattes.
*d up35.15
     7  1.25e6,1.5e6,2.0e6,3.2e6,4.0e6,5.0e6,6.3e6,8.0e6/
*ident up66
*/ heatr -- 27jan99
*/   fix a typing error when adding double constants.
*/   only affects 64-bit installations.  from mattes.
*d heatr.1598
      data step/1.5e0/
*ident up67
*/ groupr -- 27jan99
*/   fix a typing error when adding double constants.
*/   only affects 64-bit installations.  from mattes.
*d up52.99
      data up,dn/1.001e0,0.999e0/
*ident up68
*/ purr -- 27jan99
*/   fix a typing error when adding double constants.
*/   only affects 64-bit installations.  from mattes.
*d up36.15
     7  1.25e6,1.5e6,2.0e6,3.2e6,4.0e6,5.0e6,6.3e6,8.0e6/
*ident up69 
*/ groupr -- 28jan99
*/   add a missing save variable.  this variable is an
*/   index to a scratch area, and the problem can be
*/   intermittent, depending on the machine.  it only
*/   occurs for discrete levels (mt=51,52,...) when
*/   given in file 6.  noticed by frankle (lanl).
*d up25.5
      save ltt3,lttn,ifls
*ident up70
*/ dtfr -- 28jan99
*/   fix a problem with the graphics file for the
*/   neut-phot matrix from dtfr.  it shows up in viewr.
*d dtfr.1397
      write(nplot,'(1x