      SUBROUTINE PACKPOT(RCC,CGAMA,QSQ)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(5,3),B(5,3),C(5,2),RC(5),P(5),VSCF(5),VCOR(5)
      DIMENSION V(5),D(5)
      DIMENSION X(96)
      DATA KASE/0/
      DATA NC/5/,II/1/,K/2/,TMU/1.0D+00/
      DATA A(1,1)/0.172995D+02/,A(1,2)/-0.982289D+00/,
     X A(1,3)/-0.917494D-01/
      DATA A(2,1)/0.655752D+02/,A(2,2)/-0.112538D+01/,
     X A(2,3)/-0.766468D-01/
      DATA A(3,1)/0.105137D+03/,A(3,2)/-0.144877D+01/,
     X A(3,3)/-0.522417D-01/
      DATA A(4,1)/0.178531D+03/,A(4,2)/-0.190354D+01/,
     X A(4,3)/-0.239120D-01/
      DATA A(5,1)/0.314919D+03/,A(5,2)/-0.240207D+01/,
     X A(5,3)/ 0.000000D+00/
      DATA B(1,1)/0.122126D+00/,B(1,2)/-0.628088D+00/,
     X B(1,3)/-0.470607D-01/
      DATA B(2,1)/0.569431D-01/,B(2,2)/-0.663263D+00/,
     X B(2,3)/-0.399192D-01/
      DATA B(3,1)/0.252919D-01/,B(3,2)/-0.976097D+00/,
     X B(3,3)/-0.215588D-01/
      DATA B(4,1)/0.162599D-01/,B(4,2)/-0.116525D+01/,
     X B(4,3)/-0.280915D-01/
      DATA B(5,1)/0.636620D-01/,B(5,2)/-0.188254D+01/,
     X B(5,3)/ 0.000000D+00/
      DATA C(1,1)/16.75D+00/,C(1,2)/278.0D+00/
      DATA C(2,1)/3.88D+00/,C(2,2)/290.0D+00/
      DATA C(3,1)/0.0D+00/,C(3,2)/57.0D+00/
      DATA RC/5.72050D+00,6.27735D+00,6.39141D+00,99999.999D+00,99999.99
     X9D+00/
      DATA D/0.122126D+00,0.569431D-01,0.252919D-01,0.162599D-01,
     X 0.636620D-01/
      DATA XSAV/0.0D+00/
      DATA ISKIP,ISAV /2*0/
      X(1)=CGAMA
      R=RCC/0.52917706D+00
      IF(R.LT.2.D0)R=2.D0
C-------------------------------------------------------
C  READ IN POTENTIAL PARAMETERS ON FIRST CALL
C-------------------------------------------------------
      IF (ISKIP.NE.0) GO TO 50
C-------------------------------------------------------
C  WRITE OUR COEFFS IN TABLE FORM
C-------------------------------------------------------
      WRITE (6,15000)
           DO 10 J=1,3
   10      WRITE (6,20000) (A(I,J),I=1,NC)
      WRITE (6,25000)
           DO 20 J=1,3
   20      WRITE (6,20000) (B(I,J),I=1,NC)
C-------------------------------------------------------
C  I LABELS THE LEGENDRE COEFFS.  DIFFER IN SYM AND UNSYM CASES.
C   J=1,2 GO WITH THE R-6 AND R-8 TERMS IN VDW POTENTIAL
C-------------------------------------------------------
      WRITE (6,40000)
           DO 30 J=1,2
   30      WRITE (6,45000) (C(I,J),I=1,3)
      WRITE (6,45000) (RC(I),I=1,NC)
C-------------------------------------------------------
C  D(I) IS B(I,1) PRIME, THE PART OF B WHICH MATCHES VDW AT RC.  FOR OLD
C  EG SURFACES WHICH READ IN B PRIME, SET D(I)=B(I,1).
C-------------------------------------------------------
      IF (KASE.EQ.1) THEN
           DO 40 I=1,NC
   40      D(I)=B(I,1)
      END IF
      WRITE (6,55000)
      WRITE (6,20000) (D(I),I=1,NC)
      ISKIP=1
   50 CONTINUE
      IF (XSAV.EQ.X(II)) GO TO 70
C-------------------------------------------------------
C  CALCULATE THE LEGENDRE POLYS NEEDED ON THE FIRST PASS AT A GIVEN ANGL
C-------------------------------------------------------
           DO 60 J=1,NC
C-------------------------------------------------------
C  FOR SYMM MOLEC N=2*J-2.  FOR UNSYM MOLEC N=J-1
C-------------------------------------------------------
           N=2*J-2
           S=1.D+00
           IF (MOD(N,2).EQ.1) S=-S
   60      P(J)=PLM(N,0,X(II))/DSQRT(N+0.5D+00)*S
      XSAV=X(II)
   70 CONTINUE
C-------------------------------------------------------
C  CALCULATE POTENTIAL AT R FOR THE GIVEN ANGLE
C-------------------------------------------------------
      RS=R*R
      RM2=1.D+00/RS
      VT=0.0D+00
           DO 100 I=1,NC
           VSCF(I)=A(I,1)*DEXP(A(I,2)*R+A(I,3)*RS)
           IF (R.GT.RC(I)) GO TO 80
           VCOR(I)=-B(I,1)*DEXP(B(I,2)*R+B(I,3)*RS)
           GO TO 90
   80      CONTINUE
           RM6=RM2**3
           RM8=RM2*RM6
           VCOR(I)=-C(I,1)*RM6-C(I,2)*RM8
           IF (D(I).EQ.B(I,1)) GO TO 90
           VCOR(I)=VCOR(I)-(B(I,1)-D(I))*DEXP(B(I,2)*R+B(I,3)*RS)
   90      CONTINUE
           V(I)=VSCF(I)+VCOR(I)
           VT=VT+V(I)*P(I)
  100      CONTINUE
      IF (K.NE.2) GO TO 110
      QSQ=VT*TMU
      QSQ=QSQ*27210.6D+00
      RETURN
  110 CONTINUE
      QSQ=E-XL*RM2/TMU-VT
      IF (K.EQ.0) RETURN
      IF (QSQ.GE.0.D+00) GO TO 120
      WRITE (6,60000) E,XL,R
      Q=0.D+00
      RETURN
  120 CONTINUE
      Q=DSQRT(TMU*QSQ)
      RETURN
15000 FORMAT(1H0,' A COEFFS'/)
20000 FORMAT(1H ,10E13.6)
25000 FORMAT(1H0,' B COEFFS'/)
40000 FORMAT(1H0,'C COEFFS AND RC'/)
45000 FORMAT(1H ,3F14.7,4F10.2)
55000 FORMAT(1H0,' B PRIMES'/)
60000 FORMAT(1H ,'WARNING...FOR E=',F10.7,' AND (L+0.5)**2=',F10.2,' AT
     XR=',F10.5,' QSQ IS NEGATIVE, INDICATING 3 TURNING POINTS.')
      END
