c                                                       May 2006
c     PRECOM2.for contains the most interconnected subprograms
c          ALPHA                          NUTRA 
c          ASHELL                         OMEGA
c          COLLEC                         SINGLE
c          EMMY                           TCOR
c          FACT                           WELL
c     
c********************************************************************
c
      function alpha(ppi,hpi,pnu,hnu,ee,ze,gpi,gnu)
c     A simplified version of OMEGA for use in TCOR for calculating
c     the correction function to the state densities with good isospin
c
c     written in 1990 
c
c     Calculates 2-component particle-hole state densities in the
c     simple equi-spacing model
c     (no pairing, shell or finite well depth corrections).
c     
c     called from: TCOR
c
c
c     * I am alpha and omega, the beginning and the end.
c     * I will give unto him that is athirst of the fountain of the
c     * water of life freely.   (Rev. 21.6)
c
      common /state/ fac(52), zee
      pour = 1.e+35
      alpha = 0.
      ze = ee
      nppi = ppi
      nhpi = hpi
      npnu = pnu
      nhnu = hnu
      if (nppi.lt.0) go to 20
      if (nhpi.lt.0) go to 20
      if (npnu.lt.0) go to 20
      if (nhnu.lt.0) go to 20
      nppi1 = nppi + 1
      nhpi1 = nhpi + 1
      npnu1 = npnu + 1
      nhnu1 = nhnu + 1
      nph = nppi + nhpi + npnu + nhnu
      if (nph.lt.1) go to 20
      rho = gpi**ppi / fac(nppi1)
      rho = rho * gpi**hpi / fac(nhpi1)
      rho = rho * gnu**pnu / fac(npnu1)
      rho = rho * gnu**hnu / fac(nhnu1)
      ph = nph
c     --------------------------------------
c     1. Calculate pauli correction function
c     --------------------------------------
      qpi = amax1(ppi,hpi)
      epaul = ppi*(ppi+1.) + hpi*(hpi+1.)
      epaul = qpi*qpi - epaul*0.25
      epaul = epaul/gpi
      thresh = qpi*qpi/gpi
      qnu = amax1(pnu,hnu)
      epauli = pnu*(pnu+1.) + hnu*(hnu+1.)
      epauli = qnu*qnu - epauli*0.25
      epauli = epauli/gnu
      epauli = epauli + epaul
      thresh = qnu*qnu/gnu + thresh
      ze = ee - epauli
      edif = ee - thresh
      if (edif.gt.0.) then
           term3 = 0.
           act = 0.
           if (nppi.gt.0) act = 1.
           if (nhpi.gt.0) act = act+1.
           if (npnu.gt.0) act = act+1.
           if (nhnu.gt.0) act = act+1.
           act1 = 1./act
           act2 = act1 + act1
           act3 = act2 + act1
           if (qpi.gt.0.) then
                den = 5.7 + 0.6*act
                den = den * gpi * edif * act / ph
                den = 9.35 + den
                t3 = 0.
                if (nppi.gt.0) t3 = (ppi-act2)*(ppi-act3) + act1
                if (nhpi.gt.0) t3 = (hpi-act2)*(hpi-act3) + act1 + t3
                term3 = t3/(gpi*den)
           end if
           if (qnu.gt.0.) then
                den = 5.7 + 0.6*act
                den = den * gnu * edif * act / ph
                den = 9.35 + den
                t3 = 0
                if (npnu.gt.0) t3 = (pnu-act2)*(pnu-act3) + act1
                if (nhnu.gt.0) t3 = (hnu-act2)*(hnu-act3) + act1 + t3
                term3 = t3/(gnu*den) + term3
           end if
           ze = ze - term3
      end if
c     --------------------------
c     2. Calculate state density
c     --------------------------
      if (edif.le.0.) go to 20
      if (ze.le.0.) go to 20
      alpha = (ph-1.)*alog10(ze) - alog10(fac(nph))
      alpha = rho * 10.**alpha
      if (alpha.lt.0.1) alpha=0.
      if (alpha.gt.pour) alpha=0.
20    return
      end
c
c*********************************************************************
c
      function ashell(kp,api0,anu0,temp)
c     Shell corrected equilibrium level density parameters using simple 
c     staircase weighting functions.
c     Always evaluated with the maximum weight in the center of the
c     shell gap (as for closed shell nuclei).
c     For other nuclei, take a weighted average with a0.
c
c     Written in 1996
c
c     Called from EMMY
c
      dimension asp(2), asp0(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)
      asp0(1) = api0
      asp0(2) = anu0
      span = fspan * temp
      do 12 ia=1,2
      xmovf = range(ia)
c     ** Range is the number of mass units over which shell effects
c     ** wash out.
      idifa = iabs(idif(kp,ia))
      difa = idifa
      if (difa.ge.xmovf) then
           asp(ia) = asp0(ia)
           go to 12
      end if
      db2 = db(ia) * 0.5
      wtmax = span / ds(kp,ia)
      ninter = wtmax + 1.
c     ** With this ninter, the last interval will have a fractional weight
      sum = 0.
      wtsum = 0.
      weight = wtmax + 1.
      do 10 j=1,ninter
      xj = j
      ehigh = xj * ds(kp,ia)
      elow = ehigh - ds(kp,ia)
      weight = weight - 1.
      if (weight.le.0.) go to 10
      wtsum = wtsum + weight
      if (elow.gt.db2) then
           sum = sum + weight
c          ** Bins in this region have one s.p. state in them
      else if (elow.eq.db2) then
           sum = sum + (xj+1)*0.5*weight
c          ** Degen level is at low end of bin. Bin has 1 other s.p. st.
      else if (ehigh.eq.db2) then                  
           sum = sum + xj*0.5*weight
c          ** Degen level is at top of bin.  Add in half of it.
      else if (ehigh.gt.db2) then
           sum = sum + xj*weight
c          ** This bin has the degenerate level inside it; maybe 1 other
      end if
c          ** The cases where ehigh<db2 are in the gap and make no
c          ** contribution.
10    continue
      wtsum = wtsum * ds(kp,ia)
      asp(ia) = 1.645 * sum / wtsum
      if (idifa.gt.0) then
           fwash = sqrt(difa/xmovf)
           asp(ia) = asp0(ia)*fwash + asp(ia)*(1.-fwash)
      end if
12    continue
      ashell = asp(1) + asp(2)
      return
      end
c
c*********************************************************************
c
      function collec(kc,nppi,nhpi,npnu,nhnu,kp,in)
c
c     written in 1990
c
c     Calculate correction to pauli energy due to collective pairing
c     coll = ethresh(coll pair) - ethresh(esm + pair)
c
c     (epbar and ehbar are corrected for finite well depth effects)
c     (coll is to be subtracted from e in state density calcns)
c
c     called from: COMPOUND, OMEGA
c
      dimension e(2), n(2), nq(2)
      common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin,
     1 xncrit(7,2)
      common /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull
      coll = 0.
      ppi = nppi
      hpi = nhpi
      pnu = npnu
      hnu = nhnu
      n(1) = nppi + nhpi
      n(2) = npnu + nhnu
      nq(1) = max0(nppi,nhpi)
      nq(2) = max0(npnu,nhnu)
      e(1) = ppi*epbar + hpi*ehbar
      e(2) = pnu*epbar + hnu*ehbar
      do 20 ia=1,2
      if (cpre(1,kin,ia).le.0.) go to 20
      if (in.gt.0) then
           if (nq(ia).eq.2) then
                if (n(ia).eq.2) go to 20
           end if
      end if
c     ** Above tests bypass calculation of collective pairing 
c     ** contributions for direct transfer of a pair of nucleons of
c     ** type ia, since pairs will not be broken in this process.
      xn = n(ia)
      q = nq(ia)
      xnrat = xn / xncrit(kp,ia)
      ckp = cpre(1,kin,ia) * gres(kc,kp,ia) / gres(1,kin,ia)
c     ** Xncrit was calculated with the target pairing gap.
c     ** Ckp the is target pairing energy scaled by s.p. state densities.
c     ** Thus it is the pairing energy in the nucleus under study that
c     ** is based on the target delta0.
c     ** Delrat (below) is delta/delta0 in the active nucleus.
c     ** C is gres*delta*delta/4.  It is analagous to
c     ** ckp = c0 = gres*delta0*delta0/4.
c     ** For xnrat < 0.446, ephase is really the threshold energy.
c     ** For xnrat > 0.446, it is the energy of the phase transition.
c     ** In either case delta and c are set to zero for e(ia) < ephase.
c     ** This leads to a zero value for coll.
      if (xnrat.le.0.446) then
           ephase = ckp * (3.23*xnrat - 1.57*xnrat*xnrat)
      else
           ephase = ckp * (0.716 + 2.44*xnrat**2.17)
      end if
      if (e(ia).le.ephase) go to 20
      erat = e(ia) / ckp
      delrat = 0.996 - 1.76*(xnrat**1.6)/(erat**0.68)
      c = ckp * delrat * delrat
      x = q / gres(kc,kp,ia)
      y = x*x + 4.*c/gres(kc,kp,ia)
      coll = coll - c + q*sqrt(y) - q*x
20    continue
      collec = coll
      return
      end
c
c********************************************************************
c
      subroutine emmy(kc,kp,multi)
c
c     Written in 1982; revised 1990
c     Constant temperature eqb state densities at low E added 1996
c
c     Particle emission rates from particle-hole states
c     and general two-component equilibrium evaporation rates.
c     Includes sum over allowed exit channel isopins in
c     an isopin conserved calculation.
c     Constant temperature part assumed spin cutoff factor is
c     constant below the match point.
c
c     called from: PRECO
c     calls to: ASHELL, OMEGA, WELL
c
      dimension prod(51), ref(3), ur(51)
      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 /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)
      common /state/ fac(52), ze
      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 /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull
      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)
      flow = 1.e-18
      spill = 1.e+18
      umax = e - ben(kc,kp)
      if (umax.le.0.) go to 38
      mkp = kp
      if (kp.gt.3) mkp = kp-3
      jout = jpout(kp) + jnout(kp)
      iout = ipout(kc) + inout(kc)
      ap = jout
      apr = ap * (acom(kc)-ap) / acom(kc)
      sd = nsd(kp)
      rel = sd * apr * 3.72e-04
c     --------------------------------
c     1. equilibrium evaporation rates
c     --------------------------------
      eeh = e
c     ** EEH should be E-ROSE(7) but with 2ndary emiss, the 2nd CN can
c     ** have low E and need a constant temperature state density, 
c     ** which extends down to zero excitation energy.
c     ** Since composite state densities cancel out in equilibrium 
c     ** branching ratios, this makes no difference and the constant 
c     ** temperature form is not used there.  They are, however, used in 
c     ** the residual nucleus state densities.
c     ** But be sure that same CN state density is used in calculating 
c     ** the fission rates!
      api0 = 1.645 * gres(kc,kp,1)
      anu0 = 1.645 * gres(kc,kp,2)
      aspr0 = api0 + anu0
      asp = 1.645 * (gres(kc,7,1)+gres(kc,7,2))
      umatch = 2.5 + 150./(acom(kc)-ap)
      temp = sqrt(aspr0/umatch) - 1.25/umatch
      temp = 1. / temp
c     ** Matching energy determines transition to constant temp state
c     ** densities for residual nuclei in evaporation rates.
      if (ishell.gt.0) then
           aspr = ashell(kp,api0,anu0,temp)
           temp = sqrt(aspr/umatch) - 1.25/umatch
           temp = 1. / temp
      else
           aspr = aspr0
      end if
      const = 2.* sqrt(aspr*umatch)
      x = (aspr*umatch)**0.25 * umatch
      const = exp(const) / x
      x = exp((umatch+rose(kp))/temp)
      const = const / x
      ref(1) = 2. * sqrt(asp*eeh)
      imax=0
      if (iso.gt.0) then
           imax=2
           if (ct(kc,kp,2).le.0.) imax=1
           ee = eeh - esym(kc,7,1)
           ref(2) = 0.
           if (ee.gt.0.) ref(2) = 2. * sqrt(asp*ee)
           ref(3) = ref(2)
      end if
c     ** i=0 is for the T-mixed part of the evaporation rates.
c     ** i=1,2 are for the two resid. isospins in the T-conserved part
      do 10 ne=2,neps1(kc,kp)
      ee = eeh
      wrs = 0.
      wrt = 0.
      ur(ne) = umax - eps(ne)
      prod(ne) = eps(ne) * sigin(kp,ne) * rel
      do 6 i=0,imax
      u = ur(ne)
      if (i.eq.1) then
           ee = ee - esym(kc,7,1)
           if (ee.le.0.) go to 8
           u = u - esym(kc,kp,1)
      else if (i.eq.2) then
           u = u + esym(kc,kp,1) - esym(kc,kp,2)
      end if
      if (u.le.0.) go to 8
      u = u - rose(kp)
      if (u.ge.umatch) then
           quo = (u/ee)**1.25 * (aspr/asp)**0.25
           prex = prod(ne) / quo
           ex = 2. * sqrt(aspr*u)
      else
           u = u + rose(kp)
           prex = prod(ne) * ee**1.25 * asp**0.25
           prex = const * prex
           ex = u / temp
      end if
      ex = ex-ref(i+1)
      wr = prex * exp(ex)
      if (wr.lt.flow) wr=0.
      if (wr.gt.spill) wr=0.
      if (i.eq.0) then
           wrs = wr
      else
           x = wr*ct(kc,kp,i)
           wrt = wrt + x
c          ** If multiple preeq emission and T consv are both active,
c          ** then set fraction of primary emission going to lower T.
           if (multi.eq.0) go to 6
           if (kc.gt.1) go to 6
           if (kp.gt.2) go to 6
           if (i.eq.1) then
                xisoeq(kp,ne) = 1.
           else
                if (wrt.gt.0.) xisoeq(kp,ne) = 1. - x/wrt
           end if
      end if
6     continue
8     if (iso.gt.0) weist(kp,ne) = wrt
10    weiss(kp,ne) = wrs
c     --------------------------------
c     2. preequilibrium emission rates
c     --------------------------------
      if (kc.gt.1) then
           if (e.le.esym(kc,7,1)) go to 38
      end if
c     ** This bypasses 2ndary preeq. emission with T conservation
c     ** when energetically forbidden.
      npilo1 = max0(jpin-ipout(kc),jpout(kp)) + 1
      do 36 ne=2,neps1(kc,kp)
      u = ur(ne)
      i = 1
16    ee = ur(ne) - esym(kc,kp,i)
      if (ee.le.0.) go to 36
      do 32 np1=ndwn+1,nup1
c     ** np1 is the particle number index in the emitting nucleus.
c     ** np (below) is the particle number in the residual nucleus.
c     ** nh is the same in both nuclei, so nph is for the residual.
      mp1 = np1
      if (kp.gt.3) mp1 = mxdimp-np1+1
      if (iso.gt.0) then
c          ** Set indices for XISO, used in secondary particle emission.
c          ** Indices reflect configuration in primary residual states.
           jp1 = np1-1
           if (kp.eq.2) jp1 = mxdim2-np1+2
      end if
      np = np1-1-jout
      p = np
      if (np.lt.0) go to 32
      nh = np1-1 - jin+iout
      nph = np+nh
      if (nph.lt.1) go to 32
      npihi1 = np1 - max0(jnin-inout(kc),jnout(kp))
c     ** Now set indices for the packed array of emission rates.
      do 30 nppi1=npilo1,npihi1
      wr = 0.
      sig = 0.
      if (rho(np1,nppi1).le.0.) go to 30
      mppi1 = nppi1
      if (kp.gt.3) mppi1=mxdimp-nppi1+2
      if (iso.gt.0) then
           jppi1 = nppi1
           if (kp.eq.2) jppi1=mxdim2-nppi1+3
      end if
c     ** Calculate residual nucleus state density and emission rates.
      nppi = nppi1-jpout(kp)-1
      npnu = np - nppi
      nhpi = nppi1-1 - jpin + ipout(kc)
      nhnu = np1 - nppi1 - jnin + inout(kc)
      if (nh.eq.1) then
           v = ehole(nhnu+1)
      else if (nh.eq.2) then
           v = vfull
      end if
c     ** For the initial composite nucleus, allow auxilliary configurations
c     ** in the residual nucleus. These correpond to CN states where one
c     ** hole (or more than one), by chance, occupies a single particle
c     ** state next to the Fermi level and that, when the Fermi level
c     ** moves down during particle emission, is no longer a hole state 
c     ** but a particle state.  Thus the residual nucleus has one less 
c     ** degree of freedom than the "main" residual configurations.
      j2lim = nhpi
      j4lim = nhnu
      if (kc.eq.1) then
           j2lim = min0(nhpi,nppi)
           j4lim = min0(nhnu,npnu)
      end if
      do 29 j2 = j2lim,nhpi
      do 29 j4 = j4lim,nhnu
      nhj = j2 + j4
      h = nhj
      nphj = np + nhj
      rhol = 0.
      if (nphj.lt.1) go to 29
      ph = nphj
c     ** Calculate finite well depth corrections to the residual
c     ** nucleus state density.
      fwd(1) = 1.
      fwd(2) = 1.
      if (nhj.le.0) go to 28
      if (nh.gt.2) then
           go to 28
      else if (nh.eq.2) then
           if (ee.le.vfull) go to 28
      else
           if (ee.le.etest(nhnu+1)) go to 28
      end if
      fwd(1) = well(nphj,nhj,ee,v,vfull)
      fwd(2) = well(nphj+1,nhj,ee,v,vfull)
28    epbar = ee/ph
      ehbar = epbar
      if (fwd(1).gt.0.) epbar = epbar * fwd(2) / fwd(1)
      if (nhj.gt.0) ehbar = (ee - p*epbar) / h
      rhol = omega(kc,nppi,j2,npnu,j4,ee,kp,i,0) * ct(kc,kp,i) * fwd(1)
29    sig = sig + rhol
      wr = prod(ne) * sig / rho(np1,nppi1)
      if (wr.lt.flow) wr=0.
      if (wr.gt.spill) wr=0.
c     ** If multiple preeq emission and T consvervation are both active,
c     ** then set the fraction of primary emission going to the lower T.
      if (iso.eq.0) go to 30
      if (multi.eq.0) go to 30
      if (kc.gt.1) go to 30
      if (kp.gt.2) go to 30
      if (np1.gt.mxdim2+1) go to 30
      if (i.eq.1) then
           xiso(jp1,jppi1,ne) = 1.
      else
           x = wr + w(mkp,mp1,mppi1,ne)
           if (x.gt.0) xiso(jp1,jppi1,ne) = 1. - wr/x
      end if
30    w(mkp,mp1,mppi1,ne) = w(mkp,mp1,mppi1,ne) + wr
32    continue
      if (iso.le.0) go to 36
      if (i.gt.1) go to 36
      if (ct(kc,kp,2).le.0.) go to 36
      i=2
      go to 16
36    continue
c     ** Next stmts zero primary emission rates for states with nh=0
c     ** with A_a>1.  They are already zero for incident nucleons.
c     ** [Emission rates are calculated to allow generation of xiso
c     ** for secondary emission following direct stripping.]
c     ** To allow particle emission to begin at (p,h)=(Aa,0),
c     ** comment all stmts below up to (but NOT including) stmt 38.   
      if (kc.gt.1) go to 38
      if (jin.eq.1) go to 38
      do 34 ne=2,neps1(1,kp)
      mp1 = ndwn + 1
      mppi1 = npilo1
      if (kp.gt.3) then
           mp1 = mxdimp - mp1 + 1
           mppi1 = mxdimp - mppi1 + 2
      end if
34    w(mkp,mp1,mppi1,ne) = 0.
38    return
      end
c
c*******************************************************************
c
      function fact(x)
c
c     Calculates factorials using Sterling's Approximation with a
c     correction factor for better accuracy at low x
c
c     Corr'n factor changed (Feb 2007) from (12*x-1) to 12*x in denom.
c     This does slightly better at low x.
c
c     Called from: SINGLE
c
      pi = 3.14159265
      e  = 2.71828183
      fact = 1.
      if (x.le.1.) go to 10
      rat = (x/e)**x
      root = sqrt(2.*pi*x)
      corr = 1. + 1./(12*x)
      fact = rat * root * corr
10    return
      end
c
c*********************************************************************
c
      subroutine nutra(kp,multi,mxdim2)
c     
c     Written 1979; revised for shell effects 1982;
c     Revised again 1990, 2000 and 2005.
c
c     Phenomenological energy spectra for nucleon transfer.
c     Includes isospin conservation, collective pairing and finite 
c     well depth effects.
c
c     called from: PRECO
c     calls to: OMEGA, WELL
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 /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)
      common /state/ fac(52), ze
      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 /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull
      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
      spill = 1.e+18
      xin = jin
      kin = jin + jpin
      xpin = jpin
      jout = jpout(kp) + jnout(kp)
      xout = jout
      atarg = acom(1) - xin
      ecom = e - ben(1,7)
      elabi = ecom * acom(1) / atarg
      ares = acom(1) - xout
      umax = e - ben(1,kp)
      vab = 0.25 * 50. * xin
      xkab = 12.
      if (jin.eq.4) then
           if(ecom.gt.20.) xkab = 12. - 11.*(ecom-20.)/ecom
      end if
      if (jin.eq.2) xkab=1.
      if (jin.eq.3) xkab=1.
      if (jout.eq.2) xkab=1.
      if (jout.eq.3) xkab=1.
      if (jin.eq.jout) xkab=1.
c     ** xkab is an enhancement factor for (nucln,He4) and (He4,nucln)
      velo = xin / (elabi+vab)
      x = jout * nsd(kp) * xkab
      y = 6.
      if (jin.eq.1) y=2.
      if (jin.eq.4) y=4.
c     ** y is Aproj*(2s+1)
      xhold = 0.0125 * x * ct(1,7,1) / (y*ecom*xin)
      if (jin.gt.jout) xhold = xhold * sqrt(ecom) / 7.25
c      if (jin.eq.jout) xhold = xhold * sqrt(ecom) / 15.5 [WRONG may 06]
      if (jin.eq.jout) xhold = xhold * sqrt(ecom) / 14.5
c     ** trann (tranp) is number of transferred neutrons (protons)
      tranp = iabs(jpout(kp)-jpin)
      za = (rzz+rzz)/atarg
      trann = iabs(jnout(kp)-jnin)
      ph = trann + tranp
      nph = ph
      nex = 0
      if (nph.gt.0) go to 4
c     ** nex = 1(2) allows for exchange of a proton(neutron) in the
c     ** inelastic channel for projectiles with A = 2 to 4.
      nex = 1
      ph = 2.
      nph = 2
      tranp = 2.
4     jp = jpout(kp) - jpin
      jn = jnout(kp) - jnin
      jpickp = max0(jp,0)
      jstrip = jpickp - jp
      jpickn = max0(jn,0)
      jstrin = jpickn - jn
      nppi = jstrip
      if (nex.eq.1) nppi=1
      nhpi = jpickp
      if (nex.eq.1) nhpi=1
      npnu = jstrin
      if (nex.eq.2) npnu=1
      nhnu = jpickn
      if (nex.eq.2) nhnu=1
      x = xhold * (3800./ares)**nph
      if (kin.eq.1) x = xhold * (5500./ares)**nph
      x = x * velo**(2.*ph)
      pow = 2.*(2.+xpin)*nhpi + npnu+npnu
      x = x * za**pow
      v = ehole(1)
c     ** Currently the only projectile for which there is a different
c     ** ehole for p and n excitation is neutrons.  Pickup reactions
c     ** for incident neutrons must involve a picked up proton so
c     ** ehole(1) is the appropriate Veff to use. 
c     **
c     ** The factor 'reduce' (below) reduces the intensity of
c     ** configurations with an extra p-h pair in the residual nucleus.
c     ** Redux is set to this value when these configs are calculated.
      reduce = nppi*nppi + 1.5*nhpi*nhpi + npnu*npnu + nhnu*nhnu
      reduce = reduce * sqrt(elabi/xin)
      reduce = reduce * 7 / (v * atarg*atarg)
      if (nppi+nhpi.eq.0) v = v * za
c     ** Above stmt reduces the effective well depth for pure n transfer
c     ** when calculating finite well depth corrections in the
c     ** state densities.
      do 24 ne=2,neps1(1,kp)
      i = 1
      u = umax-eps(ne)
      sig = 0.
8     ee = u-esym(1,kp,i)
      if (ee.le.0.) go to 22
c     ** Next statements allow for extra configurations 
c     ** (base + 1p1h or 2p2h).
c     ** nexmax is the maximum number of extra p-h pairs allowed.
      nexmax = 3
      do 18 nexp = 0,nexmax
      kppi = nppi + nexp
      khpi = nhpi + nexp
      do 18 nexn = 0,nexmax-nexp
      kpnu = npnu + nexn
      khnu = nhnu + nexn
      kpart1 = kppi + kpnu + 1
      next = nexp + nexn
      redux = 1.
      if (next.gt.0) redux = reduce**next
c     ** Below are new loops added to allow for additional residual 
c     ** configurations corresponding to transfer at the Fermi level.
      if (next+nex.eq.0) then
           nppilo = 0
           nhpilo = 0
           npnulo = 0
           nhnulo = 0
      else
           nppilo = kppi
           nhpilo = khpi
           npnulo = kpnu
           nhnulo = khnu
      end if
      do 18 j1=nppilo,kppi
      do 18 j2=nhpilo,khpi
      do 18 j3=npnulo,kpnu
      do 18 j4=nhnulo,khnu
      np = j1 + j3
      nh = j2 + j4
      nphtemp = np + nh
      if (nphtemp.eq.0) go to 18
      p = np
      h = nh
      ph = nphtemp
c     ** End of creation of new loops.
      do 12 ia=1,2
      gpart(ia) = gres(1,kp,ia)
12    ghole(ia) = gres(1,kp,ia)
      fwd(1) = 1.
      fwd(2) = 1.
      if (nh.le.0) go to 14
c     ** Finite well depth correction results are evaluated with veff 
c     ** for the first interaction even though nh may be higher than 1.
      fwd(1) = well(nphtemp,nh,ee,v,vfull)
      fwd(2) = well(nphtemp+1,nh,ee,v,vfull)
14    epbar = ee/ph
      ehbar = epbar
      if (fwd(1).gt.0.) epbar=epbar*fwd(2)/fwd(1)
      if (nh.gt.0) ehbar = (ee-p*epbar)/h
      rhol = omega(1,j1,j2,j3,j4,ee,kp,i,1) * redux*fwd(1)*ct(1,kp,i)
      sig = rhol + sig
c     ** Next stmts store (complex,N) strength into strong(np1,nppi1,ne)
c     ** to make it available for secondary preeq emission.  Only 
c     ** strength to T< states is stored here.  
c     ** Strength to T> states (if any) is included in repent.
      if (multi.eq.0) go to 18
      if (kp.gt.2) go to 18
      if (i.gt.1) go to 18
      rhol = rhol * eps(ne) * sigin(kp,ne) * x 
      rhol = rhol / (sigcni*ct(1,7,1))
      eint = (eps(ne+1)-eps(ne-1)) * 0.5
      if (kp.eq.1) then
           strong(kpart1,kppi+1,ne)= strong(kpart1,kppi+1,ne)+rhol*eint
      else 
           np1 = mxdim2 - kpart1 + 1
           nppi1 = mxdim2 - kppi + 1
           strong(np1,nppi1,ne) = strong(np1,nppi1,ne) + rhol*eint
      end if
18    continue
c     ** Recyle to new residual isospin if appropriate.
      if (iso.le.0) go to 22
      if (i.gt.1) go to 22
      if (ct(1,kp,2).le.0.) go to 22
      i = 2
      go to 8
22    sig = sig * eps(ne) * sigin(kp,ne) * x
      if (sig.lt.flow9) sig=0.
      if (sig.gt.spill) sig=0.
24    snutra(kp,ne) = snutra(kp,ne) + sig
      if (kin.ne.kp) go to 26
c     ** Next stmts change exchanged pair for inelastic channel.
      if (nex.eq.2) go to 26
           nex = 2
           tranp = 0.
           trann = 2.
           ph = 2.
           go to 4
26    return
      end
c
c********************************************************************
c
      function omega(kc,nppi,nhpi,npnu,nhnu,ee,kp,it,in)
c
c     written in 1990 
c
c     Calculate 2-component particle-hole state densities in the
c     shell-shifted equi-spacing model
c           normal pairing is included in epauli
c           collective pairing corrections added here
c           isospin conservation corrections also added here
c           finite well depth corrections added in calling routine
c     
c     called from: COMPOUND, NUTRA, EMMY
c     calls to: COLLEC, TCOR, SINGLE
c
c
c     * I am alpha and omega, the beginning and the end.
c     * I will give unto him that is athirst of the fountain of the
c     * water of life freely.   (Rev. 21.6)
c
      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)
      common /state/ fac(52), ze
      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 /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull
      pour = 1.e+35
c     ** Set up particle and hole number related quantities 
      nppi1 = nppi+1
      nhpi1 = nhpi+1
      npnu1 = npnu+1
      nhnu1 = nhnu+1
      ppi = nppi
      hpi = nhpi
      pnu = npnu
      hnu = nhnu
      omega = 0.
      nph = nppi + nhpi + npnu + nhnu
      if (nph.lt.1) go to 20
      nq1 = max0(nppi1,nhpi1)
      nq2 = max0(npnu1,nhnu1)
      ph=nph
c     ** Main terms in the Pauli correction term
      coll = 0.
      if (icol.gt.0) coll = collec(kc,nppi,nhpi,npnu,nhnu,kp,in)
      adjust = 0.
      if (in.eq.0) go to 8
c     ** IN=1 means omega is called from nutra.  ADJUST (calculated
c     ** below) cancels out simple pairing corrections for direct
c     ** transfer (in nutra) of a pair of neutrons or protons, since
c     ** pairs are not broken and no pairing correction is needed in 
c     ** that case.
      if (nppi+nhpi.eq.2) then
           if (nppi.ne.1) adjust = cpre(kc,kp,1)
      end if
      if (npnu+nhnu.eq.2) then
           if (npnu.ne.1) adjust = adjust + cpre(kc,kp,2)
      end if
8     hold = thresh(kp,nq1,1) + thresh(kp,nq2,2) + coll - adjust
      if (ee.lt.hold) go to 20
c     ** Set single particle state densities
      if (ishell.gt.0) then
           call single(kc,nppi,nhpi,npnu,nhnu,ee,kp,coll,adjust)
      else
           do 10 ia=1,2
           gpart(ia) = gres(kc,kp,ia)
10         ghole(ia) = gres(kc,kp,ia)
      end if
c     ** Calculate non-energy dependent factors
c     ** ACT is the number of active classes of excitons.
      rho = 1.
      act=0.
      if (nppi1.gt.1) then
           rho = rho * (gpart(1)**nppi) / fac(nppi1)
           act = 1.
      end if
      if (nhpi1.gt.1) then
           rho = rho * (ghole(1)**nhpi) / fac(nhpi1)
           act = act+1.
      end if
      if (npnu1.gt.1) then
           rho = rho * (gpart(2)**npnu) / fac(npnu1)
           act = act+1.
      end if
      if (nhnu1.gt.1) then
           rho = rho * (ghole(2)**nhnu) / fac(nhnu1)
           act = act+1.
      end if
c     ** The above tests are needed in the S2-ESM at low excitation 
c     ** energies because the g's can be zero giving 0**0.
c     ** Now calculate the energy dependent part of the Pauli corr'n.
      gpi = gres(kc,kp,1)
      gnu = gres(kc,kp,2)
      qpi = nq1-1
      qnu = nq2-1
      edif = ee-hold
      term3 = 0.
      act1 = 1./act
      act2 = act1+act1
      act3 = act2+act1
      if(qpi.gt.0.) then
           den = 5.7 - 0.6*act
           den = den * gpi * edif * act / ph
           den = 9.35 + den
           t3 = 0.
           if (nppi.gt.0) t3 = (ppi-act2)*(ppi-act3) + act1
           if (nhpi.gt.0) t3 = (hpi-act2)*(hpi-act3) + act1 + t3
           term3 = t3/(gpi*den)
      end if
      if(qnu.gt.0.) then
           den = 5.7 - 0.6*act
           den = den * gnu * edif * act / ph
           den = 9.35 + den
           t3 = 0.
           if (npnu.gt.0) t3 = (pnu-act2)*(pnu-act3) + act1
           if (nhnu.gt.0) t3 = (hnu-act2)*(hnu-act3) + act1 + t3
           term3 = t3/(gnu*den) + term3
      end if
c     ** Now put it all together
      ze = ee - epauli(kp,nq1,1)-epauli(kp,nq2,2) - coll+adjust - term3
      if (ze.le.0) go to 20
      omega = (ph-1.)*alog10(ze) - alog10(fac(nph))
      omega = rho * 10.**omega
      if (omega.lt.0.1) omega=0.
      if (omega.gt.pour) omega=0.
      if (iso.gt.0) omega = omega * tcor(kc,ppi,hpi,pnu,hnu,ee,kp,it)
20    return
      end
c
c**********************************************************************
c
      subroutine single(kc,nppi,nhpi,npnu,nhnu,e,kp,coll,adjust)
c       
c     Effective single particle state densities in the S2-ESM
c
c     Written 1994 from stand-alone program shell.for
c          Replaces 1982 subroutine (revised in 1990,1992)
c
c     called from: COMPOUND and OMEGA
c     
      dimension del(2), fe(2), fs(2), ish(2), nq(2), np(2), nh(2),
     1 nhole(2), nhlo(2)
      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)
c     ------------------------------------------------------------
c     1. Initialize variables and calculate preliminary quantities
c     ------------------------------------------------------------
      do 170 ia = 1,2
      fe(ia) = 1.
      fs(ia) = 0.
      gpart(ia) = gres(kc,kp,ia)
170   ghole(ia) = gres(kc,kp,ia)
      npi = nppi + nhpi
      nnu = npnu + nhnu
      nph = npi+nnu
      ph = nph
      nqpi1 = max0(nppi,nhpi) + 1
      nqnu1 = max0(npnu,nhnu) + 1
      nq(1) = nqpi1 - 1
      nq(2) = nqnu1 - 1
      np(1) = nppi
      np(2) = npnu
      nh(1) = nhpi
      nh(2) = nhnu
      ihlotot=0
      ihtot=0
      ex = e - thresh(kp,nqpi1,1)-thresh(kp,nqnu1,2) - coll + adjust
      if (ex.lt.0.) go to 190
      exfree = ex
      phfree = ph
c     ** Exfree is the permutable excitation energy, i.e. that not 
c     ** tied up in exiting excitons across a shell gap. 
c     ** Phfree is the number of excitons not trapped between the Fermi
c     ** level and a shell gap.
c     ** Both are refined in the next loop.
      do 172 ia = 1,2
      nhole(ia) = nh(ia)
      nhlo(ia) = 0.
      del(ia) = 0.
      dbig = db(ia)
      dsmall = ds(kp,ia)
c     ** Now calculate the fraction of shell effects remaining.
c     ** If idifa.ge.range, shell effects are assumed to be washed out
      i = 0
      if (dbig.gt.dsmall) then
           xmovf = range(ia)
           idifa = abs(idif(kp,ia))
           difa = idifa
           if (difa.lt.xmovf) i=1
      end if
      ish(ia) = i
c     ** Ish(ia)=0 means no shell effects; 1 means shell effects active
      if (i.eq.0) go to 172
      if (nq(ia).le.0) go to 172
      fe(ia) = difa / xmovf
      fs(ia) = 1. - fe(ia)
c     ** Fs is the fraction of shell effects still present
c     ** Fe is the fraction washed out with distance from shell closure
      if (idif(kp,ia).gt.0) then
           ih = nh(ia)
      else
           ih = np(ia)
           nhole(ia) = ih
      end if
      iq = nq(ia)
      ihlo = min0(idifa,iq) - iq + ih
      if (ihlo.le.0) go to 172
      nhlo(ia) = ihlo
c     ** Ihlo and nhlo(ia) are the number of excitons of type ia
c     ** "trapped" between the Fermi level and a shell gap at threshold.
      ihlotot = ihlotot + ihlo
      ihtot = ihtot + ih
c     ** Ihlotot (& hlotot) are the total number of trapped excitons 
c     ** (neutron+proton) at threshold.
c     ** Ihtot is the total number of excitons on the shell gap side of 
c     ** the Fermi level when there are trapped excitons.
      hlo = ihlo
      dlim = (ex+dsmall*0.5) / dbig
      limdel = dlim
      limdel = amin0(limdel,ihlo)
      dlim = limdel
c     ** Dlim is the maximum number of holes which can be excited across 
c     ** the shell gap.
      exn = ex * hlo / ph
      deltg = amin1(dlim,exn/dbig)
      if (ih.lt.nph) then
           if (ih.gt.ihlo) then
                x = amin1(difa*dsmall/dbig,0.5)
                deltg = deltg * (1.-x)
           end if
      end if
      del(ia) = deltg
c     ** Deltg=del(ia) is the average number of holes excited across
c     ** the gap.
c     ** Now calculate corrections to exfree and phfree when excitons 
c     ** have been excited across a shell gap.
      exfree = exfree - deltg*dbig * fs(ia)
      phfree = phfree + (deltg-hlo) * fs(ia)
172   continue
      if (ish(1)+ish(2).le.0) go to 190
c     ** Ish(1)+ish(2)=0 means no shell effects active in this nucleus
      hlotot = ihlotot
      exph = 0.
      if (phfree.gt.0.002) exph = exfree / phfree
      esm = 0.
      esm1 = 0.
c     --------------------------------------------------------
c     2. Generate correction to N_{Ei} = g0exf for gap crossings
c     --------------------------------------------------------
c     ** g0exf = g_{i0} * (E-threshold) / (n-1); it is denoted N_{Ei} in
c     ** the code writeup where i is pi or nu for protons or neutrons.
c     ** Corrections are denoted Delta_{Na} and Delta_{Nb} in the 
c     ** writeup but are esm1 and esm-esm1 here.  Esm1 is stored 
c     ** separately, because it is used alone in evaluating the "highly 
c     ** constrained" cases.
      do 176 ia = 1,2
      if (nq(ia).le.0) go to 176
      ihlo = nhlo(ia)
      if (ihlo.le.0) go to 176
      if (nph.le.1) go to 176
      es1 = 0.
      esm2 = 0.
      hlo = ihlo
      ih = nhole(ia)
      h = ih
      ihhi = ih - ihlo
      hhi = ihhi
      hfree = hhi + del(ia)
      if (hfree.le.0.) go to 176
      dbig = db(ia)
      dlim = (ex+ds(kp,ia)*0.5) / dbig
c     ** Dlim is the maximum number of excitons that could cross 
c     ** the shell gap.
      scale = h/ph
      limesm = dlim * scale
      limesm = min0(limesm,ihlo)
c     ** Limesm is the maximum number scaled by this class' share 
c     ** of total permutable excitation energy
c     ** and limited by initial number available. 
      dlim = limesm
      x = 0.
      if (limesm.gt.0) then
           do 174 l=1,limesm
           y = l
174        x = x + 1./y
           x = x * dbig
      end if
      if (limesm.lt.ihlo) then
           y = (ex*scale - dlim*dbig) / (dlim+1.)
           x = x + y 
      end if
      if (ihlotot+ihhi.lt.nph) then
           if (ihlotot+ihlo.gt.3) then
                difa = abs(idif(kp,ia))
                y = 1. - difa*ds(kp,ia)*0.5/dbig
                scale = scale * y
                x = x * y
           end if
      end if
      x = x * gres(kc,kp,ia) / (hhi+1.)
      es1 = hlo * dbig / h
      es1 = amin1(es1,x)
      deltg = del(ia)
      if (hhi+deltg.gt.0.) then
           nd3 = mdfar(kp,ia)
           xnd3 = nd3
           x = xnd3 - hhi - deltg - 1.
           if (x.gt.0.) esm2 = 0.5 * x * (hlo-deltg) / hlo
           if (nd3.le.ihhi+1) esm2 = - amin1(ex/(hlo*dbig),1.)
      end if
      esm = esm + (es1+esm2) * scale * fs(ia)
      esm1 = esm1 + es1*scale * fs(ia)
176   continue
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     * Main loop to calculate xnum1 and xnum2 and the g's  *
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     -----------------------------
c     1. Set preliminary quantities
c     -----------------------------
      do 186 ia = 1,2
      g0 = gres(kc,kp,ia)
      if (nq(ia).le.0) go to 186
      dsmall = ds(kp,ia)
      ds2 = 0.5 * dsmall
      dbig = db(ia)
      db2 = 0.5 * dbig
      iq = nq(ia)
      q = iq
      ih = nhole(ia)
      ip = np(ia) + nh(ia) - ih
      p = ip
      h = ih
      iqp = iq - ip
      iqh = iq - ih
      idifa = abs(idif(kp,ia))
      difa = idifa
      mnear = mdnear(kp,ia)
      mfar = mdfar(kp,ia)
      dnear = mnear
      dfar = mfar
      dod = dodd(kp,ia)
c     ** Xnum1 (xnum2) is the number of s.p. states accessible on the
c     ** same (opposite) side of the gap as the Fermi level in the S2ESM
c     ** Nd1,nd2,nd3 are numbers of active states in the corresponding 
c     ** degenerate levels.
c     ** Xc1,xc2,xc3 are their contributions to xnum1 and xnum2.
c     ** Xc2 and xc3 both contribute to xnum2.
c     ** Eta1,eta2,eta3 are the average occupation numbers of the
c     ** degenerate levels.
c     ** Eta1t,eta2t,eta3t are the threshold occupation numbers, and
c     ** eta1l,eta2l,eta3l are the limiting values at high energy.  
c     ** Delt is the maximum number of excitons that can be excited
c     ** into or out of a degenerate level.
      if (nph.gt.1) then
c     ----------------------------------------------------------------
c     2. Accessible single particle state with shells washed out.  n>1
c     ----------------------------------------------------------------
           xeff = 1. / (ph-1.)
           g0ex = g0 * ex
           g0exf = g0ex * xeff
           xnum1e = p + g0exf - esm
           xnum2e = h + g0exf - esm
           if (ip.eq.1) xnum1e = xnum1e + esm*0.5
           if (ih.eq.1) xnum2e = xnum2e + esm*0.5
           if (xnum1e.lt.0.) xnum1e = 0.
           if (xnum2e.lt.0.) xnum2e = 0.
           if (ish(ia).le.0) then
                xnum1 = xnum1e
                xnum2 = xnum2e
                go to 182
           end if 
c     ------------------------------------------------------------------
c     3. Accessible S2ESM s.p. states on the Fermi level or particle 
c        side of gap for n>1  --  xnum1
c     ------------------------------------------------------------------
           nd1 = max0(mnear-idifa-iqp,0)
           xnd1 = nd1
           eta1t = amin0(ip,nd1)
           xnum1 = amax0(nd1,ip) + g0exf - esm
           xc1 = xnd1
           if (nd1.gt.1) then  
                exn = eta1t * exph
                g0exn = g0 * exn
                eta1l = amin1(p-1.,xnd1*0.5)
                delt = 0.
                if (ex.gt.dod-ds2) then
                     dif = amax0(ip-nd1,0)
                     b = dif + dod*g0 - 0.5
                     delt = 2.*g0exn + b*b
                     delt = sqrt(delt) - b
                     if (eta1l.lt.xnd1*0.5) then
                          if (ip.lt.nph) delt = delt*0.75
                     end if
                end if
                eta1 = amax1(eta1t-delt,eta1l)
                if (eta1.ge.1.) then
                     xcomb = fact(xnd1) / (fact(eta1) * fact(xnd1-eta1))
                     xcomb = xcomb**(1./eta1)
                     xc1 = eta1*xcomb
                     xnum1 = xnum1 - xnd1 + xc1
                else if (ip.eq.1) then
                     xnum1 = xnum1 - xnd1/(g0exf+ph)
                end if
           end if
           if (xnum1.lt.0.) xnum1=0.
           xnum1 = xnum1*fs(ia) + xnum1e*fe(ia)
c     -----------------------------------------------------------------
c     4. Accessible S2ESM s.p. states on the side of gap opposite the
c        Fermi level for n>1.  --  xnum2
c     -----------------------------------------------------------------
           if (ih.le.0) go to 182
           nd2 = 0
           nd3 = 0
           nd3t = 0
           eta2t = 0.
           eta3t = 0.
           dbe = dbig
           ihlo = nhlo(ia)
           hlo = ihlo
           ihhi = ih - ihlo
           hhi = ihhi
c          ** hhi (hlo) is num of "holes" above (below) gap at thresh.
c          -------------------------------------------------------------
c          4a. Limiting values for the occupation of degenerate levels
c          and contributions to xnum2 from any ESM states below the gap.
c          -------------------------------------------------------------
           if (iq.le.idifa-mnear) then
c               ** Last exciton at threshold is below 1st degenerate level
                xnum2 = h
                dbe = dod + (difa-dnear-1.-q) * dsmall
                dbe = amax1(dod,dbe) + dbig
c               ** dbe is an effective value of the shell gap including
c               ** effects from the odd spacing just outside it.
                etest1 = ex - dbe + dbig
                if (etest1.ge.-ds2) then
                     xnum2 = idifa - iqh
                     nd2 = mnear
                end if
                etest = etest1 - dbig
                if (etest.ge.-ds2) then
                     xnum2 = xnum2 + dfar
                     nd3 = mfar
                end if
           else if (iq.le.idifa) then
c               ** Last exciton at threshold is in 1st degenerate level
                xnum2 = idifa - iqh
                nd2 = min0(idifa-iqh,mnear)
                eta2t = ih - max0(idifa-iqh-mnear,0)
                etest = ex - dbig
                if (etest.ge.-ds2) then
                     xnum2 = xnum2 + dfar
                     nd3 = mfar
                end if
           else
c               ** Last exciton at threshold is across the shell gap
                xs = iq - idifa - mfar - 1
                if (xs.ge.0.) dbe = dbe + dod + xs*dsmall
                nd2 = min0(idifa-iqh,mnear)
                nd2 = max0(nd2,0)
                eta2t = nd2
                nd3 = mfar - max0(iqh-idifa,0)
                nd3t = nd3
                xnd3t = nd3t
                eta3t = iq - max0(idifa,iqh)
                eta3t = amin1(eta3t,xnd3t)
                xnum2 = idifa - iqh + max0(iq-idifa,mfar)
           end if
           xnd2 = nd2
           xnd3 = nd3
c          -----------------------------------------------------------
c          4b. Average occupation numbers of degenerate levels 2 and 3
c          -----------------------------------------------------------
           limdel = (ex+ds2) / dbe
           limdel = min0(limdel,ihlo)
           deltc = limdel
c          ** Deltc is the maximum number of excitons of this class 
c          ** that can be excited across the gap.  (used 2 places below)
           exn = ex * hlo  / ph
           deltg = amin1(deltc,exn/dbe)
c          ** Deltg calculated earlier assumed dbe = db(ia)
           if (ih.lt.nph) then
                if (ihhi.gt.0) then
                     x = amin1(difa*dsmall/dbig,0.5)
                     deltg = deltg * (1.-x)
                end if
           end if
           eta3l = amin1(xnd3*0.5,hhi+deltg-1.)
           eta3l = amax1(eta3l,0.)
           eta2l = amin1(h-1.-eta3l,xnd2*0.5)
           deltg = amin1(deltg,hlo-eta2l)
           xc3 = 0.
           eta3 = 0.
           if (nd3.gt.1) then
                eta3 = amin1(eta3t+deltg,xnd3)
                x = eta3t+deltg - eta3
                if (eta3.gt.eta3l) then
                     delt = 0.
                     b = 0.
                     if (ex.gt.dod-ds2) then
                          dif = amax0(iq-idifa-mfar,0)
                          exn = (eta3t+deltg) * exph
c                         ** exph assumes that deltg is not reduced 
c                         ** for eta2l.
                          exn = exn - dod*x - ds2*x*(dif+dif-1.+x)
                          exn = amax1(exn,0.)
                          b = dif+x + dod*g0 - 0.5
                          delt = 2. * g0 * exn + b*b
                          delt = sqrt(delt) - b
                          if (eta3l.lt.xnd3*0.5) then
                                if (ih.lt.nph) delt = delt*0.75
                          end if
                     end if
                     eta3 = amax1(eta3-delt,eta3l)
                end if
           end if
           xc2 = 0.
           eta2 = 0.
           if (nd2.gt.1) then
                eta2 = eta2t - deltg
                if (eta2t.lt.eta2l) then
c                    ** all "holes" are below the gap at threshold
c                    ** add'l holes available below lower degen level
                     delt = 0.
                     if (ex.ge.dod-ds2) then
                          g0exn = g0 * (h-eta2t) * exph
                          b = g0*dod - 0.5
                          delt = sqrt(2.*g0exn+b*b) -b
                     end if
                     eta2 = amin1(eta2+delt,eta2l)
                else
                     eta2 = amax1(eta2,eta2l)
                end if
           end if
c          -----------------------------------------------------
c          4c. Contributions from the degenerate levels to xnum2
c          -----------------------------------------------------
           if (nd3.gt.0) then
                xc3 = xnd3
                if (eta3.ge.1.) then
                     xcomb = fact(xnd3) / (fact(eta3)*fact(xnd3-eta3))
                     xcomb = xcomb**(1./eta3)
                     xc3 = eta3*xcomb
                end if
                if (hhi+deltg.lt.1.) then
                     xc3 = xnd3 * deltg
                     if (ex.gt.dbe+ds2) xc3 = xc3 - xc3/(g0exf+ph)
                else if (hhi+deltg.lt.2) then
                     if (ex.ge.(1.-hhi)*dbe+ds2) then
                          xc3 = xc3 - xc3*(2.-hhi-deltg)/(g0exf+ph)
                     end if
                end if
                xnum2 = xnum2 - xnd3 + xc3
           end if
           if (nd2.gt.0) then
                if (eta2.gt.1.) then
                     xcomb = fact(xnd2) / (fact(eta2)*fact(xnd2-eta2))
                     xcomb = xcomb**(1./eta2)
                     xc2 = eta2*xcomb
                     xnum2 = xnum2 - xnd2 + xc2
                else
                     xc2 = xnd2
                end if
c               ------------------------------------------------
c               4d. Including effects of "Constraints" when most
c               excitons are trapped by gaps.
c               ------------------------------------------------
                lcns1 = 0
                lcns2 = 0
c               ** lcns1 and lcns2 refer to two different types or 
c               ** levels of constraint.
                limc = (ex-ds2)/dbe
                limdel = (ex+ds2)/dbe
c               ** Limc is limdel shifted up one s.p. spacing in energy
                if (ihhi+limc.eq.0) lcns1 = 1
c               ** Lcns1=1 means there are no free excitons in this class
                if (ihtot.eq.nph) then
                     if (ihhi.lt.2) then
                          jhhi = nhole(3-ia) - nhlo(3-ia)
                          if (jhhi.lt.2) then
                               if (ihhi+jhhi.eq.2) limc=limdel
                               jtest = 2
                               if (ih.lt.nph) jtest=4
                               jtest = min0(jtest,ih+1-jhhi)
                               if (ihhi+limc.lt.jtest) lcns2 = 1
                          end if
                     end if
                end if
c               ** Lcns2=1 means most excitons are trapped.
                lcnst = lcns1 + lcns2
c               ** Lcnst=0 means there are no constraints
                if (lcnst.le.0) go to 180
                if (ex.le.ds2) go to 180
                if (limc.eq.limdel) then
                     if (lcns1.eq.0) xnum2 = xnum2 - (g0exf-esm1)*0.5
                     xnum2 = xnum2 - xc2 - xc3
                     if (lcnst.eq.2) then
                          xc2 = xc2 * 0.5
                          if (ihlotot.eq.nph) xc2 = 0.
                     else
                          ext = ex - (deltc-hhi)*dbe
                          denom = dbe + hhi*dbe
                          if (limc+ihhi.gt.1) then
                               if (limc+ihhi.eq.jtest-1) then
                                    ext = ext + dbe
                                    denom = (denom+dbe)*2.
                               end if
                          end if
                          x = 1. - 0.25*ext/denom
                          xc2 = xc2 * x
                          xc3 = xc3 * x
                     end if
                     xnum2 = xnum2 + xc2 + xc3
                else if (ihtot.eq.nph) then
                     cns = lcnst
                     if(ihlotot.gt.1) cns= cns/(hlotot*(hlotot-1.))
                     xnum2 = xnum2 + (xc2+xc3)*cns
                     if (lcnst.eq.2) then
                          xnum2 = xnum2 - (g0exf-esm1)*0.5
                     else if (lcns2.eq.1) then
                          if (ihhi+limdel.lt.jtest) then
                               xnum2= xnum2- (g0exf-esm1)*0.5
                          end if
                     end if
                end if
180        end if
           if (ihhi+limdel.gt.0) xnum2 = xnum2 + g0exf - esm
           if (xnum2.lt.0.) xnum2=0.
           xnum2 = xnum2*fs(ia) + xnum2e*fe(ia)
c     -----------------------------------------------------------
c     5. Convert xnum1 and xnum2 into effective g values; for n>1
c     -----------------------------------------------------------
182        g0exf = g0ex * xeff
           if (ish(ia).le.0) then
                if (ip.gt.0) gpart(ia) = g0 * xnum1 / (p+g0exf)
                if (ih.gt.0) ghole(ia) = g0 * xnum2 / (h+g0exf)
           else if (idif(kp,ia).gt.0) then
c               **i.e. if Z>Zmag or N>Nmag
                if (ip.gt.0) gpart(ia) = g0 * xnum1 / (p+g0exf)
                if (ih.gt.0) ghole(ia) = g0 * xnum2 / (h+g0exf)
           else                                              
                if (ih.gt.0) gpart(ia) = g0 * xnum2 / (h+g0exf)
                if (ip.gt.0) ghole(ia) = g0 * xnum1 / (p+g0exf)
           end if
      else
c     ----------------------------------------------
c     6. Calculate xnum1 or xnum2 for case where n=1
c     ----------------------------------------------
           if (ish(ia).le.0) go to 186
           if (idifa.lt.mnear) then
                if (ex.lt.dod-ds2) then
                     xnum1 = dnear - difa
                else if (ex.lt.ds2) then
                     xnum1 = dnear - difa + 1.
                else
                     xnum1 = 1.
                end if
           else
                xnum1 = 1.
           end if
           xnum1 = xnum1*fs(ia) + fe(ia)
           if (idifa.le.0) then
                xnum2 = xnum1
           else if (idifa.le.mnear) then
                if (ex.le.ds2) then
                     xnum2 = difa
                else if (ex.lt.db2+db2-ds2) then
                     xnum2 = 0.
                else if (ex.le.db2+db2+dod-ds2) then
                     xnum2 = dfar
                else if (ex.le.db2+db2+ds2) then
                     xnum2 = dfar + 1.
                else
                     xnum2 = 1.
                end if
           else
                gs = difa * dsmall
                if (ex.lt.gs-dsmall-db2) then
                     xnum2 = 1.
                else if (ex.le.gs-db2-dod) then
                     xnum2 = dnear + 1.
                else if (ex.le.gs-db2) then
                     xnum2 = dnear
                else if (ex.lt.gs-dsmall+db2) then
                     xnum2 = 0.
                else if (ex.lt.gs-dsmall+db2+dod) then
                     xnum2 = dfar
                else if (ex.le.gs+db2+dod) then
                     xnum2 = dfar + 1.
                else
                     xnum2 = 1.
                end if
           end if
           xnum2 = xnum2*fs(ia) + fe(ia)
c     -----------------------------------------------------
c     7. Convert xnum1, xnum2 to effective g value when n=1
c     -----------------------------------------------------
           if (idif(kp,ia).gt.0.) then
                if (ip.gt.0) gpart(ia) = g0*xnum1
                if (ih.gt.0) ghole(ia) = g0*xnum2
           else
                if (ih.gt.0) gpart(ia) = g0*xnum2
                if (ip.gt.0) ghole(ia) = g0*xnum1
           end if
      end if
186   continue
190   return
      end
c
c*******************************************************************
c
      function tcor(kc,tppi,thpi,tpnu,thnu,ee,kp,i)
c
c     written in 1990; revised June 1991
c     
c     MAIN TERMS use type 2 blocking and state densities with their
c     natural particle and hole numbers 
c     CORRECTION TERMS have no blocking and one state density
c
c     called form: COMPOUND, OMEGA
c     calls to: ALPHA
c
      common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin,
     1 xncrit(7,2)
      common /state/ fac(52), zee
      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
      flow = 1.e-18
      spill = 1.e+18
      tz2 = tze(kp)*2.
      e = ee + esym(kc,kp,i) - esym(kc,kp,i+1)
c     ----------------------------------------------------
c     Set particle and hole numbers, 
c     reversing the roles of protons and neutrons if Tz<0.
c     ----------------------------------------------------
      if (itc(i).gt.0) then
           ppi = tppi
           hpi = thpi
           pnu = tpnu
           hnu = thnu
           gpi = gres(kc,kp,1)
           gnu = gres(kc,kp,2)
      else
           ppi = tpnu
           hpi = thnu
           pnu = tppi
           hnu = thpi
           gpi = gres(kc,kp,2)
           gnu = gres(kc,kp,1)
      end if
      ze = 0.
      denom = alpha(ppi,hpi,pnu,hnu,ee,ze,gpi,gnu)
c     ** Alpha returns a correct ESM value for ze=E-A(p,ppi)
      gpir = gpi/gnu
      gnur = gnu/gpi
      b = ppi*gnur + hnu*gpir
      cm = pnu + hpi
      c = cm + tz2
c     ** The index i value transferred into TCOR is T-Tmin+1, 
c     ** but TCOR needs T-abs(Tz)+1.  Call this "it".
c     ** Here twot is 2Tmin
      delt = 0.5 * twot(kp) - abs(tz(kp))
      it = delt+0.5
      it = it + i
c     ----------------------------------------------------
c     Calculate blocking fraction for type 2 isospin flips
c     ----------------------------------------------------
      xi = it -1
      ph = ppi + hpi + pnu + hnu
      sym = (tze(kp)+xi) * 4. / (gpi+gnu) 
      zeb = ze - sym
      if (zeb.gt.0.) then
           f = (zeb/ze)**(ph-1.)
           if (f.lt.flow) f=0.
           if (f.gt.spill) f=0.
           f = 1.-f 
      else
           f=1.
      end if
c     -----------------------------------------------------
c     Calculate state densities according to T-Tz of states
c     -----------------------------------------------------
      zet = ze + esym(kc,kp,i) - esym(kc,kp,i+1)
      if (zet.gt.0.) fac2 = (ph-1.) * (ph-2.) / (gpi*gnu * zet*zet)
      t2 = 0.
      if (it.eq.1) then
           t1 = denom
           if (zet.gt.0.) then
                bt = tz2+2.
                x1 = b / (c+3.)
                x2 = bt * ppi*hnu * fac2 / (c+2.)
                t2 = x1 + x2
           end if
      else if (it.eq.2) then
           bt = tz2+2.
           bp0 = f*ppi
           bp1 = bp0 - f
           if (bp0.gt.bt) bp0=bt
           if (bp1.gt.bt) bp1=bt
           bh0 = f*hnu
           bh1 = bh0 - f
           if (bh0.gt.bt) bh0=bt
           if (bh1.gt.bt) bh1=bt
           b1 = bt - bp1 - bh0 + bp1*bh0/bt
           b2 = bt - bp0 - bh0 + bp0*bh0/bt
           b3 = bt - bp0 - bh1 + bp0*bh1/bt
           x1 = (pnu+1.) / (b1+cm+1.)
           x2 = 0.
           if (b2.gt.0.) x2 = b2 / (b2+cm)
           x3 = (hpi+1.) / (b3+cm+1.)
           t1 = x1 * alpha(ppi-1.,hpi,pnu+1.,hnu,ee,ze,gpi,gnu)
           t1 = t1 + x2 * denom
           t1 = t1 + x3 * alpha(ppi,hpi+1.,pnu,hnu-1.,ee,ze,gpi,gnu)
           if (zet.gt.0.) then
                x1 = x1 * ppi*gnur / (pnu+1.)
                x3 = x3 * hnu*gpir / (hpi+1.)
                bt = tz2+4.
                x11 =  x1 * (b-gnur) / (c+6.)
                x12 = x1 * bt * (ppi-1.)*hnu * fac2 / (c+5.)
                x21 = x2 * b / (c+5.)
                x22 = x2 * bt * ppi*hnu * fac2 / (c+4.)
                x31 = x3 * (b-gpir) / (c+6.)
                x32 = x3 * bt * ppi*(hnu-1.) * fac2 / (c+5.)
                t2 = x11 + x12 + x21 + x22 + x31 + x32
           end if
      else if (it.eq.3) then
           bt = tz2+4.
           bp0 = f*ppi
           bp1 = bp0 - f
           bp2 = bp1 - f
           if (bp0.gt.bt) bp0=bt
           if (bp1.gt.bt) bp1=bt
           if (bp2.gt.bt) bp2=bt
           bh0 = f*hnu
           bh1 = bh0 - f
           bh2 = bh1 - f
           if (bh0.gt.bt) bh0=bt
           if (bh1.gt.bt) bh1=bt
           if (bh2.gt.bt) bh2=bt
           b11 = bt - bp2 - bh0 + bp2*bh0/bt
           b12 = bt - bp1 - bh0 + bp1*bh0/bt
           b13 = bt - bp1 - bh1 + bp1*bh1/bt
           b22 = bt - bp0 - bh0 + bp0*bh0/bt
           b23 = bt - bp0 - bh1 + bp0*bh1/bt
           b33 = bt - bp0 - bh2 + bp0*bh2/bt
           den = b11+cm+2.
           x11 = (pnu+1.) * (pnu+2.) / (den * (den-1.))
           den = b12+cm+1.
           x12 = (pnu+1.) * b12 / (den * (den-1.))
           den = b13+cm+2.
           x13 = (pnu+1.) * (hpi+1.) / (den * (den-1.))
           den = b23+cm+1.
           x23 = b23 * (hpi+1.) / (den * (den-1.))
           den = b33+cm+2.
           x33 = (hpi+1.) * (hpi+2.) / (den * (den-1.))
           den = b22+cm
           x22 = 0.
           if (den.gt.0.) x22 = b22 * (b22-1.) / (den * (den-1.))
           t1 = x11 * alpha(ppi-2.,hpi,pnu+2.,hnu,ee,ze,gpi,gnu)
           y = alpha(ppi-1.,hpi,pnu+1.,hnu,ee,ze,gpi,gnu) 
           t1 = t1 + (x12+x12) * y
           y = alpha(ppi-1.,hpi+1.,pnu+1.,hnu-1.,ee,ze,gpi,gnu)
           t1 = t1 + (x13+x13) * y
           t1 = t1 + x22 * denom
           y = alpha(ppi,hpi+1.,pnu,hnu-1.,ee,ze,gpi,gnu)
           t1 = t1 + (x23+x23) * y
           y = alpha(ppi,hpi+2.,pnu,hnu-2.,ee,ze,gpi,gnu)
           t1 = t1 + x33 * y
           if (zet.gt.0.) then
                x11 = x11 * ppi*(ppi-1.)*gnur*gnur / ((pnu+1.)*(pnu+2.))
                x13 = x13 * ppi * hnu / ((pnu+1.)*(pnu+2.))
                x33 = x33 * hnu*(hnu-1.)*gpir*gpir / ((hpi+1.)*(hpi+2.))
                x12 = x12 * ppi*gnur / (pnu+1.)
                x23 = x23 * hnu*gpir / (hpi+1.)
                bt = tz2+6.
                x111 = x11 * (ppi-2.)*gnur / (c+9.)
                x113 = x11 * hnu*gpir / (c+9.)
                x133 = x13 * (hnu-1.)*gpir / (c+9.)
                x333 = x33 * (hnu-2.)*gpir / (c+9.)
                x = bt * fac2 / (c+8.)
                x112 = x11 * (ppi-2.)*hnu * x
                x132 = x13 * (ppi-1.)*(hnu-1.) * x
                x332 = x33 * ppi * (hnu-2.) * x
                x121 = x12 * (ppi-1.)*gnur / (c+8.)
                x123 = x12 * hnu*gpir / (c+8.)
                x233 = x23 * (hnu-1.)*gpir / (c+8.)
                x122 = x12 * bt * (ppi-1.)*hnu * fac2 / (c+7.)
                x232 = x23 * bt * ppi*(hnu-1.) * fac2 / (c+7.)
                x221 = x22 * b / (c+7.)
                x222 = x22 * bt * ppi*hnu * fac2 / (c+6.)
                t2 = x111 + x333 + x112 + x332 + x221 + x222
                t2 = t2 + 2.*(x121 + x122 + x132 + x232 + x233)
                t2 = t2 + 3.*(x113 + x133) + 4.*x123
           end if
      end if
      tcor = t1
      if (zet.gt.0.) then
      tcor = tcor - t2 * alpha(ppi,hpi,pnu,hnu,e,ze,gpi,gnu)
      end if
      if (tcor.lt.0.) tcor=0.
      if (denom.gt.0.) then
      tcor = tcor / denom
      else
      tcor = 1.
      end if
      return
      end
c
c*********************************************************************
c
      function well(n,nh,ee,v,vfull)
c
c     Finite well depth correction functions for state densities
c     This includes averaging when an effective well depth due to 
c     surface localization of the initial interaction is less than
c     the central well depth.
c
c     called from: PRECO, NUTRA, EMMY
c
c     * Therefore with joy shall ye draw waters out of the wells
c     * of salvation. (Isaiah 12.3)
c
c     * But whosoever drinketh of the water that I shall give him
c     * shall never thirst, but the water that I shall give him shall
c     * be in him a well of water springing up into everlasting life.
c     * (John 4.14)
c
      xx = 1.
      w = v * (vfull-v)/(vfull+vfull)
      if (n.eq.1) go to 8
c     ------------------------------------------------------------
c     For n>1 calculate the approximate average of the finite well 
c     depth correction over a range of well depths if needed.
c     ------------------------------------------------------------
      jmax = 5
      if (v.ge.vfull-0.5) jmax=1
      wtsum = 0.
      fwtsum = 0.
      ph1 = n-1
      h = nh
      do 7 j=1,jmax
      xj = j-1
      wt = (1.+exp(xj))*(1.+exp(-xj))
      wt = 1./wt + 0.001*xj
2     vj = v + xj*w
      if (vj.gt.vfull) go to 6
      if (vj.lt.0.) go to 6
      f = 1.
      term = 1.
      xsign = 1.
c     ** Well is only called for h=1,2 in the exciton model,
c     ** but for 3 nucleon pickup or auxilliary configurations 
c     ** in nutra, it can be called for higher n.
c     ** This programming is general.  
      do 3 ih = 1,nh
      hi = ih
      ev = ee - hi*vj
      if (ev.le.0.01) go to 4
      term = term * (h-hi+1.) / hi
      xsign = -xsign
      ev = alog10(term) + ph1*alog10(ev/ee)
      if (ev.lt.-5) go to 3
      f = f + xsign * 10.**ev
3     continue
4     wtsum=wtsum+wt
      fwtsum=fwtsum+wt*f
6     if (xj.le.0.) go to 7
      xj=-xj
      go to 2
7     continue
      if (wtsum.gt.0.) xx=fwtsum/wtsum
      go to 10
c     ------------------------------------------
c     Use a closed form average of fwd for n=h=1
c     ------------------------------------------
8     xx=0.
      if (ee.gt.1.16*vfull) go to 10
      if (ee.lt.v) xx=1.
      if (v.ge.vfull) go to 10
      xx=exp((ee-v)/w)
      xx=1./(1.+xx)
10    well=xx
      return
      end

