      SUBROUTINE TFOLD
C *** SEE R H 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) GO TO 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) GO TO 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) GO TO 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) GO TO 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.) GO TO 200
            RHO = ANORM*EXP(-PA2)
            GO TO (50,60,70,80), NNA
C      1S SHELL
 50         RHO = RHO*NA
C      1S + 1P
 60         RHO = RHO*(4.+(NA-4.)*2.*PA2/3.)
            GO TO 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.)
            GO TO 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) GO TO 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.) GO TO 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))
               GO TO 120
 110           S = SP-2.*PLAB*P*PTCOS(NCOS)
               IF (S.LT.RHL4) GO TO 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
            GO TO 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.) GO TO 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) GO TO 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) GO TO 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
