c                                                       October 2006
c     PRECOM3.for contains the more independent subprograms
c          ANGEL                          REPENT
c          BIND                           RESOUND       
c          CROSS                          SETISO 
c          INEL                           SIGPAR 
c          KNOCK                          SPEC   
c          PAULI                          SPECOUT
c
c*********************************************************************
c       
      subroutine angel(kp,length,newpg)
c     written in 1980; revised 1990, 1992
c
c     Continuum angular distributions
c     From empirical systematics based on exponentials of cos theta.
c
c     Third term in 'a' is needed for incident n,p,d but not alpha
c          Tentatively used for incident t, He-3
c     Transition energies tentatively set at
c          Et1 = 130 MeV (n,p,d,t,He3 incident) or 260 MeV (alpha)
c          Et3 = 35 MeV (n,p,d,t,He3 incident)  [not used for alpha]
c     (Only nucleon values of Et1 and n,p,d values of Et3 known.)
c
c     called from: PRECO
c     calls to: BIND
c
c     *The angel of the Lord encampeth round about them that fear him,
c     *and delivereth them.  (Psalm 34.7)
c
      dimension fmsd(51), jang(19), sigma(19), slope(51), slopef(51),
     1 slopes(51), total(51), xcos(19), xnorm(51), xnormf(51),
     2  xnorms(51)
      common /alph/ title
      common /angels/ elastf(51), elcoll(51), kang, tomsc(51),
     1 tomsd(51), xs(6)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      character*64 title
      character*1 newpg
      real*8 arg
      do 2 i=1,10
2     sigma(i) = 0
      do 4 i=1,19
      theta = i*10 - 10
      jang(i) = theta
      theta = theta * 3.1413 / 180.
4     xcos(i) = cos(theta)
      kin = jin + jpin
c     ---------------------------------------
c     calculate energy independent quantities
c     and print headings
c     ---------------------------------------
      jout = jpout(kp) + jnout(kp)
      xpin = jpin
      azcom = rzz + xpin
      esys = e - ben(1,7) + bind(acom(1),azcom,jpin,jnin)
c     ** esys is denoted e_a in the program writeup.
      bin = bind(acom(1),azcom,jpout(kp),jnout(kp))
      lpage = neps1(1,kp)-1
      lpage = 22+lpage+lpage+lpage
      if (length.ge.lpage) then
         write (iwri,18)
         write (iwri,18)
      else
         write (iwri,8) newpg,title
8        format (a1,' ',a64)
      end if
      write (iwri,12) jpout(kp),jnout(kp)
12    format (' angular distributions (mb/MeV)  z=',1i2,'  n=',1i2)
      write (iwri,14) esys,bin
14    format (17h sys. excit. en.=,1f6.2,18h;  sys. bind. en.=,1f5.2)
      er = 1.
      if (jin.gt.3.) er = 2.
c     ** er = Et1/130.
      esysr = esys / er
      e1 = esys
      if (esysr.gt.190.) then
           e1 = 130.*er
      else if (esysr.gt.80.) then
           x = (esysr-130.) / 12.
           x = 1. + exp(x)
           e1 = esys /x
           x = (136.-esysr) / 12.
           x = 1. + exp(x)
           e1 = e1 + 130.*er/x
      end if
c     ** Next, er should be reset to Et3/35 and esysr recalculated.
c     ** The value of er=1 from term 1 for jin.le.3 is appropriate for 
c     ** nucleons and is assumed here for d, t and h.
c     ** Term 3 (E3 and Et3 in it) is not used for jin=4.
c     ** Thus just use er = 1 and esysr = esys.
      er = 1
      esysr = esys
      e3 = esys
      if (esysr.gt.51.) then
           e3 = 35. * er
      else if (esysr.gt.21.) then
           x = (esysr-35.) / 3.2
           x = 1. + exp(x)
           e3 = esys / x
           x = (36.6-esysr) / 3.2
           x = 1. + exp(x)
           e3 = e3 + 35.*er/x
      end if
      xmb = 1
      if (jout.eq.4) xmb=2.
      if (jpout(kp).eq.0) xmb=0.5
c     -----------------------------------------
c     calculate and print angular distributions
c     -----------------------------------------
      write (iwri,18)
      write (iwri,16) (jang(j),j=1,10)
16    format (' eps:deg',1i3,9i8)
      write (iwri,18)
18    format (1h )
      do 28 ne=2,neps1(1,kp) 
      tot = tomsd(ne) + tomsc(ne)
      tote = tot
      if (kp.eq.kin) tote = tote + elastf(ne) + elcoll(ne)
      total(ne) = tote
      if (tote.le.0.) go to 28
      fmsd(ne) = 0.
      if (tot.gt.0.) fmsd(ne) = tomsd(ne)/tot
      epscm = eps(ne)+bin
      y = epscm*e1 / (esys*130.)
      a = 5.2*y + 4.*y*y*y
      y = epscm*e3 / (esys*35.)
      a3 = 1.9*y*y*y*y*xmb
c     ** Parameters for main components
      if (jin.lt.4) a = a + a3
c     **** Provisionally, Ma = 2 for incident neutrons; add a3 in again
      if (jin+jpin.eq.1) a = a + a3
      xnorm(ne) = a * tot / (12.5664*sinh(a))
      slope(ne) = a
      if (kp.eq.kin) then
c          ** Parameters for elastic and collective components
           acoll = 1.5 * a
           if (kp.lt.3) then
                acoll = a
           else if (kp.eq.4) then
                acoll = 1.9 * a
           else if (kp.eq.6) then
                acoll = 4.3 * a
           end if
           slopes(ne) = acoll
           aelfwd = acoll + 18.
           slopef(ne) = aelfwd
           xnorms(ne) = acoll * elcoll(ne) / (12.5664*sinh(acoll))
           xnormf(ne) = aelfwd * elastf(ne) / (12.5664*sinh(aelfwd))
      end if
c     ** Calculate double differential cross sections up to 90 deg
      do 24 i=1,10
      arg = a * xcos(i)
      sig = fmsd(ne)*dsinh(arg) + dcosh(arg)
      sig = sig * xnorm(ne)
      if (kp.eq.kin) then
c          ** Include elastic and collective components
c          ** (Prelim.  Slopes not well benchmarked for complex proj.)
           arg = acoll * xcos(i)
           sig = sig + (dsinh(arg)+dcosh(arg))*xnorms(ne)
           arg = aelfwd * xcos(i)
           sig = sig + (dsinh(arg)+dcosh(arg))*xnormf(ne)
      end if
24    sigma(i) = sig
      write (iwri,26) eps(ne),(sigma(i),i=1,10)
26    format (1f6.2,10(1pe8.2e1))
28    continue
      lpage = 9 + 2*neps1(1,kp)
      if (length.lt.lpage) then
           write (iwri,8) newpg,title
           write (iwri,12) jpout(kp),jnout(kp)
      else
           write (iwri,18)
      end if
      write (iwri,18)
      write (iwri,30)(jang(j),j=11,19)
30    format (' eps:deg',1i3,8i8,5x,'total')
      write (iwri,18)
c     ** Calculate double differential cross sections at back angles
      do 38 ne=2,neps1(1,kp)
      if (total(ne).le.0.) go to 38
      do 34 j=11,19
      arg = slope(ne) * xcos(j)
      sig = fmsd(ne)*dsinh(arg) + dcosh(arg)
      sig = sig * xnorm(ne)
      if (kp.eq.kin) then
           arg = slopes(ne) * xcos(j)
           sig = sig + (dsinh(arg)+dcosh(arg))*xnorms(ne)
           arg = slopef(ne) * xcos(j)
           sig = sig + (dsinh(arg)+dcosh(arg))*xnormf(ne)
      end if
34    sigma(j) = sig
      write (iwri,36) eps(ne), (sigma(i),i=11,19), total(ne)
c      write (iwri,36) eps(ne), (sigma(i),i=11,19), slope(ne)
c     ** Above is alternative printout giving the main slope parameter
36    format (1f6.2,10(1pe8.2e1))
38    continue
      return
      end
c
c*********************************************************************
c
      function bind(acom,azcom,jpout,jnout)
c
c     Calculate binding energies from semi-empirical masses
c     with pairing (and shell) effects not included.
c
c     called from: ANGEL
c
      common /angels/ elastf(51), elcoll(51), kang, tomsc(51),
     1 tomsd(51), xs(6)
      ancom = acom - azcom
      xout = jpout + jnout
      kp = jpout + jpout + jnout
      xpout = jpout
      ares = acom - xout
      azres = azcom - xpout
      anres = ares - azres
      acom3 = acom**0.33333
      ares3 = ares**0.33333
      bin = 15.68 * xout
      bin = bin - 18.56*(acom3*acom3-ares3*ares3)
      x = (ancom-azcom) * (ancom-azcom) / acom
      y = (anres-azres) * (anres-azres) / ares
      bin = bin - 28.07*(x-y)
      bin = bin + 33.228*(x/acom3-y/ares3)
      bin = bin - 0.717*(azcom*azcom/acom3-azres*azres/ares3)
      bin = bin + 1.211*(azcom*azcom/acom-azres*azres/ares)
      bind = bin - xs(kp)
      return
      end
c
c*********************************************************************
c
      subroutine cross(kc,kp)
c
c     written in 1982; revised 1990
c
c     Calculate optical model reaction cross sections
c     using the empirical parameterization
c     of Narasimha Murthy, Chaterjee, and Gupta
c     going over to the geometrical limit at high energy.
c
c           Proton cross sections scaled down with signor for a<100
c           (appropriate for becchetti-greenlees potential).
c           p2 reduced, and global red'n factor introduced below Bc.
c           Neutron cross sections scaled down with signor for a<40
c           Scaled up for A>210 (appropriate for Mani et al potential).
c
c     parameter values set in subroutine sigpar
c
c     called from: PRECO
c
c     *For the preaching of the cross is to them that perish 
c     *foolishness, but unto us which are saved it is the
c     *power of God. (1 Corinthians 1.18)
c
      common /direct/ betalam(10), coul(7), elam(10), elast(51), fracel,
     1 lam(10), lamax, scoll(51), sigbar(7), sigcni, siggeom(51), 
     2 snock(6,51), snutra(6,51), wcoll(10)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      common /par/ xl0(6), xl1(6), xm0(6), xm1(6) ,xn0(6), xn1(6),
     1 xn2(6), xp0(6), xp1(6), xp2(6)
      flow = 1.e-18
      spill = 1.e+18
      kin = jnin + jpin + jpin
      jout = jpout(kp) + jnout(kp)
      xout = jout
      ares = acom(kc) - xout
      athrd = ares**0.3333
      signor = 1.
c     ** signor reduces p and n result for light targs as per expt.
      if (kp.eq.1) then
           if (ares.lt.40.) signor=0.7+ares*0.0075
           if (ares.gt.210.) signor = 1. + (ares-210.)/250.
           xlamb = xl0(1)/athrd + xl1(1)
           xmu = xm0(1)*athrd + xm1(1)*athrd*athrd
           xnu = xn0(1)*athrd*ares + xn1(1)*athrd*athrd + xn2(1)
           ec = 0.5
           ecsq = 0.25
           p = xp0(1)
           xnulam = 1.
           etest = 32.
c          ** etest is the energy above which the rxn cross section is
c          ** compared with the geometrical limit and the max taken.
c          ** xnulam here is a dummy value for neutrons needed later.
           ra = 0.
      else 
           ra = 1.20
           if (kp.eq.2) then
                ra = 0.
                if (ares.le.60.) then
                     signor = 0.92
                else if (ares.lt.100.) then
                     signor = 0.8 + ares*0.002
                end if
           end if
           xpout = jpout(kp)
           rz = rzz + jpin - jpout(kp)
           if (kc.eq.3) rz = rz - 1.
           ec = 1.44 * xpout * rz / (1.5*athrd+ra)
           ecsq = ec * ec
           p = xp0(kp) + xp1(kp)/ec + xp2(kp)/ecsq
           xlamb = xl0(kp)*ares + xl1(kp)
           a = ares**xm1(kp)
           xmu = xm0(kp) * a
           xnu = a* (xn0(kp)+xn1(kp)*ec+xn2(kp)*ecsq)
           if (jout.eq.2) ra=0.8
           if (jout.eq.3) ra=0.8
c          ** New values of ra are for calculating the geometrical limit
c          ** to the cross section.
           if (kp.eq.2) then
                c = min(3.15,ec*0.5)
                w = 0.7 * c / 3.15
c               **C and w are for global corr'n factor for elab<ec.
c               ** They are scaled down for light targs where ec is low.
           end if
           xnulam = xnu / xlamb
           if (xnulam.gt.spill) xnulam=0.
           if (xnulam.lt.flow) go to 20
           if (kp.eq.2) then
                etest = sqrt(xnulam) + 7.
           else
                etest = 1.2 * sqrt(xnulam)
           end if
c          ** For xnulam.gt.0, sig reaches a maximum at sqrt(xnulam).
      end if                                           
20    a = -2.*p*ec + xlamb - xnu/ecsq
      b = p*ecsq + xmu + 2.*xnu/ec
      ecut = 0.
      cut = a*a - 4.*p*b
      if (cut.gt.0.) ecut = sqrt(cut)
      ecut = (ecut-a) / (p+p)
c     ** If cut>0, sigma=0 at elab=ecut
c     ** If cut<0, OM sigma (no global corr'n) is at a min at elab=ecut
      ecut2 = ecut
      if (cut.lt.0) ecut2 = ecut - 2.
c     ** ecut2 is the energy below which sigma is set to zero
      neps = neps1(1,kp)
      if (kc.gt.1) neps = max0(neps1(1,1),neps1(1,2))
      do 36 ne=2,neps
      elab = eps(ne) * acom(kc) / ares
      sig = 0.
      if (elab.le.ec) then
           if (elab.gt.ecut2) then
                sig = (p*elab*elab+a*elab+b) * signor
                if (kp.eq.2) then
                     signor2 = (ec-elab-c) / w
                     signor2 = 1. + exp(signor2)
                     sig = sig / signor2
c                    signor2 is empirical global corr'n at low elab
                end if
           end if
      else
           sig = (xlamb*elab+xmu+xnu/elab) * signor
           geom = 0.
           if (xnulam.lt.flow) go to 34
           if (elab.lt.etest) go to 34
           geom = sqrt(xout*eps(ne))
           if (kp.le.2) then
                geom = 1.23*athrd + ra + 4.573/geom
           else
                geom = 1.12*athrd + 0.9 + ra + 4.573/geom
           end if
           geom = 31.416 * geom * geom
           sig = amax1(geom,sig)
      end if
34    if (kc.eq.1) then
           if(kp.eq.kin) then
c               ** Calculate the geometric cross sections for collective
c               ** state excitation calculations in RESOUND
                geom = sqrt(xout*eps(ne))
                if (kp.le.2) then
                     geom = 1.23*athrd + ra + 4.573/geom
                else
                     geom = 1.12*athrd + 0.9 + ra + 4.573/geom
                end if
                siggeom(ne) = 31.416 * geom * geom
           end if
      end if
36    sigin(kp,ne) = sig
      return
      end
c
c*********************************************************************
c
      subroutine inel(kp)
c
c     Phenomenological energy spectra for inelastic scattering
c     of complex particles by n, p and alpha pair excitation
c     and of nucleons with alpha pair excitation.
c     Calculations do not include T conservation, shells, pairing, or
c     finite well depth corrections.
c
c     called from: PRECO
c     call to: RESOUND
c
      common /direct/ betalam(10), coul(7), elam(10), elast(51), fracel,
     1 lam(10), lamax, scoll(51), sigbar(7), sigcni, siggeom(51), 
     2 snock(6,51), snutra(6,51), wcoll(10)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3),
     1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff
      flow9 = 1.001e-09
      spill = 1.e+18
      xpin = jpin
      xnin = jnin
      xin = jin
      acomp = acom(1)
      atarg = acomp-xin
      rnn = atarg-rzz
c
c     * * * * * * * * * * * * * * * * * * * * * * * *
c     * Calculate normalization denominators and    *
c     * the non-energy dependent part of main terms *
c     * * * * * * * * * * * * * * * * * * * * * * * *
c
      xnor = jin * nsd(kp)
      xnorm = xnor * sigcni * ct(1,7,1) / 14.
c     ----------------------------------------------
c     1. projectile=ejectile
c        s.p. state densities & emission denominator
c     ----------------------------------------------
      xnor = xnor * sigbar(kp)/6.
      coup = coul(kp)
      emax = e - ben(1,kp)
      xnor = xnor * (emax+2.*coup) * (emax-coup) * (emax-coup)
      if (jin.gt.1) go to 2
c     ** gbn (gba) is single particle state density for a kp-type
c     ** particle in the nucleon (alpha) residual nucleus
      gba = (rzz-1.)/13.
      if (jnin.eq.1) gba=(rnn-1.)/13.
      go to 10
2     if (jin.eq.2) then
           gbn = (acomp-1.) / 52.
           gba = (acomp-4.) / 52.
      else if (jin.eq.3) then
           gbn = (acomp-1.) / 156.
           gba = (acomp-4.) / 156.
      else if (jin.eq.4) then
           gbn = (acomp-1.) / 208.
           gba = (acomp-4.) / 208.
      end if
c     -------------------------------
c     2. neutron excitation term (xn)
c     -------------------------------
      gii = (rnn+xnin-1.)/13.
      gir = rnn/13.
c     ** gii (gir) is the single particle state density for the 
c     ** excited particle and hole in their own (the inelastic)
c     ** residual nucleus.
      p = rnn/atarg
c     **  p = N/A is the probability of exciting a neutron pair
      xnori = 2. * sigbar(1) * gii*gbn/6.
      emax = e - ben(1,1)
      coup = coul(1)
      xnori = xnori * (emax+2.*coup) * (emax-coup)*(emax-coup)
      xnorb = xnor * gir*gir
      xn = gir*gir * p / (xnori+xnorb)
      thren = 1. / gir
c     ------------------------------
c     3. proton excitation term (xp)
c     ------------------------------
      gii = (rzz+xpin-1.) / 13.
      gir = rzz / 13.
      p = rzz / atarg
c     ** p = Z/A is the probability of exciting a proton pair
      xnori = 2. * sigbar(2) * gii*gbn/6.
      emax = e - ben(1,2)
      coup = coul(2)
      xnori = xnori * (emax+2.*coup) *(emax-coup)*(emax-coup)
      xnorb = xnor * gir*gir
      xp = gir*gir * p / (xnori+xnorb)
      threp = 1. / gir
c     -----------------------------
c     4. alpha excitation term (xa)
c     -----------------------------
10    gii = (acomp-4.) / 208.
      gir = atarg / 208.
      p = 0.08 * rzz / (2.*atarg)
c     ** p is the probability of exciting an alpha pair
      if (rnn.le.116.) go to 14
      if (rnn.ge.129.) go to 14
      if (rnn.le.126.) then
           y = 0.02 + 0.06*(126.-rnn)/10.
           p = p * y / 0.08
      else
           y = 0.02 + 0.06*(rnn-126.)/3.
           p = p * y / 0.08
      end if
14    xnori = 4. * sigbar(6) * gii*gba/6.
      emax = e - ben(1,6)
      coup = coul(6)
      xnori = xnori * (emax+2.*coup) * (emax-coup)*(emax-coup)
      xnorb = xnor * gir*gir
      xa = gir*gir * p / (xnori+xnorb)
      if (jin.eq.4) xa = xa * 2
c     ** With alphas incident, there is only one kind of emission from
c     ** alpha excitation, so xnori=xnorb and no sum in denominator,
c     ** because competition is only among different energies of alphas.
      threa = 1. / gir
c
c     * * * * * * * * * * * * * * *
c     * calculate energy spectrum *
c     * * * * * * * * * * * * * * *
c
      umax = e - ben(1,kp)
      do 16 ne=2,neps1(1,kp)
      u = umax - eps(ne)
c     ** note: because all residual states have a particle-hole pair
c     ** of a single kind (n, p, or alpha) the pauli correction term
c     ** is zero.
      if (u.le.0.) go to 18
      x = 0.
      if (u.gt.threa) x=xa
      if (jin.gt.1) then
           if (u.gt.thren) x = x + xn
           if (u.gt.threp) x = x + xp
      end if
      sig = xnorm * eps(ne) * sigin(kp,ne) * x * u
      if (sig.lt.flow9) sig=0.
      if (sig.gt.spill) sig=0.
16    snock(kp,ne) = sig
c     
c     * * * * * * * * * * * * * * * * * * * * * * * * * *
c     * Contribution from strong collective resonances  *
c     * * * * * * * * * * * * * * * * * * * * * * * * * *
c     
18    call resound(kp,ct(1,7,1))
      return
      end
c
c*********************************************************************
c
      subroutine knock(kp)
c
c     written in 1980, revised 1990
c
c     Phenomenological energy spectra for knockout of n, p, or alpha
c     rxn must involve at least one complex particle.
c     Calculations do not include t conservation, shells, pairing or
c     finite well depth corrections.
c
c     called from: PRECO
c
c     * Ask, and it shall be given you.
c     * Seek, and ye shall find.
c     * Knock, and it shall be opened unto you.  (Matthew 7.7)
c
      common /direct/ betalam(10), coul(7), elam(10), elast(51), fracel,
     1 lam(10), lamax, scoll(51), sigbar(7), sigcni, siggeom(51), 
     2 snock(6,51), snutra(6,51), wcoll(10)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3),
     1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff
      flow9 = 1.001e-09
      spill = 1.e+18
      xin = jin
      xpin = jpin
      xnin = jnin
      jout = jpout(kp) + jnout(kp)
      xout = jout
      xpout = jpout(kp)
      xnout = jnout(kp)
      acomp = acom(1)
      atarg = acomp - xin
      rnn = atarg - rzz
      ares = acomp - xout
c     --------------------------------------------------------
c     calculate s.p. state densities and (2s+1)*A type factors
c     --------------------------------------------------------
      if (jin.gt.1) go to 2
      rel=2.
c     ** ga is the single particle state density for the projectile
c     ** in the real residual nucleus.
      ga = (rzz+1.-xpout)/13.
      if (jnin.eq.1) ga=(rnn+1.-xnout)/13.
      go to 8
2     if (jin.eq.2) then
           rel = 6.
           ga = ares / 52.
      else if (jin.eq.3) then
           rel = 6.
           ga = ares / 156.
      else if (jin.eq.4) then
           rel = 4.
           ga = ares / 208.
      end if
8     if (jout.gt.1) go to 10
c     ** gbt (gbr) is the s.p. state density for the emitted particle
c     ** type in the inelastic (real) residual nucleus.
      gbt = rzz / 13.
      gbr = (rzz+xpin-1.) / 13.
      if (jnout(kp).eq.0) go to 12
      gbt = rnn / 13.
      gbr = (rnn+xnin-1.) / 13.
      go to 12
10    gbt = atarg / 208.
      gbr = ares / 208.
c     -----------------------
c     calculate normalization
c     -----------------------
12    thresh = 1./gbr + 1./ga
      xnorm = sigcni * ct(1,7,1) / 14.
      xnorb = jout * nsd(kp)
      xnorm = xnorm * xnorb
      xnorb = xnorb * sigbar(kp) * ga*gbr/6.
      emax = e - ben(1,kp)
      coup = coul(kp)
      xnorb = xnorb * (emax+2.*coup) * (emax-coup)*(emax-coup)
      xnora = rel*sigbar(7)*gbt*gbt/6.
      emax = e - ben(1,7)
      coup = coul(7)
      xnora = xnora * (emax+2.*coup) * (emax-coup)*(emax-coup)
      xnorm = xnorm * ga * gbr / (xnora+xnorb)
c     ** Statements below multiply the normalization by the probability 
c     ** of the projectile striking the knocked out particle.
      if (jpout(kp).eq.0) then
           xnorm = xnorm * (1.-rzz/atarg)
      else if (jpout(kp).eq.1) then
           xnorm = xnorm * rzz / atarg
      else
           xnorm = xnorm * 0.08 * rzz / (2.*atarg)
           if (rnn.le.116.) go to 20
           if (rnn.ge.129.) go to 20
           if (rnn.le.126.) then
                x = 0.02 + 0.06*(126.-rnn)/10.
                xnorm = xnorm * x / 0.08
           else
                x = 0.02 + 0.06*(rnn-126.)/3.
                xnorm = xnorm * x / 0.08
           end if
      end if
c     ------------------------
c     calculate energy spectra
c     ------------------------
20    apauli = thresh * 0.5
c     ** The pauli correction term above is based on cluster degrees
c     ** of freedom: a 1p0h configuration of the projectile type plus
c     ** 0p1h of the ejectile type.  This is consistent with the state
c     ** density calculations.
      umax = e - ben(1,kp) - apauli
      do 22 ne=2,neps1(1,kp)
      u = umax - eps(ne)
      if (u.lt.apauli) go to 24
      sig = xnorm * u * eps(ne) * sigin(kp,ne)
      if (sig.lt.flow9) sig=0.
      if (sig.gt.spill) sig=0.
22    snock(kp,ne) = sig
24    return
      end
c
c*********************************************************************
c
      subroutine pauli(kc,kp)
c
c     Calculate pauli correction function for state densities in the
c     ESM or S2-ESM and including simple pairing.
c     (Collective pairing and the third term in A are energy dependent
c     and must be included later in the calculations.)
c
c     Calculate equilibrium shell and pairing shifts.
c
c     Calculate preliminary variables used in subroutine SINGLE
c     to get shell corrected single particle state densities.
c
c     called from: PRECO
c
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin,
     1 xncrit(7,2)
      common /shell/ db(2), ds(7,2), epauli(7,14,2), dodd(7,2), fspan,
     1 ghole(2), gpart(2), idif(7,2), ishell, mdnear(7,2), mdfar(7,2),
     2 ndwn, nup2, range(2), rose(7), thresh(7,14,2)
      rose(kp) = 0.
      do 20 ia=1,2
      dbig = db(ia)
      dsmall = ds(kp,ia)
      dsx = dsmall*0.5
      phd = nphd(kc,kp,ia)
      if (phd.lt.0.) phd=-phd
      if (dbig.le.0.) go to 10
      mov = (dbig+dsmall) / (dsmall+dsmall)
      xmove = mov + mov
      dodr = (xmove*dsmall-dbig) * 0.5
c     ** dodr is the odd spacing minus dsmall/2
      dodd(kp,ia) = dodr + dsx
      if (mdnear(kp,ia).le.0) then
           mdnear(kp,ia) = mov
           mdfar(kp,ia) = mov
      end if
      mnear = mdnear(kp,ia)
      mfar = mdfar(kp,ia)
      idifa = iabs(idif(kp,ia))
      difa = idifa
      fwash = difa / range(ia)
      na = iabs(mnear-idifa)
      nc = mfar + idifa
c     ------------------------------
c     1.  Pauli Correction Functions
c     ------------------------------
10    do 16 nq1=2,nup2
      nq = nq1-1
      q = nq
c     ** Normally these calculations would be bypassed for q<phd.
c     ** However, when nucleon transfer at the Fermi level is allowed
c     ** in NUTRA, then q<phd is needed, but qh has to be set to zero.
      paul = q*q * dsmall
      if (dbig.gt.0.) then
           if (fwash.lt.1.) then
                xa = amin0(nq,na)
                pa = q-xa
                pb = nq - min0(nq,idifa)
                pc = nq - min0(nq,nc)
                pauls = (pa+pc)*dodr + (pa*pa+pc*pc)*dsx + pb*dbig
                if (idifa.gt.mnear) pauls = pauls + xa*(q+q-xa)*dsmall
                paul = (1.-fwash)*pauls + fwash*paul
           end if
      end if
      qh = amax1(q-phd,0.)
      thresh(kp,nq1,ia) = paul + cpre(kc,kp,ia)
      paul = paul - (q*(q+1.)+qh*(qh+1.))*dsmall*0.25
      epauli(kp,nq1,ia) = paul + cpre(kc,kp,ia)
16    continue
c     ---------------------------------
c     2.  Equilibrium Shell Corrections
c     ---------------------------------
      if (dbig.gt.0.) then
           if (fwash.lt.1) then
                fwash = sqrt(fwash)
                xmove = mov
c               ** Next line gives closed shell equilibrium correction.
                paul = xmove * (dbig*0.5 - xmove*dsx)
                paul = (1.-fwash) * paul
                rose(kp) = rose(kp) + paul
           end if
      end if
      rose(kp) = rose(kp) + ceq(kc,kp,ia)
20    continue
      return
      end
c
c     *****************************************************************
c
      subroutine repent(kc,neps)
c
c     Calculates recycling points in Main and sets values for needed 
c     variables when multiple preequil. particle emission is allowed.
c
c     called from: PRECO
c     Calls to: CROSS and SETISO
c
c     * Then Peter said unto them, "Repent and be baptized every one of 
c     * you in the name of Jesus Christ for the remission of sins, and
c     * ye shall receive the gift of the Holy Spirit."  (Acts 2:38)
c
      common /emit/ clos(13,13), fiss(3), fweis, fweist, fteq, 
     1 mxdime, mxdimp, mxdim2, ndwn1, nup1, preeq1(6,51), preeq2(2,51),
     2 ptot, rho(13,13), ttot(6,51), weiss(6,51), weist(6,51)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      common /repeat/ bfis(3), ecut, ifiss, irecyc, kpmax, nemin,
     1 neprev, nepsin, einit, rgamma, rzcn(3), rncn(3)
      common /transit/ gamhp(13,13), gamhn(13,13), gampn(13,13),
     1 gamnp(13,13), gh(2), gp(2)
      common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3),
     1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff
      common /wmain/ w(3,13,14,51), strong(9,10,51), seq(6,51),
     1 seq2(2,51), steq(6,51), steq2(2,51), xiso(9,10,51), xisoeq(2,51),
     2 spre(2,51)
c
      if (kc.eq.1) then
c          ** Choose first energy for emission from n-resid.,
c          ** skipping over those with sigin=0
           write(*,410)
410        format(' Emission from neutron residual')
           kc = 2
           kpmax = 2
           neprev = 1
412        neprev = neprev+1
           if (neprev.gt.neps1(1,1)) go to 414
c          ** Above test bypasses primary n emission if not
c          ** energetically allowed.  (Highly unlikely.)
           if (sigin(1,neprev).le.0.) go to 412
           nemin = neprev
           e = einit - eps(neprev) - ben(1,1)
           echek = amin1(ben(2,1),ben(2,2)) + eps(2)
           if (e.lt.echek) go to 414
           if (nsigin.eq.0) then
                call cross(2,1)
                call cross(2,2)
           end if
           do 411 ne=2,neps1(1,2)
           nee = ne
           if (sigin(2,ne)/sigin(1,ne).gt.rgamma) go to 413
411        continue  
413        ecut = amin1(ben(2,1)+eps(2),ben(2,2)+eps(nee))
c          ** ECUT cuts off secondary emission calculations if there
c          ** is not enough energy for particle emission (or fission).
c          ** Using eps(nee) instead of eps(2) for the proton 
c          ** requirement takes account of low emission energy where 
c          ** sigin is zero or so small that gamma emission will
c          ** successfully compete.
           if (ifiss.gt.0) then
                ecut = amin1(ecut,bfis(2))
           end if
           if (iso.gt.0) then
                nt(2) = ntmin(1)
                call setiso(2,rzcn(2),rncn(2))
           end if
           if (e.gt.ecut) then
                if (e.gt.esym(2,7,1)) go to 436
                if (fteq.lt.0.95) go to 436
           end if
c          ** Statement 436 begins initializing arrays for new CN.
c          ** If E<Esym for the lowest primary emission energy, then no
c          ** emission is possible from the n-residual (except 2ndary
c          ** t-mixed equilibrium emission).  Try p-residual.
414        kc = 3
           write (*,420)
           neprev = 1
416        neprev = neprev+1
           if (neprev.gt.neps1(1,2)) go to 432
           if (sigin(2,neprev).le.0.) go to 416
           nemin = neprev
           e = einit - eps(neprev) - ben(1,2)
           echek = amin1(ben(3,1),ben(3,2)) + eps(2)
           if (e.lt.echek) go to 432
           if (nsigin.eq.0) then
                call cross(3,1)
                call cross(3,2)
           end if
           do 418 ne=2,neps1(1,2)
           nee = ne
           if (sigin(2,ne)/sigin(1,ne).gt.rgamma) go to 419
418        continue  
           write (*,420)
419        ecut = amin1(ben(3,1)+eps(2),ben(3,2)+eps(nee))
           if (ifiss.gt.0) then
                ecut = amin1(ecut,bfis(3))
           end if
           if (iso.gt.0) then
                nt(3) = ntmin(2)
                call setiso(3,rzcn(3),rncn(3))
           end if
           if (e.gt.ecut) then
                if (e.gt.esym(3,7,1)) go to 436
                if (fteq.lt.0.95) go to 436
           end if
      else if (kc.eq.2) then
           if (neprev.lt.neps1(1,1)) then
c               ** Set new energy in neutron residual
                neprev = neprev+1
                e = einit - eps(neprev) - ben(1,1)
                if (e.gt.ecut) then
                     if (e.gt.esym(2,7,1)) go to 436
                     if (nt(2).eq.ntmin(1)) then
                           if (fteq.lt.0.95) go to 436
                     end if
                end if
           else if (iso.gt.0) then
                if (nt(2).lt.ntmax(1)) then
c                    ** Another T is possible in n-residual.
                     neprev = nemin
                     e = einit - eps(neprev) - ben(1,1)
                     nt(2) = nt(2) + 1
                     call setiso(2,rzcn(2),rncn(2))
                     if (e.gt.ecut) then
                          if (e.gt.esym(2,7,1)) go to 436
                     end if
                end if
           end if
c          ** Choose first energy for p-residual
           write(*,420)
420        format(' Emission from proton residual')
           kc = 3
           neprev = 1
424        neprev = neprev+1
           if (neprev.gt.neps1(1,2)) go to 432
c          ** Above test bypasses primary p emission if not energetically
c          ** allowed.  
           if (sigin(2,neprev).le.0.) go to 424
           nemin = neprev
           e = einit - eps(neprev) - ben(1,2)
           echek = amin1(ben(3,1),ben(3,2)) + eps(2)
           if (e.lt.echek) go to 432
           if (nsigin.eq.0) then
                call cross(3,1)
                call cross(3,2)
           end if
           do 426 ne=2,neps1(1,2)
           nee = ne
           if (sigin(2,ne)/sigin(1,ne).gt.rgamma) go to 428
426        continue  
428        ecut = amin1(ben(3,1)+eps(2),ben(3,2)+eps(nee))
           if (ifiss.gt.0) then
                ecut = amin1(ecut,bfis(3))
           end if
           if (iso.gt.0) then
                nt(3) = ntmin(2)
                call setiso(3,rzcn(3),rncn(3))
           end if
           if (e.gt.ecut) then
                if (e.gt.esym(3,7,1)) go to 436
                if (fteq.lt.0.95) go to 436
           end if
      else
           if (neprev.lt.neps1(1,2)) then
                neprev = neprev+1
                e = einit - eps(neprev) - ben(1,2)
                if (e.gt.ecut) then
                     if (e.gt.esym(3,7,1)) go to 436
                     if (nt(3).eq.ntmin(2)) then
                          if (fteq.lt.0.95) go to 436
                     end if
                end if
           else if (iso.gt.0) then
                if (nt(3).lt.ntmax(2)) then
                     neprev = nemin
                     e = einit - eps(neprev) - ben(1,2)
                     nt(3) = nt(3)+1
                     call setiso(3,rzcn(3),rncn(3))
                     if (e.gt.ecut) then
                          if (e.gt.esym(3,7,1)) go to 436
                     end if
                end if
           end if
      end if
c     ** Secondary emission calculations are done!  Reset to original CN
432   kc = 1
      kpmax = 6
      e = einit
      neps = nepsin
      irecyc = 2
      if (nsigin.eq.0) then
           call cross(1,1)
           call cross(1,2)
      end if
      go to 450
c     ** Do initialization of arrays for secondary emission
436   do 438 ne=1,mxdime
      do 438 kp=1,3
      weiss(kp,ne) = 0.
      weist(kp,ne) = 0.
      do 438 npp=1,mxdimp+1
      do 438 np=1,mxdimp
438   w(kp,np,npp,ne) = 0.
      do 440 npp=1,mxdimp
      do 440 np=1,mxdimp
      gamhp(np,npp) = 0.
      gamhn(np,npp) = 0.
      gampn(np,npp) = 0.
      gamnp(np,npp) = 0.
440   clos(np,npp) = 0.
      do 444 kp=1,2
      emu = e - ben(kc,kp)
      do 442 ne=2,nepsin+1
      ner = ne
      if (eps(ne).gt.emu) go to 444
442   continue
      ner = ner+1
444   neps1(kc,kp) = ner-1
      neps = max0(neps1(kc,1),neps1(kc,2)) - 1
c     ** Next, for first time thru on an emitting nucleus, irecyc=0
c     ** Otherwise, irecyc=1
      if (neprev.eq.nemin) then
           irecyc = 0
           if (iso.gt.0) then
                if (nt(kc).gt.ntmin(kc-1)) then
c                    ** 2nd T value, so STRONG needs to be renormalized
c                    ** and xisoeq and spre need to be reset.
                     do 446 ne=1,mxdime
                     x = xisoeq(kc-1,ne)
                     spre(kc-1,ne) = spre(kc-1,ne) * (1./x - 1.)
                     xisoeq(kc-1,ne) = 1. - x
                     do 446 np1=1,mxdim2
                     do 446 nppi1=1,np1
                     x = xiso(np1,nppi1,ne)
                     if(x.gt.0.) x= strong(np1,nppi1,ne) * (1./x - 1.)
                     if (x.lt.0.) x=0.
446                  strong(np1,nppi1,ne) = x
                     irecyc = 1
                end if
           end if
           if (kc.eq.3) then
                if (irecyc.eq.0) then
c                    ** First pass thru p-residual; flip packed arrays
                     do 448 np1=1,mxdim2
                     mp1 = mxdim2 - np1 + 1
                     do 448 nppi1=1,np1
                     mppi1 = mxdim2 - nppi1 + 2
                     do 448 ne=2,mxdime
                     strong(np1,nppi1,ne) = strong(mp1,mppi1,ne)
448                  xiso(np1,nppi1,ne) = xiso(mp1,mppi1,ne)
                end if
           end if
      else
           irecyc = 1
      end if
450   return
      end
c
c*********************************************************************
c
      subroutine resound(kp,ctin)
c
c     Written 1998
c
c     Calculate contributions from strong collective resonance states
c          2+, 3- and possibly 4+ spectroscopic strength
c          and giant resonance states
c
c     called from INEL
c
      common /direct/ betalam(10), coul(7), elam(10), elast(51), fracel,
     1 lam(10), lamax, scoll(51), sigbar(7), sigcni, siggeom(51), 
     2 snock(6,51), snutra(6,51), wcoll(10)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      xin = jin
      einc = e - ben(1,kp)
      chcm = 1. - xin/acom(1)
      eincm = einc * chcm
      do 20 i = 1,lamax
      ilor = 0
c     ** ilor gives lineshape: ilor=0 for Guassian, 1 for Lorentzian
      if (lam(i).lt.0) then
           ilor = 1
           lam(i) = -lam(i)
      end if
      xl = lam(i) + lam(i) + 1
      wcol = wcoll(i)
      ielast = 0
c     ** ielast is zero for normal collectives; 1 for elastic peak
c     ** Now do special calculations for elastics, if requested.
      if (i.gt.1) go to 2
      if (lam(1).eq.0) then
           if (elam(1).lt.0.0001) then
c               ** Elastic component calc'd with collective model
                ielast = 1
                xpin = jpin
                bet = 1. + xpin/(xin*xin)
                at = (acom(1)-xin) ** 0.33333333
                betalam(i) = 1.5 * bet / at
                frel1 = elam(2) + elam(3)
                frel1 = frel1 * 0.5 / (xin*eincm)
c               ** fracel is fraction of elastic following the flatter
c               ** (collective model) angular distribution 
c               ** systematics.  The rest is more forward peaked.
c               ** Frel1 is the flat fraction for the elastic strength 
c               ** calculated using the collective model.
c               ** Now calculate Coulomb elastic component.
c               ** (frel2 is the flatter fraction.)
                strel = 0.
                frel2 = 0.
                if (jpin.gt.0) then
                     strel = xpin*xpin * rzz*rzz / (eincm*eincm)
                     rnn = acom(1) - xin - rzz
                     if (kp.eq.2) then
                          strel = 400. * strel
                          frel2 = 0.020
                     else if (kp.eq.3) then
                          strel = 1000. * strel
                          frel2 = 0.004 + (rnn-rzz)/(9*acom(1))
                     else if (kp.eq.4) then
                          strel = 1000. * strel
                          frel2 = 0.040 - (rnn-rzz)/(6*acom(1))
                     else if (kp.eq.5) then
                          strel = 1400. * strel
                          frel2 = 0.006 + (rnn-rzz)/(7*acom(1))
                     else
                          strel = 600. * strel
                          frel2 = 0.032
                     end if
                end if
           end if
      end if
2     strong = 304 / (xin * eincm**1.5)
      strong = strong * betalam(i) * betalam(i) / xl
      strong = strong * acom(1)**1.3333
      strong = strong * sigcni / siggeom(1)
      strong = strong * ctin
      if (ielast.gt.0) then
           str = strong * sqrt(eincm) * sigcni / siggeom(1)
           fracel = frel1*str + frel2*strel
           fracel = fracel / (str+strel)
      end if
      do 10 ne = 2,neps1(1,kp)
      u = e - ben(1,kp) - eps(ne)
      eout = eps(ne) * chcm
      if (ilor.eq.0) then
c          ** Gaussian line shape (recommended)
           edif = abs(u-elam(i))
           if (edif.gt.5.*wcol) go to 10
           prob = 0.5 * edif * edif / (wcol*wcol)
           prob = 0.399 / (exp(prob) * wcol)
      else
c          ** Lorentz line shape for indicated GR states 
           e0sq = elam(i)*elam(i) + 0.25*wcol*wcol
           term2 = u*u*wcol
           denom = (u*u - e0sq)**2 + term2*wcol
           prob = 0.637 * term2 / denom
      end if
      prob1 = prob *sqrt(eout) *sigin(kp,ne)/siggeom(ne)
      if (ielast.eq.0) then
           scoll(ne) = scoll(ne) + strong*prob1
      else
           elast(ne) = strong*prob1 + strel*prob
      end if
10    continue
20    continue
      return
      end
c
c*********************************************************************
c
      subroutine setiso(kc,rzcn,rncn)
c
c     Extracted from PRECO in 1995
c
c     Sets values for Tz, 2T(kp), Ct^2 and Esym.
c
c     Called from PRECO, REPENT
c
      common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3),
     1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff
      if (kc.eq.1) then
           ntmin(1) = 0
           ntmin(2) = 0
      end if
      tz(7) = (rncn-rzcn) * 0.5
      tz(1) = tz(7) - 0.5
      tz(2) = tz(7) + 0.5
      tz(3) = tz(7)
      t = nt(kc)
      twot(7) = 2.* (abs(tz(7))+t)
      twot(1) = 2.* abs(tz(1))
      if (twot(1).lt.twot(7)-1.1) then
           twot(1)=twot(1)+2.
           if (kc.eq.1) ntmin(1) = 1
      end if
      twot(2) = 2.* abs(tz(2))
      if (twot(2).lt.twot(7)-1.1) then
           twot(2)=twot(2)+2.
           if (kc.eq.1) ntmin(2) = 1
      end if
      twot(3) = twot(7)
c     ** For n, p, t, and He-3, two possible isospins are considered in
c     ** the residual nucleus. 
c     ** Twot(kp) is two times the lower allowed value. 
c     ** For n and t (assuming N>Z) the higher value is always allowed.
c     ** For p and he-3 it is only allowed if t(kp) < t(7).
c     ** The quantity ct(kp,it) is the square of the isospin coupling
c     ** Clebsch-Gordan coefficient for the kp channel.  
c     ** it=1 is for the lower residual isospin value, while it=2 is 
c     ** for the upper one.
c     ** If the upper isospin value is not allowed, ct(kp,2)=0.
c     ** Programming below for Ct(1,it) and Ct(2,it) is taken from eq.
c     ** (2.82) on pg 927 of deShalit and Feshbach.
      x = twot(1)+1.
      if (twot(1).lt.twot(7)) then
           ct(kc,1,1) = 0.5 + tz(7)/x
           ct(kc,1,2) = 0.5 - tz(7)/(x+2.)
           if (kc.eq.1) ntmax(1) = ntmin(1)+1
      else
           ct(kc,1,1) = 0.5 - tz(7)/x
           ct(kc,1,2) = 0.
           if (kc.eq.1) ntmax(1) = ntmin(1)
      end if
      x = twot(2)+1.
      if (twot(2).lt.twot(7)) then
           ct(kc,2,1) = 0.5 - tz(7)/x
           ct(kc,2,2) = 0.5 + tz(7)/(x+2.)
           if (kc.eq.1) ntmax(2) = ntmin(2)+1
      else
           ct(kc,2,1) = 0.5 + tz(7)/x
           ct(kc,2,2) = 0.
           if (kc.eq.1) ntmax(2) = ntmin(2)
      end if
      ct(kc,3,1) = 1.
      ct(kc,3,2) = 0.
      do 82 kp=1,3
      tz(kp+3) = tz(kp)
      twot(kp+3) = twot(kp)
      ct(kc,kp+3,1) = ct(kc,kp,1)
82    ct(kc,kp+3,2) = ct(kc,kp,2)
c     ** A separate calculation is done for each composite isospin.
c     ** Thus ct(7,2) is always zero.
c     **
c     ** The next loop generates itc(kp) which determines whether the
c     ** N>Z or the Z>N form of the isospin correction is used. 
c     ** In the latter case ppi,hpi,gpi are interchanged with
c     ** pnu,hnu,gnu in TCOR.
c     ** It also generates tze(kp) which is the effective neutron (or       
c     ** proton, for Z>N) excess in the nucleus specified by kp.  
c     ** Most of the time tze(kp) = abs(tz(kin)).
      do 84 kp=1,7
      tzkp = tz(kp)
      if (itest.ge.0) then
         if (tzkp.ge.-0.001) then
            itc(kp) = 1
            tze(kp) = tzeff
         else
            itc(kp) = -1
            tze(kp) = 0
         end if
      else
         if (tzkp.le.0.001) then
            itc(kp) = -1
            tze(kp) = tzeff
         else
            itc(kp) = 1
            tze(kp) = 0
         end if
      end if
84    continue
      acom = rzcn + rncn
      x = 112./acom
c     ** In preco-2000, the multiplier above was (incorrectly) 110.
c     ** To use only the volume symmetry energy, comment out the 
c     ** next line and make the conforming change further down for
c     ** residual nuclei.
      x = x - 133./acom**1.3333
      t = twot(7)*0.5
      esym(kc,7,1) = x * (t-tz(7)) * (t+tz(7))
      esym(kc,7,2) = x * (t+1.-tz(7)) * (t+1.+tz(7))
      do 86 kp=1,6
      aout = 3.
      if (kp.lt.3) aout=1.
      if (kp.eq.3) aout=2.
      if (kp.eq.6) aout=4.
      x = 110. / (acom-aout)
c     ** To use only the volume symmetry energy, comment out the next 
c     ** line & make the conforming change for the composite nucleus.
      x = x - 133./(acom-aout)**1.3333
      t=twot(kp)*0.5
      esym(kc,kp,1) = x * (t-tz(kp)) * (t+tz(kp))
      esym(kc,kp,2) = x * (t+1.-tz(kp)) * (t+1.+tz(kp))
86    esym(kc,kp,3) = x * (t+2.-tz(kp)) * (t+2.+tz(kp))
      return            
      end   
c
c*********************************************************************
c
      subroutine sigpar
c 
c     Store parameters for calculating approximate optical model 
c     reaction cross sections.
c
c     called from: PRECOE
c
      common /par/ xl0(6), xl1(6), xm0(6), xm1(6) ,xn0(6), xn1(6),
     1 xn2(6), xp0(6), xp1(6), xp2(6)
c     ** Neutron parameters generated for the optimcal model potential
c     ** of Mani, Melkanoff and Iori.
      xp0(1) = -312.
      xp1(1) = 0.
      xp2(1) = 0.
      xl0(1) = 12.10
      xl1(1) = -11.27
      xm0(1) = 234.1
      xm1(1) = 38.26
      xn0(1) = 1.55
      xn1(1) = -106.1
      xn2(1) = 1280.8
c     ** Protons parameters for Becchetti-Greenlees potential (but with
c     ** a sub-barrier correction function and xp2 changed from -449)
      xp0(2) = 15.72
      xp1(2) = 9.65
      xp2(2) = -300.
      xl0(2) = 0.00437
      xl1(2) = -16.58
      xm0(2) = 244.7
      xm1(2) = 0.503
      xn0(2) = 273.1
      xn1(2) = -182.4
      xn2(2) = -1.872
c     ** Deuteron parameters for Perey and Perey potential
      xp0(3) = 0.798
      xp1(3) = 420.3
      xp2(3) = -1651.
      xl0(3) = 0.00619
      xl1(3) = -7.54
      xm0(3) = 583.5
      xm1(3) = 0.337
      xn0(3) = 421.8
      xn1(3) = -474.5
      xn2(3) = -3.592
c     ** Triton parameters for potential of Hafele, Flynn et al.
      xp0(4) = -21.45
      xp1(4) = 484.7
      xp2(4) = -1608.
      xl0(4) = 0.0186
      xl1(4) = -8.90
      xm0(4) = 686.3
      xm1(4) = 0.325
      xn0(4) = 368.9
      xn1(4) = -522.2
      xn2(4) = -4.998
c     ** He-3 parameters for potential of Gibson et al.
      xp0(5) = -2.88
      xp1(5) = 205.6
      xp2(5) = -1487.
      xl0(5) = 0.00459
      xl1(5) = -8.93
      xm0(5) = 611.2
      xm1(5) = 0.35
      xn0(5) = 473.8
      xn1(5) = -468.2
      xn2(5) = -2.225
c     ** alpha parameters for Huizenga and Igo potential
      xp0(6) = 10.95
      xp1(6) = -85.21
      xp2(6) = 1146.
      xl0(6) = 0.0643
      xl1(6) = -13.96
      xm0(6) = 781.2
      xm1(6) = 0.29
      xn0(6) = -304.7
      xn1(6) = -470.
      xn2(6) = -8.580
      return
      end
c
c*********************************************************************
c
      subroutine spec(kc,kp,multi,ihf,neprev,sigcni)
c
c     Calculate preeq energy spectra normalized to unit rxn strength
c     Store strength from primary emission for secondary emission calc
c
c     called from: PRECO
c
      common /emit/ clos(13,13), fiss(3), fweis, fweist, fteq, 
     1 mxdime, mxdimp, mxdim2, ndwn1, nup1, preeq1(6,51), preeq2(2,51),
     2 ptot, rho(13,13), ttot(6,51), weiss(6,51), weist(6,51)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      common /hauser/ csres0, csres1(9,52), csres2(3,52,52), csr(9),
     1 eres0, eres1(9,52), eres2(3,52,52)
      common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3),
     1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff
      common /wmain/ w(3,13,14,51), strong(9,10,51), seq(6,51),
     1 seq2(2,51), steq(6,51), steq2(2,51), xiso(9,10,51), xisoeq(2,51),
     2 spre(2,51)
c
c     Depletion from direct rxns is in clos(np,nppi) and fweis
c
      jout = jpout(kp) + jnout(kp)
      mkp = kp
      if (kp.gt.3) mkp=kp-3
      npmin = max0(ndwn1,jout+1)
      npilo1 = max0(jpin-ipout(kc),jpout(kp)) + 1
      do 36 ne=2,neps1(kc,kp)
      preeq = 0.
      eint = (eps(ne+1)-eps(ne-1)) * 0.5
      if (e.le.esym(kc,7,1)) go to 34
      do 32 np1=npmin,nup1
      mp1 = np1
      if (kp.gt.3) mp1 = mxdimp - np1 + 1
      np = np1-1
      npihi1 = np1 - max0(jnin-inout(kc),jnout(kp))
      do 32 nppi1=npilo1,npihi1
      mppi1 = nppi1
      if (kp.gt.3) mppi1 = mxdimp - nppi1 + 2
      x = w(mkp,mp1,mppi1,ne) * clos(np1,nppi1)
      preeq = preeq + x
      if (multi.gt.0) then
           if (kc.eq.1) then
c               ** Set strength available for secondary preeq emission
c               ** indexed for states in the residual nucleus and
c               ** normalized by the energy interval in primary emission.
                if (kp.eq.1) then
                    if (np1.le.mxdim2+1) then
                        if (iso.gt.0) x = x * xiso(np,nppi1,ne)
c                       ** xiso sets the strength going to the lowest T
c                       ** If another T is possible, strong is reset in
c                       ** REPENT before recycling.
                        strong(np,nppi1,ne) = strong(np,nppi1,ne)+x*eint
                    else
                        if (iso.gt.0) x = x * xisoeq(1,ne)
                        spre(1,ne) = spre(1,ne) + x*eint
c                       ** Preequilibrium strength from these states 
c                       ** with higher p is held for 2ndary equilibrium 
c                       ** emission.
                    end if
                else if (kp.eq.2) then
                    jp = mxdim2 - np + 1
                    jppi1 = mxdim2 - nppi1 + 3
                    if (np1.le.mxdim2+1) then
                        if (iso.gt.0) x = x * xiso(jp,jppi1,ne)
                        strong(jp,jppi1,ne) = strong(jp,jppi1,ne)+x*eint
                    else
                        if (iso.gt.0) x = x * xisoeq(2,ne)
                        spre(2,ne) = spre(2,ne) + x*eint
                    end if
                end if
           end if
      end if
32    continue
34    if (kc.eq.1) then
           preeq1(kp,ne) = preeq
           seq(kp,ne) = weiss(kp,ne) * fweis
           if (iso.gt.0) steq(kp,ne) = weist(kp,ne) * fweist
           if (ihf.eq.1) then 
                preeqc = preeq * sigcni * ct(1,7,1)
                csres0 = csres0 - preeqc*eint
                csres1(kp,ne) = csres1(kp,ne) + preeqc 
           else if (ihf.eq.2) then
                preeqc = preeq + seq(kp,ne) + steq(kp,ne)
                preeqc = preeqc * sigcni * ct(1,7,1)
                csres0 = csres0 - preeqc*eint
                csres1(kp,ne) = csres1(kp,ne) + preeqc 
           end if
c          ** Seq and steq are used to print the equilibrium components 
c          ** and to feed strength to secondary particle emission.
c          ** Division into two possible isospins is done in REPENT
      else
           preeq2(kp,ne) = preeq2(kp,ne) + preeq
           seq2(kp,ne) = seq2(kp,ne) + weiss(kp,ne) * fweis
           if (iso.gt.0) steq2(kp,ne)= steq2(kp,ne)+ weist(kp,ne)*fweist
           if (ihf.eq.1) then
                preeqc = preeq * sigcni * ct(1,7,1)
                kpp = kp + kc - 2
                csres1(kc-1,neprev) = csres1(kc-1,neprev) - preeqc
                if(csres1(kc-1,neprev).lt.0.0001) csres1(kc-1,neprev)=0.
                csres2(kpp,neprev,ne) = csres2(kpp,neprev,ne) + preeqc
           else if (ihf.eq.2) then
                wfw = weiss(kp,ne)*fweis
                if (iso.gt.0) wfw = wfw + weist(kp,ne)*fweist
                preeqc = (preeq+wfw) * sigcni * ct(1,7,1)
                kpp = kp + kc - 2
                x = preeqc * eint * 2. / (eps(neprev+1)-eps(neprev-1))
                csres1(kc-1,neprev) = csres1(kc-1,neprev) - x
                if(csres1(kc-1,neprev).lt.0.0001) csres1(kc-1,neprev)=0.
c               **above stmt supresses later printing of small residuals
                csres2(kpp,neprev,ne) = csres2(kpp,neprev,ne) + preeqc
           end if
c          ** fweis and fweist include strength from primary emission
      end if
36    continue
      return
      end
c
c*********************************************************************
c
      subroutine specout(kp,newpg,multi)
c
c     Separated from spec and revised 1995
c
c     Normalize and print energy spectra for
c           1. direct reactions (calcns called from PRECO)
c           2. preequilibrium emission
c           3. evaporation
c
c     called from: PRECO
c
      character*64 title
      character*1 newpg
      common /alph/ title
      common /angels/ elastf(51), elcoll(51), kang, tomsc(51),
     1 tomsd(51), xs(6)
      common /direct/ betalam(10), coul(7), elam(10), elast(51), fracel,
     1 lam(10), lamax, scoll(51), sigbar(7), sigcni, siggeom(51), 
     2 snock(6,51), snutra(6,51), wcoll(10)
      common /emit/ clos(13,13), fiss(3), fweis, fweist, fteq, 
     1 mxdime, mxdimp, mxdim2, ndwn1, nup1, preeq1(6,51), preeq2(2,51),
     2 ptot, rho(13,13), ttot(6,51), weiss(6,51), weist(6,51)
      common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3),
     1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6),
     2 nsigin, nphd(3,7,2), rzz, sigin(7,51)
      common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3),
     1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff
      common /wmain/ w(3,13,14,51), strong(9,10,51), seq(6,51),
     1 seq2(2,51), steq(6,51), steq2(2,51), xiso(9,10,51), xisoeq(2,51),
     2 spre(2,51)
      flow9 = 1.001e-09
      kin = jin + jpin
c     ------------------
c     1.  Print headings
c     ------------------
      write (iwri,4) newpg,title
4     format(a1,' ',a64)
      write (iwri,6) jpout(kp),jnout(kp)
6     format (' Particle Spectra (mb/MeV)  z=',1i2,'  n=',1i2)
      if (iso.gt.0) write(iwri,8) nt(1),esym(1,kp,1),esym(1,kp,2)
     1,ct(1,kp,1),ct(1,kp,2)
8     format (' T(cn)=Tz+',i1,'   Esym(b)=',f6.3,f8.3,'   Ct2(b)='
     1,f5.3,f7.3)
      write (iwri,22)
      multip = 0
      if (multi.gt.0) then
           if (kp.lt.3) multip = 1
      end if
c     ** Multip=1 means that there are secondary emission results to 
c     ** print.
      if (multip.gt.0) then
           if (nt(1).eq.0) then
                write (iwri,10)
10              format ('  eps   eps--     direct     --    exciton',
     1          5x,'--  equilibrium   --  sum')
           else
                write (iwri,12)
12              format ('  eps   eps--     direct     --    exciton',
     1          5x,'--  equilibrium   --  sum      tsum')
           end if
           write (iwri,14)
14         format ('   cm  chan-- nutra    knock --primary  second',
     1     ' --primary  second --')
      else
           if (nt(1).eq.0) then
                write (iwri,18)
18              format ('  eps   eps--     direct     --exciton--',
     1          ' equil --  sum')
           else
                write (iwri,20)
20              format ('  eps   eps--     direct     --exciton--',
     1          ' equil --  sum      tsum')
           end if
           write (iwri,24)
24         format ('   cm  chan-- nutra    knock --       --',
     1     7x,'--')
      end if
      write(iwri,22)
22    format(1h )
      sumnu = 0.
      sumke = 0.
      sumel = 0.
      sumpe = 0.
      sumpe2 = 0.
      sumw = 0.
      sumw2 = 0.
      jout = jnout(kp) + jpout(kp)
      chtocm = jout
      chtocm = 1. - chtocm/acom(1)
c     ------------------------------
c     2. Primary Emission
c     ------------------------------
      sigcni = sigcni * ct(1,7,1)
c     ** Depletion from direct rxns is in clos(np,nppi) and fweis.
      do 36 ne=2,neps1(1,kp)
      epscm = eps(ne) * chtocm
      pre1 = preeq1(kp,ne) * sigcni
      if (pre1.lt.flow9) pre1=0.
      evap1 = seq(kp,ne) * sigcni
      if (iso.gt.0) then
           evap1 = evap1 + steq(kp,ne)*sigcni
      end if
      if (evap1.lt.flow9) evap1=0.
c     -----------------------
c     3. Secondary Emission
c     -----------------------
      if (multip.gt.0) then
           pre2 = preeq2(kp,ne) * sigcni
           if (pre2.lt.flow9) pre2=0.
           evap2 = seq2(kp,ne) * sigcni
           if (iso.gt.0) then
                evap2 = evap2 + steq2(kp,ne)*sigcni
           end if
           if (evap2.lt.flow9) evap2=0.
      else
           pre2 = 0.
           evap2 = 0.
      end if
c     -------
c     4. Sums
c     -------
      tomsc(ne) = evap1 + evap2
      tomsd(ne) = pre1 + pre2 + snutra(kp,ne) + snock(kp,ne)
      tomsd(ne) = pre1 + pre2 + snutra(kp,ne) + snock(kp,ne)
      total = tomsd(ne) + tomsc(ne)
      if (kp.eq.kin) total = total + elast(ne) + scoll(ne)
      ttot(kp,ne) = ttot(kp,ne) + total
      if (ttot(kp,ne).le.0.) go to 36
      eint = (eps(ne+1) - eps(ne-1)) * 0.5
      sumnu = sumnu + snutra(kp,ne) * eint
      sumke = sumke + snock(kp,ne) * eint
      if (kp.eq.kin) then
           sumke = sumke + scoll(ne) * eint
           sumel = sumel + elast(ne) * eint
           snock(kp,ne) = snock(kp,ne) + elast(ne) + scoll(ne)
           elcoll(ne) = scoll(ne) + elast(ne) * fracel
           elastf(ne) = elast(ne) * (1.-fracel)
c          ** Snock contains angle-integrated spectra of elastic,
c          ** direct inelastic, and collective cross section.
c          ** Elcoll and elastf are used for printing angular dist'ns
c          ** where for complex projectiles, elastic and collective
c          ** processes are more forward peaked than other MSD strength.
      end if
      sumpe = sumpe + pre1 * eint
      sumw = sumw + evap1 * eint
      if (multip.gt.0) then
           sumpe2 = sumpe2 + pre2 * eint
           sumw2 = sumw2 + evap2 * eint
           if (nt(1).eq.0) then
                write(iwri,30) epscm, eps(ne), snutra(kp,ne),
     1          snock(kp,ne), pre1,pre2, evap1,evap2, total
30              format(f6.2,f5.1, 8(1pe9.2e1))
           else
                write(iwri,30) epscm, eps(ne), snutra(kp,ne),
     1          snock(kp,ne), pre1,pre2, evap1,evap2, total,ttot(kp,ne)
           end if
      else
           if (nt(1).eq.0) then
                write(iwri,30) epscm, eps(ne), snutra(kp,ne),
     1          snock(kp,ne), pre1, evap1, total
           else
                write(iwri,30) epscm, eps(ne), snutra(kp,ne),
     1          snock(kp,ne), pre1, evap1, total, ttot(kp,ne)
           end if
      end if
36    continue
      write(iwri,22)
      sumall = sumnu + sumke + sumpe + sumw
      if (multip.gt.0) then
           sumall = sumall + sumpe2 + sumw2
           write(iwri,38) sumnu,sumke, sumpe,sumpe2, sumw,sumw2, sumall
      else
           write(iwri,38) sumnu,sumke, sumpe, sumw, sumall
      end if
38    format (' sums',6x,7(1pe9.2e1))
      if (kp.eq.kin) then
           if (sumel.gt.0.) write (iwri,40) sumel
40         format (' elastic',12x,1pe9.2e1)
      end if
      if (kp.eq.1) then
           if (fiss(1).gt.0.) then
                fiss(1) = fiss(1) * sigcni
                if (multip.gt.0) then
                     fiss(2) = (fiss(2)+fiss(3)) * sigcni
                     write (iwri,42) fiss(1), fiss(2)
                else
                     write (iwri,44) fiss(1)
                end if
42              format (' fission',39x,2(1pe9.2e1))
44              format (' fission',30x,1pe9.2e1)
           end if
      end if
      sigcni = sigcni / ct(1,7,1)
      return
      end