      Subroutine tfold
c *** see RH Landau, Program lpott, 1981
c     tfold folds t matrices for this nucleus and stores them in retij,
c     imtij.nb,nbnbnbnbn, cant get to unfolded t after folding
c      nifty(5)=6   fold t and denom with e3body,only Realdenom
c         n.b. different folded amps for e regular and e3body
      Implicit Real*8(a-h,i,k,m,o-z)
      Dimension ptp(24), ptcos(16), wtp(24), wtcos(16)
      Dimension nifty(20), bnuc(20), bnucf(20), bn(16), bnf(16), retij
     1(14), imtij(14), rt(14), it(14), rtcos(14), itcos(14), kapd(100),
     2rewave(14,100), imwave(14,100), d(14,76), eta(14,76), gin(14,100),
     3 g(14), kapg(100), denom(14,25), ecmden(25)
      Common /nlspfl/ gin,kapg,denom,ecmden,fmn,rewave,imwave,kapd,
     1nkapgs,necms
      Common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,mn,nz,nes
     1,nwaves,nifty,na
      ma = na*mn
      aovera = (na-1.)/na
      twomu = 2.*(mpi+mn)*(na-1.)*mn/(mpi+ma)
      mn2 = mn*mn
      mpi2 = mpi*mpi
      rhl = mn2+mpi2
      rhl4 = (mpi+mn)**2
      rhl3 = (mpi-mn)**2
c     note% factor of 2 in rho from costh intgen,?a? indep of nucleus
      a = (1.66/hbarc)*(na/12.)**(.333333)
      If (na.eq.4) a = 1.4/hbarc
      If (na.eq.3) a = 1.4/hbarc
      anorm = 2.*pi*(a**3)/pi**1.5/na
c     set up points and wts for integration over p and cos theta
      nps = 12
      ncoss = 6
      Call Gauss (nps,021,300.d0,0.d0,ptp,wtp)
      Call Gauss (ncoss,011,-1.d0,1.d0,ptcos,wtcos)
      If (na.le.4) nna = 1
      If ((na.gt.4).and.(na.le.16)) nna = 2
      If ((na.gt.16).and.(na.le.40)) nna = 3
      If (na.gt.40) nna = 4
c     Do loop over tabulated energies
      Read (3,310) nstore
      Rewind 3
      nfold = 60
      If (nfold.lt.nstore) GoTo 10
      Write (6,290) nstore
      nfold = nstore
 10   one = 0.
      nkape = nfold
c       also fold denom If e3b is used
      If (nifty(5).gt.4) nkape = nfold+necms
      Do 230 nkap=1,nkape
         If (nkap.gt.nfold) GoTo 20
c       non-denom
         kap = kapd(nkap)
         kap2 = kap*kap
         plab = (kap/mn)*(Sqrt(mn2+kap2)+Sqrt(mpi2+kap2))
         rhl1 = 2.*Sqrt(mpi2+plab*plab)
         If (nifty(5).lt.5) GoTo 20
c      calc folding and k0 for e3b,use static-like approx to k0
         k0 = mpi2+ma*ma+rhl1*ma
         k02 = (k0-(mpi+ma)**2)*(k0-(mpi-ma)**2)/4./k0
         k0 = Sqrt(k02)
         epinuc = Sqrt(k02+mpi2)+Sqrt(k02+ma*ma)-ma+mn-5.-((k0*aovera/2.
     1   )**2)/twomu
c     Do loop over p
 20      Do 30 n=1,14
            rt(n) = 0.
            it(n) = 0.
 30      Continue
         Do 220 np=1,nps
            p = ptp(np)
            p2 = p*p
            If (nkap.gt.nfold) GoTo 40
c      sp is old folding technique
            sp = rhl+rhl1*Sqrt(mn2+p*p)
            If (nifty(5).ge.5) sp = epinuc-p2/twomu
 40         pa2 = (p*a)**2
            If (pa2.gt.150.) GoTo 200
            rho = anorm*exp(-pa2)
            GoTo (50,60,70,80), nna
c      1s shell
 50         rho = rho*na
c      1s + 1p
 60         rho = rho*(4.+(na-4.)*2.*pa2/3.)
            GoTo 90
c      1s +1p + (1d+2s)
 70         rho = rho*(4.+12.*(2./3.)*pa2+(na-16.)*(4.*pa2*pa2/15.+(3./2
     1      .)*(1.-2.*pa2/3.)**2)/2.)
            GoTo 90
c      1s +1p +1d +2s +1f
 80         rho = rho*(4.+12.*(2./3.)*pa2+12.*(4.*pa2*pa2/15.+(3./2.)*(1
     1      .-2.*pa2/3.)**2)+(na-40.)*(8./105.)*pa2**3)
 90         Continue
            If (nkap.gt.nfold) GoTo 150
c     Do loop over cos theta
            Do 100 n=1,nwaves
               rtcos(n) = 0.
               itcos(n) = 0.
 100        Continue
            If (nkap.eq.nfold) onecos = 0.
            Do 140 ncos=1,ncoss
               If (nifty(5).lt.5.) GoTo 110
c       3 body def of e, neg ke possible If use denom
               s = sp-2.*p*k0*ptcos(ncos)*aovera/twomu
               Call tnoff (kap,kap,s,k0,k0,ptcos(ncos))
               GoTo 120
 110           s = sp-2.*plab*p*ptcos(ncos)
               If (s.lt.rhl4) GoTo 140
c      s.lt.s(k. e. =0.)    get  out  of  Do  loop
               kapfld = Sqrt((s-rhl4)*(s-rhl3)/4./s)
               Call tncm (kapfld)
 120           Do 130 n=1,14
                  rtcos(n) = rtcos(n)+wtcos(ncos)*retij(n)
                  itcos(n) = itcos(n)+wtcos(ncos)*imtij(n)
 130           Continue
               If (nkap.eq.nfold) onecos = onecos+wtcos(ncos)
 140        Continue
            GoTo 180
c      fold 1/denom,stored in denom,only fold p*p/twomu(what else@)
c       set all of folded denom Real to avoid folded garbage
 150        ecm = ecmden(nkap-nfold)-p*p/twomu-5.
            Call lagrng (ecm,ecmden,g,denom,necms,14,2,25,14)
            Do 170 n=1,14
               n2 = 2*n-1
               If (g(n2).eq.0.) GoTo 160
               rtcos(n) = 2./g(n2)
 160           it(ncos) = 0.
 170        Continue
 180        rhl2 = rho*p*p*wtp(np)
            Do 190 n=1,14
               rt(n) = rt(n)+rhl2*rtcos(n)
               it(n) = it(n)+rhl2*itcos(n)
 190        Continue
            If (nkap.eq.nfold) one = one+rhl2*onecos
 200        If (np.ne.nps) GoTo 220
c     last point in Do, store folded t in dummny d/eta
            Do 210 n=1,14
               d(n,nkap) = rt(n)
               eta(n,nkap) = it(n)
 210        Continue
 220     Continue
 230  Continue
c     store up folded t*s, reset retij,imtij
      Rewind 3
      Write (3,310) nfold
      Write (3,320) (kapd(ne),(d(nwave,ne),eta(nwave,ne),nwave=1,14),ne=
     11,nfold)
      Rewind 3
      Write (6,340) one
      Do 250 ne=1,nfold
         Do 240 nwave=1,14
            rewave(nwave,ne) = d(nwave,ne)
            imwave(nwave,ne) = eta(nwave,ne)
 240     Continue
c     -----compute partial wave decomp of amplitude as check for print o
c          t here refers to the partial wave decomp
         rhl = 2./3.
         rt(1) = rewave(1,ne)
         rt(2) = imwave(1,ne)
         rt(3) = rt(1)/3.+rhl*rewave(4,ne)
         rt(4) = rt(2)/3.+rhl*imwave(4,ne)
         rt(5) = rewave(2,ne)+2.*rewave(3,ne)
         rt(6) = imwave(2,ne)+2.*imwave(3,ne)
         rt(7) = rt(5)/3.+rhl*(rewave(5,ne)+2.*rewave(6,ne))
         rt(8) = rt(6)/3.+rhl*(imwave(5,ne)+2.*imwave(6,ne))
         Write (6,300) kapd(ne),(rt(n),n=1,8)
 250  Continue
      If (nkape.ne.nfold) GoTo 260
      Return
c     put folded denom back on 4
 260  Rewind 4
      Write (4,330) necms
      Write (6,350) (ne,ecmden(ne),(denom(nwave,ne),nwave=1,6),ne=1,
     1necms)
      nint = nfold+1
      Do 280 ne=nint,nkape
         nn = ne-nfold
         Do 270 nwave=1,14
            n2 = 2*nwave
            n1 = n2-1
            If (d(nwave,ne).ne.0.) denom(n1,nn) = 1./d(nwave,ne)
            denom(n2,nn) = 0.
 270     Continue
         Write (4,360) ecmden(nn),(denom(nwave,nn),nwave=1,14)
 280  Continue
      Write (6,350) (ne,ecmden(ne),(denom(nwave,ne),nwave=1,6),ne=1,
     1necms)
      Rewind 4
      Return
 360  format(e10.3/(8e10.3))
c
 290  format (51h ***** in tfold,nstore lt 60,nfoldset= nstore ,****,i5/
     1)
 300  format (9x,f8.0,8x,8e12.4)
 310  format (i5)
 320  format (7e11.4/6e11.4/6e11.4/6e11.4/4e11.4)
 330  format (i3)
 340  format (53h1 two body t matrices now folded and stored on tape3;,
     1 22h folded denom on tape4,/,28h nb tape3 should be saved...,
     2 32h so Dont have to fold each time,/,5h one=,e15.6,
     3 38h kap,ret,imt(p,n,s,p waves) now follow)
 350  format (15h denom(nf,fold),i5,e12.4,10x,6e13.3)
      End
