C N. R. BADNELL  UoS v1.11   -   QUB vx.x                      29/04/05
C
C     PROGRAM NXHAM/STG3NX
C
C *** THIS IS THE HAMILTONIAN STAGE OF RMATRX NX
C
C     VERSION WHICH ALLOWS NRANG2 TO DIMINISH WITH INCREASE
C     IN LARGE L USING THE ARRAY NVARY
C
      PROGRAM MAIN
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C SUN
      REAL*4 TARRY(2),TIME
C
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO)
      PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM6=(IDIM4+1)/2)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ICDIM3=ICDIM2+ICDIM2-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(ILDIM4=ILDIM3+ILDIM3)
      PARAMETER(INDIM2=IDIM3*IDIM3)
      PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/ANGBUF/BUFV1(LBUFFV),BUFV2(LBUFFV),
     1              IBUFV1(LBUFFV),IBUFV2(LBUFFV)
      COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4),
     1              IRHSGL(IVDIM1)
      COMMON/ASYMPR/CRL(ITDIM1,ITDIM1,IDIM4)
      COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2),
     1              ICONCH(ITDIM3),KCONCH(ITDIM3),
     2              L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1),
     3              NSCHAN(IPDIM1),NSPREC(IPDIM1)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/HPOINT/KRECH(0:ILDIM4)
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RADBUF/BUFF1(LBUFF1),BUFF2(LBUFF1)
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON/SETL  /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2),
     1              NSETL(ISDIM1,ILDIM2)
      COMMON/SETS  /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1              NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2)
      COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2)
      COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
      COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1)
      COMMON/STG1  /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2)
C
      COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN
C
 1000 FORMAT(TR31,' END OF NXHAM'/TR31,' ------------'///)
C
      CALL RECSET(LRGLLO,LRGLUP,ICONTN)
      IF (ICONTN .EQ. 1) THEN
         CALL STGRRX
      ELSE
      CALL STGRRD
      END IF
C
      CALL STGHRD(ICONTN)
C
C      READ STG1 TAPE
C
      CALL RETAP1
C
      CALL CRLMAT
C
C      WRITE ASYMPTOTIC FILE
C
      CALL WRITE1
C
      DO 10  LRGL=LRGLLO,LRGLUP
      DO 11  NPTY0=NPTYMN,NPTYMX
      IF(NPTYIN.LT.0)THEN
      NPTY=NPTY0
      ELSE
      NPTY=MOD(LRGL,2)+NPTY0
      NPTY=MOD(NPTY,2)
      ENDIF
        CALL CUSETL
        CALL SETCHN(LRGLLO)
        CALL STHMAT(LRGLLO)
      DO 6  ISP=1,NSP
C
         IF (LRGL .EQ. LRGLUP .AND.NPTY0.EQ.NPTYMX.AND. ISP .EQ. NSP)
     1      THEN
            MORE=0
         ELSE
            MORE= 2
         END IF
C
C
C      READ THE ASYMPTOTIC COEFFICIENTS AND THE HAMILTONIAN MATRIX
C      FROM DIRECT ACCESS FILE JBUFD.
C      CALCULATE ENERGY LEVELS AND SURFACE AMPLITUDES.
C      WRITE THE ASYMPTOTIC FILE BLOCKS FOR THIS LRGL,NPTY
C      AND TARGET SPIN
C
         CALL RDHMAT(ISP,MORE)
    6 CONTINUE
   11 CONTINUE
   10 CONTINUE
C
      CLOSE (JBUFIR)
      CLOSE (JBUFFV)
      CLOSE (JBUFFZ)
      CLOSE (JBUFF1)
      CLOSE (JBUFF2)
      CLOSE (JBUFD)
      CLOSE (ITAPE1)
      ENDFILE ITAPE3
      REWIND ITAPE3
      CLOSE (ITAPE3)
C
      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
**********************************************************************
      BLOCK DATA
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
C
C     SET GLOBAL REAL CONSTANTS
C
      DATA ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS/
     1     0.0D 00,0.1D 00,0.5D 00,1.0D 00,2.0D 00,3.0D 00,4.0D 00,
     2   7.0D 00,1.1D 01,1.0D-09/
C
      END
C***********************************************************************
      SUBROUTINE CMATII(L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT)
C
C      TO CALCULATE THE CONTRIBUTION TO THE HAMILTONIAN FROM THE
C      INTERACTION OF A PAIR OF SIMILAR CONFIGURATIONS
C      THIS SUBMATRIX IS STORED IN CMAT
C      ZB IS THE ARRAY OF Z-COEFFICIENTS FOR CURRENT L,LP FOR
C      LAMBDA=LAMST TO LAMFIN
C      ICOUNT IS THE COUNTER TO LOCATE VSH FOR THIS CONFIG PAIR
C      AT EXIT ICOUNT HOLDS THE COUNT FOR THE NEXT CONFIGURATION PAIR
C      LAMCFG IS THE MAX. LAMBDA FOR WHICH THERE IS A VSH VALUE
C      FOR THIS CONFIG PAIR
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(ILDIM4=ILDIM3+ILDIM3)
      PARAMETER(INDIM2=IDIM3*IDIM3)
      PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4),
     1              IRHSGL(IVDIM1)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2)
C
      DIMENSION ZB(ILDIM1)
      DIMENSION CMAT(IDIM3,IDIM3),VRMAT(IDIM3,IDIM3),RM(INDIM2)
C
 1000 FORMAT(/,28X,'SUBROUTINE CMATII'/28X,'-----------------'/)
 1001 FORMAT('    LAMBDA VRMAT')
 1002 FORMAT('    LAMBDA  CMAT')
 1003 FORMAT(4X,I6,/2X,10E11.4,/2X,10E11.4)
C
      IF (NBUG1 .EQ. 1) THEN
      WRITE(IWRITE,1000)
      END IF
C
C    ZEROISE THE MATRIX CMAT(N,NP) READY TO ACCUMULATE THE ANGULAR TERMS
C
      DO 1     N=1,IDIM3
      DO 2    NP=1,IDIM3
          CMAT(NP,N)=ZERO
    2 CONTINUE
    1 CONTINUE
C
C      SET LIMITS FOR N AND NP
C
      NRA1=NVARY(L+1)
      NRA2=NVARY(LP+1)
C
C      SET LAMBDA LIMITS CORRESPONDING TO VSHELL FILE
C
         LAMLO=MOD (ABS (IRHSGL(ICOUNT)), 100)
         LAMUP=LAMCFG
C
      DO 3    LAM=LAMLO,LAMUP,2
           NSHSGL=ABS( IRHSGL(ICOUNT))
C
C      FIND THE NUMBER OF SHELLS CONTRIBUTING FOR THIS LAMBDA
C
              NSH=NSHSGL/10000
           IF(NSH .EQ. 0) THEN
C
C      IF NO SHELLS CONTRIBUTE HIGHER VALUES OF LAMBDA ALSO
C      CANNOT CONTRIBUTE SO RETURN
C
           ICOUNT=ICOUNT+1
             RETURN
           END IF
C
           IF(LAM .LT. LAMST  .OR.  LAM .GT. LAMFIN) THEN
C
C      SKIP TO NEXT LAMBDA LOCATION OVER THE NSH VALUES
C
           ICOUNT=ICOUNT+NSH
             GO TO 3
           END IF
C
C      ZEROISE VRMAT(N,NP) TO ACCUMULATE CONTRIBUTION FROM SHELLS
C      FOR THIS LAMBDA VALUE
C
      DO 4  N=1,IDIM3
      DO 5 NP=1,IDIM3
         VRMAT(NP,N)=ZERO
    5 CONTINUE
    4 CONTINUE
      NSHHUN=NSH*100
C
C      SUM OVER SHELLS
C
      DO 6      ISH=1,NSH
               ISIG=ABS( IRHSGL(ICOUNT))/100 - NSHHUN
                 VS=VSH(ICOUNT)
             ICOUNT=ICOUNT+1
         IF(ABS(VS) .LT. EPS) GO TO 6
               LSIG=LJCOMP(ISIG)
               NSIG=NJCOMP(ISIG)
C
C      LOCATE THE RADIAL INTEGRAL
C
      CALL LOCCCD(LSIG+1, LSIG+1, NSIG, NSIG, LAM+1, RM)
      IF (NBUG3 .EQ. 1) THEN
         WRITE(IWRITE,'(/''RK INTEGRAL FOR L='',I4,''LP='',I4,
     1         ''LAM='',I4,''L1='',I4,''L1P='',I4,''N1='',I4,''N1P='',
     2           I4)') L,LP,LAM,LSIG,LSIG,NSIG,NSIG
         WRITE(IWRITE,'(2X,8E14.7)') (RM(I),I=1,ICCLNG)
      END IF
C
C
                  I=0
      DO 7        N=1,NRA1
         IF(L .EQ. LP) THEN
           DO 71     NP=1,N
              I=I+1
              VRMAT(N,NP)=VRMAT(N,NP)+RM(I)*VS
   71      CONTINUE
C
         ELSE
C
           DO 72     NP=1,NRA2
              I=I+1
              VRMAT(NP,N)=VRMAT(NP,N)+RM(I)*VS
  72       CONTINUE
C
         END IF
C
    7 CONTINUE
C
    6 CONTINUE
      IF (NBUG1 .EQ. 1) THEN
      WRITE(IWRITE,1001)
      WRITE(IWRITE,1003)LAM,((VRMAT(NP,N),NP=1,NRA2),N=1,NRA1)
      END IF
C
C      MOVE THE ZBUFF COUNTER TO LOCATE THE Z COEFFICIENT
C      FOR THIS LAMBDA VALUE
C
                IZCONT=1+(LAM-LAMST)/2
                     Z=ZB(IZCONT)
C
C      ADD Z*VR INTO CMAT(N,NP)
C
      DO 9     N=1,NRA1
         IF(L .EQ. LP) THEN
            NPLO=N
         ELSE
            NPLO=1
         END IF
C
      DO 10   NP=NPLO,NRA2
         CMAT(NP,N)=CMAT(NP,N)+Z*VRMAT(NP,N)
   10 CONTINUE
C
    9 CONTINUE
      IF (NBUG1 .EQ. 1) THEN
      WRITE(IWRITE,1002)
      WRITE(IWRITE,1003)LAM,((CMAT(NP,N),NP=1,NRA2),N=1,NRA1)
      END IF
C
    3 CONTINUE
C
      RETURN
      END
**********************************************************************
      SUBROUTINE CMATIJ(
     1           L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT,IND)
C
C      TO CALCULATE THE CONTRIBUTION TO THE HAMILTONIAN MATRIX
C      FROM THE INTERACTION OF A PAIR OF CONFIGURATIONS WHICH ARE
C      NOT IDENTICAL
C      THIS SUB MATRIX IS STORED IN CMAT
C      ZB IS THE ARRAY OF Z-COEFFICIENTS FOR CURRENT L,LP FOR
C      LAMCFG IS THE MAX. LAMBDA ALLOWED BY TARGET CONFIGURATION
C      ANGULAR MOMENTUM
C      ICOUNT IS THE COUNTER TO LOCATE VSHELL FOR THIS CONFIG PAIR
C      AT EXIT ICOUNT HOLDS THE COUNT FOR THE NEXT CONFIGURATION PAIR
C      IND RETURNS 1 IF IC,JC COUPLED
C                  0 IF THEY DID NOT
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(ILDIM4=ILDIM3+ILDIM3)
      PARAMETER(INDIM2=IDIM3*IDIM3)
      PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4),
     1              IRHSGL(IVDIM1)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2)
C
C      IND INITIALLY SET ZERO
C
      DIMENSION ZB(ILDIM1)
      DIMENSION CMAT(IDIM3,IDIM3),RM(INDIM2)
C
 1000 FORMAT(/,28X,'SUBROUTINE CMATIJ'/28X,'-----------------'/)
 1002 FORMAT('    LAMBDA  CMAT')
 1003 FORMAT(4X,I6,/2X,10E11.4,/2X,10E11.4)
C
      IF (NBUG1 .EQ. 1) THEN
      WRITE(IWRITE,1000)
      END IF
C
      IND=0
C
      IRSL=IRHSGL(ICOUNT)
       IF(IRSL .EQ. 9999) THEN
         ICOUNT=ICOUNT+1
         RETURN
       END IF
C
C      FIND THE INTERACTING SHELLS AND THE FIRST LAMBDA FROM IRSL
C
      LAMLO=MOD(IRSL,100)
        IRS=IRSL/100
       IRHO=IRS/100
       ISIG=MOD(IRS,100)
       LRHO=LJCOMP(IRHO)
       LSIG=LJCOMP(ISIG)
      LAMUP=MIN(LRHO+LSIG,LAMCFG)
C
C      TEST WHETHER THESE LIMITS OVERLAP THE LIMITS ON LAMBDA
C      FOR THE Z COEFFICIENTS
C
       IF(LAMLO .GT. LAMFIN  .OR.
     1    LAMUP .LT. LAMST) THEN
         ICOUNT=ICOUNT+1+(LAMUP-LAMLO)/2
         RETURN
       END IF
C
       NRHO=NJCOMP(IRHO)
       NSIG=NJCOMP(ISIG)
C
C      ZEROISE THE CMAT MATRIX READY TO SUM OVER LAMBDA
C
      DO 1     N=1,IDIM3
      DO 2    NP=1,IDIM3
          CMAT(NP,N)=ZERO
    2 CONTINUE
    1 CONTINUE
C
C      SET LIMITS ON N AND NP
C
      NRA1=NVARY(L+1)
      NRA2=NVARY(LP+1)
C
C      SET IND
C
      IND=1
C
      DO 3     LAM=LAMLO,LAMUP,2

          IF (LAM .LT. LAMST  .OR.  LAM .GT. LAMFIN) THEN
C
C      MOVE ALONG VSHELL FILE
C
            ICOUNT=ICOUNT+1
            GO TO 3
          END IF
C
                VS=VSH(ICOUNT)
            ICOUNT=ICOUNT+1
          IF(ABS(VS) .LE. EPS) THEN
            GO TO 3
          END IF
C
C      MOVE ZBUFF COUNTER TO LOCATE Z FOR THIS LAMBDA
C
            IZCONT=1+(LAM-LAMST)/2
                VZ=ZB(IZCONT)*VS
C
C      LOCATE THE RADIAL INTEGRAL
C
      IF (LRHO .LT. LSIG) THEN
         CALL LOCCCD(LRHO+1, LSIG+1, NRHO, NSIG, LAM+1, RM)
      ELSE
         CALL LOCCCD(LSIG+1, LRHO+1, NSIG, NRHO, LAM+1, RM)
      END IF
C
C      RM(I) CONTAINS THE NRA1 X NRA2 MATRIX OF RADIAL
C      INTEGRALS.
C      WHEN L=LP ONLY THE LOWER HALF OF THE MATRIX IS STORED
C      RM(I) IS NOW MULTIPLIED BY THE FACTOR VZ AND STORED IN CMAT(N,NP)
C      WHEN L=LP IT IS STORED IN THE UPPER HALF OF CMAT
C
      IF (NBUG3 .EQ. 1) THEN
         WRITE(IWRITE,'(/''RK INTEGRAL FOR L='',I4,''LP='',I4,
     1         ''LAM='',I4,''L1='',I4,''L1P='',I4,''N1='',I4,''N1P='',
     2           I4)') L,LP,LAM,LRHO,LSIG,NRHO,NSIG
         WRITE(IWRITE,'(2X,8E14.7)') (RM(I),I=1,ICCLNG)
      END IF
C
                 I=0
      DO 4       N=1,NRA1
          IF(L .EQ. LP) THEN
              NPUP=N
          ELSE
              NPUP=NRA2
          END IF
C
      DO 5      NP=1,NPUP
                 I=I+1
          IF(L .EQ. LP) THEN
            CMAT(N,NP)=CMAT(N,NP)+RM(I)*VZ
          ELSE
            CMAT(NP,N)=CMAT(NP,N)+RM(I)*VZ
          END IF
C
    5 CONTINUE
C
    4 CONTINUE
C
      IF (NBUG1 .EQ. 1) THEN
      WRITE(IWRITE,1002)
      WRITE(IWRITE,1003)LAM,((CMAT(NP,N),NP=1,NRA2),N=1,NRA1)
      END IF
    3 CONTINUE
C
      RETURN
      END
**********************************************************************
      SUBROUTINE CRLMAT
C
C      CALCULATES THE REDUCED MATRIX WHOSE
C      ELEMENTS (MULTIPLIED BY A Z COEFFICIENT)
C      WILL BE USED TO CALCULATE THE ASYMPTOTIC
C      COEFFICIENTS
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM6=(IDIM4+1)/2)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(ILDIM4=ILDIM3+ILDIM3)
      PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4),
     1              IRHSGL(IVDIM1)
      COMMON/ASYMPR/CRL(ITDIM1,ITDIM1,IDIM4)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON/SETS  /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1              NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2)
      COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
      COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1)
      COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2)
C
      DIMENSION IBBPOL(IDIM1,IDIM1,KDIM6), CRLADD(IDIM4)
C
 1000 FORMAT(/,28X,'SUBROUTINE CRLMAT'/28X,'-----------------',
     1       //'     IS   JS  ICS  JCS  LAM   CRLADD')
 1001 FORMAT('  ',5I5,E14.7)
 1002 FORMAT('  CRL MATRIX  LAMBDA=',I2)
 1003 FORMAT('  ',8E14.7)
C
 3000 FORMAT('  **** STOP IN CRLMAT ****'/
     1       '  ERROR IN USING VSH FILE FOR SETS ', 2I4)
 3001 FORMAT('  **** STOP IN CRLMAT ****'/
     1       '  WRONG VSH FILE LOADED')
C
      IF (NBUG2 .GT. 1) THEN
      WRITE(IWRITE,1000)
      END IF
      LAMIND=(LAMAX+1)/2
C
C      READ FIRST RADIAL INTEGRAL POINTER BLOCK
C
      READ (ITAPE1) (((IBBPOL(I,J,K),I=1,LRANG1),J=1,LRANG1),
     1                     K=1,LAMIND),NBBPOL,(IST3(I),I=1,LRANG2)
C
C      ZEROISE CRL MATRIX
C
      DO 100  IST=1,NAST
      DO 101  JST=1,NAST
      DO 102   LM=1,LAMAX
         CRL(IST,JST,LM)=ZERO
  102 CONTINUE
  101 CONTINUE
  100 CONTINUE
C
      NRECR=0
      NB=0
C
C      READ BOUND BOUND MULTIPOLE INTEGRALS INTO RKCCD ARRAY
C      THIS SAVES SPACE
C      THESE ARE THE FIRST INTEGRALS STORED SO THEY START AT
C      THE BEGINNING OF THE SECOND RECORD
C
      ISTART=LBUFF1*2 + 1
      IFIN=ISTART + NBBPOL - 1
      CALL RDRINT(ISTART,IFIN,NB,NRECR)
C
C      READ ANGULAR INTEGRAL POINTERS AND SET
C      RECORD NUMBERS ZERO READY TO START READING
C      ANGULAR INTEGRALS
C
         READ (JBUFIR, REC=2) NS,((LAMIJ(IS,JS),IS=1,NS),JS=1,NS),
     1                        (IVCONT(ISJS),ISJS=1,(NS*(NS+1))/2 + 1)
      IF (NS .NE. NSET) THEN
         WRITE (IWRITE,3001)
         STOP
      END IF
      NRECV1=0
      NRECV2=0
C
C      LOOP OVER THE SETS
C
      DO 1  IS=1,NSET
C
      DO 2  JS=IS,NSET
C
          LAMLU=LAMIJ(IS,JS)
C
C      TEST WHETHER SETS IS,JS CAN COUPLE
C      IF THEY CAN LAMIJ CONTAINS THE LIMITS ON LAMBDA IMPOSED
C      BY SHELL COUPLING BETWEEN CONFIGURATIONS IN THE SETS
C
          IF(LAMLU .EQ. -999) GO TO 2
          IF(LAMLU .EQ. 0) GO TO 2
C
C          LAMLO=LAMLU/100
          LAMUP=MOD(LAMLU,100)
C
C      TRANSFER TARGET ANGULAR INTEGRALS FOR THIS SET
C      COMBINATION TO VSH
C
      CALL RDVSH(IS,JS,NRECV1,NRECV2,NVSHEL)
      ICOUNT=1
C
C      LOOP OVER CONFIGURATIONS IN SET IS
C
      DO 3   ICS=1,NCSET(IS)
              IF(IS .EQ. JS) THEN
                 JCSLO=ICS
              ELSE
                 JCSLO=1
              END IF
C
C      LOOP OVER CONFIGURATIONS IN SET JS
C
      DO 4   JCS=JCSLO,NCSET(JS)
C
C      THE CASE OF CONFIGURATIONS HAVING THE SAME SHELL STRUCTURE
C      IS TREATED SEPARATELY
C
          IF (IRHSGL(ICOUNT) .LT. 0) THEN
C
C      ZEROISE THE ARRAY CRLADD(LAM) READY TO ACCUMULATE TERMS
C
             DO 40  LAM=1,LAMAX
                 CRLADD(LAM)=ZERO
   40        CONTINUE
C
C      SET LAMBDA LIMITS CORRESPONDING TO VSHELL FILE
C
         LAM1=MOD (ABS (IRHSGL(ICOUNT)), 100)
         LAM2=LAMUP
C
      DO 41    LAM=LAM1,LAM2,2
           NSHSGL=ABS( IRHSGL(ICOUNT))
C
C      FIND THE NUMBER OF SHELLS CONTRIBUTING FOR THIS LAMBDA
C
              NSH=NSHSGL/10000
           IF(NSH .EQ. 0) THEN
C
C      IF NO SHELLS CONTRIBUTE HIGHER VALUES OF LAMBDA ALSO
C      CANNOT CONTRIBUTE SO JUMP OUT OF LOOP
C
              LAM2=LAM-2
           ICOUNT=ICOUNT+1
              GO TO 44
           END IF
C
           IF(LAM .EQ. 0 .OR. LAM .GT. LAMAX) THEN
C
C      SKIP TO NEXT LAMBDA LOCATION OVER THE NSH VALUES
C
           ICOUNT=ICOUNT+NSH
             GO TO 41
           END IF
      NSHHUN=NSH*100
C
C      SUM OVER SHELLS
C
      DO 42      ISH=1,NSH
               ISIG=ABS( IRHSGL(ICOUNT))/100 - NSHHUN
                 VS=VSH(ICOUNT)
             ICOUNT=ICOUNT+1
         IF(ABS(VS) .LT. EPS) GO TO 42
               LSIG=LJCOMP(ISIG)
               NSIG=NJCOMP(ISIG)
C
C      LOCATE THE RADIAL INTEGRAL
C
      CALL LOCMBB(IBBPOL, LSIG+1, LSIG+1, NSIG, NSIG, LAM, RVAL)
C
      CRLADD(LAM)=CRLADD(LAM) + VS*RVAL
      IF (NBUG2 .GT. 1) THEN
      WRITE(IWRITE,1001)IS,JS,ICS,JCS,LAM,CRLADD(LAM)
      END IF
C
   42 CONTINUE
   41 CONTINUE
C
           ELSE
C
C      NOW TREAT CASE OF CONFIGURATIONS HAVING DIFFERENT
C      SHELL STRUCTURE
C
      IRSL=IRHSGL(ICOUNT)
       IF(IRSL .EQ. 9999) THEN
         ICOUNT=ICOUNT+1
         GO TO 4
       END IF
C
C      FIND THE INTERACTING SHELLS AND THE FIRST LAMBDA FROM IRSL
C
      LAM1=MOD(IRSL,100)
        IRS=IRSL/100
       IRHO=IRS/100
       ISIG=MOD(IRS,100)
       LRHO=LJCOMP(IRHO)
       LSIG=LJCOMP(ISIG)
      LAM2=MIN(LRHO+LSIG,LAMUP)
       NRHO=NJCOMP(IRHO)
       NSIG=NJCOMP(ISIG)
C
      DO 43     LAM=LAM1,LAM2,2
C
         IF (LAM .EQ. 0 .OR. LAM .GT. LAMAX) THEN
            ICOUNT=ICOUNT+1
            GO TO 43
         END IF
C
                VS=VSH(ICOUNT)
            ICOUNT=ICOUNT+1
          IF(ABS(VS) .LE. EPS) THEN
            CRLADD(LAM)=ZERO
            GO TO 43
          END IF
C
C      LOCATE THE RADIAL INTEGRAL
C
         CALL LOCMBB(IBBPOL, LRHO+1, LSIG+1, NRHO, NSIG, LAM, RVAL)
C
      CRLADD(LAM)=VS*RVAL
      IF (NBUG2 .GT. 1) THEN
      WRITE(IWRITE,1001)IS,JS,ICS,JCS,LAM,CRLADD(LAM)
      END IF
C
   43 CONTINUE
C
      END IF
C
C      LOOP OVER ALL TARGET STATES COMPOSED FROM SET IS
C
   44 DO 5     IST=1,NSTAT(IS)
C
C      FIND THE STATE NUMBER AND THUS FIND THE COEFFICIENT AIJ
C
             ISTAT=NSETST(IS,IST)
                AI=AIJ(ISTAT,ICS)
C
C      WHEN SETS IS, JS ARE THE SAME A PAIR OF CONFIGURATIONS
C      WILL CONTRIBUTE TWICE TO COMPENSATE FOR THE SAVING MADE
C      AS A RESULT OF SYMMETRY. JCSLO=KICS IN DO LOOP 7
C
          IF(IS .EQ. JS  .AND.  ICS .NE. JCS) THEN
               AI2=AIJ(ISTAT,JCS)
          END IF
C
C      MODIFY LIMITS FOR LOOPING OVER STATES FROM SET JS WHEN IS=JS
C
           IF(IS .EQ. JS) THEN
             JSTLO=IST
           ELSE
             JSTLO=1
           END IF
C
C      LOOP OVER STATES COMPOSED FROM THE SET JS
C
      DO 6     JST=JSTLO,NSTAT(JS)
             JSTAT=NSETST(JS,JST)
                AJ=AIJ(JSTAT,JCS)
              AIAJ=AI*AJ
           IF(IS .EQ. JS  .AND.  ICS .NE. JCS) THEN
              AIAJ=AIAJ+AI2*AIJ(JSTAT,ICS)
           END IF
C
           IF(ABS(AIAJ) .LE. EPS) GO TO 6
C
C      ADD IN CONTRIBUTION TO CRL MATRIX FOR EACH LAMBDA
C
         IF (LAM1 .EQ. 0) LAM1=LAM1+2
         LAM2=MIN(LAM2,LAMAX)
C
      DO 61  LAM=LAM1,LAM2,2
         IF (JSTAT .LT. ISTAT) THEN
            CRL(JSTAT,ISTAT,LAM)=CRLADD(LAM)*AIAJ + CRL(JSTAT,ISTAT,
     1          LAM)
         ELSE
            CRL(ISTAT,JSTAT,LAM)=CRLADD(LAM)*AIAJ + CRL(ISTAT,JSTAT,
     1          LAM)
      END IF
   61 CONTINUE
C
    6 CONTINUE
C
    5 CONTINUE
C
    4 CONTINUE
C
    3 CONTINUE
C
C      TEST THAT THE VSH ARRAY HAS BEEN COMPLETELY USED FOR THIS
C      SET COMBINATION
C
      IF (ICOUNT .NE. NVSHEL + 1) THEN
         WRITE (IWRITE,3000) IS,JS
         STOP
      END IF
C
    2 CONTINUE
C
    1 CONTINUE
C
      DO 7   LAM=1,LAMAX
         IF (NBUG2 .GT. 0) THEN
         WRITE (IWRITE,1002) LAM
         WRITE (IWRITE,1003) ((CRL(ISTAT,JSTAT,LAM),ISTAT=1,NAST),
     1                        JSTAT=1,NAST)
      END IF
    7 CONTINUE
      RETURN
      END
**********************************************************************
      SUBROUTINE CUSETL
C
C      TO DETERMINE WHICH SETS COUPLE WITH EACH VALUE OF THE
C      CONTINUUM ANGULAR MOMENTUM WITHIN THE RANGE ALLOWED
C      BY LRGL.  THE INFORMATION IS STORED IN /SETL/.
C      THE COMMON BLOCK /SETL/ WILL BE STORED ON DISC FOR
C      EACH LRGL,NPTY COMBINATION.  IT IS INDEPENDENT OF THE
C      TOTAL SPIN IN THE DIRECT CASE.
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ILDIM1=MZLR3)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
C
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/SETL/LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2),
     1            NSETL(ISDIM1,ILDIM2)
      COMMON/SETS/NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1            NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2)
C
C      ARRAYS INTERMEDIATE TO /SETL/ ARRAYS
C
      DIMENSION ILVAL(ILDIM2),INSCOL(ILDIM2),INSETL(ISDIM1,ILDIM2)
C
 1000 FORMAT(////' TOTAL ANGULAR MOMENTUM ',I3,' PARITY ',I3/
     1       ' -------------------------------------')
 1001 FORMAT('      KL     L   SETS COUPLED TO L')
 1002 FORMAT(2X,15I6)
 1003 FORMAT(/'  CONTINUUM ANGULAR MOMENTUM ',
     1         'TAKES VALUES BETWEEN ',I3,' AND ',I3)
C
      WRITE(IWRITE,1000) LRGL,NPTY
      IF (NBUG8 .EQ. 1) THEN
      WRITE(IWRITE,1001)
      END IF
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 1    IS=1,NSET
            LCFG=LSET(IS)
            LMIN=MIN(LMIN,ABS(LCFG-LRGL))
            LMAX=MAX(LMAX, LCFG)
    1 CONTINUE
      LMAX=LMAX+LRGL
      ILTOT=LMAX-LMIN+1
C
      DO 2    KL=1,ILTOT
        ILVAL(KL)=LMIN+KL-1
    2 CONTINUE
C
C      LOOP OVER THIS RANGE OF CONTINUUM L
C
      DO 3    KL=1,ILTOT
               L=ILVAL(KL)
              NL=0
C
        DO 4    IS=1,NSET
              LCFG=LSET(IS)
               IPI=IPISET(IS)
C
C      PARITY TEST
C
         IF(MOD(IPI+L+NPTY,2) .NE. 0) GO TO 4
C
C      TRIANGLE RULE
C
         IF(LRGL .GT. LCFG+L  .OR.
     1      LRGL .LT. ABS(LCFG-L)) GO TO 4
C
                NL=NL+1
         INSETL(NL,KL)=IS
    4   CONTINUE
C
          INSCOL(KL)=NL
    3 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 5   IKL=1,ILTOT
         IF (INSCOL(IKL) .NE. 0) THEN
            KL=KL+1
            NSCOL(KL)=INSCOL(IKL)
            LVAL(KL)=ILVAL(IKL)
            DO 51  NL=1,NSCOL(KL)
               NSETL(NL,KL)=INSETL(NL,IKL)
   51       CONTINUE
         END IF
    5 CONTINUE
      LTOT=KL
      IF (LTOT .EQ. 0) GO TO 7
      LMIN=LVAL(1)
      LMAX=LVAL(LTOT)
      IF (NBUG8 .EQ. 1) THEN
      DO 6  KL=1,LTOT
         WRITE(IWRITE,1002)KL,LVAL(KL),(NSETL(NL,KL),NL=1,NSCOL(KL))
    6 CONTINUE
      END IF
      WRITE (IWRITE,1003) LMIN,LMAX
C
    7 RETURN
      END
**********************************************************************
      SUBROUTINE LOC1CC(L,LASTL,NRECR1,R1)
C
C      LOCATE THE ONE ELECTRON CONTINUUM CONTINUUM RADIAL INTEGRAL
C      L IS THE CONTINUUM ANGULAR MOMENTUM + 1
C      NRECR1 IS THE RECORD NUMBER OF THE BLOCK CURRENTLY IN STORE
C      R1 IS A SYMMETRIC MATRIX STORED AS A 1-DIMENSIONAL ARRAY
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM2=MZNCS)
      PARAMETER(ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(INDIM2=IDIM3*IDIM3)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
C
      DIMENSION R(IDIM3,IDIM3),R1(INDIM2)
C
      IST3L=IST3(L)
      NBUFF1=(IST3L-1)/LBUFF1      !IST3L -> IST3L-1
      IF (NBUFF1 .NE. NRECR1) THEN
         READ (JBUFF1, REC=NBUFF1)R1CC
      END IF
C
      IBUFF1=MOD(IST3L-1,LBUFF1)      !IST3L ->IST3L-1, )-1 -> )
C
      NRA=NVARY(L)
      DO  1   N=1,NRA
C
      DO  2   NP=1,N
          IBUFF1=IBUFF1+1
          IF (IBUFF1 .GT. LBUFF1) THEN
             NBUFF1=NBUFF1+1
          READ (JBUFF1, REC=NBUFF1)R1CC
          IBUFF1=IBUFF1-LBUFF1
          END IF
C
          RVALUE=R1CC(IBUFF1)
          R(NP,N)=RVALUE
          R(N,NP)=RVALUE
    2 CONTINUE
C
    1 CONTINUE
C
      I=0
C
      DO  3   N=1,NRA
C
      DO  4   NP=1,NRA
          I=I+1
          R1(I)=R(NP,N)
    4 CONTINUE
C
    3 CONTINUE
C
C      IF THE END OF THE BUFFER HAS BEEN REACHED AND
C      THE LAST LOOP ON L NOT EXECUTED READ THE NEXT
C      BLOCK OF INTEGRALS IN READINESS
C
      IF (IBUFF1+1 .GT. LBUFF1 .AND. L .LT. LASTL) THEN
         NBUFF1=NBUFF1+1
         READ (JBUFF1, REC=NBUFF1)R1CC
      END IF
C
      NRECR1=NBUFF1
C
      RETURN
      END
**********************************************************************
      SUBROUTINE LOCCCD(L1,L1P,N1,N1P,LAM,RM)
C
C      LOCATE DIRECT CONTINUUM CONTINUUM RK INTEGRALS
C      SUBROUTINE RDRK(L,LP) MUST HAVE BEEN CALLED
C      PREVIOUSLY TO SET UP THE POINTER ARRAY ICTCCD
C      AND THE STORAGE ARRAY RKCCD
C
C      IF L .NE. LP  RM IS OF LENGTH NRANG2 * NRANG2
C      IF L .EQ. LP  RM CONTAINS LOWER TRIANGLE
C      THE INTEGRALS ARE STORED IN THE ORDER:
C        L1 .LE. L1P
C        IF L1 .EQ. L1P  THEN  N1 .GE. N1P
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM2=MZNCS)
      PARAMETER(ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(INDIM2=IDIM3*IDIM3)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
C
      DIMENSION RM(INDIM2)
C
 3001 FORMAT('  **** STOP IN LOCCCD ****'/
     1       '  NO POINTER FOUND IN ICTCCD FOR'/
     2       '  L1=',I4,'L1P=',I4,'LAM=',I4)
      MAXNC1=MAXNC(L1)
      MAXNCP=MAXNC(L1P)
      NI=N1-MAXNC1
      NJ=N1P-MAXNCP
C
      IF (L1 .LT. L1P) THEN
         ICCD=ICTCCD(L1,L1P,LAM)
         NN1P=MAXNHF(L1P)-MAXNCP
         NIJ=(NI-1)*NN1P + NJ
      ELSE IF (L1 .GT. L1P) THEN
         ICCD=ICTCCD(L1P,L1,LAM)
         NN1P=MAXNHF(L1) - MAXNC1
         NIJ=(NJ-1)*NN1P + NI
      ELSE
         ICCD=ICTCCD(L1,L1P,LAM)
         IF (NI .GE. NJ) THEN
            NIJ=(NI*(NI-1))/2 + NJ
         ELSE
            NIJ=(NJ*(NJ-1))/2 + NI
         END IF
      END IF
C
      IF (ICCD .EQ. -999) THEN
         WRITE (IWRITE,3001) L1,L1P,LAM
      END IF
C
      NIJ=(NIJ-1)*ICCLNG + ICCD - ICCDIF
C
      DO  1   I=1,ICCLNG
          NIJ=NIJ+1
          RM(I)=RKCCD(NIJ)
    1 CONTINUE
C
      RETURN
      END
**********************************************************************
      SUBROUTINE LOCMBB(IBBPOL,L1,L1P,N1,N1P,LAM,RVAL)
C
C      TO LOCATE THE BOUND BOUND MULTIPOLE INTEGRAL
C
C      THE INTEGRALS ARE STORED WITH THE ORDER:
C        L1 .LE. L1P
C        IF L1 .EQ. L1P  THEN  N1 .GE. N1P
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM2=MZNCS)
      PARAMETER(ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM6=(IDIM4+1)/2)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
C
      DIMENSION IBBPOL(IDIM1,IDIM1,KDIM6)
C
 3001 FORMAT('  **** STOP IN LOCMBB ****'/
     1       '  NO POINTER FOUND IN IBBPOL FOR'/
     2       '  L1=',I4,'L1P=',I4,'LAM=',I4)
C
      LM=(LAM+1)/2
      MAXNC1=MAXNC(L1)
      MAXNCP=MAXNC(L1P)
      NI=N1-MAXNC1
      NJ=N1P-MAXNCP
C
      IF (L1 .LT. L1P) THEN
         IBB=IBBPOL(L1,L1P,LM)
         MAXN=MAXNHF(L1P)-MAXNCP
         NIJ=(NI-1)*MAXN + NJ
      ELSE IF (L1 .GT. L1P) THEN
         IBB=IBBPOL(L1P,L1,LM)
         MAXN=MAXNHF(L1)-MAXNC1
         NIJ=(NJ-1)*MAXN + NI
      ELSE
         IBB=IBBPOL(L1,L1P,LM)
         IF (NI .GE. NJ) THEN
            NIJ=(NI*(NI-1))/2 + NJ
         ELSE
            NIJ=(NJ*(NJ-1))/2 + NI
         END IF
      END IF
C
      IF (IBB .EQ. -999) THEN
         WRITE (IWRITE,3001) L1,L1P,LAM
      END IF
C
      NIJ=NIJ+IBB-ICCDIF
      RVAL=RKCCD(NIJ)
      RETURN
      END
C***********************************************************************
      SUBROUTINE ORDER(EN,NCHAN,NORDER)
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C --- DEFINE NORDER, WHICH POINTS TO VALUES OF EN IN ASCENDING ORDER.
C
C --- INPUT PARAMETERS ...
C     (EN(M),M=1,NCHAN) = ENERGIES IN ARBITRARY ORDER;
C --- OUTPUT PARAMETER ...
C     (NORDER(M),M=1,NCHAN) = THE POSITION OF THE M-TH BIGGEST ENERGY
C                             IN THE EN ARRAY.
C
      DIMENSION EN(NCHAN),NORDER(NCHAN)
C
      NORDER(1)=1
      IF(NCHAN.EQ.1) RETURN
      DO 9 M=2,NCHAN
      J=M
      X=EN(J)+1.D-7
      J1=J-1
      DO 7 I=1,J1
      IF(J.LT.M) GO TO 6
      I1=NORDER(I)
      IF(X.GT.EN(I1)) GO TO 7
    6 J=J-1
      NORDER(J+1)=NORDER(J)
    7 CONTINUE
    8 NORDER(J)=M
    9 CONTINUE
      RETURN
      END
**********************************************************************
      SUBROUTINE RDHMAT(ISP,MORE)
C
C      READS THE HAMILTONIAN MATRIX ELEMENTS STORED BY STHMAT
C      AND ARRANGES THEM IN HMAT READY FOR THE HOUSEHOLDER
C      DIAGONALISATION ROUTINE
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS)
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(INDIM2=IDIM3*IDIM3)
      PARAMETER(INDIM3=INDIM2+IDIM4)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
C
      PARAMETER (LIWORK=10*IHDIM1)                                  !LAPACK
      PARAMETER (LWORK=2*IHDIM1*IHDIM1+6*IHDIM1+1)                  !DSYEVD
C
      CHARACTER*4 PARITY(0:1)
C
      COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2),
     1              ICONCH(ITDIM3),KCONCH(ITDIM3),
     2              L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1),
     3              NSCHAN(IPDIM1),NSPREC(IPDIM1)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/HMAT(IHDIM4)
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/SETL  /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2),
     1              NSETL(ISDIM1,ILDIM2)
      COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
      COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1)
      COMMON/STG1  /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2)
C
      COMMON /LAP1/A(IHDIM1,IHDIM1),XA(IHDIM1,IHDIM1),AUX(IHDIM1,6),DUM !LAPACK
      DIMENSION ISUPP(2*IHDIM1),IWORK(LIWORK)                      !LAPACK
      DIMENSION WORK(LWORK)                                        !DSYEVD
      EQUIVALENCE (WORK(1),A(1,1))                                 !DSYEVD
C
      DIMENSION BUF1D(INDIM3),H1MAT(IDIM3,IDIM3)
      DIMENSION VALUE(IHDIM1), X(IHDIM1)
C     DIMENSION WMAT(ITDIM3)
      DIMENSION WMAT(ITDIM6,IHDIM1)
C     DIMENSION CF(ITDIM3,IDIM4)
      DIMENSION CF(ITDIM6,ITDIM6,IDIM4)
      DIMENSION ISTCHS(ITDIM6)
C
      DATA EPSI/1.0D-9/
      DATA PARITY/'EVEN','ODD'/
C
 1000 FORMAT(/,28X,'SUBROUTINE RDHMAT'/28X,'-----------------'/)
 1001 FORMAT(' ASYMPTOTIC COEFFS FOR LAMBDA = 1 TO LAMAX FOR ',
     1       'TARGET SPIN (2S+1) =',I3/)
 1002 FORMAT(' CHANNEL',I3,' (STATE',I3,', L=',I3,
     1           '), WITH CHANNEL',I3,' (STATE',I3,', L=',I3,')')
 1003 FORMAT ('  SURFACE AMPLITUDE')
 1004 FORMAT(' HAMILTONIAN MATRIX FOR TARGET SPIN (2S+1) =',I3/)
 1005 FORMAT( 2X,8F14.7)
 1006 FORMAT(16X,7F14.7)
 1007 FORMAT(30X,6F14.7)
 1008 FORMAT(44X,5F14.7)
 1009 FORMAT(58X,4F14.7)
 1010 FORMAT(72X,3F14.7)
 1011 FORMAT(86X,2F14.7)
 1012 FORMAT(100X,F14.7)
 1013 FORMAT(' CONTRIBUTION FROM CHANNEL',I3,' (STATE',I3,', L=',I3,
     1           '), WITH CHANNEL',I3,' (STATE',I3,', L=',I3,')')
 1014 FORMAT (/' EIGENVALUE=',F12.7,' RYD - RELATIVE TO TARGET GROUND',
     1        F12.7,' AU')
 1015 FORMAT (' VECTOR')
C
 2001 FORMAT (5F14.7)
C
      IF (NBUG4 .GT. 0  .OR.  NBUG9 .GT. 0) THEN
         WRITE(IWRITE,1000)
      END IF
      IF (NBUG4 .GT. 0) THEN
C
C      CONSTRUCT THE ARRAY ISTCHS(ICH), THE STATE NUMBER COUPLED TO
C      CHANNEL ICH
C
         ICH=0
         DO 10 IST=1,NAST
            DO 11 K=1,KCONAT(IST,ISP)
               ICH=ICH+1
               ISTCHS(ICH)=IST
   11       CONTINUE
   10    CONTINUE
      END IF
C
      IF (NBUG4 .GT. 1) THEN
         WRITE(IWRITE,1004) NSPVAL(ISP)
      END IF
C
C      THE RECORD NUMBERS WERE CALCULATED IN STHMAT TO BE
C      IN THE CORRECT ORDER TO FILL THE HMAT ARRAY.
C      ONLY THE UPPER HALF OF THE SUBMATRICES ON THE
C      DIAGONAL OF THE HAMILTONIAN MUST BE STORED,
C      ALL OTHERS REMAIN RECTANGULAR.
C
      NCH=NSCHAN(ISP)
      IF (NCH .EQ. 0  .AND.  MORE .NE. 0) RETURN
      IREC=NSPREC(ISP)
      IREC2=NSPREC(ISP+1)
      IF (NCH .NE. 0) THEN
         READ (JBUFD, REC=IREC) BUF1D
      END IF
      NY=LENHAM(ISP)
      NY2= (NY*(NY+1))/2
C
C      WRITE ASYMPTOTIC FILE
C
      WRITE(ITAPE3) LRGL,-NSPVAL(ISP),NPTY,NCH,NY,MORE
      WRITE(ITAPE3)(KCONAT(I,ISP),I=1,NAST)
      WRITE(ITAPE3)(L2PSPN(ICH,ISP),ICH=1,NCH)
      IF (NCH .EQ. 0) THEN
         WRITE(ITAPE3)ZERO
         WRITE(ITAPE3)ZERO
         WRITE(ITAPE3)ZERO
         RETURN
      END IF
C
C      ZEROISE HAMILTONIAN MATRIX
C
      DO 1     I=1,NY2
         HMAT(I)=ZERO
    1 CONTINUE
C
C      PREPARE TO FILL CF BUFFER WITH ASYMPTOTIC COEFFS
C      FROM THE START OF EACH HAMILTONIAN BLOCK.
C      THERE IS A COEFFICIENT FOR EACH VALUE OF LAMBDA FROM
C      ONE TO LAMAX
C
C      IBUF=0
C
C      PREPARE TO FILL HMAT
C      SET INITIAL ROW LENGTH AND STARTING POSITION FOR
C      FIRST ROW OF SUB MATRICES
C
             IROW = NY
            ISTART= 1
      DO 2    ICH=1,NCH
         L=L2PSPN(ICH,ISP)+1
         NRA1=NVARY(L)
C
C      SET STARTING POSITION FOR FIRST ROW OF SUB MATRICES
C
           JSTART=ISTART
C
      DO 3    JCH=ICH,NCH
         LP=L2PSPN(JCH,ISP)+1
         NRA2=NVARY(LP)
      DO 31   IB=1,LAMAX
C        CF(JCH,IB)=BUF1D(IB)
         CF(ICH,JCH,IB)=BUF1D(IB)
         CF(JCH,ICH,IB)=BUF1D(IB)
   31 CONTINUE
C
C      COPY THE SUB MATRIX FROM THE BUFFER TO THE
C      ARRAY H1MAT, DIMENSION NRA1 X NRA2
C
      IB=LAMAX
      DO 4      N=1,NRA1
C
      DO 5     NP=1,NRA2
           IB=IB+1
           H1MAT(NP,N)=BUF1D(IB)
    5 CONTINUE
C
    4 CONTINUE
C
           IREC = IREC+1
        IF(IREC .LT. IREC2) THEN
C
C      READ THE NEXT SUB MATRIX INTO THE BUFFER IN READINESS
C
          READ (JBUFD, REC=IREC) BUF1D
C
         END IF
C
        IF(ICH .EQ. JCH) THEN
C
C      ON THE DIAGONAL STORE UPPER HALF ONLY
C
          NSTART = JSTART
            NROW = IROW
C
         DO 61      N=1,NRA1
            NPSTAR = NSTART
C
         DO 71     NP=N,NRA2
            HMAT(NPSTAR) = H1MAT(NP,N)
            NPSTAR = NPSTAR + 1
   71    CONTINUE
C
            NSTART = NSTART + NROW
              NROW = NROW - 1
   61    CONTINUE
         JSTART=ISTART+NRA2
C
         ELSE
C
C      OFF THE DIAGONAL STORE RECTANGULAR ARRAY
C
              NROW = IROW - 1
            NSTART = JSTART
C
         DO 62      N=1,NRA1
            NPSTAR = NSTART
C
         DO 72     NP=1,NRA2
            HMAT(NPSTAR) = H1MAT(NP,N)
            NPSTAR = NPSTAR + 1
   72    CONTINUE
C
            NSTART = NSTART + NROW
              NROW = NROW - 1
   62    CONTINUE
         JSTART=JSTART+NRA2
C
        END IF
C
      IF (NBUG4 .GT. 1) THEN
         WRITE(IWRITE,1013) ICH,ISTCHS(ICH),L2PSPN(ICH,ISP),
     1                      JCH,ISTCHS(JCH),L2PSPN(JCH,ISP)
         IF (NBUG4 .EQ. 2) THEN
            NPR1=MIN(5,NRA1)
            NPR2=MIN(5,NRA2)
         ELSE
            NPR1=NRA1
            NPR2=NRA2
         END IF
C
         IF (ICH .EQ. JCH) THEN
C
C      PRINT TRIANGULAR MATRIX
C
            NBLOCK=NPR1/8+1
            NLO=1
            NUP=MIN(8,NPR1)
            DO 12 NBLOC=1,NBLOCK
               IBLANK=0
               DO 13 N=NLO,NUP
                  NPLO=N
                  NPUP=NUP
                  DO 14 NBLOCP=NBLOC,NBLOCK
                     IF (NBLOCP .EQ. NBLOC .AND. IBLANK .GT. 0) THEN
                        IF (IBLANK .EQ. 1)THEN
                        WRITE(IWRITE,1006)(H1MAT(NP,N),NP=NPLO,NPUP)
                        ELSE IF (IBLANK .EQ. 2) THEN
                        WRITE(IWRITE,1007)(H1MAT(NP,N),NP=NPLO,NPUP)
                        ELSE IF (IBLANK .EQ. 3) THEN
                        WRITE(IWRITE,1008)(H1MAT(NP,N),NP=NPLO,NPUP)
                        ELSE IF (IBLANK .EQ. 4) THEN
                        WRITE(IWRITE,1009)(H1MAT(NP,N),NP=NPLO,NPUP)
                        ELSE IF (IBLANK .EQ. 5) THEN
                        WRITE(IWRITE,1010)(H1MAT(NP,N),NP=NPLO,NPUP)
                        ELSE IF (IBLANK .EQ. 6) THEN
                        WRITE(IWRITE,1011)(H1MAT(NP,N),NP=NPLO,NPUP)
                        ELSE
                        WRITE(IWRITE,1012)(H1MAT(NP,N),NP=NPLO,NPUP)
                        END IF
                     ELSE
                        WRITE(IWRITE,1005)(H1MAT(NP,N),NP=NPLO,NPUP)
                     END IF
                     NPLO=NPUP+1
                     NPUP=MIN(NPUP+8,NPR2)
   14             CONTINUE
                  IBLANK=IBLANK+1
   13          CONTINUE
               NLO=NUP+1
               NUP=MIN(NUP+8,NPR1)
   12       CONTINUE
         ELSE
C
C      PRINT RECTANGULAR MATRIX
C
            NBLOCK=NPR2/8+1
            DO 15 N=1,NPR1
               NPLO=1
               NPUP=MIN(8,NPR2)
               DO 16 NBLOCP=1,NBLOCK
                  WRITE(IWRITE,1005)(H1MAT(NP,N),NP=NPLO,NPUP)
                  NPLO=NPLO+8
                  NPUP=MIN(NPUP+8,NPR2)
   16          CONTINUE
   15       CONTINUE
         END IF
      END IF
C
    3 CONTINUE
C
C      WRITE CF MATRIX TO ASYMPTOTIC FILE
C
C     WRITE (ITAPE3) ((CF(JCH,ILAM),ILAM=1,LAMAX),JCH=ICH,NCH)
C
            ISTART = NPSTAR
              IROW = IROW - NRA1
C
    2 CONTINUE
C
C      WRITE CF MATRIX TO ASYMPTOTIC FILE
C
      WRITE (ITAPE3) (((CF(ICH,JCH,ILAM),ICH=1,NCH),JCH=1,NCH),
     1                  ILAM=1,LAMAX)
C
      IF (NBUG4 .GT. 0) THEN
         WRITE(IWRITE,1001) NSPVAL(ISP)
         DO 17 ICH=1,NCH
            DO 18 JCH=ICH,NCH
               WRITE(IWRITE,1002)ICH,ISTCHS(ICH),L2PSPN(ICH,ISP),
     1                           JCH,ISTCHS(JCH),L2PSPN(JCH,ISP)
               WRITE(IWRITE,1005)(CF(ICH,JCH,ILAM),ILAM=1,LAMAX)
   18       CONTINUE
   17    CONTINUE
      END IF
C
C      DIAGONALIZE THE HAMILTONIAN MATRIX AND
C      CALCULATE THE SURFACE AMPLITUDES
C      WRITE THEM TO ITAPE3 FOR THE ASYMPTOTIC STAGE
C
      GROUND=ENAT(1)
C
      IF (NBUG9 .GT. 0) THEN
         WRITE(IWRITE,'(/'' EIGENVALUES, EIGENVECTORS AND SURFACE''
     1              ,'' AMPLITUDES FOR TARGET SPIN (2S+1) ='',I3)')
     2                NSPVAL(ISP)
      END IF
      DO 8 NO=1,NY
C
C      THIS DIAGONALIZATION ROUTINE RETURNS EIGENVECTORS ONE-BY-ONE.
C
CORIG         CALL HSLDR(NY,HMAT,NY2,EPSI,VALUE,X,NO)          !ORIG
C
C      LAPACK 3.0 DRIVERS
C
        IF(NO.EQ.1)THEN                                       !LAPACK START
          K=0
          DO I=1,NY
            DO J=I,NY
              K=K+1
C              A(J,I)=HMAT(K)                                  !DSYEVR
              XA(J,I)=HMAT(K)                                 !DSYEVD
            ENDDO
          ENDDO
C
C          IEEEOK=ILAENV(10,'DSYEVR','N',1,2,3,4)
C          IF(IEEEOK.NE.1)THEN
C            WRITE(6,*)'WARNING: DSYEVR USING DSTEB/DSTEIN, CALL DSYEVD!'
C            STOP 'RDHMAT: BEST TO CALL DSYEVD!'
C          ENDIF
C
C DSYEVR (REQUIRES MACHINES TO HANDLE NaNs GRACEFULLY)
C
C          CALL DSYEVR('V','A','L',NY,A,IHDIM1,VL,VU,IL,IU,-EPSI,ME,
C     X                X,XA,IHDIM1,ISUPP,HMAT,IHDIM4,IWORK,LIWORK,INFO)
C          IF(ME.NE.NY)THEN
C            WRITE(6,*)'RDHMAT: ERROR IN LAPACK DSYEVR, NOT ALL',
C     X      ' E-VALUES FOUND:',ME,NY
C            STOP 'RDHMAT: FAILURE IN LAPACK ROUTINE DSYEVR'
C          ENDIF
C          IF(HMAT(1).GT.IHDIM4)THEN
C            LWORK=NINT(HMAT(1))
C            WRITE(6,*)'***OPTIMAL USE OF DSYEVR REQUIRES LWORK=',LWORK
C          ENDIF
C
C DSYEVD (TO BE USED IF NaNs CAUSE CALCULATION TO ABORT)
C
          CALL DSYEVD('V','L',NY,XA,IHDIM1,X,WORK,LWORK,IWORK,LIWORK
     X               ,INFO)
C
          IF(INFO.NE.0)THEN
            WRITE(6,*)'RDHMAT: ERROR IN LAPACK DSYEVR/D: INFO=',INFO
            STOP 'RDHMAT: FAILURE IN LAPACK ROUTINE DSYEVR/D'
          ENDIF
C
          DO I=1,NY
            VALUE(I)=X(NY+1-I)
          ENDDO
        ENDIF
C
        J=NY+1-NO
        DO I=1,NY
          X(I)=XA(I,J)
        ENDDO                                                    !LAPACK END
C
C       DEBUG PRINTOUT OF E-SOLUTIONS
C
         IF (NBUG9 .GT. 0) THEN
         ENLEVL=(VALUE(NO)-GROUND)*TWO
         WRITE (IWRITE,1014)VALUE(NO),ENLEVL
         WRITE (IWRITE,1015)
         WRITE (IWRITE,1005)(X(I),I=1,NY)
         END IF
         NCH=NSCHAN(ISP)
         I1=0
      DO 81   ICH=1,NCH
         LP=L2PSPN(ICH,ISP)+1
         NOTERM=NVARY(LP)
         SUM=ZERO
      DO 82   N=1,NOTERM
         I1=I1+1
         SUM=SUM+X(I1)*ENDS(N,LP)
   82 CONTINUE
C        WMAT(ICH)=SUM
         WMAT(ICH,NO)=SUM
   81 CONTINUE
C        WRITE (IWRITE,1005) (WMAT(ICH), ICH=1,NCH)
C        WRITE (ITAPE3) VALUE(NO),(WMAT(ICH),ICH=1,NCH)
    8 CONTINUE
         IF (NBUG9 .GT. 0) THEN
            WRITE (IWRITE,1003)
            WRITE (IWRITE,1005) ((WMAT(ICH,NO), ICH=1,NCH),NO=1,NY)
         END IF
C
      WRITE (ITAPE3) (VALUE(I),I=1,NY)
      WRITE (ITAPE3) ((WMAT(K,I),K=1,NCH),I=1,NY)
C
      WRITE(IWRITE,'(//I7,'':  '',''HAMILTONIAN DIAGONALISED FOR L='',
     1        I3,'', PARITY='',A4,'', TARGET SPIN (2S+1)='',I2//)')
     2     NY,LRGL,PARITY(NPTY),NSPVAL(ISP)
C
      RETURN
C
      END
**********************************************************************
      SUBROUTINE RDRINT(ISTART,IFIN,NB,NRECR)
C
C      A DOUBLE BUFFERING READ ROUTINE
C      READS RADIAL INTEGRALS FROM FILE JBUFF1 INTO ARRAY
C      RKCCD STARTING AT FILE POSITION ISTART TO FILE POSITION
C      IFIN.  CALLED FROM CRLMAT TO READ BOUND BOUND MULTIPOLE
C      INTEGRALS AND FROM RDRKCC TO READ RK INTEGRALS FOR GIVEN L,LP
C
C      NRECR IS THE RECORD NUMBER OF THE BLOCK CURRENTLY IN STORE
C      NB IS THE BUFFER INDICATOR IN THE DOUBLE BUFFERING
C      NRECR AND NB ARE SET ZERO ON FIRST ENTRY
C      THEY ARE NOT SET THEREAFTER BUT CONTAIN THE VALUES
C      LEFT BY THE PREVIOUS EXECUTION
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM2=MZNCS)
      PARAMETER(ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/RADBUF/BUFF1(LBUFF1),BUFF2(LBUFF1)
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
C
 3000 FORMAT('  **** STOP IN RDRINT ****'/
     1       '  DIMENSION OF RKCCD IS NOT LARGE ENOUGH'/
     2       '  L=',I3,'LP=',I3,'IRDIM1=',I7,'ISTART=',I7,'IFIN=',I7/
     3       'TRY INCREASING MZMEG TO:',I3)
C
C
C      CHECK WHETHER DIMENSION OF RKCCD IS LARGE ENOUGH
C
      IF ( IFIN-ISTART+1 .GT. IRDIM2 ) THEN
        MXMEM=IRDIM2/1000000+1
         WRITE(IWRITE,3000) IRDIM2,ISTART,IFIN,MXMEM
         STOP
      END IF
C
C      STORE ICCDIF IN COMMON READY TO BE SUBTRACTED FROM
C      STARTING VALUES IN POINTER ARRAY
C
      ICCDIF = ISTART
C
      NBLK1 = ISTART/LBUFF1
      NBLK2 = IFIN/LBUFF1
      IB1 = MOD(ISTART,LBUFF1)
      IB2 = MOD(IFIN,LBUFF1)
      IF (IB1 .EQ. 0) THEN
         NBLK1=NBLK1-1
         IB1=LBUFF1
      END IF
      IF (IB2 .EQ. 0) THEN  ! CHECK IB2 AS WELL - NRB
         NBLK2=NBLK2-1
         IB2=LBUFF1
      END IF
C
C      READ IN THE BLOCKS OF INTEGRALS USING TWO BUFFERS
C      BUT CHECK WHETHER FIRST BLOCK IS ALREADY IN THE BUFFER
C
      IF (NRECR .NE. NBLK1) THEN
         READ (JBUFF1, REC=NBLK1) BUFF1
         NB=0
      END IF
      I=0
      IBLO=IB1
C
      DO 1     NBLK=NBLK1+1, NBLK2
         IF(NB .EQ. 0) THEN
           READ( JBUFF1, REC=NBLK) BUFF2
           DO 21     IB=IBLO,LBUFF1
              I=I+1
              RKCCD(I)=BUFF1(IB)
   21      CONTINUE
           NB=1
         ELSE
           READ( JBUFF1, REC=NBLK) BUFF1
           DO 22     IB=IBLO,LBUFF1
              I=I+1
              RKCCD(I)=BUFF2(IB)
   22      CONTINUE
           NB=0
         END IF
      IBLO=1
    1 CONTINUE
C
C      COPY THE LAST BUFFER INTO RKCCD
C
      IF (NB .EQ. 0) THEN
         DO 31   IB=IBLO,IB2
            I=I+1
            RKCCD(I)=BUFF1(IB)
   31    CONTINUE
      ELSE
         DO 32   IB=IBLO,IB2
            I=I+1
            RKCCD(I)=BUFF2(IB)
   32    CONTINUE
      END IF
      NRECR=NBLK2
C
      RETURN
      END
**********************************************************************
      SUBROUTINE RDRKCC(L,LP,NRECIC,NB,NRECRK)
C
C      FOR THIS L,LP COMBINATION STORE ALL CONTINUUM/CONTINUUM
C      RK INTEGRALS IN RKCCD
C
C      NRECIC IS THE RECORD NUMBER OF THE BLOCK CURRENTLY IN STORE
C      NRECRK IS THE RECORD NUMBER OF THE RK BLOCK CURRENTLY IN STORE
C      NB IS THE BUFFER INDICATOR IN THE DOUBLE BUFFERING OF THE RK
C      NRECIC, NRECRK AND NB ARE SET ZERO ON FIRST ENTRY
C      THEY ARE NOT SET THEREAFTER BUT CONTAIN THE VALUES
C      LEFT BY THE PREVIOUS EXECUTION
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
C
      PARAMETER(LBUFF1=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
 1001 FORMAT(' **** ERROR IN RDRKCC FOR L,LP=',2I4/
     1       ' CANNOT FIND STARTING RECORD FOR L1,L1P,LAMMIN=',3I4/
     2       ' TOTAL NUMBER OF RK INTEGRALS STORED FOR L,LP=',I4)
C
C      IF TRIANGULAR RULES FOR LAMMIN,L,LP AND LAMMIN,L1,L1P ARE
C      NOT SATISFIED RETURN
C      STORE POINTER ARRAY
C      CALCULATE START AND END POSITIONS IN THE STORAGE OF THE
C      RADIAL INTEGRALS FOR THIS L, LP
C
      IF (L .EQ. LP) THEN
         NRA=NVARY(L+1)
         ICCLNG = (NRA*(NRA+1))/2
      ELSE
         NRA1=NVARY(L+1)
         NRA2=NVARY(LP+1)
         ICCLNG = NRA1 * NRA2
      END IF
C
      LDIF=LP-L+1
      LAMMIN=LDIF
      IF (LAMMIN .GT. 2*LRANG1-1) RETURN
C
      NRECIC=L*LLPDIM + LDIF
      READ (JBUFF2, REC=NRECIC)(((ICTCCD(I,J,K),I=1,LRANG1),
     1                         J=1,LRANG1),K=1,LRANG1+LRANG1-1)
C
C      FIND START AND END POSITIONS OF RK INTEGRALS
C
      L1=1
      L1P=LAMMIN
      IF (L1P .GT. LRANG1) THEN
         L1=LAMMIN-LRANG1+1
         L1P=LRANG1
      END IF
      ICCD1=ICTCCD(L1,L1P,LAMMIN)
      NRKVAL=NCCD(LDIF,L+1)
      ICCD2=ICCD1+NRKVAL-1
C
      IF (ICCD1 .LT. LBUFF1) THEN
         WRITE(IWRITE,1001) L,LP,L1,L1P,LAMMIN,NRKVAL
        CALL EXIT (0)
      END IF
C
      CALL RDRINT(ICCD1,ICCD2,NB,NRECRK)
C
      RETURN
      END
**********************************************************************
      SUBROUTINE RDVSH(KIS,KJS,NRECV1,NRECV2,NVSHEL)
C
C      TO READ IN VSH FILE FOR SETS KIS,KJS
C      NRECV1 IS THE RECORD IN BUFV1
C      NRECV2 IS THE RECORD IN BUFV2
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM3=MZNR2)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(ILDIM4=ILDIM3+ILDIM3)
      PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/ANGBUF/BUFV1(LBUFFV),BUFV2(LBUFFV),
     1              IBUFV1(LBUFFV),IBUFV2(LBUFFV)
      COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4),
     1              IRHSGL(IVDIM1)
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/SETS  /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1              NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2)
C
 1001 FORMAT(2I5,5X,I8,5X,I8,5X,E13.5)
C
 3000 FORMAT ('  **** STOP IN RDVSH **** '/
     1        '  DIMENSION OF VSH NOT LARGE ENOUGH'/
     2        '  KIS=',I4,' KJS=',I4,' IVDIM1=',I7,' ICONT1=',I8,
     3        ' ICONT2=',I8)
C
      KISJS=NSET*(KIS-1) - (KIS*(KIS-1))/2 +KJS
      ICONT1=IVCONT(KISJS)
      ICONT2=IVCONT(KISJS+1) -1
      NVSHEL=ICONT2-ICONT1+1
      IF (ICONT2-ICONT1 .GE. IVDIM1) THEN
         WRITE(IWRITE, 3000) KIS,KJS,IVDIM1,ICONT1,ICONT2
         STOP
      END IF
C
      NBLK1=ICONT1/LBUFFV
      NBLK2=ICONT2/LBUFFV
      IB1=MOD(ICONT1,LBUFFV)
      IB2=MOD(ICONT2,LBUFFV)
      IF (IB1 .EQ. 0) THEN
         NBLK1=NBLK1-1
         IB1=LBUFFV
      END IF
C
      IF (NRECV1 .NE. NBLK1) THEN
         IF (NRECV2 .NE. NBLK1) THEN
            READ (JBUFFV, REC=NBLK1) BUFV1
            READ (JBUFIR, REC=NBLK1) IBUFV1
            NB=0
            NRECV1=NBLK1
         ELSE
            NB=1
         END IF
      ELSE
         NB=0
      END IF
C
      I=0
      IBLO=IB1
      DO  1     NBLK=NBLK1+1,NBLK2
          IF (NB .EQ. 0) THEN
             IF (NRECV2 .NE. NBLK) THEN
                READ (JBUFFV, REC=NBLK) BUFV2
                READ (JBUFIR, REC=NBLK) IBUFV2
                NRECV2=NBLK
             END IF
             NB=1
             DO  21   IB=IBLO,LBUFFV
                 I=I+1
                 VSH(I)=BUFV1(IB)
                 IRHSGL(I)=IBUFV1(IB)
   21        CONTINUE
           ELSE
             IF (NRECV1 .NE. NBLK) THEN
                READ (JBUFFV, REC=NBLK) BUFV1
                READ (JBUFIR, REC=NBLK) IBUFV1
                NRECV1=NBLK
             END IF
             NB=0
             DO 22    IB=IBLO,LBUFFV
                I=I+1
                VSH(I)=BUFV2(IB)
                IRHSGL(I)=IBUFV2(IB)
   22        CONTINUE
          END IF
          IBLO=1
    1  CONTINUE
C
C       COPY THE LAST BUFFER
C
       IF (NB .EQ. 0) THEN
          DO 31   IB=IBLO,IB2
             I=I+1
             VSH(I)=BUFV1(IB)
             IRHSGL(I)=IBUFV1(IB)
   31     CONTINUE
          NB=1
       ELSE
          DO 32   IB=IBLO,IB2
             I=I+1
             VSH(I)=BUFV2(IB)
             IRHSGL(I)=IBUFV2(IB)
   32     CONTINUE
          NB=0
        END IF
      IF (NBUG5 .EQ. 1) THEN
        IR0=0                  !CAN BE RESET TO MATCH PARTICULAR STG1NX
        WRITE(IWRITE,1001)
     X    (KIS,KJS,MOD(IR0+IR-1,LBUFFV)+1,IRHSGL(IR),VSH(IR),IR=1,I)
      ENDIF
C
      RETURN
      END
**********************************************************************
      SUBROUTINE RECOV1(I,IDAT,ICURR)
*
*      THIS ROUTINE IS CALLED ONLY IN THE CASE OF ARRAY OVERFLOW.
*      IF IPLACE=0 THE PROGRAM STOPS HERE, OTHERWISE THE PROGRAM RETURNS
*      TO THE CALLING ROUTINE AFTER SETTING IPLACE=I.
*
*      I IS THE POSITION IN THE IDRAY AND IPRAY LISTS HOLDING
*      THE NAME OF THE RELEVANT DIMENSION PARAMETER AND
*      PREPROCESSOR PARAMETER.
*      IDAT IS THE ARRAY SIZE REQUIRED BY THE DATA
*      ICURR IS THE CURRENT ARRAY SIZE
*
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
      CHARACTER*6 IDRAY
      CHARACTER*3 IPRAY
C
C      PARAMETER(IDIM1=MZLR1,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     3  LBUFF1=MZBUF)
C      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
C
      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)
*
 1001 FORMAT(/' PROGRAM TERMINATES IN RECOV1'/)
      WRITE(IWRITE,1000)IDRAY(I),IPRAY(I),ICURR,IDAT
    1 WRITE(IWRITE,1001)
      CALL EXIT (0)
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE RECSET(LLLO,LLUP,ICONTN)
C
C      READS STGANG DATA AND RECONSTRUCTS COMMON/SETS/
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
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      CHARACTER*8 FILEA,FILEB,FILEC,FILED
      CHARACTER*4 TITLE,CONT
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1)
      PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(ICDIM3=ICDIM2+ICDIM2-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
C
      PARAMETER(LBUFFV=MZBUF,LBUFFZ=MZBUF)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/SETS  /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1              NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2)
      COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2)
      COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
C
      COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN
      COMMON /NRB/MAXC
C
      NAMELIST/STGNX/MINLT,MAXLT,LRGLLO,LRGLUP,MAXC,LFIXN,MJS,NPTYIN
     X        ,NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9
C
      DIMENSION NOCCSH(ICDIM1),NOCORB(ICDIM2,ICDIM1),
     1          NELCSH(ICDIM2,ICDIM1),J1QNRD(ICDIM3,3,ICDIM1)
      DIMENSION TITLE(18)
      DIMENSION LVAL(ILDIM2),ILVAL(ILDIM2),NSCOL(ILDIM2)
 1000 FORMAT('  ',10X,'STGANG DATA'//)
 1001 FORMAT('    LSET  ISPSET  IPISET   NCSET  NCFGSE')
 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'/)
 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
      DATA CONT/'CONT'/
C
C NRB
      NPTYMN=0
      NPTYMX=1
      NPTYIN=-1
C
C
      IREAD=5
      IWRITE=6
      OPEN(UNIT=IREAD,FILE='dstgnx',STATUS='UNKNOWN')
      OPEN(UNIT=IWRITE,FILE='rout3nx',STATUS='UNKNOWN')
      JBUFIR=20
      JBUFFV=21
      JBUFFZ=22
      FILEA='ANG1.DAT'
      FILEB='ANG2.DAT'
      FILEC='ANG3.DAT'
      REWIND (IREAD)
      READ(IREAD,2001) (TITLE(I),I=1,18)
      IF (TITLE(1) .EQ. CONT) THEN
         INX2=36
         OPEN (INX2,FILE='NX2.DAT',ACCESS='SEQUENTIAL',
     1         STATUS='OLD',FORM='UNFORMATTED')
C        OPEN (INX2,FILE='tapenx2',ACCESS='SEQUENTIAL',
C    1         STATUS='OLD',FORM='UNFORMATTED')
         ICONTN=1
         READ (INX2)MAXN,MAXORB,LRANG3,NSETS1,MAXNCF,NAST,MAXNST,
     1              NCHSUM,MAXNCH,ISRAN3,NCFG
         IF (ICDIM1 .LT. NCFG) CALL RECOV1(24,NCFG,ICDIM1)
      ELSE
         ICONTN=0
         READ(IREAD, *)  IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,
     1                  IBUG8,IBUG9
      END IF
C
C      OPEN ANGULAR FILES
C
      OPEN(UNIT=JBUFIR,FILE=FILEA,ACCESS='DIRECT',STATUS='OLD',
     1     FORM='UNFORMATTED',RECL=MACDIM*LBUFFV)
      OPEN(UNIT=JBUFFV,FILE=FILEB,ACCESS='DIRECT',STATUS='OLD',
     1     FORM='UNFORMATTED',RECL=MACDIM*LBUFFV)
      OPEN(UNIT=JBUFFZ,FILE=FILEC,ACCESS='DIRECT',STATUS='OLD',
     1     FORM='UNFORMATTED',RECL=MACDIM*LBUFFZ)
C
CW    WRITE(IWRITE,1000)
C
      MAXC=1
      IF(ICONTN.EQ.1)THEN
      MAXLT=0
      MINLT=0
      LRGLLO=0
      LRGLUP=0
      NBUG1=0
      NBUG2=0
      NBUG3=0
      NBUG4=0
      NBUG5=0
      NBUG6=0
      NBUG7=0
      NBUG8=0
      NBUG9=0
C
C
      READ(IREAD,STGNX)
C
C
      IF(NPTYIN.EQ.0)NPTYMX=0
      IF(NPTYIN.EQ.1)NPTYMN=1
      IF(MAXLT.GT.0)LRGLUP=MAXLT
      IF(MINLT.GT.0)LRGLLO=MINLT
      ELSE
      READ (IREAD,*) LRGLLO,LRGLUP
      ENDIF
      LLLO=LRGLLO
      LLUP=LRGLUP
CW    WRITE (IWRITE,1005) LRGLLO,LRGLUP
      IF (ILDIM3 .LE. LRGLUP) CALL RECOV1(16,LRGLUP+1,ILDIM3)
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)
C
      ELSE
C
         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
                  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,*)
   16       CONTINUE
            IREAD1=36
            FILED='CONF.DAT'
            OPEN(UNIT=36,FILE=FILED,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
      NCFGSE(IS,ICS)=IC
         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
      NCFGSE(IS,ICS)=IC
C
    2 CONTINUE
C
      NSET=IS
C
      DO 4   I=1,NSET
            NC=NCSET(I)
CW      WRITE(IWRITE,1002)LSET(I),ISPSET(I),IPISET(I),NCSET(I),
CW   1               (NCFGSE(I,J),J=1,NC)
    4 CONTINUE
      LRANG3=LRANG3+1
      IF (ILDIM1 .LT. LRANG3) CALL RECOV1(15,LRANG3,ILDIM1)
      LLPDIM=LRANG3+LRANG3-1
C
C      FIND LRANG2 NEEDED FOR REQUIRED RANGE OF LRGL
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
      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)
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
       NSCOL(KL)=1
       GO TO 9
   10   CONTINUE
C
    9 CONTINUE
      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
      IF (LTOT .EQ. 0) GO TO 6
      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
      IF (IDIM2 .LT. LRANG2) CALL RECOV1(2,LRANG2,IDIM2)
C
      RETURN
      END
**********************************************************************
      SUBROUTINE RETAP1
C
C      READS BASIC INFORMATION FROM THE FILE RAD3
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/STG1  /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2)
C
 1000 FORMAT('  ',16I4)
 1001 FORMAT('  ',6E14.7)
      REWIND ITAPE1
      ITAPE=ITAPE1
      READ(ITAPE) RA,BSTO
      READ(ITAPE) (NVARY(L),L=1,LRANG2)
      READ(ITAPE) ((ENDS(N,L),N=1,NVARY(L)),L=1,LRANG2)
      READ(ITAPE) ((COEFF(I,L),I=1,3),L=1,LRANG2)
C
      RETURN
      END
**********************************************************************
      SUBROUTINE SETCHN(LRGLLO)
C
C      FOR GIVEN LRGL AND PARITY TO CALCULATE
C      TOTAL NUMBER OF CHANNELS, ORDER OF CHANNELS FOR L AND STATE
C      AND THE FIRST CHANNEL NUMBER FOR L AND SET
C      SET UP THE ARRAY ICONCH TO CONVERT CHANNEL NUMBERS
C      IN THE NEW NO EXCHANGE CODE TO THOSE IN THE RMATRIX
C      CODE
C      SET UP NCONAT NO. OF CHANNELS COUPLED TO EACH STATE
C      AND L2P ARRAY LVALUE OF EACH CHANNEL ORDERED AS IN THE
C      RMATRIX CODE
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(ISDIM1=MZSET)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(ILDIM4=ILDIM3+ILDIM3)
C
      COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2),
     1              ICONCH(ITDIM3),KCONCH(ITDIM3),
     2              L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1),
     3              NSCHAN(IPDIM1),NSPREC(IPDIM1)
      COMMON/HPOINT/KRECH(0:ILDIM4)
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/SETL  /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2),
     1              NSETL(ISDIM1,ILDIM2)
      COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
C
      COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN
C
      DIMENSION ICHAN(ILDIM2,ITDIM1)
C
 1000 FORMAT(/,28X,'SUBROUTINE SETCHN'/28X,'-----------------'/)
 1001 FORMAT('      IS    KL NCHSET  IST ICHAN ISTCH')
 1002 FORMAT('  ',6I6)
 1003 FORMAT('  ICONCH(NCHAN)')
 1004 FORMAT(/'  TARGET SPIN (2S+1) =',I3/'  -----------------------')
 1005 FORMAT('    NO. OF CHANNELS=',I4)
 1006 FORMAT('    L2P COUPLED TO EACH CHANNEL')
 1007 FORMAT('    NO. OF CHANNELS COUPLED TO EACH STATE')
 1008 FORMAT('  NSPREC(NSP)')
C
      IF (NBUG8 .EQ. 1) THEN
      WRITE(IWRITE,1000)
      WRITE(IWRITE,1001)
      END IF
C
C      ZEROISE ICHAN ARRAY
C
      DO 6     I=1,NAST
      DO 7    KL=1,LTOT
            ICHAN(KL,I)=0
    7 CONTINUE
    6 CONTINUE
C
C
C      COUNT THE CHANNELS IN NCHAN
C
      NCHAN=0
C
C      LOOP OVER CONTINUUM ANGULAR MOMENTUM
C
      DO 1      KL=1,LTOT
C
C      LOOP OVER SETS COUPLED TO THIS L=LVAL(KL)
C
      DO 2       I=1,NSCOL(KL)
                IS=NSETL(I,KL)
C
C      IF THERE ARE NO STATES COMPOSED FROM THIS SET GO TO NEXT SET
C
         IF(NSTAT(IS) .EQ. 0) GO TO 2
C
C      STORE FIRST CHANNEL NUMBER INVOLVING THIS SET AND L
C
          NCHSET(IS,KL)=NCHAN+1
C
C      LOOP OVER STATES COMPOSED FROM THIS SET
C
      DO 3       J=1,NSTAT(IS)
               IST=NSETST(IS,J)
             NCHAN=NCHAN+1
         IF (ITDIM3 .LT. NCHAN) CALL RECOV1(21,NCHAN,ITDIM3)
C
C      STORE THE CHANNEL NUMBER FOR EACH L AND STATE TEMPORARILY
C
            ICHAN(KL,IST)=NCHAN
C
C      STORE THE STATE NUMBER COUPLED TO EACH CHANNEL
C
            ISTCH(NCHAN)=IST
      IF (NBUG8 .EQ. 1) THEN
      WRITE(IWRITE,1002)IS,KL,NCHSET(IS,KL),IST,ICHAN(KL,IST)
     1                 ,ISTCH(NCHAN)
      END IF
    3 CONTINUE
C
    2 CONTINUE
C
    1 CONTINUE
C
C      CONSTRUCT THE CHANNEL CONVERSION ARRAY
C      FOR ORDERING THE HAMILTONIAN ON DISC
C      STORE L2P AND NCONAT ARRAYS FOR ASYMPTOTIC PROGRAM
C
C      KRECH(2*LRGL+NPTY) WILL HOLD THE STARTING RECORD NO.
C      FOR EACH LRGL, NPTY COMBINATION
C      IN THIS VERSION THE STARTING RECORD IS 3 FOR EVERY
C      LRGL,NPTY COMBINATION TO SAVE DISC SPACE.
C
C      LENHAM(ISP) WILL HOLD THE LENGTH OF THE HAMILTONIAN FOR
C      SPIN ISP
C
C      NSPREC(ISP) WILL HOLD THE STARTING RECORD NO. FOR
C      THE HAMILTONIAN FOR SPIN ISP
C
      IF(NPTYIN.LT.0)THEN
      IRECH=LRGL+LRGL+NPTY
C     IF (IRECH .EQ. LRGLLO+LRGLLO) KRECH(IRECH)=3
      ELSE
      IRECH=LRGL
C     IF (IRECH .EQ. LRGLLO) KRECH(IRECH)=3
      ENDIF
      KRECH(IRECH)=3
      NSPREC(1)=KRECH(IRECH)+1
C
      KCOUNT=0
C
      DO 8   ISP=1,NSP
         ISPIN=NSPVAL(ISP)
         LENHAM(ISP)=0
         ICOUNT=0
C
      DO 4    IST=1,NAST
C
C      IF THE STATE DOES NOT HAVE CURRENT SPIN ENTER ZERO IN KCONAT
C
         LCOUNT=0
         IF (LSPN(IST) .NE. ISPIN) GO TO 41
C
      DO 5   KL=1,LTOT
            ICH=ICHAN(KL,IST)
         IF(ICH .NE. 0) THEN
           KCOUNT=KCOUNT+1
           ICOUNT=ICOUNT+1
           LCOUNT=LCOUNT+1
           ICONCH(ICH)=ICOUNT
           KCONCH(ICH)=KCOUNT
           L2PSPN(ICOUNT,ISP)=LVAL(KL)
           L=LVAL(KL)+1
           LENHAM(ISP)=LENHAM(ISP)+NVARY(L)
         END IF
    5 CONTINUE
C
   41    KCONAT(IST,ISP)=LCOUNT
    4 CONTINUE
C
         NSCHAN(ISP)=ICOUNT
C
         IF (ITDIM6 .LT. ICOUNT) CALL RECOV1(22,ICOUNT,ITDIM6)
C
         NSPREC(ISP+1)=(ICOUNT*(ICOUNT+1))/2 + NSPREC(ISP)
    8 CONTINUE
C
      KRECH(IRECH+1)=NSPREC(NSP+1)
C
      IF (NBUG8 .EQ. 1) THEN
      WRITE(IWRITE,1003)
      WRITE(IWRITE,1002)(ICONCH(I),I=1,NCHAN)
      END IF
      DO 9   ISP=1,NSP
         WRITE(IWRITE,1004)NSPVAL(ISP)
         WRITE(IWRITE,1005)NSCHAN(ISP)
         WRITE(IWRITE,1006)
         WRITE(IWRITE,1002)(L2PSPN(I,ISP),I=1,NSCHAN(ISP))
         WRITE(IWRITE,1007)
         WRITE(IWRITE,1002)(KCONAT(I,ISP),I=1,NAST)
    9 CONTINUE
      IF (NBUG8 .EQ. 1) THEN
      WRITE(IWRITE,1008)
      WRITE(IWRITE,1002)(NSPREC(ISP),ISP=1,NSP+1)
      END IF
C
      RETURN
      END
**********************************************************************
      SUBROUTINE STGHRD(ICONTN)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
C      READ STATE DATA AND FORM COMMON/STATE/
C
      CHARACTER*8 FILEH,FILAS,FILE1,FILE2,FILE3,FILEA,FILEB,FILEC,FILED
C
      PARAMETER(MACDIM=MZMAC)
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
C
      PARAMETER(LBUFD=IDIM3*IDIM3+IDIM4)
C
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/NBUG  /NBUGI(9)
      COMMON/SETS  /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1              NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2)
      COMMON/STATE1/NAST0,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
      COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1)
C
      DIMENSION NCST(ITDIM1), NTCON(ITDIM1)
      DIMENSION EN(ITDIM1),NORDER(ITDIM1),BIJ(ITDIM1,ISDIM2),
     1          LL1(ITDIM1),LSPN1(ITDIM1),LPTY1(ITDIM1)
C
      DIMENSION TITLE(18)
      CHARACTER*4 TITLE
C NRB
      NAMELIST/STG3A/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8
     X,NBUG9,IPOLPH,RAD,ISORT
      NAMELIST/STG3B/NBUT,NCONT,IDIAG,NAST,INAST,TOLER,NEST
C
 2000 FORMAT (18A4)
C2001 FORMAT (5F14.7)
C2002 FORMAT (9I5)
C
 1000 FORMAT (TR31,'***********'/TR31,'**       **'/TR31,'** NXHAM **'/
     1        TR31,'**       **'/TR31,'***********'//)
 1002 FORMAT (////TR28,' SUBROUTINE STGHRD'/TR28,' -----------------')
 1003 FORMAT (//' OUTPUT CHANNEL                UNIT NUMBER=',I3/
     1          ' INPUT CHANNEL                 UNIT NUMBER=',I3/
     2          ' TARGET ANGULAR FILE--INTEGER  UNIT NUMBER=',I3,
     3          ' NAME= ',A8/
     3          '                      REAL     UNIT NUMBER=',I3,
     3          ' NAME= ',A8/
     4          ' RACAH ANGULAR FILE            UNIT NUMBER=',I3,
     3          ' NAME= ',A8/
     5          ' RADIAL INTEGRAL FILE          UNIT NUMBER=',I3,
     3          ' NAME= ',A8/
     6          ' RADIAL POINTER FILE           UNIT NUMBER=',I3,
     3          ' NAME= ',A8/
     7          ' SEQENTIAL FROM NXRAD          UNIT NUMBER=',I3,
     3          ' NAME= ',A8/
     8          ' UNDIAG. HAMILTONIAN FILE      UNIT NUMBER=',I3,
     3          ' NAME= ',A8/
     9          ' OUTPUT FOR ASYMPTOTIC PROGRAM UNIT NUMBER=',I3,
     3          ' NAME= ',A8)
 1004 FORMAT ('  DEBUG PARAMETERS'/12I5)
 1005 FORMAT ('  NUMBER OF ATOMIC STATES ',I4)
 1006 FORMAT ('  STATE',I3,' L=',I3,' SPIN=',I3,' PARITY=',I1)
 1007 FORMAT ('  NUMBER OF STATES WITH EACH SYMMETRY'/'  ',9 I4)
 1008 FORMAT ('  STATES INDEX FOR EACH SYMMETRY')
 1009 FORMAT ('  ',9 I4)
 1010 FORMAT ('   AIJ=',8F14.7,4(1X/6X,8F14.7))
 1011 FORMAT ('  WARNING - THE NORMALIZATION OF THE STATE IS ',
     1           F14.7, /'   THE STATE IS BEING RENORMALISED')
 1012 FORMAT ('   ENAT=',2F14.7,I8)
 1013 FORMAT ('   EXPERIMENTAL ENERGIES',4(1X,/6X,8F14.7))
C
 3001 FORMAT ('  ****  ERROR IN DATA  ****'/
     1        '  STATE ',I3,' SHOULD CONTAIN ',I3,' CONFIGURATIONS '/
     2        '  IT EXPECTS ',I3,' COEFFICIENTS')
C
      IREAD=5
      IWRITE=6
      JBUFD=26
      ITAPE3=10
      FILEA='ANG1.DAT'
      FILEB='ANG2.DAT'
      FILEC='ANG3.DAT'
      FILE1='RAD1.DAT'
      FILE2='RAD2.DAT'
      FILE3='RAD3.DAT'
      FILEH='HAM1.DAT'
C     FILAS='H.DAT'
      FILAS='H.DAT'
C
C      ZEROISE KSPIN ARRAY USED IN SORTING AND STORING
C      STATE SPINS
C
      DO 1   I=1,IPDIM1
         KSPIN(I)=0
    1 CONTINUE
C
C      ZEROISE NSTAT ARRAY USED IN COUNTING STATES BELONGING
C      TO EACH SET
C
      DO 2   IS=1,NSET
         NSTAT(IS)=0
    2 CONTINUE
C
      IF (ICONTN .EQ. 1) THEN
         INX2=36
      END IF
C
C      OPEN HAMILTONIAN AND ASYMPTOTIC FILES
C
      OPEN(UNIT=JBUFD,FILE=FILEH,ACCESS='DIRECT',STATUS='UNKNOWN',
     1     FORM='UNFORMATTED',RECL=MACDIM*LBUFD)
      OPEN(UNIT=ITAPE3,FILE=FILAS,ACCESS='SEQUENTIAL',STATUS='UNKNOWN',
     1     FORM='UNFORMATTED')
C
      WRITE (IWRITE,1000)
      WRITE (IWRITE,1002)
      WRITE (IWRITE,1003) IWRITE,IREAD,JBUFIR,FILEA,JBUFFV,FILEB,
     1                    JBUFFZ,FILEC,JBUFF1,FILE1,JBUFF2,FILE2,
     2                    ITAPE1,FILE3,JBUFD,FILEH,ITAPE3,FILAS
      IF (ICONTN .NE. 1) THEN
        READ (IREAD,*) (NBUGI(I),I=1,9)
      END IF
C
      WRITE (IWRITE,1004) (NBUGI(I),I=1,9)
      IF (ICONTN .EQ. 1) THEN
         READ (INX2) (LL1(I),LSPN1(I),LPTY1(I),I=1,NAST0)
      ELSE
      READ (IREAD,*) NAST0
      END IF
      WRITE (IWRITE,1005) NAST0
      IF (ITDIM1 .LT. NAST0) CALL RECOV1(19,NAST0,ITDIM1)
C
      DO 3   IST=1,NAST0
         IF(ICONTN .EQ. 0)READ (IREAD,*)LL1(IST),LSPN1(IST),LPTY1(IST)
         WRITE (IWRITE,1006) IST,LL1(IST),LSPN1(IST),LPTY1(IST)
         LST=LL1(IST)
         LSPNST=LSPN1(IST)
         LPTYST=LPTY1(IST)
      IF (IPDIM1 .LT. LSPNST) CALL RECOV1(23,LSPNST,IPDIM1)
C
C      RECORD THAT THIS SPIN VALUE OCCURS
C
         KSPIN(LSPNST)=1
C
C      FIND TO WHICH SET THIS STATE BELONGS
C
      DO 31   IS=1,NSET
         IF (LST .EQ. LSET(IS)  .AND.
     1       LSPNST .EQ. ISPSET(IS)  .AND.
     2       LPTYST .EQ. IPISET(IS)) THEN
            NSTAT(IS)=NSTAT(IS)+1
      IF (ITDIM2 .LT. NSTAT(IS)) CALL RECOV1(20,NSTAT(IS),ITDIM2)
            NSETST(IS,NSTAT(IS))=IST
            NCST(IST)=NCSET(IS)
         END IF
   31 CONTINUE
C
    3 CONTINUE
C
      WRITE (IWRITE,1007)(NSTAT(IS), IS=1,NSET)
      WRITE (IWRITE,1008)
C
      DO 4   IS=1,NSET
         WRITE (IWRITE,1009)(NSETST(IS,JST), JST=1,NSTAT(IS))
    4 CONTINUE
C
C      RECORD SPIN VALUES OCCURRING
C
      ICOUNT=0
C
      DO 5   ISP=1,IPDIM1
         IF (KSPIN(ISP) .NE. 0) THEN
            ICOUNT=ICOUNT+1
            NSPVAL(ICOUNT)=ISP
            KSPIN(ISP)=ICOUNT
      END IF
    5 CONTINUE
C
      NSP=ICOUNT
C
C      READ NO. OF CONFIGURATIONS IN EACH STATE AND COMPARE
C      WITH THE NO. OF CONFIGURATIONS IN THE RELEVANT SET
C      BEFORE READING THE CONFIGURATION COEFFICIENTS AND
C      ENERGY FOR EACH STATE.  TEST NORMALISATION
C
      IF (ICONTN .EQ. 1) THEN
         READ (INX2)(NTCON(IST), IST=1,NAST0)
      ELSE
         READ (IREAD,*)(NTCON(IST), IST=1,NAST0)
      END IF
      PT01=TENTH*TENTH
C
      DO 6   IST=1,NAST0
         NCONF=NCST(IST)
         IF (NCONF .NE. NTCON(IST)) THEN
            WRITE (IWRITE,3001) IST,NCONF,NTCON(IST)
            STOP
         END IF
         IF (ICONTN .EQ. 1) THEN
            READ(INX2)(BIJ(IST,IC), IC=1,NCONF),EN(IST)
         ELSE
         READ (IREAD,*) (BIJ(IST,IC), IC=1,NCONF),EN(IST)
         END IF
         WRITE (IWRITE,1010) (BIJ(IST,IC), IC=1,NCONF)
         ANORM=ZERO
      DO 61   IC=1,NCONF
         ANORM=ANORM+BIJ(IST,IC)**2
   61 CONTINUE
C
      IF (ABS(ANORM-ONE) .GE. PT01) THEN
         WRITE (IWRITE,1011) ANORM
         DO 62   IC=1,NCONF
            BIJ(IST,IC)=BIJ(IST,IC)/SQRT(ANORM)
   62    CONTINUE
         WRITE (IWRITE,1010) (BIJ(IST,IC), IC=1,NCONF)
      END IF
C
      T=TWO*(EN(IST)-EN(1))
      WRITE (IWRITE,1012) EN(IST),T
    6 CONTINUE
C
      IF (ICONTN .EQ. 1) THEN
C        FILED='DATAR3'
C
C OUR 'dstg3' DATASET, USUALLY REDEFINED ON INPUT TO UNIT 5 IN STG3
C
         FILED='dstg3'
         OPEN (37,FILE=FILED,STATUS='OLD',ACCESS='SEQUENTIAL')
         READ (37,2000) (TITLE(I),I=1,18)
C        READ (37,2002) (IDUMMY(I),I=1,9)
C        READ (37,2002) (IDUMMY(I),I=1,9)
C        READ (37,2002) (IDUMMY(I),I=1,3)
C        READ (37,2002) N1,N2,N3,N4,NEST,N5
         NEST=0
         NAST=0
         ISORT=0
C
         READ(37,STG3A)
         READ(37,STG3B)
         IF(NEST.EQ.0)NEST=NAST
      ELSE
        ISORT=0
        READ(IREAD,*) NEST
      END IF
      NEST0=IABS(NEST)
C
      IF (NEST .NE. 0) THEN
        IF (ICONTN .EQ. 1) THEN
          READ (37,*)(ENAT(I),I=1,NEST0)
        ELSE
          READ(IREAD,*)(ENAT(IST),IST=1,NEST0)
        END IF
        WRITE(IWRITE,1013)(ENAT(IST),IST=1,NEST0)
C
        IF(ISORT.GT.0)CALL ORDER(EN,NAST0,NORDER)
C
        DO I=1,NEST0
C RYDBERGS RELATIVE TO GROUND
          IF(NEST.GT.0)THEN
            ENAT(I)=HALF*ENAT(I)+EN(1)
          ELSE
            WRITE(6,*)' NAST/NEST.LT.0 (MERGE) NOT IN NX CODE YET'
            STOP' NAST/NEST.LT.0 (MERGE) NOT IN NX CODE YET'
          ENDIF
C
C B.P. CODES NOW USE NAST/NEST.LT.0 TO MERGE
C AND EST/EN(1).NE.0 FOR ABSOLUTE A.U.
C AND EST/EN(1).GT.0 RELATIVE RYD
C AND EST/EN(1).LT.0 RELATIVE CM-1
COLD      IF(NEST.LT.0)ENAT(I)=HALF*ENAT(I)+EN(1)
          I0=I
          IF(ISORT.GT.0)I0=NORDER(I)
          EN(I0)=ENAT(I)
        ENDDO
      END IF
C
C      ORDER THE STATES IN ORDER OF ASCENDING ENERGY
C
      CALL ORDER(EN,NAST0,NORDER)
C
      WRITE(IWRITE,'(/'' REORDERED BY ENERGY'',/)')
C
      DO 8   IST=1,NAST0
         JST=NORDER(IST)
         ENAT(IST)=EN(JST)
         LL(IST)=LL1(JST)
         LSPN(IST)=LSPN1(JST)
         LPTY(IST)=LPTY1(JST)
         NCONF=NCST(JST)
C
      DO 9   IC=1,NCONF
         AIJ(IST,IC)=BIJ(JST,IC)
    9 CONTINUE
C
         WRITE (IWRITE,1010) (AIJ(IST,IC), IC=1,NCONF)
    8 CONTINUE
C
      DO 10  IST=1,NAST0
         WRITE (IWRITE,1006) IST,LL(IST),LSPN(IST),LPTY(IST)
   10 CONTINUE
C
      DO 11  IST=1,NAST0
         T=TWO*(ENAT(IST)-ENAT(1))
         WRITE (IWRITE,1012) ENAT(IST),T,IST
   11 CONTINUE
C
      DO 12  IS=1,NSET
C
      DO 13  IST=1,NSTAT(IS)
         JST=NSETST(IS,IST)
C
      DO 14  KST=1,NAST0
         IF (NORDER(KST) .EQ. JST) THEN
            NSETST(IS,IST)=KST
            GO TO 13
         END IF
   14 CONTINUE
C
   13 CONTINUE
         WRITE (IWRITE,1009)(NSETST(IS,JST), JST=1,NSTAT(IS))
   12 CONTINUE
C
      IF (ICONTN .EQ. 1) THEN
         CLOSE (INX2)
         CLOSE (37)
      END IF
C
      RETURN
      END
**********************************************************************
      SUBROUTINE STGRRD

*     THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES
*     CERTAIN VARIABLES

      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
      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,IDIM5=MZNPO)
      PARAMETER(IDIM11=MZSLT)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
C
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(LBUFF2=IDIM1*IDIM1*KDIM9)
      PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
C

      DIMENSION MAXNLG(IDIM2),MAXPN(IDIM2)
      DIMENSION IRAD(IDIM11),ZE(IDIM11),C(IDIM11)

*     FORMAT STATEMENTS

2000  FORMAT(A)

1002  FORMAT(//'  ',10X,'STGRAD DATA')
1053  FORMAT(//12X,'THIS IS A NO-EXCHANGE CALCULATION'//)
1054  FORMAT(//12X,'THIS CALCULATION INCLUDES EXCHANGE'//)
1003  FORMAT(//' OUTPUT CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  SEQUENTIAL    UNIT NUMBER =',I3/
     *         ' OUTPUT SEQUENTIAL    UNIT NUMBER =',I3/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR RADIAL INT
     *EGRALS)'/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR CONTINUUM-
     *CONTINUUM POINTER ARRAYS)'/)
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,'  LAM =',I2,'  IBC =',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)

3001  FORMAT('  **** NOT YET ALLOWED IN THIS VERSION ****')


*     SET THE CHANNEL UNIT NUMBERS

      IREAD=5
      IWRITE=6
      JBUFF1=23
      JBUFF2=24
      ITAPE1=25
      FILE1='RAD1.DAT'
      FILE2='RAD2.DAT'
      FILE3='RAD3.DAT'
*     READ IN THE BASIC DATA

C
C      OPEN RADIAL FILES
C
      OPEN(UNIT=JBUFF1,FILE=FILE1,ACCESS='DIRECT',STATUS='OLD',
     1     FORM='UNFORMATTED',RECL=MACDIM*LBUFF1)
      OPEN(UNIT=JBUFF2,FILE=FILE2,ACCESS='DIRECT',STATUS='OLD',
     1     FORM='UNFORMATTED',RECL=MACDIM*LBUFF2)
      OPEN(UNIT=ITAPE1,FILE=FILE3,ACCESS='SEQUENTIAL',STATUS='OLD',
     1     FORM='UNFORMATTED')

CW    WRITE(IWRITE,1002)
CW    WRITE (IWRITE,1053)

      READ (IREAD, *)    NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,
     *                   NBUG8,NBUG9

      READ (IREAD, *) IPSEUD,NUMERC
      IF (IPSEUD.EQ.0) THEN
CW       WRITE(IWRITE,1005)
      ELSE
CW       WRITE(IWRITE,1006)
      ENDIF
      IF (NUMERC.EQ.0) THEN
CW       WRITE(IWRITE,1007)
      ELSE
CW       WRITE(IWRITE,1008)
      ENDIF

CW    WRITE(IWRITE,1009)
      READ (IREAD, *)    NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM,IBC
CW    WRITE(IWRITE, 1010)NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM,IBC
      IF (IDIM1 .LT. LRANG1) CALL RECOV1(1,LRANG1,IDIM1)
      IF (IDIM3 .LT. NRANG2) CALL RECOV1(3,NRANG2,IDIM3)
      IF (IDIM4 .LT. LAMAX) CALL RECOV1(4,LAMAX,IDIM4)
      READ (IREAD, *)    (MAXNHF(I),I=1,LRANG1)
CW    WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1)
      READ (IREAD, *)    (MAXNLG(I),I=1,LRANG1)
CW    WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1)

*     READ IN DATA DEFINING THE STATIC POTENTIAL OF THE GROUND STATE

CW    WRITE(IWRITE,1013)
      DO 10,I=1,LRANG1
         READ (IREAD, *) ISTAT
10    CONTINUE
C
      IF (NUMERC.EQ.0) THEN
*        ANALYTIC ORBITALS ORBITALS ARE TO BE USED
*        READ IN COEFFICIENTS DEFINING THE ANALYTIC ORBITALS
         DO 20,L=1,LRANG1
            DO 30,N=L,MAXNHF(L)
               READ (IREAD, *) M
               READ (IREAD, *) (IRAD(IM),IM=1,M)
               READ (IREAD, *) (ZE(IM),IM=1,M)
               READ (IREAD, *) (C(IM),IM=1,M)
30          CONTINUE
20       CONTINUE
         IF (IBC.NE.0) THEN
            READ (IREAD, *) RA,BSTO
         ENDIF
CW       WRITE(IWRITE,1018) RA,BSTO
      ENDIF
      IF (IPSEUD.EQ.0) THEN
*        THIS IS A NON-MODEL POTENTIAL CALCULATION.
         IF (IBC.GT.1) THEN
*           READ IN INTEGRATION MESH
            READ (IREAD, *) NIX
            READ (IREAD, *) HINT
            READ (IREAD, *) IHX
            READ (IREAD, *) IRX
         ENDIF
*        INITIALISE MAXNC IN NON-MODEL POTENTIAL
*        CALCULATION.
         DO 150,I=1,LRANG1
            MAXNC(I)=I-1
150      CONTINUE
C WHAT ABOUT MODEL?
      ENDIF
CW    WRITE(IWRITE,1023) (MAXNC(I),I=1,LRANG1)
C
C      READ THE MAXIMUM ENERGY FOR THE BUTTLE CORRECTION
C
      READ(IREAD, *) EK2MAX,MJS
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
C
      IF (LRANG2.GT.LRANG1) THEN
         DO 170,I=LRANG1+1,LRANG2
            MAXNHF(I)=I-1
            MAXNLG(I)=I-1
170      CONTINUE
      ENDIF



      RETURN
      END
**********************************************************************
      SUBROUTINE STGRRX
*     THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES
*     CERTAIN VARIABLES USING RMATRX OUTPUT FILE
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
      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,IDIM5=MZNPO)
      PARAMETER(IDIM6=MZNIX)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
C
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(LBUFF2=IDIM1*IDIM1*KDIM9)
      PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON /NRB/MAXC
C
      DIMENSION MAXNLG(IDIM2),MAXPN(IDIM2),IHX(IDIM6),IRX(IDIM6)
*     FORMAT STATEMENTS
2000  FORMAT(A)
1002  FORMAT(//'  ',10X,'STGRAD DATA')
1053  FORMAT(//12X,'THIS IS A NO-EXCHANGE CALCULATION'//)
1054  FORMAT(//12X,'THIS CALCULATION INCLUDES EXCHANGE'//)
1003  FORMAT(//' OUTPUT CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  CHANNEL       UNIT NUMBER =',I3/
     *         ' INPUT  SEQUENTIAL    UNIT NUMBER =',I3/
     *         ' OUTPUT SEQUENTIAL    UNIT NUMBER =',I3/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR RADIAL INT
     *EGRALS)'/
     *         ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR CONTINUUM-
     *CONTINUUM POINTER ARRAYS)'/)
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,'  LAM =',I2,'  IBC =',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)
3001  FORMAT('  **** NOT YET ALLOWED IN THIS VERSION ****')
*     SET THE CHANNEL UNIT NUMBERS
      IREAD=5
      IWRITE=6
      JBUFF1=23
      JBUFF2=24
      ITAPE1=25
      INX1=35
C
      FILE1='RAD1.DAT'
      FILE2='RAD2.DAT'
      FILE3='RAD3.DAT'
*       OPEN THE DIRECT ACCESS FILES
      OPEN(UNIT=JBUFF1,FILE=FILE1,ACCESS='DIRECT',
     *    STATUS='OLD',FORM='UNFORMATTED',RECL=MACDIM*LBUFF1)
      OPEN(UNIT=JBUFF2,FILE=FILE2,ACCESS='DIRECT',
     *    STATUS='OLD',FORM='UNFORMATTED',RECL=MACDIM*LBUFF2)
      OPEN(ITAPE1,FILE=FILE3,ACCESS='SEQUENTIAL',
     *    STATUS='OLD',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    1     FORM='UNFORMATTED')
C
CW    WRITE(IWRITE,1002)
C
CW    WRITE(IWRITE,1009)
      LAM=0
      READ (INX1) LRANG1,NRANG2,LAMAX,NIX,NPTS,NRKPTS
      READ (INX1) NELC,NZ,LFIXN,IPSEUD
      IF(MAXC.GT.1)NRANG2=MAXC
CW    WRITE(IWRITE,1010) NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM
      IF (IDIM1 .LT. LRANG1) CALL RECOV1(1,LRANG1,IDIM1)
      IF (IDIM3 .LT. NRANG2) CALL RECOV1(3,NRANG2,IDIM3)
      IF (IDIM4 .LT. LAMAX) CALL RECOV1(4,LAMAX,IDIM4)
C
      READ (INX1) (MAXNHF(I),I=1,LRANG1), (MAXNLG(I),I=1,LRANG1)
CW     WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1)
CW     WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1)
C
      DO 150  I=1,LRANG1
         MAXNC(I)=I-1
 150  CONTINUE
C SHOULD MAXNC BE REDEFINED FOR MODEL POTENTIAL?
      IF(IPSEUD.GT.0)THEN
        READ(INX1)HINT,(IHX(I),I=1,NIX),(IRX(I),I=1,NIX)
C ANSWER APPEARS TO BE NO
C       READ(INX1)(MAXNC(I),I=1,LRANG1)
      ENDIF
*     DEFINE THE MAXPN ARRAY . NOT USED?
      DO 160,I=1,LRANG1
         MAXPN(I)=MAXNHF(I)-MAXNC(I)+I-1
160   CONTINUE
*     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
*     INITIALISE THE VARIABLES LAMBC AND LAMCC (NOT USED HERE)
      IF (LAM.GE.2) THEN
         LAMBC=LAMAX
      ELSE
          LAMBC=0
      ENDIF
      IF (LAM.GE.3) THEN
         LAMCC=LAMAX
      ELSE
          LAMCC=0
      ENDIF
C
      CLOSE (INX1)
C
      RETURN
      END
*=======================================================================
      SUBROUTINE STHMAT(LRGLLO)
C
C      CALCULATES DIRECT TERMS IN THE HAMILTONIAN MATRIX
C      AND STORES THEM ON DISC IN RECORDS NRANG2 X NRANG2
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX)
      PARAMETER(ICDIM2=MZORB)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS)
C
      PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF)
C
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
      PARAMETER(KDIM9=IDIM1+IDIM1-1)
      PARAMETER(ILDIM2=ILDIM1+ILDIM1-1)
      PARAMETER(ILDIM4=ILDIM3+ILDIM3)
      PARAMETER(INDIM2=IDIM3*IDIM3)
      PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1)
      PARAMETER(ITDIM4=ITDIM2*ITDIM2)
      PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 +
     1                 ILDIM1*ISDIM2*(ICDIM2-1))
      PARAMETER(IHDIM1=ITDIM6*IDIM3)
      PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2)
      PARAMETER(IHDIM3=INDIM2*ITDIM4)
C     PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1))
      PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1)
      PARAMETER(ICDIM=(IADIM+IBDIM+1)/2)
      PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM)
     1               -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM)
      PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1)
C
      COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4),
     1              IRHSGL(IVDIM1)
      COMMON/ASYMPR/CRL(ITDIM1,ITDIM1,IDIM4)
      COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2),
     1              ICONCH(ITDIM3),KCONCH(ITDIM3),
     2              L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1),
     3              NSCHAN(IPDIM1),NSPREC(IPDIM1)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/HPOINT/KRECH(0:ILDIM4)
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1),
     1              R1CC(LBUFF1),RKCCD(IRDIM2)
      COMMON/NBUG  /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,
     1              NBUG9
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2,
     1              IDIM2),ICCDIF,ICCLNG
      COMMON/SETL  /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2),
     1              NSETL(ISDIM1,ILDIM2)
      COMMON/SETS  /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1),
     1              NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2)
      COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
      COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1)
C
      COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN
C
      PARAMETER(INDIM3=INDIM2+IDIM4)
C
      DIMENSION CMAT(IDIM3,IDIM3),HSMAT(IHDIM3)
      DIMENSION BUF1D(INDIM3),BUF2D(INDIM3)
      DIMENSION R1(INDIM2)
      DIMENSION ZB(ILDIM1)
      DIMENSION CF(ITDIM4,IDIM4)
C
C
 1000 FORMAT(/,28X,'SUBROUTINE STHMAT'/28X,'-----------------'/)
 1001 FORMAT('       L  LP  IS  JS  IC  JC IST JST')
 1002 FORMAT(4X,8I4)
 1003 FORMAT(2X,9E13.6)
 1004 FORMAT('  ONE-ELECTRON INTEGRAL')
 1005 FORMAT('  Z-COEFFICIENTS')
 1006 FORMAT('  CF MATRIX OF ASYMPTOTIC COEFFS FOR SETS ',I3,I3)
C
 3000 FORMAT('  **** STOP IN STHMAT ****'/
     1       '  ERROR IN USING VSH FILE FOR SETS ', 2I4)
C
      IF (NBUG6 .EQ. 1) THEN
         WRITE(IWRITE,1000)
         WRITE(IWRITE,1001)
      END IF
C
      IF(NPTYIN.LT.0)THEN
      LTEST=LRGLLO+LRGLLO
      LRGLPI=LRGL+LRGL+NPTY
      ELSE
      LTEST=LRGLLO
      LRGLPI=LRGL
      ENDIF
C
C      AT FIRST ENTRY READ POINTER ARRAYS FOR ANGULAR AND RK INTEGRALS
C
      IF (LRGLPI .EQ. LTEST) THEN
         READ (JBUFFZ, REC=2) L,KRECZ
         READ (ITAPE1)((NCCD(I,J),I=1,LLPDIM),J=1,LRANG2)
      END IF
C
C      WRITE POINTER BLOCK FOR THE HAMILTONIAN MATRICES FOR
C      THIS LRGL AND PARITY
C
      IRECH=LRGLPI
      NRECH=KRECH(IRECH)
      WRITE(JBUFD, REC=NRECH) NSP,(NSPVAL(I),NSCHAN(I),NSPREC(I),
     1                        LENHAM(I),I=1,NSP)
      IF (LTOT .EQ. 0) RETURN
C
C      READ FIRST BLOCK OF Z COEFFICIENTS FOR THIS LRGL AND PARITY
C
      IRECZ=LRGLPI
      NRECZ=KRECZ(IRECZ)
      READ (JBUFFZ, REC=NRECZ)BUFFZ
      IZCON=1
C
C      SET RECORD NUMBERS TO ZERO READY FOR INITIAL ENTRY
C      INTO ROUTINES TO READ RADIAL INTEGRALS FROM DISC
C
      NRECR1=0
      NRECIC=0
      NRECRK=0
      NB=0
      NRECV1=0           !NRB
      NRECV2=0           !NRB
C
      LASTL=LVAL(LTOT)+1
C
C      LOOP OVER CONTINUUM ANGULAR MOMENTUM
C
      DO 2   KL=1,LTOT
              L=LVAL(KL)
           NRA1=NVARY(L+1)
C
C      LOCATE ONE ELECTRON INTEGRAL
C
         CALL LOC1CC(L+1,LASTL,NRECR1,R1)
C
      IF (NBUG6 .EQ. 1) THEN
         WRITE(IWRITE,1004)
         NRA1SQ=NRA1*NRA1
         WRITE(IWRITE,1003)(R1(NNP),NNP=1,NRA1SQ)
      END IF
C
      DO 3  KLP=KL,LTOT
             LP=LVAL(KLP)
         NRA2=NVARY(LP+1)
         NRASQR=NRA1*NRA2
C
C      READ FROM DISC ALL CONTINUUM-CONTINUUM RK INTEGRALS
C      FOR THIS L,LP COMBINATION
C
      CALL RDRKCC(L,LP,NRECIC,NB,NRECRK)
C
C      LOOP OVER THE SETS COUPLED TO L AND LP
C
      DO 4    I=1,NSCOL(KL)
             IS=NSETL(I,KL)
         ICHAN1=NCHSET(IS,KL)
           LCFG=LSET(IS)
C
C      FIND THE SPIN, THE NO. OF CHANNELS IN THE
C      HAMILTONIAN FOR THIS SPIN AND THE RECORD
C      IN THE FILE AT WHICH IT STARTS
C
          ISPIN=ISPSET(IS)
            ISP=KSPIN(ISPIN)
            NCH=NSCHAN(ISP)
          NSREC=NSPREC(ISP)-1
C
            JLO=1
           IF(L .EQ. LP) JLO=I
C
      DO 5    J=JLO,NSCOL(KLP)
             JS=NSETL(J,KLP)
         JCHAN1=NCHSET(JS,KLP)
         LCFGP=LSET(JS)
         LAMCFG=LCFG+LCFGP
         JSPIN=ISPSET(JS)
C
C      LAMSIN IS THE INDICATOR TO SAY SETS IS, JS CAN INTERACT
C      IF THEY CANNOT IT WILL BE SET TO -999 AND USED TO SKIP
C      UNNECESSARY LOOPS AND ENSURE ZERO SUB MATRICES ARE
C      STORED IN THE HAMILTONIAN ON DISC IN THE APPROPRIATE
C      CHANNEL POSITIONS
C
         LAMSIN=0
C
C      SUBROUTINE DIRANG ASSUMES JS .GE. IS
C      TO BE CONSISTENT IF IS .GT. JS THE SETS MUST BE REVERSED
C
          IF(IS .GT. JS) THEN
            KIS=JS
            KJS=IS
          ELSE
            KIS=IS
            KJS=JS
          END IF
C
          LAMLU=LAMIJ(KIS,KJS)
C
C      TEST WHETHER SPIN ALLOWS SETS TO COUPLE. IF NOT
C      JUMP TO END OF LOOP
C
         IF (ISPIN .NE. JSPIN) GO TO 5
C
C      TEST WHETHER SETS CANNOT COUPLE FOR SOME OTHER REASON
C      IF THEY CAN LAMIJ CONTAINS THE LIMITS ON LAMBDA IMPOSED
C      BY SHELL COUPLING BETWEEN CONFIGURATIONS IN THE SETS
C
          IF(LAMLU .EQ. -999) THEN
            LAMSIN=-999
            GO TO 50
          END IF
C
          LAMLO=LAMLU/100
          LAMUP=MOD(LAMLU,100)
C
C      CALCULATE LAMST,LAMFIN THE LIMITS ON LAMBDA IN THE STORAGE
C      OF THE Z COEFFICIENTS
C
          LAMST=MAX(ABS(L-LP),LAMLO)
          LAMFIN=MIN( L+LP ,LAMUP)
          IF(LAMST .GT. LAMFIN) THEN
             LAMSIN=-999
             GO TO 50
          END IF
C
C      TRANSFER TARGET ANGULAR INTEGRALS FOR THIS SET
C      COMBINATION TO VSH
C SURELY WE MUST INITIALIZE/SAVE NRECV1,2 TO BE SAFE - NRB
C
      CALL RDVSH(KIS,KJS,NRECV1,NRECV2,NVSHEL)
      ICOUNT=1
C
C      TRANSFER Z COEFFICIENTS FOR THIS SET COMBINATION TO ZB
C
           IZ=1
      DO  51    LAM=LAMST,LAMFIN,2
               ZB(IZ)=BUFFZ(IZCON)
                 IZ=IZ+1
              IZCON=IZCON+1
          IF (IZCON .GT. LBUFFZ) THEN
             NRECZ=NRECZ+1
             READ (JBUFFZ, REC=NRECZ) BUFFZ
             IZCON=1
          END IF
   51 CONTINUE
C
      IF (NBUG6 .EQ. 1) THEN
      WRITE(IWRITE,1005)
      WRITE(IWRITE,1003)(ZB(IIZ),IIZ=1,IZ)
      END IF
C
C
C      CALCULATE THE ASYMPTOTIC COEFFICIENTS
C      CF(IJCHAN,LAMBDA) FOR CHANNELS IN THE
C      SETS IS,JS.  LAMBDA .GE. 1
C
C      ZEROISE CF ARRAY
C
   50 DO 52   ICFS=1,ITDIM4
      DO 53   LAM=1,LAMAX
         CF(ICFS,LAM)=ZERO
   53 CONTINUE
   52 CONTINUE
      ICFS=0
          IF (LAMLU .GT. 0  .AND.  LAMSIN .EQ. 0) THEN
C
C      CALCULATE LIMITS ON LAMBDA CONSISTENT WITH Z-COEFFS
C
             IF (LAMST .EQ. 0) THEN
                LAMS=2
                IZCONT=1
             ELSE
                LAMS=LAMST
                IZCONT=0
             END IF
             LAMF=MIN(LAMAX,LAMFIN)
             DO  501  IST=1,NSTAT(IS)
                 ISTAT=NSETST(IS,IST)
                 IF (ICHAN1 .EQ. JCHAN1) THEN
                    JSTLO=IST
                 ELSE
                    JSTLO=1
                 END IF
             DO  502  JST=JSTLO,NSTAT(JS)
                 JSTAT=NSETST(JS,JST)
                 IF (JSTAT .LT. ISTAT) THEN
C
C      INTERCHANGE STATE NUMBERS SINCE ONLY UPPER HALF
C      OF THE CRL MATRIX IS STORED
C
                    KJSTAT=ISTAT
                    KISTAT=JSTAT
                 ELSE
                    KJSTAT=JSTAT
                    KISTAT=ISTAT
                 END IF
                 IZ=IZCONT
                 ICFS=ICFS+1
             DO  503  LAM=LAMS,LAMF,2
                 IZ=IZ+1
                 CF(ICFS,LAM)=2*CRL(KISTAT,KJSTAT,LAM)*ZB(IZ)
  503        CONTINUE
  502        CONTINUE
  501        CONTINUE
                 IF (NBUG7 .EQ. 1) THEN
                 WRITE(IWRITE,1006) IS,JS
                 WRITE(IWRITE,1003) ((CF(IJCHAN,LAM), LAM=1,LAMAX),
     1                               IJCHAN=1,ICFS)
                 END IF
          END IF
C
C
C      COLLECT TOGETHER THE TERMS FROM PAIRING OF CONFIGURATIONS
C      WITHIN THE SETS KIS,KJS
C
C      ZEROISE HSMAT ARRAY
C
      DO 1   IHS=1,IHDIM3
         HSMAT(IHS)=ZERO
    1 CONTINUE
C
C      LOOP OVER CONFIGURATIONS IN SET KIS
C      IF THESE SETS CANNOT INTERACT PERFORM THE LOOPS ONLY
C      ONCE TO SET UP THE CHANNEL POSITIONS
C      FIRST SET THE UPPER LIMITS TO THE CONFIGURATION LOOPS
C
      IF (LAMSIN .EQ. 0) THEN
         KICSUP=NCSET(KIS)
         KJCSUP=NCSET(KJS)
      ELSE
         KICSUP=1
         KJCSUP=1
      END IF
C
      DO 6   KICS=1,KICSUP
              IF(KIS .EQ. KJS) THEN
                 JCSLO=KICS
              ELSE
                 JCSLO=1
              END IF
C
C      LOOP OVER CONFIGURATIONS IN SET KJS
C
      DO 7   KJCS=JCSLO,KJCSUP
         IF (LAMSIN .EQ. 0) THEN
C
C      THE CASE OF CONFIGURATIONS HAVING THE SAME SHELL STRUCTURE
C      IS TREATED SEPARATELY
C
          IF (IRHSGL(ICOUNT) .LT. 0) THEN
      CALL CMATII(L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT)
           ELSE
      CALL CMATIJ(L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT,IND)
C
C      THE INDICATOR IND IS SET ZERO IF IC,JC DO NOT COUPLE
C*************
C      THIS MAY CAUSE TROUBLE IF ALL PAIRS CANNOT COUPLE
C      EG. A ZERO SUB-MATRIX
C************
C
           IF(IND .EQ. 0) GO TO 7
           END IF
C
C      ADD THE CMAT, MULTIPLIED BY ITS APPROPRIATE COEFFICIENT
C      AIJ INTO THE RELEVANT CHANNEL POSITIONS IN THE DIRECT
C      PART OF THE HAMILTONIAN MATRIX DMAT
C      WHEN IS .GT. JS THE SETS HAVE BEEN INTERCHANGED.
C      THE CONFIGURATIONS MUST THEN BE CHANGED BACK TO FIND AIJ
C
           IF(IS .GT. JS) THEN
                ICS=KJCS
                JCS=KICS
           ELSE
                ICS=KICS
                JCS=KJCS
           END IF
C
         END IF
C
C      INITIALISE STORAGE PARAMETERS
C
              IHS=1
          ICHANN=ICHAN1

C
C      LOOP OVER ALL TARGET STATES COMPOSED FROM SET IS
C
      DO 8     IST=1,NSTAT(IS)
         IF (LAMSIN .EQ. 0) THEN
C
C      FIND THE STATE NUMBER AND THUS FIND THE COEFFICIENT AIJ
C
             ISTAT=NSETST(IS,IST)
                AI=AIJ(ISTAT,ICS)
C      WHEN SETS IS, JS ARE THE SAME A PAIR OF CONFIGURATIONS
C      WILL CONTRIBUTE TWICE TO COMPENSATE FOR THE SAVING MADE
C      AS A RESULT OF SYMMETRY. JCSLO=KICS IN DO LOOP 7
C
          IF(IS .EQ. JS  .AND.  ICS .NE. JCS) THEN
               AI2=AIJ(ISTAT,JCS)
          END IF
C
         END IF
C
C      MODIFY LIMITS FOR LOOPING OVER STATES FROM SET JS WHEN IS=JS
C      AND L=LP.  I.E.  ICHAN1=JCHAN1
C
           IF(ICHAN1 .EQ. JCHAN1) THEN
             JSTLO=IST
             JCHANN=ICHANN
           ELSE
             JSTLO=1
             JCHANN=JCHAN1
           END IF
C
C      LOOP OVER STATES COMPOSED FROM THE SET JS
C
      DO 9     JST=JSTLO,NSTAT(JS)
         IF (LAMSIN .EQ. -999) GO TO 91
             JSTAT=NSETST(JS,JST)
                AJ=AIJ(JSTAT,JCS)
              AIAJ=AI*AJ
           IF(IS .EQ. JS  .AND.  ICS .NE. JCS) THEN
              AIAJ=AIAJ+AI2*AIJ(JSTAT,ICS)
           END IF
C
           IF(ABS(AIAJ) .LE. EPS) THEN
             IHS=IHS+NRASQR
             GO TO 91
           END IF
C
C      STORE MATRIX AIAJ*CMAT IN HSMAT
C
C      IF THE CHANNELS ARE REVERSED WHEN CONVERTED TO THE RMATRIX
C      CONVENTION THE MATRIX MUST BE TRANSPOSED  I.E. REVERSE THE
C      INDICES
C
         IF (LP .NE. L  .AND.
     1       ICONCH(JCHANN) .LT. ICONCH(ICHANN)) THEN
            NUP=NRA2
            NPUP=NRA1
         ELSE
            NUP=NRA1
            NPUP=NRA2
         END IF
C
      DO 10      N=1,NUP
C
      DO 11     NP=1,NPUP
C
C      WHEN THE CONTINUUM ANGULAR MOMENTA ARE EQUAL CMAT MUST
C      BE FILLED UP TO A SQUARE
C
          IF(LP .EQ. L  .AND.  NP .LT. N) THEN
             CM=CMAT(N,NP)
C
C      IF THE CHANNELS ARE REVERSED WHEN CONVERTED TO THE
C      RMATRX CONVENTION THE MATRIX MUST BE TRANSPOSED
C
          ELSE IF (LP .NE. L  .AND.
     1             ICONCH(JCHANN) .LT. ICONCH(ICHANN)) THEN
             CM=CMAT(N,NP)
C
          ELSE
             CM=CMAT(NP,N)
          END IF
            HSMAT(IHS)=HSMAT(IHS)+CM*AIAJ
                  IHS=IHS+1
   11 CONTINUE
C
   10 CONTINUE
C
C
   91       JCHANN=JCHANN+1
    9 CONTINUE
C
            ICHANN=ICHANN+1
    8 CONTINUE
      IF (NBUG6 .EQ. 1) THEN
C     WRITE(IWRITE,1002)L,LP,IS,JS,IC,JC,IST,JST
      WRITE(IWRITE,1002)L,LP,IS,JS,ICS,JCS,NSTAT(IS),NSTAT(JS)
      WRITE(IWRITE,1003)(HSMAT(IHSM),IHSM=1,IHS-1)
      END IF
C
    7 CONTINUE
C
    6 CONTINUE
      IF (NBUG6 .EQ. 1) THEN
      WRITE(IWRITE,1002)L,LP,IS,JS
      WRITE(IWRITE,1003)(HSMAT(IHSM),IHSM=1,IHS-1)
      END IF
C
C      TEST THAT THE VSH ARRAY HAS BEEN COMPLETELY USED FOR THIS
C      SET COMBINATION IF THE SETS INTERACTED
C
      IF (LAMSIN .EQ. 0) THEN
      IF (ICOUNT .NE. NVSHEL + 1) THEN
         WRITE (IWRITE,3000) KIS,KJS
         STOP
      END IF
      END IF
C
C      TRANSFER HSMAT TO A BUFFER MATRIX BY MATRIX
C      INCLUDE ASYMPTOTIC COEFFICIENTS FOR EACH LAMBDA
C      AT THE BEGINNING OF THE BUFFER FOR EACH CHANNEL
C      COMBINATION
C
             ICHAN2=ICHANN-1
             JCHAN2=JCHANN-1
             IHS=0
             ICFS=0
C
C      SET BUFFER NUMBER FOR DOUBLE BUFFERING
C
              NBUF=1
C
      DO 12    ICH=ICHAN1,ICHAN2
C
C      CONVERT CHANNEL NUMBER TO RMATRIX CONVENTION
C
            ICHCON=ICONCH(ICH)
C
         IF (ICHAN1 .EQ. JCHAN1) THEN
            JCHLO=ICH
         ELSE
            JCHLO=JCHAN1
         END IF
C
      DO 13    JCH=JCHLO,JCHAN2
C
          IF(JCH .EQ. ICH) THEN
C
C      ADD IN THE ONE-ELECTRON RK INTEGRAL AND THEN
C      ADD THE STATE ENERGY INTO THE DIAGONAL ELEMENTS
C
            IHS1=IHS
            DO 131    IR=1,NRASQR
                   IHS1 =IHS1+1
               HSMAT(IHS1)=HSMAT(IHS1)+R1(IR)
  131       CONTINUE
C
            IHS1=IHS+1
            IST=ISTCH(ICH)
              E=ENAT(IST)
            DO 132    IR=1,NRA1
               HSMAT(IHS1)=HSMAT(IHS1)+E
                     IHS1 = IHS1+NRA1+1
  132       CONTINUE
C
          END IF
C
C      CALCULATE THE RECORD NUMBER SO THAT THE MATRICES ARE
C      STORED IN SEQUENCE IN THE RMATRIX CONVENTION
C
         JCHCON=ICONCH(JCH)
         IF(JCHCON .LT. ICHCON) THEN
           IJREC=NSREC+NCH*(JCHCON-1)-(JCHCON*(JCHCON-1))/2+ICHCON
         ELSE
           IJREC=NSREC+NCH*(ICHCON-1)-(ICHCON*(ICHCON-1))/2+JCHCON
         END IF
C
         IF (NBUF .EQ. 1) THEN
C
C      FILL BUF1D
C
              ICFS=ICFS+1
              DO 141  IB=1,LAMAX
                 BUF1D(IB)=CF(ICFS,IB)
 141          CONTINUE
C
           DO 142   IB=1+LAMAX,NRASQR+LAMAX
              IHS=IHS+1
              BUF1D(IB)=HSMAT(IHS)
 142       CONTINUE
C
           WRITE(JBUFD,REC=IJREC)BUF1D
           NBUF=2
C
         ELSE
C
C      FILL BUF2D
C
              ICFS=ICFS+1
              DO 143   IB=1,LAMAX
                 BUF2D(IB)=CF(ICFS,IB)
  143         CONTINUE
           DO 144   IB=1+LAMAX,NRASQR+LAMAX
              IHS=IHS+1
              BUF2D(IB)=HSMAT(IHS)
  144      CONTINUE
C
           WRITE(JBUFD,REC=IJREC) BUF2D
           NBUF=1
C
        END IF
C
   13 CONTINUE
C
   12 CONTINUE
C
    5 CONTINUE
C
    4 CONTINUE
C
    3 CONTINUE
C
    2 CONTINUE
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE WRITE1
C
C      WRITE SECTION OF ASYMPTOTIC FILE WHICH IS INDEPENDENT OF
C      LRGL, SPIN AND PARITY
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      INCLUDE 'PARAM'
C
      PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2)
      PARAMETER(IPDIM1=MZSPN)
      PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS)
      PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS)
      PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4)
      PARAMETER(IDIM2=ILDIM3+ILDIM1-1)
C
      COMMON/DISC  /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD,
     1              ITAPE1,ITAPE3
      COMMON/RAD   /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1),
     1              NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM
      COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2),
     1              LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1),
     2              NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1)
      COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1)
      COMMON/STG1  /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2)
C
C     WRITE(ITAPE3) NELC,NZ,LRANG2,LAMAX,NAST,RA,BSTO,LRANG1,
C    1              (MAXNHF(I),I=1,LRANG2),(MAXNC(I),I=1,LRANG1),NRANG2
      WRITE(ITAPE3) NELC,NZ,LRANG2,LAMAX,NAST,RA,BSTO
C
C****
C      NRANG2 SHOULD BE REPLACED BY NOTERM
C
      WRITE(ITAPE3)(ENAT(I),I=1,NAST)
      WRITE(ITAPE3)(LL(I),I=1,NAST)
      WRITE(ITAPE3)(LSPN(I),I=1,NAST)
C     WRITE(ITAPE3)(LPTY(I),I=1,NAST)
      WRITE(ITAPE3)((COEFF(I,L),I=1,3),L=1,LRANG2)
C
      RETURN
      END