C N. R. BADNELL  UoS v1.11   -   QUB vx.x                       28/11/08
C
C     PROGRAM NXRAD/STG2NX
C
C *** THIS IS THE RADIAL STAGE OF RMATRX NX
C
      PROGRAM MAIN
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C SUN
      REAL*4 TARRY(2),TIME
C
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB,IDIM14=MZFAC)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM3=IDIM9+1,KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9,
     *          KDIM6=(IDIM4+1)/2,KDIM7=1,
     *          KDIM8=1,
     *          KDIM9=IDIM1+IDIM1-1,KDIM10=IDIM7*2,
     *          KDIM11=1,KDIM12=IDIM1*IDIM12,
     *          KDIM13=1)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      COMMON /BASDER/ FM,TLC,WR,IMEMM,ITST,JR,KM,MMM,NBTP1,NG
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA
      COMMON /BNDORB/ P(IDIM9,KDIM10)
      COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /COPY  / ICOPY1,ICOPY2,ITOTAL,IPSEUD
      COMMON /CORE  / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2),
     *                 MAXNC(IDIM2),MAXPN(IDIM2),ICHECK
      COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4
      COMMON /FACT  / GAMMA(IDIM14)
      COMMON /FUNVAL/ FRH(KDIM3),U(KDIM3),X
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC
      COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI
     *                ,BVALUE
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /POTVAL/ POVALU(KDIM10)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5)
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS
      COMMON /YKSTOR/ YK(IDIM7),RK(IDIM7),RK1(IDIM7),JN2,JL2,JN4,JL4,JKM
C
 1000 FORMAT(TR30,' END OF NXRAD'/TR30,' ------------'///)
C     READ IN AND WRITE OUT THE BASIC DATA. INITIALISE CERTAIN ARRAYS
C     AND VARIABLES
      NFACT=IDIM14
      CALL SHRIEK(NFACT)
C
C      REREAD STGANG DATA AND
C      CALCULATE RANGE OF CONTINUUM ANGULAR MOMENTUM LRANGE2
C      FROM CONFIGURATION DATA
C
      CALL RESETL(MJS,MINIM,ICONTN)
C
C      FOR A CONTINUATION FROM RMATRX READ IN DATA, BOUND ORBITALS,
C      AND THE POTENTIAL FROM RMATRX FILE NX1
C
      IF (ICONTN .EQ. 1) THEN
         CALL STGRRX(MINIM,MJS)
      ELSE
         CALL STGRRD (MJS)
C     EVALUATE THE BOUND ORBITALS AT THE MESH POINTS.
         CALL EVALUE (MINIM)
C     GENERATE THE STATIC POTENTIAL
         CALL POTF
      END IF
C
C     EVALUATE THE  CONTINUUM ORBITALS
      CALL BASORB(MINIM)
C
C      CALCULATE THE BUTTLE CORRECTION
C
      CALL NEWBUT(MINIM,MJS)
C
C      WRITE SERIAL FILE FOR STGHAM
C
      CALL WRITAP
C     GENERATE THE RADIAL INTEGRALS
      CALL GENINT (MINIM)
      WRITE(IWRITE,1000)
C
C SUN
      DUM=DTIME(TARRY)
      TIME=TARRY(1)
C CRAY
CRAY  CALL SECOND(TIME)
C
      TIME=TIME/60.0
      WRITE(IWRITE,999) TIME
  999 FORMAT(//1X,9HCPU TIME=,F9.3,4H MIN)
C
C
      STOP
      END

      SUBROUTINE BASFUN(NBT,LC,NODES,RA,BSTO,WINIT,DELTA,ETA)
C
C      NEW NUMERICAL R-MATRIX ORBITAL ROUTINE
C***********************************************************************
C
C                      OPERATING INSTRUCTIONS
C
C***********************************************************************
C
C  1. THE USER MUST PROVIDE THE FOLLOWING INPUT DATA...
C
C     NBT....... THE NUMBER OF FUNCTIONS TO WHICH THE SOLUTION IS TO BE
C                ORTHOGONALIZED...... IF NBT.GT.5 IS REQUIRED THE FIRST
C                INDEX IN THE ARRAYS  US,P,U,DU,FR,FRH,FRM,ALAMDA,DELT,
C                SDELT,UNAME,ADEL,ADL,BDL SHOULD BE INCREASED.
C
C     LC........ THE  ANGULAR MOMENTUM VALUE
C
C     RA........ THE BOUNDARY RADIUS
C
C     BSTO...... THE VALUE OF THE LOGARITHMIC DERIVATIVE AT X=RA
C
C     WINIT..... THE INITIAL ENERGY OR POTENTIAL MULTIPLE VALUE
C                        .... IF ETA=0.0 OR IS LESS THAN 1.0E-8 THE WAVE
C                FUNCTION  WILL BE EVALUATED AT THIS VALUE.   IF ETA
C                IS GREATER THAN 1.0E-8 (TYPICALLY ETA=0.00001) THEN THE
C                PROGRAM WILL ITERATE FROM WR OR PR=WINIT TO THE VALUE
C                OF WR OR PR WHICH GIVES A SOLUTION SATISFYING THE
C                LOGARITHMIC BOUNDARY CONDITION AT X=RA
C
C     DELTA..... THE INCREMENT IN THE ENERGY OR POTENTIAL MULTIPLE USED
C                FOR OBTAINING THE DERIVATIVE IN NEWTONS METHOD....THIS
C                SHOULD BE OF THE SAME ORDER AS ETA
C
C     ETA....... THE PROGRAM WILL STOP ITERATING TO FIND THE EIGENVALUE
C                WHEN THE CORRECTION TO THE ENERGY OR POTENTIAL MULTIPLE
C                BECOMES LESS THAN ETA.
C
C     IREAD..... THE INPUT PERIPHERAL NUMBER
C
C     IWRITE.... THE OUTPUT PERIPHERAL NUMBER
C
C     HINT...... THE BASIC INTEGRATION STEP LENGTH
C
C     NIX....... THE NUMBER OF CHANGES OF INTEGRATION STEP OR THE NUMBER
C                OF SUBINTERVALS INTO WHICH THE INTERVAL X=0 TO RA IS
C                DIVIDED
C
C     IHX(I),I=1,NIX.... THE MULTIPLE OF THE BASIC INTEGRATION STEP IN
C                        EACH SUBINTERVAL
C
C     IRX(I),I=1,NIX.... THE TOTAL NUMBER OF INTEGRATION STEPS TO THE
C                        END OF THE ITH SUBINTERVAL
C
C
C  2. THE USER MUST PROVIDE THE POTENTIAL FUNCTION AND ORTHOGONALISATION
C      FUNCTIONS AND STORE THEM IN THE ARRAYS POVALU(1200) AND P(6,1200)
C
C     ** THESE MUST BE EVALUATED AS FOLLOWS --- MUST BEING EMPHASISED **
C
C     THE ODD ELEMENTS POVALU(2N-1) AND P(I,2N-1),N=1,IRX(NIX), SHOULD
C     CONTAIN THE FUNCTION VALUES AT THE HALF-MESH POINTS.
C
C     THE EVEN ELEMENTS POVALU(2N) AND P(I,2N),N=1,IRX(NIX), SHOULD
C     CONTAIN THE VALUES AT THE MESH POINTS.
C
C***********************************************************************
C
C                            DEBUGGING  PARAMETERS
C
C***********************************************************************
C
C      IF THESE ARE SET EQUAL TO ZERO THERE IS NO INTERMEDIATE PRINT-OUT
C
C     NBUG1.....IF THIS IS NON-ZERO THE INTERMEDIATE INTEGRATIONS AND
C               ENERGIES ARE OUTPUT
C
C     NBUG3.....IF THIS IS NON-ZERO THE ARRAYS FOR THE DETERMINATION
C                OF THE MISMATCH ARE OUTPUT
C
C***********************************************************************
C
C                             OUTPUT RESULTS
C
C***********************************************************************
C
C     ORB(K) , K=1,IRX(NIX)+1 CONTAINS THE FINAL SOLUTION
C
C     NODES....  IS THE NUMBER OF NODES IN THE FUNCTION
C
C     EIGEN  IS THE EIGENVALUE
C
C     ALAMDA(I),I=1,NBTP1 CONTAINS THE NORMALISED LAGRANGE MULTIPLIERS
C
C     ANORMI IS THE STURMIAN NORMALISATION COEFFICIENT FOR R.GT.RA
C
C     BVALUE IS THE  LOGARITHMIC DERIVATIVE VALUE AT R=RA
C
C***********************************************************************
C
C                            OTHER VARIABLE AND ARRAY DEFINITIONS
C
C***********************************************************************
C
C     SEE THE LONG WRITE-UP
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C
      PARAMETER(LARGE=MZMAX,VSMALL=1.D1**(-MZMAX),QLARGE=1.D1**MZMX2)
C
      PARAMETER(IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG)
      PARAMETER(KDIM3=IDIM9+1,KDIM10=IDIM7*2)
C
      COMMON /BASDER/ FM,TLC,WR,IMEMM,ITST,JR,KM,MMM,NBTP1,NG
      COMMON /BNDORB/ P(IDIM9,KDIM10)
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /FUNVAL/ FRH(KDIM3),U(KDIM3),X
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI
     *                ,BVALUE
C *****
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS
      COMMON /POTVAL/ POVALU(KDIM10)
C
      DIMENSION AVAL(2),BNDRY(2),DBNDRY(2)
      DIMENSION ORTH(KDIM3),YIN(KDIM3),YOUT(KDIM3),DYIN(KDIM3),
     *          DYOUT(KDIM3),US(KDIM3,IDIM7),DUS(KDIM3,IDIM7),DU(KDIM3),
     *          FR(KDIM3),FRM(KDIM3),DELT(KDIM3,KDIM3)
      DIMENSION ADL(KDIM3+1),SDELT(KDIM3+1,KDIM3+1)
      DIMENSION ADEL(KDIM3+2,KDIM3+2),BDL(KDIM3+2)
C
      CHARACTER*4 UNAME(4)
C
      DATA UNAME/'U(1)','U(2)','U(3)','U(4)'/
C     DATA UNAME/'U(1)','U(2)','U(3)','U(4)','U(5)','U(6)','U(7)'/
C
 1001 FORMAT(8E15.6)
 1002 FORMAT(/5H IEN=,I2,5X,5H JEN=,I2,5X,11H ENERGY WR=,F15.7,5H RYDS)
 1003 FORMAT(54H ****  DEBUGGING PRINT-OUT IN BASFUN FOR NBUG1=0  ****)
 1004 FORMAT(/   13H ERROR - IRX(,I1,27H)IS AN ODD NUMBER IN BASFUN/)
 1005 FORMAT(/8X,4H R  ,5(11X,A4)/)
 1006 FORMAT(/20H OUTWARD INTEGRATION)
 1007 FORMAT(/19H INWARD INTEGRATION)
 1008 FORMAT(/19H INWARD INTEGRATION,I2)
 1009 FORMAT(54H ****  DEBUGGING PRINT-OUT IN BASFUN FOR NBUG3=0  ****)
 1010 FORMAT(1H+,44X,2H /)
 1011 FORMAT(/  53H ERROR - NO CONVERGENCE IN BASFUN AFTER 99 ITERATIONS
     */)
 1012 FORMAT(/  55H ERROR - HINT,IHX,IRX AND RA ARE INCOMPATIBLE IN BASF
     *UN/)
 1013 FORMAT(20X,7F15.7)
 1014 FORMAT(/20H SOLUTION FROM MA01A/)
 1015 FORMAT(/17H ARRAYS FOR MA01A/)
 1016 FORMAT(/4H X1=,F15.8)
C
C      THE LARGEST AND SMALLEST FLOATING POINT NUMBERS ON THIS
C      MACHINE,  LARGE AND VSMALL, ARE SET IN A PARAMETER
C      STATEMENT.  QLARGE IS THE SQUARE ROOT OF LARGE.
C      THEY ARE USED IN CONSTRUCTING A SCALING FACTOR TO PREVENT
C      UNDERFLOW AND OVERFLOW  AT HIGH VALUES OF L.
C
C      CHECK COMPATIBILITY OF HINT,IHX,IRX AND RA
C
      S=ZERO
      J=0
      DO 1 I=1,NIX
      H=DBLE(IHX(I))*HINT
      S=S+H*DBLE(IRX(I)-J)
      J=IRX(I)
    1 CONTINUE
      IF(ABS(S-RA)-HALF*HINT) 3,3,2
    2 WRITE(IWRITE,1012)
      RETURN
C
C      CHECK THAT IRX(I) ARE EVEN INTEGERS
C
    3 DO 5 I=1,NIX
      J=IRX(I)
      M=MOD(J,2)
      IF(M) 4,5,4
    4 WRITE(IWRITE,1004)I
      RETURN
    5 CONTINUE
C
C      EVALUATE AND INITIALIZE SOME COMMONLY USED PARAMETERS
C
      ONEPT5=ONE+HALF
      IMATCH=IRX(NIX)-20
      IF(ETA.GT.ZERO) IMATCH=IRX(NIX)-10
C
C****  CHANGE TO ALLOW OUTWARD INTEGRATION TO REACH BOUNDARY
C      WILL HAPPEN AT HIGH L WHEN POINT OF INFLECTION IS OUTSIDE
C      RMATRIX BOUNDARY.    VMB JAN 1986
C
      IF (ETA .EQ. 0.0  .AND.  WINIT .GE. 0.0  .AND.  NBT .EQ. 0) THEN
         IMATCH=IRX(NIX)
      END IF
C
      IF(MOD(IMATCH,2).EQ.1) IMATCH=IMATCH+1
      FM=ZERO
      ITST=0
      NBTP1=NBT+1
      NG=NBT
      I9=IRX(NIX)+1
      DO 6 I=1,I9
      ORB(I)=ZERO
      DORB(I)=ZERO
    6 CONTINUE
      LANG=LC
      TLC=DBLE(LC*(LC+1))
      LCM=LC+1
      OVERLP=1.D0/DBLE(LCM)
      INST=1
      HIMT=HINT*DBLE(IHX(INST))
C      YD=HIMT**LCM
C      DYD=DBLE(LCM)*YD/HIMT
      IEN=1
      JEN=1
      WR=WINIT
      EIGEN=WR
C
C      BEGINNING OF LOOP TO FIND THE ENERGY WR, OR POTENTIAL MULTIPLE
C      PR, WHICH GIVES A FUNCTION SATISFYING THE LOGARITHMIC BOUNDARY
C      CONDITION.  THIS IS DONE BY USING NEWTONS METHOD TO FIND THE ZERO
C      OF THE FUNCTION BVAL-BSTO.  THIS LOOP IS LEFT ONLY WHEN THE
C      ENERGY, OR POTENTIAL MULTIPLE, INCREMENT DVAL IS LESS THAN ETA.
C
C      NOTE.... THIS LOOP IS ENTERED ONCE ONLY IF ETA=0.0
C
      DO 99 IEN=1,99
      IF(IEN.LT.99) GO TO 7
      WRITE(IWRITE,1011)
      RETURN
    7 WR1=WR
      DO 97 JEN=1,2
    8 IF(NBUG1.NE.0) THEN
      WRITE(IWRITE,1003)
      WRITE(IWRITE,1010)
      WRITE(IWRITE,1002)IEN,JEN,WR
      WRITE(IWRITE,1006)
      WRITE(IWRITE,1005)(UNAME(I),I=1,NBTP1)
      END IF
C
C
C  ****MODIFICATIONS FOR LARGE L
C
C      CALCULATE TWOZ
C
      TWOZ=0.0D0
      IF (POVALU(2) .GT. TINY) THEN
         NIZ=NINT(POVALU(1)/POVALU(2)-HALF)
         IF (NIZ .GE. 1) THEN
            TWOZ=POVALU(1)*HIMT*HALF
         END IF
      END IF
      ZSC=TWOZ*HALF
C
C      ESTIMATE THE POINT OF INFLECTION.  IF IT IS BEYOND
C      THE RADIUS AT WHICH R**(L+1) OVERFLOWS CALCULATE
C      SCALING FACTOR
C
      SCAST=1.0D0
      IF (ETA .EQ. 0.0D0  .AND.  WR .GT. 0.0D0  .AND.  NBT .EQ. 0
     1     .AND.  LCM .GT. 3) THEN
         IF (ZSC .NE. 0.0D0) THEN
            WRSC=WR/ZSC**2
         ELSE
            WRSC=WR
         END IF
         PTINF=(SQRT(1.D0+WRSC*TLC)-1.D0)/WRSC/ZSC
         RLIMIT=10**(LARGE*OVERLP)
         IF (PTINF .GE. RLIMIT) THEN
            SCAST1=RLIMIT/XR(I9)
            IF (SCAST1 .LT. 1.0) THEN
               SCAST=SCAST1**LCM
C              WRITE(IWRITE,'('' SCAST='',E14.5)') SCAST
            END IF
         END IF
      END IF
C
      XTEST=XR(2)
      NSTART=1
      IF (SCAST*XTEST**LCM .EQ. ZERO) THEN
         I2=-2
         DO 205  INTX=1,NIX
            I1=I2+5
            I2=IRX(INTX)-3
         DO 200 I=I1,I2
            IF (SCAST*XR(I)**LCM.GT. VSMALL)THEN
               NSTART=I-1
               INST=INTX
C      WRITE(IWRITE,'('' NSTART='',I4,''INST='',I4)') NSTART,INST
               GO TO 201
            END IF
 200     CONTINUE
 205     CONTINUE
         WRITE(IWRITE,'('' NO STARTING POSITION FOUND FOR L='',I4)')LANG
         STOP
      END IF
 201  CONTINUE
C
C      INITIALIZATION OF FUNCTION AT X=0.0 WHICH IS USED AS THE FIRST
C      POINT IN THE SIMPSONS RULE NORMALIZATION AND ORTHOGONALIZATION
C      INTEGRATIONS.
C*****INITIALISATION AT X=XR(NSTART) TO OVERCOME UNDERFLOW
C
    9 X=XR(NSTART)
      DO 10 II=1,NBTP1
      U(II) = ZERO
      DU(II) = ZERO
   10 CONTINUE
      IF(JEN.EQ.2) GO TO 12
      DO 11 I=1,NBTP1
C ****
      DO 202  III=1,NSTART
         US(I,III)=ZERO
         DUS(I,III)=ZERO
 202  CONTINUE
C
   11 CONTINUE
   12 IF(NBUG1.EQ.0) GO TO 13
      WRITE(IWRITE,1001)X,(U(NV),NV=1,NBTP1)
   13 CONTINUE
C
C      EVALUATE THE FUNCTION AND  DERIVATIVE AT HINT AND STORE THE
C      FUNCTION
C***** EVALUATE AT X=XR(NSTART+1)
C
      LP1=LCM
      LP2 = LC+2
C****
      HIMT=HINT*IHX(INST)
      X=XR(NSTART+1)
      H=HIMT
      DO 19 K=1,NBTP1
      IF(K.EQ.NBTP1) GO TO 17
C     U(K) = P(K,2)*H*H/DBLE(4*LC+6)
C***** CHANGE FOR START AT NSTART
C***** CAUTION  THIS HAS NOT BEEN TESTED BECAUSE EXPECT
C***** NO UNDERFLOW PROBLEM FOR L<LRANG1 I.E. WHEN THERE ARE
C***** BOUND TERMS
C
      U(K) = P(K,2*NSTART)*X*X/DBLE(4*LC+6)
      DU(K)=DBLE(LC+3)*U(K)*HALF
      GO TO 18
C *****
   17 U(K)=X**LP1*(ONE-TWOZ*X/(TWO*DBLE(LP1)))
      U(K)=U(K)*SCAST
      DU(K)=X**LC*(DBLE(LP1)-TWOZ*DBLE(LP2)*X/(TWO*DBLE(LP1)))
      DU(K)=DU(K)*SCAST
      DU(K)=DU(K)*H*HALF
   18 IF(JEN.NE.1) GO TO 19
C ****
      US(K,NSTART+1) = U(K)
      DUS(K,NSTART+1)=DU(K)*TWO/H
   19 CONTINUE
      IF(NBUG1.NE.0) WRITE(IWRITE,1001)X,(U(NV),NV=1,NBTP1)
C
C      EVALUATE FR AT HINT
C
C ****
      KM=2*NSTART-1
      MMM = 1
      CALL DERFUN
      DO 20 K=1,NBTP1
      FR(K)=FRH(K)*H*H/TWELVE
   20 CONTINUE
C
C      SET UP FRM AT HINT/2
C
      X=X-HALF*HIMT
C *****
      KM=2*NSTART-2
      IF(NBT.EQ.0) GO TO 22
      DO 21 K=1,NBT
      FRM(K)=U(K)
      U(K)=U(K)-DU(K)+ONEPT5*FR(K)
   21 CONTINUE
   22 FRM(NBTP1)=U(NBTP1)
      U(NBTP1)=X**LP1*(ONE-TWOZ*X/(TWO*DBLE(LP1)))
      U(NBTP1)=U(NBTP1)*SCAST
      CALL DERFUN
C ****
      KM=2*NSTART
      X=X+HALF*HIMT
      DO 23 K=1,NBTP1
      U(K)=FRM(K)
      FRM(K) = FRH(K)*H*H/TWELVE
   23 CONTINUE
C
C      STORE THE CONTRIBUTION TO THE INTEGRAL FROM THE FIRST POINT
C
      IF(NBT.EQ.0) GO TO 26
      DO 25 I=1,NBTP1
      DO 24 J=1,NBT
      DELT(I,J)=U(I)*P(J,KM)*H*FOUR
   24 CONTINUE
   25 CONTINUE
C
C      INTEGRATE OUT TO THE MATCHING POINT
C
   26 LSWT=1
C ****
      I1=NSTART+1
      DO 41 INTX=INST,NIX
      H1=H
      H=HINT*DBLE(IHX(INTX))
      HS = H*H
      I2=IRX(INTX)
      IF(IMATCH.LT.I2) I2=IMATCH
      IF(INTX.EQ.INST) GO TO 27
      LSWT = 2
      I1=IRX(INTX-1)+1
C
C      INTEGRATE OVER A RANGE OF EQUAL INTERVALS
C
   27 DO 40 IR=I1,I2
      JR = IR+1
      CALL DEVGL(NBTP1,DU,FR,FRM,H,HS,H1,LSWT)
      IF(NBUG1.EQ.1) WRITE(IWRITE,1001)X,(U(I),I=1,NBTP1)
   28 IF(NBT.EQ.0) GO TO 37
C
C      CALCULATE THE INTEGRATION FACTOR
C
      IF(IR.EQ.IMATCH) GO TO 33
      IF(IR.EQ.I2) GO TO 31
      IF(MOD(IR,2)) 29,30,29
   29 AM=FOUR
      GO TO 34
   30 AM=TWO
      GO TO 34
   31 IF(IR-IMATCH) 32,33,32
   32 AM = DBLE(IHX(INTX)+IHX(INTX+1))/DBLE(IHX(INTX))
      GO TO 34
   33 AM = ONE
C
C      ADD IN THE CONTRIBUTATION TO THE INTEGRAL FROM THE CURRENT POINT
C
   34 DO 36 I=1,NBTP1
      DO 35 J=1,NBT
      DELT(I,J)=DELT(I,J)+U(I)*P(J,KM)*H*AM
   35 CONTINUE
   36 CONTINUE
C
C      STORE THE FUNCTIONS AT EACH INTEGRATION
C
   37 IF(JEN.NE.1) GO TO 39
      DO 38 I=1,NBTP1
      US(I,JR)=U(I)
      DUS(I,JR)=DU(I)*TWO/H
   38 CONTINUE
C
   39 IF(IR.EQ.IMATCH) GO TO 42
   40 CONTINUE
   41 CONTINUE
C
C      STORE THE FUNCTIONS AND DERIVATIVES AT THE MATCHING POINT FOR THE
C      OUTWARD INTEGRATION
C
   42 DO 43 I=1,NBTP1
      YOUT(I) = U(I)
      DYOUT(I)=DU(I)*TWO/H
   43 CONTINUE
C****  VMB CHANGE JAN 1986
C
C     WRITE (IWRITE,'(''  IMATCH='',I3,'' X='',E16.5)')IMATCH,X
C
      IF (IMATCH .EQ. IRX(NIX)) THEN
C
C      NO MATCHING REQUIRED BUT FIRST CHECK IMATCH ALLOWED AT RA
C
         IF (NBT .NE. 0  .OR.  WINIT .LT. 0.0D0) THEN
            WRITE(IWRITE,'('' **** ERROR MATCH POINT ON THE BOUNDARY''/
     1                     '' NUMBER OF BOUND TERMS ='',I3,/
     2                     '' ENERGY ='',E16.5)') NBT,WINIT
            CALL EXIT (0)
          END IF
C
         ORBMAX=0.0D0
         DO 400  I=1,I9
            ORB(I)=US(1,I)
            ORBMAX=MAX(ORBMAX,ABS(ORB(I)))
            DORB(I)=DUS(1,I)
  400    CONTINUE
C
         BVALUE=DORB(I9)*RA/ORB(I9)
         GO TO 401
      END IF
C
C*****   END VMB CHANGE
C
C
C      INITIALIZE ARRAYS FOR DEVOGELAERE INTEGRATION INWARDS
C
      MMM = -1
      NBTM=NBT
      ITST=1
      NITS=1
      BNDRY(1)=ONE
      DBNDRY(1)=BSTO/RA
      IF(ABS(ETA).GT.TINY) GO TO 44
      BNDRY(1)=ONE
      DBNDRY(1)=ZERO
      BNDRY(2)=ZERO
      DBNDRY(2)=ONE
      NITS=2
C
C      LOOP OVER NUTTY .... NITS=1 FOR ITERATION ON THE MISMATCH TO FIND
C                           AN EIGENVALUE
C                           NITS=2 FOR THE SOLUTION AT A GIVEN ENERGY
C                           FOR WHICH TWO INDEPENDENT INWARD SOLUTIONS
C                           ARE NECESSARY TO OBTAIN CONTINUITY.
C
   44 DO 83 NUTTY=1,NITS
C
C      EVALUATE THE FUNCTION AND DERIVATIVES AT RA
C
      IF(NBUG1.EQ.0) GO TO 45
      IF(NITS.EQ.1) WRITE(IWRITE,1007)
      IF(NITS.EQ.2) WRITE(IWRITE,1008)NUTTY
      WRITE(IWRITE,1005)(UNAME(I),I=1,NBTP1)
   45 KM = 2*IRX(NIX)+1
      H=HINT*DBLE(IHX(NIX))
      HS = H*H
      X=RA
      IF(NUTTY.EQ.2) NBT=0
      NG=NBT
      NBTP1=NBT+1
      DO 47 I=1,NBTP1
      IF(I.EQ.NBTP1) GO TO 46
      U(I) = ZERO
      DU(I)=ZERO
      GO TO 47
   46 U(I)=BNDRY(NUTTY)
      DU(I)=-DBNDRY(NUTTY)*H*HALF
   47 CONTINUE
      IF(NBUG1.NE.0) WRITE(IWRITE,1001)X,(U(NV),NV=1,NBTP1)
C
C      EVALUATE FR AT RA
C
      CALL DERFUN
      KM=KM+1
      DO 48 K=1,NBTP1
      FR(K)=FRH(K)*HS/TWELVE
   48 CONTINUE
C
C      EVALUATE FRM AT RA+H/2
C
      X = RA+H*HALF
      DO 49 K=1,NBTP1
      FRH(K)=U(K)-DU(K)+ONEPT5*FR(K)
   49 CONTINUE
      FRHVAL=TLC/X/X-WR-ONEPT5*POVALU(KM-1)+HALF*POVALU(KM-2)
      DO 51 K=1,NBTP1
      FRM(K)=FRH(K)*FRHVAL
      IF(K.EQ.NBTP1) GO TO 50
      FRM(K)=FRM(K)+ONEPT5*P(K,KM-1)-HALF*P(K,KM-2)
   50 FRM(K)=FRM(K)*H*H*THIRD
   51 CONTINUE
      JR=I9
      X=RA
      KM = KM-1
C
C      ADD CONTRIBUTION TO ORTHOGONALITY INTEGRALS FROM FIRST POINT
C
      IF(NUTTY.EQ.1.OR.NBTM.EQ.0) GO TO 53
      DO 52 JL=1,NBTM
      ORTH(JL)=P(JL,KM)*U(1)*H
   52 CONTINUE
   53 IF(NBT.EQ.0) GO TO 56
      DO 55 I=1,NBTP1
      DO 54 J=1,NBT
      SDELT(I,J)=U(I)*P(J,KM)*H
   54 CONTINUE
   55 CONTINUE
C
C      STORE THE FUNCTION AT RA
C
   56 IF(JEN.NE.1) GO TO 59
      IF(NUTTY.EQ.2) GO TO 58
      DO 57 I=1,NBTP1
      US(I,JR) = U(I)
      DUS(I,JR)=DU(I)*TWO/H
   57 CONTINUE
      GO TO 59
   58 ORB(JR)=U(1)
      DORB(JR)=DU(1)*TWO/H
C
C      INTEGRATE IN TO THE MATCHING POINT
C
   59 LSWT = 1
      DO 80 INTX=1,NIX
      H1=H
      IMT = NIX-INTX+1
      H=-HINT*DBLE(IHX(IMT))
      HS = H*H
      IF(IMT-1) 60,61,60
   60 I1=IRX(IMT-1)
      GO TO 62
   61 I1=0
   62 I2=IRX(IMT)
      IF(IMATCH.GT.I1) I1=IMATCH
      I3=I1+1
      IF(INTX.EQ.1) GO TO 63
      LSWT = 2
C
C      INTEGRATE OVER A RANGE OF EQUAL INTEGRALS
C
   63 DO 79 IR=I3,I2
      JR = JR-1
      CALL DEVGL(NBTP1,DU,FR,FRM,H,HS,H1,LSWT)
      IF(NBUG1.EQ.1) WRITE(IWRITE,1001)X,(U(I),I=1,NBTP1)
C
C      CALCULATE THE INTEGRATION FACTOR
C
   64 IF(IR.EQ.I2) GO TO 67
      IF(MOD(IR,2)) 65,66,65
   65 AM=FOUR
      GO TO 70
   66 AM=TWO
      GO TO 70
   67 IF(I1-IMATCH) 68,69,68
   68 AM=DBLE(IHX(IMT)+IHX(IMT-1))/DBLE(IHX(IMT))
      GO TO 70
   69 AM = ONE
C
C      ADD IN THE CONTRIBUTION TO THE INTEGRAL FROM THE CURRENT POINT
C
   70 CONTINUE
   71 IF(NBT.EQ.0) GO TO 74
      DO 73 I=1,NBTP1
      DO 72 J=1,NBT
      SDELT(I,J)=SDELT(I,J)-U(I)*P(J,KM)*H*AM
   72 CONTINUE
   73 CONTINUE
C
C      STORE THE FUNCTIONS AT EACH INTERATION
C
   74 IF(JEN.NE.1) GO TO 78
      IF(NUTTY.EQ.2) GO TO 76
      DO 75 I=1,NBTP1
      US(I,JR)=U(I)
      DUS(I,JR)=DU(I)*TWO/H
   75 CONTINUE
      GO TO 78
   76 ORB(JR)=U(1)
      DORB(JR)=DU(1)*TWO/H
      IF(NBTM.EQ.0) GO TO 78
      DO 77 JL=1,NBTM
      ORTH(JL)=ORTH(JL)-U(1)*P(JL,KM)*H*AM
   77 CONTINUE
C
   78 IF(JR.EQ.IMATCH+1) GO TO 81
   79 CONTINUE
   80 CONTINUE
C
C      STORE THE FUNCTIONS AND DERIVATIVES AT THE MATCHING POINT FOR
C      THE INWARD INTEGRATION
C
   81 IF(NUTTY.EQ.2) GO TO 83
      DO 82 I=1,NBTP1
      YIN (I) = U(I)
      DYIN(I)=DU(I)*TWO/H
   82 CONTINUE
   83 CONTINUE
C
      NBT=NBTM
      NG=NBT
      NBTP1=NBT+1
      NBTP2=NBT+2
C
C      SET UP THE MATCHING EQUATIONS FOR THE ITERATION CASE
C
      DO 84 I=1,NBTP2
      BDL(I)=ZERO
   84 CONTINUE
      BDL(NBTP1) = ONE
      IF(NBT.EQ.0) GO TO 87
      DO 86 I=1,NBT
      DO 85 J=1,NBT
      ADEL(I,J)=DELT(J,I)+SDELT(J,I)
   85 CONTINUE
      ADEL(I,NBTP1) = DELT(NBTP1,I)
      ADEL(I,NBTP2) = SDELT(NBTP1,I)
   86 CONTINUE
   87 IF(ABS(ETA).LT.TINY) GO TO 100
      DO 88 I=1,NBTP1
      ADEL(NBTP1,I) = YOUT(I)
   88 CONTINUE
      ADEL(NBTP1,NBTP2) = ZERO
      IF(NBT.EQ.0) GO TO 90
      DO 89 I=1,NBT
      ADEL(NBTP2,I)=DYOUT(I)-DYIN(I)
   89 CONTINUE
   90 ADEL(NBTP2,NBTP1)=DYOUT(NBTP1)
      ADEL(NBTP2,NBTP2)=-DYIN(NBTP1)
      IF(NBUG3.EQ.0) GO TO 92
      WRITE(IWRITE,1015)
      DO 91 M=1,NBTP2
      WRITE(IWRITE,1013) (ADEL(M,N),N=1,NBTP2),BDL(M)
   91 CONTINUE
C
C      SOLVE THE MATCHING EQUATIONS
C
   92 CALL MA01A(ADEL,BDL,NBTP2,1,0,KDIM3+2,1)
      IF(NBUG1.NE.0) WRITE(IWRITE,1014)
      IF(NBUG1.NE.0) WRITE(IWRITE,1001)(BDL(J),J=1,NBTP2)
C
C      CALCULATE THE MISMATCH
C
      IF(JEN.NE.1) GO TO 94
      DO 93 NN=1,NBTP2
      ADL(NN)=BDL(NN)
   93 CONTINUE
   94 AJEN=ONE-BDL(NBTP2)*YIN(NBTP1)
      IF(NBT.EQ.0) GO TO 96
      DO 95 I=1,NBT
      AJEN  = AJEN-BDL(I)*YIN(I)
   95 CONTINUE
C
C      STORE THE MISMATCH,INCREMENT THE ENERGY AND RETURN
C
   96 AVAL(JEN)=AJEN
      WR = WR+DELTA
   97 CONTINUE
C
C      CARRY OUT ONE ITERATION USING NEWTONS METHOD AND RETURN
C
      FDASH = (AVAL(2)-AVAL(1))/DELTA
      IF(FDASH.EQ.ZERO) GO TO 106
      DVAL = -AVAL(1)/FDASH
      IF(ABS(DVAL)-ETA) 106,106,98
   98 WR = WR1+DVAL
   99 CONTINUE
C
      GO TO 106
C
C      SET UP THE MATCHING EQUATIONS FOR THE ARBITRARY ENERGY CASE
C
  100 NBTP3=NBTP2+1
      BDL(NBTP3)=ZERO
      IF(NBT.EQ.0) GO TO 102
      DO 101 I=1,NBT
      ADEL(I,NBTP3)=ORTH(I)
      ADEL(NBTP2,I)=YOUT(I)-YIN(I)
      ADEL(NBTP3,I)=DYOUT(I)-DYIN(I)
  101 CONTINUE
  102 DO 103 N=1,NBTP3
      ADEL(NBTP1,N)=ZERO
  103 CONTINUE
      ADEL(NBTP1,NBTP1)=ONE
      ADEL(NBTP2,NBTP1)=YOUT(NBTP1)
      ADEL(NBTP3,NBTP1)=DYOUT(NBTP1)
      ADEL(NBTP2,NBTP2)=-YIN(NBTP1)
      ADEL(NBTP2,NBTP3)=-U(1)
      ADEL(NBTP3,NBTP3)=-DU(1)*TWO/H
      ADEL(NBTP3,NBTP2)=-DYIN(NBTP1)
      IF(NBUG3.EQ.0) GO TO 105
      WRITE(IWRITE,1015)
      DO 104 M=1,NBTP3
      WRITE(IWRITE,1013) (ADEL(M,N),N=1,NBTP3),BDL(M)
  104 CONTINUE
  105 CALL MA01A(ADEL,BDL,NBTP3,1,0,KDIM3+2,1)
      NBTP2=NBTP3
      GO TO 108
C
C      STORE THE ENERGY EIGENVALUE AND FORM THE CONTINUOUS SOLUTION
C      IN ORB
C
  106 EIGEN=WR1
      DO 107 NN=1,NBTP2
      BDL(NN)=ADL(NN)
  107 CONTINUE
  108 IF(NBUG3.EQ.0) GO TO 109
      WRITE(IWRITE,1014)
      WRITE(IWRITE,1001)(BDL(J),J=1,NBTP2)
      WRITE(IWRITE,'(//'' MATCHING POINT='',I4//)')IMATCH
C
C      EVALUATE THE FINAL UNNORMALIZED FUNCTION AT THE MESH POINTS
C
  109 NBTN6=NBT+2
      ORBMAX=0.0D0
      DO 111 I=1,IMATCH
      ORB(I)=ZERO
      DORB(I)=ZERO
      DO 110 J=1,NBTP1
      ORB(I)=ORB(I)+US(J,I)*BDL(J)
C
C*****   CALCULATE THE MAXIMUM VALUE OF THE ORBITAL IN ORDER
C        TO SCALE TO PREVENT OVERFLOW DURING NORMLISATION
C
      ORBMAX=MAX(ORBMAX,ABS(ORB(I)))
      DORB(I)=DORB(I)+DUS(J,I)*BDL(J)
  110 CONTINUE
  111 CONTINUE
      I1=IMATCH+1
      DO 113 I=I1,I9
      ORB(I)=ORB(I)*BDL(NBTP2)+US(NBTP1,I)*BDL(NBTN6)
      ORBMAX=MAX(ORBMAX,ABS(ORB(I)))
      DORB(I)=DORB(I)*BDL(NBTP2)+DUS(NBTP1,I)*BDL(NBTN6)
      IF(NBT.EQ.0) GO TO 113
      DO 112 J=1,NBT
      ORB(I)=ORB(I)+US(J,I)*BDL(J)
      DORB(I)=DORB(I)+DUS(J,I)*BDL(J)
  112 CONTINUE
  113 CONTINUE
      BVALUE=BDL(NBTP2)*RA/ORB(I9)
C
C      NORMALIZE THE SOLUTION AND THE LAGRANGE MULTIPLIERS
C
C
C      TEST FOR OVERFLOW DURING NORMALIZATION AND SCALE THE
C      ORBITAL IF NECESSARY
C
  401 IF (ORBMAX .GT. QLARGE) THEN
         SCALEN=QLARGE/ORBMAX
      ELSE
         SCALEN=1.0D0
      END IF
C     WRITE(IWRITE,'(''   SCALEN='',E14.5)') SCALEN
C
      X1= ZERO
      I1=1
      DO 120 I=1,NIX
      H=HINT*DBLE(IHX(I))
      I2=IRX(I)+1
      IF(I.EQ.1) GO TO 114
      I1=IRX(I-1)+1
  114 DO 119 J=I1,I2
      IF(MOD(J,2)) 116,115,116
  115 AM=FOUR
      GO TO 118
  116 IF(J.EQ.I1.OR.J.EQ.I2) GO TO 117
      AM = TWO
      GO TO 118
  117 AM = ONE
  118 SCAORB=SCALEN*ORB(J)
      X1=X1+SCAORB*SCAORB*H*AM
  119 CONTINUE
  120 CONTINUE
C
      IF(NBUG3.NE.0) WRITE(IWRITE,1016) X1
      X2=SCALEN*SQRT(THREE/X1)
      NODES = 0
      X3 = X2
      IF(ORB(2).LT.ZERO) X3=-X2
C
C      EVALUATE NODES, THE NUMBER OF NODES IN THE FINAL FUNCTION
C
      ORB1=ORB(1)
      DO 125 I=2,I9
      ORB2=ORB(I)*X3
      ORB(I)=ORB2
      DORB(I)=DORB(I)*X3
      IF(ORB1) 121,124,122
  121 IF(ORB2) 124,123,123
  122 IF(ORB2) 123,123,124
  123 NODES=NODES+1
  124 ORB1=ORB2
  125 CONTINUE
      IF(NBT.EQ.0) GO TO 127
      DO 126 I=1,NBTP1
      ALAMDA(I)=BDL(I)*X3
  126 CONTINUE
C
  127 RETURN
      END
C=======================================================================
      SUBROUTINE BASORB (MINIM)
C
C      MODIFIED VERSION FOR USE WITH THE NEW FINDER ROUTINE.
C      WRITTEN AT MUENSTER.....AUGUST 1983
C


      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      LOGICAL RESTAR
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG,IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM3=IDIM9+1,KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9,
     *          KDIM10=IDIM7*2,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA
      COMMON /BNDORB/ P(IDIM9,KDIM10)
      COMMON /BUTT  / COEFF(3,IDIM2),EK2MAX,EK2MIN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI
     *                ,BVALUE
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5)
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS

1000  FORMAT(1X/TR28,' SUBROUTINE BASORB'/TR28,
     *                ' -----------------')
1001  FORMAT(/' AMPLITUDE AT RA',TR4,' EIGENVALUE',TR6,' NODES'/)
1002  FORMAT(/' AMPLITUDE AT RA',TR4,' EIGENVALUE',TR6,' NODES',TR6,
     *        ' PHASE'/)
1003  FORMAT(//' OVERLAP INTEGRALS BETWEEN THE ORBITALS FOR L =',I2/)
1004  FORMAT(TR1,10E11.2)
1005  FORMAT(TR3,F10.5,TR4,E15.7,TR4,I5,TR4,E15.7)
1006  FORMAT(//' ORBITALS FOR L =',I2)
1007  FORMAT(//' OVERLAP INTEGRALS BETWEEN THE N =',I2,' BOUND ORBITAL
     *AND EACH CONTINUUM ORBITAL SCHMIDT ORTHOGONALISATION')
1008  FORMAT(TR1,8E14.7)


C      LOOP OVER THE CONTINUUM ANGULAR MOMENTA

       WRITE(IWRITE,1000)
       NPTS=IRX(NIX)+1
       OPEN(UNIT=61,ACCESS='DIRECT',STATUS='SCRATCH',FORM='UNFORMATTED',
     *    RECL=MACDIM*NRANG2*NPTS)
       OPEN(UNIT=62,ACCESS='DIRECT',STATUS='SCRATCH',FORM='UNFORMATTED',
     *    RECL=MACDIM*(KDIM4*KDIM4+IDIM3*IDIM10))
C
C      INITIALISE CALCULATION OF THE ARRAY NVARY WHICH WILL CONTAIN
C      THE NUMBER OF CONTINUUM ORBITALS CALCULATED AT EACH VALUE OF
C      CONTINUUM ANGULAR MOMENTUM.  THE FINAL ORBITAL AT EACH L
C      WILL HAVE AN EIGENVALUE JUST GREATER THAN TWICE EK2MAX.
C      NVARY(L) WILL HAVE A MINIMUM VALUE OF 3.
C
       LPLO=MINIM
  1    DO 10,LP=LPLO,LRANG2
          NQ=IDIM13
          L1=LP-1
          MAXLG=MAXNLG(LP)
          MAXHF=MAXNHF(LP)
          ITEST=MAXHF-MAXLG
          NBT=0
          WRITE(IWRITE,1006) L1
          IF (LP.LE.LRANG1) THEN
C            CHECK IF LAGRANGE ORTHOGONALISATION IS REQUIRED
             NBT=MAXLG-LP+1
             IF (NBT.GT.0) THEN
                READ(60,REC=LP) P
             ENDIF
          ENDIF
          IF (L1.GT.2.OR.NELC.NE.NZ) THEN
             WRITE(IWRITE,1001)
          ELSE IF (L1.LE.2.AND.NELC.EQ.NZ) THEN
                  WRITE(IWRITE,1002)
          ENDIF

C         ALGORITHM FOR DETERMINING STARTING VALUES FOR FINDER

          IF (LP.EQ.LPLO) THEN
             IF (NELC.EQ.NZ) THEN
                ETRIAL=PI*PI/(RA*RA*FOUR)
             ELSE
                ETRIAL=DBLE((NELC-NZ))*PI/(RA+RA)
             ENDIF
             IF (LPLO .GT. LRANG1) THEN
                ETRIAL=0.1D0
             END IF
          ELSE
             ETRIAL=EIGENS(1,L1)
          ENDIF

C         LOOP TO GENERATE CONTINUUM ORBITALS

C
C      IF NRANG2 HAS NOT BEEN SET THEN CONTINUE ROUND LOOP UNTIL
C      THE EIGENVALUE IS TWICE THE READ IN EK2MAX
C
          DO 20,N=1,NRANG2
             IPOS(MAXHF+N,LP)=IDIM13+N
             NQ=NQ+1
             NODES=MAXHF-ITEST-LP+N
             RESTAR=.FALSE.
             CALL FINDER(NBT,L1,NODES,RESTAR,ETRIAL)

C            CALCULATE A NEW ENERGY ESTIMATE FOR THE NEXT ORBITAL

             IF (ETRIAL.LT.ZERO) THEN
                ETRIAL=ETRIAL+(PI*PI)/(RA*RA)
             ELSE
                ETRIAL=(SQRT(ETRIAL)+PI/RA)**2
             ENDIF

             EIGENS(N,LP)=EIGEN
C
C      IF L1=0 SET EK2MAX  FOR THE BUTTLE CORRECTION
C      EK2MIN IS SET ZERO IF THERE ARE NO NEGATIVE EIGENVALUES
C
             IF (N .EQ. 1) THEN
                IF (EIGEN .GE. 0.0) THEN
                   EK2MIN=0.0D0
                ELSE
                   EK2MIN=EIGEN-0.5D0
                END IF
             END IF
             IF((LP.EQ.1.OR.EK2MAX.LT.PT01).AND.N.EQ.NRANG2) THEN
                WRITE(IWRITE,'('' MAX. ENERGY READ AS EK2MAX WAS'',
     1                F14.7)') EK2MAX
                EK2MAX=EIGEN*0.5D0
                WRITE(IWRITE,'('' IT HAS BEEN RESET TO '',F14.7)')EK2MAX
             END IF
             DO 40,I=1,NPTS
                UJ(I,NQ)=ORB(I)
40           CONTINUE
             WRITE(IWRITE,1005) UJ(NPTS,NQ),EIGEN,NODES
C
C      TRUNCATE THE NUMBER OF CONTINUUM ORBITALS WHEN
C      THE EIGENVALUE REACHES TWICE EK2MAX
C
             IF (NVARY(LP) .EQ. 0  .AND.  N .GT. 2  .AND.
     1           EIGEN .GE. EK2MAX*2.0D0) THEN
                NVARY(LP)=N
                GO TO 24
             END IF
C
20        CONTINUE
C
C      CHECK THAT NVARY HAS BEEN SET
C
          IF (NVARY(LP) .EQ. 0) THEN
             WRITE (IWRITE,'(//'' ******* ERROR IN INPUT DATA'',
     1       /'' EK2MAX IS TOO LARGE FOR NRANG2'')')
             STOP
          END IF
C
C      SET NRANG2 TO THE CURRENT NUMBER OF CONTINUUM TERMS
C
 24       NRANG2=NVARY(LP)
C
         IF (NBUG7 .EQ. 1) THEN
          WRITE(IWRITE,1003) L1
          WRITE(IWRITE,'('' BEFORE ANY SCHMIDT OTHOGONALISATION'')')
          DO 21,N1=LP,MAXHF+NRANG2
             DO 22,N2=LP,N1
                N3=IPOS(N1,LP)
                N4=IPOS(N2,LP)
                RESULT=ZERO
                DO 23,I=1,NPTS
                   RESULT=RESULT+UJ(I,N3)*UJ(I,N4)*WT(I)
23              CONTINUE
                IF(N2 .EQ. N1)RESULT=RESULT-ONE
                TEMP(N2-LP+1)=RESULT
22           CONTINUE
             J=N1-LP+1
             WRITE(IWRITE,1004) (TEMP(I),I=1,J)
21        CONTINUE
         END IF

          IF (ITEST.GT.0) THEN
             DO 90,N1=1,ITEST
                N3=IPOS(MAXLG+N1,LP)
                WRITE(IWRITE,1007) MAXLG+N1
                DO 100,N=1,NRANG2
                   N4=IPOS(MAXHF+N,LP)
                   RESULT=ZERO
                   DO 110,I=1,NPTS
                      RESULT=RESULT+UJ(I,N3)*UJ(I,N4)*WT(I)
110                CONTINUE
                   OVRLAP(N,N1)=RESULT
100             CONTINUE
                WRITE(IWRITE,1008) (OVRLAP(N,N1),N=1,NRANG2)
90           CONTINUE
C
C      RESET NRANG2 TO AN ACCEPTABLY SMALLER VALUE IF PRESENTLY
C      TOO BIG FOR THE SCHMIDT ORTHOGONALISATION PROCEDURE
C      AND RESTART CALCULATION OF ORBITALS
C
             N5=LP
             CALL SCHMDT(N5)
             IF (N5 .LT. 0) THEN
                NRANG2=-N5-1
                DO 91  I=1,LRANG2
                   IF (NVARY(I) .EQ. 0) GO TO 1
                   NVARY(I)=NRANG2
 91             CONTINUE
                GO TO 1
             END IF
             WRITE(IWRITE,'(//'' SCHMIDT COEFFICIENTS''/)')
             DO 120  J=1,ITEST+NRANG2
                WRITE(IWRITE,1008)(B(J,I),I=1,ITEST+NRANG2)
120          CONTINUE
             WRITE(62,REC=LP) B,OVRLAP
          ENDIF
C
          DO 50,N=1,NRANG2
             ENDS(N,LP)=UJ(NPTS,IDIM13+N)
50        CONTINUE
C
          WRITE(61,REC=LP) ((UJ(I,IDIM13+N),I=1,NPTS),N=1,NRANG2)
C
          IF (ITEST .GT. 0  .AND.  NBUG7 .EQ. 1) THEN
          WRITE(IWRITE,1003) L1
          DO 60,N1=LP,MAXHF+NRANG2
             DO 70,N2=LP,N1
                N3=IPOS(N1,LP)
                N4=IPOS(N2,LP)
                RESULT=ZERO
                DO 80,I=1,NPTS
                   RESULT=RESULT+UJ(I,N3)*UJ(I,N4)*WT(I)
80              CONTINUE
                IF(N2 .EQ. N1)RESULT=RESULT-ONE
                TEMP(N2-LP+1)=RESULT
70           CONTINUE
             J=N1-LP+1
             WRITE(IWRITE,1004) (TEMP(I),I=1,J)
60        CONTINUE
          WRITE (IWRITE,'(/'' AMPLITUDE AT RA'',/3X,10F10.5)')
     1           (ENDS(N,LP),N=1,NRANG2)
          END IF
10     CONTINUE
       WRITE (IWRITE,'('' VARIATION OF NRANG2 WITH L''/
     1        ''   L  NRANG2'')')
C
       DO 130  L=LPLO,LRANG2
          LMINUS=L-1
          WRITE (IWRITE,'( 1X,I3,2X,I4)') LMINUS,NVARY(L)
  130     CONTINUE
C
       RETURN
       END
C=======================================================================
      BLOCK DATA
C
C      SET UP THE MOST COMMONLY USED REAL CONSTANTS IN /CONSTS/
C      THIS MAKES IT EASIER TO CHANGE THE PRECISION OF THE PROGRAM
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      COMMON/CONSTS/ZERO,ONE  ,PT01  ,PT001,PT0001,TINY  ,PI ,FSC  ,
     *              TWO ,THREE,FOUR  ,FIVE ,SIX   ,EIGHT ,TEN  ,TWELVE,
     *              HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH

      DATA ZERO,ONE  ,PT01  ,PT001 ,PT0001,TINY  ,PI  ,FSC        ,
     *     TWO ,THREE,FOUR  ,FIVE  ,SIX   ,EIGHT ,TEN   ,TWELVE,
     *     HALF,THIRD,FOURTH,FIFTH ,SIXTH ,EIGHTH,TENTH /
     *   0.0D0,1.0D0,1.0D-2,1.0D-3,1.0D-4,1.0D-6,3.14159D0,7.2973525D-3,
     *   2.0D0,3.0D0,4.0D0 ,5.0D0 ,6.0D0 ,8.0D0 ,10.0D0,12.0D0,
     *   0.5D0,0.333333333333333D0,0.25D0,0.2D0 ,0.166666666666667D0,
     *   0.125D0,0.1D0/
      END
C***********************************************************************
C
      SUBROUTINE BUTFIT(IMAX,E,F,RA,EMAX,ALPHA,BETA,NBUT,DELTA)
C
C
C  FITTING OF BUTTLE CORRECTIONS
C
C  IMAX = NUMBER OF POINTS FOR WHICH CORRECTION CALCULATED.
C  E(I) = ENERGY POINTS
C  F(I) = BUTTLE CORRECTION
C  RA = BOUNDARY RADIUS
C  EMAX = EIGENS(NPOLES,L)
C  ALPHA, BETA = FIT PARAMETERS
C  DELTA = ACCURACY ACHIEVED
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      LOGICAL POLE
      PARAMETER(IDIM3=MZNR2)
      PARAMETER(IDBUT=IDIM3+IDIM3)
      DIMENSION E(6),F(6),FKN(0:IDBUT),UK(0:IDBUT)
C
C  INITIALISATIONS
C
      RA2=RA**2
      DO 1 I=1,IMAX
    1 E(I)=RA2*E(I)
      NBUT=INT(.5D0+RA*SQRT(EMAX)/3.14159D0)
      G=-1.5707963D0
      DO 2 N=0,NBUT
      G=G+3.141592654D0
      FKN(N)=G
    2 UK(N)=G**2
      ALPHA=1.0D0
      BETA=0.0D0
C
C
C
C START ITERATIONS FOR FIT
C
      DO 1000 KK=1,5
C
      X11=0.0D0
      X12=0.0D0
      X22=0.0D0
      Y1=0.0D0
      Y2=0.0D0
      DELTA=0.0D0
C
C    START SUM OVER POINTS I
C
      DO 100 I=1,IMAX
      U=BETA+E(I)
C
C    CALCULATE FUNCTIONS B(U) AND C(U)
C
      B=0.0D0
      C=0.0D0
C    CASE OF U.GT.0.04
      IF(U.GT..04D0)THEN
        FK=SQRT(U)
        POLE=.FALSE.
        DO 10 N=0,NBUT
        IF(ABS(FK-FKN(N)).GT..3D0)THEN
          A=1.0D0/(U-UK(N))
          B=B+A
          C=C-A**2
        ELSE
          POLE=.TRUE.
          D1=FK-FKN(N)
        ENDIF
   10   CONTINUE
        IF(POLE)THEN
          D2=D1**2
          D=.33333333D0*D1*(1.D0+.066666667D0*D2*(1.D0+.0952381D0*D2))
          A=1.D0/(2.D0*FK-D1)
          BB=(D+A)/FK
          D=.33333333D0*(1.D0+D2*(.2D0+.031746032D0*D2))
          C=2.D0*C+.5D0*(D-A**2-BB)/U
          B=2.D0*B+BB
        ELSE
          T=TAN(FK)
          TK=T/FK
          B=2.D0*B+TK
          C=2.D0*C+.5D0*(1.D0+T**2-TK)/U
        ENDIF
C
C    SUM FOR U.LT..04
      ELSE
        DO 20 N=0,NBUT
        A=1.D0/(U-UK(N))
        B=B+A
   20   C=C-A**2
C
C    CASE OF U.LT..04 AND U.GT.-.04
        IF(U.GT.-.04D0)THEN
          B=2.D0*B+1.D0+.33333333D0*U*(1.D0+.4D0*U)
          C=2.D0*C+.33333333D0*(1.D0+U*(.8D0+.48571429D0*U))
C
C    CASE OF U.LT.-.04
        ELSE
          FK=SQRT(-U)
          T=TANH(FK)
          TK=T/FK
          B=2.D0*B+TK
          C=2.D0*C+.5D0*(1.D0-T**2-TK)/U
        ENDIF
C
      ENDIF
C
C    INCREMENT MATRICES X AND Y
      C=ALPHA*C
      DF=F(I)-ALPHA*B
      DD=ABS(DF)/F(I)
      IF(DELTA.LT.DD)DELTA=DD
      X11=X11+B*B
      X12=X12+B*C
      X22=X22+C*C
      Y1=Y1+DF*B
      Y2=Y2+DF*C
C
  100 CONTINUE
C
C    CHECK CONVERGENCE
C
      IF(DELTA.LT.1.D-4)RETURN
C
C  END ITERATIONS
C
C  SOLVE EQUATIONS AND INCREMENT ALPHA AND BETA
C
      DET=1.D0/(X11*X22-X12**2)
      ALPHA=ALPHA+DET*(X22*Y1-X12*Y2)
      BETA=BETA+DET*(-X12*Y1+X11*Y2)
C
 1000 CONTINUE
C
      RETURN
      END
C=====================================================================
      SUBROUTINE CORECT(N,LP,M,R0,SIGMA,C1)
C
C      TO CALCULATE THE PARAMETERS ASSOCIATED WITH THE GAUSSIAN
C      CORRECTION TO BE APPLIED TO THE BOUND ORBITAL SPECIFIED BY THE
C      QUANTUM NUMBERS (N,LP-1).
C      CORRECTION TO ORBITAL AT RADIUS R IS -C1*EXP(-((R-R0)/SIGMA)**2).
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM12=IDIM1*IDIM12)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM


C
C      DEFINITION OF S - THE GAUSSIAN IS REQUIRED TO FALL TO EXP(-S**2)
C      OF ITS MAXIMUM VALUE BY R=0.9*RA.
C
      S=TWO
C
C      CALCULATE THE FUNCTION AND ITS DERIVATIVE AT R=RA.
C
      PA=ZERO
      DP=ZERO
      DO 1 J=1,M
         TERM=C(LP,N,J)*RA**DBLE(IRAD(LP,N,J))*EXP(-ZE(LP,N,J)*RA)
         PA=PA+TERM
         DP=DP+(DBLE(IRAD(LP,N,J))/RA-ZE(LP,N,J))*TERM
1     CONTINUE
C
C      FIND PARAMETERS SIGMA, R0, C1 FROM THE LOGARITHMIC DERIVATIVE, B
C
      B=RA*DP/PA
      SIGMA=RA*(S-SQRT(S*S-FIFTH*B))/B
      R0=S*SIGMA+(ONE-TENTH)*RA
      C1=PA/EXP(-((RA-R0)/SIGMA)**2)
      RETURN
      END
C=======================================================================
      SUBROUTINE DERFUN
C
C      EVALUATES THE SECOND DERIVATIVE FUNCTION FOR THE DE VOGELAERE
C      ROUTINE DEVGL ASSOCIATED WITH THE NEW BASFUN
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG)
      PARAMETER(KDIM3=IDIM9+1,KDIM10=IDIM7*2)
C
      COMMON /BASDER/ FM,TLC,WR,IMEMM,ITST,JR,KM,MMM,NBTP1,NG
      COMMON /BNDORB/ P(IDIM9,KDIM10)
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /FUNVAL/ FRH(KDIM3),U(KDIM3),X
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /POTVAL/ POVALU(KDIM10)
C
      KM=KM+MMM
      FRHVAL=TLC/X/X-WR-POVALU(KM)
      IF(FRHVAL.LT.ZERO.AND.FM.GT.ZERO) GO TO 1
      IF(FRHVAL.GT.ZERO.AND.FM.LT.ZERO) GO TO 1
      FM=FRHVAL
      GO TO 2
    1 IF(ITST.EQ.1) GO TO 2
      IMATCH=JR
      IF(MOD(IMATCH,2).NE.0) IMATCH=IMATCH-1
      ITST=1
    2 FRH(NBTP1)=FRHVAL*U(NBTP1)
      IF(NG.EQ.0) RETURN
      DO 3 I=1,NG
      FRH(I)=FRHVAL*U(I)+P(I,KM)
    3 CONTINUE
C
      RETURN
      END
C=======================================================================
      SUBROUTINE DEVGL(M,DY,FR,FRM,H,HS,H1,LSWT)
C
C     DE VOGELAERE INTEGRATION ROUTINE
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM9=MZLAG)
      PARAMETER(KDIM3=IDIM9+1)
C
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /FUNVAL/ FRH(KDIM3),Y(KDIM3),X

      DIMENSION DY(KDIM3),FR(KDIM3),FRM(KDIM3),YR(KDIM3)
C
      HH=HALF*H
      X=X+HH
      GO TO (1,3),LSWT
    1 DO 2 I=1,M
      DY(I)=DY(I)+FR(I)
      YR(I)=Y(I)+DY(I)
      Y(I)=YR(I)+FR(I)-EIGHTH*FRM(I)
    2 CONTINUE
      GO TO 5
    3 H12=H/H1
      H12S=H12*H12
      LSWT=1
      DO 4 I=1,M
      DY(I)=H12*DY(I)+H12S*FR(I)
      YR(I)=Y(I)+DY(I)
      Y(I)=YR(I)+HALF*H12S*(FR(I)*(H12+ONE)-FOURTH*H12*FRM(I))
    4 CONTINUE
    5 CALL DERFUN
      DO 6 I=1,M
      FRM(I)=HS*FRH(I)/THREE
      DY(I)=DY(I)+FRM(I)
      Y(I)=YR(I)+DY(I)
    6 CONTINUE
      X=X+HH
      CALL DERFUN
      DO 7 I=1,M
      FR(I)=HS*FRH(I)/TWELVE
      DY(I)=DY(I)+FR(I)
    7 CONTINUE
C
      RETURN
      END
C=======================================================================
      SUBROUTINE EVALUE (MINIM)

C     THIS ROUTINE EVALUATES ALL THE BOUND ORBITALS AT THE MESH
C     POINTS


      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM10=IDIM7*2,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
C
      COMMON /BNDORB/ P(IDIM9,KDIM10)
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4
      COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS



C     OPEN THE DIRECT ACCESS SCRATCH FILE UNIT NUMBER 60 USED
C     TO STORE THE BOUND ORBITALS REQUIRED IN THE LAGRANGE
C     ORTHOGONALISATION PROCEDURE

      OPEN(UNIT=60,ACCESS='DIRECT',STATUS='SCRATCH',FORM='UNFORMATTED',
     *   RECL=MACDIM*IDIM9*KDIM10)

      NQ=0
      DO 10,L=1,LRANG1
         MAXHF=MAXNHF(L)
         MAXLG=MAXNLG(L)
         IF (MINIM .GT. L) MAXLG=0
         KT=0
         DO 20,N=L,MAXHF
C           CHECK IF THE BOUND ORBITALS ARE NUMERICAL
            IF (NUMERC.GT.0) THEN
C              THE BOUND ORBITALS ARE NUMERICAL. READ THEM IN FROM
C              THE  SEQUENTIAL FILE ITAPE1
               NQ=NQ+1
               KT=KT+1
               IPOS(N,L)=NQ
               NPTS=IRX(NIX)
               READ(ITAPE1)  UJ(1,NQ),(P(KT,I),I=1,2*NPTS)
               READ(ITAPE1) (UDP(NQ,I),I=1,NPTS+1)
               IF (L.EQ.1) THEN
                  PORI(N)=UJ(1,NQ)
                  DPORI(N)=UDP(NQ,1)
               ENDIF
               UJ(1,NQ)=ZERO
               UDP(NQ,1)=ZERO
               J=1
               DO 30,I=2,2*NPTS,2
                  J=J+1
                  UJ(J,NQ)=P(KT,I)
30             CONTINUE
               GOTO 20
            ENDIF
C           ANALYTIC BOUND ORBITALS ARE BEING USED
C           FIRST CALCULATE THE GAUSSIAN CORRECTION
            NQ=NQ+1
            M=NCO(NQ)
            IPOS(N,L)=NQ
            CALL CORECT(N,L,M,R0,SIGMA,C1)
            IF (N.LE.MAXLG) THEN
C              THIS ORBITAL IS REQUIRD FOR LAGRANGE
C              ORTHOGONALISATION
               JUP=2
               KT=KT+1
               MT=0
            ELSE
C              THIS ORBITAL IS NOT REQUIRED FOR LAGRANGE
C              ORTHOGONALISATION
               JUP=1
            ENDIF
            DO 40,I=2,IRX(NIX)+1
C              EVALUATE THE ORBITAL AT EACH MESH POINT
               H=STEP(I)*HALF
               DO 50,J=JUP,1,-1
                  R=XR(I)+(1-J)*H
                  PV=ZERO
                  ALPHA=(R-R0)/SIGMA
                  IF (-ALPHA.LT.TWELVE) PV=-C1*EXP(-ALPHA*ALPHA)
                  DO 60,K=1,M
                     PV=PV+C(L,N,K)*R**IRAD(L,N,K)*EXP(-ZE(L,N,K)*R)
60                CONTINUE

                  IF (J.EQ.1) THEN
                     UJ(I,NQ)=PV
                  ENDIF
                  IF (JUP.EQ.2) THEN
                     MT=MT+1
                     P(KT,MT)=PV
                  ENDIF
50             CONTINUE
40          CONTINUE
            UJ(1,NQ)=ZERO

C           RENORMALISE THE ORBITAL
            AN=ZERO
            DO 100,I=2,IRX(NIX)+1
               AN=AN+UJ(I,NQ)*UJ(I,NQ)*WT(I)
100         CONTINUE
            AN=SQRT(AN)
            DO 70,I=2,IRX(NIX)+1
               UJ(I,NQ)=UJ(I,NQ)/AN
70          CONTINUE

            DO 90,I=1,M
               C(L,N,I)=C(L,N,I)/AN
90          CONTINUE

               IF (N.LE.MAXLG) THEN
               DO 80,I=1,2*IRX(NIX)
                  P(KT,I)=P(KT,I)/AN
80             CONTINUE
            ENDIF
20       CONTINUE

C        WRITE THE ORBITALS TO THE DIRECT ACCESS SCRATCH FILE
C        UNIT NUMBER 60.
         WRITE(60,REC=L) P

10    CONTINUE
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE FINDER(NBT,LC,NODES,RESTAR,ETRIAL)

C     A NEW FINDER ROUTINE WRITTEN AT MUENSTER...AUGUST 1983

C     NBT..... THE NUMBER OF BOUND ORBITALS TO ORTHOGONALISE TO.
C     LC ..... THE L VALUE.
C     NODES... THE REQUIRED NUMBER OF NODES IN THE EIGENFUNCTION.
C     RESTAR.. LOGICAL VARIABLE - IF TRUE THEN THIS IS A RESTART RUN.
C                               - IF FALSE THIS IS NOT A RESTART RUN.
C     ETRIAL.. INITIAL ESTIMATE OF THE EIGENENERGY.



      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      LOGICAL RESTAR
      PARAMETER (ERROR=1.0D-6)
C
      PARAMETER(IDIM7=MZPTS,IDIM9=MZLAG)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM3=IDIM9+1)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI
     *                ,BVALUE

C
C     FORMAT STATEMENTS
C
100   FORMAT(1X,'  DEBUG FROM NEW FINDER'/
     *          '   =====================')
101   FORMAT(' NODE= ',I5,' NODES= ',I5,' BVALUE= ',E14.6,
     *       ' EIGEN= ',E14.6,' ETRIAL= ',E14.6)
102   FORMAT(1X,'  STOP IN NEW FINDER.......CANNOT FIND THE CORRECT NUMB
     *ER OF NODES')
103   FORMAT(1X,' ERROR IN ROOT...IFLAG = ',I5)
104   FORMAT(1X,' ERROR IN NEW FINDER...BASFUN ITERATION ENTERED BUT INC
     *ORRECT NUMBER OF NODES FOUND')

C
C     INITIALISE VARIABLES.
C
      IF (NBUG4.GT.0) WRITE(IWRITE,100)
      DEL=PI/RA
      DEL2=DEL*DEL
      DELL=DEL
      ISUM=0
10    CONTINUE
C
      CALL BASFUN(NBT,LC,NODE,RA,BSTO,ETRIAL,ZERO,ZERO)
C
      IF (NBUG4.GT.0) WRITE(IWRITE,101) NODE,NODES,BVALUE,EIGEN,ETRIAL
C
C     FIRST CHECK WHETHER THIS IS A RESTART RUN
C
      IF (RESTAR) RETURN
C
C     CHECK THAT THE FUNCTION HAS THE CORRECT NUMBER OF NODES.
C
      IF (NODE.NE.NODES) THEN
C
C       FUNCTION HAS NOT THE CORRECT NUMBER OF NODES,MODIFY ENERGY
C       ACCORDINGLY.
C
        IF (ISUM.GT.0) THEN
          IF ((NODE-NODES).EQ.NHOLD) THEN
             DELL=DELL*HALF
             DEL2=DEL2*HALF
             ETRIAL=EHOLD
          ELSE
             NHOLD=NODES-NODE
          ENDIF
        ELSE
          NHOLD=NODES-NODE
        ENDIF
C
        EHOLD=ETRIAL
        ISUM = ISUM+1
        REALNN=DBLE(NODES-NODE)
        IF ((ETRIAL.LT.ZERO).OR.((SQRT(ABS(ETRIAL))+REALNN*DELL).
     1                        LT.ZERO)) THEN
          ETRIAL=ETRIAL+ABS(REALNN)*REALNN*DEL2
        ELSE
           ETRIAL=(REALNN*DELL+SQRT(ABS(ETRIAL)))**2
        ENDIF
        IF (ISUM.EQ.100) THEN
          WRITE(IWRITE,102)
          STOP
        ELSE
          GOTO 10
        ENDIF
C
      ELSE
C
C        FUNCTION HAS THE CORRECT NUMBER OF NODES.
C
        IF (BVALUE.LT.BSTO) THEN
C
C         WE HAVE AN UPPER BOUND TO THE ENERGY,NOW FIND LOWER BOUND.
C
          EHIGH=ETRIAL
          BLOW=BVALUE
          DEL1=DEL/THREE
          DEL2=DEL2/THREE
          EHOLD=ETRIAL
C
20        CONTINUE
C
          IF ((ETRIAL.LT.ZERO).OR.((SQRT(ABS(ETRIAL))-DEL1).LT.ZERO))
     1        THEN
             ETRIAL=ETRIAL-DEL2
          ELSE
             ETRIAL = (SQRT(ABS(ETRIAL))-DEL1)**2
          ENDIF
C
          CALL BASFUN(NBT,LC,NODE,RA,BSTO,ETRIAL,ZERO,ZERO)
C
          IF(NBUG4.GT.0)WRITE(IWRITE,101)NODE,NODES,BVALUE,EIGEN,ETRIAL
C
          IF (NODE.NE.NODES) THEN
C
C           ENERGY DECREASED TOO FAR,DECREASE INCREMENT
C
            DEL1=DEL1/THREE
            DEL2=DEL2/THREE
            ETRIAL=EHOLD
            GOTO 20
C
          ELSE
C
            IF (BVALUE.LT.BSTO) THEN
C
C             BETTER UPPER BOUND TO THE ENERGY FOUND,TRY AGAIN
C             FOR A LOWER BOUND.
C
              EHIGH=ETRIAL
              BLOW=BVALUE
              EHOLD=ETRIAL
              GOTO 20
C
            ELSE
C
C             LOWER BOUND TO THE ENERGY NOW FOUND.
C
              ELOW=ETRIAL
              BHIGH=BVALUE
            ENDIF
C
          ENDIF
C
        ELSE
C
C         LOWER BOUND TO THE ENERGY FOUND,NOW TRY FOR AN UPPER BOUND.
C
          ELOW=ETRIAL
          BHIGH=BVALUE
          DEL1=DEL/THREE
          DEL2=DEL2/THREE
          EHOLD=ETRIAL
C
30        CONTINUE
C
          IF (ETRIAL.LT.ZERO) THEN
            ETRIAL=ETRIAL+DEL2
          ELSE
            ETRIAL = (SQRT(ABS(ETRIAL))+DEL1)**2
          ENDIF
C
          CALL BASFUN(NBT,LC,NODE,RA,BSTO,ETRIAL,ZERO,ZERO)
C
          IF(NBUG4.GT.0)WRITE(IWRITE,101)NODE,NODES,BVALUE,EIGEN,ETRIAL
C
          IF (NODE.NE.NODES) THEN
C
C           ENERGY INCREASED TOO FAR,DECREASE INCREMENT.
C
            DEL1=DEL1/THREE
            DEL2=DEL2/THREE
            ETRIAL=EHOLD
            GOTO 30
          ELSE
            IF (BVALUE.GE.BSTO) THEN
C
C             BETTER LOWER BOUND TO THE ENERGY FOUND,TRY AGAIN
C             FOR THE UPPER BOUND.
C
              ELOW=ETRIAL
              BHIGH=BVALUE
              EHOLD=ETRIAL
              GOTO 30
C
            ELSE
C
C             UPPER BOUND TO THE ENERGY FOUND.
C
              EHIGH=ETRIAL
              BLOW=BVALUE
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C
C     WE NOW HAVE UPPER AND LOWER BOUNDS TO THE EIGENENERGY,NOW
C     USE ROOT FINDING ROUTINE TO GET A BETTER ESTIMATE.
C
      IFLAG=1
      ABSERR=ZERO
      RELERR=ERROR*TENTH
      B=ELOW
      C=EHIGH
      IF (ETRIAL.EQ.ELOW) THEN
        FT = BHIGH - BSTO
      ELSE
        FT = BLOW - BSTO
      ENDIF
C
   40 CONTINUE
C
      CALL ROOT(T,FT,B,C,RELERR,ABSERR,IFLAG)
C
      IF(IFLAG.GT.1)THEN
        IF(ABS(ELOW-EHIGH).LT.ERROR/TENTH)THEN
          IFLAG=1
          B=HALF*(ELOW+EHIGH)
        ENDIF
      ENDIF
C
      IF (IFLAG.LT.0) THEN
C
C     STILL NOT CLOSE ENOUGH TO THE ROOT
C
        IF (T.EQ.ELOW) THEN
          FT=BHIGH-BSTO
          GOTO 40
C
        ELSE IF (T.EQ.EHIGH) THEN
          FT=BLOW-BSTO
          GOTO 40
C
        ELSE
          ETRIAL=T
C
          CALL BASFUN(NBT,LC,NODE,RA,BSTO,ETRIAL,ZERO,ZERO)
C
          IF(NBUG4.GT.0)WRITE(IWRITE,101)NODE,NODES,BVALUE,EIGEN,ETRIAL
C
          FT=BVALUE-BSTO
C
C         IF ABS(FT) IS LESS THAN ERROR..SOLUTION FOUND
C
          IF (ABS(FT).LT.ERROR)GOTO 50
          GOTO 40
C
        ENDIF
C
      ELSE IF (IFLAG.GT.1) THEN
C
C       ERROR IN ROOT
C
        WRITE(IWRITE,103) IFLAG
        STOP
C
      ELSE
C
C       IFLAG=1...ROOT SUCESSFULLY LOCATED TO REQUIRED ACCURACY
C
        ETRIAL=B
C
        CALL BASFUN(NBT,LC,NODE,RA,BSTO,ETRIAL,ZERO,ZERO)
C
        IF (NBUG4.GT.0)WRITE(IWRITE,101) NODE,NODES,BVALUE,EIGEN,ETRIAL
        GOTO 50
C
      ENDIF
50    IF (NODE.EQ.NODES) THEN
C
C       ALL IS WELL...EIGENVALUE FOUND WITH CORRECT NUMBER
C       OF NODES
C
        RETURN
      ELSE
C
C       INCORRECT NUMBER OF NODES.....GIVE UP
C
        WRITE(IWRITE,104)
        STOP 'INCORRECT NUMBER OF NODES'
      ENDIF
C     RETURN
      END
C=======================================================================
      SUBROUTINE GEN1CC (MINIM)
C
C      GENERATES AND STORES ALL CONTINUUM-CONTINUUM ONE ELECTRON
C      INTEGRALS
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)

      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9)

      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1,
     *                JBUFF1,JBUFF2
      DIMENSION X(3)
 1000 FORMAT ('      L   N1   N2  X(I)')
 1001 FORMAT (2X, 3I5, 3E13.5)
C
      IF (NBUG6 .GT. 0) THEN
         WRITE (IWRITE,1000)
      END IF
C

      DO 40,L=MINIM,LRANG2
         MAXHF=MAXNHF(L)
         CALL SKIPER(L,0)
         IST3(L)=IBUFF1+1+(LBUFF1*NBUFF1)
         IF (IBUFF1+1.GT.LBUFF1) THEN
            PRINT*,' ERROR IN CODE STOP IN  GEN1CC'
            STOP
         ENDIF
         DO 10,N1=1,NVARY(L)
            DO 20,N2=1,N1
               CALL ONEELE(N1+MAXHF,L,N2+MAXHF,L,X(1))
               DO 30,I=1,1
                  IBUFF1=IBUFF1+1
                  BUFF1(IBUFF1)=X(I)
                  IF (IBUFF1.EQ.LBUFF1) THEN
                     WRITE(JBUFF1,REC=NBUFF1) BUFF1
                     NBUFF1=NBUFF1+1
                     IBUFF1=0
                  ENDIF
30             CONTINUE
            IF (NBUG6 .GT. 0) THEN
               WRITE (IWRITE,1001) L,N1,N2,(X(I),I=1,1)
            END IF
20          CONTINUE
10       CONTINUE
40    CONTINUE
      RETURN
      END
C==============================================================
      SUBROUTINE GENCC(L,LP,NRKVAL)
C
C      GENERATES AND STORES ALL THE CONTINUUM-CONTINUUM RK INTEGRALS
C      WITH CONTINUUM ANGULAR MOMENTA L AND LP
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,
     *  IDIM7=MZPTS,IDIM8=MZAMP,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)

      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /CORE  / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2),
     *                 MAXNC(IDIM2),MAXPN(IDIM2),ICHECK
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1,
     *                JBUFF1,JBUFF2
C
 1000 FORMAT ('    LAM   L1  L1P   N1  N1P    N   NP  RKVAL')
 1001 FORMAT (2X, 7I5, E13.5)
C
C      ZEROIZE VARIABLES AND ARRAYS
C
      NRKVAL=0
      I1=MIN(2*LRANG1-1,L+LP+1)
C      I2=MIN(LRANG1+L,LRANG1+LP)
      DO 4 I=1,LRANG1
         DO 3 J=1,LRANG1
            DO 1 K=1,I1
               ICTCCD(I,J,K)=-999
    1       CONTINUE
    3    CONTINUE
    4 CONTINUE
      MAXHF=MAXNHF(L+1)
      MAXHFP=MAXNHF(LP+1)

C      FIRST CALCULATE THE DIRECT INTEGRALS

      IF (NBUG6 .EQ. 2) THEN
         WRITE (IWRITE,1000)
      END IF
C
      LAMUP=MIN(L+LP,LRANG1+LRANG1-2,LLPDIM-1)
      DO 10,LAM=ABS(L-LP),LAMUP,2
         DO 20,L1=0,LRANG1-1
            MAXC1=MAXNC(L1+1)
            MAXHF1=MAXNHF(L1+1)
            IF(MAXC1.EQ.MAXHF1) GOTO 20
            LLO=ABS(LAM-L1)
            IF (LLO.LT.L1) THEN
               LLO=2*((L1-LLO+1)/2)+LLO
            ENDIF
            L1PUP=MIN(LRANG1-1,LAM+L1)
            DO 30,L1P=LLO,L1PUP,2
               MAXC2=MAXNC(L1P+1)
               MAXHF2=MAXNHF(L1P+1)
               IF(MAXC2.EQ.MAXHF2) GOTO 30
               ICTCCD(L1+1,L1P+1,LAM+1)=IBUFF1+1+(LBUFF1*NBUFF1)
               IF (IBUFF1+1.GT.LBUFF1) THEN
                  PRINT*,' ERROR IN CODE STOP IN GENCC DIRECT'
                  STOP
               ENDIF
               DO 40,N1=MAXC1+1,MAXHF1
                  IF (L1P.EQ.L1) THEN
                     N1PUP=N1
                  ELSE
                     N1PUP=MAXHF2
                  ENDIF
                  DO 50,N1P=MAXC2+1,N1PUP
                     DO 60,N=1,NVARY(L+1)
                        IF (L.EQ.LP) THEN
                           NPUP=N
                        ELSE
                           NPUP=NVARY(LP+1)
                        ENDIF
                        DO 70,NP=1,NPUP
                           CALL RS(N1,L1+1,MAXHF+N,L+1,N1P,
     *                             L1P+1,MAXHFP+NP,LP+1,LAM,
     *                             RKVAL)
                           IBUFF1=IBUFF1+1
                           BUFF1(IBUFF1)=RKVAL
                           NRKVAL=NRKVAL+1
                           IF (IBUFF1.EQ.LBUFF1) THEN
                              WRITE(JBUFF1,REC=NBUFF1) BUFF1
                              NBUFF1=NBUFF1+1
                              IBUFF1=0
                           ENDIF
                        IF (NBUG6 .EQ. 2) THEN
                           WRITE (IWRITE,1001) LAM,L1,L1P,
     1                           N1,N1P,N,NP,RKVAL
                        END IF
70                      CONTINUE
60                   CONTINUE
50                CONTINUE
40             CONTINUE
30          CONTINUE
20       CONTINUE
10    CONTINUE

      RETURN
      END
C=======================================================================
      SUBROUTINE GENINT (MINIM)
C
C     THIS ROUTINE CONTROLS THE EVALUATION OF ALL THE RADIAL INTEGRALS
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM4=MZLMX,
     *  IDIM7=MZPTS,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(KDIM6=(IDIM4+1)/2,KDIM9=IDIM1+IDIM1-1,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9)
C
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON /YKSTOR/ YK(IDIM7),RK(IDIM7),RK1(IDIM7),JN2,JL2,JN4,JL4,JKM
      COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1,
     *                JBUFF1,JBUFF2
C
      DIMENSION IBBPOL(IDIM1,IDIM1,KDIM6)
C
 1000 FORMAT(' '/28X,' SUBROUTINE GENINT'/28X,' -----------------')
 1001 FORMAT(' GENERATE THE MULTIPOLE INTEGRALS')
 1002 FORMAT(' GENERATE THE BOUND-BOUND MULTIPOLE INTEGRALS')
 1003 FORMAT(' GENERATE THE BOUND-CONTINUUM MULTIPOLE INTEGRALS')
 1004 FORMAT(' GENERATE THE CONTINUUM-CONTINUUM MULTIPOLE INTEGRALS')
 1006 FORMAT(' GENERATE THE BOUND-BOUND ONE ELECTRON INTEGRALS')
 1008 FORMAT(' GENERATE THE BOUND-CONTINUUM ONE ELECTRON INTEGRALS')
 1010 FORMAT(' GENERATE CONTINUUM-CONTINUUM ONE ELECTRON INTEGRALS')
 1012 FORMAT(' GENERATE THE BOUND-BOUND RK INTEGRALS')
 1014 FORMAT(' GENERATE THE BOUND-CONTINUUM RK INTEGRALS')
 1015 FORMAT(' L =',I5,' LP =',I5)
 1016 FORMAT(' GENERATE THE CONTINUUM-CONTINUUM RK INTEGRALS')

      JKM=999
      IBUFF1=0
      NBUFF1=2
      NBUFF2=2
C
      WRITE(IWRITE,1000)
C

C     GENERATE THE MULTIPOLE RADIAL INTEGRALS

      WRITE(IWRITE,1001)
      LAMIND=(LAMAX+1)/2

      IF (LAMAX.GT.0) THEN
         WRITE(IWRITE,1002)
      ENDIF
      CALL GENMBB(NBBPOL,IBBPOL)
C
      WRITE(IWRITE,1010)
      CALL GEN1CC (MINIM)

C     WRITE THE ONE ELECTRON POINTER ARRAYS TO RAD3
C
      WRITE(ITAPE3)(((IBBPOL(I,J,K),I=1,LRANG1),J=1,LRANG1),
     1                     K=1,LAMIND),NBBPOL,(IST3(I),I=1,LRANG2)
C
C     GENERATE THE TWO ELECTRON RADIAL INTEGRALS

      WRITE(IWRITE,1016)
C
C      PREPARE FOR STORAGE OF RK POINTERS
C
      LAMDIM=LRANG1+LRANG1-1
      NBUF=(MINIM-1)*LLPDIM+1
      DO 10,L=MINIM-1,LRANG2-1
         LPUP=MIN (LRANG2-1, LLPDIM-1+L)
         LDIF=0
          NBUFF2=NBUF
         DO 20,LP=L,LPUP
            LDIF=LDIF+1
C
            WRITE(IWRITE,1015) L,LP
            CALL SKIPER(L+1,LP+1)
            CALL GENCC(L,LP,NRKVAL)
C
            WRITE(JBUFF2,REC=NBUFF2) (((ICTCCD(I,J,K),I=1,LRANG1),
     1                               J=1,LRANG1),K=1,LAMDIM)
            NCCD(LDIF,L+1)=NRKVAL
C
            NBUFF2=NBUFF2+1
20       CONTINUE
         NBUF=NBUF+LLPDIM
10    CONTINUE


C
C      WRITE UNFINISHED BLOCK TO DISC
C
      IF (IBUFF1.GT.0) THEN
         DO 30  I=IBUFF1+1,LBUFF1    !IBUFF1 -> IBUFF1+1
            BUFF1(I)=ZERO
  30     CONTINUE
         WRITE(JBUFF1,REC=NBUFF1) BUFF1
      ENDIF

C
C      WRITE NO EXCHANGE POINTERS TO RAD3
C
         WRITE(ITAPE3) ((NCCD(LDIF,L),LDIF=1,LLPDIM),L=1,LRANG2)
      RETURN
      END
C=======================================================================
      SUBROUTINE GENMBB(NBBPOL,IBBPOL)
C
C      GENERATES ALL THE BOUND-BOUND MULTIPOLE INTEGRALS
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM4=MZLMX,
     *  IDIM7=MZPTS,IDIM8=MZAMP,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM6=(IDIM4+1)/2,KDIM9=IDIM1+IDIM1-1,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9)
C
      COMMON /CORE  / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2),
     *                 MAXNC(IDIM2),MAXPN(IDIM2),ICHECK
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1,
     *                JBUFF1,JBUFF2
      DIMENSION X(2)
      DIMENSION IBBPOL(IDIM1,IDIM1,KDIM6)
C
 1000 FORMAT ('    LAM   L1  L1P   N1  N1P  RVAL')
 1001 FORMAT (2X, 5I5, E13.5)
C
      IF (NBUG6 .GT. 0) THEN
         WRITE (IWRITE,1000)
      END IF
C
C
C      ZEROIZE ARRAYS
C
      DO 3,I=1,LRANG1
         DO 2,J=1,LRANG1
            DO 1,K=1,LAMIND
               IBBPOL(I,J,K)=-999
1           CONTINUE
2        CONTINUE
3     CONTINUE

      NBBPOL=0
      LUP=LRANG1-1
      DO 10,L1=0,LUP
         MAXC1=MAXNC(L1+1)
         MAXHF1=MAXNHF(L1+1)
         IF (MAXC1.EQ.MAXHF1) GOTO 10
         DO 20,L2=L1,LUP
            MAXC2=MAXNC(L2+1)
            MAXHF2=MAXNHF(L2+1)
            IF (MAXC2.EQ.MAXHF2) GOTO 20
            LAMLO=ABS(L1-L2)
            IF (LAMLO.EQ.0) THEN
               LAMLO=2
            ENDIF
            LAMST=(LAMLO+1)/2-1
            LAMUP=MIN(L1+L2,LAMAX)
            DO 30,LAM=LAMLO,LAMUP,2
               LAMST=LAMST+1
               IBBPOL(L1+1,L2+1,LAMST)=IBUFF1+1+(LBUFF1*NBUFF1)
               IF (IBUFF1+1.GT.LBUFF1) THEN
                  PRINT*,' ERROR IN CODE STOP IN GENMBB'
                  STOP
               ENDIF
               DO 40,N1=MAXC1+1,MAXHF1
                  IF (L1.EQ.L2) THEN
                     N2UP=N1
                  ELSE
                     N2UP=MAXHF2
                  ENDIF
                  DO 50,N2=MAXC2+1,N2UP
                     CALL RADINT(N1,L1+1,N2,L2+1,LAM,X(1))
                        IUP=1
                     DO 60,I=1,IUP
                        NBBPOL=NBBPOL+1
                        IBUFF1=IBUFF1+1
                        BUFF1(IBUFF1)=X(I)
                        IF (IBUFF1.EQ.LBUFF1) THEN
                           WRITE(JBUFF1,REC=NBUFF1) BUFF1
                           NBUFF1=NBUFF1+1
                           IBUFF1=0
                        ENDIF
C     WRITE (IWRITE,1001) LAM,L1,L2,N1,N2,X(I)
                     IF (NBUG6 .GT. 0) THEN
                        WRITE (IWRITE,1001) LAM,L1,L2,N1,N2,X(I)
                     END IF
60                   CONTINUE
50                CONTINUE
40             CONTINUE
30          CONTINUE
20       CONTINUE
10    CONTINUE
      RETURN
      END
C
C**********************************************************
      SUBROUTINE LSQ(P,Q,A,IMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C
      DIMENSION P(6),Q(6),X(3,3),C(3),A(3)
C
C
  600 FORMAT(/18H COEFFICIENTS ARE ,3F14.7)
  601 FORMAT(' X(I,J) = ',F14.7)
  602 FORMAT(/'C(J) = ',F14.7)
C
C
      N=IMAX
C
      DO 100 I=1,3
         DO 101 J=1,3
            X(I,J)=0.0D0
            L=I+J-2
            IF(L.NE.0) GOTO 105
            X(I,J)=N
            GOTO 101
  105       DO 102 K=1,N
               X(I,J)=X(I,J)+P(K)**L
  102       CONTINUE
C           WRITE(6,601) X(I,J)
  101    CONTINUE
  100 CONTINUE
C
      DO 103 J=1,3
         C(J)=0.0D0
         L=J-1
         IF(L.NE.0) GOTO 106
         DO 107 K=1,N
            C(J)=C(J)+Q(K)
  107    CONTINUE
         GOTO 103
  106    DO 104 K=1,N
            C(J)=C(J)+(P(K)**L)*Q(K)
  104    CONTINUE
C        WRITE(6,602) C(J)
  103 CONTINUE
C
      CALL MA01A(X,B,3,0,1,3,1)
      CALL PMUL(X,C,A,1,3)
C
C
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE MA01A(A,B,M,N,M1,IA,IB)
C
C      SOLUTION OF SIMULTANEOUS EQUATIONS AND OR MATRIX INVERSION
C
C
C      A           THE M*M MATRIX OF LEFT HAND SIDES OR THE MATRIX BEING
C                  INVERTED. OVERWRITTEN ON EXIT BY THE INVERSE MATRIX
C
C      B           THE M*N MATRIX OF THE RIGHT HAND SIDES. OVERWRITTEN
C                  ON EXIT BY SOLUTIONS
C
C      M           THE ORDER OF THE A-MATRIX. THIS MUST NOT BE GREATER
C                  THAN THE DIMENSION OF THE C AND IND ARRAYS. THE UPPER
C                  LIMIT CAN BE EXTENDED BY RECOMPILING WITH LARGER
C                  DIMENSIONS FOR THE PRIVATE ARRAYS C AND IND
C
C      N           THE NUMBER OF THE RIGHT HAND SIDES IN THE
C                  SIMULTANEOUS EQUATIONS
C
C      M1          =0 ONLY SIMULTANEOUS EQUATIONS ARE SOLVED IF N.GT.0
C                       IF N=0 A FURTHER ENTRY TO MA01A WITH M1.LT.0
C                       REQUIRED TO OBTAIN THE INVERSE OF A
C                  .GT.0 MATRIX INVERSION. IN ADDITION SIMULTANEOUS
C                       EQUATIONS ARE SOLVED IF N.GT.0
C                  .LT.0 ONLY USED IF PREVIOUS ENTRY TO MA01A
C                       WITH M1=0. IN THIS CASE THE MATRIX INVERSION IS
C                       COMPLETED
C
C      IA          DEFINES THE DIMENSIONS OF THE ARRAY WHERE
C                  THE A-MATRIX IS STORED
C
C      IB          DEFINES THE SECOND DIMENSION OF THE ARRAY WHERE
C                  B-MATRIX IS STORED
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM9=MZLAG)
      PARAMETER(KDIM3=IDIM9+1)
C
      DIMENSION A(IA,IA),B(IA,IB),C(KDIM3+2),IND(KDIM3+2)
C
      ZERO=DBLE(0)
      ONE=DBLE(1)
      MM=M-1
      IF(M1) 32,1,1
C
C      CHECK FOR ZERO DIAGONAL ELEMENTS
C
    1 IF(MM) 55,2,2
    2 DO 3 I=1,M
      IF(A(I,I).EQ.ZERO) A(I,I)=DBLE(10)**(-35)
    3 CONTINUE
C
C      DO THE TRIVIAL CASE OF A 1*1 MATRIX
C
      IF(MM) 55,4,7
    4 IF(N-1) 5,6,7
    5 A(1,1)=ONE/A(1,1)
      GO TO 55
    6 B(1,1)=B(1,1)/A(1,1)
      GO TO 5
C
C      THIS IS NOT A TRIVIAL CASE
C
C      FIND THE FIRST PIVOTAL ELEMENT AND STORE THE CORRESPONDING ROW
C      NUMBER IN I4. IND DEFINES THE ORDER OF THE ROWS OF THE ORIGINAL
C      A-MATRIX BEFORE ROW INTERCHANGE
C
    7 AMAX=ZERO
      DO 9 I=1,M
      IND(I)=I
      IF(ABS(A(I,1))-AMAX) 9,9,8
    8 AMAX=ABS(A(I,1))
      I4=I
    9 CONTINUE
C
C      EACH TIME THROUGH THE FOLLOWING LOOP THE A-MATRIX IS
C      REDUCED BY ONE
C
      IF(MM) 24,24,10
   10 DO 23 J=1,MM
C
C      INTERCHANGE THE I4TH AND THE JTH ROWS OF THE A-MATRIX AND STORE
C      ORDER IN IND IF I4 .NE.J
C
      IF(I4-J) 15,15,11
   11 ISTO=IND(J)
      IND(J)=IND(I4)
      IND(I4)=ISTO
      DO 12 K=1,M
      STO=A(I4,K)
      A(I4,K)=A(J,K)
      A(J,K)=STO
   12 CONTINUE
C
C      INTERCHANGE THE I4TH AND THE JTH ROWS OF THE B-MATRIX IF N.GT. 0
C
      IF(N) 15,15,13
   13 DO 14 K=1,N
      STO=B(I4,K)
      B(I4,K)=B(J,K)
      B(J,K)=STO
   14 CONTINUE
C
C      THE JTH ROW NOW CONTAINS THE PIVOTAL ELEMENT IN THE JTH POSITION
C      ELIMINATE THE JTH ELEMENT FROM EACH OF THE REMAINING ROWS OF THE
C      A-MATRIX AND THE B-MATRIX AND STORE THE MULTIPLIERS IN THE LOWER
C      TRIANGLE
C
   15 AMAX=ZERO
      J1=J+1
      DO 22 I=J1,M
      A(I,J)=A(I,J)/A(J,J)
      DO 18 K=J1,M
      A(I,K)=A(I,K)-A(I,J)*A(J,K)
      IF(K-J1) 16,16,18
C
C      FIND THE NEXT PIVOTAL ELEMENT AND STORE THE CORRESPONDING ROW
C      NUMBER IN I4
C
   16 IF(ABS(A(I,K))-AMAX) 18,18,17
   17 AMAX=ABS(A(I,K))
      I4=I
   18 CONTINUE
   19 IF(N) 22,22,20
   20 DO 21 K=1,N
      B(I,K)=B(I,K)-A(I,J)*B(J,K)
   21 CONTINUE
   22 CONTINUE
   23 CONTINUE
C
C      THE ELIMINATION IS NOW COMPLETE AND THE A-MATRIX HAS BEEN
C      REDUCED TO THE PRODUCT OF AN UPPER AND LOWER TRIANGLE MATRIX
C
   24 IF(N) 31,31,25
C
C      NOW CARRY OUT THE BACK SUBSTITUTION AND STORE RESULT IN THE
C      B-MATRIX IF THERE IS AT LEAST ONE RIGHT HAND SIDE
C
   25 DO 30 I1=1,M
      I=M+1-I1
      DO 29 J=1,N
      IF(M-I) 28,28,26
   26 I2=I+1
      DO 27 K=I2,M
      B(I,J)=B(I,J)-A(I,K)*B(K,J)
   27 CONTINUE
   28 B(I,J)=B(I,J)/A(I,I)
   29 CONTINUE
   30 CONTINUE
   31 IF(M1) 55,55,32
C
C      REPLACE THE A-MATRIX BY ITS INVERSE WHEN M1.NE. ZERO
C
C      FIRST INVERT THE LOWER TRIANGLE MATRIX AND STORE ON ITSELF
C
   32 IF(MM) 40,40,33
   33 DO 39 I1=1,MM
      I=M+1-I1
      I2=I-1
      DO 37 J1=1,I2
      J=I2+1-J1
      J2=J+1
      W1=-A(I,J)
      IF(I2-J2) 36,34,34
   34 DO 35 K=J2,I2
      W1=W1-A(K,J)*C(K)
   35 CONTINUE
   36 C(J)=W1
   37 CONTINUE
      DO 38 K=1,I2
      A(I,K)=C(K)
   38 CONTINUE
   39 CONTINUE
C
C      NOW INVERT THE UPPER TRIANGLE MATRIX AND FORM THE ORIGINAL
C      A-MATRIX APART FROM COLUMN INTERCHANGE.  THIS OVERWRITES THE
C      ORIGINAL A-MATRIX
C
   40 DO 50 I1=1,M
      I=M+1-I1
      I2=I+1
      W=ONE/A(I,I)
      DO 48 J=1,M
      IF(I-J) 41,42,43
   41 W1=ZERO
      GO TO 44
   42 W1=ONE
      GO TO 44
   43 W1=A(I,J)
   44 IF(I1-1) 47,47,45
   45 DO 46 K=I2,M
      W1=W1-A(I,K)*A(K,J)
   46 CONTINUE
   47 C(J)=W1
   48 CONTINUE
      DO 49 J=1,M
      A(I,J)=C(J)*W
   49 CONTINUE
   50 CONTINUE
C
C      RE-ORDER THE COLUMNS OF THE INVERSE A-MATRIX TO COINCIDE WITH
C      THE ORDER OF THE ROWS OF THE A-MATRIX ON INPUT
C
      DO 54 I=1,M
   51 IF(IND(I)-I) 52,54,52
   52 J=IND(I)
      DO 53 K=1,M
      STO=A(K,I)
      A(K,I)=A(K,J)
      A(K,J)=STO
   53 CONTINUE
      ISTO=IND(J)
      IND(J)=J
      IND(I)=ISTO
      GO TO 51
   54 CONTINUE
C
   55 RETURN
      END
C=======================================================================
      SUBROUTINE MESH
C
C      NEW VERSION WRITTEN BY VAL BURKE 7 OCTOBER 1985
C
C      THIS ROUTINE AUTOMATICALLY GENERATES THE INTEGRATION MESH,
C      BASED ON THE NUCLEAR BEHAVIOUR OF BOUND ORBITALS, THE NUMBER
C      OF CONTINUUM ORBITAL LOOPS AND THE CURRENT ARRAY DIMENSIONS
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C
      PARAMETER(MSHDIM=16)
C
      PARAMETER(IDIM1=MZLR1,
     *  IDIM6=MZNIX,IDIM7=MZPTS,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM12=IDIM1*IDIM12)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM

C
      DATA MINFAC/1/,MAXFAC/2/
C
      NCORSE=NRANG2*MSHDIM
      IF (NCORSE .LT. 96) NCORSE=96
C
C      CALCULATE THE COARSEST MESH, HMAX, AND THE MESH REQUIRED
C      NEAR THE ORIGIN. CALCULATE HMIN = HMAX/2**M WHERE M+1 IS
C      THE NUMBER OF STEP SIZES
C
      HMAX=RA/NCORSE
      HINNER=.025D0/NZ
      DELTA=HINNER/5.0D0
      M=0
      HMIN=HMAX
1     IF (HMIN .GT. HINNER+DELTA) THEN
         M=M+1
         HMIN=HMIN/2
         GO TO 1
      END IF
C
      HINT=HMIN
      NIX=M+1
      IF (NIX .GT. IDIM6) CALL RECOV1(6,NIX,IDIM6)
C
C      SET UP THE IHX ARRAY
C
      IH=1
      DO 2   I=1,NIX
         IHX(I)=IH
         IH=IH+IH
    2 CONTINUE
C
C      CONSIDER SEPARATELY M .LT. 4  AND  M .GE. 4
C      NA IS THE NUMBER OF STEPS AT EACH STEP SIZE
C      IT IS A MULTIPLE OF 16 AND CAN TAKE VALUES
C      FROM 16*MAXFAC DOWN TO 16*MINFAC
C
      NAFAC=MAXFAC
C
      IF (M .GE. 4) THEN
    3    NA=16*NAFAC
         NTOT=NCORSE+(M-1)*NA+NA/8
         IF (NTOT .GE. IDIM7) THEN
            NAFAC=NAFAC-1
            IF (NAFAC .LT. MINFAC) CALL RECOV1(7,NTOT,IDIM7)
            GO TO 3
         END IF
C
C      SET UP IRX ARRAY
C
         IRX(1)=NA
         IRX(2)=NA+IRX(1)
         IRX(3)=NA+IRX(2)
         IRX(4)=NA+NA/8+IRX(3)
         DO 4   I=5,M
            IRX(I)=NA+IRX(I-1)
    4    CONTINUE
C
      ELSE
C
         MPOW2=2**M
    5    NA=16*NAFAC
         NTOT=NCORSE+M*NA-NA*(MPOW2-1)/MPOW2
         IF (NTOT .GE. IDIM7) THEN
            NAFAC=NAFAC-1
            IF (NAFAC .LT. MINFAC) CALL RECOV1(7,NTOT,IDIM7)
            GO TO 5
         END IF
C
C      FILL IRX ARRAY
C
         IA=0
         DO 6    I=1,M
            IA=IA+NA
            IRX(I)=IA
    6    CONTINUE
C
      END IF
C
C      NUMBER OF STEPS AT EACH STEP SIZE MUST BE EVEN
C
      IF (MOD(NTOT,2) .NE. 0) THEN
         NTOT=NTOT+1
         IRX(M)=IRX(M)+2
      END IF
C
      IRX(NIX)=NTOT
C
C      PERFORM CHECK
C
      RVAL=ZERO
      IR=0
      DO 7  I=1,M+1
         RVAL=RVAL+HINT*(IRX(I)-IR)*IHX(I)
         IR=IRX(I)
    7 CONTINUE
CW    WRITE (IWRITE,'('' RVAL='',E14.7)')RVAL
C
      RETURN
      END
C*********************************************************************
      SUBROUTINE NEWBUT(MINIM,MJS)
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM3=IDIM9+1,KDIM10=IDIM7*2,
     *          KDIM12=IDIM1*IDIM12)
C
      LOGICAL OK
C
      COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA
      COMMON/BASIC/ BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     1              LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /BNDORB/ P(IDIM9,KDIM10)
      COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN
      COMMON/INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI
     *                ,BVALUE
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
C
      DIMENSION EN(IDIM3,IDIM2),AMP(IDIM3,IDIM2)
      DIMENSION EL(6),RT(6),RCN(6),A(3),EK(6)
C
 1000 FORMAT('  BUTTLE CORRECTION COEFFICIENTS FOR EACH CHANNEL L'/
     1       '     L',12X,'A(0)',16X,'A(1)',16X,'A(2)'/)
 1001 FORMAT('  ',I4,3F20.6)
C
 6000 FORMAT(' TWO BUTTLE POLES ARE EQUAL-CONTINUE WITHOUT CORRECTION')
 6001 FORMAT(' SUBROUTINE NEWBUT - ENERGY RANGE TOO SMALL TO AVOID BUTTL
     *E POLE, L = ',I3)
 6003 FORMAT(' '/28X,' SUBROUTINE NEWBUT'/28X,' -----------------'/
     1       ' EK2MAX = ',F14.7,' EK2MIN = ',F14.7)
 6004 FORMAT(' SUBROUTINE NEWBUT - CANNOT EVALUATE BUTTLE CORRECTION FOR
     * L = ',I3,' AND ENERGY ',F14.7,'RYD.')
C
C LOOP OVER CHANNEL ANGULAR MOMENTA TO CALCULATE BUTTLE POLES
C UP TO LRANG2
C
C      IPTS=2*IRX(NIX)
      I9=IRX(NIX)+1
C
      WRITE(IWRITE,6003) EK2MAX,EK2MIN
         WRITE(IWRITE,1000)
C
C      SET COEFF(3,1) TO INDICATE TO THE ASYMPTOTIC PROGRAM
C      WHICH FITTING PROCEDURE IS USED
C
      IF (MINIM .GT. 1) THEN
         IF (MJS .EQ. 1) THEN
            COEFF(3,1)=-20000.0D0
         ELSE
            COEFF(3,1)=0.0D0
         END IF
      END IF
C
      DO 50 L=MINIM,LRANG2
         NPOLES=NVARY(L)
C
C SET UP ENERGY POINTS IN EK()
C
C     WRITE(IWRITE,6005)
      IMAX=6
      EKINT=IMAX-1
      EKINT=(EK2MAX-EK2MIN)/EKINT
      E=EK2MIN-EKINT
      DO 99 I=1,IMAX
         E=E+EKINT
         EK(I)=E
   99 CONTINUE
C
         L1=L-1
         NBT=0
         IF(L.LE.LRANG1) NBT=MAXNLG(L)-L+1
         IF (NBT .GT. 0) READ(60,REC=L) P
         ITEST=0
         IF(L.GT.LRANG1) GOTO 51
         ITEST=MAXNHF(L)-MAXNLG(L)
   51    IF (ITEST .GT. 0) GO TO 52
C
C    IF MAXNHF(L).EQ.MAXNLG(L) THE BUTTLE POLES AND AMPLITUDES ARE EQUAL
C    TO THE CONTINUUM EIGENVALUES AND FUNCTIONS EVALUATED IN STG1 AND
C    STORED IN EIGENS AND ENDS.
C
C
         DO 53 N=1,NPOLES
            AMP(N,L)=ENDS(N,L)
            EN(N,L)=EIGENS(N,L)
   53    CONTINUE
         GO TO 60
C
C    IF MAXNLG(L).LT.MAXNHF(L) EVALUATE BUTTLE POLES
C
   52    ITEST=0
         OK=.TRUE.
         DO 54 N=1,NPOLES
            NC=N-ITEST
            IF(NC) 55,55,56
   55       W=EN(N,L)
            EH=0.0D0
            GOTO 57
   56       W=EIGENS(NC,L)
            EH=0.0D0
   57       CALL BASFUN(NBT,L1,NODE,RA,BSTO,W,DELTA,EH)
            EN(N,L)=EIGEN
            AMP(N,L)=ORB(I9)
            N1=N-1
            IF(N1) 54,54,58
   58       DO 59 NLN=1,N1
               IF(ABS(EIGEN-EN(NLN,L)).LT.TINY) OK=.FALSE.
   59       CONTINUE
   54    CONTINUE
         IF (OK) GO TO 60
         WRITE(IWRITE,6000)
         RETURN
C
C CHECK THAT THE CHANNEL ENERGIES DO NOT COINCIDE WITH THE
C BUTTLE POLES AND THEN CALCULATE COEFFICIENTS FOR QUADRATIC FIT.
C
C
C    SET CHANNEL ENERGIES AND TEST AGAINST POLES
C
C
C      FROM OPACITY VERSION OF STGR
C
C      ENSURE THAT EL(I) NOT TOO CLOSE TO POLES.  EXCLUDE RANGE
C      ((1-X)*E(J)+X*E(J-1)) TO ((1-X)*E(J)+X*E(J+1))
C      WHERE E(J)=EN(J,L) AND  X=0.2
C
   60 DO 63 I=1,IMAX
         EL(I)=EK(I)
   63 CONTINUE
      X=0.2D0
      NA=1
      DO 35  I=1,IMAX
C
C      FIND FIRST POLE ABOVE EL(I)
C
      DO 20  N=NA,NPOLES
         IF (EN(N,L) .GT. EL(I)) THEN
            NN=N
            GO TO 30
         END IF
   20 CONTINUE
   30 NA=NN
      EA=EN(NA,L)
C
C      CASE OF ALL POLES ABOVE
C
      IF (NA .EQ. 1) THEN
         EP=EN(2,L)
         D=X*(EP-EA)
         IF ((EA-EL(I)) .LT. D) EL(I)=EA-D
C
C      SOME POLES BELOW
C
      ELSE
C
C      NEAREST POLE BELOW
C
         EB=EN(NA-1,L)
         D=X*(EA-EB)
         IF ((EL(I)-EB) .LT. D) EL(I)=EB+D
         IF ((EA-EL(I)) .LT. D) EL(I)=EA-D
      END IF
   35 CONTINUE
C
C      CHECK NUMBER OF INDEPENDENT VALUES
C
      INDEP=1
      DO 36  N=2,IMAX
         IF (EL(N-1) .NE. EL(N)) INDEP=INDEP+1
   36 CONTINUE
      IF (INDEP .LT. 3) THEN
         WRITE (IWRITE,6001) L1
         STOP
      END IF
C
C
C    ARRAY EL() NOW CONTAINS ACCEPTABLE ENERGIES FOR ANGULAR MOMENTUM
C    L . NOW EVALUATE THE BUTTLE CORRECTION AT THESE ENERGIES.
C
   67    DO 68 II=1,IMAX
            IF(ABS(EL(II)).LT.TINY)EL(II)=EL(II)+TWO*TINY
            EH=ZERO
            CALL BASFUN(NBT,L1,NODE,RA,BSTO,EL(II),DELTA,EH)
            RT(II)=BVALUE-BSTO
            IF(ABS(RT(II)).GT.TINY) GO TO 692
            WRITE(IWRITE,6004) L1,EL(II)
C            IFLAG=1
            STOP
  692       SUM=ZERO
            DO 69 N=1,NPOLES
               IF(AMP(N,L).EQ.0.0D0) GO TO 69
               SUM=SUM+AMP(N,L)*AMP(N,L)/(EN(N,L)-EL(II))
   69       CONTINUE
            SUM=SUM/RA
            R=ONE/RT(II)
            RCN(II)=R-SUM
   68    CONTINUE
C
C      FIT BUTTLE CORRECTION USING MJS METHOD IF MJS IS NON ZERO
C      UNLESS BSTO IS ZERO
C      USE FIT TO QUADRATIC IF MJS =0 OR BSTO IS ZERO
C
      IF (BSTO .NE. ZERO  .OR.  MJS .EQ. 0) THEN
C
         CALL LSQ(EL,RCN,A,IMAX)
         DO 691 I=1,3
            COEFF(I,L)=A(I)
  691    CONTINUE
C
           DELTB=0.0D0
           DO 693 I=1,IMAX
           DD=ABS(RCN(I)-A(1)-EL(I)*(A(2)+EL(I)*A(3)))/RCN(I)
  693      IF(DELTB.LT.DD)DELTB=DD
C
      ELSE
C
C  FOR BSTO.NE.0. USE IMPROVED FITTING PROCEDURE
        CALL BUTFIT(IMAX,EL,RCN,RA,EIGENS(NPOLES,L),ALPHA,BETA,
     +              NBUT,DELTB)
        COEFF(1,L)=ALPHA
        COEFF(2,L)=BETA
        COEFF(3,L)=-10000*NBUT
C
      ENDIF
C
         WRITE(IWRITE,1001) L1,(COEFF(I,L),I=1,3)
C
C  ACCURACY OF FIT
      WRITE(6,6010)DELTB
 6010 FORMAT(25X,'LARGEST ERROR IN FIT =',E15.6)
C
C
   50 CONTINUE
      RETURN
      END
C=======================================================================
      SUBROUTINE ONEELE(N11,L1,N12,L2,ALBVAL)
C
C      CALCULATE THE ONE-ELECTRON MATRIX ELEMENT BETWEEN ORBITALS
C      DEFINED BY THE QUANTUM NUMBERS (N1,L1-1) AND (N2,L2-1),
C      AND STORE THE RESULT IN ALBVAL.
C      WHEN BOTH ORBITALS ARE CONTINUUM ORBITALS THE SCHMIDT
C      COEFFICIENTS ARE USED TO EXPRESS ONE CONTINUUM ORBITAL IN TERMS
C      OF ORBITALS SATISFYING A DIFFERENTIAL EQUATION AND BOUND ORBITALS
C      WHEN ONE OR BOTH ORBITALS ARE BOUND THEN
C      ANALYTIC DIFFERENTIATION OF THE BOUND ORBITAL IS USED.
C      SEE COMMENTS IN EVALUE AND CORECT FOR DETAILS OF THE
C      CORRECTION APPLIED TO THE BOUND ORBITALS.
C NRB:
C      ADDED MASS-VELOCITY CONTRIBUTION IN THE HIGH-L LIMIT I.E.
C      NO BOUND-BOUND OR BOUND-CONTINUUM, JUST CONTINUUM-CONTINUUM
C      WITH NO LAGRANGE OR SCHMIDT CONTRIBUTION. VALID FOR A CONTINUATION
C      RUN WHERE (NORMALLY) THE MIN CONTINUUM L EXCEEDS THE MAX BOUND L.
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)

      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9,KDIM10=IDIM7*2,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
C
      COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA
      COMMON /COPY  / ICOPY1,ICOPY2,ITOTAL,IPSEUD
      COMMON /CORE  / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2),
     *                 MAXNC(IDIM2),MAXPN(IDIM2),ICHECK
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /POTVAL/ POVALU(KDIM10)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5)
      COMMON /NRBREL/IRELMV
C
C      RETURN IF THE ANGULAR MOMENTA ARE NOT COMPATIBLE
C
      N1=N11
      N2=N12
      ALBVAL = ZERO
      IF(L1.NE.L2) RETURN
C
C      DEFINE THE POSITION OF THE MODEL POTENTIAL IN THE POTHAM ARRAY.
C
      LLL=LPOSX(L1)
    1 ZN  =DBLE(NZ+NZ)
      ALSQ=DBLE(L1*(L1-1))
C
C      FIND IF ONE ORBITAL IS BOUND. IF SO PUT IT INTO THE SECOND
C      POSITION SO THAT THE ORBITAL DEFINED BY (N2,L2) IS BOUND.
C
      MAXHF = L1-1
      IF(L1.LE.LRANG1) MAXHF = MAXNHF(L1)
      IF(N2.LE.MAXHF) GO TO 2
      IF(N1.GT.MAXHF) GO TO 11
      N3=N1
      N1=N2
      N2=N3
    2 I1=IPOS(N1,L1)
      I2=IPOS(N2,L2)
      IF (NUMERC.EQ.0) THEN
         M=0
         DO 40,I=1,L2-1
            DO 50,K=I,MAXNHF(I)
               M=M+1
50          CONTINUE
40       CONTINUE
         M2=NCO(M+N2-L2+1)
         CALL CORECT(N2,L2,M2,R0,SIGMA,C1)
      ENDIF
      IST = 1
      R = ZERO
C
C      CARRY OUT THE INTEGRATION USING SIMPSONS RULE
C
      DO 10 I=1,NIX
      H =HINT*DBLE(IHX(I))
      MI=1
      SUM = ZERO
      IFI = IRX(I)
      DO 9 JP=IST,IFI
      R = R+H
C
C      DETERMINE THE CONTRIBUTION INVOLVING THE ANGULAR MOMENTUM
C      AND THE NUCLEAR CHARGE
C
      RINV = ONE/R
      J1=JP+1
      IF(IPSEUD.GT.0) POVAL=POTHAM(LLL,JP)
      IF(IPSEUD.LT.0)POVAL=TWO*POTHAM(LLL,JP)*RINV
      IF(IPSEUD.EQ.0) POVAL=ZN*RINV
      X1=UJ(J1,I2)*(POVAL-ALSQ*RINV*RINV)
C
C      DETERMINE THE CONTRIBUTION INVOLVING THE SECOND DERIVATIVE
C      ACTING ON THE BOUND-ORBITAL
C
      IF (NUMERC.EQ.0) THEN
C
C       ANALYTIC BOUND ORBITALS ARE BEING USED
C
         PV=ZERO
         ALPHA=(R-R0)/SIGMA
         IF(-ALPHA.LT.TWELVE)PV=-(TWO*ALPHA*ALPHA-ONE)*(TWO/(SIGMA*SIGMA
     *   ))*C1*EXP(-ALPHA*ALPHA)
         DO 3 JK=1,M2
            IR = IRAD(L2,N2,JK)
            ZEE=-ZE(L2,N2,JK)
            ATERM=DBLE(IR*(IR-1))+TWO*ZEE*DBLE(IR)*R+ZEE*ZEE*R*R
            PV=PV+C(L2,N2,JK)*EXP(ZEE*R)*ATERM*R**(IR-2)
3        CONTINUE
      ELSE
C
C        NUMERICAL BOUND ORBITALS ARE BEING USED
C
         PV=UDP(I2,J1)
      ENDIF
C
C      DETERMINE THE SIMPSONS RULE WEIGHTING FACTOR
C
      GO TO (4,5),MI
    4 MI = 2
      WEIGHT = FOUR
      GO TO 8
    5 MI=1
      IF(JP.EQ.IFI) GO TO 6
      WEIGHT = TWO
      GO TO 8
    6 IF(I.EQ.NIX) GO TO 7
      WEIGHT = THREE
      GO TO 8
    7 WEIGHT = ONE
C
C      ADD IN THE CONTRIBUTIONS
C
    8 SUM=SUM+UJ(J1,I1)*(X1+PV)*WEIGHT
    9 CONTINUE
      IST = IFI+1
      ALBVAL = ALBVAL-SUM*H/SIX
   10 CONTINUE
      RETURN
C
C      BOTH ORBITALS ARE CONTINUUM ORBITALS. THE SCHIMDT COEFFICIENTS
C      ARE USED TO EXPRESS ONE CONTINUUM ORBITAL IN TERMS OF ORBITALS
C      SATISFYING A DIFFERENTIAL EQUATION AND BOUND ORBITALS. THE
C      SECOND DIFFERENTIATION CAN THEN BE CARRIED OUT ANALYTICALLY.
C
   11 CONTINUE
      I1=IPOS(N1,L1)
      I2=IPOS(N2,L2)
C
C      ITEST DETERMINES WHETHER SCHMIDT ORTHOGONALIZATION IS USED
C
      MAXLG = MAXHF
      IF(L1.LE.LRANG1) MAXLG = MAXNLG(L1)
      ITEST = MAXHF-MAXLG
C
C      FIRST EVALUATE CONTRIBUTION FROM THE CONTINUUM ORBITALS USING
C      SIMPSONS RULE
C
   12 SUM = ZERO
      X   = ZERO
      IST = 1
      DO 17 I=1,NIX
      H=HINT*DBLE(IHX(I))
      MI=1
      ADD = ZERO
      IFI = IRX(I)
      DO 16 J=IST,IFI
      JP=J+1
      X = X+H
C WEIGHT
      GO TO (14,15),MI
   14 MI = 2
      WEIGHT=FOUR
      GO TO 13
   15 MI=1
      IF(J.EQ.IFI)GO TO 69
      WEIGHT=TWO
      GO TO 13
   69 IF(I.EQ.NIX)GO TO 68
      WEIGHT=THREE
      GO TO 13
   68 WEIGHT=ONE
C ADD
   13 IF(IPSEUD.GT.0) POVAL2=POTHAM(LLL,J)
      IF(IPSEUD.LT.0)POVAL2=TWO*POTHAM(LLL,J)/X
      IF(IPSEUD.EQ.0) POVAL2=ZN/X
      POVAL1=POVALU(2*J)
      ADD=ADD+UJ(JP,I1)*UJ(JP,I2)*WEIGHT*(POVAL2-POVAL1)
   16 CONTINUE
      IST = IFI+1
      ADD = ADD*H/THREE
      SUM=SUM+ADD
   17 CONTINUE
      ALBVAL = -SUM/TWO
C
C      NOW ADD IN ENERGY TERM
C
      N1P=N1-MAXHF+ITEST
      N2P=N2-MAXHF+ITEST
C      N3P=N1-MAXHF
      N4P=N2-MAXHF
      IF(ITEST.EQ.0) GO TO 21
      DO 18 I=1,N4P
      I3 = I+ITEST
      ALBVAL=ALBVAL+B(N1P,I3)*B(N2P,I3)*EIGENS(I,L1)/TWO
   18 CONTINUE
      DO 20 I=1,ITEST
      DO 19 J=1,N4P
      J3=J+ITEST
      ALBVAL=ALBVAL+B(N1P,I)*B(N2P,J3)*OVRLAP(J,I)*EIGENS(J,L1)/TWO
   19 CONTINUE
   20 CONTINUE
      GO TO 22
   21 IF(N1.NE.N2) GO TO 22
      ALBVAL = ALBVAL+EIGENS(N1P,L1)/TWO
C
C  NRB:
C      ADD-IN THE MASS-VELOCITY CONTRIBUTION (CONTINUUM-CONTINUUM ONLY)
C
   22 IF(IRELMV.EQ.0)GO TO 79
C TEST ZERO-OUT UNPHYSICALLY LARGE DIAGONAL ELEMENT (AS PER STG1 GEN1CC)
      IF(N1.EQ.N2)GO TO 79
      SUM = ZERO
      IST = 1
      POVAL2 = EIGENS(N1P,L1) + EIGENS(N2P,L2)
      DO 77 I=1,NIX
      H=HINT*DBLE(IHX(I))
      MI=1
      ADD = ZERO
      IFI = IRX(I)
      DO 76 J=IST,IFI
      JP=J+1
C WEIGHT
      GO TO (74,75),MI
   74 MI = 2
      WEIGHT=FOUR
      GO TO 78
   75 MI=1
      IF(J.EQ.IFI)GO TO 72
      WEIGHT=TWO
      GO TO 78
   72 IF(I.EQ.NIX)GO TO 73
      WEIGHT=THREE
      GO TO 78
   73 WEIGHT=ONE
C ADD
   78 POVAL1=POVALU(2*J)
      ADD=ADD+UJ(JP,I1)*UJ(JP,I2)*WEIGHT*(POVAL2+POVAL1)*POVAL1
   76 CONTINUE
      IST = IFI+1
      ADD = ADD*H/THREE
      SUM=SUM+ADD
   77 CONTINUE
C
C      NOW ADD IN ENERGY TERM
C
      IF(N1.EQ.N2)SUM = SUM + EIGENS(N1P,L1) * EIGENS(N2P,L1)
C
      SUM = -SUM*FSC*FSC/EIGHT
      ALBVAL = ALBVAL + SUM
C
C     END OF CONTINUUM-CONTINUUM MASS-VELOCITY CONTRIBUTION
C
   79 IF(ITEST.EQ.0) RETURN
C
C      NOW ADD IN CONTRIBUTION FROM THE BOUND ORBITALS USING
C      SIMPSONS RULE
C
      DO 31 I=1,ITEST
      N2=MAXLG+I
      I2=IPOS(N2,L2)
      M=0
      DO 80,K1=1,L2-1
         DO 90,K2=K1,MAXNHF(K1)
            M=M+1
90       CONTINUE
80    CONTINUE
      M2=NCO(M+N2-L2+1)
      CALL CORECT(N2,L2,M2,R0,SIGMA,C1)
      IST = 1
      R = ZERO
C
      DO 30 J=1,NIX
      H=HINT*DBLE(IHX(J))
      MI = 1
      SUM = ZERO
      IFI=IRX(J)
      DO 29 K=IST,IFI
      R=R+H
C
C      DETERMINE THE CONTRIBUTION INVOLVING THE ANGULAR MOMENTUM AND
C      THE POTENTIAL
C
      RINV=ONE/R
      X1=UJ(K+1,I2)*(POVALU(2*K)-ALSQ*RINV*RINV)
C
C      DETERMINE THE CONTRIBUTION INVOLVING THE SECOND DERIVATIVE
C      ACTING ON THE BOUND ORBITAL
C
      PV=ZERO
      ALPHA=(R-R0)/SIGMA
      IF(-ALPHA.LT.TWELVE)PV=-(TWO*ALPHA*ALPHA-ONE)*(TWO/(SIGMA*SIGMA))
     *   *C1*EXP(-ALPHA*ALPHA)
      DO 23 JK=1,M2
         IR=IRAD(L2,N2,JK)
         ZEE=-ZE(L2,N2,JK)
         ATERM=DBLE(IR*(IR-1))+TWO*ZEE*DBLE(IR)*R+ZEE*ZEE*R*R
         PV=PV+C(L2,N2,JK)*EXP(ZEE*R)*ATERM*R**(IR-2)
   23 CONTINUE
C
      GO TO (24,25),MI
   24 MI=2
      WEIGHT = FOUR
      GO TO 28
   25 MI=1
      IF(K.EQ.IFI) GO TO 26
      WEIGHT = TWO
      GO TO 28
   26 IF(J.EQ.NIX) GO TO 27
      WEIGHT = THREE
      GO TO 28
   27 WEIGHT = ONE
C
   28 SUM = SUM+UJ(K+1,I1)*(X1+PV)*WEIGHT
   29 CONTINUE
      IST = IFI+1
      ALBVAL=ALBVAL-SUM*H*B(N2P,I)/SIX
   30 CONTINUE
   31 CONTINUE
C
      RETURN
      END
C
C*********************************************************************
C
      SUBROUTINE PMUL(A,B,C,K,N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C
C     MULTIPLICATION OF VECTOR BY MATRIX OR MATRIX BY MATRIX
C     K=1 FOR B,C VECTORS, K=N FOR B,C MATRICES
C     OUTPUTS RESULTS FOR C(N,N) OR C(N,1)
C
C     DIMENSION A(N,N),B(N,K),C(N,K),AB(N)
      DIMENSION A(3,3),B(3,K),C(3,K),AB(3)
      DO 20 J=1,K
      DO 20 I=1,N
      DO 18 KL=1,N
   18 AB(KL)=A(I,KL)*B(KL,J)
      SUM=0.0D0
      DO 19 KL=1,N
   19 SUM=SUM+AB(KL)
   20 C(I,J)=SUM
      RETURN
      END
C=====================================================================

      SUBROUTINE POTF
C
C      THIS ROUTINE AUTOMATICALLY CALCULATES THE STATIC POTENTIAL OF THE
C      LOWEST POSSIBLE TARGET CONFIGURATION WITH THE GIVEN RADIAL
C      ORBITALS.
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)

      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM10=IDIM7*2,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3),
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(P75=0.75D0,P125=0.125D0,SIXT=16.0D0,TNINE=9.0D0)

      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /POTVAL/ POVALU(KDIM10)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS

 1001 FORMAT (///'  POVALU'/)
 1002 FORMAT ('  ', 4E16.7)
C
      DO 60,I=1,IRX(NIX)*2
         POVALU(I)=ZERO
60    CONTINUE

      DO 10,I=1,LRANG1
         DO 20,J=I,MAXNHF(I)
            NQ=IPOS(J,I)
            NUMORB=ISTAT(I,J)
            IF (NUMORB.EQ.0) THEN
               GOTO 20
            ELSE
               K=IRX(NIX)+1
               SUM1=ZERO
               SUM2=ZERO
               POVALU(IRX(NIX)*2)=POVALU(IRX(NIX)*2)+ONE*NUMORB
               DO 30,L=2*IRX(NIX)-2,2,-2
                  K=K-1
                  SUM1=SUM1+(STEP(K+1)/TWO)*(UJ(K,NQ)**2+UJ(K+1,NQ)**2)
                  SUM2=SUM2+(STEP(K+1)/TWO)*(UJ(K,NQ)**2/XR(K)+
     *                                     UJ(K+1,NQ)**2/XR(K+1))
                  POVALU(L)=POVALU(L)+(ONE-SUM1+XR(K)*SUM2)*NUMORB
30             CONTINUE
            ENDIF
20       CONTINUE
10    CONTINUE

      DO 40,I=2,2*IRX(NIX),2
         POVALU(I)=TWO*(DBLE(NZ)-POVALU(I))
40    CONTINUE


C     NOW INTERPOLATE THE MISSING POINTS USING A 4-POINT INTERPOLATION
C     FORMULA.... A 3-POINT FORMULA IS USED AT THE ORIGIN

      POVALU(1)=P75*(DBLE(NZ)+POVALU(2))-P125*POVALU(4)
      IX=5
      A=TNINE/SIXT
      B=ONE /SIXT
      POVALU(3)=A*(POVALU(4)+POVALU(2))-B*(POVALU(6)+TWO*DBLE(NZ))
      DO 50,I=1,NIX
         DO 70,J=IX,IRX(I)*2-3,2
            POVALU(J)=A*(POVALU(J+1)+POVALU(J-1))-B*(
     *                   POVALU(J+3)+POVALU(J-3))
70       CONTINUE
         IF (I.NE.NIX) THEN
C           CHANGE OF INTERVAL...DOUBLING OF STEP SIZE IS ASSUMED
            K1=2*IRX(I)+1
            K2=K1-2
            POVALU(K1)=A*(POVALU(K1+1)+POVALU(K1-1))-B*(
     *                    POVALU(K1+3)+POVALU(K1-5))
            POVALU(K2)=A*(POVALU(K2+1)+POVALU(K2-1))-B*(
     *                    POVALU(K2+2)+POVALU(K2-3))
            IX=K1+2
            GOTO 50
         ELSE
C           FINAL POINT TO BE INTERPOLATED
            K1=IRX(NIX)*2-1
            POVALU(K1)=A*(POVALU(K1+1)+POVALU(K1-1))-B*(
     *                    POVALU(K1-3)+DBLE(2*(NZ-NELC)))
         ENDIF
50    CONTINUE
      K=1
      DO 80,I=2,2*IRX(NIX),2
         K=K+1
         POVALU(I)=POVALU(I)/XR(K)
         POVALU(I-1)=POVALU(I-1)/(XR(K)-STEP(K)/TWO)
80    CONTINUE
C
      IF (NBUG5 .EQ. 1) THEN
         WRITE(IWRITE,1001)
         WRITE (IWRITE,1002) (POVALU(I), I=1,IRX(NIX)*2)
      END IF
C
      RETURN
      END
C=======================================================================
      SUBROUTINE RADINT(N1,L1,N2,L2,K,X)
C
C      EVALUATES X, THE RADIAL MULTIPOLE INTEGRAL OF ORDER K BETWEEN
C      TWO NUMERICAL ORBITALS SPECIFIED BY THE QUANTUM NUMBERS (N1,L1-1)
C      AND (N2,L2-1).
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,
     *  IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM1=IDIM13+2*IDIM3)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3),
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)

      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS

      X=ZERO
      J1=IPOS(N1,L1)
      J2=IPOS(N2,L2)
      DO 10,I=1,IRX(NIX)+1
         X=X+WT(I)*UJ(I,J1)*UJ(I,J2)*XR(I)**K
10    CONTINUE
      RETURN
      END
C=====================================================================
      SUBROUTINE RECOV1(I,IDAT,ICURR)
C
C      THIS ROUTINE IS CALLED ONLY IN THE CASE OF ARRAY OVERFLOW.
C
C      I IS THE POSITION IN THE IDRAY AND IPRAY LISTS HOLDING
C      THE NAME OF THE RELEVANT DIMENSION PARAMETER AND
C      PREPROCESSOR PARAMETER.
C      IDAT IS THE ARRAY SIZE REQUIRED BY THE DATA
C      ICURR IS THE CURRENT ARRAY SIZE
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      CHARACTER*6 IDRAY
      CHARACTER*3 IPRAY
C      PARAMETER(IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO,
C     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM,
C     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB,IDIM14=MZFAC)
C      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4,ISDIM1=MZSET,
C     1  ISDIM2=MZNCS,ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,
C     2  ITDIM6=MZCHS,IPDIM1=MZSPN,ICDIM1=MZNCF,ICDIM2=MZORB)
C      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      DIMENSION IDRAY(27),IPRAY(27)
      DATA IDRAY/'IDIM1','IDIM2','IDIM3','IDIM4','IDIM5','IDIM6',
     1        'IDIM7','IDIM8','IDIM9','IDIM10','IDIM11','IDIM12',
     2        'IDIM13','IDIM14','ILDIM1','ILDIM3','ISDIM1',
     3        'ISDIM2','ITDIM1','ITDIM2','ITDIM3','ITDIM6',
     4        'IPDIM1','ICDIM1','ICDIM2','IRDIM1','LBUFF1'/
      DATA IPRAY/'LR1','LR2','NR2','LMX','NPO','NIX',
     1           'PTS','AMP','LAG','SHM','SLT','NR1',
     2           'ORB','FAC','LR3','LR4','SET',
     3           'NCS','TAR','NSS','CHF','CHS',
     4           'SPN','NCF','ORB','RKC','BUF'/
 1000 FORMAT(/' * ARRAY OVERFLOW *'/
     1    /1X,A6,' (MZ',A3,') SHOULD BE INCREASED FROM',
     2    I7,' TO AT LEAST ',I7)
C
 1001 FORMAT(/' PROGRAM TERMINATES IN RECOV1'/)
 1002 FORMAT(/' CHECK TO SEE IF OTHER ARRAYS ARE GOING TO BE EXCEEDED')
C
      WRITE(IWRITE,1000)IDRAY(I),IPRAY(I),ICURR,IDAT
    1 WRITE(IWRITE,1001)
      CALL EXIT (0)
C
      RETURN
      END
C=======================================================================
      SUBROUTINE RESETL(MJS0,MINIM,ICONTN)
C
C      READS STGANG DATA. CALCULATES LRANG2, THE RANGE OF
C      VALUES OF CONTINUUM ANGULAR MOMENTUM +1
C      ARRANGES CONFIGURATIONS INTO SETS HAVING THE SAME L,S,PI
C      THE CONFIGURATIONS ARE ASSUMED TO BE READ IN ALREADY
C      GROUPED ACCORDING TO L,S,PI
C      FOLLOWS THE CODING OF CUSETL IN STGANG
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      CHARACTER*4 TITLE,CONT
      CHARACTER*3 RELOP
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ICDIM3=ICDIM2+ICDIM2-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *  	      NBUG9
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
C
      COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN
      COMMON /NRB/MAXC,LFIXN0
      COMMON /NRBREL/IRELMV
C
      DATA CONT/'CONT'/
      DIMENSION LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1              NCSET(ISDIM1)
      DIMENSION NJCOMP(ICDIM2),LJCOMP(ICDIM2)
      DIMENSION NOCCSH(ICDIM1),NOCORB(ICDIM2,ICDIM1),
     1          NELCSH(ICDIM2,ICDIM1),J1QNRD(ICDIM3,3,ICDIM1)
      DIMENSION LVAL(ILDIM2),ILVAL(ILDIM2),NSCOL(ILDIM2)
      DIMENSION TITLE(18)
C
      NAMELIST/STGNX/MINLT,MAXLT,LRGLLO,LRGLUP,MAXC,LFIXN,MJS,NPTYIN
     X ,RELOP,NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9
C
 1000 FORMAT(//'  ',10X,'STGANG DATA')
 1001 FORMAT('    LSET  ISPSET  IPISET   NCSET')
 1002 FORMAT(10(2X,I6))
 1005 FORMAT ('  HAMILTONIAN MATRICES ARE TO BE CALCULATED OVER VALUES'
     1       /'  OF THE TOTAL ANGULAR MOMENTUM FROM ',I2,' TO ',I2/
     2        '  AND TOTAL PARITY EVEN AND ODD IN EACH CASE'/)
C
 1007 FORMAT(///48X,'TARGET INPUT DATA')
 1008 FORMAT(///' NUMBER OF CONFIGURATIONS =',I3)
 1009 FORMAT(' NUMBER OF OCCUPIED SHELLS IN THESE CONFIGURATIONS'/30I3)
 1010 FORMAT(' CONFIGURATION',I3/6X,' OCCUPIED ORBITALS ARE',32X,10I3)
 1011 FORMAT(5X,' NUMBER OF ELECTRONS IN RESPECTIVE OCCUPIED SHELLS',
     15X,10I3)
 1012 FORMAT(5X,' COUPLING SCHEME')
 1013 FORMAT(5X,3I3,6(I10,2I3))
 1014 FORMAT(/' THERE ARE',I3,' ORBITALS WHOSE N,L VALUES ARE'/
     1       6X,10(I8,I3))
 1019 FORMAT(//' ****** LRANG2 ******'/
     1         '  L+1 TAKES VALUES BETWEEN ',I3,' AND ',I3)
C
C      THE FOLLOWING FORMAT STATEMENTS ARE TO READ THE CARD INPUT DATA.
C
 2000 FORMAT(12I5)
 2001 FORMAT(18A4)
C
C NRB
      NPTYMN=0
      NPTYMX=1
      NPTYIN=-1
      IRELMV=0
C
C
      IREAD=5
      IWRITE=6
      OPEN(UNIT=IREAD,FILE='dstgnx',STATUS='UNKNOWN')
      OPEN(UNIT=IWRITE,FILE='rout2nx',STATUS='UNKNOWN')
      REWIND (IREAD)
      READ(IREAD,2001) (TITLE(I),I=1,18)
      IF (TITLE(1) .EQ. CONT) THEN
         INX2=36
C        OPEN (INX2,FILE='tapenx2',ACCESS='SEQUENTIAL',
         OPEN (INX2,FILE='NX2.DAT',ACCESS='SEQUENTIAL',
     1         STATUS='OLD',FORM='UNFORMATTED')
         REWIND (INX2)
         ICONTN=1
         READ(INX2,ERR=301)MAXN,MAXORB,LRANG3,NSETS1,MAXNCF,NAST,MAXNST,
     1             NCHSUM,MAXNCH,ISRAN3,NCFG,IRELMV
 301     IF (ICDIM1 .LT. NCFG) CALL RECOV1(24,NCFG,ICDIM1)
      ELSE
         ICONTN=0
         READ(IREAD, *)  NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,
     1                   NBUG8,NBUG9
      END IF
C
CW    WRITE(IWRITE,1000)
C
      MAXC=1
      LFIXN0=-1
      IF(ICONTN.EQ.1)THEN
      MINLT=0
      MAXLT=0
      LRGLLO=0
      LRGLUP=0
      MJS=-1
      LFIXN=-1
      RELOP='  '
      NBUG1=0
      NBUG2=0
      NBUG3=0
      NBUG4=0
      NBUG5=0
      NBUG7=0
      NBUG8=0
      NBUG9=0
C
      READ(IREAD,STGNX)
C
      IF(RELOP.EQ.'MVD'.OR.RELOP.EQ.'YES')IRELMV=1
      IF(RELOP.EQ.'NO')IRELMV=0
      IF(NPTYIN.EQ.0)NPTYMX=0
      IF(NPTYIN.EQ.1)NPTYMN=1
      LFIXN0=LFIXN
      IF(MJS.LT.0)MJS=1
      MJS0=MJS
      IF(MINLT.GT.0)LRGLLO=MINLT
      IF(MAXLT.GT.0)LRGLUP=MAXLT
      ELSE
      READ (IREAD,*) LRGLLO,LRGLUP
      ENDIF
C
C
CW    WRITE (IWRITE,1005) LRGLLO,LRGLUP
      IF (ILDIM3 .LE. LRGLUP) CALL RECOV1(16,LRGLUP+1,ILDIM3)
C
CW    WRITE(IWRITE,1007)
C
C      READ IN THE QUANTUM NUMBERS DEFINING THE CONFIGURATIONS
C
      IF (ICONTN .EQ. 1) THEN
         READ(INX2) (NOCCSH(I),I=1,NCFG)
CW       WRITE(IWRITE,1008)NCFG
         DO 101 I1=1,NCFG
            N1=NOCCSH(I1)
            READ(INX2) (NOCORB(J,I1),J=1,N1),(NELCSH(J,I1),J=1,N1)
            M1=N1+N1-1
            READ(INX2)((J1QNRD(J,K,I1),K=1,3),J=1,M1)
  101    CONTINUE
         READ(INX2)(NJCOMP(I),LJCOMP(I),I=1,MAXORB)
         REWIND (INX2)
      ELSE
         READ(IREAD,*)  MAXORB,NELC,NSET,NKEY
         IF (ICDIM2 .LT. MAXORB) CALL RECOV1(25,MAXORB,ICDIM2)
         READ(IREAD,*)(NJCOMP(I),LJCOMP(I),I=1,MAXORB)
         IREAD1=IREAD
         IF (NKEY .NE. 2) THEN
            IF (NKEY .EQ. -2) THEN
               READ(IREAD,*) NOCC
               IF (NOCC .GT. ICDIM1) CALL RECOV1(24,NOCC,ICDIM1)
               READ(IREAD,*) (NOCCSH(I),I=1,NOCC)
               DO 91 I=1,NOCC
                  JACT=NOCCSH(I)
                  READ(IREAD,*)
                  READ(IREAD,*)
   91          CONTINUE
            ELSE
               READ(IREAD,*)  NOPTN
               READ(IREAD,*)
               READ(IREAD,*)
               DO 15 M=1,NOPTN
                  READ(IREAD,*)
   15          CONTINUE
            END IF
            DO 16 N=1,NSET
               READ(IREAD,*)  LL,LSPN,LPTY
   16       CONTINUE
            IREAD1=36
            OPEN(UNIT=36,FILE='CONF.DAT',ACCESS='SEQUENTIAL',
     1           STATUS='OLD',FORM='UNFORMATTED')
            REWIND (36)
         END IF
C
C      READ IN THE QUANTUM NUMBERS DEFINING THE CONFIGURATIONS
C
         READ(IREAD1,*)NCFG
CW       WRITE(IWRITE,1008)NCFG
         IF (ICDIM1 .LT. NCFG) CALL RECOV1(24,NCFG,ICDIM1)
         READ(IREAD1,*) (NOCCSH(I),I=1,NCFG)
CW       WRITE(IWRITE,1009) (NOCCSH(I),I=1,NCFG)
         DO 192 I=1,NCFG
            N=NOCCSH(I)
            READ(IREAD1,*) (NOCORB(J,I),J=1,N)
CW          WRITE(IWRITE,1010) I,(NOCORB(J,I),J=1,N)
            READ(IREAD1,*) (NELCSH(J,I),J=1,N)
CW          WRITE(IWRITE,1011) (NELCSH(J,I),J=1,N)
            M=N+N-1
            READ(IREAD1,*) ((J1QNRD(J,K,I),K=1,3),J=1,M)
CW          WRITE(IWRITE,1012)
CW          WRITE(IWRITE,1013)((J1QNRD(J,K,I),K=1,3),J=1,M)
  192    CONTINUE
         IF (NKEY .NE. 2) REWIND(36)
C
      END IF
C
CW    WRITE(IWRITE,1001)
C
C      IS WILL COUNT THE SETS
C      LRANG3 WILL CONTAIN MAX. CONFIGURATION ANGULAR MOMENTUM + 1
C
      IS=0
      LRANG3=0
C
C      LOOP OVER ALL CONFIGURATIONS READ
C
      DO 2     IC=1,NCFG
              ISH=NOCCSH(IC)
             ISH2=ISH+ISH-1
           IPARTY=0
C
C      LOOP OVER THE SHELLS OF THIS CONFIG. TO FIND THE PARITY
C
         DO 3    ISHEL=1,ISH
                    IL=NOCORB(ISHEL,IC)
                IPARTY=IPARTY+LJCOMP(IL)*NELCSH(ISHEL,IC)
    3    CONTINUE
C
           IPARTY=MOD(IPARTY,2)
             LCFG=(J1QNRD(ISH2,2,IC)-1)/2
            ISPIN= J1QNRD(ISH2,3,IC)
C
       IF(IS .NE. 0) THEN
         IF (LCFG .EQ. LSET(IS)  .AND.
     1   ISPIN .EQ. ISPSET(IS)  .AND.
     2  IPARTY .EQ. IPISET(IS)) THEN
C
C      FILL CURRENT SET
C
        ICS=ICS+1
         IF (ISDIM2 .LT. ICS) CALL RECOV1(18,ICS,ISDIM2)
       NCSET(IS)=ICS
         GO TO 2
         END IF
      END IF
C
C      START NEW SET
C
         IS=IS+1
         IF (ISDIM1 .LT. IS) CALL RECOV1(17,IS,ISDIM1)
        ICS=1
        LSET(IS)=LCFG
      LRANG3=MAX(LRANG3,LCFG)
      ISPSET(IS)=ISPIN
      IPISET(IS)=IPARTY
       NCSET(IS)=1
C
    2 CONTINUE
C
      NSET=IS
C
      DO 4   I=1,NSET
CW      WRITE(IWRITE,1002)LSET(I),ISPSET(I),IPISET(I),NCSET(I)
    4 CONTINUE
      LRANG3=LRANG3+1
      IF (ILDIM1 .LT. LRANG3) CALL RECOV1(15,LRANG3,ILDIM1)
      LLPDIM=LRANG3+LRANG3-1
C
C      TO DETERMINE WHICH SETS COUPLE WITH EACH VALUE OF THE
C      CONTINUUM ANGULAR MOMENTUM WITHIN THE RANGE ALLOWED
C      THUS DETERMINE THE RANGE REQUIRED FOR THE CONTINUUM
C      ANGULAR MOMENTUM
C
      MINIM=999
      LRANG2=1
      DO 5   LRGL=LRGLLO,LRGLUP
      DO 6   NPTY0=NPTYMN,NPTYMX
      IF(NPTYIN.LT.0)THEN
      NPTY=NPTY0
      ELSE
      NPTY=MOD(LRGL,2)+NPTY0
      NPTY=MOD(NPTY,2)
      ENDIF
C
C      FIND RANGE OF CONTINUUM ANGULAR MOMENTUM AND
C      SET UP ILVAL(KL) - THE POSSIBLE VALUES OF CONTINUUM
C      ANGULAR MOMENTUM
C
      LMIN=999
      LMAX=0
C
      DO 7    IS=1,NSET
            LCFG=LSET(IS)
            LMIN=MIN(LMIN,ABS(LCFG-LRGL))
            LMAX=MAX(LMAX, LCFG)
    7 CONTINUE
      LMAX=LMAX+LRGL
      ILTOT=LMAX-LMIN+1
C
      DO 8    KL=1,ILTOT
        ILVAL(KL)=LMIN+KL-1
    8 CONTINUE
C
C      LOOP OVER THIS RANGE OF CONTINUUM L
C
      DO 9    KL=1,ILTOT
               L=ILVAL(KL)
              NL=0
C
        DO 10    IS=1,NSET
              LCFG=LSET(IS)
               IPI=IPISET(IS)
C
C      PARITY TEST
C
         IF(MOD(IPI+L+NPTY,2) .NE. 0) GO TO 10
C
C      TRIANGLE RULE
C
         IF(LRGL .GT. LCFG+L  .OR.
     1      LRGL .LT. ABS(LCFG-L)) GO TO 10
C
                NL=NL+1
         NSCOL(KL)=1
         GO TO 9
   10   CONTINUE
C
    9 CONTINUE
C
C      NOW DELETE VALUES OF THE CONTINUUM ANGULAR MOMENTUM WHICH
C      CANNOT COUPLE WITH ANY SET AND CONSTRUCT /SETL/
C
      KL=0
      DO 11   IKL=1,ILTOT
         IF (NSCOL(IKL) .NE. 0) THEN
            KL=KL+1
            LVAL(KL)=ILVAL(IKL)
         END IF
   11 CONTINUE
      LTOT=KL
C
      IF (LTOT .EQ. 0) GO TO 6
C
      LMIN=LVAL(1)
      LMAX=LVAL(LTOT)
         LRANG2=MAX(LMAX+1, LRANG2)
         MINIM=MIN (LMIN+1, MINIM)
    6 CONTINUE
    5 CONTINUE
CW    WRITE (IWRITE,1019) MINIM, LRANG2
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE ROOT(T,FT,B,C,RELERR,ABSERR,IFLAG)
C
C    THIS SUBROUTINE IS TAKEN FROM THE BOOK OF: SHAMPINE AND GORDON,
C    'COMPUTATIONAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS'
C
C  ROOT COMPUTES A ROOT OF THE NONLINEAR EQUATION F(X)=0, WHERE F(X)
C  IS A CONTINUOUS REAL FUNCTION OF A SINGLE REAL VARIABLE X. THE
C  METHOF OF SOLUTION IS A COMBINATION OF BISECTION AND THE SECANT RULE.
C
C  NORMAL INPUT CONSISTS OF A CONTINUOUS FUNCTION F AND IN INTERVAL
C  (B,C) SUCH THAT F(B)*F(C).LE.0.0. EACH ITERATION FINDS NEW VALUES OF
C  B AND C SUCH THAT THE INTERVALL (B,C) IS SHRUNK AND F(B)*F(C).LE.0.0.
C
C  THE STOPPING CRITERION IS
C
C           ABS(B-C).LE.2.0*(RELERR*ABS(B)+ABSERR)
C
C  WHERE RELERR=RELATIVE ERROR AND ABSERR=ABSOLUTE ERROR ARE INPUT
C  QUANTITIES. SET THE FLAG, IFLAG, POSITIVE TO INITIALISE THE
C  COMPUTATION. AS B,C AND IFLAG ARE USED FOR BOTH INPUT AND OUTPUT,
C  THEY MUST BE VARIABLES IN THE CALLING PROGRAM.
C  IF 0.0 IS A POSSIBLE ROOT, ONE SHOULD NOT CHOOSE ABSERR=0.0.
C
C
C  THE OUTPUT VALUE OF B IS THE BETTER APPROXIMATION TO A ROOT AS
C  B AND C ARE ALWAYS REDEFINED SO THAT ABS(F(B)).LE.ABS(F(C)).
C
C  TO SOLVE THE EQUATION, ROOT MUST EVALUATE F(X) REPEATEDLY. THIS IS
C  DONE IN THE CALLING PROGRAM. WHEN AN EVALUATION OF F IS NEEDED AT T,
C  ROOT RETURNS WITH IFLAG NEGATIVE. EVALUATE FT=F(T) AND CALL ROOT
C  AGAIN. DO NOT ALTER IFLAG.
C
C  WHEN THE COMPUTATION IS COMPLETE, ROOT RETURNS TO THE CALLING
C  PROGRAM WITH IFLAG POSITIVE:
C
C     IFLAG=1  IF F(B)*F(C).LT.0 AND THE STOPPING CRITERION IS MET.
C
C          =2  IF A VALUE B IS FOUND SUCH THAT THE COMPUTED VALUE F(B)
C              IS EXACTLY ZERO. THE INTERVAL (B,C) MAY NOT SATISFY
C              THE STOPPING CRITERION.
C
C          =3  IF ABS(F(B)) EXCEEDS THE INPUT VALUES ABS(F(B)),
C              ABS(F(C)).  IN THIS CASE IT IS LIKELY THAT B IS CLOSE
C              TO A POLE OF F.
C
C          =4  IF NO ODD ORDER ROOT WAS FOUND IN THE INTERVALL.
C              A LOCAL MININMUM MAY HAVE BEEN OBTAINED.
C
C          =5  IF TOO MANY FUNCTION EVALUATIONS WERE MADE.
C              (AS PROGRAMMED, 500 ARE ALLOWED.)
C
C  THIS CODE IS A MODIFICATION OF THE CODE  Z E R O I N  WHICH IS
C  COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, NUMERICAL
C  COMPUTING: AN INTRODUCTION  BY L.F. SHAMPINE AND R.C. ALLEN
C
C***********************************************************************
C  THE ONLY MACHINE DEPENDENT CONSTANT IS BASED ON THE MACHINE UNIT    *
C  ROUNDOFF ERROR U. IT IS CALCULATED IN THE SUBROUTINE  M A C H I N   *
C  WHICH MUST HAVE BEEN CALLED BEFORE THE FIRST CALL OF  R O O T.      *
C  NOTE....IN THIS CODE MACHIN IS NOT CALLED BUT U IS SET TO 1.0E-12   *
C          IN THE PARAMETER STATEMENT BELOW.                           *
C***********************************************************************
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER (ZERO=0.0D0,ONE=1.0D0,EIGHT=8.0D0,HALF=0.5D0)
      PARAMETER (U=1.0D-12)
C
      SAVE RE,AE,IC,ACBS,A,FA,FC,FX,KOUNT                    !NRB
C
      IF (IFLAG.GE.0) GOTO 100
      IFLAG=ABS(IFLAG)
      GOTO (200,300,400), IFLAG
  100 RE=MAX(RELERR,U)
      AE=MAX(ABSERR,ZERO)
      IC=0
      ACBS=ABS(B-C)
      A=C
      T=A
      IFLAG=-1
      RETURN
  200 FA=FT
      T=B
      IFLAG=-2
      RETURN
  300 FB=FT
      FC=FA
      KOUNT=2
      FX=MAX(ABS(FB),ABS(FC))
    1 IF(ABS(FC).GE.ABS(FB)) GOTO 2
C
C  INTERCHANGE B AND C SO THAT ABS(F(B)).LE.ABS(F(C))
C
      A=B
      FA=FB
      B=C
      FB=FC
      C=A
      FC=FA
    2 CMB=HALF*(C-B)
      ACMB=ABS(CMB)
      TOL=RE*ABS(B)+AE
C
C  TEST STOPPING CRITERION AND FUNCTION COUNT
C
      IF(ACMB.LE.TOL) GOTO 8
      IF(KOUNT.GE.500) GOTO 12
C
C  CALCULATE NEW ITERATE IMPLICITLY AS B+P/Q  WHERE WE ARRANGE P.GE.0.
C  THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW.
C
      P=(B-A)*FB
      Q=FA-FB
      IF(P.GE.ZERO) GOTO 3
      P=-P
      Q=-Q
C
C  UPDATE A, CHECK REDUCTION OF THE SIZE OF BRACKETING INTERVAL IS
C  SATISFACTORY. IF NOT, BISECT UNTIL IT IS.
C
    3 A=B
      FA=FB
      IC=IC+1
      IF(IC.LT.4) GOTO 4
      IF(EIGHT*ACMB.GE.ACBS) GOTO 6
      IC=0
      ACBS=ACMB
C
C  TEST FOR TOO SMALL A CHANGE
C
    4 IF(P.GT.ABS(Q)*TOL) GOTO 5
C
C  INCREMENT BY TOLERANCE
C
      B=B+SIGN(TOL,CMB)
      GOTO 7
C
C  ROOT OUGHT TO BE BETWEEN B AND (C+B)/2
C
    5 IF(P.GE.CMB*Q) GOTO 6
C
C  USE SECANT RULE
C
      B=B+P/Q
      GOTO 7
C
C  USE BISECTION
C
    6 B=HALF*(C+B)
C
C  HAVE COMPLETED COMPUTATION OF NEW ITERATE B
C
    7 T=B
      IFLAG=-3
      RETURN
  400 FB=FT
      IF(FB.EQ.ZERO) GOTO 9
      KOUNT=KOUNT+1
      IF(SIGN(ONE,FB).NE.SIGN(ONE,FC)) GOTO 1
      C=A
      FC=FA
      GOTO 1
C
C  FINISHED.  SET IFLAG
C
    8 IF(SIGN(ONE,FB).EQ.SIGN(ONE,FC)) GOTO 11
      IF(ABS(FB).GT.FX) GOTO 10
      IFLAG=1
      RETURN
    9 IFLAG=2
      RETURN
   10 IFLAG=3
      RETURN
   11 IFLAG=4
      RETURN
   12 IFLAG=5
      RETURN
      END
C=======================================================================
      SUBROUTINE RS(N2,L2,N1,L1,N4,L4,N3,L3,K,RKVAL)
C
C      THIS SUBROUTINE EVALUATES THE FUNCTION Y(P2,P4,K/R) AND THE
C      SLATER INTEGRAL R(P1,P2,P3,P4,K), WHERE P1,P2,P3,P4 ARE NUMERICAL
C      ORBITAL FUNCTIONS DEFINED OVER A PRESCRIBED SET OF INTEGRATION
C      MESH POINTS.
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM3=MZNR2,
     *  IDIM7=MZPTS,
     *  IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM1=IDIM13+2*IDIM3)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)

      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS
      COMMON /YKSTOR/ YK(IDIM7),RK(IDIM7),RK1(IDIM7),JN2,JL2,JN4,JL4,JKM
      DIMENSION SUM1(IDIM7),SUM2(IDIM7),RINT(IDIM7),AR(IDIM7)
C
C      ********              NOTES FOR THE USER                 ********
C
C      THE FUNCTIONS P1,P2,P3,P4 ARE TO BE STORED IN THE ARRAY UJ(J,I)
C      WITH I=1,4 RESPECTIVELY.....J IS THE INTEGRATION POINT NUMBER
C      WITH J=1 CORRESPONDING TO R=0.0
C
C      ********                                                 ********
C
C      INITIALIZATION OF CONSTANTS AND ARRAY INDICES
C
C      LOCATE THE POSITIONS OF THE P1,P2,P3,P4 FUNCTIONS IN THE UJ ARRAY
C      FROM THEIR N,L VALUES SPECIFIED IN THE ARGUMENT LIST
C
      TWELTH=ONE/TWELVE
      J1=IPOS(N1,L1)
      J2=IPOS(N2,L2)
      J3=IPOS(N3,L3)
      J4=IPOS(N4,L4)
C      K1=K+1
C
C      DO NOT RECALCULATE THE YK INTEGRAL IF IT ALREADY EXISTS
C
      IF(K-JKM) 2,1,2
    1 IF(N2.NE.JN2.OR.L2.NE.JL2) GO TO 2
      IF(N4.NE.JN4.OR.L4.NE.JL4) GO TO 2
      GO TO 7
    2 JN2=N2
      JL2=L2
      JN4=N4
      JL4=L4
      IF(K.EQ.JKM) GO TO 22
      DO 20 J=2,NPTS
   20 RK(J)=XR(J)**K
      DO 21 J=2,NPTS
   21 RK1(J)=ONE/(RK(J)*XR(J))
      JKM=K
C
C      INITIALIZATION
C
   22 ARA=ZERO
      DO 3 J=2,NPTS
      SUM1(J)=UJ(J,J2)*UJ(J,J4)
    3 CONTINUE
      DO 4 J=2,NPTS
      RINT(J)=SUM1(J)*RK(J)
      AR(J)=SUM1(J)*RK1(J)
    4 CONTINUE
C
C       SIMPSONS RULE INTEGRATION FROM 0 TO RA
C
      DO 5 J=2,NPTS
      ARA=ARA+(WT(J)*AR(J))
    5 CONTINUE
C
C       SIMPSONS RULE INTEGRATIONS FROM 0 TO R
C
      YK(1)=ZERO
      RINT(1)=ZERO
      AR(1)=ZERO
      RIN=ZERO
      DO 8 J=3,NPTS,2
      SUM1(J)  =STEP(J)  *THIRD *(RINT(J-2)+FOUR*RINT(J-1)+RINT(J))
      SUM2(J)  =STEP(J)  *THIRD *(AR(J-2)+FOUR*AR(J-1)+AR(J))
    8 CONTINUE
      DO 11 J=3,NPTS,2
      SUM1(J-1)=STEP(J-1)*TWELTH*(FIVE*RINT(J-2)+EIGHT*RINT(J-1)-RINT(J)
     1)
      SUM2(J-1)=STEP(J-1)*TWELTH*(FIVE*AR(J-2)+EIGHT*AR(J-1)-AR(J))
   11 CONTINUE
      DO 9 J=3,NPTS,2
      RINT(J-1)=RIN+SUM1(J-1)
      RIN=RIN+SUM1(J)
      RINT(J)=RIN
      AR(J-1)=ARA-SUM2(J-1)
      ARA=ARA-SUM2(J)
      AR(J)=ARA
    9 CONTINUE
      DO 10 J=2,NPTS
      YK(J)=(RINT(J)*RK1(J)+AR(J)*RK(J))*WT(J)
   10 CONTINUE
C
C      EVALUATE THE SLATER INTEGRAL
C
    7 CONTINUE
      RKVAL=ZERO
      DO 12 J=2,NPTS
         RKVAL=RKVAL+UJ(J,J1)*UJ(J,J3)*YK(J)
12    CONTINUE
      RETURN
      END
C=======================================================================
      SUBROUTINE SCHMDT(LP)
C
C      SCHMIDT ORTHOGONALIZE THE CONTINUUM ORBITALS TO THE BOUND
C      ORBITALS WHICH ARE NOT INCLUDED IN THE ORTHOGONALIZATION IN
C      BASFUN FOR ANGULAR MOMENTUM LP-1.
C      THE SCHMIDT COEFFICIENTS ARE STORED IN THE B-ARRAY
C

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG,IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)

      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)

      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5)
C
C      N1 IS THE NUMBER OF SCHMIDT ORTHOGONALIZED BOUND ORBITALS AND
C      NPTS IS THE NUMBER OF INTEGRATION POINTS
C
      MAXHF=MAXNHF(LP)
      MAXLG=MAXNLG(LP)
      N1=MAXHF-MAXLG
      NPTS=IRX(NIX)+1
C
C      ZEROIZE THE B-ARRAY AND SET THE REQUIRED DIAGONAL ELEMENTS TO
C      UNITY
C
      N=NRANG2+N1
      DO 2 I=1,N
      DO 1 J=1,N
      B(I,J)=ZERO
    1 CONTINUE
    2 CONTINUE
      DO 3 I=1,N1
      B(I,I)=ONE
    3 CONTINUE
C
C      CALCULATE THE SCHMIDT COEFFICIENTS
C
      DO 8 I=1,NRANG2
      N2=N1+I-1
      ANORM=ONE
      DO 5 J=1,N2
      TEMP(J)=ZERO
      DO 4 K=1,N1
      TEMP(J)=TEMP(J)-B(J,K)*OVRLAP(I,K)
    4 CONTINUE
      ANORM=ANORM-TEMP(J)**2
    5 CONTINUE
C
C      RETURN IF NRANG2 IS TOO LARGE
C
      IF (ANORM .LE. TINY) THEN
         LP=-I
         GO TO 15
      END IF
      ANORM=ONE/SQRT(ANORM)
      DO 7 J=1,N2
      B(N2+1,J)=ZERO
      DO 6 K=1,N2
      B(N2+1,J)=B(N2+1,J)+TEMP(K)*B(K,J)
    6 CONTINUE
      B(N2+1,J)=B(N2+1,J)*ANORM
    7 CONTINUE
      B(N2+1,N2+1)=ANORM
    8 CONTINUE
C
C      SCHMIDT ORTHOGONALIZE THE CONTINUUM ORBITALS
C
      I1=IPOS(MAXHF+1,LP)-1
      DO 14 I=1,NPTS
      DO 9 J=1,NRANG2
      I2=I1+J
      TEMP(J)=UJ(I,I2)
    9 CONTINUE
      DO 13 J=1,NRANG2
      I2=I1+J
      SUM=ZERO
      N2=N1+J
      DO 12 K=1,N2
      IF(K.LE.N1)GO TO 10
      I3=K-N1
      X1=TEMP(I3)
      GO TO 11
   10 K1=K+MAXLG
      I3=IPOS(K1,LP)
      X1=UJ(I,I3)
   11 SUM=SUM+X1*B(N2,K)
   12 CONTINUE
      UJ(I,I2)=SUM
   13 CONTINUE
   14 CONTINUE
C
   15 RETURN
      END
C=======================================================================
      SUBROUTINE SHRIEK(NFACT)

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM14=MZFAC)
      COMMON /FACT  / GAMMA(IDIM14)

      GAMMA(1)=DBLE(1)
      IF (NFACT.LT.2) RETURN
      DO 1,I=2,NFACT
         II=I-1
         GAMMA(I)=DBLE(II)*GAMMA(II)
1     CONTINUE
      RETURN
      END
C=======================================================================
      SUBROUTINE SKIPER(L1,L2)
C
C      ROUTINE FOR READING THE CONTINUUM ORBITALS AND THE SCHMIDT
C      COEFFICIENTS CORRESPONDING TO THE ANGULAR MOMENTA (L1-1) AND
C      (L2-1) FROM THE SCRATCH DISC IDISC1 INTO THE UJ ARRAY. IF ONLY
C      ONE ANGULAR MOMENTUM IS REQUIRED THEN L2 MUST BE EQUAL TO ZERO
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C

      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG,IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)

      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)

      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5)

      IF (L1.LE.0) THEN
         PRINT*,' ===== ERROR IN SKIPER ===='
         STOP
      ENDIF

      NPTS=IRX(NIX)+1
      MAXHF=MAXNHF(L1)
      MAXLG=MAXNLG(L1)
      ITEST=MAXHF-MAXLG
      N1=IDIM13+1
      N2=N1+NVARY(L1)-1

      READ(61,REC=L1) ((UJ(I,N),I=1,NPTS),N=N1,N2)
      DO 10,N=1,NVARY(L1)
         IPOS(MAXHF+N,L1)=N1+N-1
10    CONTINUE
      IF (ITEST.GT.0) THEN
         READ(62,REC=L1) B,OVRLAP
      ENDIF

      IF (L2.EQ.0) THEN
         RETURN
      ENDIF
      MAXHF=MAXNHF(L2)
      MAXLG=MAXNLG(L2)
      ITEST=MAXHF-MAXLG
      N1=N2+1
      N2=N1+NVARY(L2)-1

      READ(61,REC=L2) ((UJ(I,N),I=1,NPTS),N=N1,N2)
      DO 12,N=1,NVARY(L2)
         IPOS(MAXHF+N,L2)=N1+N-1
12    CONTINUE
      IF (ITEST.GT.0) THEN
         READ(62,REC=L2) B,OVRLAP
      ENDIF

      RETURN
      END
C*********************************************************************
      SUBROUTINE STGRRD (MJS)
C     THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES
C     CERTAIN VARIABLES
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      CHARACTER*8 FILE1,FILE2,FILE3
      CHARACTER*72 LVALUE(4)*1
      DATA LVALUE/'S','P','D','F'/
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB,IDIM14=MZFAC)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM9=IDIM1+IDIM1-1,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA
      COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1,
     *                JBUFF1,JBUFF2
      COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /COPY  / ICOPY1,ICOPY2,ITOTAL,IPSEUD
      COMMON /CORE  / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2),
     *                 MAXNC(IDIM2),MAXPN(IDIM2),ICHECK
      COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4
      COMMON /FACT  / GAMMA(IDIM14)
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *                NBUG9
      COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS
C
C      SET DEFAULT VALUE OF LFIXN, THE L VALUE AT WHICH THE
C      NUMBER OF CONTINUUM TERMS IS ALLOWED TO DECREASE BELOW
C      NRANG2.
      DATA LNDEF/10/
C
C     FORMAT STATEMENTS
2000  FORMAT(A)
1000  FORMAT(TR31,'***********'/TR31,'**       **'/TR31,'** NXRAD **'/
     *       TR31,'**       **'/TR31,'***********'//)
1002  FORMAT(////TR28,' SUBROUTINE STGRRD'/TR28,' -----------------')
1003  FORMAT(//' OUTPUT CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  SEQUENTIAL    UNIT NUMBER =',I3/
     *         '  (FOR MODEL POTENTIAL)'/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,
     *         ' NAME= ',A8/
     *         '  (FOR RADIAL INTEGRALS)'/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,
     *         ' NAME= ',A8/
     *         '  (FOR RK POINTER ARRAYS)'/
     *         ' OUTPUT SEQUENTIAL    UNIT NUMBER =',I3,
     *         ' NAME= ',A8/)
1004  FORMAT(' NBUG PARAMETERS'/12I5)
1005  FORMAT(' ===   NON-MODEL POTENTIAL CALCULATION   ===')
1006  FORMAT(' ===       MODEL POTENTIAL CALCULATION   ===')
1007  FORMAT(' === ANALYTIC  BOUND ORBITALS BEING USED ===')
1008  FORMAT(' === NUMERICAL BOUND ORBITALS BEING USED ===')
1009  FORMAT(/' BASIC DATA'/)
1010  FORMAT(' NELC =',I3,'  NZ =',I3,'  LRANG1 =',I3,'  NRANG2 =',I3,
     *       '  LFIXN =',I3, ' LAMAX =',I2,'  IBC =',I2)
1011  FORMAT(' MAXNHF=',35I3/)
1012  FORMAT(' MAXNLG=',35I3/)
1013  FORMAT(/' DATA DEFINING THE STATIC POTENTIAL -',
     1       ' SHELL OCCUPANCY OF THE GROUND STATE')
6014  FORMAT(1X,35(1X,I1,A1))
1014  FORMAT(1X,35I3)
1015  FORMAT(//' THE RADIAL FUNCTIONS'/TR37,' SLATER-TYPE',TR5,
     *' POWER OF R',TR5,' EXPONENT')
1016  FORMAT(/TR13,I2,A,'  ORBITAL')
1017  FORMAT(TR37,E12.5,TR9,I3,TR9,E12.5)
1018  FORMAT(//' R-MATRIX BOUNDARY CONDITIONS,  RA =',F10.5,'   BSTO =
     *',F10.5//' AMPLITUDE OF THE FUNCTIONS AT RA'/)
1019  FORMAT(//' INTEGRATION MESH'//' NIX=',I3)
1020  FORMAT(' HINT =',E14.7)
1021  FORMAT(' IHX=',20I5)
1022  FORMAT(' IRX=',20I5//)
1023  FORMAT(' MAXNC=',20I5)
1024  FORMAT(' LPOT=',I3)
1025  FORMAT(' LPOS=',20I5)
1026  FORMAT(I5,A,TR4,E14.7)
 2005 FORMAT(' *** ANGULAR AND RADIAL DATA INCONSISTENT'/
     1       '  NELC =',I3,' IN CONFIGURATIONS AND =',I3,
     2       ' IN RADIAL DATA')
C
      ETA=0.0D0
C
C     SET THE CHANNEL UNIT NUMBERS
      IREAD=5
      IWRITE=6
      JBUFF1=23
      JBUFF2=24
      ITAPE3=25
      ITAPE1=35
C
      FILE1='RAD1.DAT'
      FILE2='RAD2.DAT'
      FILE3='RAD3.DAT'
C       OPEN THE DIRECT ACCESS FILES
      OPEN(UNIT=JBUFF1,FILE=FILE1,ACCESS='DIRECT',
     *    STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF1)
      OPEN(UNIT=JBUFF2,FILE=FILE2,ACCESS='DIRECT',
     *    STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF2)
      OPEN(ITAPE3,FILE=FILE3,ACCESS='SEQUENTIAL',
     *    STATUS='UNKNOWN',FORM='UNFORMATTED')
C     READ IN THE BASIC DATA
      WRITE(IWRITE,1000)
      WRITE(IWRITE,1002)
      WRITE(IWRITE,1003)IWRITE,IREAD,ITAPE1,JBUFF1,FILE1,
     *                  JBUFF2,FILE2,ITAPE3,FILE3
      READ (IREAD, *)    NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,
     *                   NBUG8,NBUG9
      WRITE(IWRITE,1004) NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,
     *                   NBUG8,NBUG9
      READ (IREAD, *) IPSEUD,NUMERC
      IF (IPSEUD.EQ.0) THEN
         WRITE(IWRITE,1005)
      ELSE
         WRITE(IWRITE,1006)
C EFFECTIVE POTENTIAL (2*ZA/R)
        IPSEUD=IABS(IPSEUD)
      ENDIF
      IF (NUMERC.EQ.0) THEN
         WRITE(IWRITE,1007)
      ELSE
         WRITE(IWRITE,1008)
      ENDIF
      WRITE(IWRITE,1009)
      READ (IREAD, *)   NELC1,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM,IBC
      IF (NELC1 .NE. NELC) THEN
         WRITE (IWRITE,2005) NELC, NELC1
         STOP
      END IF
      IF (LRANG1 .GT. IDIM1) CALL RECOV1(1,LRANG1,IDIM1)
      IF (NRANG2 .GT. IDIM3) CALL RECOV1(3,NRANG2,IDIM3)
      IF (LAMAX  .GT. IDIM4) CALL RECOV1(4,LAMAX,IDIM4)
C
C      CONSTRUCT NVARY(L)
C
      IF (LFIXN .LT. LNDEF) THEN
         LFIXN=LNDEF
         WRITE(IWRITE,'('' LFIXN HAS BEEN RESET TO '',I3)')LFIXN
      END IF
C
      DO 1    I=1,LFIXN
         NVARY(I)=NRANG2
    1 CONTINUE
C
      DO 2    I=LFIXN+1,LRANG2
         NVARY(I)=0
    2 CONTINUE
C
      WRITE(IWRITE,1010) NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,IBC
      READ (IREAD, *)    (MAXNHF(I),I=1,LRANG1)
      WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1)
      READ (IREAD, *)    (MAXNLG(I),I=1,LRANG1)
      WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1)
C     READ IN DATA DEFINING THE STATIC POTENTIAL OF THE GROUND STATE
      WRITE(IWRITE,1013)
      NBOUND=0
      DO 10,I=1,LRANG1
         NHF=MAXNHF(I)
         NLG=MAXNLG(I)
         IF (IDIM9 .LE. NLG-I) CALL RECOV1(9,NLG-I+1,IDIM9)
         IF (IDIM10 .LT. NHF-NLG) CALL RECOV1(10,NHF-NLG,IDIM10)
         IF (IDIM12 .LT. NHF) CALL RECOV1(12,NHF,IDIM12)
         NBOUND=NBOUND+NHF-I+1
      READ (IREAD, *)  (ISTAT(I,J),J=I,MAXNHF(I))
      WRITE(IWRITE,6014) (J,LVALUE(I),J=I,MAXNHF(I))
      WRITE(IWRITE,1014) (ISTAT(I,J),J=I,MAXNHF(I))
10    CONTINUE
      IF (IDIM13 .LT. NBOUND) CALL RECOV1(13,NBOUND,IDIM13)
      IF (NUMERC.EQ.0) THEN
C        ANALYTIC ORBITALS ARE TO BE USED
         M1=0
         ZMIN=ONE/TINY
         WRITE(IWRITE,1015)
C        READ IN COEFFICIENTS DEFINING THE ANALYTIC ORBITALS
         DO 20,L=1,LRANG1
            DO 30,N=L,MAXNHF(L)
               WRITE(IWRITE,1016) N,LVALUE(L)
               M1=M1+1
               READ (IREAD, *) M
               IF (IDIM11 .LT. M) CALL RECOV1(11,M,IDIM11)
               READ (IREAD, *) (IRAD(L,N,K),K=1,M)
               READ (IREAD, *) (ZE  (L,N,K),K=1,M)
               READ (IREAD, *) (C   (L,N,K),K=1,M)
               DO 31, K=1,M
                  IR=IRAD(L,N,K)
                  IR2=IR+IR+1
                  IF (IDIM14 .LT. IR2)CALL RECOV1(14,IR2,IDIM14)
 31            CONTINUE
               NCO(M1)=M
               X=ZERO
C              CHECK THE NORMALISATION OF EACH ANALYTIC ORBITAL
               DO 40,K=1,M
                  IF (ZE(L,N,K).LT.ZMIN) THEN
                     ZMIN=ZE(L,N,K)
                     IMIN=IRAD(L,N,K)
                  ENDIF
                  X1=ZERO
                  K2=IRAD(L,N,K)+1
                  Z1=ZE(L,N,K)
                  DO 180,K1=1,M
                     X1=X1+C(L,N,K1)*GAMMA(IRAD(L,N,K1)+K2)/
     *                  (ZE(L,N,K1)+Z1)**(K2+IRAD(L,N,K1))
180               CONTINUE
                  X=X+C(L,N,K)*X1
40             CONTINUE
               IF (ABS(X-ONE).GT.PT01) THEN
C                 THE COEFFICIENTS MAY BE CLEMENTI TYPE
C                 RENORMALISE THEM
                  DO 50,K=1,M
                     IR=IRAD(L,N,K)
                     ZR=TWO*ZE(L,N,K)
                     C(L,N,K)=C(L,N,K)*SQRT(ZR/GAMMA(IR+IR+1))*ZR**IR
50                CONTINUE
C                 CHECK NORMALISATION
                  X=ZERO
                  DO 190,K=1,M
                     X1=ZERO
                     K2=IRAD(L,N,K)+1
                     Z1=ZE(L,N,K)
                     DO 182,K1=1,M
                        X1=X1+C(L,N,K1)*GAMMA(IRAD(L,N,K1)+K2)/
     *                     (ZE(L,N,K1)+Z1)**(K2+IRAD(L,N,K1))
182                  CONTINUE
                     X=X+C(L,N,K)*X1
190               CONTINUE
                  IF (ABS(X-ONE).GT.PT01) THEN
C                    ERROR IN INPUT DATA. ORBITAL NOT NORMALISED
                     PRINT*,' ORBITAL',N,LVALUE(L),'HAS NORMALISATION',X
                     STOP
                  ENDIF
               ENDIF
               Y=SQRT(X)
               DO 60,K=1,M
                  C(L,N,K)=C(L,N,K)/Y
                  WRITE(IWRITE,1017) C(L,N,K),IRAD(L,N,K),ZE(L,N,K)
60             CONTINUE
30          CONTINUE
20       CONTINUE
         IF (IBC.EQ.0) THEN
C           ESTIMATE THE BOUNDARY RADIUS RA
            ITEST=0
            TINORB=PT001
            RA=-LOG(TINORB)/ZMIN
            RA=RA+IMIN*LOG(RA)/ZMIN
            IRA=INT(RA+HALF)
220         RA=DBLE(IRA)
            BSTO=ZERO
            M1=0
            SMALL=FIFTH
C           CHECK THAT THE BOUNDARY RADIUS IS CORRECT FOR ALL
C           ORBITALS.
            DO 70,L=1,LRANG1
               DO 80,N=L,MAXNHF(L)
                  M1=M1+1
                  M=NCO(M1)
100               X=ZERO
                  DO 90,K=1,M
                     X=X+C(L,N,K)*RA**IRAD(L,N,K)*EXP(-ZE(L,N,K)*RA)
90                CONTINUE
                  IF (X.GT.TINORB) THEN
                     RA=RA+SMALL
                     ITEST=1
                     GOTO 100
                  ENDIF
80             CONTINUE
70          CONTINUE
            IF (ITEST.EQ.0) THEN
                IF (IRA.GE.1) THEN
                   IRA=IRA-1
                   GO TO 220
                ENDIF
            ENDIF
         ELSE
            READ (IREAD, *) RA,BSTO
         ENDIF
         WRITE(IWRITE,1018) RA,BSTO
         M1=0
         DO 110,L=1,LRANG1
            DO 120,N=L,MAXNHF(L)
               M1=M1+1
               X=ZERO
               DO 130,K=1,NCO(M1)
                  X=X+C(L,N,K)*RA**IRAD(L,N,K)*EXP(-ZE(L,N,K)*RA)
130            CONTINUE
               WRITE(IWRITE,1026) N,LVALUE(L),X
120         CONTINUE
110      CONTINUE
      ENDIF
      IF (IPSEUD.EQ.0) THEN
C        THIS IS A NON-MODEL POTENTIAL CALCULATION.
         IF (IBC.GT.1) THEN
C           READ IN INTEGRATION MESH
            READ (IREAD, *) NIX
            READ (IREAD, *) HINT
            READ (IREAD, *) (IHX(I),I=1,NIX)
            READ (IREAD, *) (IRX(I),I=1,NIX)
         ELSE
C           CALCULATE THE INTEGRATION MESH
            CALL MESH
         ENDIF
C        INITIALISE LPOT,LPOSX AND MAXNC IN NON-MODEL POTENTIAL
C        CALCULATION.
         LPOT=1
         DO 150,I=1,LRANG2
            MAXNC(I)=I-1
            LPOSX(I)=1
150      CONTINUE
      ELSE
C        READ IN THE MESH AND MODEL POTENTIAL FROM THE
C        SEQUENTIAL FILE ITAPE1.
         OPEN(UNIT=ITAPE1,ACCESS='SEQUENTIAL',FORM='UNFORMATTED',
     *        FILE='STG1.POT',STATUS='OLD')
         READ(ITAPE1) LRANG1,LRANG2,(MAXNC(I),I=1,LRANG1)
         READ(ITAPE1) HINT,NIX,(IHX(I),I=1,NIX),(IRX(I),I=1,NIX)
         READ(ITAPE1) LPOT,(LPOSX(I),I=1,LPOT)
         NPTS=IRX(NIX)
         DO 140,I=1,LPOT
            READ(ITAPE1) (POTHAM(I,J),J=1,NPTS)
140      CONTINUE
      ENDIF
      WRITE(IWRITE,1019) NIX
      IF (IDIM6 .LT. NIX) CALL RECOV1(6,NIX,IDIM6)
      WRITE(IWRITE,1020) HINT
      WRITE(IWRITE,1021) (IHX(I),I=1,NIX)
      WRITE(IWRITE,1022) (IRX(I),I=1,NIX)
      IF (IDIM7 .LE. IRX(NIX))CALL RECOV1(7,IRX(NIX)+1,IDIM7)
      IF (IPSEUD .NE. 0) THEN
         WRITE(IWRITE,1023) (MAXNC(I),I=1,LRANG1)
         WRITE(IWRITE,1024) LPOT
         IF (IDIM8 .LT. LPOT) CALL RECOV1(8,LPOT,IDIM8)
         WRITE(IWRITE,1025) (LPOSX(I),I=1,LPOT)
      END IF
C      READ THE MAXIMUM ENERGY FOR THE BUTTLE CORRECTION
C      AND THE PARAMETER FOR BUTTLE FITTING
C      MJS=0 FOR QUADRATIC  MJS NOT 0 FOR MJS METHOD
C      NB  QUADRATIC IS ALWAYS USED IF BSTO IS NON ZERO
C
         READ(IREAD, *)EK2MAX,MJS
         WRITE(IWRITE,'(/'' ENERGY LIMITS FOR THE BUTTLE CORRECTION''/
     1         '' MAXIMUM ENERGY = '',F14.7)')EK2MAX
C
      IF (MJS .NE. 0) THEN
         WRITE(IWRITE,'('' MJS FITTING PROCEDURE REQUESTED'')')
         IF (BSTO .NE. 0) THEN
            WRITE(IWRITE,'(''  BUT QUADRATIC FIT WILL BE USED SINCE'',
     1             '' BSTO IS NON-ZERO'')')
            MJS=0
         END IF
      ELSE
         WRITE(IWRITE,'('' QUADRATIC FITTING REQUESTED'')')
      END IF
C
C     DEFINE THE MAXPN ARRAY .  NOT USED?
      DO 160,I=1,LRANG1
         MAXPN(I)=MAXNHF(I)-MAXNC(I)+I-1
160   CONTINUE
C     EXTEND THE MAXNHF AND MAXNLG ARRAYS
      IF (LRANG2.GT.LRANG1) THEN
         DO 170,I=LRANG1+1,LRANG2
            MAXNHF(I)=I-1
            MAXNLG(I)=I-1
170      CONTINUE
      ENDIF
C     INITIALISE THE VARIABLES LAMBC AND LAMCC
      IF (LAM.GE.2) THEN
         LAMBC=LAMAX
      ELSE
          LAMBC=0
      ENDIF
      IF (LAM.GE.3) THEN
         LAMCC=LAMAX
      ELSE
          LAMCC=0
      ENDIF
C     SET UP ARRAYS FOR USE IN SIMPSONS RULE INTEGRATION
      ONE3 =ONE /THREE
      TWO3 =TWO *ONE3
      FOUR3=FOUR*ONE3
      IFI=1
      XR(1)  =ZERO
      STEP(1)=ZERO
      DO 200,I=1,NIX
         H       =HINT*DBLE(IHX(I))
         WT(IFI) =ONE3*(H+STEP(IFI))
         RSTART=XR(IFI)
         IST   =IFI+2
         IFI   =IRX(I)+1
         IFX   =0
         IF (I.GT.1) IFX=IRX(I-1)
         DO 210,J=IST,IFI,2
            XR(J-1)  =RSTART+(J-2-IFX)*H
            XR(J  )  =RSTART+(J-1-IFX)*H
            STEP(J-1)=H
            STEP(J  )=H
            WT(J-1  )=FOUR3*H
            WT(J    )=TWO3*H
210      CONTINUE
200   CONTINUE
      WT(IFI)=ONE3*H
      NPTS=IFI
      RETURN
      END
C=======================================================================
      SUBROUTINE STGRRX (MINIM,MJS)
C
C     THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES
C     CERTAIN VARIABLES USING RMATRX OUTPUT FILE
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      CHARACTER*8 FILE1,FILE2,FILE3
C      CHARACTER*72 LVALUE(4)*1
C      DATA LVALUE/'S','P','D','F'/
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,
     *  IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM,
     *  IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM1=IDIM13+2*IDIM3,
     *          KDIM9=IDIM1+IDIM1-1,KDIM10=IDIM7*2,
     *          KDIM12=IDIM1*IDIM12)
      PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
C               KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3)
      PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9)
      PARAMETER(P75=0.75D0,P125=0.125D0,SIXT=16.0D0,TNINE=9.0D0)
      COMMON /BNDORB/ P(IDIM9,KDIM10)
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA
      COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1,
     *                JBUFF1,JBUFF2
      COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN
      COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC,
     *                TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE,
     *                HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH
      COMMON /COPY  / ICOPY1,ICOPY2,ITOTAL,IPSEUD
      COMMON /CORE  / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2),
     *                 MAXNC(IDIM2),MAXPN(IDIM2),ICHECK
      COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4
      COMMON /INFORM/ IREAD,IWRITE,IPUNCH
      COMMON /INIT  / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH
      COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM
      COMMON /NBUG  / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     *  	      NBUG9
      COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC
      COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2)
      COMMON /POTVAL/ POVALU(KDIM10)
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
      COMMON /SIMP  / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS
      COMMON /NRB/MAXC,LFIXN0
      DIMENSION PTEMP(IDIM7),UTEMP(IDIM7),IRXP(IDIM6)
C
C      SET DEFAULT VALUE OF LFIXN, THE L VALUE ABOVE WHICH THE
C      NUMBER OF CONTINUUM TERMS IS ALLOWED TO DECREASE BELOW
C      NRANG2.
      DATA LNDEF/10/
C
C     FORMAT STATEMENTS
 5002 FORMAT ('  ', 4E16.7)
2000  FORMAT(A)
1000  FORMAT(TR31,'***********'/TR31,'**       **'/TR31,'** NXRAD **'/
     *       TR31,'**       **'/TR31,'***********'//)
1002  FORMAT(////TR28,' SUBROUTINE STGRRX'/TR28,' -----------------')
1003  FORMAT(//' OUTPUT CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  SEQUENTIAL    UNIT NUMBER =',I3/
     *         '  (FOR MODEL POTENTIAL)'/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,
     *         ' NAME= ',A8/
     *         '  (FOR RADIAL INTEGRALS)'/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,
     *         ' NAME= ',A8/
     *         '  (FOR RK POINTER ARRAYS)'/
     *         ' OUTPUT SEQUENTIAL    UNIT NUMBER =',I3,
     *         ' NAME= ',A8/)
1004  FORMAT(' NBUG PARAMETERS'/12I5)
1005  FORMAT(' ===   NON-MODEL POTENTIAL CALCULATION   ===')
1006  FORMAT(' ===       MODEL POTENTIAL CALCULATION   ===')
1007  FORMAT(' === ANALYTIC  BOUND ORBITALS BEING USED ===')
1008  FORMAT(' === NUMERICAL BOUND ORBITALS BEING USED ===')
1009  FORMAT(/' BASIC DATA'/)
1010  FORMAT(' NELC =',I3,'  NZ =',I3,'  LRANG1 =',I3,'  NRANG2 =',I3,
     *       '  LFIXN =',I3, ' LAMAX =',I2)
1011  FORMAT(' MAXNHF=',35I3/)
1012  FORMAT(' MAXNLG=',35I3/)
1013  FORMAT(/' DATA DEFINING THE STATIC POTENTIAL')
1014  FORMAT(' L =',I3/(' ',35I3))
1015  FORMAT(//' THE RADIAL FUNCTIONS'/TR37,' SLATER-TYPE',TR5,
     *' POWER OF R',TR5,' EXPONENT')
1016  FORMAT(/TR13,I2,A,'  ORBITAL')
1017  FORMAT(TR37,E12.5,TR9,I3,TR9,E12.5)
1018  FORMAT(//' R-MATRIX BOUNDARY CONDITIONS,  RA =',F10.5,'   BSTO =
     *',F10.5//' AMPLITUDE OF THE FUNCTIONS AT RA'/)
1019  FORMAT(//' INTEGRATION MESH'//' NIX=',I3)
1020  FORMAT(' HINT =',E14.7)
1021  FORMAT(' IHX=',20I5)
1022  FORMAT(' IRX=',20I5//)
1023  FORMAT(' MAXNC=',20I5)
1024  FORMAT(' LPOT=',I3)
1025  FORMAT(' LPOS=',20I5)
1026  FORMAT(I5,A,TR4,E14.7)
C
      ETA=0.0D0
C
C     SET THE CHANNEL UNIT NUMBERS
      IREAD=5
      IWRITE=6
      JBUFF1=23
      JBUFF2=24
      ITAPE3=25
      INX1=35
C
      FILE1='RAD1.DAT'
      FILE2='RAD2.DAT'
      FILE3='RAD3.DAT'
C       OPEN THE DIRECT ACCESS FILES
      OPEN(UNIT=JBUFF1,FILE=FILE1,ACCESS='DIRECT',
     1    STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF1)
      OPEN(UNIT=JBUFF2,FILE=FILE2,ACCESS='DIRECT',
     1    STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF2)
      OPEN(ITAPE3,FILE=FILE3,ACCESS='SEQUENTIAL',
     1    STATUS='UNKNOWN',FORM='UNFORMATTED')
C
      OPEN (INX1,FILE='NX1.DAT',ACCESS='SEQUENTIAL',STATUS='OLD',
     1     FORM='UNFORMATTED')
C     OPEN (INX1,FILE='tapenx1',ACCESS='SEQUENTIAL',STATUS='OLD',
C
      WRITE(IWRITE,1000)
      WRITE(IWRITE,1002)
      WRITE(IWRITE,1003)IWRITE,IREAD,ITAPE1,JBUFF1,FILE1,
     *                  JBUFF2,FILE2,ITAPE3,FILE3
      WRITE (IWRITE,1004)NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,
     1                   NBUG8,NBUG9
C
      WRITE(IWRITE,1009)
      READ (INX1) LRANG1,NRANG2,LAMAX,NIX,NPTS,NRKPTS
      READ (INX1) NELC,NZ,LFIXN,IPSEUD
C EFFECTIVE CHARGE
      IPSEUD=-IABS(IPSEUD)
      IF(MAXC.GT.1)NRANG2=MAXC
      IF (LRANG1 .GT. IDIM1) CALL RECOV1(1,LRANG1,IDIM1)
      IF (NRANG2 .GT. IDIM3) CALL RECOV1(3,NRANG2,IDIM3)
      IF (LAMAX  .GT. IDIM4) CALL RECOV1(4,LAMAX,IDIM4)
         IF (MINIM .LE. LRANG1) THEN
            WRITE (IWRITE,'('' *** NOEXCHANGE PROGRAM NOT VALID'', /
     1           '' *** BOUND TERMS ARE BEING IGNORED ***'')')
         END IF
C
C      OPEN FILE TO STORE BOUND ORBITALS REQUIRED IN LAGRANGE
C      ORTHOGONALISATION
C
      IF (LRANG1 .GE. MINIM) THEN
         OPEN(60,ACCESS='DIRECT',STATUS='SCRATCH',
     1        FORM='UNFORMATTED',RECL=MACDIM*IDIM9*KDIM10)
      END IF
C
C      CONSTRUCT NVARY(L)
C
      IF (LFIXN .LT. LNDEF) THEN
         LFIXN=LNDEF
         IF(LFIXN0.GT.0)LFIXN=LFIXN0
         WRITE(IWRITE,'('' LFIXN HAS BEEN RESET TO '',I3)')LFIXN
      END IF
C
      IF(MAXC.GT.1.AND.LFIXN.LT.MINIM)LFIXN=MINIM
      IF(LFIXN0.GT.0)LFIXN=LFIXN0
C
      DO 1    I=1,LFIXN
         NVARY(I)=NRANG2
    1 CONTINUE
C
      DO 2    I=LFIXN+1,LRANG2
         NVARY(I)=0
    2 CONTINUE
C
      WRITE(IWRITE,1010) NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX
      NUMERC=1
      READ (INX1) (MAXNHF(I),I=1,LRANG1), (MAXNLG(I),I=1,LRANG1)
      WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1)
      WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1)
      NBOUND=0
      DO 10,I=1,LRANG1
         NHF=MAXNHF(I)
         NLG=MAXNLG(I)
         IF (IDIM9 .LE. NLG-I) CALL RECOV1(9,NLG-I+1,IDIM9)
         IF (IDIM10 .LT. NHF-NLG) CALL RECOV1(10,NHF-NLG,IDIM10)
         IF (IDIM12 .LT. NHF) CALL RECOV1(12,NHF,IDIM12)
         NBOUND=NBOUND+NHF-I+1
10    CONTINUE
      IF (IDIM13 .LT. NBOUND) CALL RECOV1(13,NBOUND,IDIM13)
C
      READ (INX1) HINT,(IHX(I),I=1,NIX),(IRX(I),I=1,NIX)
      WRITE(IWRITE,1019) NIX
      IF (IDIM6 .LT. NIX) CALL RECOV1(6,NIX,IDIM6)
      WRITE(IWRITE,1020) HINT
      WRITE(IWRITE,1021) (IHX(I),I=1,NIX)
      WRITE(IWRITE,1022) (IRX(I),I=1,NIX)
      IF (IDIM7 .LE. IRX(NIX))CALL RECOV1(7,IRX(NIX)+1,IDIM7)
      LPOT=1
      DO 150  I=1,LRANG2
         MAXNC(I)=I-1
        LPOSX(I)=1
 150  CONTINUE
      IF(IPSEUD.NE.0)THEN
        READ(INX1) (MAXNC(I),I=1,LRANG1)
        READ(INX1) LPOT,(LPOSX(I),I=1,LPOT)
        DO 155 I=1,LPOT
          READ(INX1)(POTHAM(I,J),J=1,IRX(NIX))
 155    CONTINUE
      ENDIF
C
         READ (INX1) RA,BSTO,EK2MAX
         WRITE (IWRITE,1018) RA,BSTO
      IF(MAXC.GT.1)EK2MAX=ZERO
         WRITE(IWRITE,'(/'' ENERGY LIMITS FOR THE BUTTLE CORRECTION''/
     1         '' MAXIMUM ENERGY = '',F14.7)')EK2MAX
C
      IF (MJS .NE. 0) THEN
         WRITE(IWRITE,'('' MJS FITTING PROCEDURE REQUESTED'')')
         IF (BSTO .NE. 0) THEN
            WRITE(IWRITE,'(''  BUT QUADRATIC FIT WILL BE USED SINCE'',
     1             '' BSTO IS NON-ZERO'')')
            MJS=0
         END IF
      ELSE
         WRITE(IWRITE,'('' QUADRATIC FITTING REQUESTED'')')
      END IF
C
C     DEFINE THE MAXPN ARRAY .   NOT USED?
      DO 160,I=1,LRANG1
         MAXPN(I)=MAXNHF(I)-MAXNC(I)+I-1
160   CONTINUE
C     EXTEND THE MAXNHF AND MAXNLG ARRAYS
      IF (LRANG2.GT.LRANG1) THEN
         DO 170,I=LRANG1+1,LRANG2
            MAXNHF(I)=I-1
            MAXNLG(I)=I-1
170      CONTINUE
      ENDIF
C
C      READ IN BOUND ORBITALS
C
      NQ=0
      NPTS=IRX(NIX)
      DO 301  L=1,LRANG1
         MAXHF=MAXNHF(L)
         DO 302  N=L,MAXHF
            NQ=NQ+1
            IPOS(N,L)=NQ
            READ (INX1) (UJ(I,NQ),I=1,NPTS+1)
 302     CONTINUE
 301  CONTINUE
C
C      READ IN POTENTIAL
C
      READ (INX1) (POVALU(I),I=1,NPTS*2)
C
C     SET UP ARRAYS FOR USE IN SIMPSONS RULE INTEGRATION
      ONE3 =ONE /THREE
      TWO3 =TWO *ONE3
      FOUR3=FOUR*ONE3
      IFI=1
      XR(1)  =ZERO
      STEP(1)=ZERO
      DO 200,I=1,NIX
         H       =HINT*DBLE(IHX(I))
         WT(IFI) =ONE3*(H+STEP(IFI))
         RSTART=XR(IFI)
         IST   =IFI+2
         IFI   =IRX(I)+1
         IFX   =0
         IF (I.GT.1) IFX=IRX(I-1)
         DO 210,J=IST,IFI,2
            XR(J-1)  =RSTART+(J-2-IFX)*H
            XR(J  )  =RSTART+(J-1-IFX)*H
            STEP(J-1)=H
            STEP(J  )=H
            WT(J-1  )=FOUR3*H
            WT(J    )=TWO3*H
210      CONTINUE
200   CONTINUE
      WT(IFI)=ONE3*H
      NPTS=IFI
C
C      INTERPOLATE THE BOUND ORBITALS AND STORE ON SCRATCH FILE
C      READY FOR LAGRANGE ORTHOGONALISATION
C
      A=TNINE/SIXT
      B=ONE /SIXT
      DO 401  L=1,LRANG1
         MAXLG=MAXNLG(L)
         IF (L .GE. MINIM) THEN
C
C      TAKE THE FACTOR R**L OUTSIDE FOR THE FIRST IFIN POINTS
C
            IF (L .EQ. 1) THEN
               NIX1=1
               IFIN=4
               IRXP(1)=3
               GO TO 403
            END IF
            XRFIN=L*(L-1)
            XRFIN=XRFIN/(2.D0*NZ)
            IR1=IRX(1)+1
            IF (XRFIN .LT. XR(IR1)) THEN
               NIX1=1
               IFIN=NINT(XRFIN/HINT)
               IF (IFIN .LT. 4) IFIN=4
               IRXP(1)=IFIN-1
               GO TO 403
            END IF
            IRXP(1)=IRX(1)
            DO 402  NI=2,NIX
               IR2=IRX(NI)+1
               IF (XRFIN .LE. XR(IR2)) THEN
                  NIX1=NI
                  IFIN=NINT((XRFIN-XR(IR1))/(HINT*IHX(NI)))
                  IF (IFIN .LT. 3) IFIN=3
                  IFIN=IFIN+IR1
                  IRXP(NI)=IFIN-1
                  GO TO 403
               END IF
               IRXP(NI)=IRX(NI)
               IR1=IR2
  402       CONTINUE
  403       KT=0
            DO 404  N=L,MAXLG
               NQ=IPOS(N,L)
               KT=KT+1
               DO 405  I=2,IFIN
                  UTEMP(I)=UJ(I,NQ)/(XR(I)**L)
  405          CONTINUE
               H=HALF*HINT
               PTEMP(1)=HALF*(3*UTEMP(2)-UTEMP(3))*H**L
               PTEMP(2)=(P75*UTEMP(3)+P125*(3*UTEMP(2)-UTEMP(4)))
     1                  *(H+HINT)**L
               IX=3
               DO 406 I=1,NIX1
                  H=HALF*HINT*IHX(I)
                  DO 407 J=IX,IRXP(I)-1
                     PTEMP(J)=(A*(UTEMP(J)+UTEMP(J+1))-B*(
     1                   UTEMP(J-1)+UTEMP(J+2)))
     2                   *(XR(J)+H)**L
 407              CONTINUE
                  IF (I.NE.NIX1) THEN
C           CHANGE OF INTERVAL...DOUBLING OF STEP SIZE IS ASSUMED
                     K2=IRX(I)
                     K1=K2+1
                     PTEMP(K1)=(A*(UTEMP(K1)+UTEMP(K1+1))-B*(
     1                    UTEMP(K1-2)+UTEMP(K1+2)))
                     PTEMP(K2)=(A*(UTEMP(K2)+UTEMP(K2+1))-B*(
     1                    PTEMP(K1)+UTEMP(K2-1)))
     2                    *(XR(K2)+H)**L
                     PTEMP(K1)=PTEMP(K1)*(XR(K1)+2*H)**L
                     IX=K1+1
                  ELSE
                     IX=IRXP(NIX1)
                  END IF
  406          CONTINUE
               DO 408  I=IX-1,IRX(NIX)+1
                  UTEMP(I)=UJ(I,NQ)
  408          CONTINUE
               DO 409  I=NIX1,NIX
                  H=HALF*HINT*IHX(I)
                  DO 410 J=IX,IRX(I)-1
                     PTEMP(J)=A*(UTEMP(J)+UTEMP(J+1))-B*(
     1                   UTEMP(J-1)+UTEMP(J+2))
 410              CONTINUE
                  IF (I.NE.NIX) THEN
C           CHANGE OF INTERVAL...DOUBLING OF STEP SIZE IS ASSUMED
                     K2=IRX(I)
                     K1=K2+1
                     PTEMP(K1)=A*(UTEMP(K1)+UTEMP(K1+1))-B*(
     1                    UTEMP(K1-2)+UTEMP(K1+2))
                     PTEMP(K2)=A*(UTEMP(K2)+UTEMP(K2+1))-B*(
     1                    PTEMP(K1)+UTEMP(K2-1))
                     IX=K1+1
                  ELSE
                     K1=IRX(NIX)
                     PTEMP(K1)=P75*UTEMP(K1) +
     1                   P125*(3*UTEMP(K1+1)-UTEMP(K1-1))
                  END IF
  409          CONTINUE
               DO 411 I=1,NPTS-1
                  P(KT,I+I-1)=PTEMP(I)
                  P(KT,I+I)=UJ(I+1,NQ)
  411          CONTINUE
  404       CONTINUE
C
            IF (KT .GT. 0) WRITE(60,REC=L)P
C
         END IF
  401 CONTINUE
      REWIND (INX1)
      RETURN
      END
C***********************************************************************
      SUBROUTINE WRITAP
C
C      WRITES BASIC INFORMATION ONTO THE STG1R PERMANENT OUTPUT FILE
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C

      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,
     *  IDIM11=MZSLT,IDIM12=MZNR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM12=IDIM1*IDIM12)
C
      COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN,
     *                LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN
      COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA
      COMMON /BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN
      COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4
      COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11),
     *                IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12),
     *                LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2),
     *                NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM
C
C     IF(ITOTAL .LE. 0) RETURN
      ITAPE=ITAPE3
      REWIND ITAPE
      WRITE(ITAPE)RA,BSTO
      WRITE(ITAPE)(NVARY(L),L=1,LRANG2)
      WRITE(ITAPE)((ENDS(N,L),N=1,NVARY(L)),L=1,LRANG2)
      WRITE(ITAPE)((COEFF(I,L),I=1,3),L=1,LRANG2)
C
      RETURN
      END