      PROGRAM CHUCK
      IMPLICIT REAL*8(A-H,O-Z)
C     PROGRAM CHUCK(INPUT,OUTPUT,PUNCH,TAPE5=INPUT,TAPE6=OUTPUT
C    1,TAPE7=PUNCH,TAPE8=2001,TAPE9,TAPE4,TAPE2)
C***********************************************************************
C
C                  *****     CHUCK-3     *****
C
C     CHANGES FROM BASIC COLORADO EXPORT VERSION,  J.R.COMFORT
C
C   1. THE NUMBER OF CONTROL INTEGERS WAS INCREASED FROM 16 TO 20.
C   2. CONTROL INTEGER 1 WAS MODIFIED.  ANGLE DATA ARE CONTROLED AS--
C          C.I. 1 = 0  REPEAT SAME SET OF ANGLES AS LAST CASE.  DO NOT READ
C                        CARD 2.  (DEFAULT SET IS 0-180 DEG. IN 5 DEG. STEPS)
C                 = 1  READ NUMBER OF ANGLES, INITIAL ANGLE, STEP SIZE.
C                 = 2  READ INITIAL ANGLE, FINAL ANGLE, STEP SIZE
C                        (OR FINAL ANGLE, INITIAL ANGLE, STEP SIZE).
C   3. CONTROL INTEGER 2 WAS ADDED TO SUPPRESS FORM-FACTOR PRINTING.
C          C.I. 2 = 0  PRINT ALL FORM-FACTOR INFORMATION.
C                 = 1  SAME.
C                 = 2  SUPPRESS SINGLE-PARTICLE WAVE FUNCTIONS, ICODE=2.
C                 = 3  SUPPRESS PARAMETER LISTINGS ALSO (ANY ICODE).
C   4. CONTROL INTEGER 3 IS NOT USED ANYMORE, POLARIZATIONS ARE ALWAYS
C        PRINTED.
C   5. CONTROL INTEGER 4 IS NOT USED ANYMORE INSTEAD CONTROL INTEGER
C        13 IS USED TO CONTROL PUNCHING OF CROSS SECTIONS AND POLARIZAT-
C        IONS.
C   6. CONTROL INTEGER 9 GIVES 3-CYCLE SEMILOG PLOTS BY DEFAULT.. PLOTS
C        SUPPRESSED IF C.I. 9 IS GREATER THAN 5.
C   7. THE PHYSICAL CONSTANTS IN THIS VERSION ARE BASED ON TAYLOR, PARKE
C        LANGENBERG (RMP 41,375,1969), EXCEPT FOR AMU WHICH IS BASED ON
C        1973 TABLE OF FUNDAMENTAL CONSTANTS (E.R.COHEN).
C   8. THE SCRATCH TAPES ARE TAPE 8 AND TAPE 9.
C   9. THE DATE AND TIME ARE CURRENT DATE AND CPU TIME OBTAINED BY CALLS TO
C        SUBROUTINE CUDATE.  THE DATA WORDS ARE DOUBLE PRECISION AND THE
C        FORMAT IS (A10,2X,F7.2).  THE FORMATS ARE 9003 IN MAIN AND 9019
C        ACHK, BCHK, SOLVE, AND XSECT.
C  10. PRINTING OF THE FIRST TITLE WAS MOVED FROM THE MAIN ROUTINE TO
C        SUBROUTINE ACHK, AFTER THE FIRST READ STATEMENT.
C  11. THE WATSON L.S POTENTIAL WAS ADDED TO POTENT.
C  12. THE FINITE-RANGE ROUTINE FNLOC WAS MODIFIED TO BE THE SAME AS DWUCK.
C  13. ALL NON-LOCALITY CORRECTIONS USE THE OLD L.E.A. METHOD (HULTHEN FORM
C        IN DWUCK-2) AS STANDARD.  THE ROST-KUNZ (EXPONENTIAL) FORM IS O???
C        BY SETTING PNLOC LESS THAN ZERO FOR EACH PROJECTILE FOR WHICH IT IS
C        TO BE APPLIED.
C  14. GAUSSIAN MICROSCOPIC INTERACTION INCLUDED IN CATHEN, WITH ICODE=4.
C  15. SUBROUTINE DWPLOT IS A MODIFIED VERSION.
C  16. THE ANGULAR MOMENTUM PACKAGE HAS BEEN MODIFIED SO AS TO WORK WITH
C        LOG FACTORIALS.  THE INTERNAL TABLES ARE EXPANDED UP TO FACTORIAL 82
C        THE HANDLING OF THE YXFCT FUNCTION CALLS IN SUBROUTINE SOLVE WAS
C        MODIFIED.
C  17. EXPLICIT DOUBLE PRECISIONING HAD TO BE DONE IN THE INTEGRATION SETUP
C        ON THE DEC-10.  THIS INVOLVED SEVERAL VARIABLES IN SUBROUTINES ???
C        AND SOLVE, AND A REWRITTING OF SUBROUTINE INVERT.  NOTE--
C           A. THE WRITING OF REAL WORDS INTO AND FROM THE BUFFER IN SUBROUTINE
C              WRITEB AND READB WAS DOUBLED IN LENGTH.
C  18. THE LIMIT ON THE NUMBER OF INTEGRATION POINTS WAS INCREASED TO 40.
C  19. THE LIMIT ON THE NUMBER OF TRANSITION AMPLITUDES WAS INCREASED FROM ???
C        TO 12800 AND THE LIMIT ON STOP CODE 111 WAS INCREASED FROM 192 TO ???
C        THE NUMBER OF SUBCHANNELS WAS INCREASED FROM 32 TO 40.  THE MAXIMUM
C        NUMBER OF PARTIAL WAVES FOR ALL CHANNELS HAS GONE FROM 400 TO 6???
C  20. AN ARRAY JPRT WAS ADDED TO CONTAIN THE PARITY OF EACH NUCLEAR STATE.
C        ODD VALUE (E.G., 1 OR -1) WILL CHANGE THE SIGN ON THE LTRT ARRAY
C        INTENDED FOR USE WITH 0- STATES.
C  21. CORRECTIONS MADE FOR THE PROPER HANDLING OF INTEGRATION COORDINATES
C        IN MICROSCOPIC INELASTIC SCATTERING.  USE CORE MASS = A-1.
C  22. MODIFICATIONS WERE MADE IN THE HANDLING OF UNBOUND FORM FACTORS,
C        FOLLOWING VINCENT AND FORTUNE.  THE HUBY-MINES DAMPING WAS CHANGED
C        TO GAUSSIAN DAMPING.
C  23. NEGATIVE RADIUS VALUES IMPLY A1**(1/3)+A2**(1/3).
C  24. LOWER-L CUTOFF FACILITY ADDED.
C  25. A FLOATING-CORE FEATURE WAS ADDED.  SEE SUBROUTINE BCHK.
C  26. MICROSCOPIC INELASTIC SCATTERING MODIFIED FOR BETTER TREATMENT.
C  27. FOLDING OPTION HAS BEEN ADDED. THE TRANSITION DENSITY IS ENTERED
C      IN A SEPARATE SUBROUTINE FMN WHICH HAS TO BE COMPILED AND LINKED
C      TO THE PROGRAM EVERY TIME A FOLDING CALCULATION IS MADE.
C  28. FORM FACTORS FOR BETA, GAMMA AND OCTUPOLE BANDS.
C
C***********************************************************************
C     DOUBLE PRECISION BETA
      COMMON THETAN,THETA1,DTHETA,DRT,RMAX,ALPHA(15),IDAT(6)
     1,ICON(20),LMAX,KT,KMAX,LLCO,JOUT,NLTR,IBUFF,IWORD
     2,KLAM( 8),LTRT( 8),JTRT( 8),NADD( 8),LAM(16),ISP(16),JSP(16)
     3,ILK(16),JPRT( 8)
      COMMON EN( 8),FM( 8),Z( 8),FMA( 8),ZA( 8),RY( 8),FS( 8),ECM( 8)
     1,FMU( 8),RC( 8),FK2( 8),FK( 8),ETA( 8),AFACT( 8),DRF( 8),FN( 8)
     2,PNLOC( 8),DUM( 8),SIG( 8)
     3,IS( 8),NS( 8),NLI( 8),LKTAB(32),ICUP(32),JCUP(32)
     4,IN( 8),INN( 8, 8),INDTAB(10)
      COMMON/LOCAL /ARRAY(4400)
      COMMON/BUFSPC/ BUFFR(800)
C
C     CALL ERRSET(10)
      KMAX=400
      IBUFF=800
      THETAN=37.D0
      THETA1=0.D0
      DTHETA=5.D0
 5001 CONTINUE
C
C     CUDATE IS DECODER SUBROUTINE FOR DATE AND TIME ROUTINES
C
      CALL CUDATE(IDAT)
C
C
C     READ IN CARD SET 1
C
      READ(5,9000,END=1)ICON,ALPHA
      IF (ICON(1).GE.8) STOP
      WRITE(6,9001)
      WRITE(6,9005)
      WRITE(6,9002)
      WRITE(6,9003) ICON,ALPHA,IDAT
C      WRITE(6,9003) ICON,ALPHA,ICON,XTIME
      REWIND 8
      REWIND 9
      CALL ACHK
      IF(JOUT.GT.1) STOP
      IF(JOUT.NE.0) GO TO 5001
      REWIND 8
      REWIND 9
      CALL BCHK
      GO TO 5001
 1    STOP
 9000 FORMAT(20I1,15A4)
 9001 FORMAT(  '1 COUPLED CHANNELS CODE--CHUCK3-VERSION 01/JULY/1978.  J
     1.R. COMFORT EXTENDED VERSION 25/JUNE/1979.'/
     2,55X,'M.N. HARAKEH EXTENDED VERSION 1/JANUARY/1982.'/
     3,55X,'M. ELVERS AND J. ENDRES MODIFIED VERSION 10/NOVEMBER/2009.')
 9002 FORMAT(41H  1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0,22H    RUN IDE
     1NTIFICATION)
 9003 FORMAT(1H ,20I2,4X,15A4,4X,I2,2(1H/,I2),I4,2(1H.,I2))
C 9003 FORMAT(1H ,20I2,4X,15A4,4X,A10)
 9005 FORMAT(18H0CONTROL INTEGERS  )
      END
      SUBROUTINE ACHK
      IMPLICIT REAL*8(A-H,O-Z)
C     ABS(XXX)=DABS(XXX)
C     FLOAT(III)=DFLOAT(III)
C     SIGN(XXX,YYY)=DSIGN(XXX,YYY)
C     DOUBLE PRECISION BETA
C     Modified by (JE and ME)
C     USE THE NEXT 1 CARDS FOR 32-BIT MACHINES
      DOUBLE PRECISION VCC
      COMMON THETAN,THETA1,DTHETA,DRT,RMAX,ALPHA(15),IDAT(6)
     1,ICON(20),LMAX,KT,KMAX,LLCO,JOUT,NLTR,IBUFF,IWORD
     2,KLAM( 8),LTRT( 8),JTRT( 8),NADD( 8),LAM(16),ISP(16),JSP(16)
     3,ILK(16),JPRT( 8)
      COMMON EN( 8),FM( 8),Z( 8),FMA( 8),ZA( 8),RY( 8),FS( 8),ECM( 8)
     1,FMU( 8),RC( 8),FK2( 8),FK( 8),ETA( 8),AFACT( 8),DRF( 8),FN( 8)
     2,PNLOC( 8),FMFCT( 8),SIG( 8)
     3,IS( 8),NS( 8),NLI( 8),LKTAB(32),ICUP(32),JCUP(32)
     4,IN( 8),INN( 8, 8),INDTAB(10)
      COMMON/LOCAL/ W(800),V(400),SAV1(800),SAV2(800),SAV3(1600)
      COMMON/BUFSPC/ X(800,1)
      COMMON /KPROJ/ MU
      COMMON /CORR/ BETAR,BETAI,VBETA,RHO0R,RHO0I,INDDEF,INDIM
      DIMENSION GAMMA(2),AKIN(32),BKIN(8,32),IJ(2),LDFRM(4),ABETA(3)
      EQUIVALENCE (BKIN(1,1),EN(1))
      DATA LCORE/600/
      DATA GAMMA/3H/2 ,3H   /
      DATA NLKMAX/16/
      DATA HBARC,AMU/197.329,931.502/
C modified by (JE, ME)
      SAVE
C
      KMAX2=KMAX+KMAX
      IF(ICON(1).EQ.0) GO TO 6
C
C     READ IN CARD SET 2
C
      READ(5,9008)THETAN,THETA1,DTHETA
      IF(THETA1.LE.0.0) THETA1=0.D0
      IF(DTHETA.LE.0.0) DTHETA=5.D0
      IF(ICON(1).EQ.1) GO TO 5
      IF(THETAN.GE.THETA1) GO TO 4
      ST=THETAN
      THETAN=THETA1
      THETA1=ST
    4 CONTINUE
      IF(THETAN.LE.0.0) THETAN=180.D0
      NTHETA=NINT((THETAN-THETA1)/DTHETA)+1
      THETAN=NTHETA
      GO TO 6
    5 CONTINUE
      IF(THETAN.LE.0.0) THETAN=37.D0
    6 CONTINUE
      WRITE(6,9010)THETAN,THETA1,DTHETA
C
C     READ IN CARD SET 3
C
      READ (5,9004)LMAX,NLTR,(LTRT(I),I=1,NLTR),(KLAM(I),I=1,NLTR)
     1    ,(JPRT(I),I=1,NLTR)
      WRITE(6,9026)LMAX,(LTRT(I),I=1,NLTR)
      WRITE(6,9027)NLTR,(KLAM(I),I=1,NLTR)
      ISUM=0
      DO 7 I=1,NLTR
      ISUM=ISUM+JPRT(I)
    7 CONTINUE
      IF(ISUM.NE.0) WRITE(6,9029) (JPRT(I),I=1,NLTR)
      LMAX=MIN0(LMAX+1,LCORE/NLTR)
C
C     READ IN CARD SET 4
C
      READ (5,9008)DRT,RMAX,FLCO
      IF(DRT.LE.0.0) DRT=0.1D0
      IF(RMAX.EQ.0.0) RMAX=20.D0
      LLCO=FLCO
      IF(LLCO.LT.0) LLCO=0
      WRITE(6,9028)DRT,RMAX,LLCO
      KT=NINT(ABS(RMAX)/DRT)
      KTMIN=KT
      CALL CUDATE(IDAT)
      WRITE(6,9019)ALPHA,IDAT
C      WRITE(6,9019)ALPHA,IDAT,XTIME
C
C     READ IN DIAGONAL CHANNEL CARD SETS
C
      DO 200 I=1,NLTR
      K=I
      JTRT(I)=IABS(LTRT(I))
      LTEMP=LTRT(I)
      TEMP=GAMMA(1)
      IF(PHASEF(LTRT(I)))86,85,85
   85 CONTINUE
      LTEMP=LTEMP/2
      TEMP=GAMMA(2)
   86 CONTINUE
      WRITE(6,9011)I,LTEMP,TEMP
C
C
      CALL KINEX(AKIN,W(1  ),ECM(1),DRT,K,KMAX,ICON)
      DO 140 M=1,19
      BKIN(I,M)=AKIN(M)
  140 CONTINUE
      IS(I)=FS(I)
      FS(I)=FS(I)/2.0
      NS(I)=IS(I)+1
C
C
      READ (5,9018)II,JJ,LAMBDA,ISPIN,JSPIN,ICODE,LDFRM(1),ICOUEX,ABETA
      LAMBDA=0
      DO 148 J=1,3
      LDFRM(J+1)=LDFRM(J)+2
      IF(ABETA(J).NE.0.0) WRITE(6,9007)LAMBDA,ABETA(J),LDFRM(J)
  148 CONTINUE
      DO 150 M=1,KMAX
      V(M  )=0.D0
  150 CONTINUE
      CALL POTENT(AKIN,W(1  ),V(1  ),KT,KMAX,ICODE,LAMBDA,LDFRM,ABETA
     1 ,JOUT,ICON,K)
      WRITE(9) W,V
      IF(ICON(11).EQ.0) GO TO 190
C
C     WRITE OUT K(I,I)**2 ELEMENTS
C
      WRITE(6,9901)I
      KT2=KT+KT
      IND=0
      R=DRF(I)
      DO 180 M=1,KT2,10
      INC=IND+1
      IND=MIN0(INC+9,KT2)
      WRITE(6,9900)R,(W(J  ),J=INC,IND)
      WRITE(7,9900)R,(W(J  ),J=INC,IND)
      R=R+5.0*DRF(I)
  180 CONTINUE
  190 CONTINUE
  200 CONTINUE
C
C
      JOUT=0
      KTFLG=0
C
C      READ IN CHANNEL COUPLING
C
      CALL CUDATE(IDAT)
      WRITE(6,9019)ALPHA,IDAT
      DO 210 I=1,8
      DO 210 J=1,8
      INN(I,J)=0
  210 CONTINUE
      LK=0
      NCUP=0
  220 CONTINUE
      READ (5,9018)II,JJ,LAMBDA,ISPIN,JSPIN,ICODE,LDFRM(1),ICOUEX,ABETA
     1,RCX,FNRNG,AGAMMA,VBETA
      ICOD=ICODE
      ICODE=IABS(ICODE)
      IF(ICOD.EQ.7.OR.(ICODE.LT.7.OR.ICODE.GT.8)) GO TO 221
      IF(ABS(VBETA).LT.1.0E-6) VBETA=1.D0
      BETAR= FNRNG/VBETA
      BETAI=AGAMMA/VBETA
  221 CONTINUE
      IF(II.EQ.0) GO TO 985
      IF(JJ.EQ.0) JJ=1
      I=IABS(II)
      J=IABS(JJ)
      MU=IABS(KLAM(I)-KLAM(J))
      WRITE(6,9017)II,JJ,ICODE
      LK=LK+1
      IF(LK.LE.NLKMAX) GO TO 232
  231 CONTINUE
      WRITE(6,9054)
      JOUT=2
      GO TO 2000
  232 CONTINUE
      NCUP=NCUP+1
      ICUP(NCUP)=J
      JCUP(NCUP)=I
      LKTAB(NCUP)=LK
      INN(I,J)=1
      IF(I.EQ.J) GO TO 260
      IF(II.LT.0) GO TO 260
      NCUP=NCUP+1
      ICUP(NCUP)=I
      JCUP(NCUP)=J
      LKTAB(NCUP)=LK
      INN(J,I)=1
  260 CONTINUE
      LAM(LK)=LAMBDA
      ISP(LK)=ISPIN
      IF(ISPIN.EQ.0) JSPIN=LAMBDA+LAMBDA
      JSP(LK)=JSPIN
C   THE FOLLOWING STATEMENT IS TO WRITE A ZERO FOR FNRNG SINCE FNRNG IS
C   GIVEN A VALUE FOR THE PURPOSE OF USING CORRECTION TERMS FOR BETA
C   AND OCTUPOLE VIBRATIONS
      IF(ICOD.EQ.-7.OR.ICOD.EQ.8)
     1WRITE(6,9037) LAM(LK),ISP(LK),JSP(LK),ICOUEX,RCX,ABETA(3)
      IF(ICOD.EQ.7.OR.(ICODE.LT.7.OR.ICODE.GT.8))
     1WRITE(6,9037) LAM(LK),ISP(LK),JSP(LK),ICOUEX,RCX,FNRNG
      IF(ICOD.EQ.7) WRITE(6,9038) AGAMMA
      IF(ICOD.EQ.-7) WRITE(6,9039) VBETA
      IF(ICOD.EQ.-7) WRITE(6,9041) FNRNG,AGAMMA
      IF(ICOD.EQ.8) WRITE(6,9040) VBETA
      IF(ICOD.EQ.8) WRITE(6,9042) FNRNG,AGAMMA
      DO 265 K=1,3
      LDFRM(K+1)=LDFRM(K)+2
      IF(ABETA(K).NE.0.0) WRITE(6,9007)LAM(LK),ABETA(K),LDFRM(K)
  265 CONTINUE
      IF(ABETA(1).NE.0.0) GO TO 270
      JOUT=1
      WRITE(6,9050)
  270 CONTINUE
C
C     CHECK COUPLING FOR CONSISTENCY
C
      IF(ISPIN.NE.0) GO TO 282
      IF((ISPIN.EQ.0).AND.(IS(I).EQ.IS(J))) GO TO 282
      JOUT = 1
      WRITE(6,9051)IS(I),IS(J),ISPIN
 282  CONTINUE
      IF(PHASEF(IS(I)+IS(J)+ISPIN).GT.0.0) GO TO 283
      JOUT=1
      WRITE(6,9051)IS(I),IS(J),ISPIN
  283 CONTINUE
      IF(PHASEF(LTRT(I)+LTRT(J)+JSPIN).GT.0.0) GO TO 284
      JOUT=1
      WRITE(6,9052)LTRT(I),LTRT(J),JSPIN
  284 CONTINUE
      IF(DFLOAT(ISIGN(1,LTRT(I))*ISIGN(1,LTRT(J)))
     X  *PHASEF(JPRT(I))*PHASEF(JPRT(J))
     1*PHASEF(LAMBDA).GT.0.0) GO TO 285
      JOUT=1
      WRITE(6,9053)LTRT(I),LTRT(J),LAMBDA
  285 CONTINUE
      IF(I.LE.NLTR) GO TO 286
      JOUT=1
      WRITE(6,9055)
  286 CONTINUE
      IF(J.LE.NLTR) GO TO 287
      JOUT=1
      WRITE(6,9055)
  287 CONTINUE
C
C
      BFACT=-ABETA(1)*2.0*AMU/HBARC**2
      ILK(LK)=ICODE
      IF(ICODE.EQ.0) GO TO 300
      IF(ICODE.EQ.1) GO TO 320
      IF(ICODE.EQ.2) GO TO 340
      IF(ICODE.EQ.3) GO TO 300
      IF(ICODE.EQ.4) GO TO 360
      IF(ICODE.EQ.5) GO TO 380
      IF(ICODE.EQ.6) GO TO 400
C
C    ICODE.EQ.7 FOR EXCITATION OF GAMMA BAND.OR. EXCITATION OF BETA BAND
C    ICODE.EQ.8 FOR EXCITATION OF OCTUPOLE BANDS
C
      IF(ICODE.EQ.7) GO TO 300
      IF(ICODE.EQ.8) GO TO 300
      WRITE(6,9056)ICODE
      JOUT=2
      GO TO 2000
C
C     COLLECTIVE FORM FACTOR
C
  300 CONTINUE
      IF(ICODE.NE.0) GO TO 305
      IF((KLAM(I).EQ.KLAM(J)).AND.PHASEF(JTRT(J)+KLAM(J)).GT.0.0)
     1  GO TO 305
      JOUT=1
      WRITE(6,9050)
  305 CONTINUE
      RCX=RCX*AFACT(J)
      IF(RCX.LT.0.0) RCX=RCX+RCX*FMFCT(J)
      RCX=ABS(RCX)
      AKIN(18)=FMFCT(J)
      AKIN(19)=1.D0
      AKIN(14)=AFACT(J)
      AKIN(15)=DRF(J)
      DO 310 M=1,KMAX
      MK=M+M-1
      X(MK  , 1)=0.D0
      X(MK+1, 1)=0.D0
  310 CONTINUE
C
      CALL POTENT(AKIN,X,SAV2,KT,KMAX,ICOD,LAMBDA
     1,LDFRM,ABETA,JOUT,ICON,0)
      TEMP = 1.D0
      IF(ICOD.EQ.7) TEMP=AGAMMA*ABETA(1)
      IF(ICOD.EQ.-7.OR.ICODE.EQ.8) TEMP=VBETA
      R = 0.D0
      SUMR1=0.D0
      SUMI1=0.D0
      SUMR2=0.D0
      SUMI2=0.D0
      DO 315 M=1,KT
      MK=M+M-1
      X(MK  ,1)=BFACT*X(MK  ,1)*TEMP
      X(MK+1,1)=BFACT*X(MK+1,1)*TEMP
      R=R+DRT
      R2=R*R
      RL=R2*R**LAMBDA
      IF(LAMBDA.EQ.0.OR.LAMBDA.EQ.1) RL=RL*R2
      SUMR2=SUMR2+X(MK  ,1)*RL
      SUMI2=SUMI2+X(MK+1,1)*RL
      IF(LAMBDA.GT.1) GO TO 315
      SUMR1=SUMR1+X(MK  ,1)*RL/R2
      SUMI1=SUMI1+X(MK+1,1)*RL/R2
  315 CONTINUE
      TEMP=ZA(I)*ABETA(1)/BFACT
      IX=JTRT(I)
      IY=JTRT(J)
      KLMM=KLAM(I)
      XIX=SQRT(DFLOAT(IX+1))
      LL2=2*LAMBDA
      IF(ICODE.EQ.0) TEMP=TEMP*XIX*VCC(IX,LL2,IY,KLMM,0)
      IF(ICODE.EQ.3) TEMP=TEMP*XIX/SQRT(DFLOAT(LL2+1))*PHASEF((IX-IY)/2)
      IF(ICODE.EQ.7.OR.ICODE.EQ.8) TEMP=TEMP*XIX*VCC(IX,LL2,IY,-KLMM,MU)
      SUMR2=TEMP*SUMR2/(12.56637*RHO0R)
      SUMI2=TEMP*SUMI2/(12.56637*RHO0I)
      SUMR1=TEMP*SUMR1/(12.56637*RHO0R)
      SUMI1=TEMP*SUMI1/(12.56637*RHO0I)
      IF(LAMBDA.NE.0) GO TO 313
C   THE MONOPOLE TRANSITION IS CALCULATED FOR THE OPERATOR (1/2)*R**2
      SUMR2=SUMR2*1.772454
      SUMI2=SUMI2*1.772454
  313 CONTINUE
      IF(LAMBDA.NE.1) GO TO 314
C   THE DIPOLE TRANSITION IS CALCULATED FOR THE OPERATOR (1/2)*R**3*Y1
      SUMR2=SUMR2*0.5
      SUMI2=SUMI2*0.5
  314 CONTINUE
      IF(LAMBDA.EQ.0) WRITE(6,9043) SUMR1,SUMI1,SUMR2,SUMI2
      IF(LAMBDA.EQ.1) WRITE(6,9044) SUMR1,SUMI1,SUMR2,SUMI2
      IF(LAMBDA.GT.1) WRITE(6,9045) LAMBDA,SUMR2,LAMBDA,SUMI2
      IF(INDDEF.EQ.1.AND.INDIM.EQ.1) THEN
        WRITE(6,9046)
        WRITE(6,9048)
        WRITE(6,9047)
      ELSEIF(INDDEF.EQ.1) THEN
        WRITE(6,9046)
        WRITE(6,9047)
      ELSEIF(INDIM.EQ.1) THEN
        WRITE(6,9048)
        WRITE(6,9047)
      ELSE
      ENDIF
      IF(ICOUEX.EQ.0) GO TO 318
C
C     CALCULATE COULOMB EXCITATION
C
      ETAK= 6.0*ETA(J)*FK(J)*RCX**LAMBDA
     1/(FMU(J)*DFLOAT(LAMBDA+LAMBDA+1))
      ETAK=-ABETA(1)*ETAK
      R=0.D0
      DO 316 M=1,KT
      R=R+DRF(I)
      IF(R.LT.RCX) GO TO 316
      MK=M+M-1
      X(MK  ,1)=X(MK  ,1)+ETAK/R**(LAMBDA+1)
  316 CONTINUE
  318 CONTINUE
      GO TO 960
C
C     SINGLE PARTICLE ORBITAL
C
  320 CONTINUE
      CALL KINEX(AKIN,SAV1,ECMT,DRT,0,KMAX,ICON)
      DO 330 M=1,KMAX
      MK=M+M-1
      SAV2(MK  )=0.D0
      SAV2(MK+1)=0.D0
  330 CONTINUE
      CALL POTENT(AKIN,SAV1,SAV2,KT,KMAX,1,0,0,0.D0,JOUT,ICON,0)
      DO 335 M=1,KMAX
      MK=M+M-1
      X(MK  ,1)=SAV1(MK)*BFACT
      X(MK+1,1)=0.D0
  335 CONTINUE
      GO TO 960
C
C     MICROSCOPIC MODEL FORM FACTOR
C
  340 CONTINUE
      DO 350 M=1,KMAX
      MK=M+M-1
      X(MK  , 1)=0.D0
      X(MK+1, 1)=0.D0
  350 CONTINUE
      DR=DRT
      IERR=0
      CALL CATHEN(LAMBDA,ISPIN,JSPIN,KT,KMAX,SAV1,SAV3,X,DR,ICON,IERR)
      IF(IERR.EQ.0) GO TO 352
      JOUT=IERR
      WRITE(6,9050)
  352 CONTINUE
      IF(KT.LT.0) KTFLG=MAX0(-KT,KTFLG)
      KT=IABS(KT)
      AKIN(15)=DR
      IF(ABS(FMA(J)-FMA(I)).GT.1.2) ILK(LK)=1
      DO 355 M=1,KMAX
      MK=M+M-1
      X(MK  , 1)=X(MK  , 1)*BFACT
      X(MK+1, 1)=X(MK+1, 1)*BFACT
  355 CONTINUE
      GO TO 960
C
C     'CLUSTER' (COLLECTIVE) COMPLEX F.F., HANDLED AS STRIPPING
C
  360 CONTINUE
      DO 365 M=1,KMAX
      MK=M+M-1
      X(MK  ,1)=0.D0
      X(MK+1,1)=0.D0
  365 CONTINUE
      CALL KINEX(AKIN,SAV1,ECMT,DRT,0,KMAX,ICON)
      AKIN(19)=1.D0
      CALL POTENT(AKIN,X,SAV2,KT,KMAX,ICODE,LAMBDA
     1,LDFRM,ABETA,JOUT,ICON,0)
      DO 368 M=1,KMAX
      MK=M+M-1
      X(MK  ,1)=BFACT*X(MK  ,1)
      X(MK+1,1)=BFACT*X(MK+1,1)
  368 CONTINUE
      GO TO 960
C
C     RESCALE LAST F.F. WITH NEW STRENGTH
C
  380 CONTINUE
      DO 390 M=1,KMAX
      MK=M+M-1
      X(MK  ,1)=ABETA(1)*X(MK  ,1)
      X(MK+1,1)=ABETA(1)*X(MK+1,1)
  390 CONTINUE
      ILK(LK)=ILK(LK-1)
      GO TO 960
C
C     STRONG COUPLING COLLECTIVE MODEL
C     SPECTROSCOPIC FACTOR ALGEBRA
C
  400 CONTINUE
      JT=JSPIN
      IF((PHASEF(JTRT(I)+KLAM(I)).GT.0.0).AND.
     1  (PHASEF(JTRT(J)+KLAM(J)).GT.0.0)) GO TO 402
      JOUT=1
      WRITE(6,9050)
  402 CONTINUE
      IF(FMA(J).GT.FMA(I)+0.2) GO TO 405
C     STRIPPING
      MK=KLAM(I)-KLAM(J)
      BFACT=BFACT*SQRT(DFLOAT(JTRT(J)+1)/DFLOAT(JTRT(I)+1))
     1*VCC(JTRT(J),JT,JTRT(I),KLAM(J),MK)
      GO TO 410
C     PICKUP
  405 CONTINUE
      MK=KLAM(J)-KLAM(I)
      BFACT=BFACT*SQRT(DFLOAT(JTRT(I)+1)/DFLOAT(JTRT(J)+1))
     1*VCC(JTRT(I),JT,JTRT(J),KLAM(I),MK)
  410 CONTINUE
      IF(MK.LT.0) BFACT=BFACT*PHASEF((JT-MK)/2)
      ILK(LK)=1
      GO TO 320
C
  960 CONTINUE
      IF(ICODE.EQ.7.OR.ICODE.EQ.8) GO TO 980
      IF((ABS(FNRNG)+ABS(PNLOC(I))+ABS(PNLOC(J))).EQ.0.0) GO TO 980
      IJ(1)=I
      IJ(2)=J
      REWIND 9
      I4=MIN0(I,J)-1
      IF(I4.EQ.0) GO TO 964
      DO 962 M=1,I4
      READ(9)
  962 CONTINUE
  964 CONTINUE
      READ(9) (SAV3(M),M=1,800)
      I4=IABS(J-I)-1
      IF(I4.LT.0) BACKSPACE 9
      IF(I4.LE.0) GO TO 968
      DO 966 M=1,I4
      READ(9)
  966 CONTINUE
  968 CONTINUE
      READ(9) (SAV3(M),M=801,1600)
      CALL FNLOC(SAV3,X,SAV2,LAMBDA,BKIN,FNRNG,KT,IJ,AKIN,ICODE,KMAX)
  980 CONTINUE
C
C     WRITE COUPLING ELEMENTS ON TAPE
C
      WRITE(8)(X(M,1),M=1,KMAX2)
C
C     WRITE OUT K(I,J)**2 ELEMENTS
C
      IF(ICON(12).EQ.0) GO TO  984
      WRITE(6,9902)I,J
      KT2=KT+KT
      DR=AKIN(15)
      IND=0
      R=DR
      DO 983 M=1,KT2,10
      INC=IND+1
      IND=MIN0(INC+9,KT2)
      WRITE(6,9900)R,(X(K,1),K=INC,IND)
      WRITE(7,9900)R,(X(K,1),K=INC,IND)
      R=R+5.0*DR
  983 CONTINUE
  984 CONTINUE
      GO TO 220
  985 CONTINUE
      IF(KTFLG.NE.0) KT=MIN0(KT,KTFLG-5)
      NLI(1)=LK
      NLI(2)=NCUP
      IF(RMAX.LT.0.0) KT=KTMIN
      L=LMAX-1
      WRITE(6,9021)L,KT
 2000 CONTINUE
      RETURN
C
 9004 FORMAT(26I3)
 9007 FORMAT(18H                  ,9H   LAM  =,I4,14H        BETA =,F9.3
     1,9H   LDFRM=,I4)
 9008 FORMAT(10F8.4)
 9009 FORMAT(18H0INPUT DATA       ,9H   DR   =,F9.4,9H   RMAX =,F9.4)
 9010 FORMAT(18H0ANGLE DATA       ,9H   THETN=,F9.4,9H   THET1=,F9.4,9H
     1  DTHET=,F9.4)
 9011 FORMAT( 8H0CHANNEL,I2,8X,9H   ISPIN=,I3,A3,2X,80(1H*))
 9017 FORMAT( 8H0CHANNEL,I2,10H - CHANNEL,I2,9H COUPLING ,5X,67(1H*)
     1,9H   ICODE=,I4 )
 9018 FORMAT(8I3,7F8.4)
 9019 FORMAT(1H1,17X,15A4,4X,I2,2(1H/,I2),I4,2(1H.,I2))
 9021 FORMAT(1H0,17X,9H   LMAX =,I4,5X,9H   NSTEP=,I4)
 9026 FORMAT(18H0ANG MOM DATA     ,9H   LMAX =,I4,14H        2*JTR=,8I4)
 9027 FORMAT(18H                  ,9H   NJTR =,I4,14H        2*KLM=,8I4)
 9028 FORMAT(18H0CARD SET 4 DATA   ,9H   DRF  =,F9.4,9H   RMAX =,F9.4,9H
     1   LLCO =,I4)
 9029 FORMAT(40X,5HJPRT=,8I4)
 9037 FORMAT(18H                  ,9H   LAM  =,I4,14H        2*STR=,I4,1
     14H        2*JTR=,I4,14H        ICOUX=,I4,9H   RCOUL=,F9.4,9H   FNR
     2NG=,F9.4)
 9038 FORMAT(1H0,20X,'FOR ICODE=7 GAMMA =',F7.4,/)
 9039 FORMAT(1H0,20X,'FOR ICODE=-7  BETA-BETA2 =',F7.4,/)
 9040 FORMAT(1H0,20X,'FOR ICODE=8 BETA3 =',F7.4,/)
 9041 FORMAT(1H0,20X,'CORRECTION TERM TO CONSERVE NUMBER OF PARTICLES FO
     1R A MONOPOLE BETA VIBRATION',
     2/21X,'BETA0R =',F7.4,'   BETA0I =',F7.4,/)
 9042 FORMAT(1H0,20X,'CORRECTION TERM TO REMOVE SPURIOUS CENTER OF MASS
     1MOTION FOR ISOSCALAR DIPOLE TRANSITIONS IN OCTUPOLE BANDS',
     2/21X,'BETA1R =',F7.4,'   BETA1I =',F7.4,/)
 9043 FORMAT(10X,'VOLUME REAL MONOPOLE MATRIX ELEMENT =',E15.6,
     110X,'VOLUME IMAGINARY MONOPOLE MATRIX ELEMENT =',E15.6,
     2/10X,'REAL MONOPOLE MATRIX ELEMENT        =',E15.6,
     310X,'IMAGINARY MONOPOLE MATRIX ELEMENT        =',E15.6)
 9044 FORMAT(10X,'CENTER OF MASS REAL MATRIX ELEMENT =',E15.6,
     110X,'CENTER OF MASS IMAGINARY MATRIX ELEMENT =',E15.6,
     2/10X,'REAL DIPOLE MATRIX ELEMENT         =',E15.6,
     310X,'IMAGINARY DIPOLE MATRIX ELEMENT         =',E15.6)
 9045 FORMAT(10X,I1,'TH POLE REAL MATRIX ELEMENT =',E15.6,
     110X,I1,'TH POLE IMAGINARY MATRIX ELEMENT =',E15.6)
 9046 FORMAT(1X,'WARNING: BECAUSE OF IMPROPER NORMALIZATION, THE MATRIX
     1 ELEMENT FOR THIS DEFORMED FORM FACTOR CAN BE QUITE OFF.')
 9047 FORMAT(1X,'         USE  B E L G E N  TO CALCULATE THE CORRECT
     1 MATRIX ELEMENTS.')
 9048 FORMAT(1X,'WARNING: IMAGINARY MATRIX ELEMENTS ARE NOT CORRECT. TH
     1EY SHOULD BE CALCULATED SEPARATELY FOR EACH IMAGINARY POTENTIAL')
 9050 FORMAT(26H0NO TRANSITION STRENGTH   ,50(1H$))
 9051 FORMAT(39H0SPINS DO NOT SATISFY PROPER STATISTICS,10H  2*IS(I)=,I2
     1,10H  2*IS(J)=,I2,8H  2*STR=,I2,3X,50(1H$))
 9052 FORMAT(39H0SPINS DO NOT SATISFY PROPER STATISTICS,10H  2*I (I)=,I2
     1,10H  2*I (J)=,I2,8H  2*JTR=,I2,3X,50(1H$))
 9053 FORMAT(22H0PARITY DOES NOT CHECK ,10H   2*J(I)=,I2,10H   2*J(J)=,
     1 I2,10H   LAMBDA=,I2,3X,60(1H$))
 9054 FORMAT(30H0TOO MANY CHANNEL COUPLINGS   ,90(1H$))
 9055 FORMAT(30H0COUPLING.GT.NO. OF CHANNELS  )
 9056 FORMAT(20H0ICODE.GT.8, ICODE= ,I2)
 9901 FORMAT(1H0,7HCHANNEL,I3,15H  DIAGONAL K**2)
 9902 FORMAT(1H0,7HCHANNEL,I3,10H - CHANNEL,I3,10H  COUPLING)
 9900 FORMAT(F6.2,10E12.3)
      END
      SUBROUTINE KINEX(AKIN,W,ECMT,DRT,I,KMAX,ICON)
      IMPLICIT REAL*8(A-H,O-Z)
C     SQRT(XXX)=DSQRT(XXX)
C     ABS(XXX)=DABS(XXX)
C     USE THE NEXT 1 CARDS FOR 32-BIT MACHINES
C modified by (JE,ME)
      DOUBLE PRECISION Q,ECM1,FM1,FM2,FMT,W1,W2,T1
      DIMENSION AKIN(32),BKIN(32),W(800),ICON(20)
      EQUIVALENCE (EN,BKIN(1)),(FM,BKIN(2)),(Z,BKIN(3)),(FMA,BKIN(4))
     1,(ZA,BKIN(5)),(RY,BKIN(6)),(FS,BKIN(7)),(ECM,BKIN(8))
     2,(FMU,BKIN(9)),(RC,BKIN(10)),(FK2,BKIN(11)),(FK,BKIN(12))
     3,(ETA,BKIN(13)),(AFACT,BKIN(14)),(DRF,BKIN(15)),(FN,BKIN(16))
     4,(PNLOC,BKIN(17)),(BFACT,BKIN(18)),(FACT,BKIN(19))
      DATA HBARC,AMU/197.329,931.502/
C modified by (JE,ME)
      SAVE
C
C
      IREL=ICON(10)
      ICON2=ICON(2)
      READ (5,9000)Q ,FM,Z,FMA,ZA,RY,AC,PNLOC,FS,FN
      IF(I.NE.1) GO TO 50
      EN=Q
      ECM1=EN*FMA/(FM+FMA)
      FMA1=FMA
      Q=0.D0
      GO TO 52
   50 CONTINUE
      Q=Q+FN
   52 CONTINUE
      FMU=FM*FMA/(FM+FMA)
      IF(I.NE.0) GO TO 55
      EN=Q
      ECM=Q
      ECM1=0.D0
      GO TO 56
   55 CONTINUE
      EN=(ECM1+Q )*FM/FMU
      ECM=EN*FMU/FM
   56 CONTINUE
      DRF=DRT*FMA1/FMA
      AFACT=FMA**.333333333
      BFACT=FM **.333333333
      RC=ABS(RY)*AFACT
      IF(RY.LT.0.0) RC=RC+ABS(RY)*BFACT
      IF(IREL.EQ.0) GO TO 75
C
C     IREL=ICON(10).NE.0 GIVES RELATIVISTIC KINEMATICS
C
      IF(I.EQ.0) GO TO 66
      IF(I.EQ.1) GO TO 66
      IF(FN .GT.0.0) GO TO 66
      EN=EN+(ECM1+Q)**2/(2.0*FMA*AMU)
   66 CONTINUE
      E=EN
      FM1=FM*AMU
      FM2=FMA*AMU
      FMT=FM1+FM2
      IF(I.EQ.0) E=EN*(0.5*EN+FMT)/FM2
      T1=SQRT(2.0*E*FM2+FMT*FMT)
      W1=(FMT*FM1+FM2*E)/T1
      W2=FM2*(E+FMT)/T1
      FMU=W1*W2/(W1+W2)
      ECM=T1-FMT
      IF(I.EQ.1) ECM1=ECM
      FACT=2.0*FMU/(HBARC*HBARC)
      FMU=FMU/AMU
      FK2=(E*E+2.0*E*FM1)*(FM2/(T1*HBARC))**2
      GO TO 76
   75 CONTINUE
      FACT=2.0*FMU*AMU/(HBARC*HBARC)
      FK2=ECM*FACT
   76 CONTINUE
      IF(FK2.EQ.0.0) GO TO 80
      FK=SQRT(ABS(FK2))
      ETAK=1.43998*Z*ZA*FACT
      ETA=0.5*ETAK/FK
      GO TO 81
   80 CONTINUE
      FK=0.D0
      ETA=0.D0
      ETAK=0.D0
   81 CONTINUE
      ECMT=ECM1
      Q=ECM-ECM1
      IF((I.EQ.0).AND.(ICON2.GE.3)) GO TO 82
      WRITE(6,9012)EN,RY,FS
      WRITE(6,9013)FM,FMA,Q
      WRITE(6,9014)Z,ZA,DRF,PNLOC
      RHO=FK*RC
      WRITE(6,9015)ECM,RC,RHO
      DRHOC=DRF*FK
      WRITE(6,9016)FK,ETA,DRHOC
   82 CONTINUE
      R=0.D0
      RC2=RC*RC
      F1=0.5*ETAK/RC
      DO 100 M=1,KMAX
      R=R+DRF
      MK=M+M-1
      IF(R.GT.RC) GO TO 90
      FC=F1*(3.0-R*R/RC2)
      GO TO 95
   90 CONTINUE
      FC=ETAK/R
   95 CONTINUE
      IF(I.EQ.0) GO TO 98
      W(MK  )=FK2-FC
      W(MK+1)=0.D0
      GO TO 99
   98 CONTINUE
      W(MK+1)=FK2-FC
      W(MK  )=0.D0
   99 CONTINUE
  100 CONTINUE
C
C
      DO 110 J=1,32
      AKIN(J)=BKIN(J)
  110 CONTINUE
      RETURN
 9000 FORMAT(10F8.4)
 9012 FORMAT(18H0INPUT DATA       ,9H   ELAB =,F9.4,9H   RC0  =,F9.4,
     19H   2*S  =,F9.4)
 9013 FORMAT(18H                  ,9H   MASSP=,F9.4,9H   MASST=,F9.4,
     19H   Q    =,F9.4)
 9014 FORMAT(18H                  ,9H   ZP   =,F9.4,9H   ZT   =,F9.4,
     19H   DR   =,F9.4,9H   PNLOC=,F9.4)
 9015 FORMAT(18H0DERIVED DATA     ,9H   ECM  =,F9.4,9H   RC   =,F9.4,
     19H   RHOC =,F9.4)
 9016 FORMAT(18H                  ,9H   K    =,F9.4,9H   ETA  =,F9.4,
     19H   DRHOC=,F9.4)
      END
      SUBROUTINE POTENT(AKIN,U,V,K,KMAX,ICODET,LAM,LDFRM,ABETA,IERR
     1 ,ICON,IPART)
      IMPLICIT REAL*8(A-H,O-Z)
C     SQRT(XXX)=DSQRT(XXX)
C     ABS(XXX)=DABS(XXX)
C     AMAX1(XXX,YYY)=DMAX1(XXX,YYY)
C     EXP(XXX)=DEXP(XXX)
C     FLOAT(III)=DFLOAT(III)
C     AMIN1(XXX,YYY)=DMIN1(XXX,YYY)
C     DOUBLE PRECISION B
      COMMON /KPROJ/ MU
      COMMON /CORR/ BETAR,BETAI,VBETA,RHO0R,RHO0I,INDDEF,INDIM
      DIMENSION U(800),V(800),UT(5),G(4),ICON(20)
     1,        AKIN(32),BKIN(32),CN(16),YLAM(16),X( 8),W( 8),B(3,17)
     2,CP(16),LDFRM(3),ABETA(3),BETA(3),FACTR(3)
      DIMENSION UB(800),VB(1600)
      EQUIVALENCE (G(1),FN),(G(2),FL),(G(3),FJ2),(G(4),FSS)
      EQUIVALENCE (EN,BKIN(1)),(FM,BKIN(2)),(Z,BKIN(3)),(FMA,BKIN(4))
     1,(ZA,BKIN(5)),(RY,BKIN(6)),(FS,BKIN(7)),(ECM,BKIN(8))
     2,(FMU,BKIN(9)),(RC,BKIN(10)),(FK2,BKIN(11)),(FK,BKIN(12))
     3,(ETA,BKIN(13)),(AFACT,BKIN(14)),(DRF,BKIN(15))
     4,(PNLOC,BKIN(17)),(BFACT,BKIN(18)),(FACT,BKIN(19))
     5,(BIF,BKIN(28)),(G,BKIN(29))
      DATA B(1,1),B(2,1),B(3,1)/6H NX=1 ,6HVOLUME,6H W-S  /
      DATA B(1,2),B(2,2),B(3,2)/6H NX=2 ,6HSURFAC,6HE W-S /
      DATA B(1,3),B(2,3),B(3,3)/6H NX=3 ,6H2ND DE,6HRIV   /
      DATA B(1,4),B(2,4),B(3,4)/6H NX=4 ,6HL.S VO,6HLUME  /
      DATA B(1,5),B(2,5),B(3,5)/6H NX=5 ,6HL.S SU,6HRFACE /
      DATA B(1,6),B(2,6),B(3,6)/6H NX=6 ,6HVOL*R*,6H*POWR /
      DATA B(1,7),B(2,7),B(3,7)/6H NX=7 ,6HSUR*R*,6H*POWR /
      DATA B(1,8),B(2,8),B(3,8)/6H NX=8 ,6H2DER*R,6H**PWR /
      DATA B(1,9),B(2,9),B(3,9)/6H NX=9 ,6HHARMON,6HIC OSC/
      DATA B(1,10),B(2,10),B(3,10)/6H NX=10,6H GAUSS,6HIAN   /
      DATA B(1,11),B(2,11),B(3,11)/6H NX=11,6H DEFOR,6HM VOL /
      DATA B(1,12),B(2,12),B(3,12)/6H NX=12,6H DEFOR,6HM SURF/
      DATA B(1,13),B(2,13),B(3,13)/6H NX=13,6H DEFOR,6HM 2DER/
      DATA B(1,14),B(2,14),B(3,14)/6H NX=14,6H EXTER,6H FORMF/
      DATA B(1,15),B(2,15),B(3,15)/6H NX=15,6H FOLD ,6HFORMF /
      DATA B(1,16),B(2,16),B(3,16)/6H NX=16,6H WATSO,6HN L.S /
      DATA B(1,17),B(2,17),B(3,17)/6H NX= O,6HUT OF ,6HRANGE /
      DATA X/0.0950125098,0.2816035507,0.4580167776,0.6178762444,
     1       0.7554044083,0.8656312023,0.9445750230,0.9894009439/
      DATA W/0.1894506104,0.1826034150,0.1691565193,0.1495959888,
     1       0.1246289712,0.0951585116,0.0622535239,0.0271524594/
      DATA SQRPI,PI/1.772453851,3.141593/
      DATA ETA4,ETA5,ETA6/ 4.D0,10.D0,60.D0/
      DATA F4PI/12.56637061/
C
      IND=0
      INDDEF=0
      INDIM=0
      RHO0R=0.D0
      RHO0I=0.D0
      ICOD=ICODET
      ICODE=IABS(ICODET)
      ICO=ICODE
      IF(ICODE.EQ.0.OR.ICODE.EQ.3.OR.ICODE.EQ.7.OR.ICODE.EQ.8) ICO=0
C
      ICON2=ICON(2)
      IF((IPART.EQ.0).AND.(ICON2.GE.3)) GO TO 40
      WRITE(6,9008)
   40 CONTINUE
      DO 45 M=1,32
      BKIN(M)=AKIN(M)
   45 CONTINUE
      DR=DRF
      RM=0.D0
C
C
C     READ IN CARD SET 5,6,OR 7   POTENTIAL CARDS
C
   70 CONTINUE
      IF(IND.EQ.0) GO TO 71
      READ(5,9000) BETAR,BETAI
      IF(ICOD.EQ.-7) WRITE(6,9002) BETAR,BETAI
      IF(ICOD.EQ.8) WRITE(6,9003) BETAR,BETAI
      BETAR=BETAR/VBETA
      BETAI=BETAI/VBETA
      IND=0
   71 CONTINUE
      READ (5,9000)FZ,VR,RY,AR,VSOR,VI,RZ,AI,VSOI,PWR
      NX=ABS(FZ)
      VRT=VR
      VIT=VI
      IF(NX.LE. 0) NX=17
      IF(NX.GT.17) NX=17
      RR=ABS(RY)*AFACT
      RI=ABS(RZ)*AFACT
      IF(RY.LT.0.0) RR=RR+ABS(RY)*BFACT
      IF(RZ.LT.0.0) RI=RI+ABS(RZ)*BFACT
   74 CONTINUE
      IF((IPART.EQ.0).AND.(ICON2.GE.3)) GO TO 75
      WRITE(6,9509)(B(J,NX),J=1,3),VR,RY,AR,RR,VSOR
      WRITE(6,9510)                VI,RZ,AI,RI,VSOI,PWR
   75 CONTINUE
      RY=ABS(RY)
      RZ=ABS(RZ)
C      NX=1 VOLUME WOODS-SAXON
C      NX=2 SURFACE WOODS-SAXON
C      NX=3 SECOND DERIVATIVE WOODS-SAXON
C      NX=4 L.S POTENTIAL FOR WOODS-SAXON VOLUME
C      NX=5 L.S POTENTIAL FOR WOODS-SAXON SURFACE
C      NX=6 WOOD-SAXON*R**PWR
C      NX=7 1ST-DER WOOD-SAXON*R**PWR
C      NX=8 2ND-DER WOOD-SAXON*R**PWR
C      NX=9 HARMONIC OSCILLATOR
C      NX=10 GAUSSIAN*R**PWR
C      NX=11 DEFORMED VOLUME
C      NX=12 DEFORMED SURFACE
C      NX=13 DEFORMED SECOND DERIVATIVE
C      NX=14 FORMF EXTERNAL FORM FACTOR
C      NX=15 FOLDING FORM FACTOR
C      NX=16 WATSON L.S POTENTIAL
      KFLAG=0
      VR=VR*FACT
      VI=VI*FACT
      IF(ICODE.EQ.1) GO TO 81
      IF(ICODE.EQ.2) GO TO 81
      KT=FK*MAX(RR,RI)+ETA5
      LPLUS=MAX0(LPLUS,KT)
   80 CONTINUE
      KT=(2.3*ETA4*MAX(AR,AI)+MAX(RR,RI))/DR
      GO TO 82
   81 CONTINUE
      RM=MAX(RR,RM)
      IF(EN.GT.0.0) GO TO 80
      ALPHB=2.3*ETA4
      TEMP=1.0+ALPHB/10.0
      KT=((ALPHB-ETA*TEMP)/FK+RR)/DR
   82 CONTINUE
      K=MIN0(MAX0(K,KT),KMAX)
   83 CONTINUE
      IF(AR.NE.0.0) GO TO 85
      F1=0.D0
      F2=0.D0
      GO TO 86
   85 CONTINUE
      F2=EXP(-DR/AR)
      F1=EXP( RR/AR)
   86 CONTINUE
      IF(AI.NE.0.0) GO TO 95
      F3=0.D0
      F4=0.D0
      GO TO 96
   95 CONTINUE
      F4=EXP(-DR/AI)
      F3=EXP( RI/AI)
   96 CONTINUE
      IF(IPART.NE.0) GO TO 98
      IF(ICO  .NE.0) GO TO 98
      IF(ICODE.EQ.7.OR.ICODE.EQ.8) GO TO 97
      IF(NX.GT.5) GO TO 98
   97 CONTINUE
      IF(ICODE.EQ.7.AND.NX.LT.12) GO TO 3010
      IF (ICODE .EQ. 7 .AND. NX .GT. 15) GOTO 3010
      IF(ICODE.EQ.8.AND.NX.LT.12) GO TO 3020
      IF(ICODE.EQ.8.AND.NX.GT.15) GO TO 3020
      IF(NX.EQ.14.OR.NX.EQ.15) GO TO 98
      PWR=PWR+1.D0
      IF(AR.NE.0.0) VR=VR*(RR/AR)**PWR
      IF(AI.NE.0.0) VI=VI*(RI/AI)**PWR
   98 CONTINUE
      IF(NX.GT.16) GO TO 2009
C  THIS SECTION IS ADDED TO CALCULATE VOLUME INTEGRALS OF THE REAL AND
C  IMAGINARY PARTS OF THE OMP TO USE FOR NORMALIZATION OF THE MULTIPOLE
C  MOMENTS OF THESE POTENTIALS AS CALCULATED IN SUBROUTINE ACHK. SINCE
C  THIS IS NOT DONE PROPERLY FOR DEFORMED POTENTIALS, IT GIVES ONLY
C  APPROXIMATE VALUES FOR THESE CASES. ALSO IN THE CASE OF MORE THAN ONE
C  IMAGINARY POTENTIAL THE MULTIPOLE MOMENTS ARE NOT CORRECT. FOR ALL
C  THESE CASES USE PROGRAM  BELGEN  TO OBTAIN THE CORRECTLY NORMALIZED
C  MULTIPOLE MOMENTS.
      IF(ICO.NE.0.OR.IPART.NE.0) GOTO 12
      IF(ABS(RHO0I).GT.1.0E-7) THEN
        INDIM=1
        GOTO 12
      ENDIF
      F5=F1
      F6=F3
      IF(LAM.NE.0) GOTO 2
      IF(ICODE.NE.3) GOTO 1
      IF(PWR.NE.0.0) GOTO 12
      MNX=NX-5
      VRT=VRT/3.
      VIT=VIT/3.
      GOTO 7
    1 IF(ICODE.NE.7) GOTO 12
      MNX=NX-11
      GOTO 7
    2 IF(LAM.NE.1) GOTO 4
      IF(ICODE.NE.3) GOTO 3
      IF(PWR.NE.1.0) GOTO 12
      MNX=NX-5
      VRT=VRT*RR/10.
      VIT=VIT*RI/10.
      GOTO 7
    3 IF(ICODE.NE.8) GOTO 12
      MNX=NX-11
      GOTO 7
    4 IF(ICODE.NE.0) GOTO 5
      IF(NX.EQ.2.OR.NX.EQ.11) MNX=1
      IF(NX.EQ.3.OR.NX.EQ.12) MNX=2
      GOTO 7
    5 IF(ICODE.NE.3) GOTO 6
      MNX=NX-1
      GOTO 7
    6 MNX=NX-11
    7 CONTINUE
      IF(NX.GE.11) INDDEF=1
      R=0.D0
      DO 11 M=1,K
      R=R+DR
      R2=R*R
      F5=F5*F2
      F6=F6*F4
      XTEMP=1.D0
      IF(M.EQ.1.OR.M.EQ.K) XTEMP=0.5D0
      IF(MNX.EQ.1) GOTO 8
      IF(MNX.EQ.2) GOTO 9
      GOTO 11
    8 IF(VR.NE.0.) RHO0R=RHO0R+XTEMP*VRT*R2*F5/(1.0+F5)
      IF(VI.NE.0.) RHO0I=RHO0I+XTEMP*VIT*R2*F6/(1.0+F6)
      GOTO 11
    9 IF(VR.NE.0.) RHO0R=RHO0R-XTEMP*VRT*R2*F5/(1.0+F5)**2
      IF(VI.NE.0.) RHO0I=RHO0I-XTEMP*VIT*R2*F6/(1.0+F6)**2
   11 CONTINUE
   12 CONTINUE
C   
      GO TO (100,200,300,400,500,600,700,800,900,1000,
     1  1100,1200,1300,1400,1500,1600), NX
C
C     VOLUME WOODS SAXON
C
  100 CONTINUE
      DO 160 M=1,K
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )-VR*F1/(1.0+F1)
      U(MK+1)=U(MK+1)-VI*F3/(1.0+F3)
  160 CONTINUE
      GO TO 2000
C
C     1ST DER WOODS SAXON
C
  200 CONTINUE
      DO 260 M=1,K
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )+VR*F1/(1.0+F1)**2
      U(MK+1)=U(MK+1)+VI*F3/(1.0+F3)**2
  260 CONTINUE
      GO TO 2000
C
C     2ND DER WOODS SAXON
C
  300 CONTINUE
      DO 360 M=1,K
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )-VR*F1*(1.0-F1)/(1.0+F1)**3
      U(MK+1)=U(MK+1)-VI*F3*(1.0-F3)/(1.0+F3)**3
  360 CONTINUE
      GO TO 2000
C
C     L.S VOLUME WOODS SAXON
C
  400 CONTINUE
      IF(AR.NE.0.0) VR=VR/AR
      R=0.D0
      DO 460 M=1,K
      R=R+DR
      F1=F1*F2
      V(M)=V(M)-VR*F1/(R*(1.0+F1)**2)
  460 CONTINUE
      GO TO 2000
C
C     L.S 1ST DER WOODS SAXON
C
  500 CONTINUE
      IF(AR.NE.0.0) VR=VR/AR
      SPO=0.D0
      IF(VSOR.GT.0.0)SPO=1.D0
      R=0.D0
      DO 560 M=1,K
      R=R+DR
      F1=F1*F2
      TEMP=1.D0+F1
      V(M)=V(M)+VR*F1*((1.-F1)/TEMP+SPO*AR/R)/(R*TEMP*TEMP)
  560 CONTINUE
      GO TO 2000
C
C     WOOD SAXON*R**PWR
C
  600 CONTINUE
      R=0.D0
      DO 660 M=1,K
      R=R+DR
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      RPWR=R**PWR
      U(MK  )=U(MK  )-VR*F1*RPWR/(1.0+F1)
      U(MK+1)=U(MK+1)-VI*F3*RPWR/(1.0+F3)
  660 CONTINUE
      GO TO 2000
C
C      1ST DER WOOD SAXON*R**PWR
C
  700 CONTINUE
      R=0.D0
      DO 760 M=1,K
      R=R+DR
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      RPWR=R**PWR
      U(MK  )=U(MK  )+VR*F1*RPWR/(1.0+F1)**2
      U(MK+1)=U(MK+1)+VI*F3*RPWR/(1.0+F3)**2
  760 CONTINUE
      GO TO 2000
C
C     2ND DER WOODS-SAXON*R**PWR
C
  800 CONTINUE
      R=0.D0
      DO 860 M=1,K
      R=R+DR
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      RPWR=R**PWR
      U(MK  )=U(MK  )-VR*RPWR*F1*(1.0-F1)/(1.0+F1)**3
      U(MK+1)=U(MK+1)-VI*RPWR*F3*(1.0-F3)/(1.0+F3)**3
  860 CONTINUE
      GO TO 2000
C
C     HARMONIC OSCILLATOR
C
  900 CONTINUE
      READ (5,9000)G
      WRITE(6,9508)G
      RY=ABS(RY)
      F1=1.0/RY**2
      F2=F1/RY
      F3=0.5D0
      F4=SQRPI*0.5
      L=FL
      IF(L.EQ.0) GO TO 931
      DO 930 I=1,L
      F3=F3+1.D0
      F4=F4*F3
      F2=F2*F1
  930 CONTINUE
  931 CONTINUE
      NN=FN
      T1=1.D0
      T2=F4
      CN(1)=(-1.0)**NN
      IF(NN.EQ.0) GO TO 941
      DO 940 I=1,NN
      F3=F3+1.D0
      T1=T1*FLOAT(I)
      T2=T2*F3
      CN(I+1)=-CN(I)*F1*FLOAT(NN+1-I)/(FLOAT(I)*F3)
  940 CONTINUE
  941 CONTINUE
      T2=SQRT(2.0*F2*T2/T1)/F4
      IF(VR.NE.0.0) T2=T2*VR
      KT=10.0*RY/DR
      KT=MIN0(KT,K)
      R=0.D0
      F1=F1/2.0
      DO 960 M=1,KT
      MK=M+M-1
      R=R+DR
      R2=R**2
      F2=CN(1)
      IF(NN.EQ.0) GO TO 951
      F3=1.D0
      DO 950 I=1,NN
      F3=F3*R2
      F2=F2+CN(I+1)*F3
  950 CONTINUE
  951 CONTINUE
      U(MK  )=U(MK  )+F2*T2*EXP(-F1*R2)*R**L
  960 CONTINUE
      GO TO 2000
C
C     EXP(-((R-R0)/A)**2)*R**POWR
C
 1000 CONTINUE
      IF(VR.EQ.0.0) GO TO 1061
      IF(AR.EQ.0.0) GO TO 1061
      M1=NINT((RR-6.0*AR)/DR)+1
      M1=MAX0(M1,1)
      M2=NINT((RR+6.0*AR)/DR)+1
      M2=MIN0(M2,K)
      R=FLOAT(M1-1)*DR
      DO 1060 M=M1,M2
      MK=M+M-1
      R=R+DR
      TEMP=((R-RR)/AR)**2
      U(MK  )=U(MK  )-VR*EXP(-TEMP)*R**PWR
 1060 CONTINUE
 1061 CONTINUE
      IF(VI.EQ.0.0) GO TO 1071
      IF(AI.EQ.0.0) GO TO 1071
      M1=NINT((RI-6.0*AI)/DR)+1
      M1=MAX0(M1,1)
      M2=NINT((RI+6.0*AI)/DR)+1
      M2=MIN0(M2,K)
      R=FLOAT(M1-1)*DR
      DO 1070 M=M1,M2
      MK=M+M-1
      R=R+DR
      TEMP=((R-RI)/AI)**2
      U(MK+1)=U(MK+1)-VI*EXP(-TEMP)*R**PWR
 1070 CONTINUE
 1071 CONTINUE
      GO TO 2000
C
C   THE FOLLOWING STATEMENTS ARE INCLUDED TO TREAT CORRECTION TERMS FOR
C   BETA AND OCTUPOLE VIBRATIONS
C
 1099 CONTINUE
      IND=1
      VR=VR*BETAR
      VI=VI*BETAI
C
C
C     DEFORMED VOLUME OR SURFACE BY YLM EXPANSION
C     ALSO DEFORMED THE POTENTIALS FOR COUPLINGS OF BETA, GAMMA
C     AND OCTUPOLE BANDS
C
 1100 CONTINUE
      IF(VBETA.EQ.1.0.AND.IND.EQ.0) GO TO 1171
      DO 1101 I=1,3
      BETA(I)=ABETA(I)
      FACTR(I)=SQRT(DFLOAT(LDFRM(I)+LDFRM(I)+1))*BETA(I)/(SQRPI*2.0)
 1101 CONTINUE
      LMAX=MAX0(LAM,LDFRM(3))+1
      LTEMP=2
      IF(ICOD.EQ.8) LTEMP=3
      IF(IND.EQ.1.AND.ICOD.EQ.-7) LTEMP=0
      IF(IND.EQ.1.AND.ICOD.EQ. 8) LTEMP=1
      IF(ICODE.EQ.7.OR.ICODE.EQ.8) LMAX=MAX0(LAM+LTEMP,LDFRM(3))+1
      DO 1103 I=1,8
      YLAM(I  )=0.D0
 1103 YLAM(I+8)=0.D0
      T2=(-1.0)**LAM
      T1=(-1.0)**LDFRM(1)
      DO 1109 I=1,8
      CN(I  )=0.D0
      CN(I+8)=0.D0
      P2=0.D0
      P3=1.D0
      Z=X(I)
      DO 1108 M=1,LMAX
      L=M-1
      FL=L-1
      IF(L.EQ.0) GO TO 1102
      P3=((2.0*FL+1.0)*Z*P2-FL*P1)/(FL+1.0)
 1102 CONTINUE
      DO 1104 J=1,3
      IF(L.NE.LDFRM(J)) GO TO 1104
      CN(I  )=CN(I  )+P3*FACTR(J)
      CN(I+8)=CN(I+8)+P3*FACTR(J)*T1
 1104 CONTINUE
      IF(ICODE.EQ.8.OR.ICODE.EQ.7) GO TO 1106
      GO TO 1105
 1106 CONTINUE
      IF(IABS(L-LAM).GT.LTEMP.OR.L+LAM.LT.LTEMP) GO TO 1107
      T3=SQRT(DFLOAT((LTEMP+LTEMP+1)*(L+L+1))/DFLOAT(LAM+LAM+1))*0.5
     1/SQRPI*VCC(L+L,LTEMP+LTEMP,LAM+LAM,0,0)*VCC(L+L,LTEMP+LTEMP,
     2LAM+LAM,0,MU)*P3*W(I)*SQRT(DFLOAT(L+L+1))*SQRPI
      YLAM(I  )=YLAM(I  )+T3
      YLAM(I+8)=YLAM(I+8)+T3
      GO TO 1107
 1105 IF(L.NE.LAM) GO TO 1107
      YLAM(I  )= P3*W(I)*SQRT(DFLOAT(L+L+1))*SQRPI
      YLAM(I+8)=YLAM(I  )*T2
 1107 CONTINUE
      P1=P2
      P2=P3
 1108 CONTINUE
 1109 CONTINUE
      DO 1170 I=1,2
      GO TO (1110,1120),I
 1110 CONTINUE
      IF(VR.EQ.0.0) GO TO 1170
      VX=VR
      RX=RR
      AX=AR
      F1=1.D0
      FTEMP=F2
      IFL=-1
      GO TO 1130
 1120 CONTINUE
      IF(VI.EQ.0.0) GO TO 1170
      VX=VI
      RX=RI
      AX=AI
      F1=1.D0
      FTEMP=F4
      IFL=0
 1130 CONTINUE
      DO 1135 J=1,16
      CP(J)=EXP((1.D0+CN(J))*RX/AX)
 1135 CONTINUE
      VX=VX/BETA(1)
      IF(LAM.EQ.0.AND.ICODE.NE.7) VX=-VX*BETA(1)/(2.0*SQRPI)
      JT=NX-10
      DO 1160 M=1,K
      MK=M+M+IFL
      VTEMP=0.D0
      F1=F1*FTEMP
      DO 1155 J=1,16
      F3=F1*CP(J)
      GO TO (1140,1145,1150), JT
 1140 CONTINUE
      VTEMP=VTEMP+YLAM(J)*VX*F3/(1.0+F3)
      GO TO 1155
 1145 CONTINUE
      VTEMP=VTEMP-YLAM(J)*VX*F3/(1.0+F3)**2
      GO TO 1155
 1150 CONTINUE
      VTEMP=VTEMP+YLAM(J)*VX*F3*(1.0-F3)/(1.0+F3)**3
 1155 CONTINUE
      U(MK  )=U(MK  )+VTEMP
 1160 CONTINUE
 1170 CONTINUE
 1171 CONTINUE
      IF(ICOD.EQ.7.OR.(ICODE.LT.7.OR.ICODE.GT.8)) GO TO 2000
      IF(IND.EQ.0) GO TO 1099
      GO TO 2000
 1200 CONTINUE
      GO TO 1100
 1300 CONTINUE
      GO TO 1100
C
C     EXTERNAL FORM FACTOR
C
 1400 CONTINUE
      READ(2,9000) F1,F2
      WRITE(6,9513) F1,F2
      F3=VR
      IGT=F2+1.D0
      IF(MOD(IGT,2).EQ.0) F3=VI
      IF(F3.EQ.0) F3=1.D0
      KT=F1
      IF(IGT.GT.2) GO TO 1440
      MK=IGT
      DO 1430 M=1,KT,5
      READ(2,9100) UT
      IF(M.GT.KMAX) GO TO 1430
      DO 1420 I=1,5
      U(MK  )=U(MK  )-UT(I)*F3
      MK=MK+2
 1420 CONTINUE
 1430 CONTINUE
      GO TO 1480
 1440 CONTINUE
      MK=IGT-2
      DO 1460 M=1,KT,5
      READ(5,9100) UT
      IF(M.GT.KMAX) GO TO 1460
      IF(MK.GE.2) GO TO 1460
      DO 1450 I=1,5
      V(MK  )=V(MK  )-UT(I)*F3
      MK=MK+2
 1450 CONTINUE
 1460 CONTINUE
 1480 CONTINUE
      K=MAX0(K,KT)
      K=MIN0(K,KMAX)
      GO TO 2000
C
C     FOLDING FORM FACTOR
C
 1500 CONTINUE
      READ(5,9000) F1,QCODE,FMUV,OPT
      F3=OPT
      KT=F1
      DRX=DRF
      VOPT=OPT*F4PI
      IQODE=QCODE
      KM=400
      MINL=LAM+1
      CALL FMN(UB,KT,DRX)
      IF(IQODE.NE.6) GO TO 1501
      DO 1502 I=1,KT
      MK=I+I-1
      U(MK)=U(MK)+UB(I)*OPT
 1502 CONTINUE
      GO TO 1580
 1501 CONTINUE
      CALL SLATR(KT,DRX,VB,MINL,FMUV,IQODE,KM)
      DO 1570 I=1,KT
      VB(I+ 800)=VB(I    )
      VB(I+1200)=VB(I+400)
 1570 CONTINUE
      WRITE(6,9001) IQODE,FMUV,OPT
      CALL RADIN(KT,DRX,VB,UB,SL,SK,VOPT,U(1),IQODE,MINL,FMUV,KM)
 1580 CONTINUE
      K=MAX0(K,KT)
      K=MIN0(K,KMAX)
      GO TO 2000
C
C     WATSON L.S POTENTIAL
C
 1600 CONTINUE
      IF(AR.NE.0.0) VR=VR/AR
      DO 1660 M=1,K
      F1=F1*F2
      V(M)=V(M)-VR*F1/(RR*(1.0+F1)**2)
 1660 CONTINUE
 2000 CONTINUE
      IF(KFLAG.NE.0) GO TO 2009
      IF(VSOR.EQ.0.0) GO TO 2009
      NX=NX+3
      IF(NX.GT.5) GO TO 2009
      KFLAG=1
      VR=VR*VSOR/45.2
      GO TO 83
 2009 CONTINUE
      IF(FZ.GT.0.0) GO TO 70
      IF(ICODE.EQ.0) GO TO 3000
      IF(ICODE.EQ.3) GO TO 3000
      IF(ICODE.EQ.4) GO TO 3000
      IF(ICODE.EQ.7) GO TO 3000
      IF(ICODE.EQ.8) GO TO 3000
C
C     SINGL PARTICLE ORBITAL
C
C
C     READ IN QUANTUM NUMBERS FO SINGLE PARTICLE ORBITAL
C
      READ(5,9000)G,VTRIAL,FISW,DAMP
      ISW=FISW
      WRITE(6,9500)
      FJ=FJ2/2.0
      IF(VTRIAL.EQ.0.0) VTRIAL=ETA6
      WRITE(6,9508)G,FISW,VTRIAL
      IF(DAMP.NE.0.0) WRITE(6,9511) DAMP
      FSS=FSS/2.0
      FACT=(FJ**2+FJ-FL**2-FL-FSS**2-FSS)*0.5
      DO 2028 M=1,KMAX
      MK=M+M-1
      U(MK  )=U(MK  )+V(M)*FACT
 2028 CONTINUE
      CALL BIND(U,V,DR,RM,FN,FL,K,FK,ETA,VTRIAL,ECM,FK2,ISW,IERR)
      BIF=RM/DR
      DO 2045 M=28,32
      AKIN(M)=BKIN(M)
 2045 CONTINUE
      FACT=PNLOC*PNLOC/4.0
C
C     NON-LOCAL CORRECTION FOR SINGLE PARTICLE FUNCTION
C
      SUM=0.D0
      R=0.D0
      DO 2075 M=1,K
      MK=M+M-1
      R=R+DR
      U(MK+1)=VTRIAL*U(MK  )+U(MK+1)
      TEMP=FACT*(U(MK+1)-FK2)
      IF(PNLOC) 2072,2075,2070
 2070 CONTINUE
      V(M)=V(M)/SQRT(1.0+TEMP)
      GO TO 2074
 2072 CONTINUE
      V(M)=V(M)*EXP(-TEMP/2.0)
 2074 CONTINUE
      SUM=SUM+(V(M)*R)**2
 2075 CONTINUE
      IF(FACT.EQ.0.0) SUM=1.D0
      IF(FACT.EQ.0.0) GO TO 2076
      SUM=1.0/SQRT(SUM*DR)
 2076 CONTINUE
      IF(DAMP.EQ.0.0) GO TO 2090
C
C     APPLY DAMPING FACTOR, EXP(-DAMP*R*R) TO FORM FACTOR
C
      R=0.D0
      H=ABS(DAMP)
      DO 2085 M=1,K
      R=R+DR
      F2=H*R*R
      F2=-MIN(F2,30.D0)
      F1=EXP(F2)
      V(M)=V(M)*F1
 2085 CONTINUE
 2090 CONTINUE
      DO 2100 M=1,KMAX
      MK=M+M-1
      IF(M.GT.K) V(M)=0.D0
      U(MK  )=V(M)*SUM
      V(M)=U(MK+1)
      U(MK+1)=0.D0
 2100 CONTINUE
 3000 CONTINUE
      IF(ABS(RHO0R).LT.1.0E-7) RHO0R=1.D0
      IF(ABS(RHO0I).LT.1.0E-7) RHO0I=1.D0
      RETURN
 3010 WRITE(6,9501)
      STOP
 3020 WRITE(6,9502)
      STOP
 9000 FORMAT(10F8.4)
 9001 FORMAT(18H0PARAMETERS       ,9H   QCODE=,F9.4,9H   RANGE=,F9.4,7H
     1  V0 =,F11.4)
 9002 FORMAT(/21X,'BETA0R =',F7.4,'   BETA0I =',F7.4,/)
 9003 FORMAT(/21X,'BETA1R =',F7.4,'   BETA1I =',F7.4,/)
 9008 FORMAT(21H0POTENTIAL PARAMETERS)
 9010 FORMAT( 9H0PARTICLE,I2,120(1H*))
 9100 FORMAT(5E16.7)
 9500 FORMAT(1H ,3A6,5(3X,A6,F9.4))
 9501 FORMAT('0  XXX  INCORRECT GAMMA.OR.BETA BAND COUPLING POTENTIAL  X
     1XX')
 9502 FORMAT('0  XXX  INCORRECT OCTUPOLE BAND COUPLING POTENTIAL  XXX')
 9508 FORMAT(18X,9H   NODES=,F9.4,9H   L    =,F9.4,9H   2*J  =,F9.4,
     19H   2*S  =,F9.4,9H   FISW =,F9.4,9H   VTRL =,F9.4)
 9509 FORMAT(3A6,9H   V RL =,F9.4,9H   R0RL =,F9.4,9H   A RL =,F9.4,
     19H   R RL =,F9.4,9H   VSOR =,F9.4)
 9510 FORMAT(18X,9H   V IM =,F9.4,9H   R0IM =,F9.4,9H   A IM =,F9.4,
     19H   R IM =,F9.4,9H   VSOI =,F9.4,9H   POWR =,F9.4)
 9511 FORMAT(18X,9H   DAMP =,F9.4)
 9513 FORMAT(18X,9H   F1   =,F9.4,9H   F2   =,F9.4)
      END
      SUBROUTINE FNLOC(U,W,VB,LTR,BKIN,FNRNG,KT,IJ,AKIN,ICODE,KMAX)
      IMPLICIT REAL*8(A-H,O-Z)
C     SIN(XXX)=DSIN(XXX)
C     COS(XXX)=DCOS(XXX)
C     EXP(XXX)=DEXP(XXX)
C     ABS(XXX)=DABS(XXX)
      DIMENSION U(800,2),W(800),VB(400),BKIN(8,18),AKIN(18),IJ(2)
     1,RT(8,7),FM(3),FMA(3),PNLOC(3),DR(3),FMU(3),FK2(3)
      DIMENSION FACTNL(2)
      DATA  RT/  0.0   ,  4.4934,  7.7253, 10.9041, 14.0662, 17.2208,
     1 20.3713, 23.5195,  2.0816,  5.9404,  9.2058, 12.4044, 15.5792,
     2 18.7426, 21.8997, 25.0528,  3.3421,  7.2899, 10.6139, 13.8461,
     3 17.0429, 20.2219, 23.3905, 26.5526,  4.5141,  8.5838, 11.9727,
     4 15.2445, 18.4681, 21.6666, 24.8501, 28.0239,  5.6467,  9.8404,
     5 13.2956, 16.6093, 19.8624, 23.0828, 26.2833, 29.4706,  6.7565,
     6 11.0702, 14.5906, 17.9472, 21.2311, 24.4748, 27.6937, 30.8960,
     7  7.8511, 12.2793, 15.8632, 19.2627, 22.5781, 25.8461, 29.0843,
     8 32.3025/
C
C     FINITE RANGE CORRECTION
C
C
C     POSITIVE FR PARAMETER CALCULATES A HULTHEN TYPE FR CORRECTION
C     NEGATIVE FR PARAMETER CALCULATES A EXPONENTIAL TYPE FR CORRECTION
C
      KMAX2=KMAX+KMAX
      K=KT
      DO 20 I=1,2
      J=IJ(I)
      FM(I)   =BKIN(J,2)
      FMA(I)  =BKIN(J,4)
      PNLOC(I)=BKIN(J,17)
      DR(I)   =BKIN(J,15)
      FK2(I)  =BKIN(J,11)
      FMU(I)  =BKIN(J,9)
   20 CONTINUE
      FM(3)   =AKIN(2)
      FMA(3)  =AKIN(4)
      PNLOC(3)=AKIN(17)
      DR(3)   =AKIN(15)
      FK2(3)  =AKIN(11)
      FMU(3)  =AKIN(9)
C
C     CHANGE CHANNELS I,J INTO 1,2
C
      IF(IJ(1).GT.IJ(2)) GO TO 35
      IJ(1)=1
      IJ(2)=2
      GO TO 40
   35 CONTINUE
      IJ(1)=2
      IJ(2)=1
   40 CONTINUE
      IF(FNRNG.EQ.0.0) GO TO 348
      IF(FM(1)-FM(2))143,348,142
  142 I1=1
      I2=2
      J1=IJ(1)
      J2=IJ(2)
      GO TO 144
  143 I1=2
      I2=1
      J1=IJ(2)
      J2=IJ(1)
C
C     I1 IS FOR LARGER MASS , I2 IS FOR SMALLER MASS
C
  144 CONTINUE
      FM1=FMU(I1)
      FM2=FMU(I2)
      TEMP=FM(I2)*FNRNG**2/FM(I1)
      FACT=TEMP*FMA(I1)/FMA(I2)
      FMX=FM(I1)-FM(I2)
      TEMP=TEMP*FMX
      IF((FMX.LT.1.5).OR.(ICODE.NE.2)) GO TO 260
      FNRNG=-ABS(FNRNG)
C
C     FORM FACTOR PART BY ROST-KUNZ METHOD
C     USED FOR MULTI-NUCLEON TRANSFER REACTIONS
C
      DR2=DR(3)**2
      FLF=LTR*(LTR+1)
      R3=DR(3)*FLOAT(K+1)
      D3=0.D0
      A3=0.D0
      KM=K+K+1
      R2=R3-DR(3)
      D2=-FLF/R3
      A2=W (KM-2)*R2
      DO 220 M=1,K
      KM=KM-2
      R1=R2-DR(3)
      A1=W (KM-2)*R1
      D1=-FLF/(R2*R2)
      IF(A3.EQ.0.0) GO TO 215
      D1=D1+(A3-2.0*A2+A1)/(A2*DR2)
      W (KM+2)=EXP(D2*FACT)*W (KM+2)
      DX2=(D3*A3-2.0*D2*A2+D1*A1)/(R2*DR2)-A2*D2*D2/R2
      W (KM+2)=W (KM+2)+0.5*DX2*FACT*FACT
  215 CONTINUE
      IF(ABS(W (KM-2)).LT.ABS(W (KM  ))) GO TO 225
      R3=R2
      R2=R1
      D3=D2
      D2=D1
      A3=A2
      A2=A1
  220 CONTINUE
  225 CONTINUE
      KM=KM/2+1
      DO 227 M=1,KM
      MK=M+M-1
      VB(M)=W(MK)
      W(MK)=0.D0
  227 CONTINUE
      RKM=FLOAT(KM)*DR(3)
      LL=MIN0(LTR,6)
      FLF=LL*(LL+1)
      DR12=DR2/12.0
      DO 250 N=1,8
      FKX=(RT(N,LL+1)/RKM)**2
      A1=0.D0
      D1=0.D0
      RF=DR(3)
      A2=DR(3)
      D2=(1.0+DR12*(FKX-FLF/(RF*RF)))*A2
      FNORM=0.D0
      AN=0.D0
      DO 230 M=1,KM
      AN=AN+(VB(M)*A2*RF)
      FNORM=FNORM+A2*A2
      VB(M+KMAX2)=A2/RF
      RF=RF+DR(3)
      D3=1.D0+DR12*(FKX-FLF/(RF*RF))
      A3=(12.0*A2-10.0*D2-D1)/D3
      D1=D2
      D2=D3*A3
      A1=A2
      A2=A3
  230 CONTINUE
      AN=AN-0.5*VB(KM)*A1*(RF-DR(3))
      FNORM=FNORM-0.5*A1*A1
      FNORM=AN*EXP(-FACT*FKX)/FNORM
      DO 240 M=1,KM
      MK=M+M-1
      W (MK)=W (MK)+VB(M+KMAX2)*FNORM
  240 CONTINUE
  250 CONTINUE
      GO TO 300
  260 CONTINUE
      IF(FNRNG.GT.0.0) GO TO 300
C
C     SINGLE PARTICLE PART, GAUSSIAN FACTOR
C
      DO 270 M=1,K
      MK=M+M-1
      FACT=EXP(-TEMP*VB(M)/FMU(3))
      W(MK  )=W(MK  )*FACT
      W(MK+1)=W(MK+1)*FACT
  270 CONTINUE
  300 CONTINUE
C
C     DISTORTED WAVE PART
C
      DO 345 M=1,K
      MK=M+M-1
      CTEMP1=TEMP*(U(MK  ,J2)/FM2-U(MK  ,J1)/FM1)
      CTEMP2=TEMP*(U(MK+1,J2)/FM2-U(MK+1,J1)/FM1)
      IF(FNRNG.LT.0.0) GO TO 320
      UT1=CTEMP1+TEMP*VB(M)/FMU(3)+1.0
      UT2=CTEMP2
      DET=UT1*UT1+UT2*UT2
      UT1= UT1/DET
      UT2= UT2/DET
      GO TO 330
  320 CONTINUE
      FACT=EXP(-CTEMP1)
      UT1=FACT*COS(CTEMP2)
      UT2=FACT*SIN(CTEMP2)
  330 CONTINUE
      UF1   =UT1*W(MK  )+UT2*W(MK+1)
      UF2   =UT1*W(MK+1)-UT2*W(MK  )
      W(MK  )=UF1
      W(MK+1)=UF2
  345 CONTINUE
  348 CONTINUE
C
C     NON LOCALITY CORRECTION FACTOR FOR DISTORTED WAVES
C
      DO 350 I=1,2
      FACTNL(I)=PNLOC(I)**2/4.0
  350 CONTINUE
      IF(FACTNL(1)+FACTNL(2).EQ.0.0) GO TO 410
      DO 400 M=1,K
      MK=M+M-1
      DO 390 I=1,2
      L=IJ(I)
      IF(FACTNL(I).EQ.0.0) GO TO 390
      CTEMP1=-FACTNL(I)*(U(MK  ,L)-FK2(I))
      CTEMP2=-FACTNL(I)* U(MK+1,L)
      IF(PNLOC(I).LT.0.0) GO TO 370
      CTEMP1=1.0-CTEMP1
      THETA=ATAN(CTEMP2/CTEMP1)/2.0
      RX=SQRT(SQRT(CTEMP1*CTEMP1+CTEMP2*CTEMP2))
      UT1= COS(THETA)/RX
      UT2= SIN(THETA)/RX
      GO TO 380
  370 CONTINUE
      CTEMP1=CTEMP1/2.0
      FACT=EXP(CTEMP1)
      CTEMP2=CTEMP2/2.0
      UT1=FACT*COS(CTEMP2)
      UT2=FACT*SIN(CTEMP2)
  380 CONTINUE
      UF1   =UT1*W(MK  )-UT2*W(MK+1)
      UF2   =UT1*W(MK+1)+UT2*W(MK  )
      W(MK  )=UF1
      W(MK+1)=UF2
  390 CONTINUE
  400 CONTINUE
  410 CONTINUE
      RETURN
      END
      SUBROUTINE CATHEN(LOPT,IOPT,JOPT,K ,KMAX,UB,VB,SI,DRT,ICON,IERR)
      IMPLICIT REAL*8(A-H,O-Z)
C     SQRT(XXX)=DSQRT(XXX)
C     FLOAT(III)=DFLOAT(III)
C     ABS(XXX)=DABS(XXX)
C modified by (JE, ME)
C     USE THE NEXT 1 CARDS FOR 32-BIT MACHINES
      DOUBLE PRECISION VCC,RACAH,WINEJ
      DIMENSION UB(800,2),VB(1600),SI(800),C(10),QNUM(4,2),G(4),IQN(3)
     1,BKIN(32),JI(2),JJ(2),FJT(2),ICON(20)
      EQUIVALENCE (EN,BKIN(1)),(FM,BKIN(2)),(Z,BKIN(3)),(FMA,BKIN(4))
     1,(ZA,BKIN(5)),(RY,BKIN(6)),(FS,BKIN(7)),(ECM,BKIN(8))
     2,(FMU,BKIN(9)),(RC,BKIN(10)),(FK2,BKIN(11)),(FK,BKIN(12))
     3,(ETA,BKIN(13)),(AFACT,BKIN(14)),(DRF,BKIN(15)),(FN,BKIN(16))
     4,(PNLOC,BKIN(17)),(BFACT,BKIN(18)),(FACT,BKIN(19))
     5,(BIF,BKIN(28)),(G,BKIN(29))
     6,(C(1),CNTROL),(C(2),QCODE),(C(3),FMUV),(C(4),OPT)
      DATA SQR4PI,SQRTEN/3.54490780,3.162277660/
C modified by (JE,ME)
      SAVE
C
      ICON2=ICON(2)
      KMAX2=KMAX +KMAX
      KMAX3=KMAX2+KMAX
      IQFLG=0
      KMFLG=0
      IQN(1)=LOPT
      IQN(2)=IOPT
      IQN(3)=JOPT
      KT=KMAX
      TOTAL=0.D0
  100 CONTINUE
      READ(5,9001)C
      CNTRL=ABS(CNTROL)
      IF(CNTRL.EQ.0.0) GO TO 3010
      IF(QCODE)2020,3000,1100
 1100 CONTINUE
C
C     LOOP FOR TWO ORBITALS
C
      ML=1
      DO 2010 I=1,2
      FJT(1)=C(5)
      FJT(2)=C(I+5)
      IF(I.EQ.1) GO TO 2004
      IF(CNTRL.EQ.1.0) GO TO 2005
 2004 CONTINUE
      CALL KINEX(BKIN,UB(1,2),ECMT,DRT,0,KMAX,ICON)
      RAM=(DRF/DRT)**3
      DO 2003 M=1,KMAX
      MK=M+M-1
      VB(MK  )=0.D0
      VB(MK+1)=0.D0
 2003 CONTINUE
C
C
      CALL POTENT(BKIN,UB(1,2),VB,K,KMAX,2,0,0,0.D0,IERR,ICON,0)
      FJ2=G(3)
      IBF=BIF
      IF(KMFLG.EQ.0) KM=IBF
      KMFLG=1
      KT=MIN0(KT,K)
 2005 CONTINUE
C
C     PRINT OUT SINGLE PARTICLE INFORMATION
C
      IF(ICON2.GE.2) GO TO 2007
      IF((I.EQ.2).AND.(CNTRL.EQ.1.0)) GO TO 2007
      WRITE(6,9501) I,CNTRL
      R=DRF
      K2=K+K
      DO 2006 M=1,K2,20
      MK4=MIN0(M+19,K2)
      WRITE(6,9052)R,(UB(N,2),N=M,MK4,2)
      R=R+DRF*10.0
 2006 CONTINUE
 2007 CONTINUE
      DO 2008 M=1,KMAX
      MK=M+M-1
      UB(ML,1)=UB(MK,2)
      ML=ML+1
 2008 CONTINUE
      DO 2009 M=1,4
      QNUM(M,I)=G(M)
 2009 CONTINUE
      JI(I)=FJT(1)
      JJ(I)=FJT(2)
 2010 CONTINUE
      SLP=0.D0
      RMSP=0.D0
      R=0.D0
      DO 2015 M=1,KT
      R=R+DRF
      RSQ=R*R
      TEMP=UB(M     ,1)*UB(M+KMAX,1)
      SLP =SLP +TEMP*RSQ
      RMSP=RMSP+TEMP*RSQ*RSQ
 2015 CONTINUE
      IF(SLP.NE.0.0) RMSP=SQRT(RMSP/SLP)
      SLP=SLP*DRF
 2020 ICODE=ABS(QCODE)
      KT=MIN0(KT,K)
      LTR1=IQN(1)+1
      JLX=QNUM(3,1)
      JPX=QNUM(3,2)
      WRITE(6,9050) JLX,JI(1),JJ(1),JPX,JI(1),JJ(2)
      WRITE(6,9051) QCODE,FMUV,OPT,SLP,RMSP
C     ENTER TIME REVERSAL PHASE
      LVR=QNUM(2,1)+QNUM(2,2)
      OPT=OPT*PHASEF((IQN(1)-LVR)/2)
      IFLAG=1
      IF(ICODE.LE.10) GO TO 2025
      IFLAG=2
      ICODE=ICODE-10
 2025 CONTINUE
C
      GO TO (2100,2100,2300,2100,2500,2100),ICODE
C
C     ICODE=1  YUKAWA POTENTIAL
C     ICODE=2  COULOMB POTENTIAL
C     ICODE=3  OPEP TENSOR POTENTIAL
C     ICODE=4  GAUSSIAN POTENTIAL
C     ICODE=5  TWO PARTICLE TRANSFER
C     ICODE=6  ZERO RANGE KNOCK OUT
C     ICODE=11-16  GIVES IMAGINARY FORM FACTOR
C
 2100 LVR=0
      TOPT=1.D0
      GO TO 2350
 2300 ICODE=1
      LVR=2
      TOPT=-SQRTEN
 2350 MNNL=IABS(LTR1-LVR-1)+1
      MXXL=LTR1+LVR
      LTR2=LTR1+LTR1-2
      LVR2=LVR+LVR
      LLX=QNUM(2,1)*2.0
      LPX=QNUM(2,2)*2.0
      JLX=QNUM(3,1)
      JPX=QNUM(3,2)
      IS1=QNUM(4,1)*2.0
      IS2=QNUM(4,2)*2.0
      OPT=OPT*PHASEF(LLX/2)
      DRF=DRT
      DO 2390 LAM =MNNL,MXXL,2
      LAMA=LAM
      LAM2 =LAM +LAM -2
      VOPT=TOPT*OPT*SQR4PI*SQRT(DFLOAT(LAM2 +1)*DFLOAT(IQN(2)+1))
     1 *VCC(LAM2 ,LVR2,LTR2 ,0,0)
     2 *RACAH(IQN(2),IQN(2),LAM2 ,LTR2 ,LVR2,IQN(3))
      RME=0.D0
      IF(IQN(2).EQ.0) RME=1.D0
      IF(IQN(2).EQ.2) RME=SQRT(DFLOAT(IS1*(IS1+2)))
      RME=RME*SQRT(DFLOAT((IS1+1)*(IQN(2)+1)))
     1*SQRT(DFLOAT((LLX+1)*(LAM2+1)))*VCC(LLX,LAM2,LPX,0,0)
     2 *SQRT(DFLOAT(JLX+1)*DFLOAT(JPX+1)*DFLOAT(IQN(3)+1))
     3 *WINEJ(LPX,LLX,LAM2 ,IS2,IS1,IQN(2),JPX,JLX,IQN(3))
       IF(JI(1).NE.0) RME=RME*PHASEF((JI(1)+IQN(3)-JLX-JJ(2))/2)
     1 *RACAH(JPX,JJ(2),JLX,JJ(1),JI(1),IQN(3))
     2 *SQRT(DFLOAT(JJ(1)+1)*DFLOAT(JJ(2)+1))
      VOPT=VOPT*RME
      VOPT=VOPT*RAM
      TOTAL=TOTAL+ABS(VOPT)
      SL=0.D0
      RMS=0.D0
      IF(VOPT.EQ.0.0) GO TO 2385
      IF(ICODE.NE.6) GO TO 2365
      SCALE=VOPT/(SQR4PI*SQR4PI)
      R=0.D0
      MK=IFLAG
      DO 2360 M=1,KT
      R=R+DRF
      RSQ=R*R
      TEMP=UB(M    ,1)*UB(M+KMAX,1)*SCALE
      SI(MK  )=SI(MK  )+TEMP
      SL      =SL      +TEMP*RSQ
      RMS     =RMS     +TEMP*RSQ*RSQ
      MK=MK+2
 2360 CONTINUE
      IF(SL.NE.0.0) RMS=SQRT(RMS/SL)
      SL=SL*DRF
      GO TO 2385
C
 2365 CONTINUE
      CALL SLATR (KT,DRF,VB(  1),LTR1,FMUV,ICODE,KMAX)
      IF(LAM.NE.LTR1) GO TO 2375
      DO 2370 I=1,KT
      VB(I+KMAX2)=VB(I     )
      VB(I+KMAX3)=VB(I+KMAX)
 2370 CONTINUE
      GO TO 2380
 2375 CONTINUE
      CALL SLATR (KT,DRF,VB(KMAX2+1),LAMA,FMUV,ICODE,KMAX)
 2380 CONTINUE
      CALL RADIN(KT,DRF,VB,UB,SL,RMS,VOPT,SI(IFLAG),ICODE,LTR1,FMUV
     1  ,KMAX)
C
 2385 CONTINUE
      I=LAM -1
      WRITE(6,9102) IQN(1),I,RME,SL,RMS
 2390 CONTINUE
      GO TO 2990
C
C     HERE FOR TWO NUCLEON TRANSFER
C
 2500 CONTINUE
      M=0
      UB(M,1)=0.0
      IF(FMUV.EQ.0.0) FMUV=1.7
      CALL DSTRIP(IQN,DRF,KT,UB,SI(IFLAG),QNUM,OPT,FMUV,KM,SL)
      TOTAL=TOTAL+ABS(SL)
      WRITE(6,9002)KM,SL
      IQFLG=1
      GO TO 2990
 2990 CONTINUE
 3000 CONTINUE
      IF(CNTROL.GT.0.0) GO TO 100
      IF(IQFLG.EQ.1) KT=-KT
      K=KT
 3010 CONTINUE
      DRT=BKIN(15)
      IF(TOTAL.EQ.0.0) IERR=1
      RETURN
 9001 FORMAT(10F8.4)
 9002 FORMAT(15H FORM FACTOR,M=,I3,E18.6)
 9050 FORMAT(33H0ANGULAR MOMENTUM COUPLING  2*J1=,I3,3X,5H2*JC=,I3,3X,
     15H2*JI=,I3,5X,5H2*J2=,I3,3X,5H2*JC=,I3,3X,5H2*JF=,I3)
 9051 FORMAT(18H0PARAMETERS       ,9H   QCODE=,F9.4,9H   RANGE=,F9.4,
     17H   V0 =,F11.4,3X,21HSING. PART. OVERLAP =,F12.4,9H   RMS R=,
     2F9.4)
 9052 FORMAT(1H ,F6.2,10E12.4)
 9102 FORMAT(18X,9H   LTR  =,I4,5X,9H   LAM  =,I4,5X,9H   RME  =,F9.4,3X
     1,21HF.F. VOL. INTEGRAL  =,F12.4,9H   RMS R=,F9.4)
 9501 FORMAT(25H0FORM FACTOR DATA   PART.,I3,3H OF,F3.0)
 9998 FORMAT(1H1)
      END
      SUBROUTINE BCHK
      IMPLICIT REAL*8(A-H,O-Z)
C***********************************************************************
C
C     THIS ROUTINE MAY BE USED IN EITHER OF TWO MODES--
C       (1) A STANDARD MODE WITH FIXED DIMENSIONS AND FIXED CORE LOAD;
C       (2) A FLOATING MODE IN WHICH THE ROUTINE WILL ACQUIRE NEEDED
C           CORE DYNAMICALLY.  THE PARTICULAR CODE FOR THIS IS DESIGNED
C           FOR OPERATION AT THE UNIV. OF PITTSBURGH DEC-10 COMPUTER.
C
C     THE STANDARD-CODE AND FLOATING-CODE SECTIONS ARE MARKED AT
C     THEIR BEGINNINGS AND ENDS.
C
C     IN THE FLOATING-CODE FORM, THE ROUTINES THAT MUST BE REPLACED ON
C     OTHER SYSTEMS ARE--
C       CORSIZ(KUSE,KLIM) -- RETURNS CURRENT CORE USAGE AND UPPER LIMIT
C             IN UNITS OF PAGES (512 WORDS).
C       HEREIS(CORE) -- INITIALIZES DYNAMIC-ALLOCATION PACKAGE.
C       ALLOC (J,LEN,KRET) -- ACQUIRES CORE OF LENGTH LEN WORDS.
C             J IS THE LOCATION IN THE DUMMY ARRAY CORE AND KRET IS
C             USED FOR ERROR RETURNS (ERROR CODE IN J).
C       FREE  (CORE(J)) -- RELEASES CORE AT LOCATION J.
C
C***********************************************************************
C
C     FOR SINGLE PRECISION VERSIONS, CHANGE THE DEFINITION OF IDBL,
C     WHICH OCCURS AT THE BEGINNING OF THE ROUTINE, TO READ
C       IDBL=1
C
C***********************************************************************
C     DOUBLE PRECISION BETA
      DIMENSION LIM(4),NEED(6)
      COMMON THETAN,THETA1,DTHETA,DRT,RMAX,ALPHA(15),IDAT(6)
     1,ICON(20),LMAX,KT,KMAX,LLCO,JOUT,NLTR,IBUFF,IWORD
     2,KLAM( 8),LTRT( 8),JTRT( 8),NADD( 8),LAM(16),ISP(16),JSP(16)
     3,ILK(16),JPRT( 8)
      COMMON EN( 8),FM( 8),Z( 8),FMA( 8),ZA( 8),RY( 8),FS( 8),ECM( 8)
     1,FMU( 8),RC( 8),FK2( 8),FK( 8),ETA( 8),AFACT( 8),DRF( 8),FN( 8)
     2,PNLOC( 8),DUM( 8),SIG( 8)
     3,IS( 8),NS( 8),NLI( 8),LKTAB(32),ICUP(32),JCUP(32)
     4,IN( 8),INN( 8, 8),INDTAB(10)
C***********************************************************************
C        STANDARD CODE
C        ======== ====
C     DOUBLE PRECISION ZI,CZ,AMAT,CMAT
      DIMENSION W(800,8),V(400,8),ZI(80,2),CZ(80,5)
     1         ,AMAT(80,40),CMAT(80,40),PLM(6400),SR(2048),SI(2048)
      COMMON/STORAG/ CORE(12800)
      COMMON/POT/ X(800,16)
      EQUIVALENCE (W(1,1),CORE(1)), (V(1,1),CORE(6401))
     1           ,(ZI(1,1),CORE(9601)), (CZ(1,1),CORE(10401))
     2           ,(AMAT(1,1),CORE(1)),  (CMAT(1,1),CORE(6401))
     3           ,(PLM(1),CORE(1)),(SR(1),CORE(6401)),(SI(1),CORE(8451))
      DATA NCHAN/   40/
      DATA MAXD /12800/
      DATA MAXT / 2048/
      DATA MAXLM/ 6400/
C        END STANDARD CODE
C***********************************************************************
C        FLOATING CODE
C        ======== ====
C     LOGICAL FIRST
C     DIMENSION CORE(1)
C     DATA FIRST/.TRUE./
C     DATA NCHAN/  100/
C        END FLOATING CODE
C***********************************************************************
C
      CALL CUDATE(IDAT)
      WRITE(6,9006) (IDAT(KJ),KJ=4,6)
      KMAX2=KMAX+KMAX
      NPAR=NLI(1)
      IDBL=2
C
C     CHECK STORAGE SIZE FOR D AMPLITUDE SPACE
C
      LK=0
      NLL=0
      DO 40 I=1,NLTR
      M=(JTRT(I)+1)*NS(I)
      LK=LK+M
      NLL=MAX0(NLL,M)
   40 CONTINUE
      M=(JTRT(1)+1)*NS(1)
      LK=(LK+LK)*M*LMAX
      NLL=(NLL+NLL)*M
      LIM(2)=LK
      LIM(3)=NLL
C
C     COUNT MAXIMUM NUMBER OF SUBCHANNELS, YLM TERMS
C
      MAXM=0
      NLL=0
      DO 60 I=1,NLTR
      NLL=NLL+((JTRT(I)+1)*NS(I)+1)/2
      MAXM=MAX0(MAXM,(JTRT(1)+JTRT(I)+IS(1)+IS(I))/2+1)
   60 CONTINUE
      LIM(1)=NLL
      LIM(4)=MAXM*LMAX
      NLL2=NLL+NLL
C
C***********************************************************************
C        STANDARD CODE
C        ======== ====
C
C     READ IN DIAGONAL AND COUPLING POTENTIALS
C
      DO 120 I=1,NLTR
      READ(9) (W(J,I),J=1,KMAX2),(V(J,I),J=1,KMAX)
  120 CONTINUE
      DO 130 I=1,NPAR
      READ(8) (X(J,I),J=1,KMAX2)
  130 CONTINUE
C
      WRITE(6,9100) LIM,NCHAN,MAXD,MAXT,MAXLM
      IF(LIM(1).GT.NCHAN) GO TO 1200
      IF(LIM(2).GT.MAXD ) GO TO 1200
      IF(LIM(3).GT.MAXT ) GO TO 1200
      IF(LIM(4).GT.MAXLM) GO TO 1200
C
      CALL CUDATE(IDAT)
      WRITE(6,9007) (IDAT(KJ),KJ=4,6)
C
      CALL INTEG(ZI,CZ,W,V,X,NLL2)
C
      CALL CUDATE(IDAT)
      WRITE(6,9008) (IDAT(KJ),KJ=4,6)
C
      CALL SOLVE(AMAT,CMAT,X,NLL2)
C
      CALL CUDATE(IDAT)
      WRITE(6,9009) (IDAT(KJ),KJ=4,6)
      WRITE(6,9019)ALPHA,IDAT
C
      CALL XSECT(X,PLM,SR,SI)
C        END STANDARD CODE
C***********************************************************************
C        FLOATING CODE
C        ======== ====
C     WRITE(6,9100) LIM,NCHAN
C     IF(LIM(1).GT.NCHAN) GO TO 1200
C     IF(FIRST) CALL CORSIZ(KORUSE,KORLIM)
C     FIRST=.FALSE.
C     CALL HEREIS(CORE)
C     NEED(1)=KMAX2*NLTR
C     NEED(2)=KMAX *NLTR
C     NEED(3)=KMAX2*NPAR
C     ISUM=NEED(1)+NEED(2)+NEED(3)
C     NEED(4)=7*IDBL*NLL2
C     NEED(5)=2*IDBL*NLL2*NLL+LIM(2)
C     NEED(6)=LIM(4)+2*LIM(3)
C     ISUM=MAX0(NEED(4)+ISUM,NEED(5),NEED(6)+LIM(2))
C     NPG=(ISUM+511)/512
C     KORNED=KORUSE+NPG+4
C     KORSAF=KORNED+1
C     WRITE(6,9112) KORSAF,KORLIM
C     IF(KORSAF.GT.KORLIM) GO TO 1000
C     CALL ALLOC(LS,2048)
C     CALL ALLOC(KOR1,ISUM)
C     CALL FREE(CORE(LS))
C     LW=KOR1
C     LV=LW+NEED(1)
C     LX=LV+NEED(2)
C     LZ=LX+NEED(3)
C     LC=LZ+ 2*IDBL*NLL2
C
C     READ IN DIAGONAL AND COUPLING POTENTIALS
C
C     KW1=LW
C     KV1=LV
C     DO 120 I=1,NLTR
C     KW2=KW1+799
C     KV2=KV1+399
C     READ(9) (CORE(J),J=KW1,KW2),(CORE(J),J=KV1,KV2)
C     KW1=KW1+800
C     KV1=KV1+400
C 120 CONTINUE
C     KX1=LX
C     DO 130 I=1,NPAR
C     KX2=KX1+799
C     READ(8) (CORE(J),J=KX1,KX2)
C     KX1=KX1+800
C 130 CONTINUE
C
C     CALL CUDATE(IDAT)
C     WRITE(6,9007) (IDAT(KJ),KJ=4,6)
C
C     CALL INTEG(CORE(LZ),CORE(LC),CORE(LW),CORE(LV),CORE(LX),NLL2)
C     LD=KOR1
C     LA=LD+LIM(2)
C     LC=LA+IDBL*NLL2*NLL
C
C     CALL CUDATE(IDAT)
C     WRITE(6,9008) (IDAT(KJ),KJ=4,6)
C
C     CALL SOLVE(CORE(LA),CORE(LC),CORE(LD),NLL2)
C     LW=LA+LIM(4)
C     LV=LW+LIM(3)
C
C     CALL (IDAT)
C     WRITE(6,9009) (IDAT(KJ),KJ=4,6)
C     WRITE(6,9019)ALPHA,IDAT
C
C     CALL XSECT(CORE(LD),CORE(LA),CORE(LW),CORE(LV))
C     CALL FREE(CORE(KOR1))
C        END FLOATING CODE
C***********************************************************************
C
 1000 CONTINUE
      RETURN
 1200 CONTINUE
      WRITE(6,9110)
      GO TO 1000
 9006 FORMAT(24H0SUBROUTINE BCHK ENTERED,10X,I4,2(1H.,I2))
 9007 FORMAT(24H0SUBROUTINE INTEG CALLED,10X,I4,2(1H.,I2))
 9008 FORMAT(24H0SUBROUTINE SOLVE CALLED,10X,I4,2(1H.,I2))
 9009 FORMAT(24H0SUBROUTINE XSECT CALLED,10X,I4,2(1H.,I2))
 9019 FORMAT(1H1,17X,15A4,4X,I2,2(1H/,I2),I4,2(1H.,I2))
 9100 FORMAT(22H0ARRAY NEEDS -- NCHAN=,I4,3X,7HNDAMPL=,I6,3X,5HNTR1=,I5,
     13X,5HNYLM=,I5,3X,12HWITH LIMITS ,I4,1H,,I6,1H,,I5,1H,,I5)
 9110 FORMAT(16H0LIMITS EXCEEDED)
 9112 FORMAT(15H0PROGRAM NEEDS ,I3,28H CORE PAGES.  YOUR LIMIT IS ,I3)
      END
      SUBROUTINE INTEG(ZI,CZ,W,V,X,NLL2)
      IMPLICIT REAL*8(A-H,O-Z)
C     SQRT(XXX)=DSQRT(XXX)
C     AMAX1(XXX,YYY)=DMAX1(XXX,YYY)
C     AMIN1(XXX,YYY)=DMIN1(XXX,YYY)
C     FLOAT(III)=DFLOAT(III)
C     DOUBLE PRECISION BETA,CZ,ZI
C modified by (JE,ME)
C     USE THE NEXT 2 CARDS FOR 32-BIT MACHINES
      DOUBLE PRECISION VCC,RACAH,WINEJ
      DOUBLE PRECISION CZ,ZI
      DIMENSION ZI(NLL2,2),CZ(NLL2,5),W(1),V(1),X(1)
      COMMON THETAN,THETA1,DTHETA,DRT,RMAX,ALPHA(15),IDAT(6)
     1,ICON(20),LMAX,KT,KMAX,LLCO,JOUT,NLTR,IBUFF,IWORD
     2,KLAM( 8),LTRT( 8),JTRT( 8),NADD( 8),LAM(16),ISP(16),JSP(16)
     3,ILK(16),JPRT( 8)
      COMMON EN( 8),FM( 8),Z( 8),FMA( 8),ZA( 8),RY( 8),FS( 8),ECM( 8)
     1,FMU( 8),RC( 8),FK2( 8),FK( 8),ETA( 8),AFACT( 8),DRF( 8),FN( 8)
     2,PNLOC( 8),DUM( 8),SIG( 8)
     3,IS( 8),NS( 8),NLI( 8),LKTAB(32),ICUP(32),JCUP(32)
     4,IN( 8),INN( 8, 8),INDTAB(10)
      COMMON/LOCAL/ FTAB(100),VTAB(100),NTAB(100),LTAB(100),JTAB(100)
     1      ,MTAB(100),NTBL(100),MI(100),IFAKE1(600),CTAB(510)
     2      ,IKTAB(510),JKTAB(510),LKTLK(510),IKTBL(510),MZTBL(512)
     3      ,ITEMP(8,8),IFAKE2(2616),DR2(8),KTT,MZ,NMM,IK,BUFC(464)
      COMMON/BUFSPC/ BUFFR(800)
      DATA MAXIJ/510/
      DATA ETA1/1.0D-08/
      DATA SQRFPI/3.5449077/
      SAVE
C
      IPWR=2*ICON(20)
      IF(IPWR.EQ.0) IPWR=8
      ETA1=10.0**(-IPWR)
      REWIND 8
      KTT=KT
      NCUP=NLI(2)
      IWORD=0
      NPAR=2
      LDEL=1
      IF(LTRT(1).NE.0) GO TO 21
      IF(  IS(1).NE.0) GO TO 21
      NPAR=1
      LDEL=2
   21 CONTINUE
      DO 40 I=1,NLTR
      DR2(I)=DRF(I)**2
   40 CONTINUE
C
C     CONSTRUCT EFFECTIVE COUPLING MATRIX
C
      DO 52 I=1,NLTR
      INN(I,I)=INN(I,I)+1
   52 CONTINUE
      DO 56 I=1,NLTR
      DO 54 J=1,NLTR
      DO 54 K=1,NLTR
      MK=0
      DO 53 M=1,NLTR
      MK=MK+INN(J,M)*INN(M,K)
   53 CONTINUE
      IF(MK.GT.1) MK=2
      ITEMP(J,K)=MK
   54 CONTINUE
      DO 55 J=1,NLTR
      DO 55 K=1,NLTR
      INN(J,K)=ITEMP(J,K)
   55 CONTINUE
   56 CONTINUE
      DO 57 I=1,NLTR
      INN(I,I)=INN(I,I)-1
   57 CONTINUE
C
C
      JX1=(1.0-PHASEF(JTRT(1)+IS(1)))/2.0
      JMAX=LMAX
      LMAX22=LMAX+LMAX-2
C
C      J LOOP
C
      DO 1000 JA=1,JMAX
      JT=JA+JA-2+JX1
      LMIN=IABS(IABS(JT-JTRT(1))-IS(1))/2
      DO 1000 JPAR=1,NPAR
C
C     SELECT SUBCHANNELS
C
      NLL=0
      MAXL=0
      DO 80 I=1,NLTR
      NLI(I)=0
      JMXX=     JT+JTRT(I)
      JMNN=IABS(JT-JTRT(I))
      IF(JMXX.LT.JMNN) GO TO 80
      XPHAS=PHASEF(JPRT(1))*PHASEF(JPRT(I))
     1*DFLOAT(ISIGN(1,LTRT(1))*ISIGN(1,LTRT(I)))
      IN(I)=SIGN(1.1D0,XPHAS)
      DO 70 J=JMNN,JMXX,2
      LMXX=MIN0(J+IS(I),LMAX22)/2
      LMNN=IABS(J-IS(I))/2
      IF(LMXX.LT.LMNN) GO TO 70
      DO 60 L=LMNN,LMXX
      PHAS=PHASEF(LMIN+JPAR-1-L)*XPHAS
      IF(PHAS.LT.0.0) GO TO 60
      NLL   =NLL   +1
      NLI(I)=NLI(I)+1
      LTAB(NLL)=L
      JTAB(NLL)=J
      NTAB(NLL)=I
      FTAB(NLL)=L*(L+1)
      VTAB(NLL)=FLOAT(J*(J+2)-IS(I)*(IS(I)+2))/4.0-FTAB(NLL)
      VTAB(NLL)=VTAB(NLL)*0.5
      MTAB(NLL)=(I-1)*KMAX
   60 CONTINUE
   70 CONTINUE
   80 CONTINUE
C modifiedy by JE,ME
      WRITE(8) NLL,NLI
     1        ,(NTAB(L),L=1,NLL),(LTAB(L),L=1,NLL),(JTAB(L),L=1,NLL)
      IF(NLI(1).EQ.0) GO TO 980
      MAXA=NLL+NLL
C
C      COEFFICENT TABLE
C
      IJ=0
      IF(NCUP.EQ.0) GO TO 115
C     I IS INITIAL CHANNEL
      DO 110 I=1,NLL
      LX=LTAB(I)*2
      JX=JTAB(I)
      II=NTAB(I)
      IF((II.NE.1).AND.(LTAB(I).LT.LLCO)) GO TO 110
      IX=JTRT(II)
      ISX=IS(II)
C     J IS FINAL CHANNEL
      DO 105 J=1,NLL
      LY=LTAB(J)*2
      JY=JTAB(J)
      JJ=NTAB(J)
      IF((JJ.NE.1).AND.(LTAB(J).LT.LLCO)) GO TO 105
      IY=JTRT(JJ)
      ISY=IS(JJ)
      DO 100 K=1,NCUP
      IF((ICUP(K).NE.II).OR.(JCUP(K).NE.JJ)) GO TO 100
      LK=LKTAB(K)
      LL2=LAM(LK)*2
      JJ2=JSP(LK)
      IS2=ISP(LK)
      TEMP=RACAH(JX,IX,JY,IY,JT,JJ2)
     1*SQRT(DFLOAT(JX+1)*DFLOAT(JY+1)*DFLOAT(JJ2+1))
     2*WINEJ(LX,ISX,JX,LL2,IS2,JJ2,LY,ISY,JY)
     3*VCC(LX,LL2,LY,0,0)*SQRT(DFLOAT(LX+1)*DFLOAT(LL2+1))
     4*PHASEF((JT-IX-JY-JJ2)/2)/SQRFPI
     5*PHASEF((LX-LL2-LY)/4)
      IF(ILK(LK).EQ.0) GO TO 82
      IF(ILK(LK).EQ.1) GO TO 84
      IF(ILK(LK).EQ.2) GO TO 86
      IF(ILK(LK).EQ.3) GO TO 88
      IF(ILK(LK).EQ.4) GO TO 84
      IF(ILK(LK).EQ.7) GO TO 90
      IF(ILK(LK).EQ.8) GO TO 90
      GO TO 92
C
C     COLLECTIVE
C
   82 CONTINUE
      KLMM=KLAM(II)
      TEMP=TEMP*SQRT(DFLOAT(IX+1)*DFLOAT(ISY+1))*VCC(IX,LL2,IY,KLMM,0)
      IPHAS=-LL2
      IF((DFLOAT(ISIGN(1,LTRT(II)))*PHASEF(JPRT(II))).LT.0.0)
     1  IPHAS=IPHAS+2
      IF((DFLOAT(ISIGN(1,LTRT(JJ)))*PHASEF(JPRT(JJ))).LT.0.0)
     1  IPHAS=IPHAS-2
      TEMP=TEMP*PHASEF(IPHAS/4)
      GO TO 92
   84 CONTINUE
C
C     PARTICLE TRANSFER
C
C     PICKUP
      IF(FMA(II).GT.FMA(JJ)+0.2) TEMP=TEMP
     1*SQRT(DFLOAT(IX+1)*DFLOAT(ISY+1))
     2*PHASEF((IX+JJ2-IY)/2)
C     STRIPPING
      IF(FMA(JJ).GT.FMA(II)+0.2) TEMP=TEMP
     1*SQRT(DFLOAT(IY+1)*DFLOAT(ISX+1))
     2*PHASEF((ISY+IS2-ISX)/2)
      GO TO 92
   86 CONTINUE
C
C     MICROSCOPIC
C
      TEMP=TEMP*PHASEF(JJ2/2)
      IF(IS2.EQ.0) TEMP=TEMP*SQRT(DFLOAT(ISY+1))
      IF(IS2.EQ.2)TEMP=TEMP*(SQRT(DFLOAT(ISX)*DFLOAT(ISY)*2.)
     1 -DFLOAT(IABS(ISX-ISY))*SQRT(3.D0))
      IF(II.LE.JJ) GO TO 92
      IPHAS=IY+JJ2-IX+ISY+IS2-ISX
      TEMP=TEMP*PHASEF(IPHAS/2)
      GO TO 92
C
C     COLLECTIVE-NON STRONG COUPLING MODEL *** K IS NOT A GOOD QUANT. NO
C
   88 CONTINUE
      TEMP=TEMP*SQRT(DFLOAT((IY+1)*(ISY+1))/DFLOAT(LL2+1))
      IPHAS=-LL2
      IF((DFLOAT(ISIGN(1,LTRT(II)))*PHASEF(JPRT(II))).LT.0.0)
     1  IPHAS=IPHAS+2
      IF((DFLOAT(ISIGN(1,LTRT(JJ)))*PHASEF(JPRT(JJ))).LT.0.0)
     1  IPHAS=IPHAS-2
      TEMP=TEMP*PHASEF(IPHAS/4)
      IF(II.LT.JJ) GO TO 92
      TEMP=TEMP*SQRT(DFLOAT(IX+1)/DFLOAT(IY+1))*PHASEF((IX-IY)/2)
   90 CONTINUE
C
C    GAMMA BAND COUPLING.OR.BETA BAND COUPLING
C    OCTUPOLE BAND COUPLING AS WELL
C
      KLMM=KLAM(II)
      MU=IABS(KLAM(II)-KLAM(JJ))
      TEMP=TEMP*SQRT(DFLOAT((IX+1)*(ISY+1)))*VCC(IX,LL2,IY,-KLMM,MU)
      IPHAS=-LL2
      IF((DFLOAT(ISIGN(1,LTRT(II)))*PHASEF(JPRT(II))).LT.0.0)
     1  IPHAS=IPHAS+2
      IF((DFLOAT(ISIGN(1,LTRT(JJ)))*PHASEF(JPRT(JJ))).LT.0.0)
     1  IPHAS=IPHAS-2
      TEMP=TEMP*PHASEF(IPHAS/4)
C   THIS PHASE IS TO TAKE CARE OF THE DEEXCITATION OF THE EXCITED BAND
C   TO THE GS BAND. IF THE EXCITED BAND HAS K=0 THIS PHASE IS ALWAYS
C   POSITIVE
      IF(KLAM(II).EQ.0) TEMP=TEMP*PHASEF((IX+LL2-IY)/2)
   92 CONTINUE
      T1=MAX(DRF(II),DRF(JJ))
      T2=MIN(DRF(II),DRF(JJ))
      TEMP=TEMP*SQRT(T1/T2)**3
      IF(TEMP.EQ.0.0) GO TO 100
      IJ=IJ+1
      JKTAB(IJ)=I
      IKTAB(IJ)=J
      LKTLK(IJ)=(LK-1)*KMAX*2
      CTAB(IJ) = -TEMP*FMU(JJ)*DRF(JJ)**3/DRF(II)
  100 CONTINUE
  105 CONTINUE
  110 CONTINUE
  115 CONTINUE
      IF(IJ.LE.MAXIJ) GO TO 120
      WRITE(6,9101) IJ,MAXIJ
      STOP 111
  120 CONTINUE
C
C     ITERATION LOOP
C
      DO 800 N1=1,NLL
      NJ=N1+N1-1
      II=NTAB(N1)
C
C      SET UP INITIAL Z  VALUES
C
      NMM=0
      DO 180 N2 = 1,NLL
      NK=N2+N2-1
      JJ=NTAB(N2)
      DO 165 M=1,5
      CZ(NK  ,M)=0.D0
      CZ(NK+1,M)=0.D0
      IF(M.GT.2) GO TO 165
      ZI(NK  ,M)=0.D0
      ZI(NK+1,M)=0.D0
  165 CONTINUE
      IF(N1.EQ.N2) GO TO 170
      IF(INN(JJ,II).EQ.0) GO TO 180
  170 CONTINUE
      NMM=NMM+1
      NTBL(NMM)=N2
      FL=LTAB(N2)-ICON(19)
      FL=MAX(FL,0.0D0)
      M1=SQRT((FL*FL+FL)/12.0)
      EXPO=1.0/(2.0*FL+1.0)
      M2=((1.0D-16)**EXPO)*(2.0*FL/(2.718*FK(JJ)*DRF(JJ)))
      M1=MAX0(M1,M2)+1
      MI(NMM)=MIN0(M1,KMAX)
      IF(N1.NE.N2) GO TO 180
      MZ=MI(NMM)
      ZI(NK  ,2)=ETA1
      CZ(NK  ,5)=ETA1
  180 CONTINUE
C
C     CONSTRUCT SECONDARY TABLES
C
      IK=0
      IF(IJ.EQ.0) GO TO 196
      DO 195 I=1,NMM
      DO 190 J=1,NMM
      DO 185 K=1,IJ
      IF(IKTAB(K).NE.NTBL(I)) GO TO 185
      IF(JKTAB(K).NE.NTBL(J)) GO TO 185
      IK=IK+1
      IKTBL(IK)=K
      MZTBL(IK)=0
      IF(MI(I).LE.MI(J)) GO TO 185
      MZTBL(IK)=MI(I)*262144+K*512+IK
      IKTBL(IK)=0
  185 CONTINUE
  190 CONTINUE
  195 CONTINUE
  196 CONTINUE
      MZTBL(IK+1)=(KMAX+1)*262144
C
C     SORT STARTING POINTS
C
      IF(IK.EQ.0) GO TO 221
      DO 220 J=1,IK
      M1=KMAX
      DO 210 K=J,IK
      I=MZTBL(K)/262144
      IF(I.GT.M1) GO TO 210
      M1=I
      K1=K
  210 CONTINUE
      N2=MZTBL(K1)
      IF(K1.EQ.J) GO TO 206
      K1=K1-1
      DO 205 K=J,K1
      M2=K1-K+J
      MZTBL(M2+1)=MZTBL(M2)
  205 CONTINUE
  206 CONTINUE
      MZTBL(J)=N2
  220 CONTINUE
  221 CONTINUE
C
C      INTEGRATION SECTION
C
      IF(ICON(18).EQ.0) GO TO 240
C
      CALL NUMERV(ZI,CZ,W,V,X,NLL2)
      GO TO 250
C
  240 CONTINUE
      CALL STORMR(ZI,CZ,W,V,X,NLL2)
C
  250 CONTINUE
C
C     WRITE LAST TWO POINTS ON TAPE
C

C modified by JE,ME
      WRITE(8) (ZI(I,1),I=1,MAXA)
      WRITE(8) (ZI(I,2),I=1,MAXA)
  800 CONTINUE
C
C     END OF ITERATION LOOP
C
  980 CONTINUE
 1000 CONTINUE
C
C     FLUSH BUFFR
C
C modified by JE,ME
      WRITE(8) NLL
      IF(IWORD.EQ.IBUFF) GO TO 1001
      IWORD=IBUFF-1

 1001 CONTINUE
      REWIND 8
      RETURN
 9101 FORMAT(28H0TOO MANY COUPLING TERMS -- ,I4,13H  WITH LIMIT ,I4)
      END
      SUBROUTINE STORMR(ZI,CZ,W,V,X,NLL2)
      IMPLICIT REAL*8(A-H,O-Z)
C     FLOAT(III)=DFLOAT(III)
      LOGICAL COFSET
C     DOUBLE PRECISION CZ,ZI
C     DOUBLE PRECISION CII1,CII2,CF1,CF2,CF3,CF4,CF5
      DIMENSION ZI(NLL2,2),CZ(NLL2,5),W(1),V(1),X(1)
      COMMON/LOCAL/ FTAB(100),VTAB(100),NTAB(100),LTAB(100),JTAB(100)
     1      ,MTAB(100),NTBL(100),MI(100),IFAKE1(600),CTAB(510)
     2      ,IKTAB(510),JKTAB(510),LKTLK(510),IKTBL(510),MZTBL(512)
     3      ,ITEMP(8,8),IFAKE2(2616),DR2(8),KT,MZ,NMM,IK,BUFC(464)
      DATA MSW,COFSET/6,.FALSE./
C
      IF(COFSET) GO TO 10
      CF1=  19.D0/240.D0
      CF2=- 96.D0/240.D0
      CF3= 194.D0/240.D0
      CF4=-176.D0/240.D0
      CF5= 299.D0/240.D0
      COFSET=.TRUE.
   10 CONTINUE
      M1=1
      MX=1
      DO 600 M=1,KT
      MK=M+M-1
  320 CONTINUE
      IF(M.LT.MZTBL(M1)/262144) GO TO 321
      N2=MZTBL(M1)-262144*(MZTBL(M1)/262144)
      K1=N2-512*(N2/512)
      IF(K1.NE.0) IKTBL(K1)=N2/512
      M1=M1+1
      GO TO 320
  321 CONTINUE
      IF(M.EQ.MZ) GO TO 530
      IF(M.LT.MZ) GO TO 600
      IF(MX.GE.MSW) GO TO 500
      MX=MX+1
C
C      INTEGRATE U(R) BY N=2 STORMER EQUATION
C
      DO 410 NY=1,NMM
      N2=NTBL(NY)
      NK=N2+N2-1
C
      CII1  =2.0*ZI(NK  ,2)-ZI(NK  ,1)+CZ(NK  ,5)
      CII2  =2.0*ZI(NK+1,2)-ZI(NK+1,1)+CZ(NK+1,5)
C
C     STEP INTEGRATION TABLES BACK ONE
C
      DO 408 K=1,4
      CZ(NK  ,K)=CZ(NK  ,K+1)
      CZ(NK+1,K)=CZ(NK+1,K+1)
  408 CONTINUE
      ZI(NK  ,1)=ZI(NK  ,2)
      ZI(NK+1,1)=ZI(NK+1,2)
      ZI(NK  ,2)=CII1
      ZI(NK+1,2)=CII2
  410 CONTINUE
      GO TO 530
C
  500 CONTINUE
C
C     INTEGRATE U(R) BY N=5 STORMER EQUATION
C
      DO 510 NY=1,NMM
      N2=NTBL(NY)
      NK=N2+N2-1
      CII1  =2.0*ZI(NK  ,2)-ZI(NK  ,1)+CF1  *CZ(NK  ,1)+CF2  *CZ(NK  ,2)
     1+CF3  *CZ(NK  ,3)+CF4  *CZ(NK  ,4)+CF5  *CZ(NK  ,5)
      CII2  =2.0*ZI(NK+1,2)-ZI(NK+1,1)+CF1  *CZ(NK+1,1)+CF2  *CZ(NK+1,2)
     1+CF3  *CZ(NK+1,3)+CF4  *CZ(NK+1,4)+CF5  *CZ(NK+1,5)
C
C     STEP INTEGRATION TABLES BACK ONE
C
      DO 508 K=1,4
      CZ(NK  ,K)=CZ(NK  ,K+1)
      CZ(NK+1,K)=CZ(NK+1,K+1)
  508 CONTINUE
      ZI(NK  ,1)=ZI(NK  ,2)
      ZI(NK+1,1)=ZI(NK+1,2)
      ZI(NK  ,2)=CII1
      ZI(NK+1,2)=CII2
  510 CONTINUE
  530 CONTINUE
C
C      COMPUTE NEW CZ VALUE
C
C
C     COMPUTE DIAGONAL PART OF CZ
C
      DO 540 NY=1,NMM
      N2=NTBL(NY)
      NK=N2+N2-1
      JJ=NTAB(N2)
      NL=MTAB(N2)+M
      NM=NL+NL-1
      CII1  =FTAB(N2)/FLOAT(M*M)-DR2(JJ)*(W(NM  )+V(NL)*VTAB(N2))
      CII2  =                   -DR2(JJ)* W(NM+1)
      CZ(NK  ,5)=CII1  *ZI(NK  ,2)-CII2  *ZI(NK+1,2)
      CZ(NK+1,5)=CII1  *ZI(NK+1,2)+CII2  *ZI(NK  ,2)
  540 CONTINUE
C
C     COMPUTE NON-DIAGONAL PART OF CZ
C
      IF(IK.EQ.0) GO TO 591
      DO 590 IM=1,IK
      I=IKTBL(IM)
      IF(I.EQ.0) GO TO 590
      NK=IKTAB(I)+IKTAB(I)-1
      NL=JKTAB(I)+JKTAB(I)-1
      LK=LKTLK(I)+MK
      CII1  = X(LK  )*CTAB(I)
      CII2  = X(LK+1)*CTAB(I)
      CZ(NK  ,5)=CZ(NK  ,5)+ CII1  *ZI(NL  ,2)-CII2  *ZI(NL+1,2)
      CZ(NK+1,5)=CZ(NK+1,5)+ CII1  *ZI(NL+1,2)+CII2  *ZI(NL  ,2)
  590 CONTINUE
  591 CONTINUE
  600 CONTINUE
      RETURN
      END
      SUBROUTINE NUMERV(ZI,CZ,W,V,X,NLL2)
      IMPLICIT REAL*8(A-H,O-Z)
C     INTEGRATION BY THE NUMEROV ALGORITHM
C     CF. A.C.ALLISON, J.COMP.PHYS. 6, 378 (1970)
C     FLOAT(III)=DFLOAT(III)
C     ABS(XXX)=DABS(XXX)
C     DOUBLE PRECISION CZ,ZI
C     DOUBLE PRECISION CII1,CII2,QR,QI,DET
      DIMENSION ZI(NLL2,2),CZ(NLL2,5),W(1),V(1),X(1)
      COMMON/LOCAL/ FTAB(100),VTAB(100),NTAB(100),LTAB(100),JTAB(100)
     1      ,MTAB(100),NTBL(100),MI(100),IFAKE1(600),CTAB(510)
     2      ,IKTAB(510),JKTAB(510),LKTLK(510),IKTBL(510),MZTBL(512)
     3      ,ITEMP(8,8),IFAKE2(2616),DR2(8),KT,MZ,NMM,IK,BUFC(464)
C
      M1=1
      DO 700 M=1,KT
      MK=M+M-1
  320 CONTINUE
      IF(M.LT.MZTBL(M1)/262144) GO TO 321
      N2=MZTBL(M1)-262144*(MZTBL(M1)/262144)
      K1=N2-512*(N2/512)
      IF(K1.NE.0) IKTBL(K1)=N2/512
      M1=M1+1
      GO TO 320
  321 CONTINUE
      IF(M.LT.MZ) GO TO 700
      IF(M.EQ.MZ) GO TO 530
C
C     INTEGRATE NUMEROV AUXILLIARY FUNCTION
C
      DO 410 NY=1,NMM
      N2=NTBL(NY)
      NK=N2+N2-1
C
      CII1  =12.0*ZI(NK  ,2)-10.0*CZ(NK  ,5)-CZ(NK  ,4)
      CII2  =12.0*ZI(NK+1,2)-10.0*CZ(NK+1,5)-CZ(NK+1,4)
C
C     STEP INTEGRATION TABLES BACK
C
      ZI(NK  ,1)=ZI(NK  ,2)
      ZI(NK+1,1)=ZI(NK+1,2)
      CZ(NK  ,4)=CZ(NK  ,5)
      CZ(NK+1,4)=CZ(NK+1,5)
      CZ(NK  ,5)=CII1
      CZ(NK+1,5)=CII2
  410 CONTINUE
  530 CONTINUE
C
C     COMPUTE NEW ZI VALUE
C
C
C     COMPUTE DIAGONAL ELEMENTS
C
      DO 540 NY=1,NMM
      N2=NTBL(NY)
      NK=N2+N2-1
      JJ=NTAB(N2)
      NL=MTAB(N2)+M
      NM=NL+NL-1
      CII1  =FTAB(N2)/FLOAT(M*M)-DR2(JJ)*(W(NM  )+V(NL)*VTAB(N2))
      CII2  =                   -DR2(JJ)* W(NM+1)
      QR=1.0-CII1/12.0
      QI=   -CII2/12.0
      DET=QR*QR+QI*QI
      QR= QR/DET
      QI=-QI/DET
      CZ(NK  ,1)=QR
      CZ(NK+1,1)=QI
      CZ(NK  ,3)=0.D0
      CZ(NK+1,3)=0.D0
      ZI(NK  ,2)=QR*CZ(NK  ,5)-QI*CZ(NK+1,5)
      ZI(NK+1,2)=QR*CZ(NK+1,5)+QI*CZ(NK  ,5)
  540 CONTINUE
C
C     COMPUTE NON-DIAGONAL COUPLINGS
C
      IF(IK.EQ.0) GO TO 600
      DO 590 IM=1,IK
      I=IKTBL(IM)
      IF(I.EQ.0) GO TO 590
      NK=IKTAB(I)+IKTAB(I)-1
      NL=JKTAB(I)+JKTAB(I)-1
      LK=LKTLK(I)+MK
      CII1  = X(LK  )*CTAB(I)
      CII2  = X(LK+1)*CTAB(I)
      CZ(NK  ,3)=CZ(NK  ,3)+(CII1  *ZI(NL  ,2)-CII2  *ZI(NL+1,2))
      CZ(NK+1,3)=CZ(NK+1,3)+(CII1  *ZI(NL+1,2)+CII2  *ZI(NL  ,2))
  590 CONTINUE
  600 CONTINUE
      DO 620 NY=1,NMM
      N2=NTBL(NY)
      NK=N2+N2-1
      CII1  =CZ(NK  ,5)+CZ(NK  ,3)/12.0
      CII2  =CZ(NK+1,5)+CZ(NK+1,3)/12.0
      ZI(NK  ,2)=CII1*CZ(NK  ,1)-CII2*CZ(NK+1,1)
      ZI(NK+1,2)=CII1*CZ(NK+1,1)+CII2*CZ(NK  ,1)
  620 CONTINUE
  700 CONTINUE
      RETURN
      END
      SUBROUTINE SOLVE(AMAT,CMAT,D,NLL2)
      IMPLICIT REAL*8(A-H,O-Z)
C     SQRT(XXX)=DSQRT(XXX)
C     SIN(XXX)=DSIN(XXX)
C     COS(XXX)=DCOS(XXX)
C     FLOAT(III)=DFLOAT(III)
C     DOUBLE PRECISION BETA,AMAT,CMAT,TEMP,D1,D2,D3,D4,D5,D6,CII,ZI
C modified by (JE,ME)
C     USE THE NEXT 3 CARDS FOR 32-BIT MACHINES
      DOUBLE PRECISION X1,X2,ETAC,ZERO,FLM,F,FP,G,GP,S,VCC,YXFCT
      DOUBLE PRECISION ARG,ATAB,BTAB,FTAB,FACT,CMR,CMI,COEFF
      DOUBLE PRECISION AMAT,CMAT,TEMP,D1,D2,D3,D4,D5,D6,CII,ZI
      DIMENSION AMAT(NLL2,1),CMAT(NLL2,1),D(1)
      COMMON THETAN,THETA1,DTHETA,DRT,RMAX,ALPHA(15),IDAT(6)
     1,ICON(20),LMAX,KT,KMAX,LLCO,JOUT,NLTR,IBUFF,IWORD
     2,KLAM( 8),LTRT( 8),JTRT( 8),NADD( 8),LAM(16),ISP(16),JSP(16)
     3,ILK(16),JPRT( 8)
      COMMON EN( 8),FM( 8),Z( 8),FMA( 8),ZA( 8),RY( 8),FS( 8),ECM( 8)
     1,FMU( 8),RC( 8),FK2( 8),FK( 8),ETA( 8),AFACT( 8),DRF( 8),FN( 8)
     2,PNLOC( 8),DUM( 8),SIG( 8)
     3,IS( 8),NS( 8),NLI( 8),LKTAB(32),ICUP(32),JCUP(32)
     4,IN( 8),INN( 8, 8),INDTAB(10)
      COMMON/LOCAL/ F(600),FP(600),G(600),GP(600),S(600),ZI(200)
     1      ,ATAB(100),BTAB(100),FTAB(100),E(200),NTAB(100),LTAB(100)
     2      ,JTAB(100),IFAKE(300),CII(2),MAXTAB(8),BLUB(394)
      COMMON/BUFSPC/ BUFFR(800)
C modified by (JE,ME)
      SAVE
C
      IWORD=IBUFF
      DELR=0.1D0
      JMAX=LMAX
      LK=0
      IX=JTRT(1)
      MAX1=JTRT(1)+1
      MXS=MAX1*NS(1)
      NPAR=2
      IF(MXS.EQ.1) NPAR=1
      INDTAB(1)=0
      MAXL=MAX0(LMAX,3)
      IND=1
      DO 1100 I=1,NLTR
      SIG(I)=0.D0
      X2=DRF(I)*FLOAT(KT)*FK(I)
      X1=X2-DRF(I)*FK(I)
      IF(FK2(I).LT.0.0) GO TO 1050
      CALL COU(X1,X2,ETA(I),MAXL,DELR,F(IND),FP(IND),G(IND),GP(IND)
     1,S(IND))
      GO TO 1051
 1050 CONTINUE
      CALL HAN(X1,X2,ETA(I),MAXL,DELR,F(IND),FP(IND),G(IND),GP(IND)
     1,S(IND))
 1051 CONTINUE
      IND=IND+LMAX
      LKTEMP=(JTRT(I)+1)*NS(I)
      LK=LK+LKTEMP*MXS*LMAX
      MAXTAB(I)=LKTEMP
      INDTAB(I+1)=LK
 1100 CONTINUE
      MAX=LK+LK
C
C     CLEAR AMPLITUDE STORAGE
C
      DO 1110 I=1,MAX
 1110 D(I)=0.D0
      IF((ICON(6)+ICON(7)).EQ.0) GO TO 1112
      CALL CUDATE(IDAT)
      WRITE(6,9019)ALPHA,IDAT
      WRITE(6,9024)
      WRITE(6,9022)
 1112 CONTINUE
      JX1=(1.0-PHASEF(JTRT(1)+IS(1)))/2.0
C
C     J LOOP
C
      DO 2000 JA=1,JMAX
      DO 1900 JPAR=1,NPAR
C modified by (JE,ME)
      READ(8) NLL,NLI
     1       ,(NTAB(L),L=1,NLL),(LTAB(L),L=1,NLL),(JTAB(L),L=1,NLL)
      K=NLI(1)
      IF(K.EQ.0) GO TO 1900
      MAXA=NLL+NLL
      DO 1120 N1=1,NLL
      I=NTAB(N1)
      L=LTAB(N1)+1
      IND=(I-1)*LMAX
      ARG=S(IND+L)-S(IND+1)
      ATAB(N1)=COS(ARG)
      BTAB(N1)=SIN(ARG)
      DET=SQRT(DFLOAT(L+L-1)*FMU(I)*FK(1)/(FK(I)*FMU(1)))
     1*SQRT(DRF(I)/DRF(1))**3/FK(I)
      FTAB(N1)=DET
 1120 CONTINUE
C
C     SET UP AND INVERT MATRIX FOR SCATTERING AMPLITUDES
C     COEFFICIENT MATRIX IS SOLVED BY PARTITIONS
C     (A2*A1(INV)*H1-H2)*SMAT=-(A2*A1(INV)*F1-F2)
C
C     FORM A1(TR,INV)*A2(TR) = AMAT
C
      DO 1140 N2=1,NLL
      NK=N2+N2-1
C modified by (JE,ME)
      READ(8) (ZI(I),I=1,MAXA)
      DO 1130 N1=1,NLL
      NJ=N1+N1-1
      CMAT(NK  ,N1)=ZI(NJ  )
      CMAT(NK+1,N1)=ZI(NJ+1)
 1130 CONTINUE
C modified by (JE,ME)
      READ(8) (ZI(I),I=1,MAXA)
      DO 1135 N1=1,NLL
      NJ=N1+N1-1
      AMAT(NK  ,N1)=ZI(NJ  )
      AMAT(NK+1,N1)=ZI(NJ+1)
 1135 CONTINUE
 1140 CONTINUE
C
      CALL INVERT(NLL,NLL,CMAT,AMAT,NLL2)
C
C
C      TRANSPOSE AMAT INTO CMAT
C
      DO 1160 N1=1,NLL
      NJ=N1+N1-1
      DO 1150 N2=1,NLL
      NK=N2+N2-1
      CMAT(NK  ,N1)=AMAT(NJ  ,N2)
      CMAT(NK+1,N1)=AMAT(NJ+1,N2)
 1150 CONTINUE
 1160 CONTINUE
      DO 1200 N2=1,NLL
      I=NTAB(N2)
      L=LTAB(N2)+1
      IND=(I-1)*LMAX+L
      CII(1)=ATAB(N2)*FTAB(N2)
      CII(2)=BTAB(N2)*FTAB(N2)
      D1           =  CII(1)*G (IND)-CII(2)*F (IND)
      D2           =  CII(1)*F (IND)+CII(2)*G (IND)
      D3           =  CII(1)*GP(IND)-CII(2)*FP(IND)
      D4           =  CII(1)*FP(IND)+CII(2)*GP(IND)
      DO 1190 N1=1,NLL
      NJ=N1+N1-1
      AMAT(NJ  ,N2)=CMAT(NJ  ,N2)*D1-CMAT(NJ+1,N2)*D2
      AMAT(NJ+1,N2)=CMAT(NJ+1,N2)*D1+CMAT(NJ  ,N2)*D2
      IF(I.NE.1) GO TO 1189
      D5=-CMAT(NJ  ,N2)*F (L)
      D6=-CMAT(NJ+1,N2)*F (L)
      CMAT(NJ  ,N2)=D5*CII(1)-D6*CII(2)
      CMAT(NJ+1,N2)=D5*CII(2)+D6*CII(1)
 1189 CONTINUE
 1190 CONTINUE
      NK=N2+N2-1
      AMAT(NK  ,N2)=AMAT(NK  ,N2)-D3
      AMAT(NK+1,N2)=AMAT(NK+1,N2)-D4
      IF(I.NE.1) GO TO 1200
      CMAT(NK  ,N2)=CMAT(NK  ,N2)+CII(1)*FP(L)
      CMAT(NK+1,N2)=CMAT(NK+1,N2)+CII(2)*FP(L)
 1200 CONTINUE
C
      CALL INVERT(NLL,K  ,AMAT,CMAT,NLL2)
C
      JT=JA+JA-2+JX1
      FJ=JT+1
      DO 1600 N1=1,K
      NN=N1+N1-1
      E(NN  )=CMAT(NN  ,N1)
      E(NN+1)=CMAT(NN+1,N1)
      JX=JTAB(N1)
      LX=LTAB(N1)+LTAB(N1)
      FL=LTAB(N1)+LTAB(N1)+1
      DO 1500 N2=1,NLL
      NK=N2+N2-1
      JJ=NTAB(N2)
      JY=JTAB(N2)
      LZ=LTAB(N2)
      LY=LTAB(N2)+LTAB(N2)
      IY=JTRT(JJ)
      MAXI=JTRT(JJ)+1
      CII(1)=(ATAB(N2)**2-BTAB(N2)**2)/FK(JJ)
      CII(2)= 2.0*ATAB(N2)*BTAB(N2)   /FK(JJ)
      D1=CMAT(NK  ,N1)*CII(1)-CMAT(NK+1,N1)*CII(2)
      D2=CMAT(NK  ,N1)*CII(2)+CMAT(NK+1,N1)*CII(1)
C
C     REACTION AND TOTAL INELASTIC XSECT
C
      FACT=FJ*FLOAT(LZ+LZ+1)/FL
      CMR=CMAT(NK  ,N1)
      CMI=CMAT(NK+1,N1)
      TEMP=0.0
      IF(ABS(CMR).GT.1.0E-18) TEMP=CMR*CMR
      IF(ABS(CMI).GT.1.0E-18) TEMP=TEMP+CMI*CMI
      IF(JJ.EQ.1) GO TO 1290
      SIG(JJ)=SIG(JJ)+TEMP*FACT
      GO TO 1299
 1290 CONTINUE
      IF(N1.NE.N2) GO TO 1299
      SIG(1)=SIG(1)+(CMAT(NK+1,N1)-TEMP)*FACT
 1299 CONTINUE
      CMAT(NK  ,N1)=D1
      CMAT(NK+1,N1)=D2
C
C     STORE S-MATRIX ELEMENTS
C
      JR=NS(1)
      JS=NS(JJ)
      PHASE=PHASEF((LX+IS(1)+IX-LY-IS(JJ)-IY)/2)
      INDEX=MAXTAB(1)*MAXTAB(JJ)*LZ+INDTAB(JJ)
      INDEX=INDEX+INDEX+1
      IND=INDEX+MAX1*MAXI*JR*JS*2-2
      MM=-JTRT(1)
      IF(MM.EQ.0) GO TO 1310
      M1=(MAX1+1)/2
      MI=MAXI
      GO TO 1315
 1310 CONTINUE
      M1=1
      MI=(MAXI+1)/2
 1315 CONTINUE
      DO 1490 MX=1,M1
      MMS=-JTRT(JJ)
      DO 1480 MY=1,MI
      MS=-IS(1)
      DO 1470 M =1,JR
      MSP=-IS(JJ)
      DO 1460 MP=1,JS
      ML=MS+MM-MSP-MMS
      MLT=IABS(ML/2)
      IF(MLT.GT.LZ) GO TO 1458
      PHAS=1.D0
      IF(ML.LT.0) PHAS=PHASEF(MLT)
C
      COEFF=VCC(LX,IS(1),JX,0,MS)*VCC(JX,IX,JT,MS,MM)
     1*VCC(LY,IS(JJ),JY,ML,MSP)*VCC(JY,IY,JT,ML+MSP,MMS)
     2*EXP(0.5*YXFCT(LZ+MLT,LZ-MLT))*PHAS*FLOAT(LY+1)
C
      D(INDEX  )=D(INDEX  )+D1*COEFF
      D(INDEX+1)=D(INDEX+1)+D2*COEFF
C
C     FILL OUT D TABLE
C
      IF((JTRT(1)+JTRT(JJ)).EQ.0) GO TO 1458
      IF((M1.NE.MAX1).AND.(MM.EQ.0)) GO TO 1458
      IF((MI.NE.MAXI).AND.(MMS.EQ.0))GO TO 1458
      PHAS=PHASE*PHAS
      IF(ML.GT.0)PHAS=PHASE*PHASEF(MLT)
      D(IND  )=D(INDEX  )*PHAS
      D(IND+1)=D(INDEX+1)*PHAS
 1458 CONTINUE
      INDEX=INDEX+2
      IND=IND-2
      MSP=MSP+2
 1460 CONTINUE
      MS=MS+2
 1470 CONTINUE
      MMS=MMS+2
 1480 CONTINUE
      MM=MM+2
 1490 CONTINUE
 1500 CONTINUE
      IF(ICON(7).EQ.0) GO TO 1600
      K1=1
      K3=1
      DO 1590 I=1,NLTR
      IF(NLI(I).EQ.0) GO TO 1590
      K2=K1+NLI(I)+NLI(I)-1
      WRITE(6,9025)I,JT,LTAB(K3),(CMAT(M,N1),M=K1,K2)
      K1=K2+1
      K3=K3+NLI(I)
 1590 CONTINUE
 1600 CONTINUE
      IF(ICON(6).EQ.0) GO TO 1900
      I=1
      K2=K+K
      WRITE(6,9023)I,JT,LTAB(1),(E(M),M=1,K2)
      IF(ICON(6).NE.2) GO TO 1900
      DO 1610 M=1,K
      MK=M+M-1
      T1=E(MK+1)
      E(MK+1)=2.0*E(MK  )
      E(MK  )=1.0-2.0*T1
 1610 CONTINUE
      WRITE(6,9023)I,JT,LTAB(1),(E(M),M=1,K2)
 1900 CONTINUE
 2000 CONTINUE
      DO 2100 I=1,NLTR
      SIG(I)=SIG(I)*12.5663706/(FLOAT((JTRT(1)+1)*NS(1))*FK(I)**2)
 2100 CONTINUE
      IF(ICON( 8).EQ.0) GO TO 2610
C
C     WRITE OUT AMPLITUDES
C
      MAX1=JTRT(1)+1
      CALL CUDATE(IDAT)
      WRITE(6,9019)ALPHA,IDAT
      IND=1
      DO 2605 I=1,NLTR
      MAXI=JTRT(I)+1
      KST=NS(1)*MAX1*NS(I)*MAXI
      KST=KST+KST
      WRITE(6,9060)I
      DO 2590 L=1,LMAX
      LTEMP=L-1
      K1=IND
      K2=IND+KST-1
      WRITE(6,9061)LTEMP,(D(J),J=K1,K2)
      IF(I.EQ.ICON(8)) WRITE(2,9062)LTEMP,(D(J),J=K1,K2)
      IND=K2+1
 2590 CONTINUE
 2605 CONTINUE
 2610 CONTINUE
 3000 CONTINUE
      RETURN
 9019 FORMAT(1H1,17X,15A4,4X,I2,2(1H/,I2),I4,2(1H.,I2))
 9022 FORMAT(12H0CHN 2*J  L1)
 9023 FORMAT(3I4,10F12.7/(12X,10F12.7))
 9024 FORMAT(22H SCATTERING AMPLITUDES )
 9025 FORMAT(3I4,10E12.4/(12X,10E12.4))
 9060 FORMAT(5H0CHN=,I2,100(1H*))
 9061 FORMAT(5H  L2=,I3,2X,10E12.4/(10X,10E12.4))
 9062 FORMAT(I4,10E12.5/(4X,10E12.5))
      END
      SUBROUTINE HAN(R ,RP,E,L,H,F,FP,G,GP,S)
C modified by (JE,ME)
C      IMPLICIT REAL*8(A-H,O-Z)
C     USE THE NEXT 1 CARDS FOR 32-BIT MACHINES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     ABS(XXX)=DABS(XXX)
      DIMENSION F(50),FP(50),G(50),GP(50),S(50)
      FF(X)=EXP(-X-E*LOG(2.0*X))
C modified by (JE,ME)
      SAVE
      DO 40 J=1,2
      FL=J-1
      T3=1.D0
      T4=1.D0
      T5=1.D0
      T6=1.D0
      G (J)=1.D0
      GP(J)=1.D0
      DO 30 I=1,20
      T3=T3*(FL+1.0-E-T4)*(FL+E+T4)/(2.0*T4)
      T5=T5/R
      T6=T6/RP
      G (J)=G (J)+T5*T3
      GP(J)=GP(J)+T6*T3
      T4=T4+1.D0
      TEST=ABS(T5*T3)
      IF(TEST.LT.1.0D-05) GO TO 35
   30 CONTINUE
      WRITE(6,9000)TEST
   35 CONTINUE
      G (J)=G (J)*FF(R )
      GP(J)=GP(J)*FF(RP)
   40 CONTINUE
      F (1)=0.D0
      F (2)=0.D0
      FP(1)=0.D0
      FP(2)=0.D0
      S(1)=0.D0
      S(2)=0.D0
      LL=L
      T1=E*E
      T2=1.D0
      T3=2.D0
      DO 50 I=3,LL
      T4=T2+T3
      T5=T2*T3
      T6=T3*SQRT(T2**2+T1)
      T7=T2*SQRT(T3**2+T1)
      G (I)=(T4*(E+T5/R )*G (I-1)+T6*G (I-2))/T7
      GP(I)=(T4*(E+T5/RP)*GP(I-1)+T6*GP(I-2))/T7
      T2=T2+1
      T3=T3+1
      F (I)=0.D0
      FP(I)=0.D0
      S(I)=0.D0
   50 CONTINUE
      RETURN
 9000 FORMAT(25H0BOUNDSTATE CONVERGES TO ,E12.3,8H IN HAN   )
      END
      SUBROUTINE XSECT(D,PLM,SR,SI)
      IMPLICIT REAL*8(A-H,O-Z)
C     SIN(XXX)=DSIN(XXX)
C     COS(XXX)=DCOS(XXX)
C     ALOG(XXX)=DLOG(XXX)
C     SQRT(XXX)=DSQRT(XXX)
C     FLOAT(III)=DFLOAT(III)
C     DOUBLE PRECISION BETA
      DIMENSION D(1),PLM(1),SR(1),SI(1)
      COMMON THETAN,THETA1,DTHETA,DRT,RMAX,ALPHA(15),IDAT(6)
     1,ICON(20),LMAX,KT,KMAX,LLCO,JOUT,NLTR,IBUFF,IWORD
     2,KLAM( 8),LTRT( 8),JTRT( 8),NADD( 8),LAM(16),ISP(16),JSP(16)
     3,ILK(16),JPRT( 8)
      COMMON EN( 8),FM( 8),Z( 8),FMA( 8),ZA( 8),RY( 8),FS( 8),ECM( 8)
     1,FMU( 8),RC( 8),FK2( 8),FK( 8),ETA( 8),AFACT( 8),DRF( 8),FN( 8)
     2,PNLOC( 8),DUM( 8),SIG( 8)
     3,IS( 8),NS( 8),NLI( 8),LKTAB(32),ICUP(32),JCUP(32)
     4,IN( 8),INN( 8, 8),INDTAB(10)
      COMMON/LOCAL/  TSIG(256),POL(9,256),SUMR(5,5),SUMI(5,5),POLX(9)
     1,BUFC(1781)
      COMMON/BUFSPC/ BUFFR(800)
      DATA NMAX/256/
C
      IF(ICON(5).EQ.0) GO TO 10
C
C     WRITE AMPLITUDES ON TAPE 4
C
      WRITE(4)THETAN,THETA1,DTHETA,ALPHA,IDAT,LMAX,NLTR,LTRT,JTRT,NADD
     1,JPRT,FS,FK,ETA,SIG,IS,NS,INDTAB
      N=INDTAB(NLTR+1)
      N=N+N
      WRITE(4)(D(I),I=1,N)
   10 CONTINUE
      JMAX=LMAX
      JMAX2=JMAX+JMAX
      NLTR1=NLTR+1
      NTHETA=NINT(THETAN)
      IF(NTHETA.EQ.0) GO TO 2100
      MAX1=JTRT(1)+1
      MAXM=0
      DO 40 I=1,NLTR
      MAXM=MAX0(MAXM,(JTRT(1)+JTRT(I)+IS(1)+IS(I))/2+1)
   40 CONTINUE
      REWIND 9
      NX=MAXM*LMAX
      IWORD=0
      THETA=THETA1
      DO 50 M=1,NTHETA
      CALL LGNDR(PLM,MAXM,LMAX,THETA)
C modified by (JE,ME)
      WRITE(9) (PLM(I),I=1,NX)
      THETA=THETA+DTHETA
   50 CONTINUE
C modified by (JE,ME)
      WRITE(9) NX
      IF(IWORD.EQ.IBUFF) GO TO 51
      IWORD=IBUFF-1
  51  CONTINUE
C
C     CHANNEL LOOP
C
      DO 1000 I=1,NLTR
      IF(FK2(I).LT.0.0) GO TO 1000
      CALL CUDATE(IDAT)
      WRITE(6,9019) ALPHA,IDAT
      WRITE(6,9022) I,FM(I),Z(I),FMA(I),ZA(I),EN(I),IS(I),LTRT(I)
      WRITE(6,9024)
      IF((ICON(13).EQ.9).OR.(ICON(13).EQ.I)) WRITE(7,9030) I,ALPHA
      REWIND 9
      IWORD=IBUFF
C
C     ANGLE LOOP
C
      THETA=THETA1
      DO 900 N=1,NTHETA
      IF(N.GT.NMAX) GO TO 910
C modified by (JE,ME)
      READ(9) (PLM(M),M=1,NX)
      IF(I.NE.1) GO TO 75
C
C     CALCULATE COULOMB AMPLITUDE FOR CHANNEL 1
C
      IF(ETA(1).EQ.0.0) GO TO 75
      SS=SIN(THETA/114.5915  )**2
      IF(SS.LE.1.0D-10) GO TO 75
      TL=LOG(SS)
      ARG=ETA(1)*TL
      CFACT=ETA(1)/(2.0*FK(1)*SS)
      CL=COS(ARG)*CFACT
      SL=SIN(ARG)*CFACT
      GO TO 76
   75 CONTINUE
      CL=0.D0
      SL=0.D0
   76 CONTINUE
      IND=INDTAB(I)
      IND=IND+IND
      MAXI=JTRT(I)+1
      KST=NS(1)*MAX1*NS(I)*MAXI
      DO 150 J=1,KST
      SR(J)=0.D0
      SI(J)=0.D0
  150 CONTINUE
      TSIG(N)=0.D0
      DO 155 M=1,9
      POL(M,N)=0.D0
  155 CONTINUE
      INC=IND+1
      JR=NS(1)
      JS=NS(I)
      DO 500 L=1,LMAX
      LZ=L-1
      IT=1
      MM=-JTRT(1)
      DO 490 MX=1,MAX1
      MMS=-JTRT(I)
      DO 480 MY=1,MAXI
      MS=-IS(1)
      DO 470 M =1,JR
      MSP=-IS(I)
      DO 460 MP=1,JS
      ML=MS+MM-MSP-MMS
      MLT=IABS(ML/2)
      INDX=MLT*LMAX+LZ+1
      SR(IT)=SR(IT)+D(INC  )*PLM(INDX)
      SI(IT)=SI(IT)+D(INC+1)*PLM(INDX)
      INC=INC+2
      IT=IT+1
      MSP=MSP+2
  460 CONTINUE
      MS=MS+2
  470 CONTINUE
      MMS=MMS+2
  480 CONTINUE
      MM=MM+2
  490 CONTINUE
  500 CONTINUE
      IT=1
      DO 790 MX=1,MAX1
      DO 780 MY=1,MAXI
      DO 580 M =1,5
      DO 580 MP=1,5
      SUMR(M,MP)=0.D0
      SUMI(M,MP)=0.D0
  580 CONTINUE
      DO 670 M =1,JR
      DO 660 MP=1,JS
      IF(I.NE.1) GO TO 610
      IF(M .NE.MP) GO TO 610
      IF(MX.NE.MY) GO TO 610
C
C     ADD COULOMB TERM TO CHANNEL 1 AMPLITUDE
C
      SR(IT)=SR(IT)-CL
      SI(IT)=SI(IT)+SL
  610 CONTINUE
      SUMR(MP,M)=SUMR(MP,M)+SR(IT)
      SUMI(MP,M)=SUMI(MP,M)+SI(IT)
      TSIG(N)=TSIG(N)+SR(IT)**2+SI(IT)**2
      IT=IT+1
  660 CONTINUE
  670 CONTINUE
C
C     CALCULATE POLARIZATION
C
      CALL POLFCT(POLX,SUMR,SUMI,JS,JR)
      DO 770 M=1,9
      POL(M,N)=POL(M,N)+POLX(M)
  770 CONTINUE
  780 CONTINUE
  790 CONTINUE
C
C
      IF(TSIG(N).EQ.0.0) GO TO 820
      DO 815 M=1,9
      POL(M,N)=POL(M,N)/TSIG(N)
  815 CONTINUE
      TSIG(N)=TSIG(N)/FLOAT(MAX1*NS(1))
      GO TO 830
  820 CONTINUE
      DO 825 M=1,9
      POL(M,N)=0.D0
  825 CONTINUE
  830 CONTINUE
      RATIO=0.D0
      IF(I.NE.1) GO TO 855
      IF(ETA(1).EQ.0.0) GO TO 855
      RATIO=1.0
      IF(SS.LT.1.0D-10) GO TO 855
      RATIO=TSIG(N)/(CFACT*CFACT)
  855 CONTINUE
      WRITE(6,9900) THETA,TSIG(N),RATIO,(POL(M,N),M=1,9)
      IF((ICON(13).EQ.9).OR.(ICON(13).EQ.I)) WRITE(7,9031) THETA,
     1  TSIG(N),(POL(M,N),M=1,9)
      IF(I.NE.1) GO TO 890
      IF(ETA(1).EQ.0.0) GO TO 890
      TEMP=TSIG(N)
      TSIG(N)=RATIO
      RATIO=TEMP
  890 CONTINUE
      THETA=THETA+DTHETA
  900 CONTINUE
  910 CONTINUE
      IF(I.EQ.1) WRITE(6,9904) SIG(I)
      IF(I.NE.1) WRITE(6,9905) SIG(I)
      IF(ICON(9).GT.5) GO TO 1000
      IF(ICON(9).EQ.0) ICON(9)=3
      N=MIN0(NTHETA,NMAX)
      CALL CUDATE(IDAT)
      WRITE(6,9019)ALPHA,IDAT
      WRITE(6,9901)I,LTRT(I)
      IF((I.NE.1).OR. (ETA(1).EQ.0.0)) WRITE(6,9028)
      IF((I.EQ.1).AND.(ETA(1).NE.0.0)) WRITE(6,9029)
      CALL DWPLOT(N,ICON(9),TSIG(1),DTHETA,THETA1)
 1000 CONTINUE
 2100 CONTINUE
      RETURN
 9019 FORMAT(1H1,17X,15A4,4X,I2,2(1H/,I2),I4,2(1H.,I2))
 9022 FORMAT(8H0CHANNEL,I3,4X,4HFM =,F8.4,5H  Z =,F5.1,3X,2HON,7H  FMA =
     1,F9.4,6H  ZA =,F5.1,4X,6HELAB =,F9.4,04X,8H2*SPIN =,I4,4X,8H2*JTRT
     2 =,I4)
 9024 FORMAT(7H0 THETA,4X,10HSIGMA F**2,4X,10HSIGMA/RUTH,7X,4HPOL.,6X,3H
     1T22,7X,3HT21,7X,3HT20,9X,4HA.P.,6X,3HT22,7X,3HT21,7X,3HT20,7X,7HSF
     2-PROB)
 9028 FORMAT(1H+,12X,10HSIGMA F**2)
 9029 FORMAT(1H+,12X,10HSIGMA/RUTH)
 9030 FORMAT(7HCHANNEL,I3,10X,15A4)
 9031 FORMAT(F6.2,1PE11.4,9(0PF7.4))
 9900 FORMAT(F7.2,2E14.4,2(2X,4F10.5),2X,F10.5)
 9901 FORMAT(8H0CHANNEL,I3,12X,10H   2*SPIN=,I3)
 9904 FORMAT(9H0REACSIG ,E12.4)
 9905 FORMAT(9H0TOTALSIG,E12.4)
      END
      SUBROUTINE INVERT(MAX,N,A,C,NLL2)
      IMPLICIT REAL*8(A-H,O-Z)
C     COMPLEX*16 A(40,40),C(40,40),B
C     COMPLEX    A(40,40),C(40,40),B
C     DO 100 I=1,MAX
C     B=1.0/A(I,I)
C     DO 10 K=1,N
C     C(I,K)=C(I,K)*B
C  10 CONTINUE
C     DO 20 J=I,MAX
C  20 A(I,J)=A(I,J)*B
C     DO 50 J=1,MAX
C     IF(I.EQ.J) GO TO 50
C     B=A(J,I)
C     DO 35 K=1,N
C     C(J,K)=C(J,K)-C(I,K)*B
C  35 CONTINUE
C     DO 40 K=I,MAX
C  40 A(J,K)=A(J,K)-A(I,K)*B
C
C     CODE TO HANDLE DOUBLE PRECISION COMPLEX ON THE DEC-10.
C     IF ABOVE CODING IS RESTORED, DIMENSIONS AND NLL2 MUST BE FIXED UP.
C
C     DOUBLE PRECISION A,C,AR,AI,BR,BI,CR,CI,TERM,FACT
      DIMENSION A(NLL2,1),C(NLL2,1)
      DO 100 I=1,MAX
      I1=I+I-1
      AR=A(I1  ,I)
      AI=A(I1+1,I)
      FACT=1.D0
      TERM=DMAX1(DABS(AR),DABS(AI))
      IF(TERM.LE.5.D18) GO TO 2
      FACT=TERM/5.D18
      GO TO 3
    2 CONTINUE
      IF(TERM.GE.1.D-10) GO TO 5
      FACT=TERM/1.D-10
    3 CONTINUE
      AR=AR/FACT
      AI=AI/FACT
    5 CONTINUE
      TERM=AR*AR+AI*AI
      BR= AR/TERM
      BI=-AI/TERM
      BR=BR/FACT
      BI=BI/FACT
      DO 10 K=1,N
      CR=C(I1  ,K)
      CI=C(I1+1,K)
      C(I1  ,K)=CR*BR-CI*BI
      C(I1+1,K)=CR*BI+CI*BR
   10 CONTINUE
      DO 20 J=I,MAX
      AR=A(I1  ,J)
      AI=A(I1+1,J)
      A(I1  ,J)=AR*BR-AI*BI
      A(I1+1,J)=AR*BI+AI*BR
   20 CONTINUE
      DO 50 J=1,MAX
      IF(I.EQ.J) GO TO 50
      J1=J+J-1
      BR=A(J1  ,I)
      BI=A(J1+1,I)
      DO 35 K=1,N
      CR=C(I1  ,K)
      CI=C(I1+1,K)
      C(J1  ,K)=C(J1  ,K)-CR*BR+CI*BI
      C(J1+1,K)=C(J1+1,K)-CR*BI-CI*BR
   35 CONTINUE
      DO 40 K=I,MAX
      AR=A(I1  ,K)
      AI=A(I1+1,K)
      A(J1  ,K)=A(J1  ,K)-AR*BR+AI*BI
      A(J1+1,K)=A(J1+1,K)-AR*BI-AI*BR
   40 CONTINUE
   50 CONTINUE
  100 CONTINUE
      RETURN
      END
